arXiv CSopen access

1-Lipschitz Neural Networks on Hadamard Manifolds

1-Lipschitz Neural Networks on Hadamard Manifolds Davide Murari∗ , Marta Ghirardelli† , Ben Adcock# , Elena Celledoni† , Brynjulf Owren† , Carola-Bibiane Schönlieb∗ Department of Applied Mathematics and Theoretical Physics, University of Cambridge [dm2011, cbs31]@cam.ac.uk †

Department of Mathematical Sciences, Norwegian University of Science and Technology (NTNU), [marta.ghirardelli, elena.celledoni, brynjulf.owren]@ntnu.no

arXiv:2607.19335v1 [math.NA] 21 Jul 2026

Department of Mathematics, Simon Fraser University, [email protected]

ABSTRACT Controlling the Lipschitz constant of a neural network is a standard way to promote robustness and stability. Most existing constraining strategies are designed for Euclidean spaces. In this work, we construct and analyze a class of 1-Lipschitz neural networks on Hadamard manifolds. Our layers are of gradient-descent type, 1-Lipschitz, and quasi-α-firmly nonexpansive. The core building blocks of the proposed architecture are Busemann functions, and we exploit the properties of Busemann gradient flows to design 1-Lipschitz geometry-preserving layers. We provide explicit constructions and examples for hyperbolic manifolds and the manifold of symmetric positive definite (SPD) matrices. We test the proposed architecture in two numerical experiments: robust classification on the Poincaré disk and masked-Wishart covariance reconstruction. On the Poincaré disk, the proposed networks yield robust classifiers under hyperbolic perturbations. On the SPD manifold, we train SPD-valued denoisers and adopt them as a Plug-and-Play prior for a maskedWishart covariance reconstruction problem. We show improved results from the nonexpansive denoiser over static, data-only, and Log-Euclidean denoising baselines, and empirically test its convergence properties. Keywords: 1-Lipschitz neural networks, Hadamard manifolds, Busemann functions, geometric deep learning, hyperbolic neural networks, Plug-and-Play MSCcodes: 68T07, 53C20, 47H09, 47H10

1

Introduction

Controlling the sensitivity of a neural network is important in several settings where stability is a modeling requirement. This is the case in adversarially-robust classification, in Wasserstein-based generative modeling, and in learned iterative methods for inverse problems, where the network often appears inside a fixed-point iteration and must interact well with the surrounding algorithmic structure [1, 25, 26, 48, 59]. In Euclidean spaces, this has led to an extensive literature on architectures with controlled Lipschitz constant. This constraint is usually enforced through spectral normalization, orthogonality constraints, monotone and gradient-type residual blocks [1,25,26,46,48,56,59]. These constructions show that stability can be encoded structurally, rather than checked a posteriori. Many learning problems are not naturally posed in flat spaces. Data may satisfy intrinsic constraints, or admit geometric representations that a curved latent space can capture efficiently. This is one of the main themes of geometric deep learning, where the goal is to design models that respect the geometry of the domain rather than forcing it into Euclidean coordinates [22]. Hyperbolic representation learning is one of the most prominent examples of this paradigm. Negatively curved spaces are well-suited to hierarchical or tree-like data and have therefore led to a rich literature on hyperbolic embeddings and hyperbolic neural networks [32, 53, 54, 58]. More broadly, manifold-valued architectures already appear in intrinsic convolutions

on Riemannian manifolds, networks on the manifold of symmetric positive definite (SPD) matrices, and constructions on noncompact symmetric spaces [39, 45, 52, 61]. Most of this literature, however, focuses on expressive geometry-aware layers rather than a systematic stability theory. This matters in concrete applications. Manifold-valued data arise in directional statistics and phase-valued imaging on circles and spheres [13, 44], in rotation-valued signals and images such as electron backscatter diffraction data on SO(3) [6, 13], and in covariance-based vision where SPD descriptors are used for recognition and classification [37, 62]. A particularly relevant example is diffusion tensor imaging. Here, each voxel is represented by an SPD matrix encoding anisotropic diffusion [36, 63]. These settings call for intrinsic neural layers that respect the geometry of the state space while retaining the stability properties to make them useful in robust learning and inverse problems. Among stability notions, plain 1-Lipschitz continuity is often only the first step. In inverse problems and learned fixed-point methods, the more relevant classes are averaged, and firmly nonexpansive operators, since these are the ones that naturally appear in proximal algorithms, forward–backward splitting, plugand-play methods, and equilibrium formulations [7, 8, 10, 27, 40, 49, 57]. The importance of this viewpoint is that it connects stable layers with convergent learned solvers. The Euclidean architectures most relevant to the present work are those derived from gradient flows. Recent papers have shown that residual blocks of the form x 7→ x − τ ∇V (x) can provide a natural route to nonexpansive (i.e., 1-Lipschitz) and averaged networks [23, 46, 59]. They have also recently been shown to generate universal 1-Lipschitz scalar function approximators [51]. For these layers, the Lipschitz constraint follows from the convexity and L-smoothness of the potential. The same formalism extends beyond Rn . On a Riemannian manifold, one can replace the Euclidean gradient by the Riemannian one, obtaining the map
Tτ (x) = expx −τ grad V (x) . Whenever the manifold geometry supports a workable convexity theory, this gives a principled mechanism for building geometry-preserving stable layers. This is where Hadamard manifolds become the natural setting. These manifolds are complete, simply connected, and nonpositively curved. They therefore combine a broad geometric scope with the global properties needed for analysis. Geodesics are unique, and geodesic convexity admits a global theory closely related to convex analysis in Hilbert spaces [5, 19, 20, 43]. The class includes the standard hyperbolic models used in representation learning, as well as other geometries of independent interest, including the SPD manifolds. Just as importantly, the fixed-point theory survives in this setting. Averagedness must be reformulated through geodesic interpolation, leading to natural extensions of firmly nonexpansive and αfirmly nonexpansive mappings on Hadamard spaces [2,3,11,12]. These results provide the right framework for studying manifold-valued layers. To the best of our knowledge, however, there is no prior work whose main objective is the design and analysis of stable neural networks on manifolds. Existing papers on manifold-valued networks explain how to build intrinsic architectures, while the Euclidean stability literature explains how robustness and convergent dynamics can be obtained from Lipschitz, averaged, or gradient-based constructions. The gap addressed in this paper is precisely between these two threads. We transport Euclidean stability constructions to the Hadamard setting and study whether they remain both analytically useful and practically effective. For concrete layer design, we focus on gradient-type updates generated by geodesically convex potentials. On some important Hadamard manifolds, Busemann functions provide a particularly natural source of such potentials and hence a convenient modeling tool for intrinsic layers. In hyperbolic learning, they already appear, explicitly or implicitly, in prototype methods, horospherical classification, dimensionality reduction, and sliced-Wasserstein constructions [17, 24, 29, 33, 34]. On the manifold of symmetric positive definite matrices, horoballs and horospheres have also been studied from a geometric viewpoint. Explicit formulae for the associated radial fields, that is, the negative gradients of Busemann functions, have recently become available [18, 31, 60]. What is missing in these lines of work is the use of Busemann functions as building

blocks for stable gradient-type layers equipped with nonexpansive or averaged guarantees. This makes them a natural bridge between existing geometric constructions and the framework developed here. 1.1

Main contributions

In this paper, we study intrinsic gradient-descent-type layers on Hadamard manifolds and show how to endow them with stability properties familiar from the Euclidean theory. Our main contributions are as follows. • We establish a bridge between 1-Lipschitz Euclidean gradient layers and intrinsic gradient steps on Hadamard manifolds. We propose Busemann gradient-descent layers and prove their nonexpansiveness under a stepsize restriction. We also show they are quasi α-firmly nonexpansive. • We construct implementable manifold-valued layers from tractable families of potentials, with particular attention to hyperbolic manifolds and the manifold of SPD matrices where Busemann functions yield explicit, computationally convenient formulae. • We validate the proposed layers on two numerical tasks, namely adversarially robust classification in the Poincaré disk and denoising on the manifold of SPD matrices, and compare them with baseline manifold models. The numerical implementation of the networks and experiments can be found at the associated GitHub repository https://github.com/davidemurari/one-lipschitz-hadamard-networks. 1.2

Outline of the paper

The paper is organized as follows. After the notation subsection, Section 2 presents some background material and introduces the property of quasi α-firm nonexpansiveness of certain nonexpansive maps. In Section 3, we introduce Busemann functions and the proposed layers, whose nonexpansiveness condition is derived in Section 4. The implementation of the architecture is discussed in Section 5, then numerical experiments are shown in Section 6, and Section 7 concludes the paper. Deferred proofs are collected in Section A. 1.3

Notation

Throughout the paper, (M, g) denotes a smooth Hadamard manifold. We write Tx M for the tangent space at x ∈ M. We write ⟨·, ·⟩x and ∥ · ∥x for the inner product and norm induced by g on Tx M. When the base point is clear from the context, we write ⟨·, ·⟩ and ∥ · ∥. The length of a differentiable curve γ : [a, b] ⊆ R → M is denoted by Z b ℓ(γ) = ∥γ̇(t)∥dt, a

while d(x, y) is the Riemannian (geodesic) distance between points x and y on M. (M, d) is therefore a metric space. The Levi-Civita connection is denoted by ∇, the Riemannian gradient by grad, and the Riemannian Hessian by Hess. For x ∈ M, we denote by expx : Tx M → M the Riemannian exponential map and by logx (y) ∈ Tx M the inverse exponential map at x, well defined for every y ∈ M since Hadamard manifolds are globally geodesically convex. For α ∈ [0, 1] and x, y ∈ M, the notation x#α y = (1 − α)x ⊕ αy := expx α logx y

denotes the point at proportion α along the geodesic from x to y. We write Fix(T ) for the set of fixed points of a map T : M → M, i.e., Fix(T ) := {x ∈ M : T (x) = x}. We denote by O(n) the set of n × n real orthogonal matrices.

2

Background definitions and results

In this section we present some definitions and results that we use in Section 3 and Section 4. We first define Jacobi fields, which are variation fields of one-parameter families of geodesics [41, Chapter 10]. Carrying information about the curvature of the space, they provide a useful tool to analyze the behavior of nearby geodesics. In particular, they can be used to obtain a condition under which geodesics do not spread apart; see Definition 4.1. The second part is dedicated to quasi-α-firmly nonexpansive maps, and concludes by deriving such a property for nonexpansive Riemannian gradient-descent maps. 2.1

Jacobi fields

Let γ : R → M be a unit-speed geodesic. A Jacobi field along γ is a vector field satisfying the Jacobi equation ∇2t J + R(J, γ̇)γ̇ = 0, where ∇t = ∇γ̇(t) . Any Jacobi field can be decomposed as J(t) = J ∥ (t) + J ⊥ (t) = (α + βt)γ̇(t) + J ⊥ (t) for some coefficients α, β ∈ R, with J ∥ (t) parallel to γ̇(t), and J ⊥ (t) orthogonal to γ̇(t), i.e. ⟨J ⊥ (t), γ̇(t)⟩ = 0

for all t ∈ R.

We refer the reader to [41, Chapter 10] for a more in-depth study. 2.2

Quasi-α-firmly nonexpansive maps

A map T : M → M is 1-Lipschitz, i.e., nonexpansive, if d(T (x), T (y)) ≤ d(x, y) for every x, y ∈ M. This constraint does not necessarily lead to a convergent fixed point iteration xk+1 = T (xk ), k ≥ 0. Banach’s fixed point Theorem ensures that if T satisfies d(T (x), T (y)) < βd(x, y), for every x, y ∈ M and a β ∈ (0, 1), then such an iteration converges to the unique fixed point of T . However, this contraction constraint can be overly restrictive in practice. In this subsection, we introduce the necessary mathematical framework to define a broader subclass of nonexpansive maps which leads to convergent fixed point iterations. We start by defining the notion of quasi-α-firmly nonexpansive map, and introducing some of the properties of this class of functions. Definition 2.1. [12, Definition 1] Let M be a Hadamard manifold and let α ∈ (0, 1). We say T : M → M is quasi-α-firmly nonexpansive if Fix(T ) is non empty, and for every y ∈ Fix(T ) one has d2 (T x, y) ≤ d2 (x, y) −

1−α 2 d (T x, x) α

