使用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)

%RESULT

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)

%RESULT

(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]]))

- andreyxdd

---

你确定结果不会因特征向量的排序/方向而有所不同吗？ - Ander Biguri

我已经添加了一个例子。当然，无论是在Python还是在Matlab中使用eig(C)，输出结果都是相同的，但我需要确切地使用eig(C,I)。 - andreyxdd

我明白，但是要让你知道的是：它们不一定输出相同，因为向量/值的顺序和方向可能不同。它们在数学上是相同的，但数字可能会有所不同。 - Ander Biguri

是的，已经理解了，谢谢。 - andreyxdd

使用 I 表示对角矩阵而非单位矩阵是非常糟糕的符号表示法。 - Rodrigo de Azevedo

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的例子。重要提示：必须使用np.array来表示I和C，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 类型。保留 HTML 标记。

- Paul Panzer

---

没错，Matrix里只有正数项。谢谢，我会试一下的！ - andreyxdd

@andreyxdd 很棒。只有一件事（很重要）：不要使用np.matrix类！它不能正常工作。请改用np.array。 - Paul Panzer

Answer 2

不同于其他答案，我假设符号代表单位矩阵，x=x。

您想要解决的问题Cx=λIx是所谓的标准特征值问题，大多数特征值求解器处理该格式描述的问题，因此Numpy函数的签名为eig（C）。

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

相反，如果您的问题确实是广义特征值问题，例如频率方程Kx=ω²Mx，则可以使用支持对称矩阵的问题语句的scipy.linalg.eigh。

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

关于特征值的差异，Numpy实现不保证它们的排序，因此可能只是不同的排序方式，但如果您的问题确实是广义问题（I不是单位矩阵...），解决方案当然是不同的，您必须使用Scipy实现的eigh。

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

---

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

- gboffi

---

@andreyxdd 很高兴能够帮到你！ - gboffi

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

请注意根是实数。

- Rodrigo de Azevedo