Lucas De Lara

Alberto González-Sanz

Jean-Michel Loubes

Diffeomorphic Registration using Sinkhorn Divergences

Keywords: Diffeomorphic Registration, Entropic Optimal Transport, Matching Estimation

The diffeomorphic registration framework enables to define an optimal matching function between two probability measures with respect to a data-fidelity loss function. The non convexity of the optimization problem renders the choice of this loss function crucial to avoid poor local minima. Recent work showed experimentally the efficiency of entropy-regularized optimal transportation costs, as they are computationally fast and differentiable while having few minima. Following this approach, we provide in this paper a new framework based on Sinkhorn divergences, unbiased entropic optimal transportation costs, and prove the statistical consistency with rate of the empirical optimal deformations.

Introduction

Diffeomorphic deformations describe a large class of computational frameworks whose goal is to find optimal deformations of the ambient space, defined as a diffeomorphisms generated through flow equations [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF]. They amount to solving an optimization problem involving two terms: an objective loss function characterizing in which sense the deformation should be optimal; a penalization over the kinetic energy spent by the transformation. The versatility of the problem formulation along with the appealing mathematical properties of diffeomorphisms made diffeomorphic deformations widely used in various application fields. In particular, they have been popularized for diffeomorphic registration in medical image analysis. This task consists of constructing diffeomorphic matching functions between shapes in order to establish spatial correspondences [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. More recently, [START_REF] Younes | Diffeomorphic learning[END_REF] proposed to apply flows of diffeomorphisms in a machine-learning context, where the optimal deformation is designed to render the data classes linearly separable.

This paper focuses on the diffeomorphic registration problem between two shapes. More specifically, we address the setting where the shapes are represented by probability measures: a formulation that has received a growing interest over the past few years to address unlabeled landmarks [START_REF] Glaunes | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Bauer | Diffeomorphic density matching by optimal information transport[END_REF][START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF][START_REF] Feydy | Global divergences between measures: from hausdorff distance to optimal transport[END_REF]. In this case, the objective loss function, referred as the data-fidelity loss, is defined as a metric between probability measures. Squares of maximum mean discrepancies (MMD), which are well-known kernel-based distances, became the canonical choice for such settings. In particular, their use for diffeomorphic registration enjoys a well-established theory [START_REF] Glaunes | Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching[END_REF][START_REF] Glaunes | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Younes | Shapes and diffeomorphisms[END_REF]. However, they also suffer from important practical drawbacks.

As pointed out by [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF], the non convexity of the optimization problem on the diffeomorphic deformation renders the choice of the loss function crucial to avoid poor local minima, whereas an MMD possesses many. This is why they proposed to use optimal transport metrics as an alternative. More precisely, they define the data-fidelity loss as the entropy-regularized optimal transportation cost between unbalanced measures, which has two critical advantages. Firstly, it benefits from the non locality of optimal transport metrics, leading to few local minima. Secondly, entropic regularization alleviates the computational burden of standard optimal transport: it allows for fast computation and differentiation of the cost through the celebrated Sinkhorn's algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. Nevertheless, while this alternative loss for diffeomorphic registration performs better experimentally, it lacks the statistical theory that was proven for squares of MMDs. Moreover, the entropic regularization induces a well-known bias making the loss not minimal between two identical measures. The latter issue motivates the employment of a Sinkhorn divergence: a symmetric unbiased version of the standard entropy-regularized optimal transportation cost. In [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF], the authors showed that Sinkhorn divergences performed significantly better than their biased counterparts for registration purpose. However, they carried out their analysis using flows of gradients (an approach reviewed by [START_REF] Santambrogio | {Euclidean, metric, and Wasserstein} gradient flows: an overview[END_REF]) instead of flows of diffeomorphisms.

This paper addresses diffeomorphic registration for Sinkhorn-divergence-based fidelity losses from both a theoretical and practical viewpoint. By leveraging some recent advances on these divergences [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF], Genevay et al., 2019], we show in a statistically-driven approach that the deformation obtained by solving the optimization problem between empirical measures converges with the parametric rate √ n to its population counterpart, where n is the sample size. Additionally, we illustrate the practicality of our method through numerical experiments. This furnishes a new theoretically and practically grounded framework for diffeomorphic matching of probability measures.

Related work

Several papers bear resemblances with our work as they combine entropic optimal transport with diffeomorphic registration at some point of their pipeline. Let us underline the major differences with our framework. The work of [START_REF] Croquet | Unsupervised diffeomorphic surface registration and non-linear modelling[END_REF] leverages a Sinkhorn divergence as the data-fidelity loss of a regularized diffeomorphic-registration engine restricted to flows induced by stationary velocity fields (SVF), which are notoriously not tailored to match significantly different shapes [START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF]. In contrast, our approach applies to the more flexible large deformation diffeomorphic metric mapping (LDDMM) framework where the flows are time dependent. In [START_REF] Shen | Accurate point cloud registration with robust optimal transport[END_REF], the authors also interface entropic optimal transport with large diffeomorphic deformations but for a different role: optimal transport computes a prior landmark alignment instead of acting as the data attachment term. The closest approach to ours is the one of [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF] who first suggested to use entropic optimal transport as the data-fidelity loss for diffeomporphic registration. However, they relied on the biased transportation cost between unbalanced measures whereas we tackle the unbiased divergence between probabilities. Additionally, their work focuses on practical applications while we also provide theoretical background. Finally, one implementation of time-variant diffeomorphic registration driven by an unbiased Sinkhorn divergence can be found in a PhD manuscript [Feydy, 2020, Figure 4.6]. In our work, we go further by filling the theoretical gap, as well as by proposing more comprehensive experiments illustrating the behaviour of these loss functions.

Outline The rest of the paper is organized as follows. In Section 2, we specify the basic mathematical notations that will be used throughout the paper. In Section 3, we set up the general problem we address by introducing the diffeomorphic registration framework for arbitrary data-fidelity losses. In Section 4, we present the necessary background on optimal transport and entropic regularization, in order to properly define Sinkhorn divergences. Additionally, we study some indispensable regularity properties of entropic optimal transport. In Section 5, we state our main results, that is the existence and statistical consistency of the optimal deformations.

In Section 6, we recall the implementation of diffeomorphic registration, and present the numerical experiments where we benchmark Sinkhorn divergences with other losses. All the proofs are deferred to Appendix B, while Appendix A recalls key mathematical tools from empirical process theory and Frechet differentiability.

Preliminaries and notations

In this section, we introduce the definitions and notations that will be used throughout the paper. The first part is dedicated to classes of smooth functions; the second one addresses probability measures.

Smooth functions

Let d 1 ≥ 1 and X be an arbitrary subset of R d1 with non-empty interior denoted by X . For p ≥ 1 and d 2 ≥ 1, we define C p (X , R d2 ) as the set of p-continuously Frechet-differentiable functions from X to R d2 . We also define L p (R d1 , R d2 ) the set of symmetric p-multilinear operators from R d1 to R d2 . The p-th derivative of some

F ∈ C p (X , R d2 ) is denoted by F (p) . It maps any point x ∈ X to F (p) (x)[•] ∈ L p (R d1 , R d2 ). By convention we set F (0) = F . For any L ∈ L p (R d1 , R d2
), we define the operator norm as

L op := sup{ L[δ 1 , . . . , δ k ] | δ i ∈ R d1 , δ i ≤ 1} where • is the Euclidean norm. For example, if F ∈ C 1 (X , R), then F (x) op = ∇F (x)
where ∇F is the gradient of F . This enables to define, for any F ∈ C p (X , R d2 ), the functional norm,

F p,∞ := max 0≤k≤p F (k) ∞ ,
where F ∞ := sup x∈X F (x) , and F (k) ∞ := sup x∈ X F (k) (x) op for k ≥ 1. In addition, for any R > 0 we denote by C p R (X , R d2 ) the class of functions F ∈ C p (X , R d2 ) such that F p,∞ ≤ R, and write B R for the centered Euclidean ball of radius R.

Actions on probability measures

We write E[X] for the expectation of any random variable X. The symbol ⊗ denotes the product of measures. For two measures µ and ν on R d , the relation µ ν means that µ is absolutely continuous with respect to ν, that is

(ν(E) = 0 =⇒ µ(E) = 0) for every measurable set E ⊆ R d .
We define two kinds of actions involving probability measures. Let µ be a probability measure on R d and f : R d → R be a measurable function. The action of µ on f defines the real number:

