How does one find the axis of rotation for a pure rotation matrix when said matrix is also symmetric

eigenvalues-eigenvectorslinear algebramatricessymmetric matrices

I'm a programmer working on an advanced C++ 3D mathematics library.

Now, things have been going well, in fact, basically everything about the library has been fully implemented, but one final bit of code still alludes me: finding the axis of rotation for a pure rotation matrix when said matrix is also symmetric.

I've got a nice bit of math going when it comes to non-symmetric matrices

Given a non-symmetric 3x3 pure rotation matrix [M]

M = { { a, b, c },
      { d, e, f },
      { g, h, i } }

det(M) = 1

M * T(M) = T(M) * M = I

M =/= T(M)

an eigenvector [u] which sits along the axis of rotation can be found 

u = { h - f,
      c - g,
      d - b }

such that its normal is axis of rotation [r] of the matrix

r = u / |u|

but this math breaks the moment you give it a symmetric matrix, as the 'h – f', 'c – g', and 'd – b' portions will all resolve to zero, which obviously isn't the normal vector I want.

Now, I do understand linear algebra, but only a little bit. I've been researching this problem for a few days now, and while there are resources that talk about this, most of them either don't address the problem I'm having, or explain it in a way that my scrub brain simply can't keep up with.

They tell me to do things like 'diagonalize M and solve for u', but I've not a clue what that actually entails doing, let alone in a generalized way, let alone (even more) teaching my C++ library to do it in a generalized way given any symmetric pure rotation matrix.

So ya, that's my plight. Hoping that one of y'all could help bail me out on this one and show me how to solve this problem. 😀

Btw, again, this is needed for writing code, so if your answer could be written in a way that addresses that need and also the fact that I'm a linear algebra noobie, that would be awesome.

Thanks in advance!

Best Answer

If the matrix is symmetric, then $M^TM=M^2=I$, so it’s a 180-degree rotation about some axis. The method you’re using is known to fail for this angle, as you’ve discovered. (Indeed, among the conditions in the gray block in your question is $M\ne M^T$.) Essentially, it examines the skew-symmetric part of $M$, which is $\sin\theta$ times the “cross-product matrix” of the rotation axis. This unfortunately vanishes when $\theta=\pi$.

This question and this one describe a different method that uses the symmetric part of $M$ instead: Compute $$T = M+M^T-(\operatorname{tr}(M)-1)I.$$ Any nonzero row or column of this matrix is a vector parallel to the rotation axis.

Addendum (2019.11.19): As you noted, “diagonalize $M$ and solve for $u$” isn’t particularly helpful advice: most rotation matrices aren’t diagonalizable over the reals in the first place. The rotation axis is fixed by the rotation, though, so what you can do for any rotation matrix is compute the eigenspace of $1$, i.e., find a null vector of $M-I$. However, the method in your question and the above method are more efficient. Because of roundoff and other computational errors, in practice computing a null vector of a matrix is done by computing its SVD and taking the singular vector that corresponds to the least singular value.

Related Question