Finding smallest eigenvectors of large sparse matrix, over 100x slower in SciPy than in Octave

Question

I am trying to compute few (5-500) eigenvectors corresponding to the smallest eigenvalues of large symmetric square sparse-matrices (up to 30000x30000) with less than 0.1% of the values being non-zero.

I am currently using scipy.sparse.linalg.eigsh in shift-invert mode (sigma=0.0), which I figured out through various posts on the topic is the prefered solution. However, it takes up to 1h to solve the problem in most cases. On the other hand the function is very fast, if I ask for the largest eigenvalues (sub seconds on my system), which was expected from the documentation.

Since I am more familiar with Matlab from work, I tried solving the problem in Octave, which gave me the same result using eigs (sigma=0) in just a few seconds (sub 10s). Since I want to do a parameter sweep of the algorithm including the eigenvector computation, that kind of time gain would be great to have in python as well.

I first changed the parameters (especially tolerance), but that didn't change much on the timescales.

I am using Anaconda on Windows, but tried to switch the LAPACK/BLAS used by scipy (which was a huge pain) from mkl (default Anaconda) to OpenBlas (used by Octave according to the documentation), but could not see a change in performance.

I was not able to figure out, whether there was something to change about the used ARPACK (and how)?

I uploaded a testcase for the code below to the following dropbox-folder: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

In Python

import numpy as np
from scipy.sparse import csr_matrix, csc_matrix, linalg, load_npz   
M = load_npz('M.npz')
evals, evecs = linalg.eigsh(M,k=6,sigma=0.0)

In Octave:

M=dlmread('M.txt');
M=spconvert(M);
[evecs,evals] = eigs(M,6,0);

Any help is appriciated!

Some additional options I tried out based on the comments and suggestions:

Octave: eigs(M,6,0) and eigs(M,6,'sm') give me the same result:

[1.8725e-05 1.0189e-05 7.5622e-06 7.5420e-07 -1.2239e-18 -2.5674e-16]

while eigs(M,6,'sa',struct('tol',2)) converges to

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

much faster, but only if the tolerance values is above 2, otherwise it does not converge at all and the values are strongly different.

Python: eigsh(M,k=6,which='SA') and eigsh(M,k=6,which='SM') both do not converge (ARPACK error on no convergence reached). Only eigsh(M,k=6,sigma=0.0) gives some eigenvalues (after almost an hour), which are different to octave for the smallest ones (even 1 additional small value is found):

[3.82923317e-17 3.32269886e-16 2.78039665e-10 7.54202273e-07 7.56251500e-06 1.01893934e-05]

If the tolerance is high enough I also get results from eigsh(M,k=6,which='SA',tol='1'), which come close the other obtained values

[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]

again with a different number of small eigenvalues. The computation time is still almost 30min. While the different very small values might be understandable, since they might represent multiples of 0, the different multiplicity baffles me.

Additionally, there seems to be some fundamental differences in SciPy and Octave, which I cannot figure out, yet.

Babe

1 - I assume you meant to put brackets around [evals,evecs] in the octave code? 2 - can you include a small example for M? or maybe a generator script for one if that's possible? – Lutetium

1 - Yes exactly, I edited my post. 2 - I tried out the performance for some sub-matrices of my data and it seems that Octave is always faster, but for smaller matrices below 5000x5000 its just a factor of 2-5 times, above that it gets really ugly. And since its "real data" I can not give a generator script. Is there a standard way to upload an example somehow? Due to the sparsity, a npz-file is reasonably small. – Babe

I guess you can share a link to any cloud storage facility. – Chacha

Thx. I included a dropbox link in the original post and updated the code to a working example. – Babe

Juts to reinforce your point, I tested on Matlab R2019b and got 84 seconds vs 36.5 mins in Python 3.7, Scipy 1.2.1 (26 times faster). – Marthmartha

