# Finite elements approximation of the boundary value problems of geodesics

Gaël Le Ruz, Damiano Lombardi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Gaël Le Ruz, Damiano Lombardi. Finite elements approximation of the boundary value problems of geodesics. 2025. hal-05354459

HAL Id: hal-05354459

https://inria.hal.science/hal-05354459v1

Preprint submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Finite elements approximation of the boundary value problems of geodesics.

Gaël, Le Ruz

LJLL, CNRS, Inria, Sorbonne Université, 75005 Paris

e-mail: gael.le-ruz@inria.fr

web page: www.lrzgael.github.io.

Damiano, Lombardi

LJLL, CNRS, Inria, Sorbonne Université, 75005 Paris

e-mail: damiano.lombardi@inria.fr

November 7, 2025

## Abstract

In this contribution we investigate the numerical approximation of the boundary value problems of geodesics, i.e. the log operation, on Riemannian manifolds. In particular, by leveraging the variational formulation of the geodesics problem, we propose a finite elements discretisation in time. We also investigate the numerical approximation of the Hessian of the squared distance with a sensitivity analysis on the first problem.

## 1 Introduction

In the following, let $\mathcal{M}$ denote a Riemannian manifold of dimension $d$. A minimizing curve on $\mathcal{M}$ is a curve that realizes the shortest path between two points. In local coordinates, such a curve $x$ can be characterized as a minimizer of the geodesic energy

$$ \mathcal{E}(x) = \frac{1}{2} \int_{0}^{1} \dot{x}(t)^{T} g(x(t)) \dot{x}(t) dt, $$

where $g$ denotes the Riemannian metric on $\mathcal{M}$. The associated Euler-Lagrange equations yield the well-known geodesic equation:

$$ \frac{d^2 x_k(t)}{dt^2} + \sum_{i,j=1}^{d} \Gamma_{ij}^{k}(x(t)) \frac{dx_i(t)}{dt} \frac{dx_j(t)}{dt} = 0, \quad k = 1, 2, \dots, d, \quad t \in [0, 1], $$

where $\Gamma_{ij}^k$ denote the Christoffel symbols of the metric $g$.

A curve satisfying the geodesic equation is not necessarily a global minimizer of $\mathcal{E}$, since it corresponds merely to a stationary point of the functional. Nevertheless, such a curve may
be interpreted as a *straight curve*, in the sense that its covariant acceleration vanishes. If $\nabla$ denotes the Levi-Civita connection — the unique torsion-free connection compatible with the metric $g$ — the geodesic equation can be written in the compact form

$$ \nabla_{\dot{x}(t)} \dot{x}(t) = 0, \quad t \in [0, 1]. $$

Thus, a curve $x$ satisfying this equation, although not necessarily minimizing the energy functional, is a curve with zero covariant acceleration. Such a curve is referred to as an *affine geodesic*, or simply a *geodesic*.

In the numerical approximation of geodesics, two distinct problems arise. The first is the *initial value problem* (IVP), which consists in finding a geodesic $x(t)$ with prescribed initial position and velocity, namely

$$ x(0) = x_0, \quad \dot{x}(0) = v_0. $$

The mapping

$$ \mathrm{Exp}_{x_0}(v_0) = x(1) $$

is called the *exponential map* at $x_0$. It is defined on the entire tangent space $T_{x_0}M$, but is bijective only within a ball centered at 0 of radius $\mathrm{inj}(x_0)$, the *injectivity radius* of $M$ at $x_0$. Within this ball, the notions of geodesics and minimizing curves coincide, and the energy functional admits a unique minimizer.

The second is the *boundary value problem* (BVP), which consists in finding a geodesic $x$ joining two prescribed points,

$$ x(0) = x_0, \quad x(1) = x_1. $$

The mapping

$$ \mathrm{Log}_{x_0}(x_1) = \dot{x}(0) $$

is called the *logarithmic map* at $x_0$. It is the local inverse of the exponential map and is therefore uniquely defined only when $x_1$ lies within the injectivity domain of the exponential map at $x_0$, i.e., when the geodesic connecting $x_0$ and $x_1$ is unique and minimizing.

When working on a surface discretized by a triangular mesh, the most natural formulation to consider is the BVP, that is, the search for the shortest path between two points on the surface. The MMP algorithm [16] and the Chen-Han algorithm [4] provide exact solutions to this problem, albeit at a high computational cost. Several approximations [20] and adaptations for large-scale data sets [23] have since been proposed. Alternatively, the Riemannian distance between two points — corresponding to the length of the minimizing curve — can be computed by solving an eikonal equation using the fast marching method [11, 15]. In this approach, geodesics are recovered by tracing the characteristics of the resulting distance function. More recently, the *heat method* [6, 7], based on Varadhan's asymptotic formula [22], has been introduced as an efficient alternative for distance computation. In contrast, the IVP involves geodesic tracing rather than distance computation. In this framework, solving the IVP typically requires dedicated geodesic tracing methods [5].

In the continuous setting, however, when an atlas — a collection of local coordinate charts covering the entire manifold — and the associated Riemannian metric $g$ are available, the
situation is reversed: the IVP is comparatively easier to handle, since it amounts to integrating the geodesic equation from prescribed initial conditions. It is typically approximated using an explicit Euler scheme or, more accurately, a symplectic integrator applied to the geodesic equation written in its Hamiltonian form. In contrast, solving the BVP — finding the geodesic connecting two fixed points — remains a challenging task. A class of shooting methods has been proposed to address this problem, which consists in minimizing the residual of the following objective function:

$$ \mathcal{J}(v) = \|\mathrm{Exp}_{x_0}(v) - x_1\|. $$

Such approaches have been proposed in [3, 10, 12], including specialized formulations for particular manifolds, such as the Stiefel manifold, where explicit Jacobians can be derived to enable efficient Newton-type shooting methods [21]. However, the minimization of $\mathcal{J}$ may lead to convergence toward local minima that do not correspond to the desired geodesic.

In contrast, our approach consists in solving the BVP directly from the geodesic equation using a fixed-point iteration that enforces the prescribed boundary conditions. A similar strategy was employed in [9], where a finite-difference approximation of the geodesic equation was used. In this work, we propose and justify a finite-element discretization in time of the geodesic equation. This method requires an initial curve that is subsequently “straightened” through the fixed-point iteration. Accordingly, we adopt the term path-straightening method, as introduced in [19] for the computation of geodesics in shape spaces.

In section 2, we propose and justify a fixed-point iteration formulated at the continuous level. We then introduce a finite-element discretization in time of the geodesic equation and study the behavior of the fixed-point iteration in the discrete setting. In section 3, we show how the proposed framework naturally extends to the approximation of the Hessian of the Riemannian squared distance—a quantity that arises in many optimization problems defined on manifolds [8, 17, 18]. Finally, in section 4, we present numerical experiments on the unit sphere embedded in $\mathbb{R}^3$, illustrating the convergence results established in section 2 and section 3.

## 2 The boundary value problem: computing the minimizing path between two points

Let $\mathcal{M}$ be a Riemannian manifold of dimension $d$. Consider $(\mathcal{U}, \varphi)$ a chart of $\mathcal{M}$, that is, $\mathcal{U}$ is an open subset of $\mathcal{M}$ and $\varphi: \mathcal{U} \to \varphi(\mathcal{U}) \subset \mathbb{R}^d$ is a homeomorphism. The elements of $\varphi(\mathcal{U})$ are referred to as the coordinates of points in $\mathcal{U}$. Let $\mathcal{B}_{\mathcal{M}}$ denote a ball in $\varphi(\mathcal{U})$. We define the set of curves on $\mathcal{B}_{\mathcal{M}}$ as $V = H^1([0, 1], \mathcal{B}_{\mathcal{M}})^{.1}$. Let $x_0, x_1 \in \mathcal{B}_{\mathcal{M}}$, a curve $x \in V$ with fixed endpoints $x(0) = x_0$ and $x(1) = x_1$ can be decomposed as

$$ x(t) = (1-t)x_0 + tx_1 + u(t), $$

where $u \in V_0 = H_0^1([0, 1], \mathcal{Y})$, and the set $\mathcal{Y}$ is defined by

$$ \mathcal{Y} = \{y \mid (1-t)x_0 + tx_1 + y \in \mathcal{B}_{\mathcal{M}}, t \in [0, 1]\}. $$

<footnote n="1">¹Since $H^1$ in one dimension is continuously embedded in $\mathcal{C}^0$, this definition is equivalent to considering the set of continuous functions.</footnote>
The velocity of the curve $x(t)$ is then given by

$$ \dot{x}(t) = (x_1 - x_0) + \dot{u}(t). $$

For the sake of compactness in the notation, we denote $\Delta = x_1 - x_0 \in \mathbb{R}^d$ and $\ell(t) = (1 - t)x_0 + tx_1$. Let $g(x(t)) \in \mathbb{R}^{d \times d}$ be the metric tensor at $x(t)$, which is a symmetric positive definite matrix. From now on, we denote by $S^{d \times d}$ the cone of symmetric positive definite matrices. A minimizing curve can be defined as a minimum of the geodesic energy $\mathcal{E}: V \to \mathbb{R}^+$:

