Parametrization of the special orthogonal group in $n$ dimensions We decompose a rotation $O$ in $n$ dimensions into a product of $\binom{n}{2}$ elementary rotations each one in a plane $(i,j)$ where $1 \le i < j \le n$. in other words we write:
\begin{equation}
O = \prod_{i=1}^n \prod_{j=i+1}^n O^{(i,j)}
\end{equation}
where
\begin{equation}
\left(
\begin{array}{cc}
O^{(i,j)}_{i,i} & O^{(i,j)}_{i,j} \\ O^{(i,j)}_{j,i} & O^{(i,j)}_{j,j}
\end{array}
\right)
=
\left(
\begin{array}{cc}
\cos(\varphi_{\ell_{i,j}}) & -\sin(\varphi_{\ell_{i,j}}) \\ \sin(\varphi_{\ell_{i,j}}) & \cos(\varphi_{\ell_{i,j}})
\end{array}
\right)
\end{equation}
The remaining diagonal and cross-diagonal elements of $O^{(i,j)}$ are equal to one and to zero respectively. Here $\ell_{i,j} := (i-1)(2 n-i)/2  + (j-i)$.
Now, with the help of Mathematica I have checked that the following equality holds:
\begin{eqnarray}
O_{i,j} = \sum\limits_{i \le i_1 \le i_2 \le \cdots \le i_{j-1} \le n}
{\mathcal A}^{(1)}_{i,i_1} {\mathcal A}^{(2)}_{i_1,i_2} \cdots {\mathcal A}^{(j)}_{i_{j-1},j}
\end{eqnarray}
subject to $2 \le i_1,  3 \le i_2, 4 \le i_3, \ldots, j \le i_{j-1}$ when $j\ge 2$ and $O_{i,j}= {\mathcal A}^{(1)}_{i,1}$ when $j=1$.
Here the matrix elements ${\mathcal A}^{(1)}$ depend only on the first $(n-1)$ angles, the matrix elements ${\mathcal A}^{(2)}$ depend only on the second $(n-2)$ angles and so on and so forth. The matrix elements ${\mathcal A}$ read:
\begin{eqnarray}
&&{\mathcal A}^{(p)}_{q,\ell} := \\[8pt] &&
\begin{cases}
\left(1_{q=p} + 1_{q>p} \sin(\varphi_{B_p + q-p-1})\right) \prod\limits_{\xi = B_p + q-p}^{E_p} \cos(\varphi_\xi) & \text{if } \ell=p \\[6pt]
0 & \text{if } p+1 \le \ell < q \\[6pt]
\cos(\varphi_{B_p + q-p-1}) & \text{if $q=\ell$ and $q > p$} \\[6pt]
-\left(1_{q=p} + 1_{q>p} \sin(\varphi_{B_p + q-p-1})\right) \prod\limits_{\xi = B_p + q-p}^{B_p+\ell-p-2} \cos(\varphi_\xi) \cdot \sin(\varphi_{B_p+\ell-p-1}) & \text{if $q+1 \le \ell \le n$}
\end{cases}
\end{eqnarray}
and $B_p := (p-1) (2 n-p)/2 + 1$ and $E_p := p (2 n - p -1)/2$ for $p=1,\ldots,n-1$. Note that $B_p + n-p-1 = E_p$.
Can anybody provide a proof for this factorization?
 A: Note that the elementary rotations are sparse matrices. Yet their product produces a matrix that is not sparse.  Interestingly enough there are certain sub-sequences in the sequences of all matrices $O^{i,j}$, sub-sequences whose product produces fairly simple matrices. To be precise the following equality holds:
\begin{equation}
\prod\limits_{j=i+1}^n O^{(i,j)} = 
\left(
\begin{array}{ll}
1 & 0 & \cdots & 0 & \cdots & 0  \\
\vdots & \ddots \\
0 & \cdots & 1 & 0 & \cdots & 0 \\
0 & \cdots & 0 & {\mathcal A}^{(i)}_{i,i} & \cdots & {\mathcal A}^{(i)}_{i,n} \\
0 & \cdots & 0 & {\mathcal A}^{(i)}_{i+1,i} & \cdots &  {\mathcal A}^{(i)}_{i+1,n} \\
\vdots &  & \vdots & \vdots & & \vdots \\
0 & \cdots & 0 & {\mathcal A}^{(i)}_{n,i} & \cdots & {\mathcal A}^{(i)}_{n,n}
\end{array}
\right)_{n \times n}
\end{equation} 
In here the element ${\mathcal A}^{(i)}_{i,i}$ sits at the $(i-1,i-1)$ position. I haven't proven this identity but I generated the matrices symbolically in Mathematica and checked that it indeed holds. Intuitively we can understand that such an equality should hold.  Indeed, multiplying $O^{i,i+1}$ by $O^{i,i+2}$ will produce a three by three matrix. Multiplying that matrix further by $O^{i,i+3}$ will produce a four by four matrix. Finally taking the whole product as on the right hand side will produces a $n-i+1$ by $n-i+1$ matrix. Interestingly enough that matrix does not contain any sums but instead each element is just a product of cosine or sine functions. Now, let us take a closer look at the matrix $\left\{{\mathcal A}^{(i)}_{q,l}\right\}_{q,l=i,i}^{n,n}$. That matrix is almost upper diagonal, except for the first column which is not zero. For example for $n=6$ and $i=3$ we have $\left(B_{3}, E_{3}\right) = \left(10,12\right)$ and that matrix reads:
\begin{equation}
\begin{array}{cccc}
 \cos (\phi_{10}) \cos (\phi_{11}) \cos (\phi_{12}) & -\sin (\phi_{10}) & -\cos (\phi_{10}) \sin (\phi_{11}) & -\cos (\phi_{10}) \cos (\phi_{11}) \sin (\phi_{12}) \\
 \cos (\phi_{11}) \cos (\phi_{12}) \sin (\phi_{10}) & \cos (\phi_{10}) & -\sin (\phi_{10}) \sin (\phi_{11}) & -\cos (\phi_{11}) \sin (\phi_{10}) \sin (\phi_{12}) \\
 \cos (\phi_{12}) \sin (\phi_{11}) & 0 & \cos (\phi_{11}) & -\sin (\phi_{11}) \sin (\phi_{12}) \\
 \sin (\phi_{12}) & 0 & 0 & \cos (\phi_{12}) \\
\end{array}
\end{equation}
So what we end up with is a product of matrices from the left to the right , such that the left most matrix is $n$ dimensional, the second matrix acts in the last $n-1$ dimensions, the third matrix acts in the last $n-2$ dimensions and so on and so forth until the last matrix is just an identity matrix. Since the first column of all the last $n-1$ matrices is $(1,0,\cdots,0)$ it is clear that the first column of the first matrix will be equal to the first column of the final product. Likewise, since the second column of the last $n-2$ matrices is $(0,1,0,\cdots,0)$ it is clear that the second column of the final product willl involve elements from the first and from the second matrix only. Likewise we easily conclude that the $i$th column of the final product involves elements from the first $i$ matrices. Now, from the very definition of matrix multiplication and taking into account the structure of our matrices  we get the formula for $O_{i,j}$. 
