algorithm - 使用伴随矩阵求根

Question

我想找到单变量多项式的所有实根。例如，我可以使用 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形式。三角矩阵的特征值列在对角线上，特征值问题得到解决。”< br/>...等等，什么？我试过:

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

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

拜托，谁能给我解释一下？

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

PPS:还有http://adorio-research.org/wordpress/?p=184 (虽然不确定它与第一种方法有何不同...)

Answer 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 算法期间保留了这种形式，您可以使用带有 Givens 旋转的 QR 分解，其中仅发生靠近对角线的旋转。

在未移位的 QR 算法中，您应该至少在右下方的 3x3 block 中观察到明显的变化。

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

只要任何次对角线元素几乎为零，就会发生矩阵 split 。如果拆分发生在最后一行，则最后一个对角线条目是移位矩阵 (A-s*I) 的特征值。如果拆分将最后一个 2x2 block 分开，那么可以很容易地确定其特征值，这又是移位矩阵的特征值。

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

• 附录一

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

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

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

较大的对角线 block 将出现在多个或聚类特征值中。

附录二

是线

   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)