Two non-co-linear vectors define a 2D subspace that passes through the origin. In 3D, you can represent a 2D subspace by its Normal. It's very easy to define the angle between the two planes: $\theta = \cos^{-1}(N_1\cdot N_2)$ (assuming unit length) and it's a simple matter to find a new plane "in between" by interpolating to find a new normal: $N_t = \frac{\sin((1-t)\theta)}{\sin(\theta)}N_1 + \frac{\sin(t\theta)}{\sin(\theta)}N_2$. Essentially, we're rotating the new Normal vector within the 2D sub-surface defined by the two original Normals.
Rotating a vector this way extends easily into higher dimensions, but a Normal defines a hyper-plane, not a linear, 2D subspace. So I can't use the same approach. If I have two pairs of orthogonal vectors defining two distinct linear 2d subspaces, is there a way to measure the angle between them or find a new subspace that is an interpolation between the two?
Edit:
Working with the problem, I realized that even in 3D, interpolating the normals wouldn't be sufficient for my problem, but I would actually need to interpolate the space defined by the vector pairs. It's still straightforward to do in 3D. Referring to the image below, if I have vectors $X$ and $Y$ and wish to interpolate to $X'$ and $Y'$, it's simple if I find vectors $Z$ and $Z'$ by the cross-product approach. Then we can find $F = \left[X\ Y\ Z\right]$ and $F'=\left[X'\ Y'\ Z'\right]$. The rotation matrix making the full transform from $F$ to $F'$ would just be $F'F^T$. From that rotation matrix, it's easy to extract the axis and angle of rotation and then simply interpolate through the angle.
In higher dimensions, cross-product doesn't work. However, we can still find a point in terms of the first frame of reference and project it into the second: $P(1) = \left[X'(X\cdot P(0) + Y'(Y\cdot P(0))\right]$. Given any pair of $P(0)$ and $P(1)$, I can use spherical interpolation find a point on the shortest arc between the two points, but that's not what I need. I need to interpolate the 2D space itself and find the corresponding point within the plane. If I need to add any number of extra axes because there are multiple equally short paths, that's fine.
Best Answer
I think this is the same answer as yasmar's, but from a more computational perspective.
It makes sense that the pair of vectors in the two subspaces forming the least angle with each other is distinguished and these vectors should be rotated into each other if the interpolation is to follow a geodesic. To find these vectors, you can find the stationary values of the dot product between arbitrary vectors in the subspaces,
$$(\cos\phi\vec{a}_1+\sin\phi\vec{a}_2)\cdot(\cos\theta\vec{b}_1+\sin\theta\vec{b}_2)\;,$$
where $\vec{a}_1$ and $\vec{a}_2$ form an orthonormal basis for the first subspace and likewise $\vec{b}_1$ and $\vec{b}_2$ for the second one. (You can easily calculate these by orthonormalization if you don't already have them.) With $g_{ij}:=\vec{a}_i\cdot\vec{b}_j$, this becomes
$$ g_{11}\cos\phi\cos\theta+ g_{12}\cos\phi\sin\theta+ g_{21}\sin\phi\cos\theta+ g_{22}\sin\phi\sin\theta\;. $$
Differentiating with respect to $\theta$ and $\phi$ yields
$$ \begin{eqnarray} -g_{11}\cos\phi\sin\theta +g_{12}\cos\phi\cos\theta -g_{21}\sin\phi\sin\theta +g_{22}\sin\phi\cos\theta &=&0\;,\\ -g_{11}\sin\phi\cos\theta -g_{12}\sin\phi\sin\theta +g_{21}\cos\phi\cos\theta +g_{22}\cos\phi\sin\theta &=&0\;, \end{eqnarray} $$
and adding and subtracting these equations leads to
$$ \begin{eqnarray} (g_{22}-g_{11})\sin(\phi+\theta)+(g_{12}+g_{21})\cos(\phi+\theta)&=&0\;,\\ (g_{22}+g_{11})\sin(\phi-\theta)+(g_{12}-g_{21})\cos(\phi-\theta)&=&0\;,\\ \end{eqnarray} $$
from which $\phi+\theta$ and $\phi-\theta$, and hence $\phi$ and $\theta$, can be determined.
That gives you the vectors $\vec{a}_<:=\cos\phi\vec{a}_1+\sin\phi\vec{a}_2$ and $\vec{b}_<:=\cos\theta\vec{b}_1+\sin\theta\vec{b}_2$ that should be rotated into each other. The vectors $\vec{a}_>$ and $\vec{b}_>$ with $\vec{a}_> \perp \vec{a}_<$ and $\vec{b}_> \perp \vec{b}_<$ also fulfill $\vec{a}_> \perp \vec{b}_<$ and $\vec{b}_> \perp \vec{a}_<$ (else the stationary point wouldn't be stationary). You can choose their signs so as to form the smaller of the two possible angles between them.
Now these pairs of vectors form two planes that are orthogonal to each other and within which the interpolation rotations should take place, with $\vec{a}_<$ being rotated into $\vec{b}_<$ and $\vec{a}_>$ being rotated into $\vec{b}_>$, both with (different) linear rates of change of the respective rotation angles.