Compute all eigenvectors and eigenvalues of small symmetric matrices

Question

My problem is to compute eigenvectors and eigenvalues of a lot of small (n < 30) symetric, positive definite matrices.

So far I am using LAPACK's DSYEV.

The priority is speed more than accuracy. Is there a faster algorithm (maybe not asymptotically, but in practice for small input)?

Answer 1

In my experience the answer to this question is not clear-cut; it is dependent on the span of the eigenvalues and relative matrix structure itself. That said, your current approach evokes an implictly shifted QR solver that is essentially the standard for this exact type of problem. With some experimentation, however, you may squeak out a modest performance increase.

A few options:

1. If your problem is well conditioned, compute using single precision.

2. DSYEVR is a LAPACK driver for real symmetric matrices that uses a MRRR algorithm to compute the eigenvalues first, then grabs the requested eigenvectors through an inverse-type problem. It is possible for your particular matrix that it could be quicker.

3. Use a parallel library. While your problem is small, I might still expect modest performance gains when using a relatively small number of cores. ScaLAPACK is always a safe bet here, though I would first recommend SLEPc. SLEPc operates through PETSc and serves as a wrapper for several EVP solvers, many in their parallel and thus, usually, scalable form. This will give you the opportunity to try several approaches on the fly and evaluate their relative efficiency for your problem. Iterative algorithms used in SLEPc can take an eigenvalue tolerance parameter that might also improve performance in your case.

Note: I realize this answer isn't especially satisfying. If you provide some details on the nature of your matrix and its construction, someone might impart some kind of specialized knowledge on you.

Answer 2

You might want to try the SelfAdjointEigenSolver class https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html in the Eigen C++ class library http://eigen.tuxfamily.org/index.php?title=Main_Page

I did some numerical experiments with a 30x30 SPD matrix with a condition number of 1000 constructed according to the procedure described by Neumaier here Generating Symmetric Positive Definite Matrices using indices

I used a Windows 7 computer with a 2.8 GHz AMD Phenom II X4 830 CPU. The compiler is VS 2013.

For the experiments with Lapack ssyev/dsyev, I used OpenBlas http://www.openblas.net/. In previous experiments I have found the BLAS operations in this library to have very good performance. I suspect (but did not verify) that their implementation of ssyev/dsyev simply uses the standard Fortran code but that the performance of that routine depends significantly on the quality of the underlying BLAS functions.

Eigen::SelfAdjointEigenSolver was significantly faster than ssyev/dsyev in this case-- ~0.25 ms per call compared with ~2 ms per call.

I saw only very small reduction of the time when I switched from double to single precision. But I noticed large errors in the eigenvalues computed. This was true for both Eigen and Lapack.