# Stabilizing a 3x3 real symmetric matrix eigenvalue calculation

I have many 3x3 real symmetric matrices for which I need to determine the eigenvalues. Wikipedia gives a nice non-iterative algorithm for this case, which I have translated into C++:

#include <iostream>
#include <vector>
#include <cmath>
#include <boost/math/constants/constants.hpp>

template<typename Real>
std::vector<Real> eigenvalues(const Real A[3][3])
{
using boost::math::constants::third;
using boost::math::constants::pi;
using boost::math::constants::half;

static_assert(std::numeric_limits<Real>::is_iec559,
"Template argument must be a floating point type.\n");

std::vector<Real> eigs(3, std::numeric_limits<Real>::quiet_NaN());
auto p1 = A[0][1]*A[0][1] + A[0][2]*A[0][2] + A[1][2]*A[1][2];
auto diag_sq = A[0][0]*A[0][0] + A[1][1]*A[1][1] + A[2][2];
if (p1 == 0 || 2*p1/(2*p1 + diag_sq) < std::numeric_limits<Real>::epsilon())
{
eigs[0] = A[0][0];
eigs[1] = A[1][1];
eigs[2] = A[2][2];
return eigs;
}

auto q = third<Real>()*(A[0][0] + A[1][1] + A[2][2]);
auto p2 = (A[0][0] - q)*(A[0][0] - q) + (A[1][1] - q)*(A[1][1] -q) + (A[2][2] -q)*(A[2][2] -q) + 2*p1;
auto p = std::sqrt(p2/6);
auto invp = 1/p;
Real B[3][3];
B[0][0] = A[0][0] - q;
B[0][1] = A[0][1];
B[0][2] = A[0][2];
B[1][1] = A[1][1] - q;
B[1][2] = A[1][2];
B[2][2] = A[2][2] - q;
auto detB = B[0][0]*(B[1][1]*B[2][2] - B[1][2]*B[1][2])
- B[0][1]*(B[0][1]*B[2][2] - B[1][2]*B[0][2])
+ B[0][2]*(B[0][1]*B[1][2] - B[1][1]*B[0][1]);
auto r = invp*invp*invp*half<Real>()*detB;
if (r >= 1)
{
eigs[0] = q + 2*p;
eigs[1] = q - p;
eigs[2] = 3*q - eigs[1] - eigs[0];
return eigs;
}

if (r <= -1)
{
eigs[0] = q + p;
eigs[1] = q - 2*p;
eigs[2] = 3*q - eigs[1] - eigs[0];
return eigs;
}

auto phi = third<Real>()*std::acos(r);
eigs[0] = q + 2*p*std::cos(phi);
eigs[1] = q + 2*p*std::cos(phi + 2*third<Real>()*pi<Real>());
eigs[2] = 3*q - eigs[0] - eigs[1];

return eigs;
}

int main()
{
float M[3][3];
M[0][0] = 1.25e6;
M[1][0] = 1.25e6;
M[0][1] = 1.25e6;
M[1][1] = 1.25e6;
M[0][2] = 0.0;
M[2][0] = 0.0;
M[1][2] = 0.0;
M[2][1] = 0.0;
M[2][2] = 0.0;

auto eigs = eigenvalues<float>(M);
std::cout << "Eigenvalues: ";
std::cout.precision(std::numeric_limits<float>::digits10 + 5);
std::cout << std::fixed << std::scientific;
for (auto eig : eigs)
{
std::cout << eig << ", ";
}
std::cout << '\n';
}


The eigenvalues of the matrix are $2.5\times10^6$, 0, and 0. However, the program returns $2.5\times10^6$, $0.0625$, and $0$. Yes, the ratio of the second to the first is roughly the float epsilon, and $q$ and $p$ are nearly equal. But is there a way to stabilize this algorithm so that the loss of precision is not so dramatic?

