Finding eigenvectors of 3x3 covariance matrices


An important step in finding the orientation for an OBB is finding the eigenvectors of a 3x3 covariance matrix. A standard method such as Jacobi iterations given in numerical recipes in C is one way to approach the problem, but this routine is rather general and rather expensive. Our efficient top-down tree building routine spends 40% of its time in this eigenvector computation routine. So, we figure we can get a speedup of as much as 1.66 if we can make the eigenvector computation more efficient.

The first thing to note is that we don't need a general routine. Our matrices will always be 3x3 real symmetric matrices. I think they will always be positive semi-definite, but I'm not sure.

The following solution method was suggested in part by leonard mcmillan.

First, the eigenvalues of a matrix will equal the roots of its characteristic equation. For an nxn matrix, the characteristic equation is an n degree polynomial. So, the characteristic equation of our 3x3 covariance matrix will be a cubic, for which the roots can be computed directly.

Now, given a matrix A, if the eigenvalues are along the diagonal of a matrix D, and the eigenvectors are the columns of a matrix E, then the following equation holds:

AE = DE

or, expressed differently,

(A-D)E = 0

If we know the eigenvalues, then D is known, and since A is known, then finding the eigenvectors, E, involves solving a 3x3 linear system. And THAT's not too hard.

Now, remember that this is a homogeneous system, so the solution set is defined to within a scale factor, only. This doesn't complicate things too, much, really.

Consider just one column of E, and call it e. Then we have the expression

(A-D)e = 0

Now set the first element of e to 1, then leave the other two elements as unknowns, and solve the linear system for those. This works great, unless the actual eigenvector has the first element to be zero (in other words, the eigenvector lies in the yz plane).

The fix for this is to examine the system to determine which of the elements of the eigenvector will be largest. THAT is the element which we want to set to 1, and the other two elements are the unknowns for which we solve the system.