Stabilizing a 3x3 real symmetric matrix eigenvalue calculation

Question

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?

Answer 1

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.

Answer 2

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];