@Bill, are the scipy-arpack eigenvalues == -2.57e-16 ... in the post ? Did you run which="LM", sigma=0, k=6, tol=0, v0=None (v0 random) ? – Academic

Were both runs with tol=0 ? Are the eigenvalues and eigenvectors np.allclose ? If you could post them, I'd like to take a look. Thanks ! – Academic

Answer 1

Added 19 May: Cholesky inner solver:

The doc for scipy eigsh says

shift-invert mode ... requires an operator to compute the solution of the linear system (A - sigma * I) x = b ... This is computed internally via a sparse LU decomposition (splu) for an explicit matrix, or via an iterative solver for a general linear operator.

ARPACK calls this "inner solver" many many times, depending on tol etc.; obviously, slow inner solver => slow eigs. For A posdef, sksparse.cholmod is waaay faster than splu.

Matlab eig also uses cholesky:

If A is Hermitian and B is Hermitian positive definite, then the default for algorithm is 'chol'

---

Fwiw, np.linalg.eigh finds all eigenvalues and eigenvectors of the 7 Gb dense matrix A.A in under an hour on my old 4-core imac -- wow. Its spectrum looks like this:

---

February 2020, TL;DR

A conjecture and some comments, since I don't have Matlab / Octave:

To find small eigenvalues of symmetric matrices with eigenvalues >= 0, like yours, the following is waaay faster than shift-invert:

# flip eigenvalues e.g.
# A:     0 0 0 ... 200 463
# Aflip: 0 163 ... 463 463 463
maxeval = eigsh( A, k=1 )[0]  # biggest, fast
Aflip = maxeval * sparse.eye(n) - A
bigevals, evecs = eigsh( Aflip, which="LM", sigma=None ... )  # biggest, near 463
evals = maxeval - bigevals  # flip back, near 463 -> near 0
# evecs are the same

eigsh( Aflip ) for big eigenpairs is faster than shift-invert for small, because A * x is faster than the solve() that shift-invert must do. Matlab / Octave could conceivably do this Aflip automatically, after a quick test for positive-definite with Cholesky.
Can you run eigsh( Aflip ) in Matlab / Octave ?

Other factors that may affect accuracy / speed:

Arpack's default for the start vector v0 is a random vector. I use v0 = np.ones(n), which may be terrible for some A but is reproducible :)

This A matrix is almost exactly sigular, A * ones ~ 0.

Multicore: scipy-arpack with openblas / Lapack uses ~ 3.9 of the 4 cores on my iMac; do Matlab / Octave use all cores ?

---

Here are the scipy-Arpack eigenvalues for several k and tol, grepped from logfiles under gist.github:

k 10  tol 1e-05:    8 sec  eigvals [0 8.5e-05 0.00043 0.0014 0.0026 0.0047 0.0071 0.0097 0.013 0.018] 
k 10  tol 1e-06:   44 sec  eigvals [0 3.4e-06 2.8e-05 8.1e-05 0.00015 0.00025 0.00044 0.00058 0.00079 0.0011] 
k 10  tol 1e-07:  348 sec  eigvals [0 3e-10 7.5e-07 7.6e-06 1.2e-05 1.9e-05 2.1e-05 4.2e-05 5.7e-05 6.4e-05] 