$$ \mathcal{E}(x) = \frac{1}{2} \int_{0}^{1} \dot{x}(t)^{T} g(x(t)) \dot{x}(t) dt. \qquad (1) $$

This equation can be cast as a variational problem for $u \in V_0 = H_0^1([0, 1], \mathcal{Y})$ by simply substituting the expression of $x$:

$$ \mathcal{E}(u) = \frac{1}{2} \int_{0}^{1} (\Delta + \dot{u}(t))^{T} g(\ell(t) + u(t)) (\Delta + \dot{u}(t)) dt. \qquad (2) $$

With a slight abuse of notation, we henceforth identify the geodesic energy with $\mathcal{E}: \mathcal{Y} \to \mathbb{R}^+$. Let us remark that an essentially bounded metric tensor, $g \in L^{\infty}(\mathcal{B}_{\mathcal{M}}, S^{d \times d})$, is sufficient to ensure that the energy is finite.

We consider the problem of finding the stationary points of the geodesic energy. To this end, we write the Euler-Lagrange equations of the functional defined in Eq. (2). For any $1 \le l \le d$, it holds that:

$$ \frac{d}{dt} \left( \sum_{j=1}^{d} g_{lj}(\ell(t) + u(t)) (\Delta_j + \dot{u}_j(t)) \right) - \frac{1}{2} \left( \sum_{i,j=1}^{d} \partial_{x_i} g_{ij}(\ell(t) + u(t)) (\Delta_i + \dot{u}_i(t)) (\Delta_j + \dot{u}_j(t)) \right) = 0. $$

This is in fact another form of the geodesic equation. Let us write the weak formulation of this equation. Let $v \in V_0$ and define $\mathcal{F}: V_0 \times V_0 \to \mathbb{R}$ by

$$ \mathcal{F}(u, v) = \int_{0}^{1} \dot{v}(t)^{T} g(\ell(t) + u(t)) (\Delta + \dot{u}(t)) dt \qquad (3) \\ + \frac{1}{2} \int_{0}^{1} (\Delta + \dot{u}(t))^{T} [\nabla_{x} g(\ell(t) + u(t)) \cdot v(t)] (\Delta + \dot{u}(t)) dt. $$

A stationary point $u_* \in V_0$ is such that, for all $v \in V_0$, $\mathcal{F}(u_*, v) = 0$. In the following, we assume that the endpoints $x_0$ and $x_1$ are chosen such that the corresponding stationary point of $\mathcal{E}$ is also a global minimum. In other words, the geodesic curve defined by the above equation coincides with the minimizing curve.
## 2.1 A fixed-point iteration

As the problem of finding a stationary point is utterly non-linear, we propose a fixed-point iteration to obtain an approximation. We define the fixed point at the continuous level. Let $u^{(0)} \in V_0$. For $k \in \mathbb{N}^*$, and with a slight abuse of notation, we define $\mathcal{F}(u^{(k-1)}) : V_0 \times V_0 \to \mathbb{R}$ as

$$ 
\begin{aligned}
\mathcal{F}(u^{(k-1)})(w, v) = & \int_0^1 \dot{v}(t)^T g(\ell(t) + u^{(k-1)}(t)) (\Delta + \dot{w}(t)) dt \quad (4) \\
& + \frac{1}{2} \int_0^1 \Delta^T [\nabla_x g(\ell(t) + u^{(k-1)}(t)) \cdot v(t)] \Delta dt \\
& + \frac{1}{2} \int_0^1 (\dot{u}^{(k-1)}(t))^T [\nabla_x g(\ell(t) + u^{(k-1)}(t)) \cdot v(t)] \dot{w}(t) dt \\
& + \int_0^1 \Delta^T [\nabla_x g(\ell(t) + u^{(k-1)}(t)) \cdot v(t)] \dot{w}(t) dt.
\end{aligned}
$$

The approximation $u^{(k)}$ is implicitly defined as the element of $V_0$ satisfying

$$ \mathcal{F}(u^{(k-1)})(u^{(k)}, v) = 0, \quad \text{for all } v \in V_0. $$

By construction, any fixed point of the above equation corresponds to a stationary point of the geodesic energy. It follows from the definition of $\mathcal{F}(u^{(k-1)})$ in (4) that the problem for $u^{(k)}$ is linear. In particular, we can define a bilinear form $A(u^{(k-1)}) : V_0 \times V_0 \to \mathbb{R}$ by

$$ 
\begin{aligned}
A(u^{(k-1)})(w, v) = & \int_0^1 \dot{v}(t)^T g(\ell(t) + u^{(k-1)}(t)) \dot{w}(t) dt \quad (5) \\
& + \int_0^1 \left( \Delta + \frac{1}{2} \dot{u}^{(k-1)}(t) \right)^T [\nabla_x g(\ell(t) + u^{(k-1)}(t)) \cdot v(t)] \dot{w}(t) dt,
\end{aligned}
$$

and a linear form $b(u^{(k-1)}) : V_0 \to \mathbb{R}$ by

$$ 
\begin{aligned}
b(u^{(k-1)})(v) = & - \int_0^1 \dot{v}(t)^T g(\ell(t) + u^{(k-1)}(t)) \Delta dt \quad (6) \\
& - \frac{1}{2} \int_0^1 \Delta^T [\nabla_x g(\ell(t) + u^{(k-1)}(t)) \cdot v(t)] \Delta dt,
\end{aligned}
$$

such that the k-th step of the fixed-point iteration consists in finding $u^{(k)} \in V_0$ such that,

$$ A(u^{(k-1)})(u^{(k)}, v) = b(u^{(k-1)})(v), \quad \forall v \in V_0. \quad (7) $$

## 2.2 Well-posedness of the fixed-point iteration.

Let us prove that the problem defined in (7) is well-posed, thereby ensuring that the fixed-point iteration is well-defined. We first present three lemmas that, under suitable assumptions, show that the linear functional on the right-hand side of (7) is continuous and that the bilinear form on the left-hand side is continuous and coercive.
**Lemma 1.** Let $g \in W^{1,\infty}(\mathcal{B}_M, \mathcal{S}^{d \times d})$. For any $u \in V_0$, the linear form $b(u) : V_0 \to \mathbb{R}$ is continuous.

*Proof.* We consider separately the different terms appearing in (6); the definition of $b(u)$. We denote by $\lambda_1(t)$ the largest eigenvalue of $g(\ell(t) + u(t))$. By the regularity assumption on the metric tensor, this quantity is bounded. Let $M(u) = \sup_{\tau \in [0,1]} \lambda_1(\tau)$. The first term of the linear form satisfies

$$ \left| \int_0^1 \dot{v}(t)^T g(\ell(t) + u(t)) \Delta t \right| \leq |\Delta| M(u) \int_0^1 |\dot{v}(t)| dt, $$

which gives, by the virtue of the Cauchy–Schwarz inequality,

$$ \left| \int_0^1 \dot{v}(t)^T g(\ell(t) + u(t)) \Delta t \right| \leq |\Delta| M(u) \|\dot{v}\|_{L^2([0,1])^d} \leq |\Delta| M(u) \|v\|_{V_0}. $$

For the second term, we denote

$$ I_2 = \frac{1}{2} \int_0^1 \sum_{i,j=1}^d \sum_{l=1}^d \partial_{x_l} g_{ij}(\ell(t) + u(t)) v_l(t) \Delta_i \Delta_j dt. $$

By hypothesis, the derivatives of the metric tensor are bounded. In particular, we have

$$ |I_2| \leq \frac{1}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \int_0^1 \sum_{i=1}^d |\Delta_i| \sum_{j=1}^d |\Delta_j| \sum_{l=1}^d |v_l(t)| dt, $$

where

$$ \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} = \sup_{\substack{1 \leq i,j,l \leq d \\ x \in \mathcal{B}_M}} |\partial_{x_l} g_{ij}(x)|. $$

Each sum can be bounded using the Cauchy–Schwarz inequality for vectors. It follows that

$$ |I_2| \leq \frac{d^3}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} |\Delta|^2 \int_0^1 |v(t)| dt. $$

Combining the estimates, we obtain

$$ |b(u)(v)| \leq \left( |\Delta| M(u) + \frac{d^3}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} |\Delta|^2 \right) \|v\|_{V_0}. $$

In the following, the continuity constant of $b(u)$ is denoted by $C_b(u)$ and defined as

$$ C_b(u) = |\Delta| M(u) + \frac{d^3}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} |\Delta|^2. \qquad (8) $$

<figure><img src="image_3.png" /></figure>

Let us now proceed with the continuity of the bilinear form $A(u)$.
**Lemma 2.** Let $g \in W^{1,\infty}(\mathcal{B}_M, \mathcal{S}^{d \times d})$ and $u$ be a continuous curve with maximal finite speed, we denote $\sup_{\tau \in [0,1]} |\Delta + \dot{u}(\tau)| = c(u)$. The bilinear form $A(u) : V_0 \times V_0 \to \mathbb{R}$ is continuous.