for every x ∈ M. The class of quasi-α-firmly nonexpansive maps is constrained enough to ensure the convergence of the fixed point iterations, whenever a fixed point exists, as formalized in the following theorem. Theorem 2.2. [11, Corollary 6] Let M be a finite-dimensional Hadamard manifold, and T : M → M a nonexpansive map. If T is quasi-α-firmly nonexpansive then the iterates xn := T (xn−1 ) for n ∈ N converge strongly to some x∗ ∈ Fix(T ). Quasi-α-firm nonexpansiveness is closed under composition under some assumptions. Closeness under composition is a condition that simplifies the design of constrained neural networks, since it allows working at the layer level, rather than considering the whole architecture. This is also true for the nonexpansiveness property. We leverage this in the network design process presented in Section 5. Proposition 2.3. [12, Corollary 13] Let M be a Hadamard manifold. Let T1 : D1 → M where D1 ⊂ M, and for j = 2, . . . , m let Tj : Dj → M where Dj := {Tj−1 (x) : x ∈ Dj−1 }. If Tj is quasi-α-firmly nonexpansive with constant αj on Dj for j = 1, . . . , m, and Fix(Tm ◦ Tm−1 ◦ · · · ◦ T1 ) ⊂ D1 is nonempty,

then the composite map T := Tm ◦ Tm−1 ◦ · · · ◦ T1 is quasi-α-firmly nonexpansive on D1 with constant given recursively by ξ m−1 + ξm αm = for m ≥ 3, ξ m−1 ξm + ξ m−1 + ξm where ξ j :=

1−αj 1−αj ξ1 +ξ2 αj for j ≥ 2, ξj := αj for j ≥ 1, and α2 := ξ1 ξ2 +ξ1 +ξ2 .

On Hadamard manifolds, we can derive a more practical construction of some quasi-α-firmly nonexpansive maps. This is formalized in the following theorem. Theorem 2.4. Let M be a Hadamard manifold and let α ∈ (0, 1). Let Tα : M → M be defined as Tα (x) = (1 − α)x ⊕ αN (x) = expx (α logx (N (x))),

x ∈ M,

for a nonexpansive map N : M → M. If Fix(N ) ̸= ∅, then Tα is quasi-α-firmly nonexpansive and Fix(Tα ) = Fix(N ). Proof. The proof is similar to the Hilbert space identity for averaged operators, see e.g. [5, Proposition 2.2]. Let p ∈ Fix(N ) and x ∈ M. On Hadamard manifolds, the distance function is uniformly convex [3, Definition 2.1], hence d2 (Tα (x), p) ≤ (1 − α)d2 (x, p) + αd2 (N (x), p) − α(1 − α)d2 (x, N (x)). Now d(N (x), p) = d(N (x), N (p)) ≤ d(x, p) by nonexpansiveness of N , and d(x, Tα (x)) = αd(x, N (x)) by construction. Therefore, d2 (Tα (x), p) ≤ d2 (x, p) −

1−α 2 d (x, Tα (x)). α

By the injectivity of the Riemannian exponential and logarithm on Hadamard manifolds, we have that x = Tα (x)

⇐⇒

α logx (N (x)) = 0

⇐⇒ N (x) = x,

and hence Fix(N ) = Fix(Tα ). This concludes the proof. The result above allows us to model quasi-α-firmly nonexpansive maps by modeling nonexpansive ones. This is why we focus on building 1-Lipschitz neural networks. Such a construction will be evident in the SPD-inverse problem of Section 6.2. Additionally, for the gradient-descent-type layers we study in the upcoming sections, we have an even more compatible structure as pointed out by Definition 2.5. Theorem 2.5. Let M be a Hadamard manifold. Let V : M → R be of class C 1 , and consider the Riemannian gradient descent step Tτ (x) = expx (−τ grad V (x)),

Tτ : M → M.

Assume Tτ is nonexpansive for τ ∈ [0, τ̄ ] and Fix(Tτ ) ̸= ∅. Then, for every τ ∈ (0, τ̄ ), Tτ is quasi-α-firmly nonexpansive for any α ∈ [τ /τ̄ , 1). Proof. Fix τ ∈ (0, τ̄ ] and α ∈ (0, 1). Define Gα,τ (x) := (1 − α)x ⊕ αTτ (x). Since Tτ is nonexpansive and Fix(Tτ ) ̸= ∅ by assumption, Definition 2.4 implies that Gα,τ is quasi-α-firmly nonexpansive. On the other hand, we also have Gα,τ (x) = Tατ (x) for every x ∈ M. Thus Tατ is quasi-α-firmly nonexpansive for any τ ∈ (0, τ̄ ] and any α ∈ (0, 1). Now let τ ′ ∈ (0, τ̄ ) and choose α ∈ [τ ′ /τ̄ , 1). Then τ :=

τ′ ∈ (0, τ̄ ] . α

By the previous argument, Tατ = Tτ ′ is quasi-α-firmly nonexpansive. As τ ′ ∈ (0, τ̄ ) was arbitrary, renaming τ ′ as τ yields the claim.

To the best of our knowledge, the literature does not provide a result showing the nonexpansiveness of the gradient-descent map Tτ for a general class of potentials V , such as geodesically convex and L-smooth potentials. For this reason, in Section 3 we derive a more specialized construction of V that allows us to prove the nonexpansiveness of the corresponding Riemannian gradient-descent step. We then use such a result to design our nonexpansive manifold-valued neural networks. The analysis of the more general case remains open, and could expand the range of 1-Lipschitz manifold-valued networks one can design following our principles.

3

Busemann gradient-descent layers

In Rn equipped with the Euclidean metric, one can consider layers of the form x 7→ x − τ ∇V (x) = x − τ A⊤ σ(Ax + b),

x, b ∈ Rn , A ∈ Rn×n , τ ≥ 0, Rτ n where σ is a scalar activation function, V (x) = 1⊤ n φ(Ax + b), φ(τ ) = 0 σ(s)ds, and 1n ∈ R is a vector of ones. Assuming σ is non-decreasing and L-Lipschitz, the above map is 1-Lipschitz and α-averaged whenever τ ∈ (0, 2/(L∥A∥22 )), with α = τ ∥A∥22 L/2 [59]. We recall that T : Rn → Rn is α-averaged if T (x) = (1 − α)x + αN (x) for a 1-Lipschitz map N : Rn → Rn , as considered for the construction in Definition 2.4. In this section, we consider Riemannian descent gradient layers on Hadamard manifolds. We propose a parametrization of the potential V based on Busemann functions, thereby obtaining a concrete class of nonexpansive models on Hadamard manifolds for which Busemann functions and their Riemannian gradients can be efficiently implemented. We first provide the basic notions and properties about such a family of functions, and refer the reader to [20, Chapter II.8] for a more in-depth study. We then present the Busemann gradient-descent layers and show that they are 1-Lipschitz under a certain condition on the stepsize (Definition 4.1). Let p ∈ M, ξ ∈ ∂M a point at infinity and γ : [0, +∞) → M the (unique) unit-speed geodesic ray from p to ξ. Let b denote the Busemann function associated to γ defined by b(x) = lim (d(γ(s), x) − s). s→+∞

(1)

When we want to remark the associated geodesic ray or point at infinity, we interchangeably use the notations bγ and bξ . The existence of the limit in (1) is guaranteed by [20, Chapter II.8, Lemma 8.18]. The level sets of b are called horospheres: for t ≥ 0 Hγ(t) := {x ∈ M : b(x) = b(γ(t))} (2) is the horosphere at γ(t), a smooth hypersurface orthogonal to γ at γ(t).

Lemma 3.1. Let (M, g) be a Hadamard manifold. Every Busemann function b on M is geodesically convex, 1-Lipschitz, and of class C 2 . The above result is quite standard, and we include a proof in Section A. Definition 3.2 below shows that Busemann flows, i.e., Riemannian gradient flows of Busemann functions, generate unit-speed geodesics, often referred to as Busemann geodesics. This property plays a central role in showing Definition 4.1. Lemma 3.2. Let (M, g) be a Hadamard manifold, and b a Busemann function on M. Let x ∈ M, and let (Φt )t≥0 denote the Busemann flow given by d Φt (x) = − grad b(Φt (x)), dt

Φ0 (x) = x.

Then Φt (x) = expx (−t grad b(x)) for every t ≥ 0, and (t, x) 7→ Φt (x) is C 1 jointly in (t, x). Moreover, b(Φt (x)) = b(x) − t, and Hess bΦt (x) (grad b(Φt (x)), grad b(Φt (x))) = 0.

(3)

Proof. The flow line t 7→ Φt (x) is a unit-speed geodesic asymptotic to ξ, see [28, Section 2.3], with initial data Φ0 (x) = x, Φ̇0 (x) = − grad b(x). On the other hand, the curve t 7−→ expx (−t grad b(x)) =: η(t) is by definition the unique geodesic with initial conditions η(0) = x, η̇(0) = − grad b(x). By uniqueness of geodesics with prescribed initial data, the two curves must coincide for every t ≥ 0. Definition 3.1 guarantees that b is C 2 , hence the flow map (t, x) 7→ Φt (x) is C 1 jointly in (t, x). Consider now the derivative of b(Φt (x)) with respect to t: d b(Φt (x)) = ⟨grad b(Φt (x)), − grad b(Φt (x))⟩ = −∥ grad b(Φt (x))∥2 = −1. dt Integrating the above yields b(Φt (x)) = b(x) − t, while deriving we find d2 b(Φt (x)) = Hess bΦt (x) (Φ̇t (x), Φ̇t (x)) dt2 = Hess bΦt (x) (grad b(Φt (x)), grad b(Φt (x))) = 0, i.e., the Hessian of b vanishes in the directions parallel to γ̇ (i.e., to grad b(x)). We now use the theory of Busemann functions and their gradient flows to design nonexpansive neural networks on Riemannian manifolds. We briefly introduce these layers now, and then study their properties in Section 4.

Definition of the Busemann gradient-descent layers Let γ be a geodesic ray on M, b the Busemann function generated by γ, λ ≥ 0, β ∈ R, φ : R → R of class C 1 and consider the potential V (x) = φ(λb(x) + β),

grad V (x) = λφ′ (λb(x) + β) grad b(x).

(4)

The Busemann (Riemannian) gradient-descent layer is defined by τ ≥ 0.

Tτ (x) = expx (−τ grad V (x)),

(5)

It follows from Definition 3.2 that, for every fixed base point x ∈ M, Tτ (x) = Φa(b(x)) (x),

(6)

where a(r) := τ λφ′ (λr + β), Φ is the Busemann flow as in (3), and the variable flow time is evaluated pointwise through a(b(x)). We can interpret this as a time-reparametrization of Busemann gradient flows.

4

Nonexpansiveness of Busemann gradient-descent layers

In this section, we present the main theoretical result of our work, Definition 4.1, which provides a condition on the stepsize ensuring that the Busemann gradient-descent layer in Equation (5) is 1-Lipschitz. We present two versions of the proof. The first one relies only on basic properties of Busemann functions on Hadamard manifolds, and provides a strategy to prove the product-manifold extension presented in Definition 4.4. The second one uses the construction in, e.g., [4, Theorems 2.2,3.1], [35, Theorems 6,9] based on variations through geodesics and Jacobi fields, providing a more geometric interpretation. Theorem 4.1. Let (M, g) be a Hadamard manifold and b a Busemann function on M. Let λ ≥ 0, τ ≥ 0, β ∈ R, and φ : R → R be a C 1 function such that φ′ is globally Lipschitz continuous, i.e. φ ∈ C 1,1 . Assume φ′ (s) ≥ 0

for all s ∈ R,

and ′′

0 ≤ τ λ2 φ′′ (s) ≤ 2

for a.e. s ∈ R, ′

(7)

where φ denotes the a.e. derivative of the Lipschitz function φ . Then, the map Tτ in (5) is nonexpansive on M.

Figure 1. Sharpness example for a single Busemann layer on D2 . For V (x) = 12 ReLU(bP (x))2 , choose two active points Pi = −ri P on the diameter opposite P . The update Tτ (x) = Φτ ReLU(bP (x)) (x) reparametrizes the Busemann flow by the value of bP , so the point with larger bP moves farther and the two trajectories cross. On this pair, d(Tτ (P1 ), Tτ (P2 ))/d(P1 , P2 ) = |1 − τ |, showing nonexpansiveness for 0 ≤ τ ≤ 2 and failure beyond the bound.

