GEORCE-FM: Joint Fréchet Mean & Geodesics

• GEORCE-FM is an efficient algorithm that jointly computes the Fréchet mean and associated geodesics by framing the problem as a discrete-time optimal control task in a local coordinate chart.
• The method circumvents repeated geodesic boundary value problems by employing closed-form joint updates, leading to global convergence and local quadratic convergence guarantees.
• Empirical results demonstrate that GEORCE-FM outperforms traditional methods in runtime and solution quality, with extensions to adaptive mini-batching and Finsler manifold settings.

The GEORCE-FM algorithm is an efficient and theoretically grounded method for simultaneously computing the Fréchet mean and associated geodesics on both Riemannian and Finsler manifolds. By framing the joint computation as a single discrete-time optimal control problem in a local coordinate chart, GEORCE-FM circumvents the traditional need to solve nested geodesic boundary value problems in each iteration of Fréchet mean estimation. The method accommodates adaptive extensions for large-scale data settings and provides rigorous guarantees of global convergence and local quadratic convergence. Empirical evidence demonstrates that GEORCE-FM dominates earlier manifold mean estimation techniques in both runtime and solution quality across classical and learned manifold models.

1. Fréchet Mean Estimation on Manifolds

The Fréchet mean $\mu$ for a set of points $a_1, \dots, a_N$ on a Riemannian manifold $(\mathcal{M}, g)$ is the minimizer of the sum of squared distances: $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$ with $\operatorname{dist}$ the geodesic distance. Computing $\operatorname{dist}(y, a_i)$ generally requires solving a geodesic boundary value problem, for example by integrating the ODE

$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$

or minimizing the energy functional

$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$

where $G$ is the metric tensor in a local chart. Most existing algorithms, therefore, embed a costly geodesic solver within each Fréchet mean update, leading to compounding computational expense.

2. Joint Energy Formulation and Optimization

GEORCE-FM addresses this computational bottleneck by jointly optimizing the mean location $y$ and all associated geodesics’ discretizations. In local coordinates, each geodesic $a_1, \dots, a_N$0 from $a_1, \dots, a_N$1 to $a_1, \dots, a_N$2 is discretized by a set of $a_1, \dots, a_N$3 points $a_1, \dots, a_N$4 with step vectors $a_1, \dots, a_N$5, subject to boundary conditions $a_1, \dots, a_N$6, $a_1, \dots, a_N$7, and state updates $a_1, \dots, a_N$8. The joint discrete energy minimized by GEORCE-FM is: $a_1, \dots, a_N$9 with Lemma 1 showing that as $(\mathcal{M}, g)$0, minimizing this energy over $(\mathcal{M}, g)$1 matches the original Fréchet mean problem. The algorithm thus eliminates the repeated inner geodesic solves typically required.

3. Algorithmic Workflow and Pseudocode

The GEORCE-FM algorithm operates in the following principal steps, all conducted within a single local chart (i.e., no need for chart transitions or global log/exponential map evaluations):

1. Initialization: Set $(\mathcal{M}, g)$2 (e.g., $(\mathcal{M}, g)$3), $(\mathcal{M}, g)$4, $(\mathcal{M}, g)$5.
2. Metric Evaluation: For all $(\mathcal{M}, g)$6, set $(\mathcal{M}, g)$7.
3. Gradient Computation: For all $(\mathcal{M}, g)$8, compute $(\mathcal{M}, g)$9.
4. Closed-Form Joint Update:
   • Compute

   $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$0

• Update the mean: $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$1.
• Compute dual variables and new controls for all $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$2:

  $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$3

  $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$4

1. Projected Line Search: Armijo backtracking line search is performed over the control updates to ensure sufficient energy decrease.
2. Mixing: Weighted averaging (with step size $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$5) of previous and proposed $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$6 ensures convergence and regularization.
3. Repeat until the decrease in sum energy $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$7 meets the tolerance.

This algorithmic procedure matches Algorithm 2 and its Riemannian update rules as given in (Rygaard et al., 6 Nov 2025). The Finslerian extension involves similar machinery but adapts $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$8 to depend on both position and direction (see Section 4).

4. Extensions: Finsler Manifolds and Adaptive Mini-Batching

Finsler Geometry: In Finsler manifolds, the metric is defined via a fundamental tensor $$\mu = \arg\min_{y \in \mathcal{M}} \sum_{i=1}^N \operatorname{dist}^2(y, a_i)$$9 that is homogeneous and strongly convex in the tangent vector $\operatorname{dist}$0. GEORCE-FM extends naturally to this setting by replacing the Riemannian $\operatorname{dist}$1 with $\operatorname{dist}$2, and updating the optimization objective accordingly: $\operatorname{dist}$3 The update rules, line-search procedure, and convergence proofs carry over nearly unchanged.