*Proof.* Similarly to what was done for the linear form, we consider the terms in (5) separately. The first one is:

$$I_1 = \int_0^1 \dot{v}(t)^T g(\ell(t) + u(t)) \dot{w}(t) dt.$$

As previously, we denote by $\lambda_1(t)$ the largest eigenvalue of the metric tensor and set $M(u) = \sup_{\tau \in [0,1]} \lambda_1(\tau)$. It holds:

$$|I_1| \le M(u) \|\dot{w}\|_{L^2([0,1])^d} \|\dot{v}\|_{L^2([0,1])^d} \le M(u) \|w\|_{V_0} \|v\|_{V_0}.$$

The second term is:

$$I_2 = \int_0^1 \sum_{i,j=1}^d \left( \Delta_i + \frac{1}{2} \dot{u}_i(t) \right) \left[ \sum_{l=1}^d \partial_{x_l} g_{ij}(\ell(t) + u(t)) v_l(t) \right] \dot{w}_j(t) dt.$$

As before, the derivatives of the metric tensor are bounded, and we have

$$|I_2| \le \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \int_0^1 \sum_{i=1}^d \left| \Delta_i + \frac{1}{2} \dot{u}_i(t) \right| \sum_{j=1}^d \left| \dot{w}_j(t) \right| \sum_{l=1}^d |v_l(t)| dt.$$

The first sum is bounded by the Euclidean norm at time $t$. By hypothesis, this is finite,

$$\sum_{i=1}^d \left| \Delta_i + \frac{1}{2} \dot{u}_i(t) \right| \le \frac{1}{2} \sum_{i=1}^d \left| \Delta_i + (\Delta_i + \dot{u}_i(t)) \right| \le \frac{d^{\frac{1}{2}}}{2} (|\Delta| + c(u)).$$

By the Cauchy-Schwarz inequality, we get

$$|I_2| \le \frac{d^{\frac{3}{2}}}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} (|\Delta| + c(u)) \|v\|_{V_0} \|w\|_{V_0}.$$

Putting it all together, we obtain

$$|A(u)(w, v)| \le \left( M(u) + \frac{d^{\frac{3}{2}}}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} (|\Delta| + c(u)) \right) \|v\|_{V_0} \|w\|_{V_0}.$$

□

We now prove the coercivity of the bilinear form $A(u)$ under suitable assumptions.

**Lemma 3.** Let the metric tensor $g \in W^{1,\infty}(\mathcal{B}_M, \mathcal{S}^{d \times d})$, and let $\lambda_d(x)$ denote the smallest eigenvalue of $g(x)$ at any point $x \in \mathcal{M}$. We define $\lambda_d^* = \inf_{x \in \mathcal{B}_M} \lambda_d(x)$. Let $u \in V_0$ be a regular curve with maximal finite speed, we denote $\sup_{\tau \in [0,1]} |\Delta + \dot{u}(\tau)| = c(u)$, and let $C_P > 0$ be the Poincaré constant for $V_0$. Moreover, we assume that

$$(|\Delta| + c(u)) \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} < \frac{4 \lambda_d^* d^{-\frac{3}{2}}}{(1 + C_P^2)}.$$

Then, the bilinear form $A(u) : V_0 \times V_0 \to \mathbb{R}$ is coercive.
Proof. We consider separately the different terms in the sum. The first one, setting $w = v$, is

$$I_1 = \int_0^1 \dot{v}(t)^T g(\ell(t) + u(t)) \dot{v}(t) dt.$$

Let $C_P > 0$ be the Poincaré constant; it holds that

$$\|v\|_{L^2([0,1])^d} \le C_P \|\dot{v}\|_{L^2([0,1])^d}.$$

Hence,

$$\|v\|_{V_0} \le (1 + C_P^2)^{1/2} \|\dot{v}\|_{L^2([0,1])^d}.$$

Since the metric tensor is positive definite, we have $\lambda_d^* > 0$. It follows that

$$|I_1| \ge \lambda_d^* \|\dot{v}\|_{L^2([0,1])^d}^2 \ge \frac{\lambda_d^*}{(1 + C_P^2)} \|v\|_{V_0}^2.$$

We now consider the second and third terms in the sum. We begin with the second one, setting $w = v$:

$$I_2 = \int_0^1 \sum_{i,j=1}^d \left( \Delta_i + \frac{\dot{u}_j(t)}{2} \right) \left[ \sum_{l=1}^d \partial_{x_l} g_{ij}(\ell(t) + u(t)) v_l(t) \right] \dot{v}_j(t) dt.$$

We can bound the absolute value of $I_2$ as

$$|I_2| \le \frac{1}{2} (|\Delta| + c(u)) d^{\frac{3}{2}} \|\nabla g\|_{L^{\infty}(B_M, S^{d \times d})} \int_0^1 |v(t)| |\dot{v}(t)| dt.$$

Using the inequality between arithmetic and geometric means, we obtain

$$|I_2| \le \frac{1}{2} (|\Delta| + c(u)) d^{\frac{3}{2}} \|\nabla_x g\|_{L^{\infty}(B_M)} \frac{1}{2} \int_0^1 (|v(t)|^2 + |\dot{v}(t)|^2) dt.$$

Hence, since $I_2 \ge -|I_2|$, we get

$$A(u)(v, v) \ge \left( \frac{\lambda_d^*}{(1 + C_P^2)} - \frac{1}{4} (|\Delta| + c(u)) d^{\frac{3}{2}} \|\nabla_x g\|_{L^\infty(B_M)} \right) \|v\|_{V_0}^2.$$

Given the hypothesis on the metric tensor and on the curve at the previous iteration of the fixed-point scheme, this concludes the proof. In the following, we denote by $c_A(u)$ the coercivity constant of $A(u)$, defined as

$$c_A(u) = \frac{\lambda_d^*}{(1 + C_P^2)} - \frac{1}{4} (|\Delta| + c(u)) d^{\frac{3}{2}} \|\nabla_x g\|_{L^\infty(B_M)} \quad (9)$$

<figure><img src="image_4.png" /></figure>
**Proposition 1.** Let the metric tensor be $g \in W^{1,\infty}(\mathcal{B}_M, \mathcal{S}^{d \times d})$. Let the curve at iteration $(k-1)$ be such that $\sup_{\tau \in [0,1]} |\dot{u}^{(k-1)}(\tau)| = c(u^{(k-1)}) < \infty$. Moreover, assume that the metric tensor satisfies the condition

$$ (|\Delta| + c(u^{(k-1)})) \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} < \frac{4 \lambda_d^* d^{-\frac{3}{2}}}{(1 + C_P^2)}, $$

where $\lambda_d^*$ denotes the minimal eigenvalue of $g$ on $\mathcal{B}_M$, and $C_P > 0$ is the Poincaré constant associated with $V_0$. Then, the $k$-th fixed-point iteration is well-posed and its solution is unique.

*Proof.* Lemma 1 establishes that the linear form $b(u^{(k-1)})$ is continuous. Moreover, Lemmas 2 and 3 show that, under the above hypotheses, the bilinear form $A(u^{(k-1)})$ is continuous and coercive. Hence, by the virtue of the Lax-Milgram theorem, the variational problem (7) is well-posed and admits a unique solution $u^{(k)} \in V_0$. $\square$

## 2.3 Convergence of the fixed-point

We prove in this section that the fixed-point iteration converges under suitable assumptions. In particular, we show that appropriate regularity conditions on the metric ensure that the mapping defined by one fixed-point iteration is a contraction.

Recall that the $k$-th fixed-point iteration can be written as the problem of finding $u^{(k)} \in V_0$ such that, for all $v \in V_0$,

$$ A(u^{(k-1)})(u^{(k)}, v) = b(u^{(k-1)})(v). $$

Note that, while solving a single iteration by fixing $u^{(k-1)}$ yields a well-posed linear problem, the overall fixed-point iteration remains non-linear, as the dependence on $u^{(k-1)}$ is itself non-linear. The solution operator is defined as $T: V_0 \to V_0$ by

$$ u^{(k-1)} \mapsto u^{(k)} = T(u^{(k-1)}). $$

Let $u_1, u_2 \in V_0$, and define $w_1 = T(u_1)$ and $w_2 = T(u_2)$. For all $v \in V_0$, it holds that

$$ A(u_1)(w_1, v) = b(u_1)(v), $$

$$ A(u_2)(w_2, v) = b(u_2)(v). $$

Taking the difference of the two equations and rearranging terms, we obtain

$$ A(u_1)(w_1 - w_2, v) = (b(u_1) - b(u_2))(v) + (A(u_2) - A(u_1))(w_2, v). \quad (10) $$

This equation is the starting point of our proof. We now state the Lipschitz continuity properties of $A$ and $b$, which will make it possible to bound the two terms on the right-hand side of (10).

**Lemma 4.** Let $g \in C^{1,1}(\mathcal{B}_M, \mathcal{S}^{d \times d})$, that is, the gradient of the metric is uniformly Lipschitz continuous on $\mathcal{B}_M$, with constant $L_g > 0$. Let $u_1, u_2 \in V_0$ be continuous curves. Then, for all $v \in V_0$,