Before proving the theorem, we provide an illustrative example of the bound and its tightness in Figure 1. More details on the calculations are in Appendix B. Proof of Definition 4.1 (version 1). If τ = 0, or λ = 0, then Tτ = Id, and there is nothing to prove. Hence, we assume τ λ > 0 and that φ is not constant. Define a(r) := τ λφ′ (λr + β). Since φ′ is Lipschitz, the scalar map a is Lipschitz on R. Moreover, a′ (r) = τ λ2 φ′′ (λr + β) for a.e. r ∈ R. Hence 0 ≤ a′ (r) ≤ 2 for a.e. r ∈ R whenever τ satisfies (7). Let γ : [0, 1] → M be the (unique) length-minimising geodesic segment from x to y, i.e., γ(0) = x, γ(1) = y, and ℓ(γ) = d(x, y). Then, since d(Tτ (x), Tτ (y)) ≤ ℓ(Tτ ◦ γ), it suffices to show that ℓ(Tτ ◦ γ) ≤ ℓ(γ). Define w(s) := b(γ(s)). Since (t, x) 7→ Φt (x) is C 1 , the curve s 7→ Φa(w(s)) (γ(s)) = Tτ (γ(s)) is absolutely continuous. Therefore, it is enough to show that ∥(Tτ ◦ γ)′ (s)∥ ≤ ∥γ̇(s)∥ for a.e. s ∈ [0, 1] so that Z 1 Z 1 ′ ℓ(Tτ ◦ γ) = ∥(Tτ ◦ γ) (s)∥ ds ≤ ∥γ̇(s)∥ds = ℓ(γ). 0

0

It follows from Definition 3.2 that Tτ (γ(s)) = Φa(b(γ(s))) (γ(s)), where (Φt )t≥0 is the Busemann flow in (3). Since w is C 1 and a is Lipschitz, the curve a ◦ w is absolutely continuous. By the one-dimensional chain rule for Lipschitz functions [42, Corollary 3.68], (a ◦ w)′ (s) = a′ (w(s))w′ (s) for a.e. s ∈ [0, 1]. Therefore, for almost every s ∈ [0, 1], one has (Tτ ◦ γ)′ (s) = DΦa(w(s)) (γ̇(s)) + a′ (w(s)) w′ (s)

∂ Φt (γ(s)) . ∂t t=a(w(s))

We have w′ (s) = ⟨grad b(γ(s)), γ̇(s)⟩, i.e. w′ (s) is the component of γ̇(s) in the direction of grad b(γ(s)). For every s ∈ [0, 1], let ξH (s) ∈ span{grad b(γ(s))}⊥ be the orthogonal component so that γ̇(s) = w′ (s) · grad b(γ(s)) + ξH (s).

Then, ∥γ̇(s)∥2 = |w′ (s)|2 + ∥ξH (s)∥2 . By the properties of pushforward maps, DΦt (grad b(·)) = grad b(Φt (·)) and hence DΦa(b(γ(s))) (grad b(γ(s))) = grad b(Φa(b(γ(s))) (γ(s))) = grad b(Tτ (γ(s))). Similarly, ∂ = − grad b(Φa(b(γ(s))) (γ(s))) = − grad b(Tτ (γ(s))). Φt (γ(s)) ∂t t=a(b(γ(s))) Putting everything together, we have (Tτ ◦ γ)′ (s) = DΦa(b(γ(s))) (γ̇(s)) − a′ (w(s))w′ (s)DΦa(b(γ(s))) (grad b(γ(s)))   = DΦa(b(γ(s))) γ̇(s) − a′ (w(s))w′ (s) grad b(γ(s)) . Since Φ̇t = X ◦ Φt with X = − grad b, and since b ∈ C 2 is convex, we have ⟨∇η X, η⟩ = −⟨∇η grad b, η⟩ = −Hess b(η, η) ≤ 0. Thus X satisfies the Riemannian monotonicity condition of [35] with constant ν = 0. Hence its exact flow is nonexpansive. Since X is C 1 , the flow is differentiable with respect to the initial condition, and the corresponding differential satisfies ∥DΦt (x)η∥ ≤ ∥η∥,

t ≥ 0, η ∈ Tx M.

Applying this estimate with x = γ(s) and t = a(b(γ(s))) gives DΦa(b(γ(s))) (γ(s))η ≤ ∥η∥,

(8)

provided a(b(γ(s))) ≥ 0. This is why we assume a ≥ 0, equivalently φ′ ≥ 0. Therefore, ∥(Tτ ◦ γ)′ (s)∥ ≤ γ̇(s) − a′ (w(s))w′ (s) grad b(γ(s)) . Using the orthogonal decomposition of γ̇(s) at the source point, γ̇(s) − a′ (w(s))w′ (s) grad b(γ(s)) = ξH (s) + (1 − a′ (w(s)))w′ (s) grad b(γ(s)). Since ξH (s) ⊥ grad b(γ(s)), we obtain γ̇(s) − a′ (w(s))w′ (s) grad b(γ(s))

2

= ∥ξH (s)∥2 + |w′ (s)|2 (1 − a′ (w(s)))2 .

Finally, since 0 ≤ a′ (r) ≤ 2 for a.e. r ∈ R, it follows that γ̇(s) − a′ (w(s))w′ (s) grad b(γ(s))

2

≤ ∥ξH (s)∥2 + |w′ (s)|2 = ∥γ̇(s)∥2 ,

which implies ∥(Tτ ◦ γ)′ (s)∥ ≤ ∥γ̇(s)∥. This concludes the proof. Remark 4.2. The first version of the proof applies under the stated lower-regularity assumptions φ ∈ C 1,1 . The second proof is included as a more geometric argument, and is written under the additional simplifying assumption φ ∈ C 2 , so that a(r) = τ λφ′ (λr + β) is C 1 and the computations involving a′ are classical. Proof of Definition 4.1 (version 2). Define a(r) := τ λφ′ (λr + β).

In this proof, we assume φ ∈ C 2 . Then a′ (r) = τ λ2 φ′′ (λr + β) for every r ∈ R. Moreover, a(r) ≥ 0 and since φ′ is Lipschitz, the scalar map a is Lipschitz on R. Fix x, y ∈ M and let σ : [0, 1] → M be the minimizing geodesic from x to y, so that σ(0) = x,

∥σ̇(s)∥ ≡ d(x, y).

σ(1) = y,

Define a two-parameter map Γ : [0, 1] × [0, 1] → M,

Γ(s, t) := Φta(b(σ(s))) σ(s) .

For each fixed s, the curve t 7→ Γ(s, t) is the Busemann geodesic through σ(s), starting at Γ(s, 0) = σ(s) and ending at Γ(s, 1) = Φa(b(σ(s))) (σ(s)) = expσ(s) (−a(b(σ(s))) grad b(σ(s))) = Tτ (σ(s)). On the other hand, for each fixed t, the curve s 7→ Γ(s, t) is a deformation of σ into the image curve η(s) := Γ(s, 1) = Tτ (σ(s)). For each s ∈ [0, 1], the variational field along t 7→ Γ(s, t) defined by J s (t) := ∂s Γ(s, t) is a Jacobi field. In particular, J s (0) = ∂s Γ(s, 0) = σ̇(s),

J s (1) = ∂s Γ(s, 1) = η̇(s).

Our goal is to show that Z 1 0

∥J s (1)∥ ds ≤

Z 1 0

∥J s (0)∥ ds,

(9)

which implies d(Tτ (x), Tτ (y)) ≤ d(x, y). Fix s = 0 and consider the t-curve γ(t) := Γ(0, t),

γ(0) = x, γ(1) = expx (−a(b(x)) grad b(x)) = expx (−τ grad V (x)).

By construction, γ is a (reparametrized) Busemann geodesic γ(t) = Φta(b(x)) (x), γ̇(t) = a(b(x))

d = −a(b(x)) grad b(γ(t)), Φs (x) ds s=ta(b(x))

with ∥γ̇∥ = |a(b(x))| independent of t, and ∇γ̇ γ̇ = 0, hence γ̇ is parallel along γ. Set u(t) := grad b(γ(t)). Since ∇grad b grad b = 0, the vector field u is parallel along γ. We decompose J(t) = J ∥ (t) + J ⊥ (t) = (c0 + c1 t)u(t) + J ⊥ (t),

⟨J ⊥ (t), u(t)⟩ = 0.

The coefficients c0 and c1 depend on the initial data at t = 0. The torsion-free commutation identity gives Dt ∂s Γ(s, t) = Ds ∂t Γ(s, t), hence Dt J(0) = Dt ∂s Γ(s, t)|t=s=0 = Ds (−a(b(σ(s))) grad b(σ(s)))|s=0   d =− a(b(σ(s))) grad b(σ(s)) − a(b(σ(s))) Ds (grad b(σ(s)))|s=0 ds s=0 = −a′ (b(x)) ⟨grad b(x), S0 ⟩ grad b(x) − a(b(x))∇S0 grad b(x).

Then c0 and c1 are given by [35, Section 3] c0 = ⟨J(0), u(0)⟩ = ⟨S0 , grad b(x)⟩

c1 = ⟨Dt J(0), u(0)⟩ = −a′ (b(x)) ⟨S0 , grad b(x)⟩ = −a′ (b(x)) c0 ,

where we have used ⟨∇S0 grad b, grad b⟩ = Hess b(S0 , grad b)

= Hess b(grad b, S0 ) = ⟨∇grad b grad b, S0 ⟩ = 0,

since ∇grad b grad b = 0. Therefore ∥J ∥ (1)∥2 − ∥J ∥ (0)∥2 = (2c0 c1 + c21 ) = c20 a′ (b(x)) (−2 + a′ (b(x))) ≤ 0, as 0 ≤ a′ (r) ≤ 2 for any r ∈ R by assumption. Since u is parallel along γ, Dt J ⊥ = Dt J(t) − c1 u(t). Using Dt J = Ds ∂t Γ s=0 and ∂t Γ(s, t) = −a(b(σ(s))) grad b(Γ(s, t)) we obtain Dt J(t) = −

d u(t) − a(b(x))∇J(t) grad b(γ(t)) a(b(σ(s))) ds s=0

= c1 u(t) − a(b(x))∇J(t) grad b(γ(t)). Consequently, Dt J ⊥ (t) = −a(b(x))∇J(t) grad b(γ(t)). Since J(t) = (c0 + c1 t) grad b(γ(t)) + J ⊥ (t) and ∇grad b grad b = 0, we have ∇J(t) grad b(γ(t)) = ∇J ⊥ (t) grad b(γ(t)). Hence Dt J ⊥ (γ(t)) = −a(b(x))∇J ⊥ (t) grad b(γ(t)). Together with the convexity of b by Definition 3.1, we then obtain d ⊥ ∥J (t)∥2 = 2 ⟨∇γ̇(t) J ⊥ (t), J ⊥ (t)⟩ dt = −2a(b(x)) Hess bγ(t) (J ⊥ (t), J ⊥ (t)) ≤ 0, since a(r) ≥ 0 for any r ∈ R by assumption. In particular, ∥J ⊥ (1)∥2 − ∥J ⊥ (0)∥2 ≤ 0. Combining the two estimates and using the orthogonality of the decomposition, we obtain  
∥J(1)∥2 − ∥J(0)∥2 = ∥J ∥ (1)∥2 − ∥J ∥ (0)∥2 + ∥J ⊥ (1)∥2 − ∥J ⊥ (0)∥2 ≤ 0. Repeating the same reasoning to all s ∈ [0, 1] gives ∥J s (1)∥2 − ∥J s (0)∥2 ≤ 0 (equivalently ∥J s (1)∥ − ∥J s (0)∥ ≤ 0 ) for any s ∈ [0, 1]. Hence
d Tτ (x), Tτ (y) ≤ and Tτ is nonexpansive.

Z 1 0

s

∥J (1)∥ ds ≤

Z 1 0

∥J s (0)∥ ds = L(σ) = d(x, y),

Remark 4.3 (Combination or composition of potentials). Let ωi ≥ 0, λi ≥ 0, and βi ∈ R. Let bi be Busemann functions on M and let φ : R → R of class C 1,1 and such that ess sups∈R φ′′ (s) ≤ M2 . Consider the maps Ti (x) = expx (−τi grad Vi (x)),

Vi (x) = wi φ(λi bi (x) + βi ),

i = 1, . . . , m.

If τi wi λ2i M2 ≤ 2 for every i then every Ti is nonexpansive, and the composition T = Tm ◦ Tm−1 ◦ · · · ◦ T1 is nonexpansive too. Consider now the simultaneous step ! m X x 7→ expx −τ grad Vi (x) . i=1

