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

Question

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

Answer 1

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

1. Простой цикл по всем строкам, используя roots на каждой строке.
2. Полностью loopless подход, основанное на идее о 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 это только: OLD. Это больше не распространяется на версии MATLAB R2009 и выше.

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

Loops ваш друг здесь^_^

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

Here - хороший обзор, а here - это более подробный обзор, а также некоторые примеры кода 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