---
title: Closed-Form Geodesics in Matrix Manifolds
url: https://www.emergentmind.com/topics/closed-form-geodesics
type: topic
---

# Closed-Form Geodesics in Matrix Manifolds

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 [2103.13327].

## 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 \)) [2103.13327].

## 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 \) [2103.13327].

## 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 [2103.13327].

## 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 [2103.13327].

## 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 [2103.13327].

## 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 [2103.13327].

Source: https://www.emergentmind.com/topics/closed-form-geodesics