5

Until now I used numpy.linalg.eigvals to calculate the eigenvalues of quadratic matrices with at least 1000 rows/columns and, for most cases, about a fifth of its entries non-zero (I don't know if that should be considered a sparse matrix). I found another topic indicating that scipy can possibly do a better job.

However, since I have to calculate the eigenvalues for hundreds of thousands of large matrices of increasing size (possibly up to 20000 rows/columns and yes, I need ALL of their eigenvalues), this will always take awfully long. If I can speed things up, even just the tiniest bit, it would most likely be worth the effort.

So my question is: Is there a faster way to calculate the eigenvalues when not restricting myself to python?

Community
  • 1
  • 1
alice
  • 254
  • 2
  • 4
  • 9
  • If python is not a must, then any other lower level language (C++ / even C#) will give you a speed boost. Only matter of a suitable implementation. – ZorleQ Jul 04 '13 at 10:54
  • 1
    Whatever you do, bear in mind that a lot of `numpy` is a Python-friendly wrapper around functionality written in languages such as C, Fortran, assembler. I see from the docs that `numpy.linalg.eigvals` i a wrapper around functions in the LINPACK library. This doesn't mean that you can't find faster solvers, but you may have to look beyond numpy, scipy and LAPACK to find them. – High Performance Mark Jul 04 '13 at 11:09
  • Do you use iterative methods ? If so maybe you can parallelise them ? – Darek Jul 25 '13 at 14:29

2 Answers2

5

@HighPerformanceMark is correct in the comments, in that the algorithms behind numpy (LAPACK and the like) are some of the best, but perhaps not state of the art, numerical algorithms out there for diagonalizing full matrices. However, you can substantially speed things up if you have:

Sparse matrices

If your matrix is sparse, i.e. the number of filled entries is k, is such that k<<N**2 then you should look at scipy.sparse.

Banded matrices

There are numerous algorithms for working with matrices of a specific banded structure. Check out the solvers in scipy.linalg.solve.banded.

Largest Eigenvalues

Most of the time, you don't really need all of the eigenvalues. In fact, most of the physical information comes from the largest eigenvalues and the rest are simply high frequency oscillations that are only transient. In that case you should look into eigenvalue solutions that quickly converge to those largest eigenvalues/vectors such as the Lanczos algorithm.

Slater Victoroff
  • 21,376
  • 21
  • 85
  • 144
Hooked
  • 84,485
  • 43
  • 192
  • 261
  • OP explicitly says that they need all of the eigenvalues, and that the matrices are about 80% sparse. I don't know if 80% is sufficient sparsity for sparse eigenvalue algorithms to outperform dense ones, but it's worth trying. – Danica Jul 25 '13 at 14:37
  • 2
    @Dougal I know he thinks he needs all of the eigenvalues, but I've learned that often a great approximation can be made with only the largest eigenvalues (for obvious reasons!). Lancozs algorithms will eventually converge onto the smaller and smaller eigenvalues, and this information is **certainly** better than having no eigenvalues at all! – Hooked Jul 25 '13 at 14:39
  • Yeah, sorry, I read that section of your answer too quickly and assumed you hadn't really read the question - sorry. :) But a protip: the username "alice" and the pronoun "he" typically don't go together. Probably better to stick to gender-neutral [singular they](http://en.wikipedia.org/wiki/Singular_they). – Danica Jul 25 '13 at 14:41
  • 1
    @Dougal who am I to enforce gender onto a username? I think however, that OP is the proper pronoun for all Stack Exchange quieres. Perhaps a good question for http://english.stackexchange.com? – Hooked Jul 25 '13 at 14:44
1

An easy way to maybe get a decent speedup with no code changes (especially on a many-core machine) is to link numpy to a faster linear algebra library, like MKL, ACML, or OpenBLAS. If you're associated with an academic institution, the excellent Anaconda python distribution will let you easily link to MKL for free; otherwise, you can shell out $30 (in which case you should try the 30-day trial of the optimizations first) or do it yourself (a mildly annoying process but definitely doable).

I'd definitely try a sparse eigenvalue solver as well, though.

Danica
  • 28,423
  • 6
  • 90
  • 122