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

$p$0

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

2. Geodesic ODE and Closed-Form Solutions

The geodesic $p$2 satisfies the second-order ODE: $p$3 Initial conditions: $p$4, $p$5, and the decomposition: $p$6 With $p$7 and matrix

$p$8

the geodesic admits two equivalent closed forms:

• Block-exponential (dimension $p$9):

$K^n$0

• Reduced (dimension $K^n$1):

$K^n$2

with

$K^n$3

Both representations maintain computational complexity at $K^n$4 per matrix exponential, avoiding dependence on the ambient dimension $K^n$5 when $K^n$6 (Nguyen, 2021).

3. Fréchet Derivatives and Gradient Computation

For the Riemannian logarithm, the problem reduces to finding $K^n$7 such that $K^n$8, minimizing

$K^n$9

The gradient is obtained via Fréchet derivatives: $Y$0 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: $Y$1 The gradient splits into blocks: $Y$2 with $Y$3, enabling efficient exact gradient evaluation (Nguyen, 2021).

4. Trust-Region Minimization and the Riemannian Logarithm

The Riemannian logarithm is computed by minimizing $Y$4 over antisymmetric $Y$5 and arbitrary $Y$6. Standard interior-point or trust-region solvers (e.g., SciPy’s “trust-krylov”) are employed, leveraging the closed-form gradient. When $Y$7 is within the injectivity radius of $Y$8, the Hessian of $Y$9 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 $\eta$0, for a subgroup $\eta$1 (e.g., Grassmann or flag), restrict $\eta$2 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 $\eta$3, the horizontal space is $\eta$4. For a purely horizontal geodesic ($\eta$5): $\eta$6 Given $\eta$7, compute the SVD: $\eta$8 The unique minimizing logarithm for all $\eta$9: $Y^t \eta$0 The geodesic distance is

$Y^t \eta$1

These are analogous to the $Y^t \eta$2, $Y^t \eta$3 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 $Y^t \eta$4) plus Fréchet derivatives; cost per exponential is $Y^t \eta$5. For low-rank problems ($Y^t \eta$6), all computations occur in small dimension ($Y^t \eta$7). 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 $Y^t \eta$8 by gradient descent, where each step requires $Y^t \eta$9 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).