We have seen the existence of several decompositions that apply in any dimension, namely independent planes, sequential angles, and nested dimensions. In all these cases we can either decompose a matrix or construct one. We have also given special attention to 3 × 3 rotation matrices, and these warrant further attention, in both directions (Stuelpnagel 1964).

Quaternion

Given the unit quaternion $q = w + x i + y j + z k$ , the equivalent pre-multiplied (to be used with column vectors) 3 × 3 rotation matrix is ^[19]

Q={\begin{bmatrix}1-2y^{2}-2z^{2}&2xy-2zw&2xz+2yw\\2xy+2zw&1-2x^{2}-2z^{2}&2yz-2xw\\2xz-2yw&2yz+2xw&1-2x^{2}-2y^{2}\end{bmatrix}}.

Now every quaternion component appears multiplied by two in a term of degree two, and if all such terms are zero what is left is an identity matrix. This leads to an efficient, robust conversion from any quaternion – whether unit or non-unit – to a 3 × 3 rotation matrix. Given:

{\begin{aligned}n&=w\times w+x\times x+y\times y+z\times z\\s&={\begin{cases}0&{\text{if }}n=0\\{\frac {2}{n}}&{\text{otherwise}}\end{cases}}\\\end{aligned}}

we can calculate

Q={\begin{bmatrix}1-s(yy+zz)&s(xy-wz)&s(xz+wy)\\s(xy+wz)&1-s(xx+zz)&s(yz-wx)\\s(xz-wy)&s(yz+wx)&1-s(xx+yy)\end{bmatrix}}

Freed from the demand for a unit quaternion, we find that nonzero quaternions act as homogeneous coordinates for 3 × 3 rotation matrices. The Cayley transform, discussed earlier, is obtained by scaling the quaternion so that its $w$ component is 1. For a 180° rotation around any axis, $w$ will be zero, which explains the Cayley limitation.

The sum of the entries along the main diagonal (the trace), plus one, equals $4 - 4(x 2 + y 2 + z 2)$ , which is $4 w 2$ . Thus we can write the trace itself as $2 w 2 + 2 w 2 - 1$ ; and from the previous version of the matrix we see that the diagonal entries themselves have the same form: $2 x 2 + 2 w 2 - 1$ , $2 y 2 + 2 w 2 - 1$ , and $2 z 2 + 2 w 2 - 1$ . So we can easily compare the magnitudes of all four quaternion components using the matrix diagonal. We can, in fact, obtain all four magnitudes using sums and square roots, and choose consistent signs using the skew-symmetric part of the off-diagonal entries:

{\begin{aligned}t&=\operatorname {tr} Q=Q_{xx}+Q_{yy}+Q_{zz}\quad ({\text{the trace of }}Q)\\r&={\sqrt {1+t}}\\w&={\tfrac {1}{2}}r\\x&=\operatorname {sgn} \left(Q_{zy}-Q_{yz}\right)\left|{\tfrac {1}{2}}{\sqrt {1+Q_{xx}-Q_{yy}-Q_{zz}}}\right|\\y&=\operatorname {sgn} \left(Q_{xz}-Q_{zx}\right)\left|{\tfrac {1}{2}}{\sqrt {1-Q_{xx}+Q_{yy}-Q_{zz}}}\right|\\z&=\operatorname {sgn} \left(Q_{yx}-Q_{xy}\right)\left|{\tfrac {1}{2}}{\sqrt {1-Q_{xx}-Q_{yy}+Q_{zz}}}\right|\end{aligned}}

Alternatively, use a single square root and division

{\begin{aligned}t&=\operatorname {tr} Q=Q_{xx}+Q_{yy}+Q_{zz}\\r&={\sqrt {1+t}}\\s&={\tfrac {1}{2r}}\\w&={\tfrac {1}{2}}r\\x&=\left(Q_{zy}-Q_{yz}\right)s\\y&=\left(Q_{xz}-Q_{zx}\right)s\\z&=\left(Q_{yx}-Q_{xy}\right)s\end{aligned}}

This is numerically stable so long as the trace, $t$ , is not negative; otherwise, we risk dividing by (nearly) zero. In that case, suppose $Q xx$ is the largest diagonal entry, so $x$ will have the largest magnitude (the other cases are derived by cyclic permutation); then the following is safe.

{\begin{aligned}r&={\sqrt {1+Q_{xx}-Q_{yy}-Q_{zz}}}\\s&={\tfrac {1}{2r}}\\w&=\left(Q_{zy}-Q_{yz}\right)s\\x&={\tfrac {1}{2}}r\\y&=\left(Q_{xy}+Q_{yx}\right)s\\z&=\left(Q_{zx}+Q_{xz}\right)s\end{aligned}}

If the matrix contains significant error, such as accumulated numerical error, we may construct a symmetric 4 × 4 matrix,

