Tell me more ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.
up vote 1 down vote favorite

Given N, 3D position vectors, I want to find the best-fit plane of those vectors.

I found a suggested answer here that said:

  1. Calculate the average of the points (easy)

  2. Calculate the Covariance matrix for points minus the average (thanks)

  3. Find the two largest eigenvectors of that matrix (bwhuh?)

  4. Calculate the Cross product of those and you get the normal of you plane (I can understand that)

So, AFAIK, this is what I think I should be doing:

average=sum(positions)/N
covariant=numpy.cov(positions-average)
eigenvectors,eigenvalues = numpy.linalg.eig(covariant)
sorted_evecs = eigenvectors[eigenvalues.argsort()]
#Take the last EV as per comment instead of crossing largest
normal = sorted_evecs[-1]

But it is plainly wrong, because the 'normal' vector that we get at the end is a 4-vector. Can anyone offer any directions / explanations please?

Sorry to mix SO and Math but thought it would be more appropriate here.

share|improve this question
3  
For 2. you subtract the average from each point before calculating the covariance matrix. This is making a translation of your space to put the average at the origin. It reduces the calculation error. If all your points have coordinates like $10^9\pm 5$, the result will be numbers like $\pm 5$, so when you square them and subtract you won't lose so much precision. – Ross Millikan Jul 21 '12 at 2:47
1  
The covariance matrix is symmetric, so it has three mutually orthogonal eigenvectors; thus taking the cross product of two of them will yield the direction of the third, so you might as well take the eigenvector for the least eigenvalue directly. (By the way, these are not the two "largest eigenvectors", they're the eigenvectors for the two largest eigenvalues.) Other than that, everything seems fine, but I don't know Numpy, so I don't know whether your library calls are right. – joriki Jul 21 '12 at 7:37
1  
The new version of the question makes it entirely a question about the libraries you're using. It seems you've got the math sorted out now; I don't think this is a good place to find out how to get a three-vector from a particular library; though you might get lucky. My wild guess would be that the library is returning vectors that contain both the eigenvector itself and, in an additional component, the eigenvalue; in that case you'd just have to figure out which three of the four components form the eigenvector. – joriki Jul 22 '12 at 13:16
By the way, I suspect that "eigenvalues" and "evals" were supposed to be the same variable? – joriki Jul 22 '12 at 13:17

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.