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.

Answer 1

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

Answer 2

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 ?

---

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

---

AA

Hope this helps.

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]