4
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • $\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
    Commented Sep 13, 2019 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
    Commented Sep 13, 2019 at 13:52
  • 2
    $\begingroup$ For computing small eigenvalues, Jacobi algorithm is more accurate than QR, DC, and bisection. See lawn15 $\endgroup$
    – arbitUser1401
    Commented Sep 13, 2019 at 17:36
  • 5
    $\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ý
    Commented Sep 13, 2019 at 18:45

3 Answers 3

Reset to default
6
$\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.

Improve this answer
$\endgroup$
6
  • $\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
    Commented Sep 13, 2019 at 12:59
  • 3
    $\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
    Commented Sep 13, 2019 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
    Commented Sep 13, 2019 at 15:54
  • 1
    $\begingroup$ Of course, if you would write the WP in a compiled language, it would probably be faster... $\endgroup$
    – GertVdE
    Commented Sep 16, 2019 at 8:47
  • 1
    $\begingroup$ I implemented the Wikipedia algorithm in Fortran, taking half time in respect to LAPACK over 10^6 3x3 matrices. $\endgroup$
    – John Snow
    Commented Sep 16, 2019 at 12:09
0
$\begingroup$

This paper gives the equations for eigenvalues and eigenvectors of 3x3 posdef matrices. But they do not discuss the numeric accuracy of these equations.

Improve this answer
$\endgroup$
-1
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • 4
    $\begingroup$ 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 ). $\endgroup$
    – Federico Poloni
    Commented May 29, 2020 at 10:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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