Search code examples
pythonmatrixsympysymbolic-matheigenvalue

SymPy could not compute the eigenvalues of this matrix


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:

  1. To reduce the precision of the floating point Float(tmp[i][j], 3), but it didn't help.
  2. I have tried to convert floats to Rational list(map(nsimplify, tmp[i])), but it didn't help.
  3. 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.


Solution

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