Decompose a 4D rotation as a pair of quaternions I'm trying to rotate points in 4D. In fact I have been able to rotate using $6$ rotation matrices. I've managed to perform 3D rotations using rotation quaternions. I understand that I can use plane angle rotations using a pair of quaternions and the vector (used as a quaternion).
$$Z_lVZ_r$$
I understand that $Z_l$ and $Z_r$ represent the $6$ rotational planes but how do I construct $Z_l$ and $Z_r$? And can I simply add the rotations together like quaternion 3D rotations?
 A: Suppose $R$ is a rotation of $\mathbb{H}\cong\mathbb{R}^4$ in the oriented planes with orthonormal bases $\{a,b\}$ and $\{c,d\}$ by angles $\theta$ and $\phi$ respectively. (Assume the ordered basis $\{a,b,c,d\}$ induces the same orientation of space as does $\{1,\mathbf{i},\mathbf{j},\mathbf{k}\}$.) There exist unit quaternions $u,v\in S^3$ such that $R(x)=uxv$.
First, solve the system of congruences
$$ \begin{cases} \alpha+\beta\equiv \theta \\ \alpha-\beta\equiv\phi \end{cases} \pmod{2\pi}$$
Up to integer multiples of $2\pi$ in each coordinate, the solution is $(\frac{1}{2}(\theta+\phi),\frac{1}{2}(\theta-\phi))$.
Then $R=ST$ where $S$ is the left isoclinic rotation by $\theta$ and $T$ is the right isoclinic rotation by $\phi$ in the $\{a,b\}$- and $\{c,d\}$-planes. That is, they may be represented by matrices
$$ \begin{array}{l} S=\begin{bmatrix}\cos(\theta) & -\sin(\theta) & 0 & 0 \\ \sin(\theta) & \phantom{-}\cos(\theta) & 0 & 0 \\ 0 & 0 & \cos(\phi) & -\sin(\phi) \\ 0 & 0 & \sin(\phi) & \phantom{-}\cos(\phi) \end{bmatrix}, \\[7pt] T=\begin{bmatrix}\cos(\theta) & -\sin(\theta) & 0 & 0 \\ \sin(\theta) & \phantom{-}\cos(\theta) & 0 & 0 \\ 0 & 0 & \cos(-\phi) & -\sin(-\phi) \\ 0 & 0 & \sin(-\phi) & \phantom{-}\cos(-\phi) \end{bmatrix} \end{array}$$
with respect to the orthonormal basis $\{a,b,c,d\}$ (not with respect to the usual basis $\{1,\mathbf{i},\mathbf{j},\mathbf{k}\}$).
Note $ST=TS$. Define $u$ and $v$ such that $ua=b$ and $cv=d$ (these are easy to solve: $u=b\overline{a}$ and $v=\overline{c}d$). Then $S(x)=\exp(\alpha u)x$ and $T(x)=x\exp(\beta v)$ so we conclude that
$$ R(x)=\exp\left(\frac{\theta+\phi}{2}b\overline{a}\right)\,x\,\exp\left(\frac{\theta-\phi}{2}\overline{c}d\right). $$
A: This is the direct method for 4D rotation factorisation given in Perez-Gracia and Thomas's paper. Given the rotation matrix $R$ compute
$$S=\frac{R-A_1RA_1-A_2RA_2-A_3RA_3}4$$
where
$$A_1=\begin{bmatrix}0&0&0&-1\\0&0&-1&0\\0&1&0&0\\1&0&0&0\end{bmatrix}$$
$$A_2=\begin{bmatrix}0&0&1&0\\0&0&0&-1\\-1&0&0&0\\0&1&0&0\end{bmatrix}$$
$$A_3=\begin{bmatrix}0&-1&0&0\\1&0&0&0\\0&0&0&-1\\0&0&1&0\end{bmatrix}$$
Then the decomposition's right factor is
$Y=\frac S{(\det S)^{1/4}}$, while the left factor is $X=RY^T$. $X$ and $Y$ are left- and right-isoclinic rotations respectively, with corresponding quaternions $x$ and $y$ given by their left-most columns; $R$ may then be effected on a quaternion $q$ as $q'=xqy$.
