财富算法特征求解

Question

是否存在 Fortunes eigensolve 算法的实现？

目前我正在尝试将其移植到Python，但到目前为止它没有提供正确的结果。 我唯一找到的是他描述算法的论文here，但我不明白我的算法哪里出了问题。

import numpy as np
from sympy import Poly, symbols
lamda = symbols('lamda')

#Creates generalized companion matrix
def gen_companion(p,S):
    #create lagrange coeffs
    L_i = []
    for i in range(len(S)):
        den_s = 1
        for j in range(len(S)):
            if (not (i == j)):
                den_s *= p(S[i]-S[j])
        p_s = p(S[i])/den_s
        L_i.append(p_s)
    #create lagrange matrix L 
    dim = S.shape[0]
    L = np.zeros((dim,dim),np.complex64)
    for i in range(dim):
        for j in range(dim):
            L[j,i] = np.abs(L_i[i])
    B = np.zeros((dim,dim),np.complex64)
    #Create Matrix S
    for i in range(dim):
        B[i,i] = S[i]
    return B-L
    


def eigensolve(p,n_iterations):
    A = np.polynomial.polynomial.polycompanion(coeffs).astype(np.complex64)
    S = qr_algorithm(A,n_iterations)
    for _ in range(n_iterations):
        L = gen_companion(p,S)
        S = qr_algorithm(L,n_iterations)
    return S[0]


def qr_algorithm(A,n_iterations):
    for _ in range(n_iterations):
        Q,R = np.linalg.qr(A)
        A = R @ Q
    return np.array([A[i,i] for i in range(A.shape[0])])

p = Poly((lamda - 10)*(lamda-20)*(lamda-4)*(lamda-5),lamda)
coeffs = np.array(list(reversed(p.all_coeffs())))


n_iterations = 10
v = eigensolve(p,n_iterations)
print("")

我只得到了与算法第二部分中输入的相同的特征值：

(20.153511+0j)
(20.153511+0j)
(20.153511+0j)
(20.153511+0j)
(20.153511+0j)

Answer 1

论文中的公式是

l[i] = P(s[i]) / QS'(s[i])

这可以实现为

QS = np.poly(S)
dQS = np.polyder(QS)
l = np.polyval(P,S)/np.polyval(dQS,S)

请注意，论文中的垂直线是错误的，或者至少不应该表示绝对值。

正如评论中所指出的，QR 算法的这种变体是 pre-alpha 版本。即使建议实现的第一个变体也使用移位来加速最小特征值的收敛，从而加速矩阵的分裂，压缩到与最小特征空间正交的子空间。

如果不太关心计算时间，也可以使用这种原始变体。它将收敛到允许读取特征值的范式。该范式的特征值按绝对值降序排列。特别是在寻根的后期阶段，伴随矩阵将已经接近对角线。如果对角线条目的顺序不正确，QR 算法将使用小旋转来交换它们，从而进行大量旋转。这种开销对根/特征值的精度没有贡献，因此应该避免。您的实现会自动执行此操作。

我想请您注意

np.diag

函数和包含常量元素的列表。然后可以将伴随矩阵构造为

np.diag(S) - np.array(dim*[l])