$$ |b(u_1)(v) - b(u_2)(v)| \le \left( d \|\nabla g\|_{L^\infty(\mathcal{B}_M, \mathcal{S}^{d \times d})} |\Delta| + \frac{d^3}{2} L_g |\Delta|^2 \right) \|v\|_{V_0} \|u_1 - u_2\|_{V_0}. \quad (11) $$
Proof. Let us denote $x^{(1)}(t) = \ell(t) + u_1(t)$ and $x^{(2)}(t) = \ell(t) + u_2(t)$. It holds

$$b(u_1)(v) - b(u_2)(v) = \int_0^1 \dot{v}(t)^T [g(x^{(1)}(t)) - g(x^{(2)}(t))] \Delta t \quad (12) \\ + \frac{1}{2} \int_0^1 \Delta^T \left( [\nabla_x g(x^{(1)}(t)) - \nabla_x g(x^{(2)}(t))] \cdot v(t) \right) \Delta t.$$

We first consider the first term:

$$I_1 = \int_0^1 \sum_{i,j=1}^d \dot{v}_i(t) [g_{ij}(x^{(1)}(t)) - g_{ij}(x^{(2)}(t))] \Delta_j dt.$$

Since $g \in W^{1,\infty}(\mathcal{B}_M, S^{d \times d})$, the following Lipschitz estimate holds for almost all $(x, y) \in (\mathcal{B}_M)^2$ and for all $1 \le i, j \le d$:

$$|g_{ij}(x) - g_{ij}(y)| \le \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} |x - y|.$$

Therefore,

$$|I_1| \le \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \int_0^1 |u_1(t) - u_2(t)| \sum_{i=1}^d |\dot{v}_i(t)| \sum_{j=1}^d |\Delta_j| dt,$$

and using the Cauchy–Schwarz inequality, first on vectors and then on functions, we obtain

$$|I_1| \le d \|\nabla g\|_{L^\infty(\mathcal{B}_M)} |\Delta| \|u_1 - u_2\|_{V_0} \|v\|_{V_0}.$$

Concerning the second term in Eq. (12), we denote

$$I_2 = \frac{1}{2} \int_0^1 \sum_{i,j=1}^d \sum_{l=1}^d \left( \partial_{x_l} g_{ij}(x^{(1)}(t)) - \partial_{x_l} g_{ij}(x^{(2)}(t)) \right) v_l(t) \Delta_i \Delta_j dt.$$

By hypothesis, the gradient of the metric is Lipschitz continuous with constant $L_g$, hence

$$|\partial_{x_l} g_{ij}(x^{(1)}(t)) - \partial_{x_l} g_{ij}(x^{(2)}(t))| \le L_g |u_1(t) - u_2(t)|.$$

By virtue of the Cauchy–Schwarz inequality on vectors, we get

$$|I_2| \le \frac{d^{\frac{3}{2}}}{2} L_g |\Delta|^2 \int_0^1 |u_1(t) - u_2(t)| |v(t)| dt.$$

Therefore,

$$|I_2| \le \frac{d^{\frac{3}{2}}}{2} L_g |\Delta|^2 \|u_1 - u_2\|_{V_0} \|v\|_{V_0},$$

which concludes the proof. $\square$
**Lemma 5.** Let $g \in C^{1,1}(\mathcal{B}_M, \mathcal{S}^{d \times d})$, namely the gradient of the metric is uniformly Lipschitz continuous on $\mathcal{B}_M$, with constant $L_g > 0$. Let $u_1, u_2 \in V_0$ be curves with maximal finite speed less or equal, we denote $\sup_{\tau \in [0,1]} |\Delta + \dot{u}_{1,2}(\tau)| = c$. Then, for all $(w, v) \in V_0 \times V_0$:

$$|(A(u_1) - A(u_2))(w, v)| \le \left( \left(d + \frac{d^3}{2}\right) \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} + \frac{d^3}{2} L_g (|\Delta| + c) \right) \|w\|_{V_0} \|v\|_{V_0} \|u_2 - u_1\|_{V_0}.$$

*Proof.* Let us denote $x_1(t) = \ell(t) + u_1(t)$ and $x_2(t) = \ell(t) + u_2(t)$, it holds

$$(A(u_1) - A(u_2))(w, v) = \int_0^1 \dot{v}(t)^T [g(x_1(t)) - g(x_2(t))] \dot{w}(t) dt + \int_0^1 \left( \left[ \left( \Delta + \frac{1}{2}\dot{u}_1(t) \right)^T \nabla_x g(x_1(t)) - \left( \Delta + \frac{1}{2}\dot{u}_2(t) \right)^T \nabla_x g(x_2(t)) \right] \cdot v(t) \right) \dot{w}(t) dt.$$

For the first term we proceed in the same way as in the proof of Lemma 4, it holds

$$|I_1| = \int_0^1 |\dot{v}(t)^T [g(x_1(t)) - g(x_2(t))] \dot{w}(t)| dt \le d \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \int_0^1 |u_1(t) - u_2(t)| |\dot{v}(t)| |\dot{w}(t)| dt.$$

We can bound $|u_1(t) - u_2(t)|$ with the $L^\infty$-norm, it holds

$$|I_1| \le d \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \|u_1 - u_2\|_{L^\infty([0,1])^d} \int_0^1 |\dot{v}(t)| |\dot{w}(t)| dt.$$

By the one-dimensional Sobolev embedding, $\|u_1 - u_2\|_{L^\infty([0,1])^d} \le \|u_1 - u_2\|_{V_0}$, hence

$$|I_1| \le d \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \|u_1 - u_2\|_{V_0} \|v\|_{V_0} \|w\|_{V_0}.$$

For the second term we denote

$$I_2 = \int_0^1 \left( \left[ \left( \Delta + \frac{1}{2} \dot{u}_1(t) \right)^T \nabla_x g(x_1(t)) - \left( \Delta + \frac{1}{2} \dot{u}_2(t) \right)^T \nabla_x g(x_2(t)) \right] \cdot v(t) \right) \dot{w}(t) dt,$$

and we compute

$$I_2 = \int_0^1 \left( \left[ \left( \Delta + \frac{1}{2} \dot{u}_1(t) \right)^T (\nabla_x g(x_1(t)) - \nabla_x g(x_2(t))) \right] \cdot v(t) \right) \dot{w}(t) dt \quad (13) \\ + \int_0^1 \left( \left[ \frac{1}{2} (\dot{u}_1(t) - \dot{u}_2(t))^T \nabla_x g(x_2(t)) \right] \cdot v(t) \right) \dot{w}(t) dt.$$

For the first term of (13) we denote

$$I_2^{(1)} = \int_0^1 \sum_{i,j=1}^d \sum_{l=1}^d (\partial_{x_l} g_{ij}(x_1(t)) - \partial_{x_l} g_{ij}(x_2(t))) \left( \Delta_i + \frac{1}{2}[\dot{u}_1]_i(t) \right) v_l(t) \dot{w}_j(t) dt.$$
Using the Lipschitz continuity of the gradient of the metric we get

$$|I_2^{(1)}| \le L_g \int_0^1 |u_1(t) - u_2(t)| \sum_{i=1}^d \left| \Delta_i + \frac{1}{2}[\dot{u}_1]_i(t) \right| \sum_{j=1}^d |\dot{w}_j(t)| \sum_{l=1}^d |v_l(t)| dt.$$

Let us remind that

$$\sum_{i=1}^{d} \left| \Delta_i + \frac{1}{2}[\dot{u}_1]_i(t) \right| \le \frac{1}{2} \sum_{i=1}^{d} \left| \Delta_i + (\Delta_i + [\dot{u}_1]_i(t)) \right| \le \frac{d^{\frac{1}{2}}}{2} (|\Delta| + c),$$

and we also bound $|u_1(t) - u_2(t)|$ with the $L^\infty$-norm, we get

$$|I_2^{(1)}| \le \frac{d^{\frac{3}{2}}}{2} L_g (|\Delta| + c) \|u_1 - u_2\|_{L^\infty([0,1])^d} \|v\|_{V_0} \|w\|_{V_0} \le \frac{d^{\frac{3}{2}}}{2} L_g (|\Delta| + c) \|u_1 - u_2\|_{V_0} \|v\|_{V_0} \|w\|_{V_0}$$

For the second term of (13) we denote

$$I_2^{(2)} = \frac{1}{2} \int_0^1 \sum_{i,j=1}^d \sum_{l=1}^d \partial_{x_l} g_{ij}(x_2(t)) ([\dot{u}_1]_i(t) - [\dot{u}_2]_i(t)) v_l(t) \dot{w}_j(t) dt.$$

In a similar way, it holds

$$|I_2^{(2)}| \le \frac{d^{\frac{3}{2}}}{2} \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \|u_1 - u_2\|_{V_0} \|v\|_{V_0} \|w\|_{V_0},$$

and this concludes the proof.

<figure><img src="image_5.png" /></figure>

We have now all the elements to prove the conditions under which the solution operator $T$ is a contraction.

