Finding eigenvectors of covariance matrix to create 3D bounding sphere

Question

I'm currently in the process of writing a function to find an "exact" bounding-sphere for a set of points in 3D space. I think I have a decent understanding of the process so far, but I've gotten stuck.

Here's what I'm working with: A) Points in 3D space B) 3x3 covariance matrix stored in a 4x4 matrix class (referenced by cells m0,m1,m2,m3,m4,ect; instead of rows and cols)

I've found the 3 eigenvalues for the covariance matrix of the points, and I've set up a function to convert a matrix to reduced row echelon form (rref) via Gaussian elimination.

I've tested both of those functions against figures in examples I've found online, and they appear to be working correctly.

The next step is to find the eigenvectors using the equation: (M - λ*I)*V

... where M is the covariance matrix, λ is one of the eigenvalues, I is the identity matrix, and V is the eigenvector.

However, I don't seem to be constructing the 4x3 matrix correctly before rref'ing it, as the far right column where the eigenvector components should be calculated are 0 before and after running rref. I understand why they are zero after (without any constants, the simplest solution to a linear system of equations is all coefficients of zero), but I'm at a loss as to what to put there.

Here's the function so far:

Vect eigenVector(const Matrix &M, const float eval) {
   Matrix A = Matrix(M);
   A -= Matrix(IDENTITY)*eval;
   A.rref();
   return Vect(A[m3],A[m7],A[m11]);
}

The 3x3 covariance matrix is passed as M, and the eigenvalue as eval. Matrix(IDENTITY) returns an identity matrix. m3,m7, and m11 correspond to the far-right column of a 4x3 matrix.

Here's the example 3x3 matrix (stored in a 4x4 matrix class) I'm using to test the functions:

Matrix(1.5f, 0.5f, 0.75f, 0,
       0.5f, 0.5f, 0.25f, 0,
      0.75f, 0.25f, 0.5f, 0,
          0,     0,    0, 0);

I'm correctly (?) getting the eigenvalues of 2.097, 0.3055, 0.09756 from my other function.

eigenVector() above correctly subtracts the passed eigenvalue from the diagonal (0,0 1,1 2,2)

Matrix A after rref():

[(1, 0, 0, -0),
(-0, 1, 0, -0),
(-0, -0, 1, -0),
(0, 0, 0, -2.09694)]

For the rref() function, I'm using a translated python function found here: http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html

What should the matrix I pass to rref() look like to get an eigenvector out?

Thanks

Answer 1

(M - λI)V is not an equation, it's just an expression. However, (M - λI)V = 0 is. And it's the equation that relates eigenvectors to eigenvalues.

So assuming your rref function works, I would imagine that you create an augmented matrix as [(M - λI) | 0], where 0 denotes a zero-vector. This sounds like what you're doing already, so I would have to assume that your rref function is broken. Or alternatively, it doesn't know how to handle 4x4 matrices (as opposed to 4x3 matrices, which is what it would expect for an augmented matrix).

Answer 2

Ah, with a few more hours of grueling research, I've managed to solve my problem.

The issue is that there is no "one" set of eigenvectors but rather an infinite number with varying magnitudes.

The method I chose was to use a REF (row echelon form) instead of RREF, leaving enough information in the matrix to allow me to substitute in an arbitrary value for z, and work backwards to solve for y and x. I then normalized the vector to get a unit eigenvector, which should work for my purposes.

My final code:

Vect eigenVector(const Matrix &M, const float eVal) {
   Matrix A = Matrix(M);
   A -= Matrix(IDENTITY)*eVal;
   A.ref();
   float K = 16; // Arbitrary value
   float J = -K*A[m6]; // Substitute in K to find J
   float I = -K*A[m2]-J*A[m1]; // Substitute in K and J to find I

   Vect eVec = Vect(I,J,K);
   eVec.norm(); // Normalize eigenvector

   return eVec;
}

The only oddity is that the eigenvectors come out facing in the opposite direction than I expected (they were negated!), but that's a moot problem.