What is the state of the art algorithm for diagonalizing real symmetric matrices?

Question

There are many methods for diagonalizing matrices, probably the most widely used is the combination of Householder transformations and the QR algorithm.

Is there any superior method for diagonalizing (large, non-sparse)real symmetric matrices? Superiority can be a bit muddy, so I define it as fast, numerically stable, does not require large amounts of extra memory and lends itself to parallelization and vectorization.

Meta: I have been a bit torn about the correct place to ask this, another candidate would have been the Mathemathics exchange, they also have a numerical linear algebra tag.

Answer 1

The QR/Francis algorithm is the go-to choice for dense eigenproblems, but there are a few competitors around:

1. The Jacobi algorithm (like QR, another algorithm with an unfortunate name, which can be confused with a method to solve linear systems...). Main idea: choose a nondiagonal element of $A_k$; apply a suitable 2x2 Givens rotations $A_{k+1}=Q_kA_kQ_k^*$ to zero it out; then repeat. In each new iteration the entries that you have zeroed out at the previous steps become nonzero again, so the matrix is always full; nevertheless, one can prove that the off-diagonal entries converge to zero under some assumptions on how selected them in the different step. Literature pointer studying stability: Demmel, Veselic 1992, Jacobi's method is more accurate than QR. Essentially, they show that when the eigenvalues span several orders of magnitude, this method gets better accuracy for the smaller ones. Note that between 1992 and now there has been another major improvement to the QR algorithm (the MRRR method by Parlett and Dhillon); I am not expert enough in the area to tell if this has bridged the gap.

2. Spectral divide-and-conquer algorithms based on doubling iterations (which are not the classical divide and conquer algorithms mentioned in a comment). These algorithms rely on BLAS level-3 operations only, so (given an implementation of products, inversions, projections etc., which pays proper attention to communication and caching issues) they should scale better to matrices of dimension larger than a few thousands. Main idea: apply an iteration such as Newton for the matrix sign $M_0=A$, $M_{k+1}=\frac12 (M_k + M_k^{-1})$ under which the eigenvectors are preserved and each eigenvalue converges to either $-1$ or $1$; then compute kernels of $(\lim M_k)\pm I$ to determine two complementary invariant subspaces of $A$; project $A$ onto these two subspaces, and repeat for the two smaller problems. Literature pointer studying stability: Nakatsukasa and Higham, 2012, Stable and Efficient Spectral Divide and Conquer Algorithms for the Symmetric Eigenvalue Decomposition and the SVD. They construct a variant of the iteration which requires no matrix inverses and converges extremely fast, and prove the stability of the resulting method.

Both methods can be faster / more accurate than QR under certain conditions; there is research interest in the numerical linear algebra community, but they are less mature than QR, which has been the main horse in the race for decades, used by virtually everyone in numerical computing. I would not recommend their use instead of QR to a practitioner in a different field who just want to get the job done; but if you are doing research on eigenvalue methods or you have a problem for which QR encounters trouble you may want to experiment with them.

As for implementations, LAPACK includes a Jacobi routine to compute the SVD, *GESVJ, and a Matlab implementation of QDWH (the algorithm described in the Nakatsukasa and Higham's paper) is on Matlab file exchange.