**Proposition 2.** Let $g \in C^{1,1}(\mathcal{B}_M, S^{d \times d})$, namely the gradient of the metric is uniformly Lipschitz continuous on $\mathcal{B}_M$, with constant $L_g > 0$. Let $c^* > 0$, and consider the subset of $V_0$ defined as

$$\mathcal{V} = \left\{ u \in V_0 \mid \sup_{t \in [0,1]} |\Delta + \dot{u}(t)| \le c^* \right\}.$$

We also consider

$$M^* = \sup_{u \in \mathcal{V}} \frac{C_b(u)}{c_A(u)},$$

where $C_b(u)$ is the continuity constant of $b(u)$ defined by (8), and $c_A(u)$ the coercivity constant of $A(u)$ defined by (9). Let $\lambda_d^*$ denote the minimal eigenvalue of $g$ over $\mathcal{B}_M$. Assume that

$$\frac{\lambda_d^*}{(1 + C_P^2)} \ge \left( d|\Delta| + \left(d + \frac{d^{\frac{3}{2}}}{2}\right) M^* + \frac{d^{\frac{3}{2}}}{4} (|\Delta| + c^*) \right) \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \\ + \frac{d^{\frac{3}{2}}}{2} L_g (|\Delta|^2 + (|\Delta| + c^*) M^*).$$

Then the fixed-point operator is a contraction.
Proof. Observe that the coercivity constant in Lemma 3 admits the lower bound

$$c_A(u) \geq c_A^* = \frac{\lambda_d^*}{(1 + C_P^2)} - \frac{1}{4} (|\Delta| + c^*) d^3 \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} \quad (14)$$

Moreover, from Lemma 4 we have the Lipschitz continuity of $b$ on $\mathcal{V}$, namely

$$|b(u_1)(v) - b(u_2)(v)| \leq L_b \|v\|_{V_0} \|u_1 - u_2\|_{V_0},$$

where

$$L_b = d \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} |\Delta| + \frac{d^3}{2} L_g |\Delta|^2. \quad (15)$$

Finally, from Lemma 5 we have the Lipschitz continuity of $A$ on $\mathcal{V}$, in the sense that

$$|(A(u_1) - A(u_2))(w, v)| \leq L_A \|w\|_{V_0} \|v\|_{V_0} \|u_1 - u_2\|_{V_0},$$

where

$$L_A = \left(d + \frac{d^3}{2}\right) \|\nabla_x g\|_{L^\infty(\mathcal{B}_M)} + \frac{d^3}{2} L_g (|\Delta| + c^*). \quad (16)$$

Using these estimates in (10) with $v = w_1 - w_2$, we get

$$c_A^* \|w_1 - w_2\|_{V_0}^2 \leq L_b \|u_1 - u_2\|_{V_0} \|w_1 - w_2\|_{V_0} + L_A \|u_1 - u_2\|_{V_0} \|w_2\|_{V_0} \|w_1 - w_2\|_{V_0}.$$

Cancel $\|w_1 - w_2\|_{V_0}$ (the inequality is trivial if it vanishes) to obtain

$$\|w_1 - w_2\|_{V_0} \leq \frac{1}{c_A^*} (L_b + L_A \|w_2\|_{V_0}) \|u_1 - u_2\|_{V_0}.$$

By definition,

$$A(u_2)(w_2, w_2) = b(u_2)(w_2).$$

From the coercivity of $A$ and the continuity of $b(u_2)$ we get

$$c_A(u_2) \|w_2\|_{V_0}^2 \leq C_b(u_2) \|w_2\|_{V_0},$$

and therefore $\|w_2\|_{V_0} \leq M^*$. Finally,

$$\|w_1 - w_2\|_{V_0} \leq \frac{1}{c_A^*} (L_b + L_A M^*) \|u_1 - u_2\|_{V_0}.$$

Thus $T$ is a contraction if and only if $L_b + L_A M^* \leq c_A^*$, which is precisely the explicit hypothesis stated above, and this concludes the proof. □

**Proposition 3.** Under the same assumptions as in Proposition 2, we set $u^{(0)} = 0 \in \mathcal{V}$ and define the sequence by

$$u^{(k+1)} = T(u^{(k)}), \text{ for } k \geq 0.$$

Then $\lim_{k \to \infty} u^{(k)} = u^*$, where $u^*$ is the unique fixed point of $T$ and a stationary point of the geodesic energy. Moreover, the rate of convergence is described by the inequality

$$\|u^{(k)} - u^*\|_{V_0} \leq q^k \|u^{(0)} - u^*\|_{V_0},$$

where $q \in (0, 1)$ is the contraction constant of $T$.
*Proof.* This follows directly from the Banach fixed-point theorem, since Proposition 2 ensures that $T$ is a strict contraction on the complete metric space $(\mathcal{V}, \|\cdot\|_{V_0})$. It remains to show that the iterates $u^{(k)}$ stay within in $\mathcal{V}$.

Let $u^{(k-1)} \in \mathcal{V}$ be given, and consider $u^{(k)} = T(u^{(k-1)})$. From the defining equation satisfied by $u^{(k)}$, we get

$$ \frac{d}{dt} \left( \sum_{j=1}^{d} g_{lj}(\ell(t) + u^{(k-1)}(t)) (\Delta_j + \dot{u}_j^{(k)}(t)) \right) - \frac{1}{2} \sum_{i,j=1}^{d} \partial_{x_i} g_{ij}(\ell(t) + u^{(k-1)}(t)) \Delta_i \Delta_j \\ - \sum_{i,j=1}^{d} \partial_{x_i} g_{ij}(\ell(t) + u^{(k-1)}(t)) (\Delta_i + \frac{1}{2} \dot{u}_i^{(k-1)}(t)) \dot{u}_j^{(k)}(t) = 0. $$

From this we infer that $\ddot{u}^{(k)} \in L^2([0, 1]; \mathcal{B}_M)$, since all coefficients are in $L^2([0, 1]; \mathcal{B}_M)$. By the Sobolev embedding theorem, this in turn yields that $\dot{u}^{(k)}$ is bounded almost everywhere, i.e.,

$$ \|\dot{u}^{(k)}\|_{L^{\infty}([0,1]; \mathcal{B}_M)} < \infty. $$

Hence, there exists a constant $c^* > 0$ such that $\sup_{t \in [0,1]} |\Delta + \dot{u}^{(k)}(t)| \le c^*$ for all $k \ge 0$, and consequently $u^{(k)} \in \mathcal{V}$. It concludes the proof. $\square$

**Remark:** From the analysis of the proposed fixed-point iteration, we can highlight a few crucial aspects. First, let us note that the hypotheses required for the iteration to be well-posed and convergent are quite mild: we assume that the metric tensor belongs to $C^{1,1}(\mathcal{B}_X, S^{d \times d})$, i.e., we require the Lipschitz continuity of the metric gradient in the ball surrounding the endpoints. Furthermore, we only need to assume that the velocity of the curve obtained at the previous iteration is finite. This assumption is reasonable and can be verified numerically. Moreover, as the iteration converges, the velocity is expected to progressively decrease until it matches the logarithmic map $\log_{x_0} x_1$.

Concerning the contractivity of the fixed point, the closer the endpoints are (that is, for small $|\Delta|$), the more contractive the fixed-point operator becomes, as expected. Moreover, in this regime, we can also expect that the constant $c^*$ can be chosen small. However, contractivity can be compromised in situations where the dimensionality of the manifold is large (due to the $d^{3/2}$ dependence), when the gradient of the metric is large, or when the metric becomes nearly degenerate (that is, when the smallest singular value of the metric tensor tends to zero). We therefore expect the proposed iteration to perform effectively in cases where the metric does not degenerate and its gradient remains sufficiently small, in the sense specified by the conditions of the above proposition, within the ball surrounding the endpoints.

## 2.4 Finite element discretization

Let $I_h$ be a partition of $[0, 1]$ with maximum step size $h$. Let $V_h \subset V_0$ be a conforming finite element space. For instance, let $(\phi_n)_{1 \le n \le N_h}$ be a finite element scalar basis of $H_0^1([0, 1])$.
Then, we can construct the corresponding vector basis as

$$ \Phi_{j,n} = \phi_n e_j, \quad \text{for } 1 \le j \le d \text{ and } 1 \le n \le N_h, $$

where $(e_j)_{1 \le j \le d}$ denotes the canonical basis of $\mathbb{R}^d$. Therefore, any function $u_h^{(k)} \in V_h$ can be written as

$$ u_h^{(k)}(t) = \sum_{j=1}^{d} \sum_{n=1}^{N_h} U_{j,n}^{(k)} \Phi_{j,n}(t), \quad \text{for } t \in [0, 1]. $$

In this discrete framework, the fixed-point iteration takes the form

$$ A(u_h^{(k-1)})(u_h^{(k)}, v_h) = b(u_h^{(k-1)})(v_h), \quad \text{for all } v_h \in V_h. \qquad (17) $$

We now express this linear problem in matrix form. It holds that

