I have a $4\times 4$ matrix $M$ which represents a general 4-dimensional rotation. $$ M = \pmatrix{a_{11} &a_{12} &a_{13} &a_{14}\\a_{21} &a_{22} &a_{23} &a_{24}\\a_{31} &a_{32} &a_{33} &a_{34}\\a_{41} &a_{42} &a_{43} &a_{44}} $$ There are also six 4D basis rotation matrices $M_{XY}$, $M_{YZ}$, $M_{ZX}$, $M_{XW}$, $M_{YW}$, $M_{ZW}$.
How to decompose 4D rotation $M$ into basis rotations (i.e., to find corresponding angles $T_{XY}$, $T_{YZ}$, $T_{ZX}$, $T_{XW}$, $T_{YW}$, $T_{ZW}$)?
Wikipedia says that any 4D rotation given by its matrix can be decomposed into two isoclinic rotations. However, I do not know if it is possible to decompose them further into basis rotations. Also, there may be another (simpler) way.
You should be able to build your basis matrices by doing something similar to the Gram-Schmidt process: by multiplying by suitably-chosen rotation matrices you can set each of the off-diagonal coefficients in turn equal to 0. For instance, consider applying a $M_{ZW}$ rotation by some angle $\theta$ to the matrix: $$M' = M_{ZW}(\theta)\ M = \pmatrix{1&0&0&0\\0&1&0&0\\0&0&\cos\theta&-\sin\theta\\0&0&\sin\theta&\cos\theta}\pmatrix{a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\a_{31}&a_{32}&a_{33}&a_{34}\\a_{41}&a_{42}&a_{43}&a_{44}} =\pmatrix{a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\*&*&a_{33}\cos\theta-a_{43}\sin\theta&a_{34}\cos\theta-a_{44}\sin\theta\\*&*&a_{33}\sin\theta+a_{43}\cos\theta&a_{34}\sin\theta+a_{44}\cos\theta}$$ (Here the *s are 'don't-cares' - they have a value, but it's essentially unimportant.) Now, by choosing $\theta$ correctly, you can make $M'_{34} = a_{34}\cos\theta-a_{44}\sin\theta$ equal to zero without affecting any other upper-triangular elements. Once you've got that matrix, you can continue in a similar fashion by finding a suitable $M_{YW}(\phi)$ to zero the $M_{24}$ element, then an $M_{XW}$ to zero the $M_{14}$ element, etc.; repeating this procedure for each combination of unequal coordinates (in some canonically-defined order) will eventually lead you to an orthogonal (by hypothesis) lower-triangular matrix, which must therefore be the identity.