µ(f ) := f dµ = E X∼µ [f (X)].
Now, consider a measurable function F : R d → R d . The action of F on µ defines a probability measure called the push-forward measure, defined as:

F µ := µ • F -1 (•).
If a random variable X follows the law µ, then the image variable F (X) follows the law F µ. The push-forward operation enables to write changes of variables. Formally,

f d(F µ) = (f • F )dµ.

Diffeomorphic measure transportation

In this section we present the necessary background on diffeomorphic registration of probability measures. We refer to [START_REF] Younes | Shapes and diffeomorphisms[END_REF] for a complete and precise treatment of this topic. Firstly, we recall how to define diffeomorphisms through flow equations. Secondly, we introduce the diffeomorphic measure transportation problem for arbitrary data-fidelity losses.

Generating diffeomorphic deformations

The diffeomorphic deformation framework can be framed as a fluid mechanics problem, where points in R d are transported by a vector field representing a stream varying across time in the ambient space. We begin by reviewing the corresponding formalism and theory.

For an integer p ≥ 1 let B p be the space of functions in C p (R d , R d ) whose derivatives up to order p vanish to zero at infinity. This together with the norm • p,∞ is a Banach space. Next, denote by V a Hilbert space with inner product •, • V and norm • V , and assume that V is continuously embedded in B p . This corresponds to the hypothesis below.

Assumption 3.1. The space V is included in B p , and there exists a constant c V > 0 such that for any v ∈ V ,

v p,∞ ≤ c V v V .
Physically, a function v ∈ V represents a stationary vector field in the ambient space, specifying the speed vector v(x) ∈ R d of the stream running at every position x ∈ R d . Then, define the class L 2

V of vector fields t ∈ [0, 1] → v t ∈ V indexed by time and space satisfying

1 0 v t 2 V dt < ∞, which is a Hilbert space endowed with the inner product, v, u L 2 V := 1 0 v t , u t V dt.
We recall that a sequence

{v n } n∈N in L 2 V converges weakly to v if for any u ∈ L 2 V , v n , u L 2 V -----→ n→+∞ v, u L 2 V . (1)
The associated norm in L 2 V is given by

v L 2 V := 1 0 v t 2 V dt,
and we use the notation

L 2 V,M := {v ∈ L 2 V | v L 2 V ≤ M }
for the centered ball of radius M > 0 in L 2 V . We can now turn to the definition of diffeomorphic deformations. Any vector field v ∈ L 2

V generates a deformation φ v := (φ v t ) t∈[0,1] , function of both time and space variables, defined as the unique solution to the following flow equation,

∀x ∈ R d , ∀t ∈ [0, 1], φ t (x) = x + t 0 v s φ s (x) ds.
(2)

The parametric curve (φ v t (x)) t∈[0,1] represents the trajectory across time of a point initially located at φ 0 (x) = x ∈ R d . Remarkably, for every t ∈ [0, 1] the transformation φ v t is a p-continuously differentiable diffeomorphism. Moreover, as a direct consequence of [Glaunes, 2005, Theorem 5], these diffeomorphic transformations are smooth over compact sets.

Lemma 3.1 (Smoothness of diffeomorphic deformations). Suppose that Assumption 3.1 holds. Then for any radius M > 0 and any compact set

K ⊂ R d , there exists a constant R = R((K, d); (V, p); M ) > 0 such that for any v ∈ L 2 V,M , max 0≤k≤p sup t∈[0,1],x∈K (φ v t ) (k) (x) op ≤ R.
In particular, sup x∈K x ≤ R.

In practice, the space of vector fields V is constructed through the choice of a kernel function. This is enabled by Assumption 3.1 which entails that V is a reproducing kernel Hilbert space (RKHS), characterized by a unique non-negative symmetric matrix-valued kernel function Ker : R d × R d → R d×d . In particular, the choice of the kernel function sets the order of regularity p of the vector fields. For instance, the typical choice of a Gaussian kernel, that is

Ker(x, y) := 1 √ 2πσ 2 exp - x -y 2 2σ 2 I d (3)
where σ > 0 is the bandwidth parameter and I d the identity matrix, leads to p = +∞.

Diffeomorphic matching of distributions

In general, diffeomorphic deformation frameworks amount to finding solutions to Equation (2) that are optimal in some sense. In this work, we focus on the diffeomorphic measure transportation framework, which aims at matching two probability measures. Formally, let Λ be a positive loss function between probability measures, and set α and β two probabilities on the ambient space R d . For a given regularization weight λ > 0, an optimal matching function between α and β is a diffeomorphism φ v solution to (2) where v minimizes

J λ (v) := Λ(φ v 1 α, β) + λ v 2 L 2 V . (4)
The first term of the objective function ( 4) is the data-fidelity loss, which tends to match φ v 1 α with β, while the second term is the regularizer, which penalizes the kinetic energy spent by the trajectories (φ v t ) t∈[0,1] , keeping them as close as possible to the identity function. The parameter λ governs the trade-off between the two contributions. The objective J λ always admits minimizers provided that the term v ∈ L 2 V → Λ(φ v 1 α, β) ∈ R + is weakly continuous. For a minimizer v * , the function φ v * 1 is an optimal matching between α and β, and the family (φ v * t ) t∈[0,1] provides an approximated interpolation between the two measures. In practical settings, one typically does not have access to the full probability measures α and β but to empirical observations. This naturally raises the question of estimating an optimal matching function between α and β on the basis of independent samples. Concretely, let x 1 , . . . , x n ∼ α and y 1 , . . . , y n ∼ β be independent samples, and define the empirical probability measures α n := n -1 n i=1 δ xi and β n := n -1 n j=1 δ yj . Plugging these discrete measures in the original objective function (4) leads to the following empirical objective function:

J λ,n (v) := Λ(φ v 1 α n , β n ) + λ v 2 L 2 V . (5)
In Theorem 5.1 we prove under some assumptions that if the data-fidelity loss Λ is a Sinkhorn divergence, a divergence derived from entropic optimal transport, then any sequence of minimizers {v n } n∈N of the empirical problem (5) converges up to the extraction of a subsequence to a minimizer of the population problem (4) as the sample size n increases to infinity.

Entropic optimal transport

In this section, we first briefly present the necessary background on optimal transport and entropic regularization, in order to properly define Sinkhorn divergences. We refer to [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal Transport: Old and New[END_REF], Peyré et al., 2019] for further insight on these topics. Then, we introduce some properties of these divergences, which will be useful to later demonstrate the main results of this paper.

Transportation costs and Sinkhorn divergences

Let α and β be two probability measures on X a subset of R d , and C : R d × R d → R + a positive ground cost function. Typically, C(x, y) := x -y 2 . The optimal transportation cost with respect to C between α and β is defined as,

T C (α, β) := min π∈Π(α,β) X ×X C(x, y)dπ(x, y), (6)
where Π(α, β) is the set of couplings admitting α as first marginal and β as second marginal. In particular, for an integer k ≥ 1 and D a distance on X , the quantity (T D k )

1 k yields a distance between measures referred as the Wasserstein distance of order k. Transportation costs and optimal transport distances became popular in many machine-learning-related problems for their appealing geometric properties, but suffer from being computationally challenging in practice. This triggered a growing literature on fast approximations of ( 6), the most popular being entropy-regularized versions, which can be computed through the Sinkhorn algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. For ε > 0, the entropy-regularized transportation cost w.r.t. C is defined as

T C,ε (α, β) := min π∈Π(α,β) X ×X C(x, y)dπ(x, y) + εKL(π|α ⊗ β), (7)
where KL(µ|ν) denotes the Kullback-Leibler divergence between probability measures µ and ν given by log dµ dν (z) dµ(z) if µ ν, and +∞ otherwise. Critically, the entropic transportation cost T C,ε suffers from the so-called entropic bias, that is T C,ε (α, α) = 0 in general. As illustrated in [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF], this entails that the minimum of T C,ε (α, •) is not reached at α but at a shrunken version of α with smaller support, making the entropic cost an unreliable loss function. The Sinkhorn divergence was originally introduced to fix this undesirable effect. It is formally defined as

