APPLICATION OF QUATERNIONS TO COM-PUTATION WITH ROTATIONS

Working Paper, Stanford AI Lab, 1979^{1}

by Eugene Salamin

Computer programs which operate on rotations may profitably represent these rotations by quaternions. The orthogonal matrix which performs a rotation by angle θ about axis n is derived. After developing some elementary properties of quaternions, it is shown how the same rotation can be represented by a unit quaternion. Formulas for interconversion between the two representations are derived. For computation with rotations, quaternions offer the advantage of requiring only 4 numbers of storage, compared with 9 numbers for orthogonal matrices. Composition of rotations requires 16 multiplications and 12 additions in quaternion representation, but 27 multiplications and 18 additions in matrix representation. Rotating a vector, with the rotation matrix in hand, requires 9 multiplications and 6 additions. However, if the matrix must be calculated from a quaternion, then this calculation needs 10 multiplications and 21 additions. The quaternion representation is more immune to accumulated computational error. A quaternion^{2} Q which deviates from unicity can be fixed by Q ← Q/ |Q|, however a matrix R which deviates from orthogonality must be fixed by the more involved calculation R ← R(R^{T} R)^{−1/2}, or some approximation thereto. One can mentally picture the rotation group manifold as the 3-sphere with diametrically opposite points identified. The intrinsic metric of the manifold is the same as the metric resulting from the imbedding of the 3sphere in Euclidean 4-space. This makes well-defined the concept of randomly chosen rotations. Finally, it is shown how Lorentz transformations in relativity can be implemented with complex quaternions.

1. Rotations and Orthogonal Matrices

Let a vector x be rotated by angle θ about the unit vector n to obtain the vector x′. Write

x = (x · n) n + (x − (x · n) n).

The first term on the right is the component of x parallel to n; it is unchanged by the rotation. The second term is the component of x perpendicular to n. It is rotated into

(x − (x · n) n) cos θ + n × (x − (x · n) n) sin θ.

