计算多个多项式的根

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()，这是一次计算所有根的好而简单的方法。当然还有更多计算多项式根的方法，每种方法都有自己的优点/缺点。有些速度要快得多，但当一些根很复杂时就不那么好了。其他的要快得多，但需要对每个根等进行相当好的初始估计。大多数其他寻根方法通常要复杂得多，这就是我在这里省略这些的原因。

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