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

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.

• Welcome to SciComp.SE! I'd say your question is on topic here, but as far as I know, the QR algorithm (not factorization!) is pretty much the only game in town for computing all eigenvectors of a dense matrix. Mar 10 '16 at 23:49
• How large is large? I can think of some really big numbers. Mar 11 '16 at 1:35
• See the discussion of Francis's algorithm in David Watkins numerical linear algebra textbook. This algorithm is commonly known as the QR algorithm, but Watkins makes a good case for not using that name, since it isn't very descriptive of how the method is actually implemented in practice in its implicitly shifted version. Mar 11 '16 at 2:03
• @BrianBorchers Good point, also since it cuts down on students' inevitable confusion between the QR factorization (used, e.g., for least squares) and the QR algorithm. I'll probably adopt this nomenclature in my class next semester. Mar 11 '16 at 8:18
• For larger matrices, you should also consider the divide and conquer algorithm. If you need all eigenvectors, div-and-conquer might well outperform (shifted) QR. You may also want to look at some of the references that discuss the latest LAPACK eigen-solver implementations, for example here, here or here. Mar 11 '16 at 20:36

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.

• Thanks for your detailed answer. I do computational chemistry and eigenproblems of large symmetric matrices take up a significant percentage of total CPU-time, sometimes the vast majority, so I am interested in any development that may help me, when it trickles down from mathematical research to a usable implementation in a library. Mar 14 '16 at 16:02
• Jacobi algorithm is definitely not faster than QR algorithm, except where matrix is already very close to diagonal. Mar 16 '16 at 16:35
• @zimbra314 I agree - both are niche algorithms. I hope I did not give the wrong impression in my answer. Mar 16 '16 at 19:00