Decompose rotation matrix to plane of rotation and angle I would like to decompose an $n$-dimensional orthogonal rotation matrix (restricting to simple rotation with a single plane of rotation) to the two basis vectors of the plane of rotation, and an angle of rotation.
The common method is decomposing the rotation matrix to an axis and angle, but this doesn't work in higher dimensions.
For example in $\mathbb{R}^3$ given the rotation matrix 
$R_{xy}=\begin{bmatrix}\cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}$
it's obvious that the plane of rotation is the $xy$-plane spanned by basis vectors $b_0 = (1,0,0)$ and $b_1=(0,1,0)$ and the angle of rotation is $\theta$.
However decomposing it mathematically is rather challenging.
What is the solution for a general (restricted to a single, but arbitrary plane) rotation matrix in $\mathbb{R}^n$?

This is a surprisingly interesting question. The solution in $\mathbb{R}^3$ is known (for example, by converting the rotation matrix $\mathbf{R}$ to a versors AKA unit quaternions to obtain the axis and the angle, followed by negating and permuting the axis vector components to get an auxiliary vector which can always (assuming a nonzero axis vector) be normalized to get the second basis vector; the third basis vector being the cross product of the other two) but in higher dimensions sounds hard.
It turns out that if $\mathbf{R}$ represents a rotation on a plane ($\mathbb{R}^2$ subspaces) in $\mathbb{R}^n$, then
$$\operatorname{Tr}\mathbf{R} = n - 2 + 2 \cos \theta$$
where $\operatorname{Tr}\mathbf{R} $ refers to trace or matrix $\mathbf{R}$, the sum of elements on the main diagonal on the $n\times n$ matrix $\mathbf{R}$. Note that $\mathbf{R}$ must be an orthogonal matrix (orthonormal matrix), i.e. $\mathbf{R}^T\mathbf{R}=\mathbf{I}$. In other words, the rotation angle $\theta$ fulfills
$$\cos \theta = 1 - \frac{n - \operatorname{Tr}\mathbf{R}}{2}$$
and since $\sin\theta = \pm\sqrt{1-\cos^2\theta}$, sign depending on which half of full rotation $\theta$ is in,
$$\sin\theta = \pm\frac{1}{2}\sqrt{\left(n - \operatorname{Tr}\mathbf{R}\right)\left(4 - n + \operatorname{Tr}\mathbf{R}\right)}$$
(I found it extremely curious and funny that the rotation angle $\theta$ would be so easy to find; it is the main reason I wrote this answer, in the hopes that this joggles some mathematicians mind enough to point us towards an algorithm that can be used to obtain the solution to the problem at hand.)
The problem left is how to find the bivectors defining the plane of rotation.
The Wikipedia article on plane of rotation mentions that the complex roots of the characteristic polynomial of matrix $\mathbf{R}$ correspond to the planes of rotation -- here, we expect only one complex conjugate root pair --; i.e. eigenplanes of the matrix $\mathbf{R}$.
Unfortunately, my math-fu fails me here -- I've never before encountered this type of problem, and I am not a mathematician --, so I am hoping a real mathematician would complete (or fix, if I am in error, as that tends to occur a bit too often) the path to an actual solution algorithm.
