计算多个多项式的根

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中避免循环的旧建议只是:OLD.它不再适用于MATLAB版本R2009及更高版本.

尽管如此,矢量化仍然是一件好事,但肯定并非总是如此.在这种情况下:分析将告诉您大部分时间花在计算块对角矩阵的特征值上.基于eig比例的算法N³(是,即三是),加上它不能以任何方式利用稀疏矩阵(比如这个块对角线),使得该方法在这个特定的上下文中是一个糟糕的选择.

循环是你的朋友^ _ ^

现在,这当然是基于eig()伴随矩阵,这是一次计算所有根的一个很好而且简单的方法.当然还有更多的方法来计算多项式的根,每个方法都有自己的优点/缺点.有些速度要快得多,但是当一些根很复杂时,它们就不那么好了.其他的更快,但需要对每个根进行相当好的初始估计等.大多数其他的根寻找方法通常要复杂得多,这就是为什么我把它们留在这里.

这是一个很好的概述,这里有一个更深入的概述,以及一些MATLAB代码示例.

如果你很聪明,如果你需要在接下来的几周内每天进行数百万次这样的计算,你应该只深入研究这些材料,否则,这不值得投资.

如果你更聪明,你会发现这无疑会在某些时候回到你身边,所以无论如何都要做到这一点.

如果你是一名学者,你掌握了所有的寻根方法,这样你就拥有了一个巨大的工具箱,这样你就可以在新工作出现时为工作挑选最好的工具.或者甚至梦想你自己的方法:)