Вычислить корни кратных многочленов

Question

Учитывая матрицу A который представляет полиномы в каждом столбце. Как эффективно вычислить корни каждого многочлена без петель?

Answer 1

Вот сравнение между 3 методами:

1. Простой цикл по всем строкам, используя roots на каждом ряду.
2. Полностью петлевой подход, основанный на идее YBE об использовании блочно-диагональной матрицы с использованием sparse в качестве промежуточного
3. Простой цикл по всем строкам, но на этот раз с использованием "встроенного" кода из roots,

Код:

%// The polynomials
m = 15;
n = 8;
N = 1e3;

X = rand(m,n);


%// Simplest approach
tic
for mm = 1:N

    R = zeros(n-1,m);
    for ii = 1:m
        R(:,ii) = roots(X(ii,:));
    end

end
toc


%// Completely loopless approach
tic
for mm = 1:N

    %// Indices for the scaled coefficients
    ii = repmat(1:n-1:m*(n-1), n-1,1);
    jj = 1:m*(n-1);

    %// Indices for the ones
    kk = bsxfun(@plus, repmat(2:n-1, m,1), (n-1)*(0:m-1).');  %'
    ll = kk-1;

    %// The block diagonal matrix
    coefs = -bsxfun(@rdivide, X(:,2:end), X(:,1)).';  %'
    one   = ones(n-2,m);
    C = full(sparse([ii(:); kk(:)], [jj(:); ll(:)],...
        [coefs(:); one(:)]));

    %// The roots
    R = reshape(eig(C), n-1,m);

end
toc


%// Simple loop, roots() "inlined"
tic    
R = zeros(n-1,m);
for mm = 1:N

    for ii = 1:m            
        A = zeros(n-1);
        A(1,:) = -X(ii,2:end)/X(ii,1);
        A(2:n:end) = 1;
        R(:,ii) = eig(A);            
    end

end
toc

Результаты, достижения:

%// m=15, n=8, N=1e3:
Elapsed time is 0.780930 seconds. %// loop using roots()
Elapsed time is 1.959419 seconds. %// Loopless
Elapsed time is 0.326140 seconds. %// loop over inlined roots()

%// m=150, n=18, N=1e2:
Elapsed time is 1.785438 seconds. %// loop using roots()
Elapsed time is 110.1645 seconds. %// Loopless
Elapsed time is 1.326355 seconds. %// loop over inlined roots()

Конечно, ваш пробег может отличаться, но общее сообщение должно быть ясным: старый совет избегать циклов в MATLAB - это просто: СТАРЫЙ. Это больше не относится к версиям MATLAB R2009 и выше.

Хотя векторизация может быть хорошей вещью, но не всегда. Как и в этом случае: профилирование скажет вам, что большая часть времени уходит на вычисление собственных значений для блочно-диагональной матрицы. Алгоритм, лежащий в основе eig масштабируется как N³ (да, это три), плюс он никоим образом не может использовать разреженные матрицы (например, эту блочно-диагональную), что делает этот подход неудачным выбором в данном конкретном контексте.

Петли твой друг здесь ^_^

Теперь это, конечно, основано на eig() сопутствующей матрицы, которая является хорошим и простым методом для вычисления всех корней за один раз. Есть, конечно, еще много способов вычисления корней полиномов, каждый из которых имеет свои преимущества / недостатки. Некоторые из них намного быстрее, но не так хороши, когда некоторые корни сложны. Другие намного быстрее, но требуют достаточно хорошей начальной оценки для каждого корня и т. Д. Большинство других методов поиска корней обычно намного сложнее, поэтому я оставил их здесь.

Вот хороший обзор, более подробный обзор, а также примеры кода MATLAB.

Если вы сообразительны, вам следует погрузиться в этот материал только в том случае, если вам нужно делать это вычисления миллионы раз в день, по крайней мере, в течение следующих нескольких недель, иначе это просто не стоит вложений.

Если вы будете умнее, вы поймете, что это, несомненно, вернется к вам в какой-то момент, так что в любом случае это стоит сделать.

И если вы академик, вы овладеете всеми методами поиска корней, и у вас будет огромный набор инструментов, чтобы вы могли выбрать лучший инструмент для работы, когда появится новая работа. Или даже придумайте свой собственный метод:)

Answer 2

Ты можешь использовать arrayfun в комбинации с roots, который даст вам результаты с точки зрения клеточных массивов.

n = size(A,2);
t = arrayfun(@(x)roots(A(:,x)), 1:n, 'UniformOutput', 0);

Вы можете использовать cell2mat преобразовать его в матрицу. Или: r = cell2mat(t), или же

r = cell2mat(arrayfun(@(x)roots(A(:,x)), 1:n, 'UniformOutput', 0));

Answer 3

Практически что roots делает, чтобы найти собственные значения матрицы компаньона.

roots(p) = eig(compan(p))

Итак, вот мой пример, который строит блочно-диагональную матрицу из сопутствующих матриц каждого полинома, а затем находит собственные значения блочно-диагональной матрицы.

>> p1=[2 3 5 7];
>> roots(p1)

ans =

  -0.0272 + 1.5558i
  -0.0272 - 1.5558i
  -1.4455          

>> eig(compan(p1))

ans =

  -0.0272 + 1.5558i
  -0.0272 - 1.5558i
  -1.4455          

>> p2=[1 2 9 5];
>> roots(p2)

ans =

  -0.6932 + 2.7693i
  -0.6932 - 2.7693i
  -0.6135          

>> p3=[5 1 4 7];
>> roots(p3)

ans =

   0.3690 + 1.1646i
   0.3690 - 1.1646i
  -0.9381          

>> A=blkdiag(compan(p1),compan(p2),compan(p3))

A =

   -1.5000   -2.5000   -3.5000         0         0         0         0         0         0
    1.0000         0         0         0         0         0         0         0         0
         0    1.0000         0         0         0         0         0         0         0
         0         0         0   -2.0000   -9.0000   -5.0000         0         0         0
         0         0         0    1.0000         0         0         0         0         0
         0         0         0         0    1.0000         0         0         0         0
         0         0         0         0         0         0   -0.2000   -0.8000   -1.4000
         0         0         0         0         0         0    1.0000         0         0
         0         0         0         0         0         0         0    1.0000         0

>> eig(A)

ans =

  -0.0272 + 1.5558i
  -0.0272 - 1.5558i
  -1.4455          
  -0.6932 + 2.7693i
  -0.6932 - 2.7693i
  -0.6135          
   0.3690 + 1.1646i
   0.3690 - 1.1646i
  -0.9381