FAQ -- How to solve the quaternion equation AZ + ZC  = E for Z.

Copyright (c) 1996 by Henry G. Baker.

When working with quaternions, we must often solve 'linear' equations.
The equation AZ=E is trivial to solve, yielding the solution Z=(1/A)E.
Likewise, the equation ZC=E is trivial to solve, yielding the solution
Z=E(1/C).

Linear quaternion equations with _two_ terms on the left hand side,
however, are a bit harder.  Equations like AZ+CZ=E succumb via
distribution: (A+C)Z=E.  Similarly, ZA+ZC=E yields via distribution:
Z(A+C)=E.

However, equations like AZB+DZC=E are quite difficult.  We can
simplify this to ZB(1/C)+(1/A)DZ=(1/A)E(1/C), or to
(1/D)AZ+ZC(1/B)=(1/D)E(1/B), but then we are stuck, because the
non-commutivity of quaternion multiplication makes it difficult to
simplify any further.

We can, of course, always solve the quaternion equation AZ+ZC=E by
converting the quaternions into matrices, and then solving the matrix
equation (see References).  While this is not terribly expensive, it
doesn't provide as much insight as a more direct solution.

The orthogonal 3-D unit vectors i, j and k span the 3-D space, and
also compose via (noncommutative) quaternion multiplication like so:

i^2 = j^2 = k^2 = -1, ij = k, jk = i, ki = j.

The characteristic equation, whose roots are the eigenvalues, for a
quaternion Q is:

z^2 - 2 Re(Q) z + |Q|^2,

where Re(Q) is the scalar part of Q, and where
|Q|^2 = QQ' = norm(Q) is the square of the 4-D 'length' of Q.

We know that _every_ quaternion A can be expressed in the form:

A = alpha + beta m,

for some real alpha, beta and some 3-D unit vector m.

Therefore, rotate the coordinate system until m coincides with the 3-D
unit vector i -- i.e.,

A = alpha + beta i.

(This rotation can be done explicitly by using the substitution Z = QZ"Q',
where Q is the quaternion that rotates the 3-D unit vector m into the
3-D unit vector i, then solving for Z", and then substituting into the
above equation for Z.  Alternatively, we can simply find a vector n
which is orthogonal to m, and define i=m, j=n.)

We now express all of the constants A,C,E and the variable Z as pairs
of _complex_ numbers, where the unit vector i becomes the complex
number i.  Thus,

A = a + 0j, C = c + dj, E = e + fj, Z = x + yj, a,c,d,e,f,x,y all complex.

Note that by construction, the 'j-part' of A is zero, and
a = alpha + beta i.  Note also that a' denotes the complex conjugate of
a -- e.g.,  a' = alpha - beta i.  Finally, note that

ja = j(alpha + beta i) = alpha j - beta i j = a'j,

so that interchanging j with a complex number conjugates that complex
number.

Given such orthogonal unit 3-D vectors i,j, we can decompose a given
quaternion Z=X+Yj into its two 'halves' by using the following
equations:

X = (Z - iZi)/2, Y = -(Z + iZi)j/2

We can now expand our equation:

AZ+ZC=E => (a+0j)(x+yj)+(x+yj)(c+dj)=(e+fj)

        => ax + ayj + xc + xdj + yjc + yjdj = e + fj

        => ax + xc + yjdj + xdj + ayj + yjc = e + fj

        => ax + cx + yd'j^2 + dxj + ayj + yc'j = e + fj

        => (a+c)x - d'y + dxj + (a+c')yj = e + fj

        => [(a+c)x-d'y] + [dx+(a+c')y]j = e + fj

But we can now _separate_ the single equation into two equations:

(a+c)x - d'y = e

dx + (a+c')y = f

These equations are easily solved via Cramer's rule (in complex numbers):

    | e   -d'|
    | f  a+c'|
x = ----------
    |a+c  -d'|
    | d  a+c'|

    |a+c   e |
    | d    f |
y = ----------
    |a+c  -d'|
    | d  a+c'|

So, since the

denominator = a^2 + (c+c')a + cc' + dd'

            = a^2 + 2Re(C)a + |C|^2

            = substitution of A (= a) into characteristic eqn for -C !

Thus, the denominator becomes zero iff the eigenvalues of A & -C overlap.

    | e   -d'|   |a+c  e |
    | f  a+c'| + | d   f | j     (a+c')e+d'f + [(a+c)f-de] j
Z = ------------------------  =  ---------------------------
          |a+c  -d'|               A^2 + 2 Re(C) A + |C|^2
          | d  a+c'|

For quaternions, if the eigenvalues overlap, then they are identical.
If the eigenvalues are identical, then the two quaternions are simple
rotations of one another, and they can be brought into conjunction via
a rotation U.  I.e., if quaternions G,H have identical eigenvalues,
then they are related by a third unit quaterion U, such that G = UHU'.

Going back to our equation AX+XC=E, if C = -U'AU, for some unit
quaternion U, then our equation becomes

AZ+ZC = AZ-ZU'AU = E.

Multiplying on the right by U', we get

(AZ+ZC)U' = AZU'+ZCU' = AZU'-ZU'AUU' = A(ZU')-(ZU')A = EU',

so we get the simpler equation:

AW-WA = E", where W=ZU', and E""=EU'.

Taking A = a+0j, as before, W = v+wj, and E" = e"+f"j, we get:

AW-WA=E" => (a+0j)(v+wj)-(v+wj)(a+0j)=(e"+f"j)

         => av + awj -va -wja = e" + f"j

         => av - av + awj - wa'j = e" + f"j

         => 0 + (a-a')wj = e" + f"j

We can now separate the two parts into two equations:

e" = 0

(a-a')w = f" => (alpha + beta i - alpha + beta i) w = f"

             => 2 beta i w = f"

             => 2 beta w = - i f"

             => w = -if"/(2 beta)

Another way to look at this same equation is to rewrite

A = |A| exp(i theta), so that

(AW-WA)/|A| = (|A| exp(i theta) W - |A| W exp(i theta))/|A|

            = exp(i theta) W - W exp(i theta) =

            = exp(i theta) (v+wj) - (v+wj) exp(i theta)

            = exp(i theta) v + exp(i theta) wj - v exp(i theta) - w j exp(i theta)

            = exp(i theta) w j - w j exp(i theta)

            = exp(i theta) w j - exp(-i theta) w j

            = [exp(i theta) - exp(-i theta)] wj

            = 2 i sin(theta) wj = E"/|A| = (e"+f"j)/|A|

Therefore, w = -if"/(2 |A| sin(theta)).

But since a = alpha + beta i = |A| (cos(theta) + i sin(theta)),
beta = |A| sin(theta), so

w = -if"/(2 beta), as before.

References:

Bartels, R.H., and Stewart, G.W.  "Algorithm 432: Solution of the
Matrix Equation AX+XB=C".  Comm. of the ACM 15, 9 (Sept. 1972),
820-826.

Gardiner, J.D., Laub, A.J., Amoto, J.J., and Moler, C.B.  "Solution of
the Sylvester Matrix Equation AXB^T + CXD^T = E".  ACM Trans. Math.
Soft. 18, 2 (June 1992), 223-231.