1

我想找到一个单变量多项式的所有实根。例如，我可以使用 Jenkins-Traub 算法，但我想学习如何使用伴随矩阵来解决它。

我知道如何将多项式转换为伴随矩阵，并且我找到了一个执行 QR 分解的脚本：http: //quantstart.com/articles/QR-Decomposition-with-Python-and-NumPy

这就是我迷失的地方：下一步该做什么？我想我必须计算多个分解，但是当我这样做时，我总是得到相同的结果（显然）。我还读到，首先将伴随矩阵转换为 Hessenberg 形式可能很有用 - 但是如何？然后是“轮班”——它们是什么？

我还找到了http://www.nr.com/webnotes/nr3web17.pdf但由于我什么都不懂，我想知道是否有更简单的方法（即使速度较慢或不太稳定）。

换句话说：阅读http://en.wikipedia.org/wiki/QR_algorithm

  • “让 A 成为我们想要计算其特征值的实矩阵”
    好的，那是我的伴生矩阵，对吗？

  • “我们计算来自第一个链接的 QR 分解 Ak=QkRk”，对吗
    ？Q, R = householder(A)

  • “然后我们形成 Ak+1 = RkQk”
    很容易，只需将 R 和 Q 相乘

  • “在一定条件下，[2] 矩阵 Ak 收敛到一个三角矩阵，即 A 的 Schur 形式。三角矩阵的特征值列在对角线上，特征值问题就解决了。”
    ……等等，什么？我试过了：

    for i in range(100):
    Q, R = householder(A)
    A = mult_matrix(R, Q)

但似乎没有任何进展，我什至看不到任何接近正确根源的数字。

请问，谁能给我解释一下？

PS：我不想盲目地使用 LAPACK 或类似的东西，因为我想了解它是如何工作的，至少在非常简化的方面。

PPS：还有http://adorio-research.org/wordpress/?p=184（不知道它与第一种方法有什么不同，虽然......）

4

1 回答 1

1

如果您的伴随矩阵是映射q(x)到的线性多项式运算的系数变换矩阵x*q(x) mod p(x)

哪里p(x)=x^(n+1)+p_n*x^n+...+p_1*x+p_0。

明确地，A 的形状为

0 0 0 ... 0 -p_0
1 0 0 ... 0 -p_1
0 1 0 ... 0 -p_2
. . . ... . . . 
0 0 0 ... 1 -p_n

已经是 Hessenberg 形式了。由于在 QR 算法期间保留了这种形式，因此您可以将 QR 分解与 Givens 旋转一起使用，其中仅发生接近对角线的旋转。

在未移位的 QR 算法中，您至少应该在右下角的 3x3 块中观察到明显的发展。

如果你取右下角 2x2 块的一个特征值并从每个对角线元素中减去它，那么你会得到按照 Francis 的移位 QR 算法。位移 s，减去的数字，是最小特征值的当前最佳估计。请注意，很可能，您现在离开了实域并且必须使用复杂的矩阵条目进行计算。您必须将 shift s 保存在内存中，并在后续步骤中添加任何新的 shift，并将组合的 shift 添加回找到的任何特征值。

只要任何对角线条目实际上为零，就会发生矩阵拆分。如果拆分发生在最后一行，则最后一个对角线条目是移位矩阵 (As*I) 的特征值。如果拆分将最后一个 2x2 块分开，则可以轻松确定其特征值，这些特征值又是移位矩阵的特征值。

如果拆分发生在其他任何地方，则 QR 算法将递归地分别应用于对角块。

  • 附录一

算法所有变体的收敛性由对角线下方的条目收敛到零来衡量。

基本算法在那些次对角条目中具有几何收敛性，(i,i-1) 处的条目的几何因子是位置 i-1 和 i 处特征值大小的分数。只要有大的跳跃，就会很快达到零。

共轭复特征值大小相同，因此算法会在对角线上产生一个 2x2 的块，求解该块的二次特征方程得到对应的特征值。

对于多个或聚集的特征值，会出现更大的对角块。

附录二

这是线

   alpha = -cmp(x[0],0) * norm(x)

在户主程序中。在大多数情况下，x[0]不会完全为 0，因此这会产生一个符号。然而，在伴随矩阵的情况下，x[0]==0通过构造，因此没有产生符号，alpha 得到错误的值 0。将其更改为更行人

    alpha = norm(x)
    if x[0] < 0: alpha = -alpha

它运作良好。

def companion_matrix(p):
    n=len(p)-1
    C=[[float(i+1 == j) for i in xrange(n-1)] for j in xrange(n)]
    for j in range(n): C[j].append(-p[n-j]/p[0])
    return C


def QR_roots(p):
    n=len(p)-1
    A=companion_matrix(p)
    for k in range(10+n):
        print "step: ",k+1," after ",(k+1)*(5+n), "iterations"
        for j in range(5+n):
            Q,R=householder(A)
            A=mult_matrix(R,Q)
        print("below diagonal")
        pprint([ A[i+1][i] for i in range(n-1) ])
        print("diagonal")
        pprint([ A[i][i] for i in range(n) ])
        print("above diagonal")
        pprint([ A[i][i+1] for i in range(n-1) ])


p=[ 1, 2, 5, 3, 6, 8, 6, 4, 3, 2, 7]
QR_roots(p)

#for a case with multiple roots at 1 and 4
#p= [1,-11,43,-73,56,-16]
#QR_roots(p)
于 2013-12-11T09:22:28.820 回答