S C,ε (α, β) := T C,ε (α, β) - 1 2 T C,ε (α, α) - 1 2 T C,ε (β, β).
As aforementioned, using a non-local similarity measure such as an entropic-optimal-transport cost instead of a local similarity measure such as a squared MMD leads to fewer local solutions when minimizing (5). Moreover, it does not suffer from the computational burden of standard optimal transport. This is why [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF] advocated the use of the entropy-regularized transportation cost (7) for diffeomorphic registration, providing empirical evidences of the benefits of this approach. However, they did not rely on the unbiased Sinkhorn divergences, for which little was known until [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF] that demonstrated several key properties. In particular, if C is continuous, e -C ε defines a positive universal kernel, and X is compact, then S C,ε is symmetric positive definite, smooth and convex in each of its input distributions. Additionally, in contrast to the standard regularized transportation cost, it metrizes the convergence in law. In particular, these properties hold for the classical cost functions C(x, y) := x -y and C(x, y) := x -y 2 defined on compact domains. The goal of this paper is precisely to use a Sinkhorn divergence for the data-fidelity loss, while providing statistical guarantees. The demonstrations are based on the dual formulation of entropic optimal transport for which we derive some important results next.

Regularity of the dual formulation

The minimization problem (7) has the following dual formulation,

T C,ε (α, β) = sup f,g∈C(X ,R) X f (x)dα(x) + X g(y)dβ(y) -ε X ×X e f (x)+g(y)-C(x,y) ε dα(x)dβ(y) + ε. (8)
The functions f and g are referred as potentials. Note that Equation ( 8) can also be compactly written as,

T C,ε (α, β) = sup f,g∈C(X ,R) (α ⊗ β) h f,g C,ε , where h f,g C,ε (x, y) := f (x) + g(y) -εe f (x)+g(y)-C(x,y) ε + ε. ( 9
)
We call the function h f,g C,ε the global potential. It will play a key role in the proofs. A remarkable property of entropic optimal transport, investigated in [START_REF] Genevay | Sample complexity of sinkhorn divergences[END_REF], Feydy et al., 2019], is that the potentials of the dual formulation inherit the regularity of the ground cost function C if the measures α and β are compactly supported. This setting will be useful to derive statistical guarantees. More specifically, it allows to restrict the set of feasible potentials to smooth functions regardless of the involved probability measures, as stated in the next lemma which readily follows from [Genevay et al., 2019, Proposition 1] (see also [del Barrio et al., 2022, Lemma 4.1] for the particular case of the quadratic ground cost).

Lemma 4.1 (Smoothness of the optimal potentials). Let µ and ν be two measures on a compact set K ⊂ R d , and suppose that the ground cost function

C belongs to C q (R d × R d , R + ) with q ≥ 1. Then, there exists a constant m = m((K, d); (C, q); ε) > 0 such that T C,ε (µ, ν) = sup f,g∈C(K,R) (µ ⊗ ν) h f,g C,ε = sup f,g∈C q m (K,R) (µ ⊗ ν) h f,g C,ε .
Naturally, the smoothness of f , g and C renders the global potential h f,g C,ε smooth as well. Combining Lemma 4.1 with the following result ensures the smoothness of the optimal global potential under smooth data-processing transformations, such as diffeomorphic transformations.

Proposition 4.1 (Smoothness of the optimal global potential). Let X be a compact subset of R d , suppose that the ground cost function C belongs to C q (R d × R d , R + ) with q ≥ 1, set p ≥ 1 and write κ := min{p, q}. Then for any m > 0 and R > 0, there exists a constant H = H(m; R; (C, q); ε; p) > 0 such that for any f, g

∈ C q m (B R , R) and T 1 , T 2 ∈ C p R (X , R d ), h f,g C,ε • (T 1 , T 2 ) ∈ C κ H (X × X , R).
We are now ready to state and prove our main results.

Main results

This section focuses on the main theoretical contributions of the paper, namely the existence and statistical consistency of the empirical optimal matching function between α and β when using a Sinkhorn divergence. Firstly, we show that the objective functions J λ and J λ,n with Λ = S C,ε admit minimizers. We recall that a function Ψ :

L 2 V → R is weakly continuous if for any sequence {v n } n∈N weakly converging to some v ∈ L 2 V (see (1)), we have Ψ(v n ) -----→ n→+∞ Ψ(v)
. [Glaunes, 2005, Theorem 7] states that J λ admits a minimum if

v ∈ L 2 V → Λ(φ v 1 α, β
) is weakly continuous and non negative while [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF], Theorem 1] guarantees the non negativeness of Sinkhorn divergences when e -C ε defines positive universal kernel. Therefore, existence of an optimal matching directly follows from the proposition below.

Proposition 5.1 (Existence of the optimal vector fields). Let α and β be two probability measures on X a compact subset of R d , suppose that the ground cost function C belongs to C 1 (R d × R d , R + ), and assume that Assumption 3.1 holds. Then the function

v ∈ L 2 V → S C,ε (φ v 1 α, β) is weakly continuous. If additionally e -C ε
defines a positive universal kernel, then J λ for Λ = S C,ε admit minimizers.

The minimizer is not unique in general due to the non convexity of the data-fidelity loss with respect to v. Uniqueness could be artificially achieved by choosing λ very large, thereby rendering the objective function strictly convex, but this would make the purpose of the regularization meaningless.

We now turn to our main theorem, which is divided in two items. The first one ensures the convergences of the empirical solutions to their population counterparts; the second one specifies the speed of this convergence.

Theorem 5.1 (Consistency of the optimal vector fields). Let α n and β n be empirical measures corresponding respectively to α and β, two probability measures on X a compact subset of R d , suppose that the ground cost function C belongs to C q (R d × R d , R + ) with q ≥ 1 and induces a positive universal kernel e -C ε . Finally, assume that Assumption 3.1 holds. If, for any n ∈ N * , v n denotes a minimizer of J λ,n for Λ = S C,ε , then the following results hold.

(i) There exists a minimizer of J λ denoted by v * such that up to the extraction of a subsequence

v n -v * L 2 V a.s. ----→ n→∞ 0 and sup t∈[0,1] φ v n t -φ v * t ∞ + (φ v n t ) -1 -(φ v * t ) -1 ∞ a.s. ----→ n→∞ 0. (ii) If κ := min{p, q} > d, then there exists a constant A = A(λ; (X , d); (C, q); ε; (V, p)) > 0 such that E [|J λ (v n ) -J λ (v * )|] ≤ A √ n .
Note that [START_REF] Glaunes | Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching[END_REF] proved a similar consistency result when the data-fidelity loss is the square of an MMD, but did not determine the speed of convergence as in (ii). The demonstration of (i) follows the steps of their proof (see [Glaunes, 2005, Theorem 16]). The idea is to show the convergence of sup v∈L 2 V,M |J λ,n (v) -J λ (v)| as n increases to infinity, where L 2 V,M contains all the minimizers independently of n. The main challenge when addressing an entropic optimal transport cost comes from the fact that it does not satisfy a triangle inequality, nor a data-processing inequality, and is hence harder to control. We remedy to this issue by proving and applying the following intermediary result: Proposition 5.2 (Uniform consistency of entropic optimal transport up to smooth data-processing transformations). Let α n and β n be empirical measures corresponding respectively to α and β, two probability measures on X a compact subset of R d , and suppose that the ground cost function C belongs to C q (R d × R d , R + ) with q ≥ 1. Set p ≥ 1 and write κ := min{p, q}. Then, the following results hold:

(i) For any R > 0 sup T1,T2∈C p R (X ,R d ) |T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β)| a.s.
-----→ n→+∞ 0.

(ii) If κ > d, then for any R > 0 there exists a constant

A = A(R; (C, q); ε; (X , d); p) > 0 such that E sup T1,T2∈C p R (X ,R d ) |T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β)| ≤ A √ n .
Notice that as a direct consequence of the triangle inequality, a similar result holds for S C,ε . Hence, as diffeomorphisms are smooth on compact sets according to Lemma 3.1, we can apply Proposition 5.2 to control sup v∈L 2 V,M |J λ,n (v) -J λ (v)|. Although Proposition 5.2 is motivated by diffeomorphic registration, we believe it has further interest. Remark in particular that the objective (4) shares similarities with generative modelling [START_REF] Goodfellow | Generative adversarial nets[END_REF]; an input distribution α is passed through a parametric function φ v 1 meant to generate a target distribution β by minimizing a certain loss Λ. In particular, generative modelling using the Wasserstein-1 distance or a Sinkhorn divergence has proved to be efficient for diverse applications [START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF][START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF]. The main difference in (4) comes from the parameter v being infinitely dimensional, and characterizing a diffeomorphism instead of a neural network. However, Proposition 5.2 is general enough to be applied in the context of generative modelling with Sinkhorn divergences, in order to derive statistical guarantees for smooth generators.