$$ A(u_h^{(k-1)})(u_h^{(k)}, v) = \sum_{\alpha=1}^{d} \sum_{\beta=1}^{d} \sum_{n=1}^{N_h} U_{\alpha,n}^{(k)} \int_{0}^{1} g_{\alpha\beta}(\ell(t) + u_h^{(k-1)}(t)) \dot{\phi}_n(t) \dot{v}_\beta(t) dt \\ + \sum_{\alpha=1}^{d} \sum_{\beta=1}^{d} \sum_{n=1}^{N_h} U_{\alpha,n}^{(k)} \int_{0}^{1} \left( \sum_{l=1}^{d} \partial_{x_l} g_{\alpha\beta}(\ell(t) + u_h^{(k-1)}(t)) v_l(t) \right) \dot{\phi}_n(t) \left( \Delta_\beta + \frac{1}{2} [\dot{u}_h^{(k-1)}(t)]_\beta \right) dt. $$

Taking $v = \Phi_{i,m}$ and renaming $\alpha$ as $j$, we obtain

$$ A(u_h^{(k-1)})(u_h^{(k)}, \Phi_{i,m}) = \sum_{j=1}^{d} \sum_{n=1}^{N_h} U_{j,n}^{(k)} \int_{0}^{1} g_{ij}(\ell(t) + u_h^{(k-1)}(t)) \dot{\phi}_n(t) \dot{\phi}_m(t) dt \\ + \sum_{j=1}^{d} \sum_{n=1}^{N_h} U_{j,n}^{(k)} dt \int_{0}^{1} \sum_{\beta=1}^{d} \left( \partial_{x_i} g_{j\beta}(\ell(t) + u_h^{(k-1)}(t)) \phi_m(t) \right) \dot{\phi}_n(t) \left( \Delta_\beta + \frac{1}{2} [\dot{u}_h^{(k-1)}(t)]_\beta \right) dt. $$

Similarly, we have

$$ b(u_h^{(k-1)})(\Phi_{i,m}) = \sum_{\alpha=1}^{d} \int_{0}^{1} g_{i\alpha}(\ell(t) + u_h^{(k-1)}(t)) \Delta_{\alpha} \dot{\phi}_m(t) dt \\ + \frac{1}{2} \sum_{\alpha=1}^{d} \int_{0}^{1} \sum_{\beta=1}^{d} \left( \partial_{x_i} g_{\alpha\beta}(\ell(t) + u_h^{(k-1)}(t)) \phi_m(t) \right) \Delta_{\alpha} \Delta_{\beta} dt. $$

Let us collect the unknown coefficients into a vector

$$ \mathbf{U}^{(k)} = (U_{1,1}^{(k)}, \dots, U_{1,N_h}^{(k)}, U_{2,1}^{(k)}, \dots, U_{d,N_h}^{(k)}) \in \mathbb{R}^{dN_h}. $$

Since the above relations hold for all $\Phi_{i,n}$, we obtain the linear system

$$ \mathbf{K}(\mathbf{U}^{(k-1)})\mathbf{U}^{(k)} = \mathbf{B}(\mathbf{U}^{(k-1)}), \qquad (18) $$
where $\mathbf{K}(\mathbf{U}^{(k-1)})$ is a $dN_h \times dN_h$ matrix with block entries $\mathbf{K}^{(n,m)}(\mathbf{U}^{(k-1)}) \in \mathbb{R}^{d \times d}$, such that for all $1 \le i, j \le d$ and $1 \le m, n \le N_h$:

$$ \begin{aligned} \left[ \mathbf{K}^{(n,m)}(\mathbf{U}^{(k-1)}) \right]_{ij} = & \int_0^1 g_{ij}(\ell(t) + u_h^{(k-1)}(t)) \dot{\phi}_n(t) \dot{\phi}_m(t) \, dt \\ & + \int_0^1 \sum_{\beta=1}^d \left( \partial_{x_i} g_{j\beta}(\ell(t) + u_h^{(k-1)}(t)) \phi_m(t) \right) \dot{\phi}_n(t) \left( \Delta_{\beta} + \frac{1}{2} [\dot{u}_h^{(k+1)}(t)]_{\beta} \right) \, dt, \end{aligned} $$

and $\mathbf{B}(\mathbf{U}^{(k-1)})$ is a vector in $\mathbb{R}^{dN_h}$ with components $\mathbf{B}^{(m)}(\mathbf{U}^{(k-1)}) \in \mathbb{R}^d$, such that for $1 \le i \le d$ and $1 \le m \le N_h$:

$$ \begin{aligned} \left[ \mathbf{B}^{(m)}(\mathbf{U}^{(k-1)}) \right]_i = & \sum_{\alpha=1}^{d} \int_0^1 g_{i\alpha}(\ell(t) + u_h^{(k-1)}(t)) \Delta_{\alpha} \dot{\phi}_m(t) \, dt \\ & + \frac{1}{2} \sum_{\alpha=1}^{d} \int_0^1 \sum_{\beta=1}^{d} \left( \partial_{x_i} g_{\alpha\beta}(\ell(t) + u_h^{(k-1)}(t)) \phi_m(t) \right) \Delta_{\alpha} \Delta_{\beta} \, dt. \end{aligned} $$

The well-posedness and convergence of the discretized fixed-point iteration follow directly from the corresponding continuous results.

**Proposition 4 (Convergence of the discrete fixed-point iteration).** Let $u_h^{(0)} = 0 \in \mathcal{V} \cap V_h$ and define the sequence

$$ u_h^{(k+1)} = T_h(u_h^{(k)}), \quad k \ge 0, $$

where $T_h$ is the discrete solution operator associated with the fixed-point iteration (17). Under the same assumptions as in Proposition 3, the discrete fixed-point iteration is well-posed and convergent. Denote by $u_h^* = \lim_{k \to \infty} u_h^{(k)}$ the unique fixed point of $T_h$. Moreover, for all $k \ge 0$ and all $v_h \in V_h$, it holds that

$$ \|u_h^{(k)} - u^*\|_{V_0} \le q^k \|u_h^*\|_{V_0} + C^{st} \|u^* - v_h\|_{V_0}, $$

where $q \in (0,1)$ is the contraction constant of $T$ stated in Proposition 2, and $C^{st} > 0$ is a constant independent of $h$. In particular,

$$ \|u_h^{(k)} - u^*\|_{V_0} \le q^k \|u_h^*\|_{V_0} + C^{st} \inf_{v_h \in V_h} \|u^* - v_h\|_{V_0}. $$

*Proof.* Well-posedness and convergence follow from the Banach fixed-point theorem together with Proposition 2, applied on the finite-dimensional space $V_h$. The discrete iteration satisfies

$$ \|u_h^{(k)} - u_h^*\|_{V_0} \le q^k \|u_h^* - u_h^{(0)}\|_{V_0} = q^k \|u_h^*\|_{V_0}. $$

By the triangle inequality,

$$ \|u_h^{(k)} - u^*\|_{V_0} \le \|u_h^{(k)} - u_h^*\|_{V_0} + \|u_h^* - u^*\|_{V_0}. \quad (19) $$
To estimate the second term, introduce $\tilde{u}_h \in V_h$ defined by

$$A(u^*)(\tilde{u}_h, v_h) = b(u^*)(v_h), \quad \forall v_h \in V_h.$$

Then,

$$\|u_h^* - u^*\|_{V_0} \le \|\tilde{u}_h - u^*\|_{V_0} + \|u_h^* - \tilde{u}_h\|_{V_0}. \quad (20)$$

The first term is bounded by Céa's lemma:

$$\|\tilde{u}_h - u^*\|_{V_0} \le C_{C\text{éa}} \inf_{v_h \in V_h} \|u^* - v_h\|_{V_0},$$

where $C_{C\text{éa}} > 0$. From the definitions of $u_h^*$ and $\tilde{u}_h$, we infer that

$$A(u_h^*)(u_h^* - \tilde{u}_h, v_h) = b(u_h^*)(v_h) - A(u_h^*)(\tilde{u}_h, v_h) + (A(u^*)(\tilde{u}_h, v_h) - b(u^*)(v_h)), \quad \forall v_h \in V_h.$$

Setting $e_h = u_h^* - \tilde{u}_h$ and $v_h = e_h$ yields

$$A(u_h^*)(e_h, e_h) = (b(u_h^*) - b(u^*)) (e_h) + (A(u_h^*) - A(u^*)) (\tilde{u}_h, e_h). \quad (21)$$

We now reintroduce the notations used in the proof of Proposition 2. Namely, $c_A^*$ denotes the uniform coercivity constant of A defined in (14), $L_A$ the Lipschitz continuity constant of A defined in (16), and $L_b$ the Lipschitz continuity constant of b defined in (15). We also recall the notation

$$M^* = \sup_{u \in V} \frac{C_b(u)}{c_A(u)},$$

where $C_b(u)$ is the continuity constant of $b(u)$ defined by (8), and $c_A(u)$ the coercivity constant of $A(u)$ defined by (9). From this definition, we infer the bound

$$\|\tilde{u}_h\|_{V_0} \le M^*.$$

From (21), it holds that

