2
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ 40k is huge for trying to find all eigenvalues. Sparse eigenvalue solvers are built to get a few selected eigenvalues effectively, not all of them. A zero eigenvalue doesn't affect this, though. (Side note: Make sure you have a decent optimized BLAS installed. The difference between scipy linked against a stock BLAS and, say, OpenBLAS is staggering.) $\endgroup$ – Christian Clason Jun 30 '16 at 21:49
  • $\begingroup$ Want to write the equations for your system? What is the potential? $\endgroup$ – nicoguaro Jun 30 '16 at 21:53
  • $\begingroup$ @ChristianClason I do only seek out the smallest eigenvalues. Sorry if that wasn't explicit enough. I fixed it. $\endgroup$ – Aurey Jun 30 '16 at 21:58
  • $\begingroup$ @nicoguaro It's not a Schrodinger-like equation where $H=T+V$, rather, it's in second quantized form and is not quadratic, meaning the size of the Hilbert space scales as $2^N$. Not sure if it'll be too helpful. $\endgroup$ – Aurey Jul 1 '16 at 0:13
  • $\begingroup$ Dirac's equation-like? There should be a form you write your problem down $\endgroup$ – nicoguaro Jul 1 '16 at 0:20
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Nice, I agree. I am more interested in the smallest algebraic values, though (which = 'SA') since it's a quantum problem, and we almost always are most interested in the lowest-lying states. $\endgroup$ – Aurey Jul 1 '16 at 0:08
  • $\begingroup$ Also, I took your advice and linked up numpy and scipy with OpenBLAS. You're right, I'm actually seeing the program finish! $\endgroup$ – Aurey Jul 1 '16 at 0:09
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.