SymPy 无法计算此矩阵的特征值答案

【问题标题】：SymPy could not compute the eigenvalues of this matrixSymPy 无法计算此矩阵的特征值
【发布时间】：2018-10-25 15:57:17
【问题描述】：

我想计算一个拉普拉斯矩阵的第二个特征值来检查对应的图是否连通，但是当我尝试使用 SymPy 的eigenvals 时，很多时候它会抛出一个错误

MatrixError: Could not compute eigenvalues for 
Matrix([[1.00000000000000, 0.0, 0.0, 0.0, -1.00000000000000, 0.0, 0.0, 0.0, 0.0, 0.0], 
        [0.0, 1.00000000000000, 0.0, 0.0, 0.0, -1.00000000000000, 0.0, 0.0, 0.0, 0.0], 
        [0.0, 0.0, 1.00000000000000, 0.0, 0.0, 0.0, 0.0, 0.0, -1.00000000000000, 0.0], 
        [0.0, 0.0, 0.0, 1.00000000000000, 0.0, 0.0, 0.0, 0.0, -1.00000000000000, 0.0], 
        [-1.00000000000000, 0.0, 0.0, 0.0, 1.00000000000000, 0.0, 0.0, 0.0, 0.0, 0.0], 
        [0.0, -1.00000000000000, 0.0, 0.0, 0.0, 3.00000000000000, 0.0, 0.0, -1.00000000000000, -1.00000000000000], 
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.00000000000000, 0.0, -1.00000000000000], 
        [0.0, 0.0, -1.00000000000000, -1.00000000000000, 0.0, -1.00000000000000, 0.0, 0.0, 3.00000000000000, 0.0], 
        [0.0, 0.0, 0.0, 0.0, 0.0, -1.00000000000000, 0.0, -1.00000000000000, 0.0, 2.00000000000000]])

环顾四周，我发现由于 SymPy 进行符号计算，浮点数可能会成为问题。所以我尝试了：

降低浮点Float(tmp[i][j], 3)的精度，但没有帮助。
我尝试将浮点数转换为 Rational list(map(nsimplify, tmp[i]))，但没有帮助。
我尝试将浮点数转换为 int list(map(int, tmp[i]))，但也没有用。

我真的不明白为什么它不起作用，即使我将每个元素都转换为int。

【问题讨论】：

你确定矩阵有真正的特征值吗？
如果你对特征值的浮点计算没问题，你应该更喜欢 numpy 或 scipy 例程。例如，scipy.linalg.eig（或者，更好的是，scipy.linalg.eigh，因为矩阵是厄米特矩阵）识别特征值没有问题
@internet_user 我确信矩阵有真正的特征值。 @Stelios 我已经看过scipy.linalg.eigh，但我不确定特征值的精度，因为我需要准确知道第二小的特征值是否大于或等于零。你认为scipy.linalg.eigh 执行的数值计算总能让我检索到这些正确信息吗？
从拉普拉斯矩阵的第 7 行/第 7 列只有零这一事实可以得出图不连通的结论。

标签： python matrix sympy symbolic-math eigenvalue

【解决方案1】：

由于拉普拉斯算子是整数矩阵，所以我们使用整数：

L = Matrix([[ 1,  0,  0,  0, -1,  0, 0,  0,  0,  0],
            [ 0,  1,  0,  0,  0, -1, 0,  0,  0,  0],
            [ 0,  0,  1,  0,  0,  0, 0,  0, -1,  0],
            [ 0,  0,  0,  1,  0,  0, 0,  0, -1,  0],
            [-1,  0,  0,  0,  1,  0, 0,  0,  0,  0],
            [ 0, -1,  0,  0,  0,  3, 0,  0, -1, -1],
            [ 0,  0,  0,  0,  0,  0, 0,  0,  0,  0],
            [ 0,  0,  0,  0,  0,  0, 0,  1,  0, -1],
            [ 0,  0, -1, -1,  0, -1, 0,  0,  3,  0],
            [ 0,  0,  0,  0,  0, -1, 0, -1,  0,  2]])

计算特征值：

>>> L.eigenvals()
{0: 3, 1: 1, 2: 1}

这很奇怪，因为矩阵是 10×10，而不是 5×5。

我尝试计算 Jordan 范式，但无法完成，因为函数 jordan_form 产生了错误消息 IndexError: list index out of range。

计算特征多项式：

>>> s = Symbol('s')
>>> p = (s * eye(10) - L).det()
>>> p
s**10 - 14*s**9 + 77*s**8 - 214*s**7 + 321*s**6 - 256*s**5 + 99*s**4 - 14*s**3

请注意，最低阶的单项式是三次的。这使我们可以得出结论，特征值 0 的重数为 3，因此，图 未连接。

让我们尝试找到特征多项式的根：

>>> solve(p,s)
[0, 0, 0, 1, 2, CRootOf(s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7, 0), CRootOf(s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7, 1), CRootOf(s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7, 2), CRootOf(s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7, 3), CRootOf(s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7, 4)]

请注意，实际上只找到了 5 个根（eigenvals 也只产生了 5 个特征值）。 5 个缺失的根是五次元 s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7 的根。

自 19 世纪以来，人们就知道并非所有 5 次（或更高）的多项式都有可以使用算术运算和根式表示的根。因此，我们可能会要求 SymPy 完成不可能。最好使用 NumPy 计算 10 个特征值的近似值。

【讨论】：

【解决方案2】：

在为它增加maxsteps 参数后，您可以使用nroots 获得特征多项式的所有10 个根的数值近似值：

>>> p = s**10 - 14*s**9 + 77*s**8 - 214*s**7 + 321*s**6 - 256*s**5 + 99*s**4 - 14*s**3
>>> [i.n(2) for i in nroots(eq,maxsteps=100)]
[0, 0, 0, 0.32, 0.68, 1.0, 2.0, 2.1, 3.2, 4.6]

solve(p, s) 中的 CRootOf 实例实际上也是解决方案，可以进行数值评估：

>>> CRootOf(s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7, 0).n(2)
0.32

【讨论】：