Remark 5.1. Proposition 5.1 and Theorem 5.1 do not hold for

T C,ε instead of S C,ε because v → T C,ε (φ v 1 α, β) is not lower bounded on L 2 V .
We also emphasize that it is preferable to use a Sinkhorn divergence in practice, since it does not suffer from the aforementioned entropic bias. In particular, the experiments from the next section illustrate that debiasing leads to more accurate registrations.

Remark 5.2. Item (ii) in Proposition 5.2 resembles classical sampling complexity bounds of entropic optimal transport such as [Genevay et al., 2019, Theorem 3], [Séjourné et al., 2019, Theorem 7] and [Mena and Niles-Weed, 2019, Corollary 1]. Our result differs critically by handling a supremum over a class of smooth push-forward maps within the expectation, which enables to prove item (ii) in Theorem 5.1.

Implementation

This section addresses the practical aspects of diffeomorphic registration through Sinkhorn divergence. Firstly, we briefly recall how to compute a minimizer of J λ,n for an arbitrary loss Λ. Then, we illustrate the procedure for Sinkhorn divergences on numerical experiments.

Resolution procedure

This subsection introduces the basic knowledge for solving a diffeomorphic registration problem. It is meant to keep the paper as self-contained as possible. Several minimization strategies coexist, corresponding to different parametrizations of the optimization problem 5. We refer to [Younes, 2010, Section 10.6] for a complete overview of the resolution procedures.

Gradient descent over the time-dependent momentum

To practically minimize J λ,n , one must first write the optimal vector fields v in a finite parametric form, and then perform a gradient descent on the coefficients of this decomposition. Recall that Assumption 3.1 implies that V is a RKHS, thereby characterized by a unique matrix-valued symmetric positive kernel function Ker : R d × R d → R d×d . For simplicity, we address the case of the Gaussian kernel defined in 3. Statistically, the bandwidth parameter σ represents the correlation between the morphed points; physically, it quantifies the fluid viscosity. When σ is small, the points have independent trajectories; when it is large, the points move as a whole.

The RKHS viewpoint enables to parametrize the optimal vector fields through a kernel trick. Firstly, note that the minimization of J λ,n can be formulated as an optimal control problem. It amounts to solving

min v∈L 2 V Λ(α n (1), β n ) + λ v 2 L 2 V ; subject to α n (t) = φ v t α n for any t ∈ [0, 1]. (10)
Then, since the constraint involves a finite number n of trajectories, the so-called reduction principle (see [Glaunes, 2005, Theorem 14]) entails that any solution to problem 10, that is any minimizer of J λ,n , can be written as,

v n t (x) = n i=1 Ker (x, z a i (t)) a i (t),
where the momentum a := (a 1 , . . . , a n ) denotes n unspecified time functions of L 2 ([0, 1], R d ), and the control trajectories z a := (z a 1 , . . . , z a n ) are defined by

z a i (t) = x i + t 0 n j=1
Ker(z a i (s), z a j (s))a j (s)ds.

(11)

This enables to recast (10) as minimizing,

E λ,n (a) := Λ 1 n n k=1 δ z a k (1) , β + λ 1 0 n i,j=1 a i (t) • Ker z a i (t), z a j (t) a j (t)dt, (12)
where • denotes the Euclidean inner product. The gradient of E λ,n was originally derived in [START_REF] Glaunes | Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching[END_REF] for the MMD case, and re-expressed in [START_REF] Glaunes | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Younes | Diffeomorphic learning[END_REF] for more general settings. It can be written as ∇E λ,n (a) = 2λa -p a where p a := (p a 1 , . . . , p a n ) denotes n functions of L 2 ([0, 1], R d ) satisfying for any i ∈ {1, . . . , n} and t ∈ [0, 1],

p a i (t) := ∇ z a i (1) Λ 1 n n k=1 δ z a k (1) , β - 1 σ 2 1 t n j=1 Ker z a i (t), z a j (t) a i (t) • p a j (t) + a j (t) • p a i (t) -2λa i (t) • a j (t) z a i (t) -z a j (t) . ( 13
)
In order to practically track all the functions of the continuous time variable, one must discretize the time scale [0, 1] into τ sub-intervals of equal sizes, which recasts a, z a and p a as (τ + 1) × n × d tensors. Then, equations ( 11) and ( 13) are successively solved at each iteration of the gradient descent by solving the associated discrete dynamical systems. By plugging the solutions z a and p a into the formula of ∇E λ,n (a) one can update the variable a with a ← -a -ξ × (2λa -p a ) where ξ denotes the step size. The computational complexity of an iteration is in O(n 2 dτ ). However, the dynamical systems can be parallelized in the number of points and the dimension. At the end of the process, we obtain the following deformation,

φ a,τ t (x) := x + 1 τ t-1 s=0 n j=1 Ker(x, z a j (s))a j (s). ( 14
)
This approach handles any data-fidelity loss Λ as long as it is differentiable with respect to the data points of the discrete distributions. Both Sinkhorn divergences and squares of MMDs satisfy this property.

Geodesic shooting of the initial momentum

A widely used variant of the above approach is the geodesic shooting of the initial momentum which relies on the equations satisfied at the minimum to uniquely constrain the time-dependent solution a(•) by its initial value, allowing for optimizing solely over a(0). More specifically, as demonstrated in [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF], the Hamiltonian viewpoint of the control problem yields the following joint dynamic of the optimal control trajectories and momentum:

z a i (t) = x i + t 0 n j=1
Ker z a i (s), z a j (s) a j (s)ds,

a i (t) = a(0) - 1 2 ∇ z a i (t) t 0   n j=1 a i (s) • Ker z a i (s), z a j (s) a j (s)   ds. ( 15
)
This entails that both the control trajectories and the momentum at any instant t are fully characterized by a(0). Slightly abusing notations we write z a = z a(0) . Additionally, the kinetic energy remains constant along optimal solutions, implying that

1 0 n i,j=1 a i (t) • Ker z a(0) i (t), z a(0) j (t) a j (t)dt = n i,j=1 a i (0) • Ker (x i , x j ) a j (0). (16)
Therefore, ( 16) together with (15) enable to recast the functional (12) to minimize as

E 0 λ,n (a(0)) := Λ 1 n n k=1 δ z a(0) k (1) , β + λ n i,j=1 a i (0) • Ker (x i , x j ) a j (0), (17)
which is a well-defined function of the time-invariant parameter a(0) ∈ R n×d only. After minimizing (17) using a gradient-descent-based method, one can shoot the obtained a(0) along the discretized system (15) to generate the optimal control trajectories z a (•) and time-dependent momentum a(•). Then, the trajectory of any new point x ∈ R d at any time t ∈ [0, 1] can be computed by integrating the flow equation as in ( 14).

Naturally, for a non-convex program such as (5) the quality of the output solution may heavily depend on the chosen resolution procedure. In the coming experiments, we compare the deformations obtained with both solving strategies.

Numerical experiments

We present a series of numerical experiments on synthetic and real 2-D and 3-D shapes. The objective is to illustrate the practical benefits of using a Sinkhorn divergence as the data-fidelity loss. Our Python code1 operates with the GeomLoss package [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF] to compute the losses and their gradients by automatic differentiation, and the KeOps package [START_REF] Charlier | Kernel operations on the gpu, with autodiff, without memory overflows[END_REF] to handle kernel-reduction operations. It is largely inspired by the example codes from these librairies' websites.2

2-D dataset

In [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF], the authors proposed an alternative measure registration framework based on the gradient flow of the data-fidelity loss. It amounts to updating the source distribution α n := n -1 n i=1 δ xi by carrying out a gradient descent on Λ(α n , β n ) with respect to the positions x 1 , . . . , x n . This model-free method enables to faithfully match one distribution to another, even when the supports have irregularities such as holes. In this section, we firstly adapt their experiments, more precisely the ones from the example section of the GeomLoss package website, by using diffeomorphic deformations instead of gradient flows. The objective is matching two blob-like point clouds in dimension 2. We proceed as follows. Firstly, we learn the optimal matching between two samples of size n = 1, 000 using each of the two previously described procedures. Secondly, we display the obtained time interpolation between two new independent samples of size m = 2, 000. In order to benchmark the influence of the data-fidelity loss, we consider a fixed setting where V is defined through a Gaussian kernel with bandwidth σ = 0.175, the regularization has weight λ = 10 -8 , and the time scale is uniformly divided into τ = 16 intervals. Then, we compare the results for different losses: (unbiased) Sinkhorn divergences, biased entropic transportation costs, and squared Gaussian maximum mean discrepancies. Recall that the squared Gaussian MMD with bandwidth parameter θ > 0 is defined as,