$$c_A^* \|e_h\|_{V_0}^2 \le (L_b + L_A M^*) \|u_h^* - u^*\|_{V_0} \|e_h\|_{V_0}.$$

Dividing by $\|e_h\|_{V_0}$, as the inequality is trivial if it is zero, gives

$$\|e_h\|_{V_0} \le \frac{L_b + L_A M^*}{c_A^*} \|u_h^* - u^*\|_{V_0}.$$

Substituting this estimate into (20) yields

$$\left(1 - \frac{L_b + L_A M^*}{c_A^*}\right) \|u_h^* - u^*\|_{V_0} \le C_{C\text{éa}} \inf_{v_h \in V_h} \|u^* - v_h\|_{V_0}.$$

We recognize here the contractivity condition of the operator T, namely

$$q = \frac{L_b + L_A M^*}{c_A^*}.$$

Since $q < 1$, the prefactor on the left-hand side equals $(1 - q) > 0$. Combining this estimate with (19) gives the desired bound and concludes the proof. □
### 3 The Hessian of the squared distance

In this section, we investigate the efficient computation of the Hessian of the squared distance function. Let $x_0, x_1 \in \mathcal{B}_M$ be such that there exists a unique geodesic $x \in V$ with fixed endpoints $x(0) = x_0$ and $x(1) = x_1$. Define

$$ \rho_{x_1}(x_0) = \frac{1}{2} \mathrm{dist}(x_0, x_1)^2. $$

It is well known that the Riemannian gradient of $\rho_{x_1}$ satisfies

$$ \mathrm{grad} \rho_{x_1}(x_0) = - \log_{x_0}(x_1) = - \dot{x}(0). $$

To compute the Hessian, we consider a variation $x_{s,j}$ defined as the geodesic between $x_0 + se_j$ and $x_1$, where $(e_j)_{1 \le j \le d}$ denotes the canonical basis of $\mathbb{R}^d$. For any $1 \le i, j \le d$, we denote

$$ z_{i;j} = \left[ \frac{d}{ds} x_{s,j} \Big|_{s=0} \right]_i, \quad \dot{z}_{i;j} = \left[ \frac{d}{dt} \frac{d}{ds} x_{s,j} \Big|_{s=0} \right]_i = \left[ \frac{d}{ds} \dot{x}_{s,j} \Big|_{s=0} \right]_i. $$

Under the Levi-Civita connection — the unique torsion-free connection compatible with the metric $g$ — the Hessian of $\rho_{x_1}$ at $x_0$ is given in coordinates by

$$ [\mathrm{Hess} \, \rho_{x_1}(x_0)]_{ij} = - \left( \dot{z}_{i;j}(0) + \sum_{k=1}^{d} \Gamma_{ij}^{k} \dot{x}_{k}(0) \right), \qquad (22) $$

where $\Gamma_{ij}^k$ are the Christoffel symbols of the second kind on $M$. This construction, often called the *geometric Hessian* [1, 2], represents the covariant derivative of the gradient and hence defines an endomorphism of the tangent space. It differs from the *covariant Hessian* [13, 14], denoted $\nabla^2 \rho_{x_1}(x_0)$, which is the second-order covariant derivative and thus a bilinear form on the tangent space. Both notions are equivalent under the *musical isomorphism*.

#### 3.1 Sensitivity analysis of the geodesic equation

From the previous section, the geodesic $x$ satisfies, for all $v \in V$,

$$ \sum_{i,j=1}^{d} \int_{0}^{1} g_{ij}(x(t)) \dot{x}_{i}(t) \dot{v}_{j}(t) + \frac{1}{2} \sum_{l=1}^{d} \partial_{x_{l}} g_{ij}(x(t)) v_{l}(t) \dot{x}_{i}(t) \dot{x}_{j}(t) dt = 0, \qquad (23) $$

Differentiating (23) with respect to the perturbation parameter yields the sensitivity equations satisfied by $z_{i;k}$:

$$ \begin{aligned} & \sum_{i,j=1}^{d} \int_{0}^{1} g_{ij}(x(t)) \dot{z}_{i;k}(t) \dot{v}_{j}(t) + \sum_{l=1}^{d} \partial_{x_{l}} g_{ij}(x(t)) z_{l;k}(t) \dot{x}_{i}(t) \dot{v}_{j}(t) dt \\ & + \frac{1}{2} \int_{0}^{1} \sum_{l=1}^{d} \partial_{x_{l}} g_{ij}(x(t)) v_{l}(t) (\dot{z}_{i;k}(t) \dot{x}_{j}(t) + \dot{x}_{i}(t) \dot{z}_{j;k}(t)) dt \qquad (24) \\ & + \frac{1}{2} \int_{0}^{1} \sum_{l,m=1}^{d} \partial_{x_{m}} \partial_{x_{l}} g_{ij}(x(t)) z_{m;k}(t) v_{l}(t) \dot{x}_{i}(t) \dot{x}_{j}(t) dt = 0, \end{aligned} $$
subject to the boundary conditions

$$z_{i;k}(0) = \delta_{ik}, \quad z_{i;k}(1) = 0.$$

The problem (24) is linear with respect to $z_{i;k} = (z_{1;k}, z_{2;k}, \dots, z_{d;k})$, and its structure closely resembles that of the problem (7). If the metric tensor possesses sufficient regularity, for instance

$$g \in W^{2,\infty}(\mathcal{B}_M, \mathcal{S}^{d \times d}),$$

then we can show that the bilinear form associated with (24) is continuous and coercive on $V$ under suitable conditions, by following the same steps as in the proof of Proposition 1. Hence, by the Lax-Milgram theorem, the problem admits a unique solution.

## 3.2 Finite elements discretization

We now consider the finite element framework introduced in section 2.4. We seek an approximate solution of the form

$$z_{i;k}(t) = (1-t) \delta_{ik} + \sum_{n=1}^{N_h} Z_{(i,n);k} \phi_n(t), \quad t \in [0,1],$$

where $\{\phi_n\}_{n=1}^{N_h}$ denotes the chosen scalar finite element basis on the unit interval. The unknown coefficients are collected into the vector

$$\mathbf{Z}_{;k} = (Z_{(1,1);k}, \dots, Z_{(1,N_h);k}, Z_{(2,1);k}, \dots, Z_{(d,N_h);k}) \in \mathbb{R}^{dN_h}.$$

The sensitivity equations (24) are linear with respect to $\mathbf{Z}_{;k}$; hence, the coefficients $\mathbf{Z}_{;k}$ can be recovered by solving the linear system

$$\mathbf{A} \mathbf{Z}_{;k} = \mathbf{F}_{;k}.$$

The matrix $\mathbf{A} \in \mathbb{R}^{(dN_h) \times (dN_h)}$ has block entries $\mathbf{A}^{(n,m)} \in \mathbb{R}^{d \times d}$, such that for all $1 \le i, j \le d$ and $1 \le m, n \le N_h$:

$$\begin{aligned} [\mathbf{A}^{(n,m)}]_{ij} = & \int_0^1 g_{ij}(x(t)) \dot{\phi}_n(t) \dot{\phi}_m(t) dt \\ & + \sum_{\alpha=1}^d \int_0^1 \partial_{x_i} g_{\alpha j}(x(t)) \phi_n(t) \dot{x}_\alpha(t) \dot{\phi}_m(t) dt \\ & + \frac{1}{2} \int_0^1 \partial_{x_j} g_{ij}(x(t)) \dot{\phi}_n(t) (\dot{x}_j(t) + \dot{x}_i(t) \delta_{ij}) \phi_m(t) dt \\ & + \frac{1}{2} \sum_{\alpha=1}^d \int_0^1 \partial_{x_i} \partial_{x_j} g_{\alpha j}(x(t)) \phi_n(t) \dot{x}_\alpha(t) \dot{x}_j(t) \phi_m(t) dt. \end{aligned}$$
The corresponding right-hand side $\mathbf{F}_{;k} \in \mathbb{R}^{dN_h}$ collects the contribution of the inhomogeneous boundary terms, i.e., the parts proportional to $(1-t)\delta_{ik}$ and its derivative $-\delta_{ik}$. It has block components $\mathbf{F}_{;k}^{(m)} \in \mathbb{R}^d$, such that for all $1 \le j \le d$ and $1 \le m \le N_h$:

$$[\mathbf{F}_{;k}^{(m)}]_i = -\sum_{\alpha=1}^{d} \int_{0}^{1} g_{\alpha i}(x(t)) (-\delta_{\alpha k}) \dot{\phi}_m(t) dt \\ -\sum_{\alpha=1}^{d} \sum_{l=1}^{d} \int_{0}^{1} \partial_{x_l} g_{\alpha i}(x(t)) ((1-t)\delta_{lk}) \dot{x}_{\alpha}(t) \dot{\phi}_m(t) dt \\ -\frac{1}{2} \sum_{\alpha=1}^{d} \int_{0}^{1} \partial_{x_i} g_{\alpha i}(x(t)) [(-\delta_{\alpha k}) \dot{x}_i(t) + \dot{x}_{\alpha}(t) (-\delta_{ik})] \phi_m(t) dt \\ -\frac{1}{2} \sum_{\alpha=1}^{d} \sum_{l=1}^{d} \int_{0}^{1} \partial_{x_l} \partial_{x_i} g_{\alpha i}(x(t)) ((1-t)\delta_{lk}) \dot{x}_{\alpha}(t) \dot{x}_i(t) \phi_m(t) dt.$$

