[Python] Fastest way to find eigenvalues and Eigenvectors for VERY large matrices

I am trying to find the eigenvectors and eigenvalues for very large square symmetric matrices. They are Laplacian matrices on the vicsek fractal. The matrices are in the scale of (5n ⋅ 4 + 1) where i am looking at n=6, 7, 8 ..., so in the scale
of 62 501, 312 501, 1 562 501 ... .

So far i have tried:


Since they are symmetric i use numpy.linalg.eigh for which we get the runtime of

|n value|# of col. & rows|Time in Seconds|
|-------|----------------|---------------|
| 4 | 2 501 | ~14 |
| 5 | 12 501 | ~442 |


I have also tried numpy for n = 7, but ended up running out of memory.

Next i tried scipy.scarse.linalg.eigsh since the laplacian matrices are sparse (there are at max 5 entries pr. column that are non-zero) but this took ~72 seconds to run for n = 4 and took well over 600 seconds for n = 5 (i ended up stopping the program for taking too long).

I have also tried the package mpmath but that ended up being much slower.

I don't think it is my own code that is slowing this down, since i only call the functions from either numpy or scipy.

def laplacianOperatorMatrix(allPointsDict: dict):
return [len(allPointsDict["neighbours"]) if i == j else -1 if j in allPointsDict["neighbours"] else 0 for j in allPointsDict]

NpEigvec, NpEigMatrix = np.linalg.eigh(laplacianOperatorMatrix(allPointsDict="some dict"))

SciEigVec, SciEigMatrix = linalg.eigsh(sparse.csc_matrix(laplacianOperatorMatrix(allPointsDict="some dict")), k=pow(5, n)*4 + 1 - 1)


Is there way to optimize the finding of the eigenvectors? Or maybe some packages/method for finding eigenvectors that i don't know of?

Continue reading...