MMD 2 θ (µ, ν) := R d ×R d exp - x -y 2 2θ 2 d(µ -ν)(x)d(µ -ν)(y).
The ground cost function for the Sinkhorn divergences is always C(x, y) := x -y 2 throughout the experiments.

Figures 1 to 3 compare the optimal matchings obtained with respectively the gradient descent on the momentum (GDM) and geodesic shooting (GS) for different values of the relevant parameters ε and θ. Note that whatever the minimization strategy, we used a fixed number of iterations with a constant learning rate, and initialized the momentum with the zero tensor. Also, while we programmed a standard gradient descent for GDM, we relied on the PyTorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] in-built L-BFGS solver for the geodesic shooting. The results are arranged as follows: Figure 1 shows the deformations for both Sinkhorn divergences and (biased) entropic transportation costs optimized with GDM; Figure 2 is the counterpart of Figure 1 for GS; Figure 3 displays the deformations generated by Gaussian maximum mean discrepancies for both resolution procedures. Firstly, we observe from Figures 1 and 2 that entropic optimal-transport metrics yield consistent results across minimization strategies. In contrast, the registration for maximum mean discrepancies depicted in Figure 3 varies with the chosen methods. This instability of the optimization problem underlines that MMDs give more local minima.

Secondly, Figures 1 and2 clearly exhibit the entropic bias: in contrast to Sinkhorn divergences, standard entropic transportation costs shrink the morphed distribution for large values of the regularization parameter ε, leading to unacceptable registrations. However, choosing a too large ε for the unbiased divergence yields a blurry, poorly accurate solution. As expected, debiasing becomes less critical as the regularization diminishes, and both entropic losses provide sharp matchings for small values of ε. Note also that there is no need to decrease ε below a certain threshold to ensure accurate deformations.

Finally, Figure 3 indicates that the consistency of the results between resolution procedures weakens as the bandwidth of the Gaussian kernel decreases. This is due to Gaussian maximum mean discrepancies ignoring disparities smaller than θ. As such, setting a large bandwidth facilitates the registration but degrades the quality of the matching. In contrast, a small bandwidth allows for sharper registration but induces more local minima. This aspect is epitomized for θ = 0.1 in the experiments: with the gradient descent on the time-dependent momentum, the morphed points end up diverging, trapped into minimizing the auto-correlation contribution of the MMD, while geodesic shooting produces a fine matching.

All in all, our experimental observations about the role of the losses are similar to the ones made by [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF] in the context of gradient flows. Critically, compared to their approach, we work with a transformation that is smooth at any time. This regularity constraint reduces the flexibility of the matching, which leads to a less accurate fitting than gradient flows. This affects particularly the anomalous parts of the targeted support, namely the holes and the tail. In contrast, regularity enables the deformation to generalize to any new out-of-sample observations. Additionally, it prevents from tearing the mass apart. The color map on the distribution α m (t) enables to track the location of the moved points through time. Notice that, as a direct consequence of the smoothness, the chromatic continuity between morphed points is preserved throughout the process.

Before turning to more complex 3-D shapes, let us push further the quality analysis of local minima on this illustrative dataset. In the sequel, we consider the same setting as before, and focus on the optimal matchings obtained by geodesic shooting for Sinkhorn divergences and Gaussian maximum mean discrepancies with different parameter values. However, instead of initializing the optimized variable a(0) to zero, we now study the stability and accuracy of the solutions over various initial values. More specifically, we rely on a warm-start strategy: solutions from the above experiments are reused as starting points in the solver. The results are gathered in Figure 4, which reports the final matchings obtained with different initializations along with their associated loss values.

Let us firstly analyze the results for the losses that previously gave the finest registrations: the Sinkhorn divergence with ε = 10 -4 (rows 1 and 5) and the Gaussian MMD with θ = 0.1 (rows 3 and 7). As anticipated, the matchings vary with the initialization. Visually, this phenomenon is stronger for the MMD than for the Sinkhorn divergence and the quality of the final matchings remains quite accurate for the optimal-transport loss. By checking the loss values, we note that the warm start downgrades the solutions for both losses, except for the MMD using initialization via S C,ε with ε = 10 -4 which gets significantly closer to the global minimum. In sum, it seems that the entropic divergence induces fewer or better local minima. Regarding the Sinkhorn divergence with ε = 1 (rows 2 and 6) and and the MMD with θ = 0.5 (rows 4 and 8), which previously yielded imprecise matchings, they have analogous behaviours with respect to warm start. We observe that the results are less robust to initialization and can be significantly improved by using already accurate solutions as starting points, underlining that the registrations obtained with the initialization to zero corresponded to bad local minima.

3-D surfaces

In a second time, we implement the diffeomorphic matching of two shapes embedded in R3 : the source is the unit sphere while the target is the centered scaled Stanford bunny, 3 both encoded through the associated uniform distributions. Similarly to the above experiments, we firstly learn the diffeomorphism on a training set of size n = 5, 000 using geodesic shooting for various losses, and then display the final matching on a testing set of size m = 10, 000. The setup is characterized by σ = 0.05, λ = 10 -8 , and τ = 16. The results can be found in Figure 1. We make comparable observations to before. Powering diffeomorphic registration with a Sinkhorn divergence instead of the biased regularized cost avoids the shrinkage effect of the entropic bias for large values of ε, and the matchings are accurate for both losses when ε is small. The Gaussian MMD requires a small bandwidth θ to potentially fit the bunny, but the solution falls into a poor local minima where several morphed points are not attracted by the target. Note also that, due to their regularity, the deformations tend to smooth the sharpest edges of the target bunny.

Conclusion

We proposed to use Sinkhorn divergences as the fidelity loss in diffeomorphic registration problems. We derived the statistical theory, and illustrated the efficiency of this method compared to past approaches based on MMDs or biased entropic transportation costs. As such, this paper paves way for accurate and smooth measure registration with certifiable asymptotic guarantees. Moreover, carrying out this work led us to further investigate the dual formulation of entropic optimal transport, complementing recent papers on the subject. A first avenue for extension could be to consider the registration of unbalanced measures using Sinkhorn divergences, which would align with the work of [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF]. A second one could be to derive sharper rates of convergences. Notably, [START_REF] Del Barrio | An improved central limit theorem and fast convergence rates for entropic transportation costs[END_REF] which demonstrates faster convergence rates for the empirical entropic transportation potentials and [START_REF] Chizat | Faster wasserstein distance estimation with the sinkhorn divergence[END_REF] which shows that debiasing decreases the approximation error of optimal transport induced by entropic regularization could serve as inspirations.

A Preliminary results

This section recalls some useful results. Section A.1 contains a brief reminder on entropy numbers of classes of functions, in order to derive an upper bound on empirical processes; Section A.2 focuses on the chain rule for composite Frechet derivatives up to arbitrary high orders.

A.1 Empirical processes

In the proof of Proposition 5.2, we will bound the sampling error between the empirical entropic transportation cost and its population counterpart by a centered empirical process indexed by a class of smooth functions.

Recalling the theory introduced in [Van Der [START_REF] Van Der | Weak convergence[END_REF]Wellner, 1996, Koltchinskii, 2011], we present in this subsection intermediary results on such processes. Let X be a compact convex subset of R d . For any probability measure µ on X and r ≥ 1, we define the L r (µ)-norm on C(X , R) as h r,µ := |h| r dµ 1/r . In empirical process theory, the complexity of classes of functions is commonly evaluated through the so-called covering and bracketing numbers. Let H be a class of function included in C(X , R), and > 0 a constant. The covering number N ( , H, L r (µ)) is defined as the minimal number of L r (µ)-balls of radius needed to cover the class of functions H. The center of the balls need not belong to H, but must have finite norm. Additionally, given two functions l and u with finite norm but not necessarily in H, the bracket [l, u] is the set of all functions h such that

