Closed-Form Geodesics in Matrix Manifolds

• Closed-form geodesics are explicit analytical solutions for distance‐minimizing curves on Riemannian manifolds, particularly in matrix manifolds like the Stiefel, Grassmann, and flag manifolds.
• They employ a two-parameter family of Riemannian metrics and offer both block-exponential and reduced formulations to efficiently compute geodesic paths, exponential maps, and logarithm maps.
• The framework integrates Fréchet derivatives with trust-region minimization, yielding robust convergence in optimization and statistical applications on low-dimensional settings.

Closed-form geodesics are explicit analytical solutions for geodesics—locally distance-minimizing curves—on Riemannian manifolds, with particular focus on matrix manifolds such as the Stiefel, Grassmann, and flag manifolds. These formulas enable efficient computation of geodesics, exponential and logarithm maps for a two-parameter family of Riemannian metrics, crucial for algorithms in optimization and statistics on these spaces. The development of closed-form geodesics, together with Fréchet derivatives and trust-region minimization, forms a unified framework for calculating Riemannian logarithms and geodesic distances even when closed-form logarithm maps are unavailable (Nguyen, 2021).

1. Stiefel Manifold and Metric Family

The Stiefel manifold $$St_K(n,p) = \{ Y \in K^{n \times p} : Y^t Y = I_p \}$$ consists of orthonormal $p$-frames in $K^n$ (real or complex). The tangent space at $Y$ comprises matrices $\eta$ such that $Y^t \eta$ is antisymmetric. The Riemannian metric is parameterized by two positive scalars $\alpha_0, \alpha_1$: $g_Y(\omega, \xi) = \Tr_R[ \omega^t (\alpha_0 I_n + (\alpha_1 - \alpha_0) Y Y^t) \xi ]$ The metric interpolates between the embedded metric ($\alpha_1 = \alpha_0$) and canonical metric ($\alpha_1 = \frac12 \alpha_0$). The associated raising operator has the form

$$g_Y(\omega) = \alpha_0 \omega + (\alpha_1-\alpha_0) Y Y^t \omega$$

This parameterization enables closed-form geodesic solutions that work efficiently for low-rank settings ($p \ll n$) (Nguyen, 2021).

2. Geodesic ODE and Closed-Form Solutions

The geodesic $Y(t)$ satisfies the second-order ODE: $$\ddot Y + Y\,\dot Y^t \dot Y + 2(1-\alpha)\,(I - Y Y^t)\,(\dot Y\,\dot Y^t)\,Y = 0, \quad \alpha = \frac{\alpha_1}{\alpha_0}$$ Initial conditions: $Y(0) = Y$, $\dot Y(0) = \eta$, and the decomposition: $$A = Y^t \eta, \qquad (I - Y Y^t)\eta = Q R,\quad Q^t Q = I,\ Q^t Y = 0,\ R \in K^{p \times k}$$ With $S_0=\eta^t\eta = -A^2 + R^t R$ and matrix

$$\hat A = \begin{bmatrix} 2\alpha A & -R^t \ R & 0 \end{bmatrix}$$

the geodesic admits two equivalent closed forms:

• Block-exponential (dimension $(p+k)$):

$$Y(t) = [Y\ Q] \exp(t \hat A) \begin{bmatrix} I_p \ 0 \end{bmatrix} \exp( (1 - 2\alpha)tA )$$

• Reduced (dimension $2p$):

$Y(t) = (Y M(t) + Q N(t))\ \exp( (1 - 2\alpha)tA )$

with

$$[M(t)\ N(t)] = \exp \left( t \begin{bmatrix} 2\alpha A & -R^t \ R & 0 \end{bmatrix} \right) \begin{bmatrix} I_p \ 0 \end{bmatrix}$$

Both representations maintain computational complexity at $O((p+k)^3)$ per matrix exponential, avoiding dependence on the ambient dimension $n$ when $k \ll p$ (Nguyen, 2021).

3. Fréchet Derivatives and Gradient Computation

For the Riemannian logarithm, the problem reduces to finding $(A,R)$ such that $Y(1) \approx Z$, minimizing

$F(A, R) = \frac12 \| Y(A, R;1) - Z \|_F^2$

