Parallelised Differentiable Straightest Geodesics

• Parallelised differentiable straightest geodesics are GPU-accelerated methods that compute geodesic curves (with zero covariant acceleration) on both smooth and discrete manifolds.
• They employ batch ODE solvers, symbolic differentiation, and specialized GPU kernels to efficiently implement the exponential map and gradient updates across meshes and grids.
• Applications range from geometry processing and shape analysis to generative modeling and optimization, delivering significant speedups and high accuracy in complex computational tasks.

Parallelised differentiable straightest geodesics provide a foundation for scalable, GPU-accelerated computation of geodesics and their derivatives on both smooth and discretised manifolds, with applications spanning geometry processing, shape analysis, generative modeling, and variational optimization. Central to these developments is the ability to implement the exponential map and related operations in a manner that is both batch-parallel and amenable to gradient-based optimization, regardless of whether the underlying geometry is defined analytically or via discrete structures such as meshes or grids.

1. Mathematical Foundations

The straightest geodesic is defined as a curve $\gamma(t)$ on a Riemannian manifold $(M,g)$ such that the covariant acceleration vanishes: $\nabla_{\dot\gamma}\dot\gamma = 0$ In coordinates,

$$\ddot{x}^k + \Gamma^k_{ij}(x)\,\dot{x}^i\,\dot{x}^j = 0$$

where $\Gamma^k_{ij}$ are the Christoffel symbols derived from the metric $g$. The exponential map $\mathrm{Exp}_p:T_pM \rightarrow M$ maps a tangent vector $v$ at base point $p$ to the endpoint $\gamma(1)$ of the geodesic starting at $p$ with initial velocity $v$. On smooth manifolds, these ODEs are numerically integrated using symbolic differentiation and ODE solvers (Kühnel et al., 2017).

On triangular meshes, straightest geodesics are realized as piecewise-linear curves that preserve the incident angle at each crossing—using barycentric coordinates within faces and explicit angle balancing as in the SG-Step routine. At a point $p=(f,b)$ (face and barycentric coordinates), geodesic iteration proceeds by either transporting within a face, across an edge, or through a vertex, maintaining the angle condition described by Polthier & Preuß. The discrete exponential map is then the endpoint after following a tangent vector of prescribed length, normalised in the tangent fan (Verninas et al., 16 Mar 2026).

2. Parallel GPU Algorithms

Batch-parallelisation is achieved by assigning each geodesic or grid-point update to an independent thread. For smooth manifolds, symbolic computation frameworks like Theano allow vectorization over batches by lifting point and velocity inputs to matrices; all tensor contractions and ODE steps are broadcast, yielding lockstep execution on GPU. Empirically, this enables marching $10^4$–$10^5$ geodesics in parallel, with tens of thousands of dimensionality processed in seconds (Kühnel et al., 2017).

For discrete meshes, the core algorithm executes the SG-Step for each geodesic independently in parallel, utilizing GPU kernels to iterate until the cumulative geodesic length is traversed. Each thread maintains local state and updates geometry, with mesh connectivity and face normals preloaded into GPU memory. As mesh size and batch count grow, complexity scales linearly, and benchmarks reach $10^5$ geodesics on $10^6$-face meshes in under $10$ ms on an RTX 4090 (Verninas et al., 16 Mar 2026).

On Cartesian grids with (possibly anisotropic) Randers-Finsler metrics, parallel Eikonal solvers use column-row fast sweeping: each sweep processes a line (row or column) in parallel, using only upwind (already swept) neighbors for updates, with optimal data flow on GPU. Each grid node update involves either one- or two-point upwind updates with closed-form solutions (Gahtan et al., 7 Feb 2026).

3. Differentiability and Gradient Computation

Native automatic differentiation fails for geodesics on discrete meshes due to non-differentiability at face transitions and combinatorial update steps. Two primary schemes are introduced (Verninas et al., 16 Mar 2026):

• Extrinsic Proxy (EP): The output point is approximated by an extrinsic rigid transformation proxy $\phi(p,v)=R_{\text{fix}}(p+v) + t_{\text{fix}}$, where $R_{\text{fix}}$ aligns tangent bases and $t_{\text{fix}}$ is a translation. The Jacobian $\partial\phi/\partial v = R_{\text{fix}}$ is used in backpropagation; derivatives w.r.t. $p$ are set to zero.
• Geodesic Finite Differences (GFD): Forward finite differences are computed for both $p$ and $v$ in local tangent frames. Small perturbations are parallel transported, Exp maps are recomputed, and the resulting differences yield $J_{Exp}^v$ and $J_{Exp}^p$. This supports full optimization w.r.t. both base and tangent.

