Define a metric on Stiefel manifold $V_{n,p}$ as
$$\langle \Delta_1,\Delta_2 \rangle = \operatorname{tr} \left( \Delta_1^T\left(I-\frac{1}{2}YY^T\right)\Delta_2 \right)$$
for all $\Delta_1, \Delta_2 \in T_Y V_{n,p}$. How to calculate the geodesic through the variational problem
$$ \min\limits_{Y(t)} \int \langle\dot{Y},\dot{Y}\rangle^{\frac12} {\rm d}t $$
where $Y(t)$ is the curve in $V_{n,p}$? Any advice is helpful.
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:
I also recommend the this recent paper which derives an efficient retraction and vector transport map for the Stiefel manifold based upon Cayley transform.
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