• What happens if you switch to double precision floats? How precise do you need the eigenvalues to be? Apr 14, 2016 at 3:21
• In double precision, the error is still much larger than the double epsilon; same for long double. The issue is not really that the eigenvalues are not precise, it's that the subsequent eigenvector calculations become unstable. This routine is at the very bottom of a deep function call stack, and I simply cannot reason clearly about how the error will propagate. Hence I would like the result to be accurate to within the <Real>::epsilon(). Apr 14, 2016 at 13:51
• What do you need the eigenvalues for, if I may ask? Apr 16, 2016 at 14:20
• @FedericoPoloni: This is for subsurface imaging. At each point in the ground, you get a different Hooke's law (81 component symmetric rank-3 tensor) then do a tensor contraction with the direction you are interested in to create the 3x3 Christoffel matrix, whose eigenvalues are the squares of the phase velocity of the waves (qP, qSH, qSV) in that particular direction. So I'm doing lots of these calculations. Apr 17, 2016 at 22:51
• Ok, so your final outcome are the eigenvalues? In that case I am afraid there aren't many shortcuts. If you needed the eigenvalues to do something else with them, then maybe there would have been a different way to do it without them... Apr 18, 2016 at 7:13

This is trying to compute the eigenvalues by computing the roots of the characteristic polynomial. In this case, the characteristic polynomial is $p(t) = t^3-2t^2x$, $x=1.25\times 10^6$, and zero is a root of multiplicity two.

There is generally no good way of computing double roots of a polynomial to precision better than $O(\sqrt{\epsilon})$, because if the polynomial is specified by its coefficients (as here) there exists a perturbation in the coefficients of magnitude $O(\epsilon)$ that causes the double root to split into two roots separated by $O(\sqrt{\epsilon})$. The same should apply to perturbations of $O(\epsilon)$ in the input matrix too.

The usual advice is that having a backward-stable algorithm (one which computes the exact result for a perturbed input) is the goal, and if you have a problem that is so sensitive to tiny perturbations in the input then you are probably asking the wrong question, although there are some exceptions. So I'd say the important question here is why it is important for you to accurately compute double eigenvalues of a matrix.

• I was hoping for an analogue to Kahan summation; and I think there is still hope for such a thing. I do indeed have a good reason for wanting the double eigenvalues to high precision (the physics I'm modelling changes at the double root); $\mathcal{O}(\epsilon)$ error is fine but $\mathcal{O}(\sqrt{\epsilon})$ is too much . . . Apr 14, 2016 at 15:40
• @NickThompson The worst you can have is a triple root, so try boost's mpfr bindings (boost.org/doc/libs/1_60_0/libs/multiprecision/doc/html/…) and use three times as much precision (or some kind of double-double arithmetic library). That will probably give you a perfectly rounded answer, saving effort and time. In the end, Kahan summation is a little bit like emulation of higher precision arithmetic. Apr 14, 2016 at 16:33
• @NickThompson My point about $O(\sqrt{\epsilon})$ is that the kind of answer you are looking for doesn't really exist: it is not just sensitive to roundoff, it is extremely sensitive to the problem's definition itself. See for example Trefethen&Bau or Higham's Accuracy and Stability of Numerical Algorithms for discussions of this. Apr 14, 2016 at 16:35
• I still feel like the judicious use of a temporary could dramatically improve the result; especially if there's a way to make it emulate higher precision arithmetic. However, if I can't get it to work in a couple days I'll concede defeat and accept this answer. Apr 14, 2016 at 16:52
• If I get really ambitious I'll even read your references . . . Apr 14, 2016 at 16:54

Please check this line, the A[2][2] should be squared: auto diag_sq = A[0][0]*A[0][0] + A[1][1]*A[1][1] + A[2][2];

• Unfortunately, this does not provide an answer to the question being asked. Though your note is correct, in this particular case, A[2][2]=0., so it won't influence the final result in in a given case. I suggest moving it to a comment (a very good one!) as opposed to an answer. Jun 6, 2018 at 0:13
• I don't have enough reputation to comment ;-) Jun 6, 2018 at 0:42
• If I'm not mistaken, the 3rd line of the determinant should be + B[0][2]*(B[0][1]*B[1][2] - B[1][1]*B[0][2]); Jun 6, 2018 at 0:43