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)$.
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)\ . $$