If you have ever played with the RotationTransform command in Mathematica, you may have at some point become impressed at how quickly it returns answers in an arbitrary number of dimensions. After all, the details section of the documentation for this command says explicitly that it “can effectively specify any element of the n-dimensional rotation group SO(n).” Since elements of the group can be evaluated by exponentiating the generator of the element, in this case an orthogonal matrix, it appears at first sight that Mathematica knows how to exponentiate an orthogonal matrix in an arbitrary number of dimensions, and very quickly at that.

It turns out that there is a simpler way to evaluate arbitrary elements of the group than by explicitly exponentiating each individual matrix. To show how this works, first quickly review the relationship between a rotation matrix and its generator. For the simplest case of rotation in a two-dimensional plane, the rotated vector is related to the initial vector by

$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}cosa& -sina\\ sina& cosa\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$

for a right-hand rotation through the angle a. The generator of this rotation is represented by the matrix $\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)$ . The square of this matrix is the negative of the identity matrix, its cube is its own negative, and the full exponentiation is

$\begin{array}{l}exp\left[a\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\right]=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+a\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)-\frac{a^2}{2!}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)-\frac{a^3}{3!}\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)+\cdots \\ \phantom{exp\left[a\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\right]}=cosa\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+sina\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\\ exp\left[a\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\right]=\left(\begin{array}{cc}cosa& -sina\\ sina& cosa\end{array}\right)\end{array}$

leading to the rotation matrix given above. The generator of the rotation can be written as an outer product of the two unit vectors along the x-axis and the y-axis,

$\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{cc}1& 0\end{array}\right)-\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{cc}0& 1\end{array}\right)=\hat{\mathbf{y}}{\hat{\mathbf{x}}}^{\mathrm{T}}-\hat{\mathbf{x}}{\hat{\mathbf{y}}}^{\mathrm{T}}$

where the transposed vector on the right of each term is multiplied leftward onto each component of the vector to produce matrices that are added together. Since the final result is now a general vector statement, it can be extended immediately to higher dimensions. The equivalent three-dimensional generator is

$L_{12}=\hat{\mathbf{y}}{\hat{\mathbf{x}}}^{\mathrm{T}}-\hat{\mathbf{x}}{\hat{\mathbf{y}}}^{\mathrm{T}}=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\end{array}\right)-\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)\left(\begin{array}{ccc}0& 1& 0\end{array}\right)=\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)$

where the subscript indicates that this is a right-hand rotation from the x-axis towards the y-axis. Note that in forming this generator, the first index appears with a negative sign as a non-transposed vector. The two remaining generators of right-hand rotations are

$\begin{array}{l}{L}_{23}=\hat{\mathbf{z}}{\hat{\mathbf{y}}}^{\mathrm{T}}-\hat{\mathbf{y}}{\hat{\mathbf{z}}}^{\mathrm{T}}=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\left(\begin{array}{ccc}0& 1& 0\end{array}\right)-\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)\left(\begin{array}{ccc}0& 0& 1\end{array}\right)=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -1\\ 0& 1& 0\end{array}\right)\\ L_{31}=\hat{\mathbf{x}}{\hat{\mathbf{z}}}^{\mathrm{T}}-\hat{\mathbf{z}}{\hat{\mathbf{x}}}^{\mathrm{T}}=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)\left(\begin{array}{ccc}0& 0& 1\end{array}\right)-\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\end{array}\right)=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right)\end{array}$

Because there are three coordinate axes, the number of pairwise combinations is also three, so that in three dimensions one can make the identifications

$L_1=L_{23}{L}_2=L_{31}{L}_3=L_{12}$

and describe rotations as being about the corresponding axes. The rotation matrices can be evaluated separately about each axis just as for the rotation in the two-dimensional plane, with the results

$\begin{array}{l}exp\left[a_1{L}_1\right]=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos{a}_1& -sin{a}_1\\ 0& sin{a}_1& cos{a}_1\end{array}\right)\\ exp\left[a_2{L}_2\right]=\left(\begin{array}{ccc}cos{a}_2& 0& sin{a}_2\\ 0& 1& 0\\ -sin{a}_2& 0& cos{a}_2\end{array}\right)\\ exp\left[a_3{L}_3\right]=\left(\begin{array}{ccc}cos{a}_3& -sin{a}_3& 0\\ sin{a}_3& cos{a}_3& 0\\ 0& 0& 1\end{array}\right)\end{array}$

