计算多项式的根 [英] Calculate roots of multiple polynomials

问题描述

给出一个表示每列中多项式的矩阵 A .如何无循环高效地计算每个多项式的根?

Given a matrix A that represents polynomials in each column. How can the roots of each polynomial be calculated efficiently without loops?

推荐答案

以下是3种方法之间的比较:

Here's a comparison between 3 methods:

1. 在所有行中进行简单循环，在每行上使用 roots .
2. 基于YBE使用块对角矩阵的想法，这是一种完全无环的方法，使用 sparse 作为中间
3. 一个简单的循环遍历所有行，但这次使用 roots 中的内联"代码.

1. A simple loop through all the rows, using roots on each row.
2. A completely loopless approach, based on YBE's idea of using a block-diagonal matrix, using sparse as an intermediate
3. A simple loop through all the rows, but this time using "inlined" code from 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及更高版本.

Of course, your mileage may vary, but the general message should be clear: The old advice of avoiding loops in MATLAB is just that: OLD. It just no longer applies to MATLAB versions R2009 and up.

虽然矢量化仍然可以是一件好事，但并非总是如此.就像在这种情况下:剖析将告诉您，大部分时间都花在计算块对角矩阵的特征值上. eig 底层的算法可缩放为N³(是的，这是一个 3 )，而且它不能以任何方式利用稀疏矩阵(例如这种对角线的方法)，因此在特定情况下该方法不是一个好的选择.

Vectorization can still be a good thing though, but certainly not always. As in this case: profiling will tell you that most time is spent on computing the eigenvalues for the block-diagonal matrix. The algorithm underlying eig scales as N³ (yes, that is a three), plus it cannot take advantage of sparse matrices in any way (like this block-diagonal one), making the approach a poor choice in this particular context.

循环是您的朋友在这里^ _ ^

Loops are your friend here ^_^

现在，这当然是基于伴随矩阵的 eig()，这是一种一次性计算所有根的好方法.当然，还有更多种计算多项式根的方法，每种方法都有其自身的优点/缺点.有些速度更快，但是当一些根源很复杂时就不好了.其他方法则要快得多，但是每个根等都需要一个相当好的初始估计值.其他大多数根查找方法通常要复杂得多，这就是为什么我在这里省略了这些原因.

Now, this is of course based on eig() of the companion matrix, which is a nice and simple method to compute all roots in one go. There are of course many more methods to compute roots of polynomials, each with their own advantages/disadvantages. Some are a lot faster, but aren't so good when a few of the roots are complex. Others are a lot faster, but need a fairly good initial estimate for every root, etc. Most other rootfinding methods are usually a lot more complicated, which is why I left these out here.

此处是一个很好的概述，而

Here is a nice overview, and here is a more in-depth overview, along with some MATLAB code examples.

如果您很聪明，则至少在接下来的几周内每天需要进行数百万次此计算时，才应深入研究此材料，否则，这是不值得的投资.

If you're smart, you should only dive into this material if you need to do this computation millions of times on a daily basis for at least the next few weeks, otherwise, it's just not worth the investment.

如果您更聪明，您会意识到这无疑会在某个时候再次出现，因此还是值得做的.

If you're smarter, you'll recognize that this will undoubtedly come back to you at some point, so it's worthwhile to do anyway.

如果您是一名学者，那么您将掌握所有寻根方法，因此您将拥有一个庞大的工具箱，因此，每当有新工作出现时，您都可以为该工作选择最佳工具.甚至幻想自己的方法:)

And if you're an academic, you master all the root-finding methods so you'll have a giant toolbox, so you can pick the best tool for the job whenever a new job comes along. Or even dream up your own method :)

这篇关于计算多项式的根的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！