What is the fastest way to compute all eigenvalues of a very big and sparse adjacency matrix in python?

Question

I'm trying to figure out if there is a faster way to compute all the eigenvalues and eigenvectors of a very big and sparse adjacency matrix than using scipy.sparse.linalg.eigsh As far as I know, this methods only uses the sparseness and symmetry attributes of the matrix. An adjacency matrix is also binary, what makes me think there is a faster way to do it.

I created a random 1000x1000 sparse adjacency matrix, and compared between several methods on my x230 ubuntu 13.04 laptop:

• scipy.sparse.linalg.eigs: 0.65 seconds
• scipy.sparse.linalg.eigsh: 0.44 seconds
• scipy.linalg.eig: 6.09 seconds
• scipy.linalg.eigh: 1.60 seconds

With the sparse eigs and eigsh, I set k, the number of the desired eigenvalues and eigenvectors, to be the rank of the matrix.

The problem starts with bigger matrices - on a 9000x9000 matrix, it took scipy.sparse.linalg.eigsh 45 minutes!

Answer 1

FILTLAN is a C++ library for computing interior eigenvalues of sparse symmetric matrices. The fact that there is a whole package devoted to just this should tell you that it's a pretty hard problem. Finding the largest or smallest few eigenvalues of a symmetric matrix can be done by shifting/inverting and using the Lanczos algorithm, but the middle of the spectrum is another matter. If you do want to use this, you can use SWIG to call a C++ program from python.

If your end goal is to compute large powers of the matrix, you could just compute eigenvectors corresponding to the largest eigenvalues, content in the knowledge that the smaller modes will be less important as you take large powers.

That said, you may indeed be better off computing the powers directly. They'll grow less and less sparse as you compute higher powers, which will mean taking up more memory; depending on how high $k$ is, you may eventually want to switch to a dense matrix.

Forgive me if these are already obvious to you: you can exploit the binary nature of the matrix by telling numpy that it consists of integers instead of floats, say by using

a = np.zeros(100,dtype=np.uint)

This will (hopefully) save some space. You can save time (but not memory) by blocking the matrix multiplications. Say you want to compute $A^{16}$; you compute $A^2$, then square this to get $A^4$, square this to get $A^8$, and so on. That way, you do $\approx\log_2 k$ matrix multiplications instead of $k$ multiplications.

You can also explore calling a parallel sparse linear algebra library like CUSP or cuSPARSE from Python if speed is your concern and you have an NVIDIA GPU.

Answer 2

I would like to comment on Daniel Shapero's answer but I don't have enough SE reputation.

The accepted answer confuses me a lot. I think shift-invert mode can be readily used to compute interior eigenvalues. See: https://docs.scipy.org/doc/scipy/tutorial/arpack.html

To answer the original question: it's rarely the case that you want all eigenvalues of a large sparse matrix. Usually, you want extrems or some cluster of interior values. In that case, for a Hermitian matrix eigsh is faster. For non-Hermitian, you will have to go with eigs. And they are much faster than numpy eig or eigh.

Answer 3

Also as an update to Daniel Shaperos answer:

If you have a GPU (doesn't matter if NVidia or AMD) available, you can try cupy. Many numpy and scipy functions are available and can be used by replacing numpy/scipy with cupy without changing the rest of the code. Arrays can be directly initialized on the GPU reducing the overhead. After your critical calculations, the cupy array can be loaded into your RAM as a numpy array.