l ≤ h ≤ u. An ( , L r (µ))-bracket is a bracket [l, u] such that l -u r,µ ≤ . Then, the bracketing number N [ ] ( , H, L r (µ))
is the minimal number of ( , L r (µ))-bracket needed to cover H. These numbers have essential applications in statistics. The supremum of a centered empirical process indexed by a class of functions with a finite bracketing number converges uniformly almost-surely to zero. Moreover, with a sharper control on the bracketing number, one can derive the following convergence rate: Proposition A.1. Let µ n be an empirical measure of a probability measure µ corresponding to a compact convex subset X of R d , and set H > 0 a constant. Consider the class of functions H := C κ H (X , R) for some integer κ ≥ 0. If κ > d/2, then there exists a constant A = A((H, κ); (X , d)) such that,

E sup h∈H |µ n (h) -µ(h)| ≤ A √ n .
Proof. Combining [Koltchinskii, 2011, Theorem 2.1] with [Koltchinskii, 2011, Theorem 3.11], we directly have that,

E sup h∈H |µ n (h) -µ(h)| ≤ 2 × c √ n E 2σn 0 log N ( , H, L 2 (µ n ))d ,
where c > 0 is some constant and σ n := sup h∈H µ n (h 2 ). By definition of H, it follows that σ n ≤ H 2 . Besides, we can upper bound the covering number in the right term by the bracketing number [Van Der Vaart and Wellner, 1996, page 84]). In addition, according to [Van Der Vaart and Wellner, 1996, Corollary 2.7.2], there exists a constant ρ = ρ((H, κ);

N [ ] (2 , H, L 2 (µ n )) (see
(X , d)) > 0 such that, log N [ ] (2 , H, L 2 (µ n )) ≤ ρ(2 ) -d/κ .
Note that the right term does not depend on µ n . All in all,

E sup h∈H |µ n (h) -µ(h)| ≤ 2c √ n E 2H 2 0 ρ(2 ) -d/κ d .
The integral is finite as κ > d/2. Consequently, the upper bound defines a constant A = A((H, κ); (X , d)). This concludes the proof.

Remark that the convexity assumption on the compact domain X is not restrictive, as it suffices to extend the probability measure µ on the convex hull of X .

A.2 Frechet derivative

The proof of Proposition 4.1 requires bounding the Frechet derivatives of arbitrary high orders of composite functions. We rely on the generalization of Faà di Bruno's formula proposed by [START_REF] Clark | Faa di bruno's formula for variational calculus[END_REF] to carry out the computation.

Let F : R d2 → R d3 and G : R d1 → R d2 be two differentiable functions up to order k ≥ 1. Denote by Ω(k) the set of partitions of {1, . . . , k}, and write |•| for the cardinality of a set. For any δ :

= (δ 1 , . . . , δ k ) ∈ (R d1 ) k , x ∈ R d1 , and ω := {ω 1 , . . . , ω |ω| } ∈ Ω(k), we define δ G ωi (x) := G (|ωi|) [(δ j ) j∈ωi ] for every 1 ≤ i ≤ |ω|.
Then, according to [Clark and Houssineau, 2013, Theorem 2],

(F • G) (k) (x)[δ 1 , . . . , δ k ] = ω∈Ω(k) F (|ω|) (G(x)) δ G ω1 (x), . . . , δ G ω |ω| (x) . (18)
This results implies a chain rule on the operator norms of derivatives of composite functions, which will greatly simplify the computations of later proofs.

Proposition A.2. Let F : R d2 → R d3 and G : R d1 → R d2 be two differentiable functions up to order k ≥ 1.

Then, for any x ∈ R d1 ,

F • G (k) (x) op ≤ ω∈Ω(k) F (|ω|) G(x) op × 1≤i≤|ω| G (|ωi|) (x) op .
Proof. According to the triangle inequality and ( 18)

(F • G) (k) (x) op ≤ ω∈Ω(k) sup δ1 ,..., δ k ≤1 F (|ω|) (G(x)) δ G ω1 (x), . . . , δ G ω |ω| (x) .
Then, we can bound the right term of this inequality by,

ω∈Ω(k) F (|ω|) (G(x)) op × 1≤i≤|ω| sup δ1 ,..., δ k ≤1 δ G ωi (x) op .
In addition, note that for any

x ∈ R d1 , sup δ1 ,..., δ k ≤1 δ G ωi (x) ≤ G (|ωi|) (x) op . Therefore, F • G (k) (x) op ≤ ω∈Ω(k) F (|ω|) G(x) op × 1≤i≤|ω| G (|ωi|) (x) op .

B Proofs of the main results

This sections details all the mathematical proofs of the paper.

Proof of Lemma 3.2. Let us start with a preliminary remark. For any v ∈ L 2 V , it follows from Assumption 3.1 that

1 0 v t p,∞ dt ≤ c V 1 0 v t V dt. Besides, by Cauchy-Schwarz inequality 1 0 v t V dt ≤ v L 2 V , leading to 1 0 v t p,∞ dt ≤ c V v L 2 V .
We now turn to the proof. Recall that by definition φ v t (x) = x + t 0 v s • φ v s (x)ds. Consequently, by the triangle inequality we have for any compact set

K ⊂ R d that sup t∈[0,1],x∈K φ v t (x) ≤ sup x∈K x + 1 0 v s ∞ ds ≤ sup x∈K x + c V v L 2 V .
Therefore, sup

v∈L 2 V,M ,t∈[0,1],x∈K φ v t (x) ≤ sup x∈K x + c V M.
Moreover, combining [Glaunes, 2005, Theorem 5] with the preliminary remark, we know that for any 1 ≤ k ≤ p, there exist two positive constants c k and c k such that for any

v ∈ L 2 V , sup t∈[0,1] (φ v t ) (k) ∞ ≤ c k exp c k v L 2 V .
Hence, sup Proof of Lemma 4.1. Let µ and ν be probability measures on a compact set K ⊂ R d . In a first time, let us show that optimal potentials (f, g) ∈ C(K, R) × C(K, R) for T C,ε (µ, ν) can be chosen as universally-bounded Lipschitz functions. The optimality condition on the potentials (see for instance [START_REF] Genevay | Entropy-Regularized Optimal Transport for Machine Learning[END_REF]) can be written as,

v∈L 2 V,M ,t∈[0,1] (φ v t ) (k) ∞ ≤ c k exp (c k M ) .
exp - f (x) ε = K exp g(y) -C(x, y) ε dν(y).
Remark that since C is continuously differentiable, f is therefore continuously differentiable. Differentiating both sides of this expression leads to,

∇f (x) = K ∇ 1 C(x, y) exp f (x) + g(y) -C(x, y) ε dν(y),
where ∇ 1 denotes the gradient with respect to x, the first variable of C. Let us define Γ f,g C,ε (x, y) := exp f (x)+g(y)-C(x,y) ε . According to the primal-dual relationship [Genevay, 2019, Proposition 7], an optimal solution π to the primal problem has the expression, dπ(x, y) = Γ f,g C,ε (x, y)dµ(x)dν(y).

Since by definition π ∈ Π(µ, ν), we consequently obtain that K Γ f,g C,ε (x, y)dν(y) = 1. Therefore, ∇f ∞ ≤ sup

x,y∈K ∇ 1 C(x, y) .

A similar argument can be made for g. This shows that f and g are -Lipschitz with = ((K, d); C) > 0. Now, note that for any constant c ∈ R, the pair (f + c, g -c) is still a pair of optimal potentials. As a consequence, they can be chosen without loss of generality such that f (x 0 ) = 0 for a given x 0 ∈ K. Thus, using the Lipschitz property we get f (x) ≤ x -x 0 , hence f ∞ ≤ diam(K). To bound g, we use [Genevay et al., 2019, Proposition 1] which states that inf x∈K {f (x) -C(x, y)} ≤ g(y) ≤ sup x∈K {f (x) -C(x, y)}. This entails that g ∞ ≤ f ∞ + sup x,y∈K |C(x, y)| ≤ diam(K) + sup x,y∈K |C(x, y)|. All in all, there exists a constant 1 = 1 ((K, d); C) such that f and g are 1 -bounded and 1 -Lipschitz continuous.

Analogously, one can bound the successive derivatives of f and g up to order q, the maximum order or differentiability of C, using [Genevay et al., 2019, Proposition 1]. In particular, this result ensures that for any

1 ≤ k ≤ q, both f (k)
∞ and g (k) ∞ are bounded by a polynomial in ε -1 whose coefficients depend only on C and K. This implies that there exists a constant m = m((K, d); (C, q); ε) > 0 such that f and g belong to C q m (K, R).

Proof of Proposition 4.2. Let m > 0 and R > 0. Set f, g ∈ C q m (B R , R). Note that the function h f,g C,ε belongs to C q (B R × B R , R). In a first time, we do not focus on any data processing operations, and show that h f,g C,ε and its derivatives up to order q are uniformly bounded. By definition,

h f,g C,ε (x, y) = f (x) + g(y) -ε exp f (x) + g(y) -C x, y ε + ε.
Before going further, we define the constant

C ∞ (R) := max 0≤k≤q sup (x,y)∈B R ×B R C (k) (x, y) op (19)
Then, using the triangle inequality and the bounds on f, g and C we obtain,

h f,g C,ε ∞ ≤ f ∞ + g ∞ + ε exp f ∞ + g ∞ + C ∞ ε + ε ≤ 2m + ε exp 2m + C ∞ ε + ε.
Notice that the upper bound does not depend on the choice of f and g. We prove similar bounds for arbitrary high orders of derivatives using the chain rule. We divide the problem by studying the function,

Γ f,g C,ε : (x, y) ∈ B R × B R → exp f (x) + g(y) -C(x, y) ε ,
which is κ-continuously differentiable. Using Proposition A.2 with F = exp, we obtain for any

1 ≤ k ≤ q, Γ f,g C,ε (k) (x, y) op ≤ Γ f,g C,ε (x, y) ω∈Ω(k) 1≤i≤|ω| ε -1 f (|ωi|) (x) + g (|ωi|) (y) -C (|ωi|) (x, y) op , Then, Γ f,g C,ε (k)
∞ ≤ exp 2m + C ∞ (R) ε ω∈Ω(k) ε -|ω| (2m + C ∞ (R)) |ω| . ( 20
)
We now turn back to h f,g C,ε . Since h f,g C,ε

(k) = f (k) + g (k) -ε Γ f,g C,ε (k)
we finally have

h f,g C,ε (k)
∞ ≤ 2m + exp 2m + C ∞ (R) ε ω∈Ω(k) ε -|ω|+1 (2m + C ∞ (R)) |ω| .
By defining,

H 0 (m; R; (C, q); ε) := (2m + ε) + ε exp 2m + C ∞ (R) ε × max 0≤k≤q    ω∈Ω(k) ε -|ω| (2m + C ∞ (R)) |ω|    , we conclude that h f,g C,ε (k)
q,∞ ≤ H 0 .

We now include data processing transformations. Set

T 1 , T 2 ∈ C p R (X , R d ). It follows from the regularity of h f,g C,ε that h f,g C,ε • (T 1 , T 2 ) ∈ C κ (X × X , R). Since T 1 (x), T 2 (y) ∈ B R , and because h f,g C,ε is bounded by H 0 on B R × B R , the function h f,g C,ε • (T 1 , T 2
) is bounded on X × X regardless of the choice of f, g, T 1 and T 2 . Here again, we use the chain rule to build higher-order bounds. From Proposition A.2 applied with

F = h f,g C,ε and G = (T 1 , T 2 ) it follows that for any 1 ≤ k ≤ κ, h f,g C,ε • (T 1 , T 2 ) (k) (x, y) op ≤ ω∈Ω(k) h f,g C,ε (|ω|) • (T 1 , T 2 )(x, y) op × 1≤i≤|ω| (T 1 , T 2 ) (|ωi|) (x, y) op . (21)
Then, remark that for any 1 ≤ k ≤ κ,

(T 1 , T 2 ) (k) (x, y) 2 op = sup δi ≤1 (T 1 , T 2 ) (k) (x, y)(δ 1 , . . . , δ k )) 2 ≤ sup δi ≤1 T (k) 1 (x)(δ 1 , . . . , δ k ) 2 + sup δi ≤1 T (k) 2 (y)(δ 1 , . . . , δ k ) 2 = T (k) 1 (x) 2 op + T (k) 2 (y) 2 op ≤ 2R 2 .
We can therefore bound the right term of (21), leading to