Here, several Busemann directions are mixed inside one exponential, destroying the one-dimensional horospherical structure used above. The split layer keeps each factor in the single-Busemann regime, where the gradient step is a variable-time gradient flow and can be analyzed more easily. However, with τ1 = · · · = τm = τ small, the split version provides an approximation of the simultaneous step (Lie-Trotter splitting, [16]), making the split version as expressive as the more complex simultaneous-step one, while still allowing nonexpansiveness certificates. The map Tτ is a nonexpansive map on the finite-dimensional Hadamard manifold M. Moreover, Definition 2.5 with τ̄ given by (7) guarantees Tτ is quasi-α-firmly nonexpansive on M, and Definition 2.2 implies that the iterates xk+1 = Tτ (xk ), x0 ∈ M, converge strongly to a fixed point of Tτ in M, provided Fix Tτ ̸= ∅. Furthermore, this property is preserved under composition, provided we assume non-trivial intersections of the fixed sets; see Definition 2.3. This allows us to recover quasi-α-firmly nonexpansive networks in two different ways: by composing layers as Tτ with suitably constrained stepsizes and ensuring nontrivial intersections of their fixed sets, or defining a modified network Nθ (x) = x#α Dθ (x) for a nonexpansive network Dθ (x) = Tτm ◦ · · · ◦ Tτ1 and α ∈ (0, 1). In the inverse problem considered in Section 6.2, we opt for the second option since it is not immediate to ensure the desired condition on the fixed sets of the layers. 4.1

Generalization to product manifolds

The construction seen so far can be extended to product manifolds for which the base manifold admits efficiently computable Busemann functions and their gradients. This extension can be used to apply the proposed theory to hyperbolic-valued images and graphs, or for considering tasks on product manifolds, such as in Diffusion Tensor Imaging (DTI). Such applications are outside the scope of this paper and will be considered in future extensions. Let N = Mm be endowed with the product Riemannian metric. We adopt the notation dM and dN for the Riemannian distances induced on M and N , respectively. Let X = (x1 , . . . , xm ), Y = (y 1 , . . . , y m ) ∈ N denote a generic pair of points on the product manifold. The product distance is dN (X, Y )2 =

m X k=1

dM (xk , y k )2 .

For a tangent vector (v1 , . . . , vm ) ∈ Tx1 M × · · · × Txm M ≃ TX N , the product exponential is componentwise:
expX (v 1 , . . . , v m ) = expx1 (v 1 ), . . . , expxm (v m ) . Fix one Busemann function b : M → R and define the map B : N → Rm ,

B(X) = b(x1 ), . . . , b(xm ) .

Let U : Rm → R and define V (X) = U (B(X)). Assuming U ∈ C 1 , the product gradient has components gradxk V (X) = ∂k U (B(X)) grad b(xk ),

k = 1, . . . , m.

Hence the explicit Riemannian gradient step of V is T (X) = expX (−τ grad V (X)), with components
T (X)k = expxk −τ ∂k U (B(X)) grad b(xk ) = Φτ ∂k U (B(X)) (xk ),

(10)

where we have used Definition 3.2. Note that the flow times for each component depend on the full vector B(X), so the components of T are coupled through U ◦ B. We study the nonexpansiveness of this gradient step with the two results below. Theorem 4.4. Let (M, g) be a Hadamard manifold and N = Mm the product manifold equipped with the product metric. Let a : Rm → Rm be Lipschitz continuous as a vector-valued function. Assume ak (z) ≥ 0,

k = 1, . . . , m,

z ∈ Rm ,

and assume that R(z) := z − a(z) is nonexpansive in the Euclidean norm on Rm . Define T : N → N,

T (X)k = Φak (B(X)) (xk ),

k = 1, . . . , m.

Then T is nonexpansive on N . The above Definition 4.4 can be proved by adapting version 1 of the proof of Definition 4.1. We provide the complete proof in Section A. Remark 4.5. For m = 1, nonexpansiveness of r 7→ r − a(r) is equivalent to |1 − a′ (r)| ≤ 1 a.e., hence to the one-dimensional criterion 0 ≤ a′ (r) ≤ 2 a.e. in Definition 4.1. Therefore, Definition 4.4 provides a strict generalization of Definition 4.1. Corollary 4.6. Let (M, g) be a Hadamard manifold and N = Mm the product manifold equipped with the product metric. Let U : Rm → R be C 1 , convex, L-smooth, and nondecreasing in every component, i.e., ∇U (z) ≥ 0 componentwise. Define V (X) = U (B(X)). Then the explicit Riemannian gradient step T (X) = expX (−τ grad V (X)) is nonexpansive on N for every τ ∈ [0, 2/L]. We recall that a C 1 function U : Rm → R is L−smooth if its gradient ∇U : Rm → Rm is L−Lipschitz. Proof. For convex L-smooth U , the Euclidean map z 7→ z − τ ∇U (z) is nonexpansive for τ ∈ [0, 2/L]. We conclude by applying Definition 4.4 with a(z) = τ ∇U (z). 5

Implementation of the Busemann gradient-descent layers

In this section, we provide some practical implementations of the proposed Busemann gradient-descent layers. We provide explicit expressions of Busemann functions in the three settings M = Rn , the n−dimensional Poincaré ball M = Dn , and the manifold S++ (n) of symmetric positive definite (SPD) matrices. We then list some admissible activation functions (i.e., those that satisfy the hypothesis of Definition 4.1) and compute the corresponding bound on the stepsize. The proposed architectures on the Poincaré ball and SPD manifold are then empirically tested in Section 6. 5.1 5.1.1

Examples of Busemann functions Euclidean space

For every u ∈ Rn with ∥u∥2 = 1, the Busemann function associated with the geodesic ray γ(t) = −tu is b(x) = u⊤ x (see [20, Example II.8.24(1)]). The norm constraint can be removed by absorbing the norm in

λ: if ξ ∈ Rn \ {0}, then ξ ⊤ x = ∥ξ∥2 bξ/∥ξ∥2 (x) =: λbξ/∥ξ∥2 (x). Let β ∈ R and φ : R → R in C 1,1 be a convex function. We consider the potential V (x) = φ(ξ ⊤ x + β). Set ∆ := ess sups∈R φ′′ (s), which is finite given that φ′ is Lipschitz continuous. The function V defined above is L-smooth with L = ∥ξ∥22 ∆. The gradient-descent map associated to it is therefore 1-Lipschitz if τ ∈ [0, 2/L] and it writes Tτ (x) = x − τ ∇V (x) = x − τ φ′ (ξ ⊤ x + β)ξ.

(11)

In this context, θ = (ξ, β) are the trainable weights of the layer. Layers defined as (11) can recover Householder reflection layers [47, Eq. 6] and the MaxMin and GroupSort activations [1, 51]. We remark that, in the case of Rn , the potential obtained by summing finitely many such potentials also leads to nonexpansive gradient steps: V (x) =

k X i=1

φ(ξi⊤ x + βi ), Tτ (x) = x − τ Ξ⊤ φ′ (Ξx + v),

 ⊤ ξ1  ..  Ξ =  . , ξk⊤



 β1   v =  ...  , βk

 τ ∈ 0,

 2 . ∆∥Ξ∥22

This is a layer of a 1-Lipschitz ResNet as considered in [46, 59]. Remark 5.1. For M = Rn , the considered Busemann functions are linear. Therefore, here, the nondecreasing assumption on φ can be dropped. It suffices to assume that φ ∈ C 1,1 and convex. 5.1.2

Hyperbolic manifolds

In the Poincaré ball Dn = {x ∈ Rn : ∥x∥2 < 1}, the points at infinity are the ones lying on the hypersphere S n−1 = ∂Dn . Let ξ ∈ S n−1 and consider a geodesic ray γ, parametrized by arc length, and tending to ξ. The Busemann function associated with γ is   ∥ξ − x∥22 , (12) b(x) = lim (d(γ(t), x) − t) = log t→∞ 1 − ∥x∥22 where ∥ξ∥2 = 1 and x ∈ Dn (see [33]). The Euclidean gradient of b is ∇b(x) =

2 2 x+ (x − ξ). 1 − ∥x∥22 ∥x − ξ∥22

The Riemannian metric on Dn is defined as  D gx =

2 1 − ∥x∥22

2

g E =: δx2 g E

where g E is the Euclidean one. The Riemannian gradient grad b is then grad b(x) =

1 1 − ∥x∥22 (1 − ∥x∥22 )2 ∇b(x) = x+ (x − ξ). 2 δx 2 2∥x − ξ∥22

Using b and grad b as above yields an implementable gradient-descent map, with the freedom to choose ξ ∈ S n−1 , λ > 0, and β ∈ R. We provide a more implementation-oriented description of this layer in Algorithm 2.

5.1.3

SPD matrices

Consider the manifold of symmetric positive definite matrices S++ (n) = {X ∈ Rn×n : X ⊤ = X, X ≻ 0} equipped with the affine-invariant Riemannian distance [55] 

dAI (X, Y ) = Log X −1/2 Y X −1/2

 F

n X

=

!1/2 log2 λi

,

i=1

where λ1 , . . . , λn are the eigenvalues of X −1/2 Y X −1/2 . Then S++ (n) is a Hadamard manifold, and TX S++ (n) ≃ S(n), the space of n × n symmetric matrices. For the split SPD layers, we use the explicit Cholesky-based expression of Busemann functions on S++ (n). This is the formula appearing, for example, in [30, Example 4]. Let U ∈ O(n), let d ∈ Rn be a vector with ∥d∥2 = 1 and whose entries are sorted in ascending order, fixing the ordering convention used in the Cholesky factorization. We write U ⊤ XU = LL⊤ , where L is lower triangular factor in the Cholesky factorization of U ⊤ XU , which has positive diagonal. We define the Busemann function bU,d (X) = −2

n X

di log Lii .

i=1

A geodesic ray starting from the identity is determined by a unit tangent direction A ∈ S(n): γA (t) = expI (tA) = exp(tA),

∥A∥F = 1,

where exp without a subscript denotes the matrix exponential. We parametrize such directions as A = U diag(d)U ⊤ ,

U ∈ O(n),

∥d∥2 = 1, d1 < d2 < · · · < dn .

In the implementation, U is learned through an orthogonal parametrization. The P vector d is obtained by sorting a raw trainable vector, subtracting its mean, and normalizing it. Thus i di = 0 and ∥d∥2 = 1. This restricts the associated direction A to be trace-free. This centering is a modeling choice for the SPD denoiser, and it is not required by the general Busemann formula on full S++ (n). The affine-invariant Riemannian gradient of bU,d is grad bU,d (X) = −U LDL⊤ U ⊤ ,

D = diag(d),

see [30, Equation 23]. For a scalar potential V (X) = φ(λbU,d (X) + β),

λ > 0,

a negative-gradient step has tangent direction −τ grad V (X) = s U LDL⊤ U ⊤ ,

s = τ λ φ′ (λbU,d (X) + β).

In our experiments φ(t) = 21 ReLU(t)2 , so φ′ (t) = ReLU(t). Using the affine-invariant exponential map and congruence invariance, expLL⊤ (sLDL⊤ ) = L exp(sD)L⊤ . Equivalently, the Cholesky factor is updated as L 7−→ L exp

s  D , 2

which gives the closed-form split update   s    s ⊤ expX (−τ grad V (X)) = U L exp D L exp D U ⊤. (13) 2 2 This is the update implemented in the split Busemann SPD denoiser, which is more computationally efficient than directly computing the affine-invariant Riemannian exponential at X. The trainable parameters of each single-Busemann step are the orthogonal matrix U , the vector d, the scalar offset β, the positive scale λ, and the bounded stepsize τ . 5.2

Examples of activations

Our theory is compatible with activation functions φ : R → R that are C 1 , have a globally Lipschitz derivative, are convex, and are non-decreasing. Some admissible ones are: • Softplus. φ′ (r) =

φ(r) = log(1 + er ),

er ≥ 0, 1 + er

and condition (7) becomes 0 ≤ τ λ2 ≤ 8 ⇐⇒ 0 ≤ τ ≤

φ′′ (r) ≤

1 , 4

8 . λ2

• Squared ReLU. 1 2 φ(r) = ReLU(r) , 2