On regular grids, differentiation is handled via implicit differentiation. After forward solution of the Eikonal, the backward pass solves the adjoint equation by reverse-time back-substitution, using only the final arrival-time field and active-stencil map—resulting in $O(N^2)$ memory and computational cost. All partial derivatives with respect to parameters (metric and drift) are given in closed-form for each update type (Gahtan et al., 7 Feb 2026).

4. Applications in Learning and Optimization

Differentiable, parallelised geodesics enable new classes of geometric ML and optimization pipelines:

• Adaptive Geodesic Convolution (AGC): Exp maps provide differentiable, learnable patch radii on meshes, improving segmentation accuracy over fixed-radius GCNNs. Filters are constructed by Exp mapping a local grid of directions at each vertex, interpolating features, and learning filter weights and patch size (Verninas et al., 16 Mar 2026).
• MeshFlow (Flow Matching): Geodesic-based generative models use stationary vector fields and iterative Exp mapping to sample paths; matching is performed using optimal transport couplings and geodesic losses, yielding improved likelihood and metric precision along with substantial GPU memory and runtime savings over projection-augmented Riemannian flow matching (Verninas et al., 16 Mar 2026).
• Second-Order Mesh Optimization (Mesh-LBFGS): Centroidal Voronoi Tessellation and similar manifold optimization tasks exploit parallel Exp for manifold-constrained LBFGS with vector transport and adjoint recursion; this achieves significant reductions in iteration count and function calls relative to first-order (Lloyd) methods (Verninas et al., 16 Mar 2026).
• Inverse Problem Solvers: Differentiable Eikonal solvers for Randers-Finsler metrics support metric and drift field optimization from sparse supervision (e.g., wildfire perimeter matching), with provable gradient accuracy and stability (Gahtan et al., 7 Feb 2026).

5. Implementation and Practical Considerations

For smooth manifolds, Theano (and similarly PyTorch/JAX) provides the primitives needed for batch ODE solvers, automatic Christoffel computation, and gradient through the entire geodesic computation graph. Batch support is a matter of tensor shape lifting; parallel GPU execution is activated via function compilation in FAST_RUN mode (Kühnel et al., 2017).

On meshes, the digeo library implements the MeshExponential interface with CUDA acceleration. The user chooses between EP and GFD schemes depending on task requirements—EP is preferred for fixed-point convolutions, GFD for full manifold optimization. Preprocessing steps include face normal and connectivity precomputation. Mixed precision (float32) delivers $<10^{-4}$ error. Integration into PyTorch pipelines relies on batching all calls to the Exp map and, for GFD, batching all perturbation directions (Verninas et al., 16 Mar 2026).

In Eikonal settings, metric feasibility is maintained through eigenvalue and drift clamping, ensuring positive-definite metrics and subunit drift norm (Randers bound). These projections are a.e. differentiable and included within the computational graph (Gahtan et al., 7 Feb 2026).

6. Accuracy, Stability, and Benchmark Performance

On spherical and toroidal test surfaces, straightest geodesic methods match or exceed the accuracy of geometry-central/potpourri3d up to $2 \times 10^{-5}$ (float32), with errors diminishing as mesh resolution increases (Verninas et al., 16 Mar 2026). Sphere and Torus canonical shapes yield errors $<10^{-3}$. For mesh geodesics, the EP backward scheme offers $\sim 4\times$ speedup over GFD, and GFD recovers both Jacobians with cosine similarity $\sim0.99$ to analytical solutions on the sphere.

Eikonal solver differentiation accuracy is machine-precision ($<10^{-6}$ error vs. central finite differences) away from nondifferentiable stencil boundaries, and gradient fields exhibit Lipschitz-type stability under parameter perturbations. Measured, a $1\%$ change in metric yields $0.48\%$ change in gradient norm. Optimization remains stable even when boundaries are crossed (Gahtan et al., 7 Feb 2026).

Remarkable parallel performance is observed: $2,500$ (5,000-dim) geodesic landmarks matched on CPU in seconds, $20,000$ (40,000-dim) on GPU in under a second (Kühnel et al., 2017); mesh methods compute $10^5$ geodesics on $10^6$-face meshes in sub-10ms timeframes, scaling linearly in mesh and batch size (Verninas et al., 16 Mar 2026).

---

Key References:

• (Kühnel et al., 2017) Kühnel et al., "Differential geometry and stochastic dynamics with deep learning numerics"
• (Verninas et al., 16 Mar 2026) "Parallelised Differentiable Straightest Geodesics for 3D Meshes"
• (Gahtan et al., 7 Feb 2026) "Wildfire Simulation with Differentiable Randers-Finsler Eikonal Solvers"