-- SamPreston - 08 Nov 2007

Closed-Form Computation of Eigenvalues of a 3x3 matrix

Calculating Eigenvalues of the Hessian

The motivation for this code was that I needed to calculate the eigenvalues and eigenvectors of a Hessian matrix. Since it's only a 3x3 matrix, I though 'how hard can it be'? Well, it turns out to be harder than one might think -- I actually only got around to calculating the eigenvalues before I gave up and used LAPACK to solve the problem (which of course posed its own set of challenges).

General Strategy

After looking around online, I found that the eigenvalues of a matrix could be solved for using the definition of an eigenvalue / eigenvector. If A is a matrix, any non-zero x solving:

Ax = ex

for any non-zero scalar e is an eigenvector of matrix A, and e is its associated eigenvalue. Given this, we can rearrange such that:

(A-Ie)x=0

where I is the identity matrix. We can actually solve for e's using the characteristic polynomial of A, which for a 3x3 matrix is a cubic equation. Looking on wikipedia, there is a closed-form solution for the roots of a cubic equation, but it turns out it's not a complete snap to code -- there are intermediate complex values. In fact, the cube root of a complex number must be calculated. It turns out that this is not as hard as one might imagine, though, because it reduces to a simple division when the complex number is put in polar form, see here. In any case, I whipped up some code to compute the eigenvalues of a matrix, including a stand-alone solver for the roots of cubic equations. However, it's not the best code I've ever written, and I would still need to write code for solving a system of three variables and three equations to find the eigenvectors for the given eigenvalues (solving (A-Ie)x=0 for each e). As I mentioned before, I just decided to use LAPACK, but didn't want the code to get lost in case I (or someone else) needed it in the future, so here it is, just cut directly from the file it was in:

code

// Calculate the coefficients of the characteristic polynomial and sent them to the solver.
// hess is the 3x3 hessian matrix
// eig is the vector that will hold the computed eigenvalues
void NDInterp::getEigenvalues(double hess[3][3], double eig[3]){
  eig[0] = eig[1] = eig[2] = 0;
  // create the characteristic polynomial of the hessian matrix
  double poly[4] = {1, 0, 0, 0};
  poly[1] = -1*(hess[1][1] + hess[2][2] + hess[3][3]);
  poly[2] = 
    hess[1][1]*hess[1][2]
    + hess[1][1]*hess[3][3]
    + hess[2][2]*hess[3][3]
    - hess[2][1]*hess[1][2]
    - hess[1][3]*hess[3][1]
    - hess[2][3]*hess[3][2];
  poly[3] =
    hess[1][1]*hess[2][3]*hess[3][2]
    + hess[1][2]*hess[2][1]*hess[3][3]
    + hess[1][3]*hess[2][2]*hess[3][1]
    - hess[1][1]*hess[2][2]*hess[3][3]
    - hess[1][2]*hess[2][3]*hess[3][1]
    - hess[1][3]*hess[2][1]*hess[3][2];
  // solve for the real roots (which are the eigenvalues)
  double real[3], imag[3];
  int realRoots, numRoots;
  getCubicRoots(poly, real, imag, realRoots, numRoots);
  for(int i=0;i<realRoots;i++){
    eig[i] = real[i];
  }
}