( ′

φ (r) = ReLU(r),

′′

φ (r) =

0, 1,

r < 0, r > 0,

with φ′′ defined for almost every r since φ′ is Lipschitz continuous. Condition (7) is then 2 0 ≤ τ λ2 ≤ 2 ⇐⇒ 0 ≤ τ ≤ 2 . λ • Smooth bounded-slope activations. Any convex nondecreasing C 1,1 activation with ess supr∈R φ′′ (r) ≤ M2

fits the theorem and yields the stepsize restriction 0 ≤ τ λ2 M2 ≤ 2 ⇐⇒ 0 ≤ τ ≤

2 λ2 M

. 2

In our numerical experiments, we consider ReLU2 /2, which has empirically performed slightly better than softplus. We plot the two activations side-by-side in Figure 2 together with their first two derivatives. Remark 5.2. The above are admissible activation functions under the C 1,1 -regularity assumption. If C 2 regularity is required, then ReLU2 (x)/2 would not be allowed.

6

Numerical experiments

This section evaluates the gradient-descent-type nonexpansive architectures we propose in two complementary regimes. The first experiment is a supervised classification problem on the Poincaré disk D2 . It is designed to test whether the constrained Busemann layers can fit non-linear hyperbolic decision boundaries while being robust under adversarial perturbations. The second experiment is an inverse problem on S++ (10), the manifold of 10 × 10 symmetric positive definite matrices, where a trained Busemann denoiser is inserted into a Plug-and-Play (PnP) reconstruction loop for masked Wishart observations. This second experiment tests whether the nonexpansive denoiser can act as an effective learned prior beyond a likelihood-only reconstruction. In both cases, we compare against baseline models that preserve the underlying geometry but lack Lipschitz constraints. We provide more details on each of the two problems and the considered architectures in the dedicated subsections below. The code accompanying the paper can be found at the GitHub repository https://github.com/davidemurari/one-lipschitz-hadamard-networks.

ϕ0

ϕ

ϕ00

value

ReLU2 /2

softplus

1.0

1.0

0.5

0.5

0.0

0.0 −2.5

0.0

2.5

−2.5

0.0

2.5

Figure 2. Two activation functions that satisfy the assumptions of our theory, along with their first two derivatives.

6.1

Hyperbolic classification on D2

We use two synthetic datasets on the Poincaré disk D2 . We sample the points first in the Euclidean plane and then map them to D2 through exp0 . The first class is sampled uniformly in B0.45,ℓ2 (0). The second class is sampled in B0.95,ℓ2 (0) \ B0.62,ℓ2 (0), with radius distributed as a Gaussian centered at 0.78 and standard deviation 0.08, clipped to [0.62, 0.95]. Here, we write Br,ℓ2 (0) for the closed Euclidean ball in R2 of radius r > 0 and center at the origin. The second dataset consists of k = 12 angular classes arranged in radial sectors. Half of the classes concentrate near the Euclidean disk radius rin = 0.7 and the remaining classes near rout = 0.8. Each point is obtained by perturbing the associated reference radius and angle with Gaussian additive noise with zero mean and fixed radial standard deviation σrad = 0.02 and angular standard deviation σang = 0.16. Both datasets contain 4000 points, split into an 80/20 train–test split. These tasks are intentionally simple enough to visualize, and they allow us to inspect the mathematical behavior of the models we are considering. We show the plot of the two datasets in Figure 3.

Figure 3. Datasets used for the classification problem.

Architectures We consider three models for this task. Two are 1-Lipschitz before the classification head, and the other is unconstrained. We now provide a precise description of how they are assembled. An algorithmic description of the Busemann and ResNet architectures is provided in Algorithm 1. The isometric model follows a similar logic. All the models start with a trivial embedding of the data points from D2 to D3 defined as x 7→ (x, 0). D3 is the manifold where the hidden maps are defined. The classifier head of all the models is a hyperbolic prototype head with C learned prototypes, where C is the number of classes. More precisely, we learn C

Algorithm 1 Schematic comparison of the two main hyperbolic classifiers. Require: Input x ∈ D2 Nonexpansive Busemann model 3

1:

Hyperbolic ResNet baseline

z ← (x, 0) ∈ D trivial isometric embedding

z ← (x, 0) ∈ D3 trivial isometric embedding

for r = 1, 2 do

for r = 1, 2 do

z ← Ar (z) Learnable isometry on D3

z ← Ar (z) Learnable Möbius affine map

for j = 1, . . . , 5 do

z ← expz (fr (z))

z ← Gbr,j (z) Busemann step, Algorithm 2

unconstrained residual step (15)

end for end for

end for

z ← A3 (z) Learnable isometry on D3

z ← A3 (z) Learnable Möbius affine map

ℓ ← H(z) ∈ RC classification head (14)

ℓ ← H(z) ∈ RC classification head (14)

return ℓ

return ℓ

target points t1 , . . . , tC ∈ D3 and return a vector in RC defined as  D3 ∋ h 7→ −d(h, t1 )

...

⊤ −d(h, tC ) ,

(14)

where h is the output of the last hidden layer. This head is designed so that the highest output entry corresponds to the closest prototype. 1-Lipschitz Busemann model. The core constrained model uses a split Busemann feature map. We always use ReLU2 /2 as the activation for the Busemann layers. Each hidden split block is a composition of singleBusemann updates as defined in Section 5.1.2, and in the implementation every such block is preceded by its own learnable hyperbolic isometry. The isometries are defined as D3 ∋ y 7→ b ⊕D3 Qy ∈ D3 where Q ∈ O(3) is a learnable orthogonal matrix, b ∈ D3 is a learnable bias, and ⊕D3 stands for the Möbius addition in D3 , see [32]. The learnable isometries are included because they preserve the network 1-Lipschitz structure and introduce additional flexibility into the model. The layer parameters are constrained so that the layer remains compatible with the nonexpansive construction. We consider two blocks of Busemann steps, each containing 5 gradient-descent steps and one isometry. 1-Lipschitz isometric model. This model is constrained as well, and provides a 1-Lipschitz baseline useful to judge the complexity of the task. It is obtained using the same structure as the Busemann model, but replacing the blocks of Busemann gradient-descent steps with an additional learnable isometry of D3 . Unconstrained hyperbolic ResNet. As a second baseline, we use an unconstrained hyperbolic residual network with the same structure as the other two models. The isometric layers are replaced by unconstrained Möbius affine maps, see [32, Equations (27)–(28)]. The Busemann blocks are replaced by unconstrained Poincaré residual steps defined as x 7→ expx (τ PT0→x (W2 σ(W1 log0 x + b1 ) + b2 )) , τ > 0,

(15)

where PT0→x is the parallel transport from the origin to x, log0 and expx are the logarithmic and exponential maps at the origin and x, respectively, and σ is a scalar nonlinearity. We use σ = ReLU and learn matrices W1 , W2 ∈ R3×3 .

Algorithm 2 One hyperbolic single-Busemann step Require: x ∈ Dd , raw direction praw , parameters λ, β, τ , scalar potential φ 1: p ← praw /∥praw ∥2 ∥x−p∥2 2: b ← bp (x) = log 1−∥x∥22 2 3: g ← ∇bp (x) 4: s ← τ λ φ′ (λb + β) 5: return expx (−sg)

Training and latent dynamics All models are trained by minimizing the multiclass cross-entropy loss, which, on one batch B of size |B|, writes L(θ) =

1 |B|

X

CE(fθ (xi ), yi ) ,

(xi ,yi )∈B

where fθ (x) denotes the vector of class scores obtained through the classifier head. We train all the models for 200 epochs with Adam, batch size 256, a one-cycle learning-rate schedule from 5 · 10−4 to 5 · 10−3 , no weight decay, and a validation set consisting of 10% of the training set. We also clip the gradient norm to 1 while training. For each model seed, the checkpoint with the best validation accuracy is retained and then evaluated on the held-out test set. The reported comparison uses the three model seeds 7, 11, 17. We report in Figure 4 the ensemble decision regions of the three models over the two classification problems.

Figure 4. For each model and each dataset, we have three seeds. We present the decision regions derived from the probability vector obtained by averaging the predicted probabilities across the three seeds. To extract probabilities from the logits, we use the softmax function. For the annular dataset, we include the circular separatrix obtained as exp0 (∂B(r1 +r2 )/2,ℓ2 (0)), which is the ideal decision boundary given our dataset.

For each model trained with seed 7, Figure 5 shows the evolution of the data points across the successive feature map layers. The plots make the different geometric mechanisms explicit. The isometric model can only relocate the data by constrained hyperbolic isometries. The Busemann model progressively pulls selected regions through nonexpansive Busemann-gradient steps, and the ResNet baseline uses unconstrained residual deformations. The end figure for each model represents, with squares, the learned prototype points t1 , t2 ∈ D3 (see (14)). Adversarial attacks and robustness evaluation Adversarial examples are generated by Riemannian projected gradient ascent in geodesic balls around each test point, using the projected-gradient attack paradigm

Figure 5. Layer-wise evolution of the annular test data through the feature maps for model seed 7. Rows correspond to the Isometry, Busemann, and ResNet models, while columns show the initial embedding, the representation after the first block, and the final representation. The Isometry model can only move the data via hyperbolic isometries, whereas the Busemann model performs structured nonexpansive deformations through compositions of Busemann-gradient steps. The ResNet baseline is unconstrained and can produce sharper residual deformations. The square markers in the last column show the pairs of learned prototypes.

and the multiclass cross-entropy loss. Each ascent step maps the Euclidean gradient to the Poincaré Riemannian gradient, moves with the exponential map, and projects back to the geodesic ball if necessary. For each test point, we run the projected-gradient ascent procedure several times from independently sampled initial points in the geodesic ball B(x0 , ε) := {x ∈ D2 : d(x0 , x) ≤ ε}. We refer to these independent initializations as restarts. A point is counted as robustly classified at radius ε if all restarts fail to find a misclassified point in B(x0 , ε). We use perturbation radii ε ∈ {0, 0.1, 0.2, 0.3, 0.4}, with 200 attack iterations, 20 random restarts, and stepsize 3ε/200. The robustness curves in Figure 6 measure worst-case behavior under repeated local searches in hyperbolic geodesic balls.

Interpretation The two datasets highlight complementary aspects of the constructions. On the annular task, the isometric baseline is not expressive enough: its mean clean accuracy is only about 60.8%. In contrast, both the Busemann model and the ResNet reach essentially perfect clean accuracy. Their robustness curves are also close to the ideal classifier provided by the intermediate circle, which we call oracle. Their normalized Areas Under the Curves (AUCs)1 are 88.6 for the Busemann model, 88.2 for the ResNet, and 90.0 for the oracle. Thus, in this task, the Busemann layers add geometric expressivity beyond isometries while preserving the nonexpansive structure. In other words, the Busemann layers play a fundamental role in achieving the final test accuracy. The angular-sector task gives a different message. Here the isometric baseline is already competitive with 1 This is a common quantity used to measure the tradeoff between accuracy and robustness. It is the area under the curve computed using the trapezoidal formula, and divided by 0.4, the total length of the epsilon range.

(a) Annular dataset

(b) Angular dataset

Figure 6. Robust accuracy as the adversarial perturbation radius ε increases. Solid curves show the mean over model seeds 7, 11, 17, and the shaded bands show the corresponding min–max envelope across the three seeds. Semantic oracle denotes the reference circle representing the ideal classification boundary in the annular dataset.

the Busemann model, indicating that this dataset can largely be solved by moving the prototype geometry with hyperbolic isometries. Nevertheless, the constrained models are more stable than the unconstrained ResNet under attack. The ResNet obtains higher clean accuracy, about 92.7%, but its robust accuracy drops to about 38.7% at ε = 0.4. The constrained models have lower clean accuracy, about 88.7%, but retain about 62.6% robust accuracy at the same radius. This supports the intended role of the nonexpansive constraint. It does not by itself guarantee higher clean accuracy, and it does not always make the Busemann model more expressive than isometries. However, it biases the learned feature map toward more stable geometry under hyperbolic perturbations. 6.2

Masked-Wishart inverse problem on S++ (10)

The second experiment studies the reconstruction of an unknown covariance matrix X⋆ ∈ S++ (n) from noisy masked principal-block observations. We consider the case n = 10. Each of the L masks ℓ = 1, . . . , L selects a coordinate set Iℓ ⊂ {1, . . . , n} of p = |Iℓ | coordinates. We denote by Pℓ ∈ Rp×n the corresponding coordinate-selection matrix, so that Pℓ x = xIℓ ∈ Rp . If x ∼ N (0, X⋆ ), then the observed subvector Pℓ x is Gaussian with covariance Cℓ (X⋆ ) = Pℓ X⋆ Pℓ⊤ . For each mask, we draw m = 20 independent samples yℓ,1 , . . . , yℓ,m ∼ N (0, Cℓ (X⋆ )) and observe the empirical covariance m 1 X ⊤ Sℓ = yℓ,r yℓ,r . m r=1 Equivalently, mSℓ follows the Wishart law Wp (Cℓ (X⋆ ), m) [50, 64]. For a candidate covariance X, the likelihood measures how compatible the predicted masked covariance Cℓ (X) is with the observed empirical covariance Sℓ . Since mSℓ ∼ Wp (Cℓ (X), m), the negative log-likelihood of the observed sample covariances, up to additive constants independent of X, is

mX log det Cℓ (X) + tr Cℓ (X)−1 Sℓ . (16) F (X) = 2 ℓ

This objective is the data-consistency term. It favors covariances whose observed principal blocks explain the empirical covariances Sℓ . The masks are fixed overlapping principal blocks. In the reported experiment, we set n = 10, block size p = 5, and L = 3, with coordinate sets I1 = {1, . . . , 5},

I2 = {3, . . . , 7},

I3 = {6, . . . , 10}.

Observed covariance entries 0 1 2

row

3 4 5 6 7 8 9 0

1

2

3

4

5

6

7

8

9

column

Figure 7. Observed-entry pattern induced by the three overlapping principal blocks. Colored cells are entries that appear in at least one observed 5 × 5 block.

Thus Pℓ XPℓ⊤ extracts a 5 × 5 principal submatrix of X. The induced observed-entry pattern is shown in Figure 7. The target covariances are sampled from a synthetic AR(1)-ensemble family. We use the stationary Gaussian AR(1) covariance model [21]: for autocorrelation parameter ρ, the covariance matrix has Toeplitz form Σ(ρ)ij = ρ|i−j| , i, j = 1, . . . , n. Each target is generated by drawing i.i.d. ρ1 , ρ2 , ρ3 ∼ U(0.2, 0.95) and setting 3

X⋆ =

1X Σ(ρq ). 3 q=1

Since only masked blocks are observed and the observations are noisy, Riemannian gradient descent with objective (16) alone does not encode the global covariance structure of the AR(1)-ensemble family. We therefore solve the reconstruction problem with a Plug-and-Play (PnP) scheme that alternates an affineinvariant gradient step for F with a learned SPD-valued denoiser, described below.

Denoisers We train two SPD-valued denoisers on full Wishart sample-covariance observations generated from the same AR(1)-ensemble family. The supervised training set contains 500 clean covariances. For each one, a full Wishart sample covariance is generated and used as the noisy input. The training task maps this full noisy covariance estimate to its clean covariance target. More precisely, for each training target (j) (j) X⋆ , j = 1, . . . , 500, we draw m = 20 independent samples zj,1 , . . . , zj,m ∼ N (0, X⋆ ) and form m

S (j) =

1 X ⊤ zj,i zj,i . m i=1

(j)

The denoiser is trained to map S (j) to X⋆ . The first model is the constrained split Busemann denoiser Dθ : S++ (10) → S++ (10). It comprises six blocks, each consisting of a learnable congruence isometry X 7→ QXQ⊤ , with Q ∈ O(10), followed by nine affine-invariant Busemann gradient steps. The Busemann potentials use the scalar nonlinearity ReLU2 /2. For more details, see Section 5.1.3. We enforce the stepsize constraint in (7) to ensure that the denoising map is nonexpansive in the affine-invariant metric. The second model is a Log-Euclidean residual denoiser Dθ (X) = exp (log X + sym(Rθ (log X))) ,

sym(A) =

A + A⊤ . 2

Here Rθ is a feedforward ReLU network applied in logarithmic coordinates. In the implementation, log X is flattened into an n2 -dimensional vector, passed through the network, reshaped back into an n × n matrix, and then symmetrized before applying the matrix exponential. For n = 10, the network has width 100 and

two affine-ReLU hidden stages before the output layer. This baseline is SPD-valued by construction, but it is not constrained to be nonexpansive in the affine-invariant geometry. Both denoisers are trained for 100 epochs in float64, then frozen for the Plug-and-Play reconstruction.

Plug-and-Play reconstruction Given an iterate Xk , the PnP update first applies one affine-invariant likelihood step for the masked-Wishart objective, Zk = expXk (−τ grad F (Xk )) , where the affine-invariant Riemannian gradient is computed by converting the closed-form expression of the Euclidean gradient:
mX ⊤ ∇F (X) = Pℓ Cℓ (X)−1 − Cℓ (X)−1 Sℓ Cℓ (X)−1 Pℓ , 2 ℓ

grad F (X) = X∇F (X)X. The matrix Zk is then averaged with the denoiser output along the affine-invariant geodesic, Xk+1 = Zk #α Dθ (Zk ), where A#α B = A1/2 A−1/2 BA−1/2

α

(17) A1/2 .

This is the weighted geometric mean or geodesic convex combination of positive definite matrices [15, Eq. (6.18)]. If α = 1, the map relies entirely on the learned denoiser. We remark that when Dθ is nonexpansive, Definition 2.4 applies, and ensures that the map in (17) is quasi-α-firmly nonexpansive for α ∈ (0, 1), and hence iterating only such steps without the data-fidelity step would ensure convergence to a fixed point, in case it exists (see Definition 2.2). This is confirmed empirically in Figure 11. The full PnP loop additionally contains the data step, so the convergence theory is not used as a proof for the inverse-problem iteration. Still, the empirical results support the validity of the proposed strategy. We use 32 validation observation-target pairs and 200 held-out test pairs using the masked-Wishart model above. The parameters τ , α, and the stopping iterate T are selected on validation by minimizing mean dAI (XT , X⋆ ), and then fixed for test evaluation. Figure 8 shows the validation grid used to select the PnP parameters for the two denoisers.

Baselines and metrics In addition to comparing the two denoisers, we benchmark different variations of the procedure. We use two types of baselines. First, we report static reconstructions that do not use the masked observations: the identity matrix, the Euclidean and Log-Euclidean means of the 500 training covariances used to train the denoisers. Second, we report a data-only baseline, obtained by affine-invariant gradient descent on the masked-Wishart objective F without any denoising step. The likelihood-only baseline selects its own data stepsize and stopping time using the same validation protocol as the PnP methods. The primary metric is the affine-invariant distance dAI (X, X⋆ ). We also record the masked-Wishart objective F (X), log-Euclidean error, relative Frobenius error, observed-entry error, and unobserved-entry error. Qualitative examples are selected using relative Frobenius error, because those panels visualize entrywise matrix errors. All quantitative tables and tuning decisions use the affine invariant distance. For the different procedures, we also test different initial guesses for the iterations. We start either from the
−1 PNtr identity matrix, the training Euclidean mean, or the training Log-Euclidean mean X LE = exp Ntr i=1 log Xi . Reconstruction results Table 1 shows that the training mean is already a strong prior baseline for this concentrated AR(1)-ensemble. The Euclidean training mean gives mean test error 0.845, and data-only descent from that mean improves it only to 0.789. The split Busemann PnP reconstruction reduces the error to about 0.368 from either training-mean initialization. This shows that this denoiser is stable to

Validation-selected mean dAI

Validation-selected mean dAI

0.2

0.2

0.15

0.15 3.0

0.12

0.12

0.1

3.5

0.1

0.08

0.08

2.5

3.0

0.01 0.005

1.5

0.003

0.03

2.5

0.02 0.01

2.0

0.005 0.003

0.002

1.5

0.002 1.0

0.001

0.001

1.0

0.0003 0.0001

averaging α

1

99

9

95

0.

0.

85

0.

0.

5

75 0.

25

0.

0.

1

15

05

0.5

0.

1

99 0.

9

95 0.

85

0.

0.

5

75 0.

25

0.

0.

1

15 0.

0.

0.

05

0.5

0.

0.0001

0.

0.0003

averaging α

Split Busemann

Log-Euclidean

Figure 8. Validation tuning for the PnP reconstruction. For each denoiser, the full validation grid is run over the tested initializations, data steps τ , averaging parameters α, and stopping times. The heatmap then shows the (τ, α) slice corresponding to the initialization with the best overall validation result for that denoiser. The star marks this global validation minimizer. Table 1 also reports the best validationselected parameters separately for each listed initialization. White cells correspond to non-finite or failed runs, which occur only for overly aggressive data steps and are not selected.

method identity Euclidean training mean Log-Euclidean training mean data-only descent split Busemann PnP split Busemann PnP split Busemann PnP Log-Euclidean PnP Log-Euclidean PnP Log-Euclidean PnP

initialization identity – – Euclidean mean identity Euclidean mean Log-Euclidean mean identity Euclidean mean Log-Euclidean mean

τ – – – 10−4 0.08 0.10 0.10 0.12 0.15 0.15

α – – – – 0.99 0.99 0.99 0.95 1.00 0.99

mean dAI ↓ 2.873 0.845 0.859 0.789 0.369 0.368 0.368 0.506 0.480 0.475

Table 1. Main test-set reconstruction results. All data-only and PnP parameters, including the stopping time, are selected based on the validation set and then frozen for test evaluation. The selected stopping times are T = 150 for data-only descent from the Euclidean mean, T = 5 for split PnP from identity, T = 50 for split PnP from either training mean, and T = 10 for the Log-Euclidean PnP rows. We tune α ∈ (0, 1], where α = 1 corresponds to applying the raw denoiser without geodesic relaxation. The averaged nonexpansive interpretation applies only when α ∈ (0, 1) and the denoiser is nonexpansive.

mean dAI

2.0

0.02

data step τ

0.05

0.03

mean dAI

0.05

data step τ

4.0

PnP model

mean ↑ improvement 0.477 0.491 0.422 0.364 0.379 0.309

baseline

split Busemann split Busemann split Busemann Log-Euclidean Log-Euclidean Log-Euclidean

Euclidean training mean Log-Euclidean training mean data-only, Euclidean mean init Euclidean training mean Log-Euclidean training mean data-only, Euclidean mean init

fraction ↑ improved 0.75 0.78 0.855 0.76 0.74 0.78

7

7

6

6

5

5

4

4

test dAI

test dAI

Table 2. Paired per-sample improvements on the test set. For each test covariance, we compute dAI (Xbaseline , X⋆ ) − dAI (XPnP , X⋆ ) using the same ground-truth target X⋆ . The “mean improvement” column averages this quantity over the test set. The “fraction improved” column is the fraction of test samples for which this paired difference is positive, i.e., for which the selected PnP reconstruction has smaller affine-invariant error than the baseline. For every row, XPnP is obtained by initializing the corresponding PnP iteration at the Euclidean training mean; its parameters τ , α, and T are those selected on validation for that initialization.

3

3

2

2

1

1

0

0

ty nti ide

an

euc

l

an ide

me

e

an ide ucl

log

an me

ly on ta nit Da ean i m

Split Busemann

id

P it Pn ty in i ent

P Pn init an me

ty nti

an

ide

ean lid

euc

me

log

l euc

an

ide

an me

ly on ta nit Da ean i m

P it Pn ty in nti

ide

Log-Euclidean

Figure 9. Per-sample test reconstruction errors corresponding to the methods and validation-selected configurations reported in Table 1. For each denoiser and initialization, the PnP parameters (τ, α, T ) are selected independently on the validation set and then fixed for test evaluation; the data-only parameters (τ, T ) are selected analogously. The boxplots summarize dAI (XT , X⋆ ) over the 200 held-out test instances, while the jittered points show the individual errors. The horizontal jitter is used only to improve visibility and has no quantitative meaning. The identity and training-mean entries are static baselines and require no parameter selection.

P Pn init an me

initialization. We conjecture that this is a consequence of the denoiser’s nonexpansiveness. The LogEuclidean denoiser also improves over data-only reconstruction, but its best PnP row is 0.475, which is weaker than the split Busemann result. Figure 9 shows the full test-error distributions, while Table 2 reports paired improvements on the same test instances. The split model improves over the data-only Euclideanmean baseline on 85.5% of test samples. Figure 12 gives representative correlation matrices for a qualitative view of these reconstruction differences. data-only identity data-only mean PnP identity

PnP mean Euclidean mean log-Euclidean mean

2.0

1.5

data-only identity data-only mean PnP identity

2.5

mean dAI (Xk , X? )

mean dAI (Xk , X? )

2.5

1.0

PnP mean Euclidean mean log-Euclidean mean

2.0

1.5

1.0

0.5 0.5 0

25

50

75

100

125

150

175

200

0

25

50

75

iteration

100

125

150

175

200

iteration

Split Busemann

Log-Euclidean

Figure 10. Mean target-error trajectories on the test set using the validation-selected stepsizes. At (b) (b) −1 PNtest iteration k, the vertical axis is Ntest b=1 dAI (Xk , X⋆ ). These curves are diagnostic rollouts. The reported table values use validation-selected stopping times, while the plots in this figure show how the same dynamics behave when continued longer. The keywords identity and mean refer to the initial guess X0 , which is either set to the identity or the Euclidean mean, respectively.

Averaged denoiser successive-iterate distance

Averaged denoiser successive-iterate distance 10−1

mean dAI (Xk+1 , Xk )

10−3

mean dAI (Xk+1 , Xk )

100

α = 0.25 α = 0.5 α = 0.75 α = 0.95

10−5 10−7 10−9

α = 0.25 α = 0.5 α = 0.75 α = 0.95

10−1

10−11 10−13

10−2 0

25

50

75

100

125

iteration

Split Busemann

150

175

200

0

25

50

75

100

125

150

175

200

iteration

Log-Euclidean

Figure 11. Averaged denoiser-only stability diagnostic. Each curve plots the mean successive-iterate distance dAI (Xk+1 , Xk ) for Xk+1 = Xk #α Dθ (Xk ). This isolates the averaged denoiser from the masked-Wishart data step, which is the component covered by the nonexpansive fixed-point theory.

Stability interpretation Figure 10 explains why early stopping is part of the reconstruction protocol. We can see that the validation-selected dynamics can initially reduce the test error but may deteriorate when continued beyond the selected stopping time. This is visible in the tuning heatmaps as failed cells for large τ , see Figure 8. The oscillations visible in the Log-Euclidean trajectory indicate that the Log-Euclidean PnP dynamics are more sensitive to early stopping, whereas the split Busemann dynamics remain more stable in this long-run diagnostic. Figure 11 isolates the denoiser-only averaged map X 7→ X#α Dθ (X). This is the part of the method for which the nonexpansive theory applies. In a Hadamard space, averaged nonexpansive iterations have fixed-point convergence guarantees when the fixed-point set is nonempty. The

decay of dAI (Xk+1 , Xk ) in the split-Busemann panel is therefore a diagnostic support for the intended averaged-map behavior. It is still not a convergence proof for the full PnP loop, because the full loop also contains the masked-Wishart gradient step. Reconstruction examples (correlation; rows selected by relative Frobenius error)

Reconstruction examples (correlation; rows selected by relative Frobenius error) |PnP − target|

1.0 0.5

0.10

0.0 0.05

best vs mean

−0.5

train mean

data-only

PnP

|PnP − target|

1.0 0.5

0.15

0.10

0.0 0.05

−0.5

−1.0

0.00

−1.0

0.00

1.0

0.15

1.0

0.15

0.5

0.10

0.0 0.05

−0.5

worst vs mean

target

0.15

median PnP

PnP

best vs mean

data-only

0.5

0.10

0.0 0.05

−0.5

−1.0

0.00

−1.0

0.00

1.0

0.15

1.0

0.15

0.5

0.10

0.0 0.05

−0.5 −1.0

worst vs mean

train mean

median PnP

target

0.5

0.10

0.0 0.05

−0.5 −1.0

0.00

Split Busemann

0.00

Log-Euclidean

Figure 12. Representative correlation reconstructions. Rows are selected by relative Frobenius error to make the entrywise visual comparison interpretable: median PnP error, best PnP improvement over the training mean, and worst PnP change relative to the training mean. Error panels use a common color scale [0, 0.15] across both denoisers.

observed unobserved

1.5

1.0

0.5

observed unobserved

2.0

relative Frobenius error

relative Frobenius error

2.0

1.5

1.0

0.5

0.0

0.0

ean

gm inin Tra

n mea only data

Split Busemann

se

PnP init ed lect

ean

m ing

in Tra

data

only

mea

n

PnP init cted

sele

Log-Euclidean

Figure 13. Observed-entry and unobserved-entry reconstruction errors for the two denoisers. The measured principal blocks determine only part of the covariance matrix. These diagnostics check whether the learned prior also improves entries not directly covered by the masks, rather than only fitting the observed blocks.

The observed/unobserved diagnostics in Figure 13 check that the improvement is not limited to entries directly included in the measured principal blocks. This is important because the masks underdetermine the inverse problem. Therefore, good performance on unobserved entries indicates that the learned denoiser is adding global covariance structure rather than only fitting the available blocks.

Conclusion for the SPD experiment The masked-Wishart experiment shows that a constrained split Busemann denoiser can be used as an effective learned prior in an SPD inverse problem. Under validationselected PnP inference, the split Busemann denoiser substantially improves over both baselines and over the Log-Euclidean SPD-valued denoiser. The theory motivates the averaged denoising component, while the complete data-plus-denoiser reconstruction is evaluated empirically using validation tuning and held-out test instances.

7

Conclusion

In this paper, we introduced a new neural network architecture on Hadamard manifolds. This model is practically implementable for manifolds where a computationally efficient closed-form expression for Busemann functions is available. We theoretically analyzed such an architecture, proving that it is 1-Lipschitz when considering small enough step sizes in the manifold-valued residual layers. The derived stepsize restriction is cheap to evaluate and practically implementable. We also extended such a theory to product manifolds, providing a viable theoretical framework for future extensions to manifold-valued graphs and images. The proposed methodology was then empirically validated by two numerical experiments. The first considers the task of robustly classifying points in the Poincaré disk D2 . The second considers an underdetermined covariance-matrix recovery problem. The core focus of this paper has been on developing the methodology and comprehensively studying its theoretical aspects from the perspective of stability and convergence. Future extensions will consider imageand graph-focused applications, such as for diffusion tensor imaging. Our architecture relies on the concept of Busemann functions. A promising further extension of our methodology is to consider alternative building blocks which allow for a similar design principle in a broader class of Riemannian manifolds, extending the availability of 1-Lipschitz manifold-valued neural networks to a wider range of Riemannian manifolds.

Acknowledgements DM acknowledges support from the EPSRC programme grant in ‘The Mathematics of Deep Learning’, under the project EP/V026259/1. CBS acknowledges support from the Royal Society Wolfson Fellowship, the EPSRC advanced career fellowship EP/V029428/1, the EPSRC programme grant EP/V026259/1, the Wellcome Innovator Awards 215733/Z/19/Z and 221633/Z/20/Z, the EPSRC funded ProbAI hub EP/Y028783/1. BA acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through grant RGPIN/2026-04531. All the authors acknowledge support from the EU through the Marie Skłodowska-Curie Actions Staff Exchanges project REMODEL, grant agreement 101131557.

REFERENCES [1] Cem Anil, James Lucas, and Roger Grosse. Sorting out lipschitz function approximation. In International conference on machine learning, pages 291–301. PMLR, 2019. [2] David Ariza-Ruiz, Laurenţiu Leuştean, and Genaro López-Acedo. Firmly nonexpansive mappings in classes of geodesic spaces. Transactions of the American Mathematical Society, 366(8):4299–4322, 2014. [3] David Ariza-Ruiz, Genaro López-Acedo, and Adriana Nicolae. The asymptotic behavior of the composition of firmly nonexpansive mappings. Journal of Optimization Theory and Applications, 167(2):409–429, 2015. [4] M. Arnold, E. Celledoni, E. Çokaj, B. Owren, and D. Tumiotto. B-stability of numerical integrators on Riemannian manifolds. Journal of Computational Dynamics, 11(1):92–107, 2024. [5] Miroslav Bacák. Convex analysis and optimization in Hadamard spaces, volume 22. Walter de Gruyter GmbH & Co KG, 2014. [6] Florian Bachmann, Ralf Hielscher, and Helmut Schaeben. Grain detection from 2d and 3d EBSD data—specification of the MTEX algorithm. Ultramicroscopy, 111(12):1720–1733, 2011. [7] Shaojie Bai, J. Zico Kolter, and Vladlen Koltun. Deep equilibrium models. In Advances in Neural Information Processing Systems, volume 32, 2019. [8] Justin Baker, Qingsong Wang, Cory D. Hauck, and Bao Wang. Implicit graph neural networks: A monotone operator viewpoint. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 1521–1548, 2023.

[9] Werner Ballmann. Lectures on Spaces of Nonpositive Curvature, volume 25 of Oberwolfach Seminars. Birkhäuser, Basel, 1995. With an appendix by Misha Brin. [10] H. Bauschke and P. Combettes. Convex analysis and monotone operator theory in Hilbert spaces. Springer, 2017. [11] Arian Bërdëllima. On a notion of averaged mappings in CAT(0) spaces. Functional Analysis and Its Applications, 56(1):27–36, 2022. [12] Arian Bërdëllima, Florian Lauster, and D Russell Luke. α-firmly nonexpansive operators on metric spaces. Journal of Fixed Point Theory and Applications, 24(1):14, 2022. [13] Ronny Bergmann, Friederike Laus, Johannes Persch, and Gabriele Steidl. Recent advances in denoising of manifold-valued images. In Handbook of Numerical Analysis, volume 20, pages 553–578. 2019. [14] G Pacelli Bessa, Jorge H de Lira, Stefano Pigola, and Alberto G Setti. Curvature estimates for submanifolds immersed into horoballs and horocylinders. Journal of Mathematical Analysis and Applications, 431(2):1000–1007, 2015. [15] Rajendra Bhatia. Positive definite matrices. Princeton Series in Applied Mathematics. Princeton University Press, Princeton, NJ, 2007. [16] Sergio Blanes, Fernando Casas, and Ander Murua. Splitting methods for differential equations. arXiv preprint arXiv:2401.01722, 2024. [17] Clément Bonet, Laetitia Chapel, Lucas Drumetz, and Nicolas Courty. Hyperbolic sliced-wasserstein via geodesic and horospherical projections. In Proceedings of 2nd Annual Workshop on Topology, Algebra, and Geometry in Machine Learning (TAG-ML), volume 221 of Proceedings of Machine Learning Research, pages 334–370, 2023. [18] Clément Bonet, Lucas Drumetz, and Nicolas Courty. Sliced-wasserstein distances and flows on cartanhadamard manifolds. Journal of Machine Learning Research, 26(32):1–76, 2025. [19] Nicolas Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, 2023. [20] Martin R Bridson and André Haefliger. Metric spaces of non-positive curvature, volume 319. Springer Science & Business Media, 2013. [21] Peter J. Brockwell and Richard A. Davis. Time Series: Theory and Methods. Springer, 1991. [22] Michael M. Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2104.13478, 2021. [23] E. Celledoni, D. Murari, B. Owren, C. Schönlieb, and F. Sherry. Dynamical systems-based neural networks. SIAM J. Sci. Comput., 45(6):A3071–A3094, 2023. [24] Ines Chami, Albert Gu, Dat P. Nguyen, and Christopher Re. Horopca: Hyperbolic dimensionality reduction via horospherical projections. In Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 1419–1429, 2021. [25] Moustapha Cissé, Piotr Bojanowski, Edouard Grave, Yann N. Dauphin, and Nicolas Usunier. Parseval networks: Improving robustness to adversarial examples. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 854–863, 2017. [26] Patrick L. Combettes and Jean-Christophe Pesquet. Lipschitz certificates for layered network structures driven by averaged activation operators. SIAM Journal on Mathematics of Data Science, 2(2):529–557, 2020. [27] Patrick L. Combettes and Jean-Christophe Pesquet. Fixed point strategies in data science. IEEE Transactions on Signal Processing, 69:3878–3905, 2021.

[28] Christopher Criscitiello and Jungbin Kim. Horospherically convex optimization on Hadamard manifolds part I: Analysis and algorithms. arXiv preprint arXiv:2505.16970, 2025. [29] Xiran Fan, Chun-Hao Yang, and Baba C. Vemuri. Horospherical decision boundaries for large margin classification in hyperbolic space. In Advances in Neural Information Processing Systems, volume 36, 2023. [30] OP Ferreira, DS Gonçalves, MS Louzeiro, SZ Németh, and J Zhu. A subdifferential characterization via busemann functions and applications to dc optimization on hadamard manifolds. arXiv preprint arXiv:2602.20931, 2026. [31] P. Thomas Fletcher, John Moeller, Jeff M. Phillips, and Suresh Venkatasubramanian. Horoball hulls and extents in positive definite space. In Algorithms and Data Structures, volume 6844 of Lecture notes in Computer Science, pages 401–412, 2011. [32] Octavian Ganea, Gary Bécigneul, and Thomas Hofmann. Hyperbolic neural networks. Advances in neural information processing systems, 31, 2018. [33] Mina Ghadimi Atigh, Martin Keller-Ressel, and Pascal Mettes. Hyperbolic busemann learning with ideal prototypes. Advances in neural information processing systems, 34:103–115, 2021. [34] Mina Ghadimi Atigh, Max van Spengler, Teng Long, Melika Ayoughi, Tejaswi Kasarla, and Pascal Mettes. Hyperbolic learning with supervision from any granularity, 2026. [35] Marta Ghirardelli, Brynjulf Owren, and Elena Celledoni. Conditional Stability of the Euler Method on Riemannian Manifolds. arXiv preprint arXiv:2503.09434, 2025. [36] Artur Gorokh, Yury Korolev, and Tuomo Valkonen. Diffusion tensor imaging with deterministic error bounds. Journal of Mathematical Imaging and Vision, 56(1):137–157, 2016. [37] Mehrtash T. Harandi, Mathieu Salzmann, and Richard Hartley. From manifold to manifold: Geometryaware dimensionality reduction for SPD matrices. In Computer Vision – ECCV 2014, pages 17–32, 2014. [38] Ernst Heintze and Hans-Christoph Im Hof. Geometry of horospheres. Journal of Differential geometry, 12(4):481–491, 1977. [39] Zhiwu Huang and Luc Van Gool. A riemannian network for spd matrix learning. In Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, pages 2036–2042, 2017. [40] Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Proximal denoiser for convergent plug-andplay optimization with nonconvex regularization. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 9483–9505, 2022. [41] J. M. Lee. Introduction to Riemannian manifolds, volume 176 of Graduate Texts in Mathematics. Springer, Cham, 2018. [42] Giovanni Leoni. A first course in Sobolev spaces. American Mathematical Soc., 2017. [43] Chong Li, Genaro López, and Victoria Martín-Márquez. Monotone vector fields and the proximal point algorithm on Hadamard manifolds. Journal of the London Mathematical Society, 79(3):663– 683, 2009. [44] Kanti V. Mardia and Peter E. Jupp. Directional Statistics. John Wiley & Sons, 2000. [45] Jonathan Masci, Davide Boscaini, Michael M. Bronstein, and Pierre Vandergheynst. Geodesic convolutional neural networks on riemannian manifolds. In 2015 IEEE International Conference on Computer Vision Workshop (ICCVW), pages 832–840, 2015. [46] Laurent Meunier, Blaise J Delattre, Alexandre Araujo, and Alexandre Allauzen. A Dynamical System Perspective for Lipschitz Neural Networks. In International Conference on Machine Learning, pages 15484–15500. PMLR, 2022.

[47] Zakaria Mhammedi, Andrew Hellicar, Ashfaqur Rahman, and James Bailey. Efficient orthogonal parametrisation of recurrent neural networks using householder reflections. In International Conference on Machine Learning, pages 2401–2409. PMLR, 2017. [48] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral normalization for generative adversarial networks. In International Conference on Learning Representations, 2018. [49] Jean-Jacques Moreau. Proximité et dualité dans un espace hilbertien. Bulletin de la Société mathématique de France, 93:273–299, 1965. [50] Robb J. Muirhead. Aspects of Multivariate Statistical Theory. Wiley, 1982. [51] Davide Murari, Takashi Furuya, and Carola-Bibiane Schönlieb. Approximation theory for 1-Lipschitz ResNets. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. [52] Xuan Son Nguyen and Shuo Yang. Building neural networks on matrix manifolds: A gyrovector space approach. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 26031–26062, 2023. [53] Maximillian Nickel and Douwe Kiela. Poincaré embeddings for learning hierarchical representations. Advances in neural information processing systems, 30, 2017. [54] Wei Peng, Tuomas Varanka, Abdelrahman Mostafa, Henglin Shi, and Guoying Zhao. Hyperbolic deep neural networks: A survey. IEEE Transactions on pattern analysis and machine intelligence, 44(12):10023–10044, 2021. [55] Xavier Pennec, Pierre Fillard, and Nicholas Ayache. A riemannian framework for tensor computing. International Journal of computer vision, 66(1):41–66, 2006. [56] Haifeng Qian and Mark N. Wegman. L2-nonexpansive neural networks. In International Conference on Learning Representations, 2019. [57] Ernest Ryu, Jialin Liu, Sicheng Wang, Xiaohan Chen, Zhangyang Wang, and Wotao Yin. Plug-andplay methods provably converge with properly trained denoisers. In International Conference on Machine Learning, pages 5546–5557. PMLR, 2019. [58] Frederic Sala, Chris De Sa, Albert Gu, and Christopher Ré. Representation tradeoffs for hyperbolic embeddings. In International conference on machine learning, pages 4460–4469. PMLR, 2018. [59] F. Sherry, E. Celledoni, M.J. Ehrhardt, D. Murari, B. Owren, and C. Schönlieb. Designing stable neural networks using convex analysis and ODEs. Phys. D, 463:Paper No. 134159, 13, 2024. [60] Ha-Young Shin. Radial fields on the manifolds of symmetric positive definite matrices. SIAM Journal on Applied Algebra and Geometry, 10(1):1–13, 2026. [61] Sho Sonoda, Isao Ishikawa, and Masahiro Ikeda. Fully-connected network on noncompact symmetric space and ridgelet transform based on helgason-fourier analysis. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 20405–20422, 2022. [62] Oncel Tuzel, Fatih Porikli, and Peter Meer. Region covariance: A fast descriptor for detection and classification. In Computer Vision – ECCV 2006, pages 589–600, 2006. [63] Andreas Weinmann, Laurent Demaret, and Martin Storath. Total variation regularization for manifoldvalued data. SIAM Journal on Imaging Sciences, 7(4):2226–2257, 2014. [64] John Wishart. The generalised product moment distribution in samples from a normal multivariate population. Biometrika, 20A(1/2):32–52, 1928.

A

Proofs omitted in the main document

We now provide a proof of Definition 3.1 after restating it. L EMMA 1. Let (M, g) be a Hadamard manifold. Every Busemann function b on M is geodesically convex, 1-Lipschitz, and of class C 2 . Proof of Definition 3.1. Let γ : [0, +∞) → M be a geodesic ray and bγ (x) its associated Busemann function given by (1). The characterization of horofunctions [20, Proposition 8.22] guarantees that the function bγ (x) is convex on Hadamard manifolds. By the triangle inequality, bγ is 1-Lipschitz. In fact, ∥∇bγ ∥ = 1 by the Gauss Lemma. In particular, bγ is differentiable a.e. and, on a Hadamard manifold, it is of class C 2 [9, 14, 38].

We now provide a proof of Definition 4.4 after restating it. T HEOREM 1. Let (M, g) be a Hadamard manifold and N = Mm the product manifold equipped with the product metric. Let a : Rm → Rm be Lipschitz continuous as a vector-valued function. Assume ak (z) ≥ 0,

z ∈ Rm ,

k = 1, . . . , m,

and assume that R(z) := z − a(z) is nonexpansive in the Euclidean norm on Rm . Define T (X)k = Φak (B(X)) (xk ),

T : N → N,

k = 1, . . . , m.

Then T is nonexpansive on N . Proof. Let X, Y ∈ N , and let Γ : [0, 1] → N be the minimizing geodesic joining them. Write Γ(s) = (γ 1 (s), . . . , γ m (s)). It is enough to show L(T ◦ Γ) ≤ L(Γ). For each k = 1, . . . , m, set rk (s) := b(γ k (s)),

r(s) := B(Γ(s)) = (r1 (s), . . . , rm (s)) ∈ Rm .

Since Γ is a smooth geodesic and b ∈ C 2 , the curve r is C 1 , hence Lipschitz on [0, 1]. Since R := I − a is nonexpansive, R ◦ r is Lipschitz. Therefore R ◦ r is differentiable for a.e. s ∈ [0, 1], and at such points ∥r(s + h) − r(s)∥2 ∥R(r(s + h)) − R(r(s))∥2 ≤ lim = ∥r′ (s)∥2 . h→0 h→0 |h| |h|

∥(R ◦ r)′ (s)∥2 = lim

Moreover, since a = I − R, the curve a ◦ r is Lipschitz as well. Hence a ◦ r is differentiable for a.e. s, and at points where both derivatives exist, (R ◦ r)′ (s) = r′ (s) − (a ◦ r)′ (s). For a.e. s, decompose each component velocity as k γ̇ k (s) = rk′ (s) grad b(γ k (s)) + ξH (s),

k ξH (s) ⊥ grad b(γ k (s)).

Since ∥ grad b∥ = 1, this gives ∥Γ̇(s)∥2 =

m X k=1

∥γ̇ k (s)∥2 = ∥r′ (s)∥22 +

m X k=1

k ∥ξH (s)∥2 .

Now set Ak (s) := ak (r(s)) = ak (B(Γ(s))). Then T (Γ(s))k = ΦAk (s) (γ k (s)). At a.e. s, using ∂t Φt = − grad b ◦ Φt , we obtain d T (Γ(s))k = DΦAk (s) (γ̇ k (s)) − A′k (s) grad b(T (Γ(s))k ). ds Since DΦt (grad b) = grad b ◦ Φt , we can rewrite this as
d T (Γ(s))k = DΦAk (s) γ̇ k (s) − A′k (s) grad b(γ k (s)) . ds Using the decomposition k γ̇ k (s) = rk′ (s) grad b(γ k (s)) + ξH (s)

and the identity rk′ (s) − A′k (s) = (R ◦ r)′k (s), we get
d k T (Γ(s))k = DΦAk (s) ξH (s) + (R ◦ r)′k (s) grad b(γ k (s)) . ds Since Ak (s) = ak (r(s)) ≥ 0, the differential of the Busemann flow ΦAk (s) is nonexpansive. Therefore d T (Γ(s))k ds

2 k ≤ ξH (s) + (R ◦ r)′k (s) grad b(γ k (s))

2

.

k (s) ⊥ grad b(γ k (s)) and ∥ grad b∥ = 1, this gives Because ξH 2

d T (Γ(s))k ds Summing over k gives d T (Γ(s)) ds

2

k ≤ ∥ξH (s)∥2 + |(R ◦ r)′k (s)| .

2

≤

m X k=1

k ∥ξH (s)∥2 + ∥(R ◦ r)′ (s)∥22 .

Using ∥(R ◦ r)′ (s)∥2 ≤ ∥r′ (s)∥2 , we obtain d T (Γ(s)) ds Thus L(T ◦ Γ) =

Z 1 0

2

≤

m X k=1

k ∥ξH (s)∥2 + ∥r′ (s)∥22 = ∥Γ̇(s)∥2 .

d T (Γ(s)) ds ≤ ds

Z 1 0

∥Γ̇(s)∥ ds = L(Γ).

Since Γ is minimizing, dN (T (X), T (Y )) ≤ L(T ◦ Γ) ≤ L(Γ) = dN (X, Y ). Therefore T is nonexpansive.

B

Calculations for tightness of the nonexpansiveness bound

We now provide an explicit example on D2 demonstrating that the bound in (7) is tight. Fix P ∈ S 1 , i.e. P ∈ R2 with ∥P ∥2 = 1. Let us consider two points Pi = −ri P,

0 < ri < 1, i = 1, 2,

and

 bP (x) = log

∥x − P ∥2 1 − ∥x∥2

 .

On the diameter through P , write any point as x = ρP , where ρ ∈ (−1, 1). Then     (1 − ρ)2 1−ρ bP (ρP ) = log = log = −2 tanh−1 (ρ). 1 − ρ2 1+ρ For Pi = −ri P , this gives

bP (Pi ) = bP (−ri P ) = 2 tanh−1 (ri ) > 0.

We then consider the potential

1 ReLU(bP (x))2 , 2 i.e., set β = 0 and λ = 1. The update it defines is V (x) =

(18)

Tτ (x) = Φτ ReLU(bP (x)) (x). At P1 and P2 , ReLU is active, so the flow time is si (τ ) = τ bP (Pi ) = 2τ tanh−1 (ri ). Let us now use the defining property of the Busemann flow bP (Φs (x)) = bP (x) − s. Let γi (τ ) := Tτ (Pi ). Since the flow stays on the diameter, we write γi (τ ) = ρi (τ )P . Then bP (γi (τ )) = bP (Pi ) − si (τ ) = 2 tanh−1 (ri ) − 2τ tanh−1 (ri ) = 2(1 − τ ) tanh−1 (ri ).
But also bP (ρi (τ )P ) = −2 tanh−1 (ρi (τ )). Therefore ρi (τ ) = tanh (τ − 1) tanh−1 (ri ) . We conclude that
γi (τ ) = Tτ (Pi ) = tanh (τ − 1) tanh−1 (ri ) P. The expansivity ratio to compute is now ℓ(τ ) := d(γ1 (τ ), γ2 (τ ))/d(P1 , P2 ). Since both Riemannian distances are between points of the form λ1 P and λ2 P , λ1 , λ2 ∈ (−1, 1), we remark that   λ1 P − λ 2 P −1 d(λ1 P, λ2 P ) = 2 tanh . 1 − λ1 λ2 ∥P ∥2 Since ∥P ∥2 = 1, we recover that −1

d(λ1 P, λ2 P ) = 2 tanh



λ1 − λ 2 1 − λ 1 λ2



= 2| tanh−1 (λ1 ) − tanh−1 (λ2 )|.

This allows us to recover d(γ1 (τ ), γ2 (τ )) = 2|τ − 1| · | tanh−1 (r1 ) − tanh−1 (r2 )| = |1 − τ |d(P1 , P2 ), and that ℓ(τ ) = |1 − τ |. Therefore, ℓ(τ ) ≤ 1 if and only if 0 ≤ τ ≤ τmax = 2, which is the theoretical bound predicted by (7) for the potential in (18).