k 20  tol 1e-05:   18 sec  eigvals [0 5.1e-06 4.5e-05 0.00014 0.00023 0.00042 0.00056 0.00079 0.0011 0.0015 0.0017 0.0021 0.0026 0.003 0.0037 0.0042 0.0047 0.0054 0.006
k 20  tol 1e-06:   73 sec  eigvals [0 5.5e-07 7.4e-06 2e-05 3.5e-05 5.1e-05 6.8e-05 0.00011 0.00014 0.00016 0.0002 0.00025 0.00027 0.0004 0.00045 0.00051 0.00057 0.00066
k 20  tol 1e-07:  267 sec  eigvals [-4.8e-11 0 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 5.8e-05 6.4e-05 6.9e-05 8.3e-05 0.00011 0.00012 0.00013 0.00015

k 50  tol 1e-05:   82 sec  eigvals [-4e-13 9.7e-07 1e-05 2.8e-05 5.9e-05 0.00011 0.00015 0.00019 0.00026 0.00039 ... 0.0079 0.0083 0.0087 0.0092 0.0096 0.01 0.011 0.011 0.012
k 50  tol 1e-06:  432 sec  eigvals [-1.4e-11 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00081 0.00087 0.00089 0.00096 0.001 0.001 0.0011 0.0011
k 50  tol 1e-07: 3711 sec  eigvals [-5.2e-10 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00058 0.0006 0.00063 0.00066 0.00069 0.00071 0.00075

versions: numpy 1.18.1  scipy 1.4.1  umfpack 0.3.2  python 3.7.6  mac 10.10.5 

Are Matlab / Octave about the same ? If not, all bets are off -- first check correctness, then speed.

Why do the eigenvalues wobble so much ? Tiny < 0 for a supposedly non-negative-definite matrix are a sign of roundoff error, but the usual trick of a tiny shift, A += n * eps * sparse.eye(n), doesn't help.

---

Where does this A come from, what problem area ? Can you generate similar A, smaller or sparser ?

Hope this helps.

Academic

Thanks for your input and sorry for the (very) late reply. The project I used this for is already complete, but I'm still curious, so I checked. Sadly, the eigenvalues in Ocatve turn out different, for k = 10 I find [-2.5673e-16 -1.2239e-18 7.5420e-07 7.5622e-06 1.0189e-05 1.8725e-05 2.0265e-05 2.1568e-05 4.2458e-05 5.1030e-05] which is also independent of the tolerance value in the range 1e-5 to 1e-7. So there is another difference here. Don't you think it is strange that scipy (including your suggestion) yield different small values dependent on the number of values queried? – Babe

Sorry, than I don't get it. You specifiy scipy 1.4.1 in your original answer and still the number of very small eigenvalues changes with k. And the eigenvalues I posted in my answer above are in Octave. They are obviously different from the ones you got with scipy. So the problem persists, or am I missing something here? – Babe

Answer 2

I know this is old now, but I had the same problem. Did you review here (https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html)?

It seems like when you set sigma to a low number (0) you should set which='LM', even though you are wanting to low values. This is because setting sigma transforms the values you want (low in this case) to appear to be high and so you still are able to take advantage of the 'LM' methods, which are much faster to get what you want (the low eigenvalues).

Conakry

Did this actually change performance for you? That would be a suprise for me. I knew the link you postet and I also implicitly specified which='LM' in my example. Because the default value for an unset which is 'LM'. I still checked though, and the performance is unchanged for my example. – Babe

Indeed, seems to have a similar difference as from Python to octave for you. I also had a large matrix that I was trying to decompose and ended up using eigsh(matrix, k=7, which='LM', sigma=1e-10). Originally I was incorrectly specifying which='SM' thinking I needed to do that to get the smallest eigenvalues and it was taking forever to give me an answer. Then, I found that article and realized that you just needed to specify it to the faster 'LM', and set k to be whatever you wanted and it would speed things up. Is your matrix actually hermitian? – Conakry

Okay, but thats a different matter. In fact if you specify sigma close to zero, than you trigger a transformation where the LM of the transformed matrix are the SM of the original one. So as I said there is no performance difference between default value of which and which='LM'. BUT a speficied sigma changes the meaning of which. Thus, which='SM' with a small sigma is simply a different query, explaining your performance difference but you would also get wrong values out of this. – Babe

Answer 3

I want to say first that I have no idea why the results you and @Bill reported are the way they are. I simply wonder if eigs(M,6,0) in Octave corresponds to k=6 & sigma=0, or perhaps it is something else?

With scipy, if I do not provide sigma, I can get a result in a decent time this way.

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigsh
from time import perf_counter
M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, k=50, which='SA', tol=1e1)
print(perf_counter() - t)
print(b)

I am entirely not sure if this makes sense though.

0.4332823531003669
[4.99011753e-03 3.32467891e-02 8.81752215e-02 1.70463893e-01
 2.80811313e-01 4.14752072e-01 5.71103821e-01 7.53593653e-01
 9.79938915e-01 1.14003837e+00 1.40442848e+00 1.66899183e+00
 1.96461415e+00 2.29252666e+00 2.63050114e+00 2.98443218e+00
 3.38439528e+00 3.81181747e+00 4.26309942e+00 4.69832271e+00
 5.22864462e+00 5.74498014e+00 6.22743988e+00 6.83904055e+00
 7.42379697e+00 7.97206446e+00 8.62281827e+00 9.26615266e+00
 9.85483434e+00 1.05915030e+01 1.11986296e+01 1.18934953e+01
 1.26811461e+01 1.33727614e+01 1.41794599e+01 1.47585155e+01
 1.55702295e+01 1.63066947e+01 1.71564622e+01 1.78260727e+01
 1.85693454e+01 1.95125277e+01 2.01847294e+01 2.09302671e+01
 2.18860389e+01 2.25424795e+01 2.32907153e+01 2.37425085e+01
 2.50784800e+01 2.55119112e+01]

The only way I found to use sigma and to get a result in a decent time is to provide M as a LinearOperator. I am not too familiar with this thing, but from what I understood my implementation represents an identity matrix, which is what M should be if not specified in the call. The reason for this is that instead of performing a direct solve (LU decomposition), scipy will use an iterative solver, which is potentially better suited. As a comparison, if you provide M = np.identity(a.shape[0]), which should be the exact same, then eigsh takes forever to yield a result. Note that this approach does not work if sigma=0 is provided. But I am not sure if sigma=0 is really that useful?

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs, eigsh, LinearOperator
from time import perf_counter


def mv(v):
    return v


M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, M=LinearOperator(shape=a.shape, matvec=mv, dtype=np.float64),
             sigma=5, k=50, which='SA', tol=1e1, mode='cayley')
print(perf_counter() - t)
print(np.sort(-5 * (1 + b) / (1 - b)))

Again, no idea if it is correct but definitely different from before. That would be great to have the input of someone from scipy.

1.4079377939924598
[3.34420263 3.47938816 3.53019328 3.57981026 3.60457277 3.63996294
 3.66791416 3.68391585 3.69223712 3.7082205  3.7496456  3.76170023
 3.76923989 3.80811939 3.81337342 3.82848729 3.84137264 3.85648208
 3.88110869 3.91286153 3.9271108  3.94444577 3.97580798 3.98868207
 4.01677424 4.04341426 4.05915855 4.08910692 4.12238969 4.15283192
 4.16871081 4.1990492  4.21792125 4.24509036 4.26892806 4.29603036
 4.32282475 4.35839271 4.37934257 4.40343219 4.42782208 4.4477206
 4.47635849 4.51594603 4.54294049 4.56689989 4.58804775 4.59919363
 4.63700551 4.66638214]

Chacha

Thank you for your input and feedback. I tried out some things to give a decent answer to your points. 1. My task at hand requires finding the k smallest eigenvalues/vectors. Therefor the approach using sigma=0 is even given in the SciPy docs: docs.scipy.org/doc/scipy/reference/tutorial/arpack.html 2. I tried out some more options, which I edited into the original question. 3. As I understand the documentaries of Octave and SciPy, eigs(M,6,0) and k=6, simga=0 should be the same. – Babe

4. Since my matrix is real and square, I thought there should be no different between SA and SM as an option, but there obviously is, at least in the computation. Am I on a wrong track here? Overall that means more questions and but no real answers or solutions from me. – Babe