// calculate the roots of the cubic equation with coefficients coeff, where 
// coeff[0] multiplies the x^3 term and coeff[3] is the constant
// rootReal and rootImag returns the real and imaginary parts of each
// root, respectively (that is, rootReal[i] and rootImag[i] form the real and
// imaginary part of a single root), and numRealRoots and numDistinctRoots
// return the number of real roots and the number of distinct roots,
// respectively.  If there are less than 3 distinct roots, it means there
// is a double (or tripple) root, but I don't give any way to figure out
// which root might be the double.  Also, real roots are returned before
// complex roots in the rootReal and rootImag arrays.
void NDInterp::getCubicRoots(double coeffs[4], double rootReal[3], double rootImag[3], int &numRealRoots, int &numDistinctRoots){
  double discriminant;
  double a, b, c, p, q;

  discriminant = 4.0*a*a*a*c - a*a*b*b + 4.0*b*b*b - 18.0*a*b*c + 27.0*c*c;

  a = coeffs[1]/coeffs[0];
  b = coeffs[2]/coeffs[0];
  c = coeffs[3]/coeffs[0];
  p = b-((a*a)/3.0);
  q = c + (2*a*a*a - 9*a*b)/27.0;

  double uReal[6];
  double uImag[6];

  for(int i=0;i<6;i++){
    uReal[i] = uImag[i] = 0.0;
  }

  double uSqrtPart = q*q/4.0 + p*p*p/27.0;

  if(uSqrtPart > 0){
    uReal[0] = sqrt(uSqrtPart);
    uReal[1] = -sqrt(uSqrtPart);
  }else if(uSqrtPart < 0){
    uImag[0] = sqrt(-uSqrtPart);
    uImag[1] = -sqrt(-uSqrtPart);
  }

  uReal[0] += q/2.0;
  uReal[1] += q/2.0;

  getComplexRoot(uReal[0], uImag[0], 3, uReal[0], uImag[0]);
  getComplexRoot(uReal[1], uImag[1], 3, uReal[1], uImag[1]);

  
  double cmplxRootFactor[2];
  cmplxRootFactor[0] = -1.0/2.0;
  cmplxRootFactor[1] = sqrt(3)/2.0;

  complexMult(uReal[0], uImag[0], cmplxRootFactor[0], cmplxRootFactor[1], uReal[2], uImag[2]);
  complexMult(uReal[0], uImag[0], cmplxRootFactor[0], -cmplxRootFactor[1], uReal[3], uImag[3]);
  complexMult(uReal[1], uImag[1], cmplxRootFactor[0], cmplxRootFactor[1], uReal[4], uImag[4]);
  complexMult(uReal[1], uImag[1], cmplxRootFactor[0], -cmplxRootFactor[1], uReal[5], uImag[5]);

  numRealRoots = 0;
  int numRoots = 0;
  int idx=0;
  double ans[2];
  while(numRoots < 3 && idx < 6){
    complexMult(3.0, 0.0, uReal[idx], uImag[idx], ans[0], ans[1]);
    complexDiv(p, 0.0, ans[0], ans[1], ans[0], ans[1]);
    ans[0] -= uReal[idx];
    ans[1] -= uImag[idx];
    ans[0] -= a/3.0;
    std::cout << "found answer " << ans[0] << " + " << ans[1] << "i" << std::endl;
    // add this if it doesn't exist yet
    bool matches = false;
    for(int n=0;n<numRoots;n++){
      if(complexEquals(rootReal[n], rootImag[n], ans[0], ans[1], 0.001)){
   matches = true;
      }
    }
    if(!matches){
      // no match found yet, add this root
      if(fabs(ans[1]) < 0.001){
   std::cout << "adding real root" << std::endl;
   // found a real root, add it to beginning of list
   numRealRoots++;
   for(int n=numRoots;n>0;n--){
     rootReal[n] = rootReal[n-1];
     rootImag[n] = rootImag[n-1];
   }
   rootReal[0] = ans[0];
   rootImag[0] = ans[1];
      }else{
   std::cout << "adding root" << std::endl;
   // found an imaginary root, add it to end of list
   rootReal[numRoots] = ans[0];
   rootImag[numRoots] = ans[1];
      }
      numRoots++;
    }else{
      std::cout << "root matched, skipping" << std::endl;
    }
    idx++;
  }// end while
  numDistinctRoots = numRoots;
  
}

// calculate the 'root'th root of the complex number 'real + imag*i'
void NDInterp::getComplexRoot(double real, double imag, int root, double &rtnReal, double &rtnImag){
  double r = sqrt(real*real+imag*imag);
  double gamma = atan2(imag, real);
  r = pow(r, 1.0/((double)root));
  gamma /= ((double)root);
  rtnReal = r*cos(gamma);
  rtnImag = r*sin(gamma);
}

void NDInterp::complexMult(double r1, double i1, double r2, double i2, double &rtnReal, double &rtnImag){
  double tmpReal, tmpImag;
  tmpReal = r1*r2 - i1*i2;
  tmpImag = r1*i2 + r2*i1;
  rtnReal = tmpReal;
  rtnImag = tmpImag;
}

void NDInterp::complexDiv(double numReal, double numImag, double denReal, double denImag, double &rtnReal, double &rtnImag){
  double tmpReal, tmpImag;
  tmpReal = (numReal*denReal + numImag*denImag)/(denReal*denReal + denImag*denImag);
  tmpImag = (numImag*denReal - numReal*denImag)/(denReal*denReal + denImag*denImag);
  rtnReal = tmpReal;
  rtnImag = tmpImag;
}

bool NDInterp::complexEquals(double r1, double i1, double r2, double i2, double eps){
  if(fabs(r1-r2) < eps && fabs(i1-i2) < eps){
    return true;
  }
  return false;
}