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).
