I've been fiddling with Sylvester's Equation in quaternions: ax-xb=c, for a,b, s.t. aa'=bb' and a+a'=b+b'. Consider the linear function f(x)=ax-xb. Note that (a-b)+(a-b)' = a - b + a' - b' = a+a' - (b+b') = a+a' - (a+a') = 0 So (a-b) = -(a'-b'). I contend that (a-b) is an *eigenvector* of f(x), i.e., f(a-b) = lambda (a-b) for some 'eigenvalue' lambda. Let's try (a-b) out: f(a-b) = a(a-b) - (a-b)b = a(a-b) - (-(a'-b'))b = a(a-b) + (a'-b')b = a(a-b) + a'b - b'b = a(a-b) + a'b - a'a = a(a-b) + a'(b-a) = a(a-b) - a'(a-b) = (a-a')(a-b) ; lambda = (a-a') So, (a-b) is an eigenvector for f(x) for the 'eigenvalue' (a-a'). But (a-a') is a (pure) quaternion *vector* -- i.e., it isn't a real or complex number at all ! Pure vectors play the part in quaternions of pure imaginaries in complex numbers. For example, the square of a pure vector v is -N(v), so pure vectors are the square roots of negative real numbers. We notice that if this Sylvester expression ax-xb is represented in *matrix* form, it is a skew symmetric matrix whose eigenvalues are (surprise, surprise) *pure imaginaries*. Indeed, the absolute value of these imaginaries are identical to the absolute value of (a-a'). So if we 'unlift' (a-a') and -(a-a') back into the complex numbers, they are indeed the eigenvalues of the matrix form of the Sylvester equation. What is the other eigenvector for ax-xb ? We consider f(f(x)): a(a-a')(a-b) - (a-a')(a-b)b = (a-a')a(a-b) - (a-a')(a-b)b = (a-a') (a(a-b) - (a-b)b) = (a-a') (a-a')(a-b) = -N(a-a') (a-b) So -N(a-a') is indeed a *real number* eigenvalue for f(f(x)) with eigenvector (a-b), which follows from the corresponding theorem that the *square* of a skew symmetric matrix is a symmetric negative semidefinite matrix. We have also proved that (a-a')(a-b) is also an eigenvector, so we can span our 2D range space. As a bonus, our eigenvectors are *orthogonal*. dot(q,r) = (qr'+rq')/2, so q,r are orthogonal iff qr'+rq' = 0. dot((a-a')(a-b),(a-b)) = (a-b)[(a-a')(a-b)]'+(a-a')(a-b)(a-b)' = (a-b)(a-b)'(a-a')' + (a-a')(a-b)(a-b)' = N(a-b)(a-a')' + (a-a')N(a-b) = N(a-b) ((a-a')' + (a-a')) = N(a-b) (a' - a + a - a') = N(a-b) 0 = 0 So, (a-b) and (a-a')(a-b) form our orthogonal eigenbasis. Finally, to solve ax-xb=c, we have to express c as a linear combination of our eigenbase vectors, which task is simplified due to their orthogonality.