# Fast and accurate eigenvalue computation for 3x3 posdef matrices

I am looking for a very fast and efficient algorithm for the computation of the eigenvalues of a $$3\times 3$$ symmetric positive definite matrix. The algorithm will be part of a massive computational kernel, thus it is required to be very efficient.

I am aware of the algorithm suggested by Wikipedia but I found this strategy not sufficiently robust. In particular, the Wikipedia algorithm often finds slightly negative eigenvalues even if the matrix is positive definite. Suggestions?

• If your matrices are positive semidefinite but singular, then any floating-point computation of the eigenvalues is likely to produce small negative eigenvalues that are effectively 0. You should be looking for ways to make the higher level computation deal with this eventuality. – Brian Borchers Sep 13 '19 at 13:51
• If the algorithm you're using is producing negative eigenvalues for matrices that are in fact strictly positive definite, then clearly that algorithm is broken. – Brian Borchers Sep 13 '19 at 13:52
• For computing small eigenvalues, Jacobi algorithm is more accurate than QR, DC, and bisection. See lawn15 – piyush_sao Sep 13 '19 at 17:36
• Can you give an example for which the cited algorithm does not give the expected result? A symmetric 3x3 matrix should be sufficiently small to be posted here... – Jakub Klinkovský Sep 13 '19 at 18:45

For a symmetric 3x3 matrix, one Householder transformation will bring your matrix in tridiagonal form. The required algorithm is given (for general $$n\times n$$ matrices) on page 459 of Matrix Computations, 4th edition, Algorithm 8.3.1. For a $$3\times 3$$ matrix, it's just one Householder reduction instead of a loop.

For the subsequent tridiagonal matrix, you can apply the implicit shift symmetric QR algorithm (see Algorithm 8.3.3 p. 463, Matrix Computations, 4th edition) which again you could unroll for $$n=3$$.

But before you unroll and code your own, just benchmark the dsyev routine from LAPACK, preferably through an optimized library as MKL.

• Sorry @GertVdE. I tried to look at LAPACK. The DSYEV routine works very well but it is too general and too computationally expensive for the case. Do you have an smaller example code? – John Snow Sep 13 '19 at 12:59
• Is your comment above based on actual benchmarking or gut-feeling? I don't want to sound negative, but rolling your own code is going to take time (program/debug) and beating MKL is hard. I don't have any code for smaller sizes. The best you can do is write down the algorithms I mentioned for the specific $3\times 3$ case (and test/compare to LAPACK). By the way, is there any chance that next year you will have to revisit your code because the application requires a $5\times 5$ matrix? If so, I wish you good luck in re-interpreting the code you wrote for $3\times 3$ and "expand" it. – GertVdE Sep 13 '19 at 13:58
• I tested the DSYEV and compared to the Wikipedia algorithm. The latter takes exactly half of the time. However, thank a lot @GertVdE. – John Snow Sep 13 '19 at 15:54
• For what it's worth: I implemented the Wikipedia code in Python and compared to scipy.linalg.eigh which directly calls LAPACK. The WP took 239 time units, eigh took 181 time units (ran 10000 random 3x3 matrices). – GertVdE Sep 16 '19 at 8:28
• I implemented the Wikipedia algorithm in Fortran, taking half time in respect to LAPACK over 10^6 3x3 matrices. – John Snow Sep 16 '19 at 12:09

For a 3x3 matrix, the analytical solution can easily be computed with some symbolic math software like wmaxima.

• The algorithm OP found on Wikipedia is, essentially, an analytical solution, and yet numerically it is not satisfying. Do you have reasons to think that the analytical solution computed by wxmaxima does not have this issue? Just because it is a closed-form formula instead of an iterative algorithm, it does not make it automatically more accurate when using machine precision; this is well-known also in simpler examples such as solving quadratic equations (see e.g. people.csail.mit.edu/bkph/articles/Quadratics.pdf ). – Federico Poloni May 29 at 10:13