计算多项式的根

Question

给定一个矩阵 A，每列代表一个多项式。如何高效地计算每个多项式的根，不使用循环？

- Razer

Answer 1

这里有三种方法的比较：

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中避免使用循环的旧建议只不过是老掉牙。它已经不再适用于R2009及以上版本的MATLAB。

向量化仍然可能是一个好东西，但并非总是如此。就像这种情况一样：分析将告诉您，大部分时间都花费在计算块对角矩阵的特征值上。eig底层算法的规模为N³（是的，是三），而且它无法以任何方式利用稀疏矩阵（例如这个块对角矩阵），使得这种方法在这种特定情况下成为一个糟糕的选择。

在这里，循环是你的朋友 ^_^

现在，这当然是基于伴随矩阵的eig()，这是一种计算所有根的简单有效的方法。当然，还有许多其他计算多项式根的方法，每种方法都有其自己的优缺点。有些方法速度更快，但在某些根为复数时效果不佳。其他方法速度更快，但需要对每个根进行相当好的初始估计等等。大多数其他求根方法通常更加复杂，这就是为什么我在这里没有涉及它们的原因。 这里提供了一个很好的概述，这里提供了更深入的概述，以及一些MATLAB代码示例。

如果你聪明的话，只有在未来几周内每天需要进行数百万次计算时才应该深入研究这些内容，否则这并不值得投资。

如果你更聪明，你会意识到这无疑会在某个时候回到你身上，所以做这件事是值得的。

如果你是一名学者，你将掌握所有的根查找方法，这样你就拥有了一个巨大的工具箱，可以在每次出现新任务时选择最好的工具。甚至可以自己构思新的方法 :)

- Rody Oldenhuis

---

非常好的分析！+1。 - jkt

我想你刚刚为我节省了几天的计算时间:D 你知道如何加快找到多项式导数根以便找到多项式最大/最小值的组合任务吗？大多数情况下，只有一个根是重要的 :/ - geometrikal

@geometrikal：这样做的方式是一样的，但是导数的系数矩阵将会是 Cprime = (numel(C)-1:-1:1) .* C(1:end-1);，其中 C 是你原始多项式的系数。顺便问一下，你所说的“组合任务”具体指什么？ - Rody Oldenhuis

@RodyOldenhuis 我需要找到一个复多项式的最大值。所以我要求导，找到它的根，然后在原始多项式中评估每个根来找到最大值。我只是想知道是否有一些加速方法，可以在不浪费时间在其他根上的情况下找到对应于最大值的根。或者你知道是否有一种方法可以通过仅评估几个有效数字的根来加速eig吗？如果是这样，我会提出这个问题。 - geometrikal

@geometrikal：听起来像是优化理论 :) 据我所知，多项式没有已知的方法可以在事先确定最小/最大值。一些想法：1）奇次多项式或正实数首项系数的多项式没有最大值，2）二阶导数测试可能会节省您一些计算，3）eigs()允许更好地控制精度，但仅适用于大型稀疏问题--进行分析以查看它是否适合您的需求。 - Rody Oldenhuis

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));

- Stewie Griffin

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     

- jkt

---

如果有20列或100列怎么办？使用A=blkdiag(compan(p1),compan(p2),compan(p3),compain(p4),...)不是一个好的方法。 - Stewie Griffin

不错的开始，不过你还需要想出一个能处理任意大小 A 的版本。此外，我会使用稀疏矩阵，并重新整理输出，使得至少有一个维度等于多项式的数量。 - Rody Oldenhuis

@RodyOldenhuis 谢谢。我只是试图阐述这个概念。 - jkt

@RobertP. 你说得对，对于大量列来说它并不是很高效。但是对于少量列来说，它可能比 arrayfun 更快。此外，我不确定，但 MATLAB 可能也会利用矩阵的分块对角形式，这可能会更快。不过这只是一个假设。 - jkt

也许使用这种方法处理10个多项式的批次是最快的。 - Daniel