In three dimensions one can also consider a rotation about a general axis by a single angle $a=\sqrt{a_1^2+a_2^2+a_3^2}$ . The generator of this rotation is the linear combination $\mathbf{a}·\mathbf{L}=a\mathbf{n}·\mathbf{L}$ , where n is a unit vector in the direction of the axis of rotation. The rotation matrix for this general rotation is found by exponentiating the matrix

$A=\left(\begin{array}{ccc}0& -a_3& a_2\\ a_3& 0& -a_1\\ -a_2& a_1& 0\end{array}\right)=a\left(\begin{array}{ccc}0& -n_3& n_2\\ n_3& 0& -n_1\\ -n_2& n_1& 0\end{array}\right)$

where $n_1^2+n_2^2+n_3^2=1$ because these numbers are direction cosines giving the inclination of the axis of rotation to the three coordinate axes. The square and cube of this matrix are

$\begin{array}{l}{A}^2=-a^{2}I+\left(\begin{array}{ccc}{a}_1^2& a_1{a}_2& a_1{a}_3\\ a_1{a}_2& a_2^2& a_2{a}_3\\ a_1{a}_3& a_2{a}_3& a_3^2\end{array}\right)\\ \\ A^3=-(a_1^2+a_2^2+a_3^2)A=-a^{2}A\end{array}$

and the exponentiated matrix can be written as

$e^A=I+\left(\frac{sina}{a}\right)A+\left(\frac{1-cosa}{a^2}\right)A^2$

where the denominators normalize powers of the coefficient matrix. When applied to a vector, the result is Rodrigues’ rotation formula:

$\begin{array}{l}{e}^A\mathbf{r}=\mathbf{r}cosa+sina\left(\begin{array}{ccc}0& -n_3& n_2\\ n_3& 0& -n_1\\ -n_2& n_1& 0\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)\\ +(1-cosa)\left(\begin{array}{ccc}{n}_1^2& n_1{n}_2& n_1{n}_3\\ n_1{n}_2& n_2^2& n_2{n}_3\\ n_1{n}_3& n_2{n}_3& n_3^2\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)\\ e^A\mathbf{r}=\mathbf{r}cosa+(\mathbf{n}\times \mathbf{r})sina+\mathbf{n}(\mathbf{n}·\mathbf{r})(1-cosa)\end{array}$

In four dimensions one can no longer describe a rotation as being about an axis, because there are two axes perpendicular to every plane and there would be an ambiguity in the description. In higher dimensions there are even more axes perpendicular to each plane, so a rotation in n dimensions is best described as being in an (n−1)-dimensional hyperplane in a direction from one unit vector towards another.

The definitions above of generators in terms of outer products are in vector notation, and so can be extended immediately to describe rotation in any hyperplane defined by two n-dimensional vectors. Given any two orthogonal unit vectors n₁ and n₂, which means

${\mathbf{n}}_1^{\mathrm{T}}{\mathbf{n}}_1=1{\mathbf{n}}_1^{\mathrm{T}}{\mathbf{n}}_2={\mathbf{n}}_2^{\mathrm{T}}{\mathbf{n}}_1=0{\mathbf{n}}_2^{\mathrm{T}}{\mathbf{n}}_2=1$

the generator of rotations in the hyperplane spanned by the two vectors is

$L_{{\mathbf{n}}_1{\mathbf{n}}_2}={\mathbf{n}}_2{\mathbf{n}}_1^{\mathrm{T}}-{\mathbf{n}}_1{\mathbf{n}}_2^{\mathrm{T}}$

Forming powers of this generator,

$\begin{array}{l}{L}_{{\mathbf{n}}_1{\mathbf{n}}_2}^2=-({\mathbf{n}}_1{\mathbf{n}}_1^{\mathrm{T}}+{\mathbf{n}}_2{\mathbf{n}}_2^{\mathrm{T}})\\ L_{{\mathbf{n}}_1{\mathbf{n}}_2}^3=-({\mathbf{n}}_2{\mathbf{n}}_1^{\mathrm{T}}-{\mathbf{n}}_1{\mathbf{n}}_2^{\mathrm{T}})\end{array}$

the n-dimensional rotation matrix is simply