h f,g C,ε • (T 1 , T 2 ) (k) ∞ ≤ ω∈Ω(k) H 0 × 1≤i≤|ω| √ 2R = H 0 ω∈Ω(k) ( √ 2R) |ω| .
We conclude by defining

H(m; R; (C, q); ε, p) := H 0 (m; R; (C, q); ε) × max 0≤k≤κ    ω∈Ω(k) ( √ 2R) |ω|    , which leads to, h f,g C,ε • (T 1 , T 2 ) κ,∞ ≤ H.
Proof of Proposition 5.1. Let {v n } n∈N be a sequence of vector fields in L 2 V weakly converging to some v ∈ L 2 V . [Glaunes, 2005, Proposition 4] implies that for every x ∈ X ,

φ v n 1 (x) -φ v 1 (x) -----→ n→+∞ 0. (22)
Next, we aim at showing that this entails φ

v n 1 α w -----→ n→+∞ φ v 1 α
, where w denotes the weak* convergence of probability measures. Firstly, note that as a consequence of the uniform-boundedness principle [Rudin, 1991, Theorem 2.5], the weak convergence of {v n } n∈N to v implies that there exists M > 0 such that {v n } n∈N ∪ {v} ⊂ L 2 V,M . Hence, according Lemma 3.1, there exists some R = R((X , d); (V, p); M ) > 0 such that the measures {φ v n 1 α} n∈N , φ v 1 α, and β are all probabilities on B R . Secondly, recall that showing the weak* convergence of amounts to check that for any bounded test functions h ∈ C(B R , R) we have that hd(φ v n 1 α) -----→ n→+∞ hd(φ v 1 α). Let h ∈ C(B R , R) be a bounded function and use the push-forward change-of-variable formula to write hd(φ v n 1 α) = (h • φ v n 1 )dα. By continuity of h and according to ( 22), the sequence of functions {h • φ v n 1 } n∈N converges point-wise to h • φ v 1 . In addition, as h is bounded, this sequence is dominated by a constant. We can therefore apply the dominated convergence theorem to obtain that φ v n 1 α w -----→ n→+∞ φ v 1 α. We conclude the proof using [Feydy et al., 2019, Proposition 13], which states that T C,ε (and consequently S C,ε ) is weak* continuous w.r.t. each of its input measures, provided that the ground cost function C is Lipschitz on their compact domains. This condition readily follows from the continuity of the derivative of C on the

compact set B R × B R . Therefore, v → S C,ε (φ v 1 α, β) is weakly continuous on L 2 V . If additionally e -C ε defines a positive universal kernel, then v → S C,ε (φ v 1 α, β
) is non negative according to [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF] Theorem 1], which implies through [Glaunes, 2005, Theorem 7] that J λ for Λ = S C,ε admits minimizers.

Proof of Proposition 5.3. Let R > 0. In a first time, we demonstrate the following Glivenko-Cantelli theorem

(i): sup T1,T2∈C p R (X ,R d ) T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β) a.s.
-----→ n→+∞ 0.

In a second time, when κ = min{p, q} ≥ d, we show the following rate of convergence (ii):

E sup T1,T2∈C p R (X ,R d ) T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β) ≤ A √ n ,
where A > 0 is a constant. In both cases, the key idea of the proof is to note that the quantity sup

T1,T2∈C p R (X ,R d ) T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β)
is the supremum of a centered empirical process indexed by a class of smooth functions, and as such can be controlled via classical results from empirical process theory (see Section A.1).

Let T 1 and T 2 be two arbitrary functions in C p R (X , R d ). By definition, the image sets T 1 (X ) and T 2 (X ) are contained in B R . Thus, using the dual formulation, the entropic transportation costs can be written as,

T C,ε (T 1 α n , T 2 β n ) = sup f,g∈C(B R ,R) (T 1 α n ⊗ T 2 β n )(h f,g C,ε ), T C,ε (T 1 α, T 2 β) = sup f,g∈C(B R ,R) (T 1 α ⊗ T 2 β)(h f,g C,ε ).
We apply Lemma 4.1 with µ = T 1 α and ν = T 2 β which are probability measures on B R . This implies that there exists a constant m = m(B R ; (C, q), ε) > 0 such that

T C,ε (T 1 α n , T 2 β n ) = sup f,g∈C q m (B R ,R) (T 1 α n ⊗ T 2 β n )(h f,g C,ε ) = sup f,g∈C q m (B R ,R) (α n ⊗ β n )(h f,g C,ε • (T 1 , T 2 )),
where we used the push-forward change-of-variable formula. Proceeding similarly with the empirical measures we get,

T C,ε (T 1 α, T 2 β) = sup f,g∈C q m (B R ,R) (T 1 α ⊗ T 2 β)(h f,g C,ε ) = sup f,g∈C q m (B R ,R) (α ⊗ β)(h f,g C,ε • (T 1 , T 2 )),
Then, by using a classical error decomposition, we can control the difference between these two terms as follows,

