How to decompose a rotation matrix in $SO(n)$ to a product of rotation matrices acting on $2D$ subspaces. I read in a paper (end of page 3) that a rotation matrix in dimension $n$ can be parametrized as a product of $n(n-1)/2$ rotation matrices of $2D$ subspaces. Each subspace is generated by a pair of canonical basis vectors:
$$
R = \prod\limits_{i=1}^{n-1}\prod\limits_{j=i+1}^n R_{ij}(\theta_{ij})
$$
The rotation matrices in the product do not commute, so if I take the product of $2$ such decompositions, the decomposition of the product is not trivially obtained from that of the factors.
My question is if there is an algorithm to obtain such a decomposition from an arbitrary rotation matrix in $SO(n)$.
 A: You can do this with a relatively minor modification of the Gaussian elimination procedure for inverting matrices.  Let $\ A=\big(a_{ij}\big)\  $ be any non-singular $\ n\times n\ $ matrix with either $\ a_{ii}\ne0\ $ or $\ a_{ji}\ne0\ $.  Let $\ R_{ij}\big(\theta_{ij}\big)\ $ be the rotation matrix
$$
\begin{array}{cccccccc}
&&&&&\text{column }\ i&&\text{column }\ j\\
&&&&&\downarrow&&\downarrow
\end{array}\\
\begin{matrix}
\ \\
\ \\
\ \\
\text{row }\ i\rightarrow\\
\ \\
\text{row }\ j\rightarrow\\
\ \\
\ \\
\end{matrix}\pmatrix{1&0&\dots&0&\dots&0&\dots&0\\
0&1&\dots&0&\dots&0&\dots&0\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\vdots\\
0&0&\dots&\cos\theta_{ij}&\dots&\sin\theta_{ij}&\dots&0\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\
0&0&\dots&-\sin\theta_{ij}&\dots&\cos\theta_{ij}&\dots&0\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\
0&0&\dots&0&\dots&0&\dots&1}\ ,
$$
where $\ \sin\theta_{ij}=\frac{a_{ji}}{a_{ii}^2+a_{ji}^2}\ $, and $\ \cos\theta_{ij}=\frac{a_{ii}}{a_{ii}^2+a_{ji}^2}\ $ if $\ a_{ii}\ne0\ $, or $\ \theta_{ij}=\frac{\pi}{2}\ $ otherwise (when $\ a_{ji}\ $ will be non-zero, by assumption). If $\ B=R_{ij}\big(\theta_{ij}\big)A\ $, and $\  b_{rs}\ $ is the entry in the $\ r^\text{th}\ $ row and $\ s^\text{th}\ $ column of $\ B\ $ then
\begin{align}
b_{ji}&=-a_{ii}\sin\theta_{ij}+a_{ji}\cos\theta_{ij}\\
&=0\,
\end{align}
and
\begin{align}
b_{ii}&=a_{ii}\cos\theta_{ij}+a_{ji}\sin\theta_{ij}\\
&=1\,
\end{align}
when $\ a_{ii}\ne0\ $, or $\ b_{ji}=a_{ii}=0\ $ and $\ b_{ii}=a_{ji}\ne0\ $ otherwise.
If  $\ a_{i_11}, a_{i_21}, \dots, a_{i_r1}\ $ are the non-zero entries in the first column of $\ A\ $, then by applying a sequence $\ R_{i_r1}\big(\theta_{i_r1}\big),$$\,R_{i_{r-1}1}\big(\theta_{i_{r-1}1}\big),$$\, \dots,$$R_{i_11}\big(\theta_{i_11}\big)\ $ of rotations of the above form to $\ A\ $, we can obtain a matrix
$$
C=\left(\prod_{k=1}^rR_{i_k1}\big(\theta_{i_k1}\big)\right)A
$$
whose first column is the first column of the $\ n\times n\ $ identity matrix.  For notational convenience, I'll fill in the gaps in the above product by setting $\ \theta_{i1}=0 $ for any $\ i\ $ such that $\ a_{i1}=0\ $.  Since $\ R_{i1}\big(\theta_{i1}\big)\ $ is then the identity matrix, we can write the above product as
$$
C=\left(\prod_{i=2}^nR_{i1}\big(\theta_{i1}\big)\right)A\ .
$$
Now applying this procedure successively to the $\ c\times c\ $ submatrices formed by the last $\ c\ $ rows and columns of $\ C\ $, and the subsequent matrices derived from it, for $\ c=n-1,n-2,\dots,2\ $, we can obtain a matrix
$$
D=\left(\prod_{j=1}^{n-1}\prod_{i=n+1-j}^nR_{ij}\big(\theta_{ij}\big)\right)A
$$
which is upper triangular with entries all $\ 1\ $ along the main diagonal.  If $\ A\ $ is a rotation matrix, however, then so is $\ D\ $, and the only upper triangular rotation matrix with entries all $\ 1\ $ along the main diagonal is the identity matrix.  So we have
$$
I=\left(\prod_{j=1}^{n-1}\prod_{i=n+1-j}^nR_{ij}\big(\theta_{ij}\big)\right)A\ ,
$$
and multiplying this equation by the inverses of the rotation matrices, we get
$$
A=\prod_{j=1}^{n-1}\prod_{i=0}^{j-1}R_{n-i\,n-j}\big({-}\theta_{n-i\,n-j}\big)\ .
$$
