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

Вот сравнение между 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.

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

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

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

Вы можете использовать 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));

Практически то, что делает 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