Determining left-eigenvectors and eigenvalues of 100 x 100 matrix

Am trying to determine the (real) left-eigenvectors and eigenvalues of a 100 by 100 matrix. Elements of the matrix, T, are probabilities from 0 to 1. Each row of the matrix, T, sums to 1.0 (100% probability). Columns do not have to sum to anything in particular.

The scipy and numpy libraries return complex numbers - which is apparently due to not being able to solve the characteristic polynomial as real numbers or the algorithms being optimised to do this. I turned to the sympy library, which also returned complex numbers as the solution.

Using sympy as described here makes sense:
https://stackoverflow.com/questions/65170764/is-it-possible-somehow-to-find-complex-eigenvalues-using-sympy

My understanding is that this will force answers as real numbers.

Here is my code:

import numpy as np
from sympy import *

T = np.loadtxt('rep10_T_ij.dat', delimiter=' ')

T = np.asmatrix(T)
T = Matrix(T)

lam = symbols('\lambda')
char_poly = det(lam * eye(T.rows) - T)
roots_char_poly = solve(char_poly, lam)
display(char_poly)
display(roots(poly(char_poly, domain=CC)))
display([N(r) for r in roots_char_poly]) 

However this does not output an answer in 12 hours with any amount of hardware thrown at it.

I read this here:
https://github.com/sympy/sympy/issues/7200

Is a 100 by 100 matrix too expensive?
I can throw more hardware at it on a cluster, but have not been able to try this yet. Is there a solution?

The advantage of using SymPy rather than NumPy or SciPy in this context is just if you want to perform the calculation exactly or symbolically. Determining that the eigenvalues are precisely real will require exact arithmetic in general. Exact or symbolic calculations are a lot slower than fixed precision floating point though so you should expect things to slow down a lot compared to using NumPy or SciPy. If your initial matrix is a matrix of floats though then it is already approximate in which case there is a good chance that there will be no benefit in using SymPy.

In any case here is some code that can show how to get the exact real eigenvalues of a matrix more efficiently using SymPy. This is usually only worth doing if your original matrix has exact rational numbers (not floats) but it should obtain precisely the real eigenvalues:

from sympy import *

def make_matrix(N):
    """Random matrix of rational numbers with all integer eigenvalues"""
    V = randMatrix(N)._rep.to_field()
    eigs = randMatrix(1, N)
    D = diag(*eigs)._rep
    Vinv = V.inv()
    M = V*D*Vinv
    return M.to_Matrix()

def eigenvals(M):
    x = symbols('x')
    p = M._rep.charpoly()
    return Poly(p, x).real_roots()

This uses the Berkowicz algorithm to compute the characteristic polynomial which is the slowest part of the operation (the charpoly method). This is much more efficient than using determinants to find the characteristic polynomial but still slow for large matrices. Finding the characteristic polynomial this way is roughly O(N**4) for an N x N dense matrix if the underlying coefficient field has O(1) arithmetic cost. Unfortunately exact rational numbers do not have O(1) arithmetic cost so in fact the complexity is significantly worse than O(N**4) as you scale up to larger matrices.

You can get a sense for the timings involved for different sizes:

In [80]: M = make_matrix(5)

In [81]: M
Out[81]: 
⎡-23268248970   17481165584  6611878976  21688404260  -332363048 ⎤
⎢─────────────  ───────────  ──────────  ───────────  ───────────⎥
⎢   6173287       6173287     6173287      6173287      6173287  ⎥
⎢                                                                ⎥
⎢-11620321902   8914715596   3291699676  10654796030  -168333588 ⎥
⎢─────────────  ──────────   ──────────  ───────────  ───────────⎥
⎢   6173287      6173287      6173287      6173287      6173287  ⎥
⎢                                                                ⎥
⎢-17978224844   13433529668  5243135950  16686586148  -354021332 ⎥
⎢─────────────  ───────────  ──────────  ───────────  ───────────⎥
⎢   6173287       6173287     6173287      6173287      6173287  ⎥
⎢                                                                ⎥
⎢-10824267464   8003053808   3040726208  10253515930  -120066032 ⎥
⎢─────────────  ──────────   ──────────  ───────────  ───────────⎥
⎢   6173287      6173287      6173287      6173287      6173287  ⎥
⎢                                                                ⎥
⎢-2998114302    2232013094   828083788   2727588566    239697782 ⎥
⎢────────────   ──────────   ─────────   ──────────    ───────── ⎥
⎣  6173287       6173287      6173287     6173287       6173287  ⎦

In [82]: %time eigenvals(M)
CPU times: user 43.8 ms, sys: 4 ms, total: 47.8 ms
Wall time: 46.6 ms
Out[82]: [14, 26, 46, 56, 82]

In [83]: M = make_matrix(10)

In [84]: %time eigenvals(M)
CPU times: user 82.7 ms, sys: 10 µs, total: 82.7 ms
Wall time: 81.8 ms
Out[84]: [8, 18, 27, 39, 58, 67, 85, 95, 96, 98]

In [85]: M = make_matrix(20)

In [86]: %time eigenvals(M)
CPU times: user 652 ms, sys: 2 µs, total: 652 ms
Wall time: 651 ms
Out[86]: [9, 15, 21, 22, 26, 26, 36, 37, 38, 40, 52, 57, 63, 64, 65, 66, 69, 70, 80, 96]

In [87]: M = make_matrix(30)

In [88]: %time eigenvals(M)
CPU times: user 6.1 s, sys: 0 ns, total: 6.1 s
Wall time: 6.1 s
Out[88]: 
[7, 8, 10, 13, 15, 18, 26, 33, 34, 35, 36, 42, 42, 43, 43, 47, 49, 50, 52, 54, 56, 62, 72, 79, 85, 
87, 89, 92, 98, 98]

In [89]: M = make_matrix(40)

In [90]: %time eigenvals(M)
CPU times: user 37.8 s, sys: 12 µs, total: 37.8 s
Wall time: 37.9 s
Out[90]: 
[1, 7, 13, 14, 15, 17, 17, 20, 26, 26, 27, 28, 29, 32, 33, 35, 37, 43, 43, 44, 44, 46, 50, 50, 52, 
52, 53, 55, 58, 64, 67, 70, 71, 72, 79, 83, 83, 88, 88, 99]

So already finding the exact integer eigenvalues of a 40x40 matrix takes 30 seconds. Scaling up to 100x100 will be slow but might be possible. I recommend having gmpy2 installed if you want to do this because it makes exact rational numbers in SymPy a lot faster.

Most likely though computing exact eigenvalues like this is not what you want to do and you should instead use NumPy/SciPy. If they give eigenvalues with small imaginary parts then you can discard those if you know that the eigenvalues should be real. Otherwise as long as the initial matrix is only approximate (i.e. floats) then the kind of slow exact computing shown above is probably not worthwhile.