如何解决Python中的多项式特征值？

Question

在我的Python代码中，我想解决多项式特征值问题：

A0 + lambda*A1 + lambda^2*A2 + lambda^3*A3 + .... = 0

An

是密集的矩阵，

lambda

是一个常数。在MATLAB中，可以使用

polyeig函数解决此问题。 Scipy似乎没有等效功能。到目前为止，我能想到的唯一方法就是形成相应的伴侣矩阵。这会产生同等的线性特征值问题，可以给出现有的Scipy求解器，但是它更大，我相信它可能会非常不适。 是否有人建议可以解决此问题的现有，开源或自由储存库？我会对可以通过F2PY或C/C ++库链接到通过Cython链接的Fortran库感到满意。

Edit：对于有兴趣解决python中非线性特征值问题的任何人，我为解决这个问题而写的代码都可以找到。请注意，我处理非线性特征值问题的更一般情况（从某种意义上说，它对lambda具有非线性依赖性）。要了解该方法，请阅读代码注释中提到的论文。

在这里，我对python的实施的实施。作为测试案例，我使用了与mathworks的相同示例。 Scaling

：特征向量与MATLAB返回的特征向量不同（因为将其定义为缩放因子）。我不确定MATLAB使用了哪种归一化，因此在这种情况下，我只是按每个向量的最大值进行标准化。

Answer 1

sorting

：MATLAB不会对特征值进行排序。在那个示例中，

polyeig

似乎以与MATLAB相同的顺序返回特征值。您可以在脚本中取消分类以对特征值进行排序。 scipy

这个讨论指出了将多项式特征值问题转变为广义特征值问题的一般方法，以后可以使用

Scipy的线性代数函数来解决该问题。希望这有帮助！

P：请注意，在上述消息的答案之一中，作者提到了错字。此外，有人补充说，这也是

import numpy as np
from scipy import linalg

def polyeig(*A):
    """
    Solve the polynomial eigenvalue problem:
        (A0 + e A1 +...+  e**p Ap)x=0 

    Return the eigenvectors [x_i] and eigenvalues [e_i] that are solutions.

    Usage:
        X,e = polyeig(A0,A1,..,Ap)

    Most common usage, to solve a second order system: (K + C e + M e**2) x =0
        X,e = polyeig(K,C,M)

    """
    if len(A)<=0:
        raise Exception('Provide at least one matrix')
    for Ai in A:
        if Ai.shape[0] != Ai.shape[1]:
            raise Exception('Matrices must be square')
        if Ai.shape != A[0].shape:
            raise Exception('All matrices must have the same shapes');

    n = A[0].shape[0]
    l = len(A)-1 
    # Assemble matrices for generalized problem
    C = np.block([
        [np.zeros((n*(l-1),n)), np.eye(n*(l-1))],
        [-np.column_stack( A[0:-1])]
        ])
    D = np.block([
        [np.eye(n*(l-1)), np.zeros((n*(l-1), n))],
        [np.zeros((n, n*(l-1))), A[-1]          ]
        ]);
    # Solve generalized eigenvalue problem
    e, X = linalg.eig(C, D);
    if np.all(np.isreal(e)):
        e=np.real(e)
    X=X[:n,:]

    # Sort eigenvalues/vectors
    #I = np.argsort(e)
    #X = X[:,I]
    #e = e[I]

    # Scaling each mode by max
    X /= np.tile(np.max(np.abs(X),axis=0), (n,1))

    return X, e

if __name__=='__main__':
    M = np.diag([3,1,3,1])
    C = np.array([[0.4 , 0 , -0.3 , 0], [0 , 0  , 0 , 0], [-0.3 , 0 , 0.5 , -0.2 ], [ 0 , 0 , -0.2 , 0.2]])
    K = np.array([[-7  , 2 , 4    , 0], [2 , -4 , 2 , 0], [4    , 2 , -9  , 3    ], [ 0 , 0 , 3    , -3]])
    X,e = polyeig(K,C,M)
    print('X:\n',X)
    print('e:\n',e)
    # Test that first eigenvector and value satisfy eigenvalue problem:
    s = e[0];
    x = X[:,0];
    res = (M*s**2 + C*s + K).dot(x) # residuals
    assert(np.all(np.abs(res)<1e-12))

也将在引擎盖下做的事情。