K={\frac {1}{3}}{\begin{bmatrix}Q_{xx}-Q_{yy}-Q_{zz}&Q_{yx}+Q_{xy}&Q_{zx}+Q_{xz}&Q_{zy}-Q_{yz}\\Q_{yx}+Q_{xy}&Q_{yy}-Q_{xx}-Q_{zz}&Q_{zy}+Q_{yz}&Q_{xz}-Q_{zx}\\Q_{zx}+Q_{xz}&Q_{zy}+Q_{yz}&Q_{zz}-Q_{xx}-Q_{yy}&Q_{yx}-Q_{xy}\\Q_{zy}-Q_{yz}&Q_{xz}-Q_{zx}&Q_{yx}-Q_{xy}&Q_{xx}+Q_{yy}+Q_{zz}\end{bmatrix}},

and find the eigenvector, $(x, y, z, w)$ , of its largest magnitude eigenvalue. (If $Q$ is truly a rotation matrix, that value will be 1.) The quaternion so obtained will correspond to the rotation matrix closest to the given matrix (Bar-Itzhack 2000) (Note: formulation of the cited article is post-multiplied, works with row vectors).

Polar decomposition

If the $n \times n$ matrix $M$ is nonsingular, its columns are linearly independent vectors; thus the Gram–Schmidt process can adjust them to be an orthonormal basis. Stated in terms of numerical linear algebra, we convert $M$ to an orthogonal matrix, $Q$ , using QR decomposition. However, we often prefer a $Q$ closest to $M$ , which this method does not accomplish. For that, the tool we want is the polar decomposition (Fan & Hoffman 1955; Higham 1989).

To measure closeness, we may use any matrix norm invariant under orthogonal transformations. A convenient choice is the Frobenius norm, $‖ Q - M ‖ F$ , squared, which is the sum of the squares of the element differences. Writing this in terms of the trace, $Tr$ , our goal is,

Find

Q

minimizing

Tr( (Q - M) T (Q - M) )

, subject to

Q T Q = I

.

Though written in matrix terms, the objective function is just a quadratic polynomial. We can minimize it in the usual way, by finding where its derivative is zero. For a 3 × 3 matrix, the orthogonality constraint implies six scalar equalities that the entries of $Q$ must satisfy. To incorporate the constraint(s), we may employ a standard technique, Lagrange multipliers, assembled as a symmetric matrix, $Y$ . Thus our method is:

Differentiate

Tr( (Q - M) T (Q - M) + (Q T Q - I) Y)

with respect to (the entries of)

Q

, and equate to zero.

Consider a 2 × 2 example. Including constraints, we seek to minimize

{\begin{aligned}&\left(Q_{xx}-M_{xx}\right)^{2}+\left(Q_{xy}-M_{xy}\right)^{2}+\left(Q_{yx}-M_{yx}\right)^{2}+\left(Q_{yy}-M_{yy}\right)^{2}\\&\quad {}+\left(Q_{xx}^{2}+Q_{yx}^{2}-1\right)Y_{xx}+\left(Q_{xy}^{2}+Q_{yy}^{2}-1\right)Y_{yy}+2\left(Q_{xx}Q_{xy}+Q_{yx}Q_{yy}\right)Y_{xy}.\end{aligned}}

Taking the derivative with respect to $Q xx$ , $Q xy$ , $Q yx$ , $Q yy$ in turn, we assemble a matrix.

2{\begin{bmatrix}Q_{xx}-M_{xx}+Q_{xx}Y_{xx}+Q_{xy}Y_{xy}&Q_{xy}-M_{xy}+Q_{xx}Y_{xy}+Q_{xy}Y_{yy}\\Q_{yx}-M_{yx}+Q_{yx}Y_{xx}+Q_{yy}Y_{xy}&Q_{yy}-M_{yy}+Q_{yx}Y_{xy}+Q_{yy}Y_{yy}\end{bmatrix}}

In general, we obtain the equation

0=2(Q-M)+2QY,

so that

M=Q(I+Y)=QS,

where $Q$ is orthogonal and $S$ is symmetric. To ensure a minimum, the $Y$ matrix (and hence $S$ ) must be positive definite. Linear algebra calls $QS$ the polar decomposition of $M$ , with $S$ the positive square root of $S 2 = M T M$ .

S^{2}=\left(Q^{\mathsf {T}}M\right)^{\mathsf {T}}\left(Q^{\mathsf {T}}M\right)=M^{\mathsf {T}}QQ^{\mathsf {T}}M=M^{\mathsf {T}}M

When $M$ is non-singular, the $Q$ and $S$ factors of the polar decomposition are uniquely determined. However, the determinant of $S$ is positive because $S$ is positive definite, so $Q$ inherits the sign of the determinant of $M$ . That is, $Q$ is only guaranteed to be orthogonal, not a rotation matrix. This is unavoidable; an $M$ with negative determinant has no uniquely defined closest rotation matrix.

Axis and angle

