用Python查找矩阵特征多项式的零点

Question

给定一个N x N 对称矩阵C 和一个N x N 对角矩阵I，求方程det(λI-C)=0 的解。换句话说，要找到C 的（广义）特征值。

我知道如何在 MATLAB 中使用内置函数解决这个问题的几种方法：

第一种方式：

function lambdas=eigenValues(C,I)
    syms x;
    lambdas=sort(roots(double(fliplr(coeffs(det(C-I*x))))));

第二种方式：

[V,D]=eig(C,I);

但是，我需要使用 Python。 NumPy 和 SymPy 中也有类似的功能，但是根据文档（numpy、sympy），它们只采用一个矩阵 C 作为输入。虽然，结果与 Matlab 产生的结果不同。此外，SymPy 生成的符号解决方案也无济于事。也许我做错了什么？如何解决？

---

示例

• MATLAB：

%输入

I =

     2     0     0
     0     6     0
     0     0     5
C =

     4     7     0
     7     8    -4
     0    -4     1

[v,d]=eig(C,I)

%结果

v =

    -0.3558   -0.3109   -0.5261
     0.2778    0.1344   -0.2673
     0.2383   -0.3737    0.0598


d =

       -0.7327    0         0
        0    0.4876         0
        0         0    3.7784

• Python 3.5：

%输入

I=np.matrix([[2,0,0],
                 [0,6,0],
                 [0,0,5]])

C=np.matrix([[4,7,0],[7,8,-4],[0,-4,1]])


np.linalg.eigh(C)

%结果

(array([-3., 1.91723747, 14.08276253]),


matrix(

       [[-0.57735027,  0.60061066, -0.55311256],

        [ 0.57735027, -0.1787042 , -0.79670037],

        [ 0.57735027,  0.77931486,  0.24358781]]))

Answer 1

至少如果I 具有正对角线条目，您可以简单地求解转换后的系统：

# example problem
>>> A = np.random.random((3, 3))
>>> A = A.T @ A
>>> I = np.identity(3) * np.random.random((3,))

# transform 
>>> J = np.sqrt(np.einsum('ii->i', I))
>>> B = A / np.outer(J, J)

# solve
>>> eval_, evec = np.linalg.eigh(B)

# back transform result
>>> evec /= J[:, None]

# check
>>> A @ evec
array([[ -1.43653725e-02,   4.14643550e-01,  -2.42340866e+00],
       [ -1.75615960e-03,  -4.17347693e-01,  -8.19546081e-01],
       [  1.90178603e-02,   1.34837899e-01,  -1.69999003e+00]])
>>> eval_ * (I @ evec)
array([[ -1.43653725e-02,   4.14643550e-01,  -2.42340866e+00],
       [ -1.75615960e-03,  -4.17347693e-01,  -8.19546081e-01],
       [  1.90178603e-02,   1.34837899e-01,  -1.69999003e+00]])

OP 的例子。重要提示：I 和 C 必须使用 np.arrays，np.matrix 将不起作用。

>>> I=np.array([[2,0,0],[0,6,0],[0,0,5]])
>>> C=np.array([[4,7,0],[7,8,-4],[0,-4,1]])
>>> 
>>> J = np.sqrt(np.einsum('ii->i', I))
>>> B = C / np.outer(J, J)
>>> eval_, evec = np.linalg.eigh(B)
>>> evec /= J[:, None]
>>> 
>>> evec
array([[-0.35578356, -0.31094779, -0.52605088],
       [ 0.27778714,  0.1343625 , -0.267297  ],
       [ 0.23826117, -0.37371199,  0.05975754]])
>>> eval_
array([-0.73271478,  0.48762792,  3.7784202 ])

如果I 有正数和负数条目，请使用eig 而不是eigh，然后再将平方根转换为complex dtype。

Answer 2

与其他答案不同，我假设符号 I 表示单位矩阵，Ix=x。

你要解决的，Cx=λIx，就是所谓的标准特征值问题， 并且大多数特征值求解器都解决了以该格式描述的问题，因此 Numpy 函数的签名为 eig(C)。

如果您的 C 矩阵是对称矩阵并且您的问题确实是标准特征值问题，我建议使用 numpy.linalg.eigh，它已针对此类问题进行了优化。

相反，如果您的问题确实是一个广义特征值问题，例如频率方程 Kx=ω²Mx，您可以使用 scipy.linalg.eigh，它支持对称矩阵的此类问题陈述.

eigvals, eigvecs = scipy.linalg.eigh(C, I)

关于特征值的差异，Numpy 实现不保证它们的顺序，所以它可能只是一个不同的顺序，但如果你的问题确实是一个普遍的问题（我 不是单位矩阵...) 解决方案当然不同，您必须使用 eigh 的 Scipy 实现。

如果差异在特征向量内，请记住特征向量在任意比例因子内是已知的，并且再次，排序可能是未定义的（但是，当然，它们的顺序与您拥有特征值的顺序相同) — scipy.linalg.eigh 的情况略有不同，因为在这种情况下，特征值已排序，特征向量相对于第二个矩阵参数（在您的示例中为 I）进行归一化。

---

Ps：scipy.linalg.eigh 行为（即，已排序的特征值和归一化的特征向量）对于 我的 用例来说非常方便，我也使用它来解决标准特征值问题。

Answer 3

使用SymPy:

>>> from sympy import *
>>> t = Symbol('t')
>>> D = diag(2,6,5)
>>> S = Matrix([[ 4, 7, 0],
                [ 7, 8,-4],
                [ 0,-4, 1]])
>>> (t*D - S).det()
60*t**3 - 212*t**2 - 77*t + 81

计算精确根：

>>> roots = solve(60*t**3 - 212*t**2 - 77*t + 81,t)
>>> roots
[53/45 + (-1/2 - sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3) + 14701/(8100*(-1/2 - sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3)), 53/45 + 14701/(8100*(-1/2 + sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3), 53/45 + 14701/(8100*(312469/182250 + sqrt(797521629)*I/16200)**(1/3)) + (312469/182250 + sqrt(797521629)*I/16200)**(1/3)]

计算根的浮点近似值：

>>> for r in roots:
...     r.evalf()
... 
0.487627918145732 + 0.e-22*I
-0.73271478047926 - 0.e-22*I
3.77842019566686 - 0.e-21*I

请注意，根是真实的。