# Geodesic of Stiefel manifold

Define a metric on Stiefel manifold $V_{n,p}$ as $$\left<\Delta_1,\Delta_2\right>=\text{tr}\Delta_1^T\left(I-\frac{1}{2}YY^T\right)\Delta_2$$ $\forall \Delta_1,\Delta_2\in T_YV_{n,p}$

• how to calculate geodesic through the variation problem

$$\min\limits_{Y(t)}\int\left<\dot{Y},\dot{Y}\right>^{\frac{1}{2}}dt$$ $Y(t)$ is the curve in $V_{n,p}$.

I guess this is called the canonical inner product. According to [1], and taking $[Q]$, the equivalence class of all $n$-by-$n$ orthogonal matrices to be a point on the Stiefel manifold, we can re-write $\left<\Delta_1,\Delta_2\right>$ as:

\begin{align*} \left<\Delta_1,\Delta_2\right> &= \frac{1}{2} \text{tr} \Bigg( Q \begin{bmatrix} A_1 & -B_1^T \\ B_1 & 0 \end{bmatrix} \Bigg)^T Q \begin{bmatrix} A_2 & -B_2^T \\ B_2 & 0 \end{bmatrix}\Bigg)\\ &=\frac{1}{2} \text{tr} A_1^T A_2 + \text{tr} B_1^T B_2 \end{align*}

I follows from [1] that the orthogonal group geodesic is given by: $$Q(t) = Q(0) \text{exp}\,\,t\begin{bmatrix}A & -B^T \\ B & 0\end{bmatrix}.$$ these are the curves of shortest length in the quotient space.

Minimizing the path length $L$ where $$L=\int\left<\dot{Y},\dot{Y}\right>^{\frac{1}{2}}dt$$ using the calculus of variations as given in the question yields the following geodesic equation: $$\ddot{Y}+\dot{Y}\dot{Y}^T+Y((Y^T\dot{Y})^2+\dot{Y}^T\dot{Y}) = 0.$$ It is reported in [1] that this is tedious to derive though. With appropriate substitutions, the paths of the form $$Y(t) = Q e^{Xt}I_{n,p}$$ where $X=\begin{bmatrix} A & -B^T \\ B & 0\end{bmatrix}$, are found to be the closed form solutions to the geodesic equation for the canonical metric.

The paper also gives a computationally more reasonable algorithm to compute that. First, authors parameterize the same geodesic by an initial position $Y(0)=Y$ and direction $\dot{Y}(0)=H$. By formulating a quadratic eigenvalue problem, they show that the geodesic is given by the following curve: $$Y(t) = Y M(t) + Q N(t)$$ where $QR:=K=(I-YY^T)H$ is the QR-decomposition of $K$ and $M(t)$ and $N(t)$ are given by the matrix exponential:

$$\begin{pmatrix}M(t)\\N(t)\end{pmatrix} = \text{exp}\, t\begin{pmatrix}A & -R^T \\ R & 0\end{pmatrix}\begin{pmatrix}I_p \\ 0\end{pmatrix}$$ $I_p$ is the identity. This can be computed easily by solving a $2p$-by-$2p$ skew-symmetric eigenvalue problem. Hint: Use SVD.

$$\\$$ $$\\$$

Side notes:

• Equivalence class encapsulates all orthogonal matrices with identical first $p$ columns.

[1] Alan Edelman, Tomas A. Arias, And Steven T. Smith, The Geometry Of Algorithms With Orthogonality Constraints, SIAM J. on Matrix Analysis and Applications, 1998