Closed form solution to rotation in arbitrary many planes in arbitrary dimensions. In each coordinate space $V$ with dimension $\dim(V)$, we can describe any rotation operator $R : V \to V$ as a product of rotations in as few as $\dim(V) - 1$ orthogonal planes in the space.
Let's say that given the canonical basis $\mathcal B$ of $V$, we take the canonical ordering of the basis vectors $\vec b_i$ and consider only rotations in the planes $\mathcal P = \{\ \mathcal P_{i, i+1} = \text{span} (b_i,b_{i+1})\ |\ b_i, b_{i+1} \in \mathcal B\ \}$. Is there a closed form solution for the matrix of the final rotation operator?
Consider an example in three dimensions. Let $R_{i, j}(\theta_k)$ be the rotation by $\theta_k$ in the plane $\mathcal P_{i, j}$.
\begin{align}
R(\theta_1, \theta_2) = R_{1,2}(\theta_1)R_{2,3}(\theta_2) =& \left[\begin{array}{ccc} \cos(\theta_1) & -\sin(\theta_1) & 0 \\ \sin(\theta_1) & \cos(\theta_1) & 0 \\ 0 & 0 & 1 \end{array}\right]
\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos(\theta_2) & -\sin(\theta_2) \\ 0 & \sin(\theta_2) & \cos(\theta_2) \end{array}\right] \\
=& \left[\begin{array}{ccc} \cos(\theta_1) & -\sin(\theta_1)\cos(\theta_2) & \sin(\theta_1)\sin(\theta_2) \\ \sin(\theta_1) & \cos(\theta_1)\cos(\theta_2) & -\cos(\theta_1)\sin(\theta_2) \\ 0 & \sin(\theta_2) & \cos(\theta_2) \end{array}\right]
\end{align}
I can see the beginnings of a pattern, but I feel as though the pattern will break shortly after increasing the number of dimensions.
Is there already a well-studied closed form for this? Is it a fool's errand?
 A: Maybe this ?


It seems that you get a Hessenberg matrix with all sines on the subdiagonal, cosine pairs on the main diagonal, superdiagonals with alternating signs and products on increasing lengths, not so easy to describe.

