3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Brian Borchers Sep 13 at 13:51
  • 2
    $\begingroup$ 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. $\endgroup$ – Brian Borchers Sep 13 at 13:52
  • $\begingroup$ For computing small eigenvalues, Jacobi algorithm is more accurate than QR, DC, and bisection. See lawn15 $\endgroup$ – piyush_sao Sep 13 at 17:36
  • 1
    $\begingroup$ 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... $\endgroup$ – Jakub Klinkovský Sep 13 at 18:45
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – John Snow Sep 13 at 12:59
  • 2
    $\begingroup$ 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. $\endgroup$ – GertVdE Sep 13 at 13:58
  • $\begingroup$ I tested the DSYEV and compared to the Wikipedia algorithm. The latter takes exactly half of the time. However, thank a lot @GertVdE. $\endgroup$ – John Snow Sep 13 at 15:54
  • $\begingroup$ 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). $\endgroup$ – GertVdE Sep 16 at 8:28
  • 1
    $\begingroup$ I implemented the Wikipedia algorithm in Fortran, taking half time in respect to LAPACK over 10^6 3x3 matrices. $\endgroup$ – John Snow Sep 16 at 12:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.