Rotation along geodesic in n-dimensional sphere

Question

Suppose I have two unit-length vectors $a, b$ in n-dimensional Euclidean space - so they are two points on the n-dimensional unit sphere.

Now there is a geodesic between these two points, which extends to a great circle around the unit sphere.

How can I calculate the point at a specified angle $\theta$ from $a$ along this geodesic defined by the two points?

---

Clearly for $\theta = 0$ we have $a$, and for $\theta = \cos^{-1}(\left<a, b\right>)$ we have $b$.

There must be some rotation matrix $R$ which represents the rotation from $a$ to $b$. Then my target point would be $R^{\frac {\theta} {\cos^{-1}(\left<a, b\right>)}}a$. How can I derive this rotation matrix?

Perhaps I can represent the great circle as sweeping a rotation parallel to the plane defined by $a$ and $b$, as in Rodrigues' formula in 2d or 3d (generalised here to higher dimensions?).

Answer 1

There’s no such thing as “the rotation from $a$ to $b$”; infinitely many different rotations take $a$ into $b$. For the geodesic, you want the one that stays in the plane spanned by $a$ and $b$. This is $a\cos\theta+c\sin\theta$, where

$$c=\frac{b-(a\cdot b)a}{|b-(a\cdot b)a|}$$

is the unit vector in the plane spanned by $a$ and $b$ that’s orthogonal to $a$ and has positive scalar product with $b$.

For $\theta=\arccos(a\cdot b)$, this yields

\begin{eqnarray} a\cos\theta+c\sin\theta&=&a\cos\arccos(a\cdot b)+\frac{b-(a\cdot b)a}{|b-(a\cdot b)a|}\sin\arccos(a\cdot b)\\ &=&a(a\cdot b)+\frac{b-(a\cdot b)a}{|b-(a\cdot b)a|}\sqrt{1-(a\cdot b)^2}\\ &=&a(a\cdot b)+b-(a\cdot b)a\\ &=&b\;. \end{eqnarray}

Answer 2

The approach sketched below the line in the question does indeed work.

The geodesic great circle connecting $a$ and $b$ can be parameterised by angle $\theta$ as a rotation of the point $a$ around the origin through the plane containing the points $a$, $b$, and the origin.

$a$ and $b$ are unit vectors together describing a plane of rotation. We can orthogonalise them (using e.g. Gram-Schmidt procedure), obtaining $n_1, n_2$ orthonormal vectors in the same hyperplane. Then the generator of rotations in this hyperplane is

$$ r_{n_1, n_2} = n_2n_1^T - n_1n_2^T $$

Now rotating by angle $\theta$ is just $e^{\theta r}$ which gives a generalisation of Rodrigues' formula for rotation

$$ R_{n_1, n_2, \theta} = I + (n_2n_1^T - n_1n_2^T) \sin \theta + (n_1n_1^T + n_2n_2^T) (\cos \theta - 1) $$

Finally the desired point is the application of this rotation to $a$, namely $R_{n_1, n_2, \theta} a$.