$exp\left[a{L}_{{\mathbf{n}}_1{\mathbf{n}}_2}\right]=I+({\mathbf{n}}_2{\mathbf{n}}_1^{\mathrm{T}}-{\mathbf{n}}_1{\mathbf{n}}_2^{\mathrm{T}})sina+({\mathbf{n}}_1{\mathbf{n}}_1^{\mathrm{T}}+{\mathbf{n}}_2{\mathbf{n}}_2^{\mathrm{T}})(cosa-1)$

The difference of sign in the last term compared to Rodrigues’ rotation formula is due to an extra minus sign from the neighboring outer product. This n×n matrix can now be applied to an n-dimensional vector to find its final value after rotation in the specified hyperplane.

Given two vectors that are not orthogonal, one can apply the Gram-Schmidt orthogonality process, renormalize the second vector and form the rotation matrix with the two newly orthogonal unit vectors. And that is precisely what Mathematica does in producing answers so quickly in any number of dimensions. The exponentiation has already been carried out in the general formula and does not need to be done explicitly for every particular case, and forming the n×n matrix is trivial for Mathematica.

To demonstrate the consistency of this general formula with the result above in three dimensions, one need merely choose any vector satisfying v₁ · n = 0 with respect to the three-dimensional axis of rotation, form its orthogonal vector v₂ = v₁×n , normalize both vectors and do a wee bit of algebra. Since there is an entire plane of rotation in which to choose the vectors, any initial choice leads by this method to the same result.

In four dimensions, exponentiation of a general orthogonal matrix is capable of explicit evaluation using a method based on the local isomorphism SO(4) $\sim$ SO(3)$\otimes$SO(3) . The general coefficient matrix is

$A=\left(\begin{array}{cccc}0& a_{12}& a_{13}& a_{14}\\ -a_{12}& 0& a_{23}& a_{24}\\ -a_{13}& -a_{23}& 0& a_{34}\\ -a_{14}& -a_{24}& -a_{34}& 0\end{array}\right)$

where the minus signs have all been kept below the diagonal for convenience. Squaring this matrix itself produces a symmetric matrix whose interpretation is not immediately clear, but if one first writes

$\begin{array}{l}A=\frac{1}{2}(A_{+}+A_{-})\\ A_{+}=\left(\begin{array}{cccc}0& a_{12}+a_{34}& a_{13}-a_{24}& a_{14}+a_{23}\\ -a_{12}-a_{34}& 0& a_{14}+a_{23}& -a_{13}+a_{24}\\ -a_{13}+a_{24}& -a_{14}-a_{23}& 0& a_{12}+a_{34}\\ -a_{14}-a_{23}& a_{13}-a_{24}& -a_{12}-a_{34}& 0\end{array}\right)\\ A_{-}=\left(\begin{array}{cccc}0& a_{12}-a_{34}& a_{13}+a_{24}& -a_{14}+a_{23}\\ -a_{12}+a_{34}& 0& a_{14}-a_{23}& a_{13}+a_{24}\\ -a_{13}-a_{24}& -a_{14}+a_{23}& 0& -a_{12}+a_{34}\\ a_{14}-a_{23}& -a_{13}-a_{24}& a_{12}-a_{34}& 0\end{array}\right)\end{array}$

then it is straightforward to verify that these two matrices commute: $[A_{+},A_{-}]=0$ . Forming their squares gives

$\begin{array}{l}{A}_{+}^2=-\left[\right(a_{12}+a_{34}{)}^2+(a_{13}-a_{24}{)}^2+(a_{14}+a_{23}{)}^2]I=-a_{+}^{2}I\\ A_{-}^2=-\left[\right(a_{12}-a_{34}{)}^2+(a_{13}+a_{24}{)}^2+(a_{14}-a_{23}{)}^2]I=-a_{-}^{2}I\end{array}$

so that the complete exponentiated matrix is

$\begin{array}{l}{e}^A=e^{A_{+}/2}{e}^{A_{-}/2}=e^{A_{-}/2}{e}^{A_{+}/2}\\ \phantom{e^A}=[cos(a_{+}/2)I+\frac{sin(a_{+}/2)}{a_{+}}{A}_{+}]\\ \times [cos(a_{-}/2)I+\frac{sin(a_{-}/2)}{a_{-}}{A}_{-}]\end{array}$

One can describe this method as arranging the six parameters of SO(4) into two vectors, remembering that one is working with four-dimensional representations of SO(3).

Uploaded 2013.01.15 — Updated 2017.12.17 analyticphysics.com