是否存在财富特征解算算法的实现?

目前,我正在try 将其移植到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)

推荐答案

本文中的公式为

    l[i] = P(s[i]) / Q_S'(s[i])  =  P(s[i]) / prod( (s[i]-s[j]) for i!=j)

这可以实现为

Q = np.poly(S)
dQ = np.polyder(Q)
l = np.polyval(P,S)/np.polyval(dQ,S)

注意,纸上的垂直线是错误的,或者至少不应该表示绝对值.

正如我在 comments 中指出的,QR算法的这个变体是前阿尔法版本.即使Francis推荐用于实现的第一个变型也使用移位来加速最小特征值的收敛,从而使矩阵分裂,收缩到与最小特征空间垂直的子空间.

如果不太关心计算时间,可以直接使用这个原始变量.它将收敛到允许读出本征值的正规形式.这个标准形的本征值按绝对值递减的顺序排列.特别是在寻根器的后期阶段,伴随矩阵已经接近对角线.如果对角线条目的顺序不正确,QR算法将使用较小的旋转来交换它们,从而进行大量旋转.这种开销不会影响根/本征值的精度,因此应该避免.这是您的实现自动完成的.

我想让您注意np.diag函数和带有常量元素的列表.然后可以将伴随矩阵构造为

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

进行这些更改并始终返回完整结果时,返回值为eigensolve,与预期值相同

[20.      +0.j 10.000002+0.j  5.000003+0.j  4.000001+0.j]

Python相关问答推荐

Python无法在已导入的目录中看到新模块

查找下一个值=实际值加上使用极点的50%

从DataFrame.apply创建DataFrame

剧作家Python:expect(locator).to_be_visible()vs locator.wait_for()

Odoo -无法比较使用@api.depends设置计算字段的日期

GL pygame无法让缓冲区与vertextPointer和colorPointer一起可靠地工作

根据条件将新值添加到下面的行或下面新创建的行中

如何让剧作家等待Python中出现特定cookie(然后返回它)?

无法定位元素错误404

给定高度约束的旋转角解析求解

在matplotlib中删除子图之间的间隙_mosaic

使用BeautifulSoup抓取所有链接

如何在PySide/Qt QColumbnView中删除列

Python—压缩叶 map html作为邮箱附件并通过sendgrid发送

如何从pandas DataFrame中获取. groupby()和. agg()之后的子列?

不允许 Select 北极滚动?

在numpy数组中寻找楼梯状 struct

使用tqdm的进度条

如何反转一个框架中列的值?

如何提高Pandas DataFrame中随机列 Select 和分配的效率?