I need to find the smallest few eigenvalues of a Hamiltonian (exact diagonalization) I use Python, and SciPy's built-in sparse eigenvalue solver. I notice, however, that for my small system (only a 40k x 40k matrix) the program takes hours, maybe even days. I began to suspect something was up.
One thing I note is that my matrix should have a zero eigenvalue. Will that cause havoc in the eigenvalue solvers? If so, what is the way to remedy this?
The SciPy tutorial explicitly states
Note that ARPACK is generally better at finding extremal eigenvalues: that is, eigenvalues with large magnitudes. In particular, using
which = 'SM'may lead to slow execution time and/or anomalous results. A better approach is to use shift-invert mode.
and goes on to describe that. Basically, if $\lambda$ is the smallest by magnitude eigenvalue of $A$, then $\nu = \frac{1}{\lambda}$ is the largest by magnitude eigenvalue of $A^{-1}$. (With the same trick, you can get the eigenvalue closest to a given $\sigma$ by looking at the largest for $(A-\sigma I)^{-1}$.) Hence, to compute, say, the three smallest by magnitude eigenvalues, you can instead use
eigenval, eigenvec = eigsh(A, 3, sigma=0, which='LM')
If you already know that $A$ has a zero eigenvalue and is therefore not invertible (and you have some idea where the next closest are), you can instead use a nonzero shift, e.g., sigma=0.01 to shift your matrix away from singular. If sigma is small enough, that should still get you the ones closest to zero.
(Also, as noted in my comment, make sure you have a proper optimized BLAS such as OpenBLAS installed -- for sparse eigenvalues, SciPy calls ARPACK (written in Fortran) which in turn makes heavy use of BLAS.)
EDIT: If it's the smallest algebraic eigenvalue you're after, you can use this neat trick: If you have an estimate $\sigma$ of the largest eigenvalue (say, by computing it using which='LM'), you can use that to shift all eigenvalues to be negative -- so that the smallest algebraic eigenvalues are also the largest magnitude ones. Then simply use which='LM' on $A-\sigma I$ and shift back.
While the smallest eigenvalues of a sparse 40k x 40k matrix might be obtained in a reasonable amount of time using ARPACK and some handy tricks through the SciPy wrapper, as mentioned here, the problem will scale at best as $N^2$.
For any larger eigenvalue problems I almost always recommend SLEPc, which also has a Python wrapper (through perhaps not as natural as the C versions) and in many cases calls parallel versions of the routines SciPy would.
