No, not *that* (politician) Hamilton! [Although the following derivation requires only undergraduate level linear algebra, I've never seen this derivation before; I've Google'd and ScholarGoogle'd, but nothing of this sort shows up.] SVD's Solve Linear Equations in Quaternions The linear mapping captured by a 4x4 matrix M of real numbers can be computed by a 4-term linear equation in quaternions: M(Z) = AZB + CZD + EZF + GZH But how can we constructively compute A,B,C,D,E,F,G,H from this 4x4 matrix M? The most significant step is Beltrami's: compute the SVD of M: M = USV' where U,V are unitary, and S is a non-negative diagonal matrix. Actually, we want to tweak our SVD slightly by making |U|=|V|=+1, i.e., U,V are *special unitary*, and allow S to have negative numbers. The next step is to factor both U & V' into left and right quaternions [Van Elfrinkhof & 'Isoclinic decomposition' -- Wikipedia]: If Z is interpreted as a column vector, then UZ is the column vector rotated by U. If Z is interpreted as a quaternion, then we can find quaternions Ul & Ur, such that the quaternion (Ul)Z(Ur) is the same quaternion as (UZ), interpreted as a quaternion. Similarly, we can find Vl & Vr, such that (Vl)Z(Vr) is the same as (V'Z). So, once we have tweaked our SVD to utilize only special unitary matrices, the implementation of the pre and post rotations by means of quaternions is straight-forward. So now comes the harder part: how to implement the equivalent of a diagonal matrix using quaternions? The answer is Hadamard's 4x4 'H2' matrix: [1 1 1 1] [1 -1 1 -1] H2 = (1/2) [1 1 -1 -1] [1 -1 -1 1] If the diagonal entries of S are [s1,s2,s3,s4], then [s1 s2 s3 s4] H2 = [a b c d] and our linear quaternion expression implementing (SZ) is: f(Z) = (a*Z - c*iZi - b*jZj -d*kZk)/2 Here i,j,k are Hamilton's basic quaternions: i^2 = j^2 = k^2 = -1 = ijk. Some examples. Extracting the 'scalar part' of Z: s1=1,s2=s3=s4=0, so [a b c d] = [1 -1 -1 -1]/2 s(Z) = (Z - iZi - jZj - kZk)/4 Extracting the 'vector part' of Z: s1=0,s2=s3=s4=1, so [a b c d] = [3 1 1 1]/2 V(Z) = (3Z + iZi + jZj + kZk)/4 Computing the 'conjugate' of Z: s1=1,s2=-1,s3=-1,s4=-1, so [a b c d] = [-1 -1 -1 -1] K(Z) = -(Z + iZi + jZj + kZk)/2 But here's the really cool part: so long as diagonal matrices have no zero diagonal elements, then they are trivially invertible with diagonal elements: [1/s1 1/s2 1/s3 1/s4] So we can trivially compute the quaternion representation for (1/M) by using 1/s1, 1/s2, 1/s3, 1/s4 instead of s1, s2, s3, s4. But this means we can find an explicit representation for the inverse of the quaternion linear function: f(Z) = (a*Z - c*iZi - b*jZj -d*kZk)/2 If g(f(Z))=Z=f(g(Z)), then we have g(Z) = (e*Z - f*iZi - g*jZj - h*kZk)*2, where e = (a^3 - a(b^2+c^2+d^2) + 2bcd)/P f = (c^3 - c(a^2+b^2+d^2) + 2abd)/P g = (b^3 - b(a^2+c^2+d^2) + 2acd)/P h = (d^3 - d(a^2+b^2+c^2) + 2abc)/P and P = (a+b+c+d)(a-b+c-d)(a+b-c-d)(a-b-c+d) -------------------------------------------- To summarize, given an arbitrary 4x4 matrix M of real numbers, we can produce a quaternion linear expression with 4 terms which represents the same mapping, thus: MZ ~ M(Z) = AZB+CZD+EZF+GZH M(Z) = (Vl) f((Ul)Z(Ur)) (Vr) = (Vl)[a*(Ul)Z(Ur) - c*i(Ul)Z(Ur)i - b*j(Ul)Z(Ur)j - d*k(Ul)Z(Ur)k](Vr)/2 = [a*(VlUl)Z(UrVr) - c*(VliUl)Z(UriVr) - b*(VljUl)Z(UrjVr) - d*(VlkUl)Z(UrkVr)]/2 i.e., A = a*(Vl)(Ul)/2 B = (Ur)(Vr) C = -c*(Vl)i(Ul)/2 D = (Ur)i(Vr) E = -b*(Vl)j(Ul)/2 F = (Ur)j(Vr) G = -d*(Vl)k(Ul)/2 H = (Ur)k(Vr) If we can explicitly represent the inverse of f(Z), then we can trivially solve equations f(Z)=R by computing Z0 = f^-1(R).