1Edited and TeX-formatted by Henry G. Baker (hbaker1@pipeline.com), Nov., 1995. Postings  (http://home.pipeline.com/~hbaker1/quaternion/stanfordaiwp79-salamin.ps.gz and .dvi.gz) by permission of Eugene Salamin.

2Editor’s note: to minimize confusion, we use capital letters for quaternions, lower-case letters for scalars, bold capital letters for matrices, and bold lower-case letters for vectors.

Collecting all terms, we see that

x′ = (cos θ) x + (sin θ) n × x + (1 − cos θ) n (n · x).  (1.1)

We want to write (1.1) in matrix form

x′ = R(θ, n)x,  (1.2)

where R(θ, n) denotes the orthogonal matrix which performs a rotation by angle θ about axis n. For this purpose, we set up a correspondence between vectors and antisymmetric matrices.

n = ( \begin{matrix}n_{1} \\ n_{2} \\ n_{3} \end{matrix} )  ↔  N = ( \begin{matrix}0 & −n_{3} & n_{2} \\ n_{3} & 0 & −n_{1} \\ −n_{2} & n_{1} & 0 \end{matrix} )  (1.3)

This correspondence is useful because Nx = n × x.

The characteristic equation of N is

det(λI − N) = λ^{3} + (n_{1}^{2} + n_{2}^{2} + n_{3}^{2}) λ = 0.

Since n is a unit vector, the eigenvalues are distinct, and so the characteristic equation is also the minimal equation. Thus (Cayley-Hamilton theorem)

N^{3} = −N.  (1.4)

From (1.1) and the definition of N, the rotation matrix for rotation by a small angle ϵ is

R(ϵ, n) = I + ϵN + O(ϵ^{2}).

For a finite rotation,

R(θ, n) = R(θ/k, n)^{k} = (I + (θ/k)N + O(k^{−2}))^{k}.

Now let k → ∞, and we get

R(θ, n) = exp(θN) = e^{θN}.  (1.5)

Expand e^{θN} into a power series. Use (1.4) to reduce all powers of N to N or N^{2}. Finally, identify the power series coefficients of N and N^{2} with their trigonometric sums. The result is

e^{θN} = I + (sin θ) N + (1 − cos θ) N^{2}.  (1.6)

By matrix multiplication,

N^{2} = ( \begin{matrix}−n_{2}^{2} − n_{3}^{2} & n_{1}n_{2} & n_{1}n_{3} \\ n_{2}n_{1} & −n_{1}^{2} − n_{3}^{2} & n_{2}n_{3} \\ n_{3}n_{1} & n_{3}n_{2} & −n_{1}^{2} − n_{2}^{2} \end{matrix} ) = nn − I,  (1.7)

in which nn denotes the dyadic product.

By substituting (1.7) into (1.6), and taking the product x′ = R(θ, n)x, we find agreement with (1.1), which was derived directly from the definition of rotation.

2. Elementary Properties of Quaternions

Quaternions are elements of a certain 4-dimensional algebra, i.e. a vector space endowed with multiplication. For the moment, we make no assumption about the ground field of the vector space.

It is convenient to think of a quaternion, Q = q_{0} + q, as the sum of a scalar, q_{0}, and an ordinary 3-vector, q. Then the quaternion product

R = P Q

is defined by

r_{0} = p_{0}q_{0} − p · q,  (2.1)
r = p_{0} q + p q_{0} + p × q.  (2.2)

We will identify a scalar (field element) with a quaternion having vanishing vector part: s = s + 0. Likewise, we identify a 3-vector with a quaternion having vanishing scalar part: v = 0 + v. If Q = q_{0} + q, then by Q^{∗} is meant the quaternion with the sign of its vector part changed, i.e. Q^{∗} = q_{0} − q.

The quaternion algebra is also given a dot product.^{3}

P · Q = p_{0}q_{0} + p · q.  (2.3)

A unit quaternion is a quaternion Q such that Q · Q = 1. A null quaternion is a quaternion Q such that Q · Q = 0. Non-zero null quaternions exist. Examples are Q = (2, 1, 0, 0) over Z_{5}, and Q = (1, i, 0, 0) over the complex numbers. For this reason, the dot product is, in general, not an inner product.

The following properties are easily established, and are left as an exercise for the reader. Remember that the ground field is arbitrary in the current context.

1. The quaternion 1 is the multiplicative identity.

1Q = Q1 = Q.

2. Quaternion multiplication is associative.

P (QR) = (P Q)R = P QR.

Furthermore,^{4}

P QR = p_{0}q_{0}r_{0} − p_{0} q · r − q_{0} p · r − r_{0} p · q − [p, q, r]
+ q_{0}r_{0} p + p_{0}r_{0} q + p_{0}q_{0} r
+ p_{0} q × r + q_{0} p × r + r_{0} p × q
− (q · r) p + (p · r) q − (p · q) r.

(2.4)

[p, q, r] = p · (q × r) = | \begin{matrix}p_{1} & p_{2} & p_{3} \\ q_{1} & q_{2} & q_{3} \\ r_{1} & r_{2} & r_{3} \end{matrix} | = −Scalar(pqr) = − \frac{pqr + (pqr)^{∗}}{2} .

3Editor’s note: P · Q = Scalar(P^{∗}Q) = Scalar(P Q^{∗}) = \frac{P Q^{∗}+(P Q^{∗})^{∗}}{2} .

4Editor’s note: [p, q, r] is the “scalar triple product” — i.e.,

3. Quaternion multiplication is non-commutative. There exist P and Q such that P Q ≠ QP. However, multiplication by a scalar is always commutative.

4. The ^{∗} operator is an anti-homomorphism.

(P Q)^{∗} = Q^{∗}P^{∗}.

5.

Q · Q = QQ^{∗} = Q^{∗}Q.

6.

(P Q) · (P Q) = (P Q) · (QP) = (P · P)(Q · Q).  (2.5)

7.

(P Q) · (P R) = (P · P)(Q · R).  (2.6)

3. Quaternion Description of Rotations

We now assume we are dealing with real quaternions. The dot product is now an inner product, and we can define the length of a quaternion Q as |Q| = \sqrt{Q · Q}. It will be shown that the rotation matrices are in a 1-2 correspondence with the unit quaternions, i.e. those quaternions, Q, for which |Q| = 1.

Consider the transformation

X′ = QXQ^{∗},  (3.1)

where Q is a unit quaternion. The scalar part, x_{0}, of X is transformed into x′_{0} = (QQ^{∗})x_{0} = x_{0}. From (2.5), we see that X′ · X′ = X · X. Since the scalar part is left unchanged, we also have x′ · x′ = x · x, which is sufficient to prove that the transformation is a rotation. By specialization of (2.4), we see that the vector part, x, of X is transformed into

x′ = (q_{0}^{2} − q · q) x + 2q_{0} q × x + 2q (q · x).  (3.2)

Note the similarity in form between (1.1) and (3.2). It is reminiscent of the trigonometric double angle formulas. By choosing

Q = Q(θ, n) = cos \frac{θ}{2} + sin \frac{θ}{2} n,

(which is not much of a choice, since every unit quaternion is of this form) we make the two equations identical.

This can also be written in a form similar to (1.5). First, note that, as a quaternion product, n^{2} = −1. By manipulation of the power series, just as was done for rotation matrices, we find

Q(θ, n) = e^{(θ/2)n} = cos \frac{θ}{2} + sin \frac{θ}{2} n.  (3.3)

We have now exhibited the correspondence between unit quaternions and proper orthogonal matrices as representations of rotations.

Q(θ, n) ↔ R(θ, n)  (3.4)

This correspondence is a group homomorphism because P (QXQ^{∗})P^{∗} = P QXQ^{∗}P^{∗} = (P Q)X(P Q)^{∗}. It is not quite an isomorphism, since both Q(θ, n) and −Q(θ, n) correspond to the same orthogonal matrix.

Given the elements of a unit quaternion, Q, we can use (3.2) to obtain the elements of the corresponding rotation matrix, R.

R = ( \begin{matrix}q_{0}^{2} + q_{1}^{2} − q_{2}^{2} − q_{3}^{2} & 2(−q_{0}q_{3} + q_{1}q_{2}) & 2(q_{0}q_{2} + q_{1}q_{3}) \\ 2(q_{0}q_{3} + q_{2}q_{1}) & q_{0}^{2} − q_{1}^{2} + q_{2}^{2} − q_{3}^{2} & 2(−q_{0}q_{1} + q_{2}q_{3}) \\ 2(−q_{0}q_{2} + q_{3}q_{1}) & 2(q_{0}q_{1} + q_{3}q_{2}) & q_{0}^{2} − q_{1}^{2} − q_{2}^{2} + q_{3}^{2} \end{matrix} ) . (3.5)

Conversely, given the elements of a proper orthogonal matrix, R, we want to compute the elements of a corresponding unit quaternion, Q. From the diagonal elements of (3.5) we get

q_{0}^{2} = \frac{1}{4}(1 + R_{11} + R_{22} + R_{33}),
q_{1}^{2} = \frac{1}{4}(1 + R_{11} − R_{22} − R_{33}),
q_{2}^{2} = \frac{1}{4}(1 − R_{11} + R_{22} − R_{33}),
q_{3}^{2} = \frac{1}{4}(1 − R_{11} − R_{22} + R_{33}).

(3.6)

From the off-diagonal elements of (3.5) we get

q_{0}q_{1} = \frac{1}{4}(R_{32} − R_{23}),
q_{0}q_{2} = \frac{1}{4}(R_{13} − R_{31}),
q_{0}q_{3} = \frac{1}{4}(R_{21} − R_{12}),
q_{1}q_{2} = \frac{1}{4}(R_{12} + R_{21}),
q_{1}q_{3} = \frac{1}{4}(R_{13} + R_{31}),
q_{2}q_{3} = \frac{1}{4}(R_{23} + R_{32}).

(3.7)

Choose the equation from (3.6) which yields the q_{k} of largest magnitude. In particular, if we find some q_{k}^{2} ≥ 1/4, then no other element of Q can be larger. Since both Q and −Q are solutions to the same R, we can take either sign of the square root. Once we have solved for one q_{k}, the three equations in (3.7) which contain q_{k} on the left hand side can be solved for the remaining elements of Q.

4. Comparison Between Quaternions and Orthogonal Matrices

If rotations are represented by 3 × 3 orthogonal matrices, then it is necessary to store all 9 elements. Composing rotations requires multiplying the two matrices. For each product matrix element, we must do 3 multiplications and 2 additions, for a total of 27 multiplications and 18 additions.

If rotations are represented by quaternions, then only 4 elements need to be stored. Composition of rotations corresponds to quaternion multiplication. When the quaternion product, (2.1) and (2.2) is written out in terms of actual elements, we see that computing each product element requires 4 multiplications and 3 additions, for a total of 16 multiplications and 12 additions.

When rotating a vector, if we already possess the rotation matrix, then only 9 multiplications and 6 additions are needed. If we have to compute the rotation matrix from the quaternion, then we can follow the procedure of (3.5). First compute all the quadratic monomials q_{j}q_{k}, which takes 10 multiplications. Next, each of the three diagonal terms can be calculated with 3 additions. Finally, each of the six off-diagonal terms can be calculated with 2 additions, provided multiplication by 2 is considered as an addition. The total computation in obtaining the matrix is 10 multiplications and 21 additions.

Buildup of computational error will cause the quaternion to become of nonunit length. This can be fixed by dividing the quaternion by |Q|. In contrast, fixing a non-orthogonal matrix is much more complicated (see Appendix A). If A is an almost orthogonal matrix, then A can be re-orthogonalized by replacing it by R = A(A^{T} A)^{−1/2}.

But we can even avoid having to rescale the quaternion provided we rotate vectors according to

x′ = \frac{QxQ^{∗}}{Q · Q} .

No such option exists when orthogonal matrices are used.

5. The Rotation Group Manifold

We already know what the rotation group manifold looks like topologically. It is a 3-sphere, i.e. the 3-dimensional surface of a 4-dimensional ball, with diametrically opposite points identified, since Q and −Q represent the same rotation. This space is usually called a projective 3-space, and denoted by P^{3}. Equivalently, the points of P^{3} can be taken to be the lines through the origin in Euclidean 4-space.

The group manifold inherits a metric, the dot product (2.3), by virtue of its imbedding in Euclidean space. A group manifold also has an intrinsic metric, which is defined by the requirement that it be invariant under group multiplication. If P, Q, and R are group elements, i.e. unit quaternions, then the invariant metric must satisfy (P Q) · (P R) = (QP) · (RP) = Q · R. But, by (2.6), the inherited metric already satisfies this condition. This validates in total the picture of the rotation group manifold as a 3-sphere (with opposite points identified).

The invariant metric implies a measure on the group manifold, and thus allows probability distributions to be defined. This is useful for Monte Carlo simulations where we need to choose “random rotations”. Random rotations are to be chosen uniformly distributed over the 3-sphere.

In order to carry out integrations over the group manifold, we can use hyperspherical coordinates.

q_{0} = cos α,
q_{1} = sin α cos β,
q_{2} = sin α sin β cos γ,
q_{3} = sin α sin β sin γ.

(5.1)

To cover the P^{3}, we only need the region q_{0} ≥ 0. The range of the angles is

\begin{aligned}0 & ≤ α ≤ π/2, \\ 0 & ≤ β ≤ π, \\ 0 & ≤ γ ≤ 2π. \end{aligned}

The length element is

ds^{2} = dα^{2} + sin^{2}α dβ^{2} + sin^{2}α sin^{2}β dγ^{2},  (5.2)

and the volume element is

dΩ = sin^{2}α sin β dα dβ dγ.  (5.3)

The total group volume is

Ω = ∫_{0}^{π/2} dα ∫_{0}^{π} dβ ∫_{0}^{2π} dγ sin^{2}α sin β = π^{2}.  (5.4)

As an application, let’s calculate the mean value of the rotation angle, θ, of randomly chosen rotations. Since q_{0} = cos(θ/2) = cos α, we have

⟨θ⟩ = 2⟨α⟩ = \frac{1}{Ω} ∫ 2α dΩ = \frac{8}{π} ∫_{0}^{π/2} α sin^{2}α dα = \frac{π}{2} + \frac{2}{π} = 126.5 deg.

Appendix A. Re-orthogonaliztion of Matrices

Let A be a matrix which is supposed to be orthogonal, but isn’t, because of accumulated computational error. We want to replace A by the “nearest” orthogonal matrix. Let X be the unknown replacement matrix. X can be found by minimizing the quantity

f = ∑_{i,j} ((X_{ij} − A_{ij})^{2} − M_{ij}(X^{T} X − I)_{ij}) .

The first term on the right hand side assures that X is close to A. In the second term, M is a Lagrange multiplier matrix. We minimize f by variation of X and M. Differentiation of f with respect to the elements of M gives the condition, X^{T} X = I, that X be orthogonal. Since X^{T} X − I is always symmetric, regardless of X, there are only 6 constraints, not 9. So we need only 6 independent Lagrange multipliers, i.e. we assume that M is symmetric.

Then, differentiating f with respect to the elements of X results in

X − A − XM = 0.  (A.1)

The solution of (A.1) proceeds as follows.

X = A(I − M)^{−1}
X^{T} = (I − M)^{−1}A^{T}
X^{T} X = (I − M)^{−1}A^{T} A(I − M)^{−1} = I
A^{T} A = (I − M)^{2}
I − M = (A^{T} A)^{1/2}
X = A(A^{T} A)^{−1/2}

Note that since A is almost orthogonal, A^{T} A is almost the identity, and so exactly one of its square roots is also almost the identity. This is obviously the square root to be used.

Appendix B. Quaternions and Relativity

This appendix will mainly be of interest to the reader acquainted with the special theory of relativity. Quaternions are now complex, and the dot product defined in (2.3) and the definition of unit quaternion are carried over exactly as stated, i.e. there is no complex conjugation. When we need complex conjugation, we denote it by an overline.

We will consider the following two transformations.

U′ = QU\overline{Q}^{∗},  (B.1)
U′ = QUQ^{∗},  (B.2)

in which Q is a complex unit quaternion. (B.1) is the Lorentz transformation of 4-vectors, e.g. space-time displacements, or momentum-energy. (B.2) is the Lorentz transformation of the electromagnetic field.

Starting with (B.1), we see from (2.5) that U′ · U′ = U · U. Suppose the vector part of U is real and the scalar part of U is imaginary, i.e. we have U = it + x, where t and x are real. Also let U′ = it′ + x′. Then the invariance of the dot product is reformulated as

−t^{′2} + x′ · x′ = −t^{2} + x · x,  (B.3)

which is the familiar condition defining the Lorentz transformation.

We proceed to prove that (B.1) is a Lorentz transformation. We must show that t′ and x′ are real, for otherwise (B.3) would not be correctly interpreted. Note that t and x are real if and only if U = it + x satisfies \overline{U} = −U^{∗}, and likewise for U′. So, to prove that (B.1) is a Lorentz transformation, we must show that \overline{U} = −U^{∗} implies \overline{U}′ = −U^{′∗}.

\overline{U}′ = \overline{QU}Q^{∗} = −\overline{Q}U^{∗}Q^{∗} = −(QU\overline{Q}^{∗})^{∗} = −U^{′∗}.

As Q varies over all unit quaternions, the (B.1) transformations can, at most, generate the proper Lorentz transformations, i.e. those which involve no space or time reflections. To show that all the proper Lorentz transformations are obtained, we only need to obtain a generating set. Such a generating set consists of the rotations and the Lorentz boosts, or pure velocity transformations.

The rotations are obtained by choosing Q real. Then Q necessarily has the canonical form (3.3) with θ and n real, and we get

Q (it + x) \overline{Q}^{∗} = Q (it + x) Q^{∗} = it + QxQ^{∗}.

The time part is unchanged, and the space part is rotated, exactly as in section 3.

The pure velocity boosts are obtained by choosing Q to have a real scalar part and imaginary vector part. Then Q necessarily has the canonical form

Q = cosh \frac{α}{2} + i sinh \frac{α}{2} m,  (B.4)

where α is real and m is a real unit vector. We find, by using (2.4),

U′ = QU\overline{Q}^{∗} = QUQ
= (q_{0}^{2} − q · q)u_{0} − 2q_{0} q · u + u − 2q (q · u) + 2q_{0} q u_{0}.

(B.5)

Substitute (B.4) into (B.5), let U = it + x, apply the hyperbolic double angle formulas, and finally, separate the imaginary scalar part from the real vector part.

t′ = (cosh α) t − (sinh α) m · x
x′ = (cosh α) mm · x + (1 − mm) · x − (sinh α) mt

(B.6)

Let’s put (B.6) into its more familiar form. Let m be along the z-axis, let the velocity be v/c = tanh α, and note that cosh α = 1/\sqrt{1 − v^{2}/c^{2}} = γ. Then (B.6) becomes

\begin{aligned}t′ & = γ(t − (v/c)z), \\ x′ & = x, \\ y′ & = y, \\ z′ & = γ(z − (v/c)t). \end{aligned}

In general, (B.6) is the Lorentz transformation for velocity v/c = (tanh α)m.

Next, we examine the transformation (B.2). Note that (B.2) is algebraically identical to (3.1), except that now the quaternions are all complex. Just as in section 3, U′ · U′ = U · U, and u′_{0} = u_{0}, so u′ · u′ = u · u. Anticipating the type of tensor whose transformation is effected by (B.2), we let u = E + iB, where E and B are real.^{5} Then we see that E · E − B · B and E · B are invariants of the transformation.

5Editor’s note: we break our convention here and allow bold capitals to denote the familiar electric and magnetic vectors.

When Q is a pure rotation, and has canonical form (3.3), E and B are independently rotated. When Q is a pure velocity boost, substitute the canonical form (B.4) into (3.2) and apply the double angle formulas to get

u′ = mm · u + (cosh α) (1 − mm) · u + 2i(sinh α) m × u.

With separation into real and imaginary parts, this becomes

E′ = mm · E + (cosh α) (1 − mm) · E − (sinh α) m × B,
B′ = mm · b + (cosh α) (1 − mm) · B + (sinh α) m × E.

(B.7)

To put (B.7) into the familiar form of the Lorentz transformation of the electric and magnetic fields, again choose m to lie along the z-axis.

E′_{x} = γ(E_{x} + (v/c)B_{y}),
B′_{x} = γ(B_{x} − (v/c)E_{y}),
E′_{y} = γ(E_{y} − (v/c)B_{x}),
B′_{y} = γ(B_{y} + (v/c)E_{x}),
E′_{z} = E_{z},
B′_{z} = B_{z}.