Adaptive (Stochastic) GEORCE-FM: For large $\operatorname{dist}$4, computational complexity is reduced by employing mini-batch stochastic updates: - At each outer iteration, a random subset $\operatorname{dist}$5 of size $\operatorname{dist}$6 is sampled. - GEORCE-FM is run for several steps within $\operatorname{dist}$7 yielding estimates $\operatorname{dist}$8. - The global mean iterate is formed using an exponential moving average in $\operatorname{dist}$9, with decaying step size $\operatorname{dist}(y, a_i)$0. - Under mild unbiasedness and correlation assumptions, the resulting sequence converges in expectation to a local minimizer (Proposition 7).

5. Theoretical Guarantees

Global Convergence: If the initial iterates satisfy the boundary constraints, and with the Armijo line search, GEORCE-FM generates a non-increasing sequence of objective values; energy is bounded below, so accumulation points are local minima (Prop. 5).

Local Quadratic Convergence: When the discrete energy $\operatorname{dist}(y, a_i)$1 is strongly convex in a neighborhood of a unique minimizer $\operatorname{dist}(y, a_i)$2, and the line search accepts full steps ($\operatorname{dist}(y, a_i)$3), quadratic convergence holds: $\operatorname{dist}(y, a_i)$4 for some $\operatorname{dist}(y, a_i)$5, by a Taylor expansion argument under smoothness and convexity assumptions.

Stochastic Convergence (Adaptive Mini-Batching): When mini-batch updates are unbiased and satisfy positivity (Eq. (26)), and each mini-batch uses a line search, convergence in expectation is assured (Prop. 7).

6. Computational Complexity and Empirical Performance

Theoretical Complexity

• Per Iteration: Each GEORCE-FM iteration is $\operatorname{dist}(y, a_i)$6 due to $\operatorname{dist}(y, a_i)$7 matrix inversions in $\operatorname{dist}(y, a_i)$8-dimensional charts.
• Comparison: This is comparable to standard Fréchet mean algorithms (which require $\operatorname{dist}(y, a_i)$9 geodesic solves per iteration at a similar cost) but benefits from amortized reuse of chart structures and avoids repeated boundary value problem solves.

Empirical Results

Benchmark	Baseline Iterations (ADAM/SGD)	GEORCE-FM Iterations	Acceleration (Runtime)	Accuracy
$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$0, $$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$1	$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$2–$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$3	$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$4–$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$5	$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$6–$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$7	Comparable
Fisher–Rao statistics	$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$8	$$\ddot\gamma^k + \Gamma^k_{ij}(\gamma)\dot\gamma^i\dot\gamma^j = 0, \qquad \gamma(0)=y,\, \gamma(1)=a_i$$9–$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$0	$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$1	Comparable
Learned VAE manifolds	Minutes	Tens of seconds	$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$2–$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$3	Comparable
Adaptive ($\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$4)	Full batch	$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$5 geodesics	$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$6	$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$7Few %

• GEORCE-FM converges rapidly on synthetic spheres, ellipsoids, and information geometry manifolds.
• On real-data latent spaces (e.g., MNIST, CelebA decoder manifolds), GEORCE-FM computes both geodesics and Fréchet mean on the order of seconds.
• In Finsler and Randers-type geometries, GEORCE-FM attains lower energy solutions and $\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$8–$\E_\gamma(y, a_i) = \frac12 \int_0^1 \dot\gamma(t)^\top G(\gamma(t)) \dot\gamma(t) dt$9 improved runtimes versus BFGS, ADAM, and Newton-type baselines.

7. Significance, Limitations, and Directions

GEORCE-FM’s principal innovation is the joint optimization over mean location and geodesic fields, facilitated by closed-form local updates within a single chart. This avoids compounding the computational and numerical overhead of embedded geodesic solvers typical of prior approaches. Its design allows extensibility to asymmetric and velocity-dependent metrics (Finsler), as well as natural stochastic mini-batching for large-scale data.

Limitations include the reliance on local charts (potentially problematic on highly nontrivial manifold topologies) and cubic scaling in intrinsic dimension $G$0 due to matrix inversions, although $G$1 is typically modest.

• A plausible implication is that further acceleration may be possible using approximate inversion strategies or manifold preconditioning for very high-dimensional problems.
• The uniform chart assumption could be relaxed for more general manifold geometries via atlas-based or manifold learning techniques.

Recent empirical benchmarking on classical and learned manifolds supports the algorithmic dominance of GEORCE-FM in both runtime and accuracy for Fréchet mean estimation, with strong theoretical backing for its global and local convergence properties (Rygaard et al., 6 Nov 2025).