Update:
With better alignment and pseudo-factors $c_{\bar1},c_4$, the pattern becomes much clearer.
The subdiagonal has all sines. The upper triangle has pairs of cosines (extended row and column indexes), and suffix products of the sines. 
$$\left(\begin{matrix}
&\color{green}{c_{\bar1}}\color{blue}{c_0}&-\color{green}{c_{\bar1}}\color{blue}{c_1}s_0&\ \ \ \color{green}{c_{\bar1}}\color{blue}{c_2}s_0s_1&-\color{green}{c_{\bar1}}\color{blue}{c_3}s_0s_1s_2&\ \ \ \color{green}{c_{\bar1}}\color{blue}{c_4}s_0s_1s_2s_3\\
&\color{magenta}{s_0}&\ \ \ \color{green}{c_0}\color{blue}{c_1}\ \ \ &-\color{green}{c_0}\color{blue}{c_2}\ \ \ \ s_1&\ \ \ \color{green}{c_0}\color{blue}{c_3}\ \ \ \ s_1s_2&-\color{green}{c_0}\color{blue}{c_4}\ \ \ \ s_1s_2s_3\\
&0&\color{magenta}{s_1}&\ \ \ \color{green}{c_1}\color{blue}{c_2}\ \ \ \ \ \ \ \ &-\color{green}{c_1}\color{blue}{c_3}\ \ \ \ \ \ \ \ s_2&\ \ \ \color{green}{c_1}\color{blue}{c_4}\ \ \ \ \ \ \ s_2s_3\\
&0&0&\color{magenta}{s_2}&\ \ \ \color{green}{c_2}\color{blue}{c_3}\ \ \ \ \ \ \ \ \ \ \ &-\color{green}{c_2}\color{blue}{c_4}\ \ \ \ \ \ \ \ \ \ \ s_3\\
&0&0&0&\color{magenta}{s_3}&\color{green}{c_3}\color{blue}{c_4}\ \ \ \ \ \ \ \ \ \ \ \\
\end{matrix}\right)$$
For the upper triangular elements,
$$R_{rc}=(-1)^{r+c}\color{green}{\cos(\theta_{r-1})}\color{blue}{\cos(\theta_{c})}\prod_{k=r}^{c-1}\sin(\theta_k).$$
By computing the products incrementally, the whole matrix can be obtained in $\approx\dfrac32d^2$ multiplies.
A: To represent rotation matrices nicely, we can use the Jordan normal form.
A rotation operator can be represented by an orthogonal matrix (a matrix $R$ satisfying $R^T R = I = R R^T$).  For these matrices, we have the result that their (complex) eigenvalues satisfy $|\lambda|^2 = 1$ (you can prove this by calculating $||\lambda v||^2$ for some eigenvector $v$).  This means that either the eigenvalues are $1$ or $-1$, or they are complex.  For complex eigenvalues we can write $\lambda_i = e^{i\theta_i} = \cos\theta_i + i \sin\theta_i$.
We can then use the real Jordan form to see that $R$ is similar to:
$$
RS = SJ\\[1cm]
J = 
\left[
\begin{array}{ccccc|ccc|ccc}
\cos\theta_1 & \sin\theta_1 \\
-\sin\theta_1 & \cos\theta_1 \\
& & \ddots \\
& & & \cos \theta_m & \sin\theta_m\\
& & & -\sin\theta_m & \cos\theta_m\\\hline
&&&&& -1\\
&&&&&& \ddots\\
&&&&&&& -1\\\hline
&&&&&&&& 1\\
&&&&&&&&& \ddots\\
&&&&&&&&&& 1
\end{array}
\right]
$$
(The fact that no identity matrices occur on the upper diagonal is due to the fact that an orthogonal matrix can be diagonalized.)  Here $2m$ is the number of non-real eigenvalues (note that these always occur in complex conjugate pairs, which allows this representation).  The matrix that transforms $R$ into this form has the real and imaginary parts of the eigenvectors of $R$ as its columns (in the same order as the eigenvalues in $J$).  For a proper rotation, the number of $-1$ eigenvalues will be even.
Hence, we can see that a general rotation operator can be represented as $m$ rotations that influence only a 2-dimensional space (spanned by the real and complex part of an eigenvector corresponding to the complex eigenvalue), and an even number of reflections in some hyperplanes (orthogonal to the eigenvector corresponding to the $-1$ eigenvalues).
Clarification
In response to the comments I would like give some additional clarifications to my answer, hoping to shed light on the difference between the two Jordan forms and the interpretation of the result.
The Jordan Normal Form
Let $M \in \mathbb{R}^{n \times n}$ be a matrix with real entries.  In general it will then have $m$ distinct eigenvalues $\lambda_i \in \mathbb{C}$.  The matrix might have enough eigenvectors to form a basis $S = [\vec{v_1^{(1)}},\dots,\vec{v_{k_1}^{(1)}},\dots,\vec{v_1^{(m)}},\dots,\vec{v_{k_m}^{(m)}}]$.  In this case you can diagonalize the matrix by $D = S^{-1} M S$.  However, in general this is not possible.  In this case you have to build a number of Jordan chains for each eigenvalue.
We let $k_i$ denote the number of Jordan chains corresponding to $\lambda_i$ and $l_{ij}$ the length of the $j$th Jordan chain of $\lambda_i$.  We then have a basis of generalized eigenvectors (for $j = 1,\dots,k_i$):
$$
(A - \lambda_i I) \vec{v_{j, p}^{(i)}} = \vec{v_{j,p-1}^{(i)}}, \qquad
p = l_{ij},\dots,2\\
(A - \lambda_i I) \vec{v_{j, 1}^{(i)}} = 0
$$
We can now arrange these in a matrix:
$$
S = [
\vec{v_{1,1}^{(1)}},\dots,\vec{v_{1,l_{11}}^{(1)}},\dots,
\vec{v_{k_1,1}^{(1)}},\dots,\vec{v_{k_1,l_{k_1 1}}^{(1)}},\dots,
\vec{v_{1,1}^{(m)}},\dots,\vec{v_{1,l_{1 m}}^{(m)}},\dots,
\vec{v_{k_m,1}^{(m)}},\dots,\vec{v_{k_m,l_{m k_m}}^{(m)}}
]
$$
We then see that $M$ takes on the following form:
$$
J = S^{-1} M S =
\left[\begin{array}{}
J_1 \\
& J_2\\
&& \ddots\\
&&& J_m
\end{array}\right],
\quad\text{where:}
\\
J_i = \left[\begin{array}{}
J_{i1}\\
& \ddots\\
&& J_{ik_i}
\end{array}\right]
\quad\text{and}\quad
J_{ij} = \left[\begin{array}{}
\lambda_i & 1\\
& \lambda_i & \ddots \\
& & \ddots & 1 \\
& & & \lambda_i
\end{array}\right]
\in \mathbb{C}^{l_{ij} \times l_{ij}}.
$$
This is then finally called the Jordan normal form.
Real Jordan form
We now make the powerful observation that if $\lambda$ is an eigenvalue of $M$ with a Jordan chain $\vec{v_1},\dots,\vec{v_l}$, then the complex conjugate $\lambda^*$ is also an eigenvalue with Jordan chain $\vec{v_1}^*,\dots,\vec{v_l}^*$. Hence, for non-real eigenvalues, we can take a linear combination of these chains to form the real and imaginary parts.
To this end, let us denote the number of non-real eigenvalues bij $2\tilde{m}$ (this is always even since they occur in pairs) and the number of real eigenvalues by $m$.  We then reorder the eigenvalues and Jordan chains above so that no pairs of complex conjugate eigenvectors occur as the first $\tilde{m}$, and define:
$$
a_i = \Re(\lambda_i) \quad
b_i = \Im(\lambda_i) \\
\vec{v^{(i)R}_{jp}} = \Re(\vec{v^{(i)}_{jp}}) \quad
\vec{v^{(i)I}_{jp}} = \Im(\vec{v^{(i)}_{jp}})
$$
If you then work out the formulas, you get:
$$
M\vec{v^{(i)R}_{jp}} = a_i \vec{v^{(i)R}_{jp}} - b_i \vec{v^{(i)I}_{jp}} + \vec{v^{(i)R}_{j (p-1)}}\\
M\vec{v^{(i)I}_{jp}} = b_i \vec{v^{(i)R}_{jp}} + a_i \vec{v^{(i)I}_{jp}} + \vec{v^{(i)I}_{j (p-1)}}
$$
Hence, setting the real eigenvalues to $c_i$ with Jordan chains $w^{(i)}_{jp}$ and redefining the matrix S to:
\begin{multline}
S = [
\vec{v^{(1)R}_{11}},\vec{v^{(1)I}_{11}},\dots,
\vec{v^{(1)R}_{1 l_{11}}},\vec{v^{(1)I}_{1 l_{11}}},\dots,
\vec{v^{(\tilde{m})R}_{k_\tilde{m} 1}},\vec{v^{(\tilde{m})I}_{k_\tilde{m} 1}},\dots,
\vec{v^{(\tilde{m})R}_{k_\tilde{m} l_{\tilde{m} k_\tilde{m}}}},\vec{v^{(\tilde{m})I}_{k_\tilde{m} l_{\tilde{m} k_\tilde{m}}}},
\\\dots,
\vec{w^{(1)}_{11}},\dots,\vec{w^{(m)}_{k_m l_{m k_m}}}
]
\end{multline}
we now see that we get the following matrix:
$$
J = S^{-1} M S = \left[\begin{array}{ccc|ccc}
K_1\\
& \ddots\\
&& K_\tilde{m}\\
\hline
&&& J_1\\
&&&& \ddots\\
&&&&& J_m
\end{array}\right]
\text{, where:}\quad
K_i = \left[\begin{array}{}
K_{i1}\\
& \ddots\\
&& K_{i k_i}
\end{array}\right]
\text{ and}\\
K_{ij} = \left[\begin{array}{cc|cc|cc|cc}
a_i & b_i & 1 & 0\\
-b_i & a_i & 0 & 1\\
\hline
&& a_i & b_i & \ddots\\
&& -b_i & a_i && \ddots\\
\hline
&&&& \ddots && 1 & 0\\
&&&&& \ddots & 0 & 1\\
\hline
&&&&&& a_i & b_i\\
&&&&&& -b_i & a_i
\end{array}\right]
\in \mathbb{R}^{2l_{ij} \times 2l_{ij}}
$$
and the $J_i$ are as above.
Orthogonal matrices
As stated above, the nice thing about orthogonal matrices is that their eigenvalues satisfy $|\lambda|^2 = 1$ and that all Jordan chains are of length 1.  Hence you can always write $\lambda = e^{i\theta} = \cos\theta + i \sin\theta$ for some $\theta \in \mathbb{R}$ and find a basis of eigenvectors.  This allows you to interpret an orthogonal matrix in a very straightforward way, namely as a series of operations on independent 1- or 2-dimensional subspaces: 


*

*For a non-real eigenvalue the real and imaginary part of an eigenvector span a 2-dimensional space on which the matrix operates as a rotation of an angle equal to the argument of the eigenvalue.

*For a real eigenvalue the matrix acts either as the identity (1) or a reflection (-1) in the hyperplane normal to the corresponding eigenvector.