To efficiently construct a rotation matrix $Q$ from an angle $θ$ and a unit axis $u$ , we can take advantage of symmetry and skew-symmetry within the entries. If $x$ , $y$ , and $z$ are the components of the unit vector representing the axis, and

{\begin{aligned}c&=\cos \theta \\s&=\sin \theta \\C&=1-c\end{aligned}}

then

Q(\theta )={\begin{bmatrix}xxC+c&xyC-zs&xzC+ys\\yxC+zs&yyC+c&yzC-xs\\zxC-ys&zyC+xs&zzC+c\end{bmatrix}}

Determining an axis and angle, like determining a quaternion, is only possible up to the sign; that is, $(u, θ)$ and $(- u, - θ)$ correspond to the same rotation matrix, just like $q$ and $- q$ . Additionally, axis–angle extraction presents additional difficulties. The angle can be restricted to be from 0° to 180°, but angles are formally ambiguous by multiples of 360°. When the angle is zero, the axis is undefined. When the angle is 180°, the matrix becomes symmetric, which has implications in extracting the axis. Near multiples of 180°, care is needed to avoid numerical problems: in extracting the angle, a two-argument arctangent with $atan2 (sin θ, cos θ)$ equal to $θ$ avoids the insensitivity of arccos; and in computing the axis magnitude in order to force unit magnitude, a brute-force approach can lose accuracy through underflow (Moler & Morrison 1983).

A partial approach is as follows:

{\begin{aligned}x&=Q_{zy}-Q_{yz}\\y&=Q_{xz}-Q_{zx}\\z&=Q_{yx}-Q_{xy}\\r&={\sqrt {x^{2}+y^{2}+z^{2}}}\\t&=Q_{xx}+Q_{yy}+Q_{zz}\\\theta &=\operatorname {atan2} (r,t-1)\end{aligned}}

The $x$ -, $y$ -, and $z$ -components of the axis would then be divided by $r$ . A fully robust approach will use a different algorithm when $t$ , the trace of the matrix $Q$ , is negative, as with quaternion extraction. When $r$ is zero because the angle is zero, an axis must be provided from some source other than the matrix.

Euler angles

Complexity of conversion escalates with Euler angles (used here in the broad sense). The first difficulty is to establish which of the twenty-four variations of Cartesian axis order we will use. Suppose the three angles are $θ 1$ , $θ 2$ , $θ 3$ ; physics and chemistry may interpret these as

Q(\theta _{1},\theta _{2},\theta _{3})=Q_{\mathbf {z} }(\theta _{1})Q_{\mathbf {y} }(\theta _{2})Q_{\mathbf {z} }(\theta _{3}),

while aircraft dynamics may use

Q(\theta _{1},\theta _{2},\theta _{3})=Q_{\mathbf {z} }(\theta _{3})Q_{\mathbf {y} }(\theta _{2})Q_{\mathbf {x} }(\theta _{1}).

One systematic approach begins with choosing the rightmost axis. Among all permutations of $(x, y, z)$ , only two place that axis first; one is an even permutation and the other odd. Choosing parity thus establishes the middle axis. That leaves two choices for the left-most axis, either duplicating the first or not. These three choices gives us 3 × 2 × 2 = 12 variations; we double that to 24 by choosing static or rotating axes.

This is enough to construct a matrix from angles, but triples differing in many ways can give the same rotation matrix. For example, suppose we use the $zyz$ convention above; then we have the following equivalent pairs:

(90°,	45°,	−105°)	≡	(−270°,	−315°,	255°)	multiples of 360°
(72°,	0°,	0°)	≡	(40°,	0°,	32°)	singular alignment
(45°,	60°,	−30°)	≡	(−135°,	−60°,	150°)	bistable flip

Angles for any order can be found using a concise common routine (Herter & Lott 1993; Shoemake 1994).

The problem of singular alignment, the mathematical analog of physical gimbal lock, occurs when the middle rotation aligns the axes of the first and last rotations. It afflicts every axis order at either even or odd multiples of 90°. These singularities are not characteristic of the rotation matrix as such, and only occur with the usage of Euler angles.

The singularities are avoided when considering and manipulating the rotation matrix as orthonormal row vectors (in 3D applications often named the right-vector, up-vector and out-vector) instead of as angles. The singularities are also avoided when working with quaternions.

Vector to vector formulation

In some instances it is interesting to describe a rotation by specifying how a vector is mapped into another through the shortest path (smallest angle). In $\mathbb {R} ^{3}$ this completely describes the associated rotation matrix. In general, given $\mathbb {S}$ , the matrix

R:=I+yx^{\mathsf {T}}-xy^{\mathsf {T}}+{\frac {1}{1+\langle x,y\rangle }}\left(yx^{\mathsf {T}}-xy^{\mathsf {T}}\right)^{2}

belongs to $SO(n + 1)$ and maps $x$ to $y$ .^[20]