Thus, assembling $\mathbf{A}$ and $\mathbf{F}_{;k}$ yields a linear system of dimension $dN_h$, whose solution provides the coefficients $\mathbf{Z}_{;k}$. Solving this system for each $k = 1, \dots, d$ gives access to the quantities $\dot{z}_{i;k}(0)$, and hence to the Hessian entries $[\text{Hess } \rho_{x_1}(x_0)]_{ik}$ via (22). Moreover, the matrix $\mathbf{A}$ is independent of $k$, and therefore needs to be assembled only once. For comparison, a finite-difference approach would necessitate solving $\frac{d(d+1)}{2}$ nonlinear systems of the type introduced in section 2.

# 4 Numerical experiments

In this section, we present numerical experiments on the unit sphere embedded in $\mathbb{R}^3$,

$$\mathcal{M} = \{ \mathbf{x} \in \mathbb{R}^3 \mid |\mathbf{x}| = 1 \}.$$

We consider the geodesic $\mathbf{x}(t)$ between the points

$$\mathbf{x}(0) = (0.6573, 0.0268, 0.7532) \quad \text{and} \quad \mathbf{x}(1) = (0.6453, -0.7138, -0.2722).$$

Interpreting $\mathcal{M}$ as an idealized Earth, with the North pole at $N = (0, 0, 1)$ and the South pole at $S = (0, 0, -1)$, these two points correspond approximately to the locations of Paris and Brasilia, respectively. We work in a chart containing both points, given by stereographic projection from the North pole. Namely, we define

$$\mathcal{U} = \mathbb{S}^2 \setminus \{N\}, \\ \varphi(\mathbf{x}) = \frac{(\mathbf{x}_1, \mathbf{x}_2)}{1 - \mathbf{x}_3} \in \mathbb{R}^2.$$

In these coordinates, the metric inherited from $\mathbb{R}^3$ is

$$g(x) = \frac{4}{(1 + (x_1)^2 + (x_2)^2)^2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}.$$
We discretize the problem using a $\mathbb{P}_1$ finite element basis on $[0, 1]$ and set $N_h = 1000$ degrees of freedom. The integrals involved are computed using a standard three-point Gauss-Legendre quadrature formula. The fixed-point iteration used to solve the discrete boundary value problem terminates once the residual of the numerical scheme (18) falls below the prescribed tolerance $\epsilon = 10^{-10}$ in the Euclidean 2-norm. To better observe the convergence towards the geodesic curve, we deliberately slowed down the iteration by introducing a relaxation

$$ \mathbf{U}^{(k)} = \omega \tilde{\mathbf{U}}^{(k)} + (1 - \omega) \mathbf{U}^{(k-1)}, \quad \text{where} \quad \mathbf{K}(\mathbf{U}^{(k-1)}) \tilde{\mathbf{U}}^{(k)} = \mathbf{B}(\mathbf{U}^{(k-1)}). $$

In this experiment, we set $\omega = 1/2$ and the relaxed fixed-point iteration converges in 48 iterations (see Figure 1).

<figure><img src="image_6.png" /><figcaption>Figure 1: $\mathbb{P}_1$ Finite element approximation of the boundary value problem, with $N_h = 1000$ degrees of freedom and tolerance $\epsilon = 10^{-10}$. Convergence of the relaxed fixed-point iteration is achieved in 48 iterations.</figcaption></figure>

We next investigate the convergence of the method with respect to the mesh size $h$. Since geodesics on $\mathcal{M}$ and the Hessian of the squared Riemannian distance are available in closed form [17, 18], we can directly compare our numerical approximations with their analytical counterparts. The approximation error for the boundary value problem is measured in the norm induced by the metric $g$, i.e., the Riemannian norm. Because the error from the fixed-point iteration is driven to a negligible residual, the dominant source of error arises from the finite element discretization. As expected for a $\mathbb{P}_1$ discretization, we observe first-order convergence, $O(h)$, in the Riemannian norm (Figure 2a).

We further assess the approximation of the Hessian of the squared distance, using the Frobenius norm to quantify the error. Here again, the convergence rate is of order $O(h)$ (Figure 2b).

# 5 Conclusions and perspectives

In the present work, we have proposed a continuous fixed-point iteration for approximating the boundary value problem of geodesics in Riemannian manifolds. We have proved that,
<figure><img src="image_7.png" /><figcaption>Figure 2: Convergence of the $P1$ finite element approximation for (a) the boundary value problem and (b) the Hessian of the squared distance, with tolerance $\varepsilon = 10^{-10}$. The reference slope 1 is shown for comparison.</figcaption></figure>

under suitable conditions, the iteration is well-posed and that the fixed point converges to a stationary point of the geodesic energy. The resulting iteration was discretized by means of finite elements, and the corresponding discrete scheme was shown to be convergent under the same assumptions. Moreover, we have introduced a method based on the sensitivity equations associated with the Euler-Lagrange equations of the geodesic energy to estimate the Hessian of the squared distance. These methods were assessed on the unit sphere embedded in $\mathbb{R}^3$.

Several perspectives arise from this work, in particular regarding alternative designs of the fixed-point iteration (4) aimed at improving the convergence speed. In addition, the application of these methods to several applications will be investigated in future work.

## References

[1] P.-A. Absil, R. Mahony, and R. Sepulchre. *Optimization algorithms on matrix manifolds*. Princeton University Press, 2008.

[2] N. Boumal. *An Introduction to Optimization on Smooth Manifolds*. Cambridge University Press, 2023.

[3] Alice Le Brigant. Computing distances and geodesics between manifold-valued curves in the srv framework. *arXiv preprint arXiv:1601.02358*, 2016.

[4] Jindong Chen and Yijie Han. Shortest paths on a polyhedron. In *Proceedings of the sixth annual symposium on Computational geometry*, pages 360–369, 1990.

[5] Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. A survey of algorithms for geodesic paths and distances. *arXiv preprint arXiv:2007.10430*, 2020.
[6] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. Geodesics in heat: A new approach to computing distance based on heat flow. *ACM Transactions on Graphics* (TOG), 32(5):1–11, 2013.

[7] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. The heat method for distance computation. *Communications of the ACM*, 60(11):90–99, 2017.

[8] Ricardo Ferreira and Joao Xavier. Hessian of the riemannian squared distance function on connected locally symmetric spaces with applications. In *Controlo 2006, 7th Portuguese conference on automatic control*, volume 2, 2006.

[9] Emin Kasap, Mustafa Yapici, and F Talay Akyildiz. A numerical study for computation of geodesic curves. *Applied Mathematics and Computation*, 171(2):1206–1213, 2005.

[10] Herbert B Keller. *Numerical solution of two point boundary value problems*. SIAM, 1976.

[11] Ron Kimmel. Fast marching methods on triangulated domains. *Proc. Nat. Acad. Sci.*, 95:8431–8435, 1998.

[12] Line Kühnel, Stefan Sommer, and Alexis Arnaudon. Differential geometry and stochastic dynamics with deep learning numerics. *Applied Mathematics and Computation*, 356:411–437, 2019.

[13] John M. Lee. *Introduction to Smooth Manifolds*. Springer New York, 2012.

[14] John M. Lee. *Introduction to Riemannian Manifolds*. Springer New York, 2018.

[15] Facundo Mémoli and Guillermo Sapiro. Fast computation of weighted distance functions and geodesics on implicit hyper-surfaces. *Journal of computational Physics*, 173(2):730–764, 2001.

[16] Joseph SB Mitchell, David M Mount, and Christos H Papadimitriou. The discrete geodesic problem. *SIAM Journal on Computing*, 16(4):647–668, 1987.

[17] Xavier Pennec. Hessian of the riemannian squared distance. 2017.

[18] Xavier Pennec. Barycentric subspace analysis on manifolds. 2018.

[19] Anuj Srivastava and Eric P Klassen. *Functional and shape data analysis*, volume 1. Springer, 2016.

[20] Vitaly Surazhsky, Tatiana Surazhsky, Danil Kirsanov, Steven J Gortler, and Hugues Hoppe. Fast exact and approximate geodesics on meshes. *ACM transactions on graphics (TOG)*, 24(3):553–560, 2005.

[21] Marco Sutti. Shooting methods for computing geodesics on the stiefel manifold. *arXiv preprint arXiv:2309.03585*, 2023.
[22] Sathamangalam R Srinivasa Varadhan. On the behavior of the fundamental solution of the heat equation with variable coefficients. *Communications on Pure and Applied Mathematics*, 20(2):431–455, 1967.

[23] Xiang Ying, Shi-Qing Xin, and Ying He. Parallel chen-han (pch) algorithm for discrete geodesics. *ACM Transactions on Graphics (TOG)*, 33(1):1–11, 2014.