The gradient is obtained via Fréchet derivatives: $$\mathcal L_f(A, E) = \lim_{h \to 0} \frac{f(A + hE) - f(A)}{h} = \sum_{i=1}^\infty f_i \sum_{a+b=i-1} A^a E A^b$$ This yields gradients with computational cost comparable to evaluating one exponential plus its Fréchet derivative (approximately three times the cost of a single exponential). The trace–Fréchet adjoint identity is applied: $\Tr[ C\,\mathcal L_f(A, E)\,D ] = \Tr[ \mathcal L_f(A, D C)\,E ]$ The gradient splits into blocks: $$\nabla_A F = (1-2\alpha)\,\mathrm{skew}(\mathring L) + 2\alpha\,\mathrm{skew}(\hat L_{11}),\qquad \nabla_R F = \hat L_{21} - \hat L_{12}^t$$ with $\mathrm{skew}(M) = \frac12(M-M^t)$, enabling efficient exact gradient evaluation (Nguyen, 2021).

4. Trust-Region Minimization and the Riemannian Logarithm

The Riemannian logarithm is computed by minimizing $F(A, R)$ over antisymmetric $A$ and arbitrary $R$. Standard interior-point or trust-region solvers (e.g., SciPy’s “trust-krylov”) are employed, leveraging the closed-form gradient. When $Z$ is within the injectivity radius of $Y$, the Hessian of $F$ is positive definite and convergence to the unique minimizing logarithm is rapid. For moderately large distances, the algorithm often still recovers nearly minimal “shooting” geodesics.

This approach extends to quotient manifolds:

• On $St(n,p)/H$, for a subgroup $H \subset O(p)$ (e.g., Grassmann or flag), restrict $A$ to the horizontal subspace and adapt cost functions to measure ambient-coordinate distance.
• The trust-region and Fréchet framework applies identically for these cases (Nguyen, 2021).

5. Grassmann and Flag Manifolds: Trigonometric Formulas

For the Grassmann manifold $Gr(p, n) = St(n, p) / O(p)$, the horizontal space is $\{\eta : Y^t\eta = 0\}$. For a purely horizontal geodesic ($A=0$): $\Exp^{Gr}_Y(\eta) = \Exp^{St}_Y(\eta) = Y\,\cos(S^{1/2}) + \eta\,S^{-1/2} \sin(S^{1/2}), \quad S = \eta^t\eta$ Given $Y,Z \in St(n,p)$, compute the SVD: $$Y^t Z = U\,\Sigma\,V^t, \quad \Sigma = \mathrm{diag}(\sigma_i) \in [0,1]^{p \times p}$$ The unique minimizing logarithm for all $\sigma_i \in (0,1)$: $\Log_Y^{Gr}(Z) = (Z - Y Y^t Z)\; V\; (I - \Sigma^2)^{-1/2} \arccos(\Sigma) U^t$ The geodesic distance is

$d_{Gr}(Y, Z) = \sqrt{\alpha_0\,\Tr[\arccos^2(\Sigma)]}$

These are analogous to the $\cos\theta$, $\sin\theta$ formulas on the sphere (Nguyen, 2021).

6. Algorithmic Complexity, Stability, and Applications

Each trust-region iteration requires one or two matrix exponentials (size at most $(p+k)$) plus Fréchet derivatives; cost per exponential is $O((p+k)^3)$. For low-rank problems ($p \ll n$), all computations occur in small dimension ($k \le p$). Fréchet-derivative routines (e.g., SciPy's expm_frechet) compute both the exponential and its Fréchet derivative at threefold the cost of the exponential alone.

The approach avoids integration along curves or multi-step shooting, instead solving a small-dimensional nonlinear least-squares problem. Numerical experiments on Stiefel and flag manifolds demonstrate robust convergence for distances up to and beyond half the injectivity radius. The Hessian remains well-conditioned for points near one another, deteriorating only as one approaches the cut locus.

One immediate application is computation of the Riemannian center of mass (Karcher mean): minimizing $\frac1N \sum_i d^2(X, Q_i)$ by gradient descent, where each step requires $N$ logarithm computations—now feasible in low-dimension settings. The methods are conceptually and algorithmically unified, supporting practical statistical and optimization tasks on matrix manifolds (Nguyen, 2021).