|T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β)| ≤ sup f,g∈C q m (B R ,R) |(α n ⊗ β n )(h f,g C,ε • (T 1 , T 2 )) -(α ⊗ β)(h f,g C,ε • (T 1 , T 2 ))|.
After taking the supremum in C p R (X , R d ) on both sides of this inequality we get,

sup T1,T2∈C p R (X ,R d ) T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β) ≤ sup T1,T2∈C p R (X ,R d );f,g∈C q m (B R ,R) (α n ⊗ β n )(h f,g C,ε • (T 1 , T 2 )) -(α ⊗ β)(h f,g C,ε • (T 1 , T 2 )) .
The right term of this inequality can be seen as a centered empirical process indexed by the class of functions

{h f,g C,ε • (T 1 , T 2 ) | T 1 , T 2 ∈ C p R (X , R d ); f, g ∈ C q m (B R , R)}.
Empirical process theory provides convergence guarantees when the index class is regular enough. Besides, we know from Proposition 4.1 that there exists a constant H := H(R; (C, q); ε, p) > 0 such that this class is included in C κ H (X × X , R). Therefore, Let us set H := C κ H (X × X , R). According to [Van Der Vaart and Wellner, 1996, Corollary 2.7.2] and [Van Der Vaart and Wellner, 1996, Theorem 2.4.1], H is a so-called (α ⊗ β)-Glivenko-Cantelli class of functions, meaning that

sup h∈H |(α n ⊗ β n )(h) -(α ⊗ β)(h)| a.s.
-----→ n→+∞ 0.

This implies (i). In addition, by Proposition A.1, if κ ≥ (2d)/2 then there exists a positive constant A := A(R; (C, q); ε; (X , d); p) such that,

E sup h∈H |(α n ⊗ β n )(h) -(α ⊗ β)(h)| ≤ A √ n .
This proves (ii).

v∈L 2 V J λ (v) ⊆ L 2 V,M .
Proof. The proof generalizes an argument made for a squared MMD in [Glaunes, 2005, Theorem 16] to a Sinkhorn divergence. Let n ∈ N and set v n a minimizer of J λ,n . Notice that the vector flow uniformly equal to zero generates the identity function, that is φ 0 t = I for any t ∈ [0, 1]. Thus, by definition of a minimizer and by non negativity of the Sinkhorn divergence, we readily have that

λ v n 2 L 2 V ≤ J λ,n (v n ) ≤ J λ,n (0) = S C,ε (α n , β n ).
Therefore, v n 2 L 2 V ≤ λ -1 S C,ε (α n , β n ). To conclude, let us bound uniformly the right-term of this inequality. According to Lemma 4.1 applied with α n and β n there exists a constant m = m((X , d); (C, q); ε) such that,

T C,ε (α n , β n ) = sup f,g∈C q m (X ,R) (α n ⊗ β n ) h f,g C,ε .
Moreover, for any x, y ∈ X , h f,g C,ε (x, y) = f (x) + g(x) -εe f (x)+g(y)-C(x,y) ε

+ ε ≤ m + m + 0 + ε.
Thus, T C,ε (α n , β n ) ≤ 2m + ε.

The same bound holds for the two auto-correlation terms of the Sinkhorn divergence, namely T C,ε (α n , α n ) and T C,ε (β n , β n ). Therefore, the triangle inequality leads to

S C,ε (α n , β n ) ≤ 4m + 2ε.
Consequently,

v n 2 L 2 V ≤ 4m + 2ε λ .
To conclude, we set M (λ; (X , d); (C, q); ε) := 4m+2ε λ . Note that this bound does not depend on n. As such, the minima {v n } n∈N all belong to L 2 V,M . A similar reasoning for v * a minimizer of J λ shows that all the minimizers of J λ also belong to L 2 V,M .

Proof of Theorem 5.2. Let M > 0 be arbitrary (for now). Set v ∈ L 2 V,M and compute

|J λ,n (v) -J λ (v)| = |S C,ε (φ v 1 α n , β n ) -S C,ε (φ v 1 α, β)|.
According to Lemma 3.1, there exists a constant R = R((X , d); (V, p); M ) such that for any φ ∈ {φ v t | t ∈ [0, 1], v ∈ L 2 V,M }, the restriction φ| X and the identity function I both belong to C p R (X , R d ). This leads to

sup v∈L 2 V,M |J λ,n (v) -J λ (v)| ≤ sup T1,T2∈C p R (X ,R d ) |S C,ε (T 1 α n , T 2 β n ) -S C,ε (T 1 α, T 2 β)|. ( 23
)
From here, let us demonstrate the convergence of the minima, that is item (i). According to Lemma B.1, there exists M = M (λ; (X , d); (C, q); ε) > 0 such that all the minimizers of J λ,n belong to L 2 V,M . Next, we show that any weakly-converging subsequences of {v n } n∈N tend to a minimizer of J λ . Set v * a minimizer of J λ , and let {u n } n∈N be a subsequence with limit u. First, let's show that lim n→+∞ J λ,n (u n ) = J λ (u). By the triangle inequality, |J λ,n (u n ) -J λ (u)| ≤ |J λ,n (u n ) -J λ (u n )| + |J λ (u n ) -J λ (u)|. The first term tends to zero by Proposition 5.2 and ( 23) specified with M (λ; (X , d); (C, q); ε), while the second term tends to zero according to Proposition 5.1 which ensures the weak continuity of J λ . Second, note that the optimality condition entails that J λ,n (u n ) ≤ J λ,n (v * ), and that lim n→+∞ J λ,n (v * ) = J λ (v * ). Then, at the limit J λ (u) ≤ J λ (v * ), meaning that u is a minimizer of J λ . Therefore, any weakly-converging subsequence {u n } n∈N of {v n } n∈N tends to a minimizer u of J λ .

To conclude on the convergence of the generated diffeomorphisms, we rely on [START_REF] Glaunes | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF]

φ u n t -φ u t ∞ + (φ u n t ) -1 -(φ u t ) -1 ∞ ≤ 2c V u n -u L 2 V exp c V u L 2 V .
We showed that u n -u L 2 V ----→ n→∞ 0. Consequently, the upper bound tends to zero as n increases to infinity.

This completes the proof of (i).

Item (ii) readily follows from Proposition 5.2 stating that if κ ≥ d, then there exists for any M > 0 a constant A = A(λ; (X , d); (C, q); ε; (V, p); M ) > 0 such that

E sup v∈L 2 V,M |J λ,n (v) -J λ (v)| ≤ A √ n .
To conclude, recall that both {v n } n∈N and v * belong to L 2 V,M for the constant M from Lemma B.1, and apply the classical deviation inequality

J λ (v n ) -J λ (v * ) ≤ 2 sup v∈L 2 V,M |J λ,n (v) -J λ (v)|.

  Figure 1: 3-D diffeomorphic matchings. Both shapes (a) and (b) are centered scaled.

  d); (V, p); M ) := max max 1≤k≤p {c k exp(c k M )}, sup x∈K x + c V M concludes the proof.

  ,R d ) |T C,ε (T 1 α n , T 2 β n ) -T C,ε (T 1 α, T 2 β)| ≤ sup h∈C κ H (X ×X ,R) |(α n ⊗ β n )(h) -(α ⊗ β)(h)|.

  Before proving Theorem 5.1, we need the next intermediary result: Lemma B.1. Under the assumptions of Theorem 5.1, there exists a positive constant M = M (λ; (X , d); (C, q); ε) such that n∈N arg min v∈L 2 V J λ,n (v) ∪ arg min

https://github.com/lucasdelara/lddmm-sinkhorn/

https://www.kernel-operations.io/geomloss/ and https://www.kernel-operations.io/keops/

http://graphics.stanford.edu/data/3Dscanrep/

Table 2: Optimal-transport-driven 2-D diffeomorphic registration optimized by GS. The colored distribution is α m (t), while the blue distribution is β m .

GDM (+19, 8%)

Table 4: 2-D diffeomorphic matching optimized by GS with warm start. A row specifies the studied loss while each column refers to an initialization. Except for the first column which indicates the initialization to zero, the columns correspond to a solution from the previous experiments. The first four rows show the final matchings while the last four rows give the associated loss values with their growth rates compared to the initialization via zero; row-wise-minimal loss values are written in bold.