Matrix in some orthonormal basis: $\frac{1}{6}\begin{bmatrix} 3 &-1 &-1 &-5\\ 1 &3 & -5 &1\\ 1 &5 &3 &-1\\ 5 &-1& 1 &3 \end{bmatrix}$
I know how to do it for $\mathbb{R}^3$, but I don't know how rotation matrices should look for $\mathbb{R}^4$.
Matrix in some orthonormal basis: $\frac{1}{6}\begin{bmatrix} 3 &-1 &-1 &-5\\ 1 &3 & -5 &1\\ 1 &5 &3 &-1\\ 5 &-1& 1 &3 \end{bmatrix}$
I know how to do it for $\mathbb{R}^3$, but I don't know how rotation matrices should look for $\mathbb{R}^4$.
Following Claim 3 in this answer, start by finding an invariant plane of the transformation. Its characteristic polynomial is $(\lambda^2-\lambda+1)^2$, so its eigenvalues are $\frac12(1\pm\sqrt3i$), each with multiplicity $2$. Taking the real and imaginary parts of one of the corresponding complex eigenvectors, say $(0,-26,1+15\sqrt3i,5-3\sqrt3i)^T$ gives us a basis for an invariant plane, namely, $(0,-26,1,5)^T$ and $(0,0,15,-3)^T$ (I’ve dropped the common factor of $\sqrt3$ from the latter to reduce clutter.) These two vectors are already orthogonal, so normalize them and use Gram-Schmidt to extend this basis to an orthonormal basis of $\mathbb R^4$. After applying this change of basis to your original matrix, you should end up with one of the form $$\begin{bmatrix}R_1&0\\0&R_2\end{bmatrix},$$ where $R_1$ and $R_2$ are $2\times2$ rotation matrices. To recover the two rotations relative to the original basis, set each $R_k$ block in turn to the identity matrix and apply the inverse change of basis.