SymPy could not compute the eigenvalues of this matrix

Question

I want to compute the second eigenvalue of a Laplacian matrix to check if the corresponding graph is connected or not, but when I try to use SymPy's eigenvals, a lot of times it happens that it throws an error

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

Looking around I found out that since SymPy does symbolic computation, floating points can be a problem for it. So I tried:

To reduce the precision of the floating point Float(tmp[i][j], 3), but it didn't help.
I have tried to convert floats to Rational list(map(nsimplify, tmp[i])), but it didn't help.
I have tried to convert floats to int list(map(int, tmp[i])), but it didn't help neither.

I really can't understand why it doesn't work out, even if I convert every element to int.

If you are fine with a floating point evaluation of eigenvalues, you should prefer either numpy or scipy routines for that purpose. For example, `scipy.linalg.eig` (or, even better, `scipy.linalg.eigh` as the matrix is hermitian) has no problem identifies the eigenvalues — Stelios, May 15 '18 at 18:18
@internet_user I am sure the matrix has real eigenvalues. @Stelios I have already seen `scipy.linalg.eigh`, but I wasn't sure about the precision of the eigenvalue, since I need to exactly know if the second smallest eigenvalue is greater or equal than zero. Do you think the numerical computation performed by `scipy.linalg.eigh` will always allow me to retrieve this correct information? — BalrogOfMoria, May 15 '18 at 20:21
You can conclude that the graph is not connected from the fact that the 7th row/column of the Laplacian matrix has only zeros. — Rodrigo de Azevedo, May 18 '18 at 13:23

Rodrigo de Azevedo · Accepted Answer · 2018-05-16T15:31:41.297

Since the Laplacian is an integer matrix, let us use integers:

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

Computing the eigenvalues:

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

This is very strange, as the matrix is 10-by-10, not 5-by-5.

I tried to compute the Jordan normal form, but could not do it, as function jordan_form produced the error message IndexError: list index out of range.

Computing the characteristic polynomial:

>>> 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

Note that the monomial of lowest degree is cubic. This allows us to conclude that the multiplicity of eigenvalue 0 is 3 and, thus, the graph is not connected.

Let us try to find the roots of the characteristic polynomial:

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

Note that only 5 roots were actually found (eigenvals produced only 5 eigenvalues, too). The 5 missing roots are the roots of the quintic s**5 - 11*s**4 + 42*s**3 - 66*s**2 + 39*s - 7.

It has been known since the 19th century that not all polynomials of degree 5 (or higher) have roots that can be expressed using arithmetic operations and radicals. Hence, we may be asking SymPy to do the impossible. Better use NumPy to compute approximate values of the 10 eigenvalues.

score 1 · Answer 2 · answered Mar 26 '21 at 14:46

You can get numerical approximations of all 10 roots of the characteristic polynomial using nroots after bumping up the maxsteps parameter for it:

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

The CRootOf instances from solve(p, s) are actually solutions, too, that can be numerically evaluated:

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

SymPy could not compute the eigenvalues of this matrix

2 Answers2

Linked