Fortune eigensolve算法Python移植异常求排查及现有实现问询

问题描述 #

是否存在Fortune eigensolve算法的现成Python实现？我尝试手动移植该算法，但始终无法得到正确结果，仅找到描述该算法的论文。运行自己的代码后，输出全是重复的特征值，不清楚问题出在哪里。

我的实现代码 #

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

# 创建广义伴随矩阵
def gen_companion(p,S):
    # 生成拉格朗日系数
    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)
    # 生成拉格朗日矩阵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)
    # 创建矩阵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)

内容的提问来源于stack exchange，提问作者Ragon