GEORCE: A Fast New Control Algorithm for Computing Geodesics

TL;DR

GEORCE reframes the geodesic boundary-value problem as a discrete optimal-control task, enabling a primal-dual update scheme that achieves global convergence and quadratic local convergence without requiring second-order derivatives. The method scales as $\mathcal{O}(Td^{3})$, with linear dependence on the number of grid points $T$, and extends from Riemannian to Finslerian geometry by incorporating a velocity-dependent fundamental tensor $G(x,v)$. Empirically, GEORCE delivers faster convergence and shorter geodesics than standard solvers across diverse manifolds, including information-geometric Fisher-Rao spaces and VAE-learned latent spaces. The findings suggest GEORCE as a robust, scalable tool for geodesic computation with broad applicability in geometry-driven data analysis and generative modeling.

Abstract

Computing geodesics for Riemannian manifolds is a difficult task that often relies on numerical approximations. However, these approximations tend to be either numerically unstable, have slow convergence, or scale poorly with manifold dimension and number of grid points. We introduce a new algorithm called GEORCE that computes geodesics via a transformation into a discrete control problem. We show that GEORCE has global convergence and quadratic local convergence. In addition, we show that it extends to Finsler manifolds. For both Finslerian and Riemannian manifolds, we thoroughly benchmark GEORCE against several alternative optimization algorithms and show empirically that it has a much faster and more accurate performance for a variety of manifolds, including key manifolds from information theory and manifolds that are learned using generative models.

Figures (34)

• Figure 1: Comparison between GEORCE and baseline methods for computing locally length minimizing geodesics between given point pairs on four different manifolds. All algorithms are terminated if the $\ell^{2}$-norm of the gradient of the discretized energy functional \ref{['eq:disc_const_energy']} is less than $10^{-4}$, or the number of iterations exceeds $100$. Note that sometimes the respective solutions overlap, making it difficult to distinguish them. Also note, that in the case of the torus the initial curve determines the direction around the hole so that the obtained geodesics follow this direction and are therefore clearly not the globally shortest geodesics between the given start and end point. Experimental details and benchmark data are in Appendix \ref{['ap:experiments']}.
• Figure 2: The estimates of bvp-solvers with different integration methods compared to GEORCE on $\mathbb{S}^{n}$ for $n=2,3,5,10$. The bvp-solvers minimize the squared error to the boundary condition using BFGSbroyden_bfgsfletcher_bfgsGoldfarb1970AFOshanno_bfgs. All methods are terminated if they take more than 24 hours or use more than 10 GB of memory on a CPU. Details on the experiment is found in Appendix \ref{['ap:hyper_parameters']}.
• Figure 3: The figure shows the $\ell^{2}$-norm of the gradient of eq. \ref{['eq:disc_const_energy']} for each iteration for GEORCE and alternative methods applied to construct geodesics for $100$ iterations for four different manifolds corresponding to Fig. \ref{['fig:synthetic_riemannian_geodesics']}. We see that the $\ell^{2}$-norm of the gradient of eq. \ref{['eq:disc_const_energy']} for GEORCE converges considerably faster than for alternative algorithms.
• Figure 4: The first two figures display the length and runtime of geodesics for $\mathbb{S}^{n}$ estimated using different algorithms from Table \ref{['tab:riemmannian_comparison_table']} and other methods in Appendix \ref{['ap:additional_experiments']}. The latter row shows the length and runtime, respectively for the $n$-Ellipsoid with half axes of $n$ equally spaced points between $0.5$ and $1.0$. The optimization solvers are described in Appendix \ref{['ap:hyper_parameters']}.
• Figure 5: The length of the estimated geodesics for the BFGS-algorithm, ADAM and GEORCE on a CPU. The methods were terminated if the $\ell^{2}$-norm of the gradient was less than $10^{-4}$ or after $1,000$ iterations. $\mathrm{E}(n)$ denotes an Ellipsoid of dimension $n$ with half axes of $n$ equally spaced points between $0.5$ and $1.0$, while $\mathcal{P}(n)$ denotes the space of $n \times n$ symmetric positive definite matrices. When the computational time was longer than $24$ hours, the value is set to $-$.

Theorems & Definitions (22)

• Proposition 1
• proof
• Proposition 2
• proof
• Proposition 3: Global convergence
• proof
• Proposition 4: Local quadratic convergence
• proof
• Proposition 6
• proof