home

Eigendecomposition for a 2x2 matrix

Believe it or not, doing the eigen-decomposition of a 2x2 symmetric matrix is not as trivial as it seems. While I was doing my M.Eng thesis, at a certain point I was implementing two algorithms:

• the eig_symm_* algorithm of the GNU Scientific Library.

• the closed-form algorithm described in:

  Jim Blinn, ”Consider the Lowly 2x2 Matrix”, IEEE Computer Graphics and Applications, March 1996, pp. 82-88.

And I found bugs in both! GSL would fail on simple cases: for diagonal matrices it would return two equal eigen-vectors. Blinn’s algorithm has a couple of minor bugs.

So, without further ado, here is the algorithm I came up with. It’s Blinn’s algorithm with a couple of corrections.

void solve_eig(double A, double B, double C, double D, 
    double* lambda1, double *v1x, double*v1y, 
    double* lambda2, double *v2x, double*v2y ) 
{
    if(B*C <= 0.1e-20 ) {
        *lambda1 = A; *v1x = 1; *v1y = 0;
        *lambda2 = D; *v2x = 0; *v2y = 1;
        return;
    }

    double tr = A + D;
    double det = A * D - B * C;
    double S = sqrt( square(tr/2) - det );
    *lambda1 = tr/2 + S;
    *lambda2 = tr/2 - S;

    double SS = sqrt( std::max(square((A-D)/2) + B * C, 0.0) );
    if( A - D < 0 ) {
        *v1x = C;
        *v1y = - (A-D)/2 + SS;
        *v2x = + (A-D)/2 - SS;
        *v2y = B;
    } else {
        *v2x = C;
        *v2y = - (A-D)/2 - SS;
        *v1x = + (A-D)/2 + SS;
        *v1y = B;
    }

    double n1 = sqrt(square(*v1x)+square(*v1y));
    *v1x /= n1; *v1y /= n1;
    double n2 = sqrt(square(*v2x)+square(*v2y));
    *v2x /= n2; *v2y /= n2;
}