# A fast and efficient algorithm for eigenvalues computation of a symmetric positive definite matrix

I am looking for a very fast and efficient algorithm for the computation of the eigenvalues of a 3x3 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

https://en.wikipedia.org/wiki/Eigenvalue_algorithm

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?

Thanks.

• 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 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 at 13:52
• For computing small eigenvalues, Jacobi algorithm is more accurate than QR, DC, and bisection. See lawn15 – piyush_sao Sep 13 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 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$$.
• 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 at 13:58
• 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 at 8:28