Jun Kitagawa

Quentin Mérigot

Boris Thibert

AND Quentin M Érigot

Convergence of a Newton algorithm for semi-discrete optimal transport

published or not. The documents may come

Introduction

Many problems in geometric optics or convex geometry can be recast as optimal transport problems: this includes the far-field reflector problem, Alexandrov's curvature prescription problem, etc. A popular way to solve these problems numerically is to assume that the source probability measure is absolutely continuous while the target measure is finitely supported. We refer to this setting as semi-discrete optimal transport. Among the several algorithms proposed to solve semi-discrete optimal transport problems, one currently needs to choose between algorithms that are slow but come with a convergence speed analysis [START_REF] Oliker | On the numerical solution of the equation and its discretizations, i[END_REF][START_REF] Luis A Caffarelli | Problem of reflector design with given far-field scattering data, Monge Ampère Equation: Applications to Geometry and Optimization[END_REF][START_REF] Kitagawa | An iterative scheme for solving the optimal transportation problem[END_REF] or algorithms that are much faster in practice but which come with no convergence guarantees [START_REF] Aurenhammer | Minkowski-type theorems and leastsquares clustering[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Fernando De Goes | Blue noise through optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3d[END_REF][START_REF] Machado | Intersection of paraboloids and application to minkowski-type problems[END_REF]. Algorithms of the first kind rely on coordinate-wise increments and the number of iterations required to reach the solution up to an error of ε is of order N 3 /ε, where N is the number of Dirac masses in the target measure. On the other hand, algorithms of the second kind typically rely on the formulation of the semi-discrete optimal transport problem as an unconstrained convex optimization problem which is solved using a Newton or quasi-Newton method.

The purpose of this article is to bridge this gap between theory and practice by introducing a damped Newton's algorithm which is experimentally efficient and by proving the global convergence of this algorithm with optimal rates. The main assumptions is that the cost function satisfies a condition that appears in the regularity theory for optimal transport (the Ma-Trudinger-Wang condition) and that the support of the source density is connected in a quantitative way (it must satisfy a weighted Poincaré-Wirtinger inequality). In §1.7, we compare this algorithm and the convergence theorem to previous computational approaches to optimal transport. 1.1. Semi-discrete optimal transport. The source space is an open domain Ω of a d-dimensional Riemannian manifold, which we endow with the measure H d g induced by the Riemannian metric g on the manifold. The target space is an (abstract) finite set Y . We are given a cost function c on the product space Ω × Y , or equivalently a collection (c y ) y∈Y of functions on Ω. We assume that the functions (c y ) are of class C 1,1 on Ω:

∀y ∈ Y, c y ∈ C 1,1 (Ω). (Reg)
We consider a compact subset X of Ω and ρ a probability density on X, i.e. such that ρdH d is a probability measure. By an abuse of notation, we will often conflate the density ρ with the measure ρdH d itself. Note that the support of ρ is contained in X, but we do not assume that it is actually equal to X. The push-forward of ρ by a measurable map T : X → Y is the finitely supported measure T # ρ = y∈Y ρ(T -1 (y))δ y . The map T is called a transport map between ρ and a probability measure µ on Y if T # ρ = µ. The semi-discrete optimal transport problem consists in minimizing the transport cost over all transport maps between ρ and µ, that is min

X c(x, T (x))ρ(x)dH d g (x) | T : X → Y s.t. T # ρ = µ . (M)
This problem is an instance of Monge's optimal transport problem, where the target measure is finitely supported. Kantorovich proposed a relaxed version of the problem (M) as an infinite dimensional linear programming problem over the space of probability measures with marginals ρ and µ.

1.2. Laguerre tessellation and economic interpretation. In the semidiscrete setting, the dual of Kantorovich's relaxation can be conveniently phrased using the notion of Laguerre tessellation. We start with an economic metaphor. Assume that the probability density ρ describes the population distribution over a large city X, and that the finite set Y describes the location of bakeries in the city. Customers leaving at a location x in X try to minimize the walking cost c(x, y), resulting in a decomposition of the space called a Voronoi tessellation. The number of customers received by a bakery is equal to the integral of ρ over its Voronoi cell,

Vor y := {x ∈ Ω | ∀z ∈ Y, c(x, y) c(x, z)}.

If the price of bread is given by a function ψ : Y → R, customers leaving at location x in X make a compromise between walking cost and price by

When the sets X and Y are contained in R d and the cost is the squared Euclidean distance, the computation of the Laguerre tessellation is a classical problem of computational geometry, for which there exists very efficient software. The shape of the Voronoi and Laguerre tessellations is depicted in Figure 1.

We want the Laguerre cells to form a partition of Ω up to a negligible set. By the implicit function theorem, this will be the case if the following twist condition holds, ∀x ∈ X, y ∈ Y → Dc y (x) ∈ T *

x Ω is injective, (Twist) where D is the differential of the function c y . The twist condition implies that for any prices ψ on Y , the transport map induced by the Laguerre tessellation T ψ (x) := arg min y∈Y (c(x, y) + ψ(y)),

is uniquely defined almost everywhere. It is easy to see (see Proposition 2.2), that T ψ is an optimal transport map between ρ and the pushforward measure T ψ# ρ = y∈Y ρ(Lag y (ψ))δ y .

1.3. Kantorovich's functional. Conversely, Theorem 1.1 below ensures that any semi-discrete optimal transport problem admits such a solution. In other words, for any probability density ρ on X and any probability measures µ on Y there exists a function (price) ψ on Y such that T ψ# ρ = ν. The proof of this theorem is an easy generalization of the proof given in [START_REF] Aurenhammer | Minkowski-type theorems and leastsquares clustering[END_REF] for the quadratic cost, but it is nonetheless included in Section 2 for the sake of completeness.

Here and after, we denote (1 y ) y∈Y the canonical basis of R Y , and • the Euclidean norm induced by this basis, while • g will denote the norm induced by the Riemannian metric g on either T x Ω or T *

x Ω (which will be clear from context). We will, slightly abusively, consider the space of probability measures P(Y ) as a subset of R Y .

Theorem 1.1. Assume (Reg) and (Twist), let ρ be a bounded probability density on Ω and ν = y∈Y ν y δ y in P(Y ). Then, the functional Φ Φ(ψ) := is concave, C 1 -smooth, and its gradient is

∇Φ(ψ) = y∈Y (ρ(Lag y (ψ))) -ν y )1 y . (1.4)
Corollary 1.2. The following statements are equivalent: (i) ψ : Y → R is a global maximizer of Φ ;

(ii) T ψ is an optimal transport map between ρ and ν ;

(iii) T ψ# ρ = ν, or equivalently, ∀y ∈ Y, ρ(Lag y (ψ)) = ν y (MA)

We call Kantorovich's functional the function Φ introduced in (1.3). Note that both this functional and its gradient are invariant by addition of a constant. The non-linear equation (MA) can be considered as a discrete version of the generalized Monge-Ampère equation that characterizes the solutions to optimal transport problems (see for instance Chapter 12 in [START_REF] Villani | Optimal transport: old and new[END_REF]).

1.4. Damped Newton algorithm. We consider a simple damped Newton's algorithm to solve semi-discrete optimal transport problem. This algorithm is very close to the one used in [START_REF] Mirebeau | Discretization of the 3d monge-ampere operator, between wide stencils and power diagrams[END_REF] and has been suggested to us by Mirebeau. To phrase this algorithm in a more general way, we introduce a notation for the measure of Laguerre cells: for ψ ∈ R Y we set

G(ψ) := y∈Y G y (ψ)1 y where G y (ψ) = ρ(Lag y (ψ)), (1.5)
so that ∇Φ(ψ) = G(ψ)µ. In the algorithm (Algorithm 1), we denote by A + the pseudo-inverse of the matrix A.

The goal of this article is to prove the global convergence of this damped Newton algorithm and to establish estimates on the speed of convergence. As shown in Proposition 6.1, the convergence of Algorithm 1 depends on the regularity and strong monotonicity of the map G = ∇Φ. As we will see, the regularity of G will depend mostly on the geometry of the cost function and the regularity of the density. On the other hand, the strong monotonicity of G will require a strong connectedness assumption on the support of ρ, in the form of a weighted Poincaré-Wirtinger inequality. Before stating our main theorem we give some indication about these intermediate regularity and monotonicity results and their assumptions.

Input: A tolerance η > 0 and an initial ψ 0 ∈ R Y such that

ε 0 := 1 2 min min y∈Y G y (ψ 0 ), min y∈Y µ y > 0. (1.6) While: G y (ψ k ) -µ y η Step 1: Compute d k = -DG(ψ k ) + (G(ψ k ) -µ)
Step 2: Determine the minimum ∈ N such that

ψ k := ψ k + 2 -d k satisfies      min y∈Y G y (ψ k )) ε 0 G(ψ k ) -µ (1 -2 -( +1) ) G(ψ k ) -µ
Step 3:

Set ψ k+1 = ψ k + 2 -d k and k ← k + 1.
Algorithm 1: Simple damped Newton's algorithm 1.5. Regularity of Kantorovich's functional and MTW condition.

In order to establish the convergence of a damped Newton algorithm for (MA), we need to study the C 2,α regularity of Kantorovich's functional Φ. However, while C 1 regularity of Φ follows rather easily from the (Twist) hypothesis (or even from weaker hypothesis, see Theorem 2.1), higher order regularity seems to depend on the geometry of the cost function in a more subtle manner. We found that a sufficient condition for the regularity of Φ is the Ma-Trudinger-Wang condition [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF], which appeared naturally in the study of the regularity of optimal transport maps. We use a suitable discretization of Loeper's geometric reformulation of the Ma-Trundinger-Wang condition [START_REF] Loeper | On the regularity of solutions of optimal transportation problems[END_REF]. are quasi-convex for all z in Y . The map exp c y is called the c-exponential with respect to y, and the domain Ω y is an exponential chart. Definition 1.2 (c-Convexity). Assuming Loeper's condition, a subset X of Ω is c-convex if for every y in Y , the inverse image (exp c y ) -1 (X) is convex. Note that by assumption, the domain Ω itself is c-convex. The connection between this discrete version of Loeper's condition and the conditions used in the regularity theory for optimal transport is detailed in Remark 1.1. The (QC) condition implies the convexity of Laguerre cells in the exponential charts, which plays a crucial role in the regularity of Kantorovich's functional.

Theorem 1.3. Assume (Reg), (Twist), and (QC). Let X be a compact, c-convex subset of Ω and let ρ in P ac (X) ∩ C α (X) for α in (0, 1]. Then, Kantorovich's functional is of class C 2,α loc on the set

K + := {ψ : Y → R | ∀y ∈ Y, ρ(Lag y (ψ)) > 0}, (1.7)
and its Hessian is given by

(z = y) ∂ 2 Φ ∂1 y ∂1 z (ψ) = Lag y (ψ)∩Lag z (ψ) ρ(x) Dc y (x) -Dc z (x) g dH d-1 g (x), ∂ 2 Φ ∂1 2 y (ψ) = - z∈Y \{y} ∂ 2 Φ ∂1 y ∂1 z . (1.8)
The proof of this theorem and a more precise statement are given in Section 4 (Theorem 4.1), showing that the C 2,α estimate can be made uniform when the mass of the Laguerre cells is bounded from below by a positive constant.

Remark 1.1. We remark that under certain assumptions on the cost c, our (QC) condition is implied by classical conditions introduced in a smooth setting by X.-N. Ma, N. Trudinger, and X.-J. Wang [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF], which include the well known (MTW) or (A3) condition. See Remark 4.3 for more specifics.

There are a wide variety of known examples satisfying these conditions. Aside from the canonical example of the inner product on R n × R n , and other costs on Euclidean spaces mentioned in [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF][START_REF] Trudinger | On the second boundary value problem for Monge-Ampère type equations and optimal transportation[END_REF], there are the nonflat examples of Riemannian distance squared andlog xy R n+1 on (a subset of) S n × S n (see [START_REF]Regularity of optimal maps on the sphere: the quadratic cost and the reflector antenna[END_REF]). The last cost is associated to the far-field reflector antenna problem. We refer the reader to [15, p. 1331] for a (more) comprehensive list of such costs. 1.6. Strong concavity of Kantorovich's functional. As noted earlier, Kantorovich's functional Φ cannot be strictly concave, since it is invariant under addition of a constant. This implies that the Hessian D 2 Φ has a zero eigenvalue corresponding to the constants. A more serious obstruction to the strict concavity of Φ at a point ψ arises when the discrete graph induced by the Hessian (where two points are connected iff ∂ 2 Φ/∂1 y ∂1 z (ψ) = 0) is not connected. This can happen either because one of the Laguerre cells is empty (hence not connected to any neighbor) or if the support of the probability density ρ is itself disconnected. In order to avoid the latter phenomena, we will require that (X, ρ) satisfies a weighted L 1 Poincaré-Wirtinger inequality. Definition 1.3 (weighted Poincaré-Wirtinger). A continuous probability density ρ on a compact set X ⊆ Ω satisfies a weighted Poincaré-Wirtinger inequality with constant

C pw > 0 if for every C 1 function f on X, f -E ρ (f ) L 1 (ρ) C pw ∇f L 1 (ρ) , (PW)
where

h L 1 (ρ) := X |h(x)|ρ(x)dH d g (x) and E ρ (f ) := X f (x)ρ(x)dH d g (x). We denote E Y the orthogonal complement (in R Y ) of the space of constant functions on Y , that is E Y := {ψ ∈ R Y | y ψ(y) = 0}.
As before, K + is the set of functions ψ whose Laguerre cells all have positive mass.

Theorem 1.4. Assume (Reg), (Twist), (QC). Let X be a compact, c-convex subset of Ω, and ρ be a continuous probability density on X satisfying (PW). Then, Kantorovich's functional Φ is strictly concave on

E Y ∩ K + .
As before, a more quantitative statement is proven in Section 5 (Theorem 5.1), establishing strong concavity of Φ assuming that the mass of the Laguerre cells is bounded from below by a positive constant. 1.7. Convergence result. Putting Proposition 6.1, Theorem 1.3 and Theorem 1.4 together, we can prove the global convergence of the damped Newton algorithm for semi-discrete optimal transport (Algorithm 1) together with optimal convergence rates.

Theorem 1.5. Assume (Reg), (Twist) and (QC) and also that (i) The support of the probability density ρ is included in a compact, c-convex subset X of Ω, and ρ ∈ C α (X) for α in (0, 1]. (ii) ρ has positive Poincaré-Wirtinger (PW) constant.

Then the damped Newton algorithm for semi-discrete optimal transport (Algorithm 1) converges globally with linear rate and locally with rate 1 + α.

Comparison to previous work. There exist a few other numerical methods relying on Newton's algorithm for the resolution of the standard Monge-Ampère equation or for the quadratic optimal transport problem. Here, we highlight some of the differences between Algorithm 1 and Theorem 1.5 and these existing results. First, we note that many authors have reported the good behavior in practice of Newton's or quasi-Newton's methods for solving discretized Monge-Ampère equations or optimal transport problems [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Fernando De Goes | Blue noise through optimal transport[END_REF][START_REF] Benamou | Numerical solution of the optimal transportation problem using the monge-ampere equation[END_REF]. Note however that none of these works contain convergence proofs for the Newton's algorithm.

Loeper and Rappetti [START_REF] Loeper | Numerical solution of the monge-ampère equation by a newton's algorithm[END_REF] (refined by Saumier, Agueh, and Khouider [START_REF] Saumier | An efficient numerical algorithm for the l2 optimal transport problem with periodic densities[END_REF]) establish the global convergence of a damped Newton's method for solving quadratic optimal transport on the torus, relying heavily on Caffarelli's regularity theory. In particular, the convergence of the algorithm requires a positive lower bound on the probability densities, while this condition is not necessary for Theorem 1.5 (see Section 5 and Appendix A where we construct explicitly probability densities with non-convex support that still satisfy the hypothesis of Theorem 1.5). A second drawback on relying on the regularity theory for optimal transport is that the damping parameter, which is an input parameter of the algorithm used in [START_REF] Loeper | Numerical solution of the monge-ampère equation by a newton's algorithm[END_REF], cannot be determined explicitly from the data. Third, the convergence proof is for continuous densities, and it seems difficult to adapt it to the space-discretized problem. On the positive side, it seems likely that the convergence proof of [START_REF] Loeper | Numerical solution of the monge-ampère equation by a newton's algorithm[END_REF][26] can be adapted to cost functions satisfying the Ma-Trudinger-Wang condition (which is equivalent to Loeper's condition (QC) that we also require).

Finally Oliker and Prussner prove the local convergence of Newton's method for finding Alexandrov's solutions to the Monge-Ampère equation det D 2 u = ν with Dirichlet boundary conditions, where ν is a finitely supported measure [START_REF] Oliker | On the numerical solution of the equation and its discretizations, i[END_REF]. Global convergence for a damped Newton's algorithm is established by Mirebeau [START_REF] Mirebeau | Discretization of the 3d monge-ampere operator, between wide stencils and power diagrams[END_REF] for a variant of Oliker and Prussner's scheme, but without convergence rates. Theorem 1.5 article can be seen as an extension of the strategy used by Mirebeau to optimal transport problems, which amounts to (a) replacing the Dirichlet boundary conditions with the second boundary value conditions from optimal transport (b) replacing the Lebesgue measure by more general probability densities and (c) changing the Monge-Ampère equation itself in order to deal with more general cost functions.

Outline. In Section 2, we establish the differentiability of Kantorovich's functional Φ, adapting arguments from [START_REF] Aurenhammer | Minkowski-type theorems and leastsquares clustering[END_REF]. In Sections 3 and 4, we prove the (uniform) second-differentiability of Kantorovich's functional when the cost function satisfies Loeper's (QC) condition. Section 5 is devoted to the proof of uniform concavity of Kantorovich's functional, when the probability density satisfies a Poincaré-Wirtinger inequality (PW). In Section 6, we combine these intermediate results to prove the convergence of the damped Newton's algorithm (Theorem 1.5), and we present a numerical illustration. Appendix A presents an explicit construction of a probability density with non-convex support over R d which satisfies the assumptions of Theorem 1.5. Appendix B contains an elementary proof of a convex geometry result used in the regularity proof. Finally, Appendix C contains an alternate proof of a crucial transversality condition, under the assumption that the target set Y is sampled from some continuous domain.

Kantorovich's functional

The purpose of this section is to present the variational formulation introduced in [START_REF] Aurenhammer | Minkowski-type theorems and leastsquares clustering[END_REF] for the semi-discrete optimal transport problem, adapting the arguments presented for the squared Euclidean cost in [START_REF] Aurenhammer | Minkowski-type theorems and leastsquares clustering[END_REF] to cost functions satisfying (Reg') and (Twist'), which are weaker than the conditions (Reg) and (Twist) presented in the introduction:

∀y ∈ Y, c y ∈ C 0 (Ω) (Reg') ∀y = z ∈ Y, ∀t ∈ R, H d g ((c y -c z ) -1 (t)) = 0 ( Twist
') Note that under (Twist'), the map T ϕ : X → Y defined by (1.2) is uniquelydefined H d g -almost everywhere. Most of the results presented here are well known in the optimal transport literature, we include proofs for completeness.

Theorem 2.1. Assume (Reg') and (Twist'), and let ρ be a bounded probability density on Ω and ν = y∈Y ν y δ y be a probability measure over Y . Then, the functional Φ defined by (1.3) is concave, C 1 -smooth, and its gradient is given by (1.4).

The proof of Theorem 2.1 relies on Propositions 2.2 and 2.3. Proposition 2.2. For any ψ : Y → R, the map T ψ is an optimal transport map for the cost c between any probability density ρ on Ω and the pushforward measure ν := T ψ# ρ.

Proof. Assume that ν = S # ρ where S is a measurable map between X and Y . Then, by definition of T ψ one has ∀x ∈ X, c(x, T ψ (x)) + ψ(T ψ (x)) c(x, S(x)) + ψ(S(x)).

Multiplying this inequality by ρ and integrating it over X gives

X (c(x, T ψ (x))+ψ(T ψ (x)))ρ(x)dH d g (x) X (c(x, S(x))+ψ(S(x)))ρ(x)dH d g (x)
Since ν = S # ρ = T ψ# ρ, the change of variable formula gives

X ψ(S(x))ρ(x)dH d g (x) = Y ϕ(y)dν = X ψ(T ψ (x))ρ(x)dH d g (x)
Substracting this equality from the inequality above shows that T ψ is optimal:

X c(x, T ψ (x))ρ(x)dH d g (x) X c(x, S(x))ρ(x)dH d g (x)
Proposition 2.3. Assume (Twist') and (Reg'). Let ρ be a probability density over a compact subset X of Ω. Then, the map G : R Y → R Y is continuous:

G(ψ) = ρ(Lag y (ψ)) y∈Y (2.9)
Lemma 2.4. Let ρ be a probability density over a compact subset X of Ω, and let f in C 0 (X) be such that ρ(f -1 (t)) = 0 for all t ∈ R. Then, the function

g : t → ρ(f -1 ((-∞, t]) is continuous.
Proof. We consider the function

h(t) = ρ(f -1 ((-∞, t)). By hypothesis, g(t) -h(t) = ρ(f -1 (t)) = 0.
H z (ψ) := {x ∈ X | c(x, y) + ψ(y) c(x, z) + ψ(z)}.
Denoting A∆B the symmetric difference of two sets, we have the following inequalities

|G y (ψ) -G y (ϕ)| ρ(Lag y (ψ)∆ Lag y (ϕ)) z =y∈Y ρ(H z (ψ)∆H z (ϕ)).
(2.10) Fix z = y ∈ Y , and denote f = c(•, y)c(•, z). Then,

H z (ψ)∆H z (ϕ) ⊆ f -1 ([ψ(z) -ψ(y), ϕ(z) -ϕ(y)]).
By (Twist') and Lemma 2.4 we know that lim ϕ→ψ ρ(H z (ψ)∆H z (ϕ)) = 0, which with (2.10) concludes the proof.

2.1. Proof of Theorem 1.1. We simultaneously show that the functional is concave and compute its gradient. For any function ψ on Y and any measurable map T : X → Y , one has min y∈Y (c(x, y) + ψ(y)) c(y, T (y)) + ψ(T (y)), which by integration gives Φ(ψ)

X (c(x, T (x)) + ψ(T (x)))ρ(x)dH d g (x) - y∈Y ψ(y)ν y . (2.11)
Moreover, equality holds when T = T ψ . Taking another function ϕ on Y and setting T = T ϕ in Equation (2.11) gives

Φ(ψ) Φ(ϕ) + G(ϕ) -ν | ψ -ϕ ,
where G is defined as in the statement of Proposition 2.3. This proves that the superdifferential of Φ is never empty, thus establishing the concavity. By Proposition 2.3, the map G is continuous, meaning that we have constructed a continuous selection of a supergradient in the superdifferential of the concave function Φ. By standard arguments from convex analysis, this proves that Φ is C 1 , and that ∇Φ(ψ) = G(ψ)ν.

Local regularity in a c-exponential chart

The results presented in this section constitute an intermediate step in the proof of C 2,α regularity of Kantorovich's functional. Let X be a compact, convex subset of R d and f 1 , . . . , f N be C 1,1 functions on X which are quasi-convex, meaning that for any scalar κ ∈ R the closed sublevel sets

K i (κ) := f -1 i ([-∞, κ]) are convex.
Let ρ be a continuous probability density over X. The purpose of this section is to give sufficient conditions to ensure the regularity of the following function Ĝ near the origin of

R N : Ĝ : λ ∈ R N → K(λ) ρ(x)dH d (x), (3.12)
where

K(λ) := X ∩ N i=1 K i (λ i ) = {x ∈ X | ∀i ∈ {1, . . . , N }, f i (x) λ i }.
3.1. Assumptions and statement of the theorem. We will impose two conditions on the functions (f i ) 1 i N . As we will see in Section 4, both conditions are satisfied when these functions (f i ) are constructed from a semi-discrete optimal transport transport problem whose cost function satisfy Loeper's condition (see Definition 1.1).

Non-degeneracy. The functions (f i ) satisfy the nondegeneracy condition if the norm of their gradients is bounded from below:

ε nd := min 1 i N min X ∇f i > 0. (ND)
This condition is necessary for the continuity of the map Ĝ even when N = 1.

Transversality. The boundary of the convex set K(λ) can be decomposed into

N + 1 facets, namely (K(λ) ∩ ∂K i (λ i )) 1 i N and K(λ) ∩ ∂ X.
The purpose of the transversality condition we consider is to ensure that adjacent facets make a minimum angle, when λ remains close to some fixed vector λ 0 . More precisely, we say that the family of functions (f i ) satisfy the transversality condition near λ 0 if there exists positive constants ε tr and T tr 1 such that for every λ in R N satisfying λλ 0 ∞ T tr for the usual ∞ norm on R N and every point x in ∂K(λ) one has,

if ∃i = j ∈ {1, . . . , N } s.t. f i (x) = λ i and f j (x) = λ j , then, ∇f i (x) | ∇f j (x) ∇f i (x) ∇f j (x) 2 1 -ε 2 tr if ∃i ∈ {1, . . . , N } s.t. f i (x) = λ i and x ∈ ∂ X, then, ∀u ∈ N x X, u | ∇f j (x) u ∇f j (x) 2 1 -ε 2 tr , (T)
where N x X is the normal cone to the convex set X at x (see (B.68) for a definition). When ∂ X is smooth at x, N x X is the ray spanned by the exterior normal to X at x. We denote by Σ(λ) the set of points x in the boundary of K(λ) satisfying one of the assumptions in (T) or, equivalently, which belong to two (or more) distinct facets of K(λ):

Σ(λ) = 1 i N (K(λ) ∩ ∂ X ∩ ∂K i (λ i )) ∪ 1 i<j N (K(λ) ∩ ∂K i (λ i ) ∩ ∂K j (λ j )).
We will see that the transversality condition (T) and the quasi-convexity of the functions (f i ) imply a uniform upper bound on the (d -2)-Hausdorff measure of Σ(λ), which is essential in establishing the smoothness of Ĝ.

Theorem 3.1. Assume that the functions satisfy the non-degeneracy condition (ND) and the transversality condition (T) near λ 0 . Let ρ be a C α probability density on X. Then, the map Ĝ defined in (3.12) is of class C 1,α on the cube Q := λ 0 + [-T tr , T tr ] N and has partial derivatives given by

∂ Ĝ ∂λ i (λ) = K(λ)∩∂K i (λ) ρ(x) ∇f i (x) dH d-1 (x). (3.13)
In addition, the norm Ĝ

C 1,α (Q)
is bounded by a constant depending only on ε tr , ε nd , ρ C α ( X) , on the diameter of X and on

C M := max 1 i N ∇f i ∞ C L := max 1 i N ∇f i Lip( X) .
Note that the C 1,α constant of Ĝ depends on the transversality constant ε tr but that it does not depend on T tr .

3.2.

Existence of partial derivatives. Without loss of generality, we assume that λ 0 = 0. We start the proof of Theorem 3.1 by showing the existence of partial derivatives for the map Ĝ. In this section, we denote e 1 , . . . , e N the canonical basis of R N . We start by rewriting the finite difference defining the partial derivative of Ĝ in direction e i using the coarea formula. Fix λ < T tr . For t > 0, one has:

1 t ( Ĝ(λ + te i ) -Ĝ(λ)) = 1 t K(λ+te i )\K(λ) ρ(x)dH d (x) = 1 t λ i +t λ i ĝ(s)ds, (3.14)
where the function ĝ is defined by ĝ(s) :=

∩ j =i K j (λ j )∩f -1 i (s) ρ(x) ∇f i (x) dH d-1 (x). (3.15)
The same reasoning also holds for t < 0. We now claim that ĝ is continuous on some interval around λ i , which by (3.14) and the Fundamental Theorem of Calculus will imply that the limit as t → 0 in (3.14) exists and is equal to ĝ(λ), thus establishing the formula (3.13). The continuity of ĝ follows from the next proposition, which is formulated in a slightly more general way.

Proposition 3.2. Let σ be a continuous non-negative function on X and let ω be the modulus of continuity of σ. Given any vector λ in R N with λ ∞ T tr , consider the function

h : s ∈ R → L∩Ss σ(x)dH d-1 (x),
where L := j =i K j (λ j ) and S s := f -1 i (s). Then h is uniformly continuous on [-T tr , T tr ] and has modulus of continuity

ω h (δ) = const • (ω(Cδ) + |δ|), (3.16)
where the constants only depend on

f i C 1,1 , diam( X), ε nd , ε tr and σ ∞ .
Taking σ = ρ/ ∇f i in the previous proposition, which is continuous using the non-degeneracy condition (ND) and the assumption f i ∈ C 1,1 ( X), we see that the function ĝ defined by (3.15) is continuous. This implies the existence of partial derivatives and establishes formula (3.13). The proof of Proposition 3.2 requires two lemmas. Lemma 3.3. Assume that the functions f i : X → R satisfy (ND). Then, for every i ∈ {1, . . . , N }, there exists a map

Φ i : X × R → R d such that: (i) For any (x, t) in X × R such that the curve Φ i (x, [0, t]) remains in X, one has f i (Φ i (x, t)) = f i (x) + t. (ii)
The map Φ i satisfies the following inequalities for every x, y ∈ X, t ∈ R:

Φ i (x, t) -Φ i (x, s) |t -s| ε nd (3.17) Φ i (x, t) -Φ i (y, t) exp(C Φ |t|) x -y , (3.18)
where C Φ := 3C L /ε 2 nd . Proof. We consider the vector field V 0 i (x) = ∇f i (x)/ ∇f i (x) 2 on X, which satisfies V 0 i ∞ 1/ε nd and whose Lipschitz constant is bounded by C Φ . This vector field is extended to R d using the orthogonal projection on X, denoted p

X , ∀x ∈ R d , V i (x) := V 0 i (p X (x)))
. By convexity of X, the map p X is 1-Lipschitz. This implies that the Lipschitz constant of V i is also bounded by C Φ . We let Φ i be the flow induced by this vector field, which exists and is for all time since V i is bounded and uniformly Lipschitz on all of R d . The inequality (3.17) follows from the definition of integral curves and the bound on

V i . Any integral curve γ : [0, T ] → R d of V i which remains in X satisfies f i (γ(t)) = f i (γ(0)) + t 0 γ (s) | ∇f i (γ(s)) ds = f i (γ(0)) + t 0 V i (γ(s)) | ∇f i (γ(s)) ds = f i (γ(0)) + t,
thus establishing (i). The inequality (3.18) follows from the bound on the Lipschitz constant of V i and from Gronwall's lemma. Lemma 3.4. Assuming the transversality condition (T), there exists a constant depending only on d and diam( X) such that for every λ ∞ T tr ,

H d-2 (Σ(λ)) const(d, diam( X)) • 1 ε tr . S t S s x F (x) A t A s B t B s L ∂K 1 (λ 1 ) Λ ∂X ∂K 2 (λ 2 )
Figure 2. Illustration of the proof of Proposition 3.2.

Proof. Given λ ∞ T tr , the transversality condition (T) implies

∀x ∈ Σ(λ), ∃u, v ∈ N x K(λ), u | v u v 2 1 -ε 2 tr ,
where N x K(λ) is the normal cone to the convex set K(λ) at x (see (B.68)). This implies that Σ(λ) is included in the set Sing(K(λ), ε tr ) of τ -singular points defined in (B.69), with τ = ε tr . The lemma then follows from Proposition B.1.

Proof of Proposition 3.2. Let t, s be small enough so that the transversality condition (T) holds (that is t, s ∈ [-T tr , T tr ]). We assume that t < s so as to fix the signs of some expressions. We consider the following partition of the facet S t ∩ L, whose geometric meaning is illustrated in Figure 2:

A t = {x ∈ S t ∩ L | Φ i (x, [0, s -t]) ⊆ L} B t = {x ∈ S t ∩ L | ∃u ∈ [0, s -t), Φ i (x, u) ∈ ∂L}.
Similarly, we define

A s = {x ∈ S s ∩ L | Φ i (x, [t -s, 0]) ⊆ L} B s = {x ∈ S s ∩ L | ∃u ∈ (t -s, 0], Φ i (x, u) ∈ ∂L}.
Recall that by definition,

h(t) = At σ(x)dH d-1 (x) + Bt σ(x)dH d-1 (x), (3.19)
where the integral is with respect to the (d -1)-dimensional Hausdorff measure. Our strategy to show the continuity of h is to prove that the first term in the sums defining h(t) and h(s) in (3.19) are close, namely

At σ(x)dH d-1 (x) - As σ(x)dH d-1 (x) const • (|s -t| + ω(C |s -t|))
(3.20) and then that the terms involving B t , B s are small (recall that both sets depend on t and s):

Bt σ(x)dH d-1 (x) + Bs σ(x)dH d-1 (x) const • |s -t| (3.21)
The combination of the estimates (3.20) and (3.21) implies the desired inequality (3.16). We now turn to the proof of these estimates, and that the constant in these estimates depend on

f i C 1,1 , diam( X), ε nd , ε tr and σ ∞ .
Proof of Estimate (3.20). By Lemma 3.3.(i), for any point x in A t one has f i (Φ i (x, st)) = s, so that the map F (x) := Φ i (x, st) induces a bijection between the sets A t and A s . As a consequence of (3.18), the restriction of F to A t is a bi-Lipschitz bijection between the sets A t and A s , with Lipschitz constant max{

F -1 Lip(As) , F Lip(At) } exp(C Φ |s -t|). Using a Lipschitz change of variable formula, we get At σ(x)dH d-1 (x) = F -1 (As) σ(x)dH d-1 (x) F -1 d-1 Lip(As) As σ(F -1 (x))dH d-1 (x) exp(C Φ (d -1) |s -t|) As σ(F -1 (x))dH d-1 (x). (3.22)
By definition of the modulus of continuity and thanks to (3.17),

σ(F -1 (x)) -σ(x) ω( F -1 (x) -x ) = ω( Φ(x, s -t) -x ) ω(|s -t| /ε nd )
Integrating this inequality, we get

As σ(F -1 (x))dH d-1 (x) As σ(x)dH d-1 (x) + H d-1 (A s )ω(|s -t| /ε nd ) As σ(x)dH d-1 (x) + H d-1 ( X)ω(|s -t| /ε nd ) (3.23)
where the second inequality uses the monotonicity of the (d -

σ(x)dH d-1 (x) exp(C Φ (d-1) |s -t|) As σ(x)dH d-1 (x) + H d-1 ( X)ω(|s -t| /ε nd ) so that At σ(x)dH d-1 (x) - As σ(x)dH d-1 (x) (exp(C Φ (d -1) |s -t|) -1) σ ∞ H d-1 ( X)+ exp(C Φ (d -1) |s -t|)H d-1 ( X)ω(|s -t| /ε nd ) const • (|s -t| + ω(|s -t| /ε nd )),
where the constant depends on C L , ε nd , ε tr , σ ∞ and diam( X). Exchanging the role of s and t completes the proof of (3.20).

Proof of (3.21). By definition, for every point x in the set B t , the curve Φ i (x, [0, st]) must cross the boundary of L at some point, so that

u(x) := min{v ∈ [0, s -t] | Φ i (x, v) ∈ ∂L}
is well defined. We write P (x) := Φ i (x, u(x)) for the corresponding point on the boundary of L. By definition of u(x), the curve Φ(x, [0, u(x)]) is included in L, so that by Lemma 3.3.(i) we have f i (P (x)) = t + u(x). This shows

P (B t ) ⊆ Λ := ∂L ∩ f -1 i ([t, s]). (3.24)
We now prove that the map P satisfies a reverse-Lipschitz inequality. Note that for any point x in B t ,

x = Φ i (P (x), -u(x)) = Φ i (P (x), t -f i (P (x))),
Using the bounds (3.18) and (3.17), we get that for any x, y in B t ,

xy Φ i (P (x), tf i (P (x))) -Φ i (P (y), tf i (P (y)))

Φ i (P (x), tf i (P (x))) -Φ i (P (y), tf i (P (x)))

+ Φ i (P (y), t -f i (P (x))) -Φ i (P (y), t -f i (P (y))) exp(C Φ T tr ) P (x) -P (y) + |f i (P (x)) -f i (P (y))| /ε nd C P (x) -P (y)
where C := exp(C Φ ) + C L /ε nd ; we have used that T tr 1. We can now bound the (d-1)-Hausdorff measure of B t from that of Λ using this Lipschitz bound and the inclusion (3.24):

H d-1 (B t ) H d-1 (P -1 (P (B t ))) C d-1 H d-1 (Λ). (3.25)
What remains to be done is to prove that the (d -1)-Hausdorff measure of Λ behaves like O(|s -t|), and this is where the transversality condition will enter.

Let us write

F j := f -1 j (λ j ), j = 0, i, ∂ X ∩ ∂L, j = 0.
Then ∂L can be partitioned (up to a H d-1 -negligible set) into faces ∂L = ∪ j =i (F j ∩ L) and using the coarea formula on each of the facets we get

(writing B := f -1 i ([t, s])) H d-1 (Λ) = j =i H d-1 (B ∩ (F j ∩ L)) = j =i B∩(F j ∩L) dH d-1 (x) = j =i s t Su∩(F j ∩L) 1 J ij (x) dH d-2 (x)du, (3.26)
where and J ij (x) is the Jacobian of the restriction of the function f i to the hypersurface F j . More precisely,

J ij (x) = ∇f i (x) -∇f i (x) | ∇f j (x) ∇f j (x) ∇f j (x) 2
if j = 0, i, and

J i0 (x) = ∇f i (x) -∇f i (x) | v 0 (x) v 0 (x)
where v 0 (x) ∈ N x X is a unit vector. Since X is convex, for H d-1 a.e. x ∈ ∂ X, the normal cone N x X consists of only one direction, thus for such x there is a unique choice of v 0 (x). Let us write v i = ∇f i (x)/ ∇f i (x) and v j for either ∇f j (x)/ ∇f j (x) or v 0 (x), we then have using (T)

J ij (x) 2 = ∇f i (x) 2 v i -v i | v j v j 2 ∇f i (x) 2 (1 -v i | v j 2 ) ε 2 nd ε 2 tr . ( 3
H d-1 (Λ) 1 ε nd ε tr j =i s t H d-2 (S u ∩ (F j ∩ L))du = 1 ε nd ε tr s t H d-2 (S u ∩ ∂L)du. (3.28)
By definition, a point belongs to the intersection S u ∩ ∂L if it lies in the singularity set Σ(λ(u)) where λ(u) = (λ 1 , . . . , λ i-1 , u, λ i+1 , . . . , λ N ). By Lemma 3.4, 3.3. Continuity of partial derivatives. We prove that the function Ĝ defined in (3.12) is continuously differentiable by controlling the modulus of continuity of its partial derivatives given in (3.13). Again, we start with a slightly more general proposition. Proposition 3.5. Let σ be a continuous function on X with modulus of continuity ω. Consider the following function on the cube Q := [-T tr , T tr ] N :

H d-2 (S u ∩ ∂L) H d-2 (Σ(λ(u))) const(d, diam( X)) • 1 ε tr . ( 3
H(λ) := K(λ)∩f -1 i (λ i ) σ(x)dH d-1 (x).
Then H is uniformly continuous on Q with modulus of continuity

ω H (δ) = const • (ω(Cδ) + |δ|),
where the constants only depend on f i C 1,1 ( X) , diam( X), ε nd , ε tr , and σ ∞ .

Proof. Proposition 3.2 yields that the function H is uniformly continuous with respect to changes of the ith variable. Let us now consider variations with respect to the jth variable, with j = i by introducing

h : s ∈ [-T tr , T tr ] → K(λ 1 ,...,λ j-1 ,s,λ j+1 ,...,λ N )∩f -1 i (λ i ) σ(x)dH d-1 (x).
for some fixed λ ∈ [-T tr , T tr ] N . We can rewrite the difference between two values of h using the coarea formula. As before, we assume s > t to fix the signs and introduce

L := X ∩ k ∈{i,j} K k (λ k ) and S := f -1 i (λ i ). We have h(s) -h(t) = L ∩K j (s)∩S σ(x)dH d-1 (x) - L ∩K j (t)∩S σ(x)dH d-1 (x) = s t L ∩S∩f -1 j (u) σ(x) J ij (x) dH d-2 (x)du,
where the Jacobian factor J ij ε nd ε tr from (3.27). Therefore,

h(s) h(t) + σ ∞ ε nd ε tr s t H d-2 (L ∩ S ∩ f -1 (u))du.
(3.30)

Just as in the proof of Proposition 3.2, the set L ∩ S ∩ f -1 (u) is included in the set Σ(λ 1 , . . . , λ j-1 , u, λ j , . . . , λ N ). Thus, by Lemma 3.4,

H d-2 (L ∩ S ∩ f -1 (u)) const(d, X) ε tr . ( 3
C h := const(d, X) σ ∞ ε nd ε 2 tr . Finally, |H(µ) -H(λ)| N j=1 |H(λ 1 , . . . , λ k-1 , µ j , . . . , µ N ) -H(λ 1 , . . . , λ j , µ j+1 , . . . , µ N )| ω h (|µ i -λ i |) + j =i C h |µ j -λ j | ω h ( µ -λ ∞ ) + (N -1)C h µ -λ ∞ ,
where ω h is the modulus of continuity defined in Proposition 3.2. This establishes the uniform continuity of the function H, with the desired modulus of continuity.

3.4. Proof of Theorem 3.1. Proposition 3.2 shows that the partial derivative Ĝ with respect to the variable λ i exists and that its expression is given by (3.15). Applying Proposition 3.5 with σ(x) = ρ(x)/ ∇f i (x) , we obtain C α regularity for each of the partial derivatives of Ĝ on the cube Q := [-T tr , T tr ] N from the C α regularity of ρ. Moreover, the C α constant of each partial derivative over Q is controlled by const(diam( X), ε nd , ε tr , ∇f i C 1,1 (X) , ρ C α (X) ).

C 2,α regularity of Kantorovich's functional

This section is devoted to the proof of the following regularity result.

Theorem 4.1. Assume (Reg), (Twist), and (QC). Let X be a compact, c-convex subset of Ω and ρ in P ac (X) ∩ C α (X) for α in (0, 1]. Then, the Kantorovich's functional Φ is uniformly C 2,α on the set

K ε := {ψ : Y → R | ∀y ∈ Y, ρ(Lag y (ψ)) > ε} (4.32)
and its Hessian is given by (1.8). In addition, the C 2,α norm of the restriction of Φ to K ε depends only on ρ ∞ , ε, diam(X), and the constants defined in Remark 4.1 below.

For the remainder of the section, for any point y in Y , we will denote X y = (exp c y ) -1 (X) ⊆ R d the inverse image of the domain X in the exponential chart at y. The set X y is convex by c-concavity of X. We consider the functions f z,y : p ∈ X y → c(exp c y (p), y)c(exp c y (p), z), which are quasi-concave by (QC). The main difficulty in deducing Theorem 4.1 from Theorem 3.1 is in establishing the quantitative transversality condition (T) for the family of functions (f z,y ) z∈Y \{y} .

Remark 4.1 (Constants). The C 2,α norm of the restriction of Φ to K ε explicitely depends on the following constants, whose finiteness (or positivity) follows from the compactness of the domain X, from the finiteness of the set Y and from the conditions (Reg), (Twist), and (QC):

ε tw := min x∈X min y,z∈Y y =z Dc y (x) -Dc z (x) g > 0 C ∇ := max (x,y)∈X×Y Dc y (x) g < +∞ C exp := max y∈Y max exp c y Lip(Xy) , (exp c y ) -1 Lip(X) < +∞, (4.33)
where we recall that c y (x) = c(x, y) and X y := exp -1 y (X). Our estimates will also rely on the following constants involving the differential of the exponential maps. As before, the tangent spaces T x Ω are endowed with the Riemannian metric g from Ω. We set: Outline. In Section 4.1, we establish a part of the transversality condition using elementary properties of convex sets (Proposition 4.2). We establish in Section 4.2 a second transversality condition using additional assumptions and proceed in Section 4.3 to the proof of Theorem 4.1. In Section 4.4, we propose an alternative transversality estimate when Y is a sample subset of a target domain Ω (Proposition 4.8).

C cond := max
4.1. Lower transversality estimates. Next, we undertake a series of proofs to obtain explicit constants in the transversality estimate (T), which depend on the choices of cost, domains, and dimension. Consider the Laguerre cell of a point y in Y in its own exponential chart, that is

L y (ψ) := (exp c y ) -1 (Lag y (ψ)) = {p ∈ X y | f z,y (p) ψ(z) -ψ(y)}.
The set L y (ψ) is the intersection of sublevel sets of the functions f z,y , and is therefore a convex subset of X y by condition (QC). The first proposition establishes that two unit outer normals to L y with the same basepoint cannot be near-opposite. Recall the definition of the normal cone from (B.68). Proposition 4.2. Assume that ψ lies in K ε/2 (see (4.32)). For any y in Y , any point p in ∂L y (ψ) and any unit normal vectors v, w ∈ N p L y (ψ) one has

v | w -1 + δ 2 0 , (4.34)
where

δ 0 := ε/(2 d-1 ρ ∞ C exp 2d diam (X) d ) 1.
The proof of this proposition follows from a general lemma about convex sets. By convexity (QC), the set L y (ψ) is contained in an intersection of two half-spaces with outward normals v and w at p, giving an upper bound on its volume in term of its diameter and the angle between v and w (see Figure 3). On the other hand, we know that the volume of L y (ψ) is bounded from below by a constant depending on ε. Comparing these bounds will give us the one-sided estimate (4.34).

Lemma 4.3. Let K be a bounded convex set of R d , let p be a boundary point of K and v, w be two unit (outward) normal vectors to K at p. Then, Proof. The left-hand side of the inequality is non-positive, so the inequality needs only to be proven when v | w 0, which we assume from now on. Making a rotation of axes and a translation if necessary, we assume that p lies at the orign and that the unit vectors span the first two coordinates of R d . Then, letting

-1 + δ 2 K v | w where δ K = H d (K) 2 d-2 diam(K) d 1.
H := {p | p | v 0}, H := {p | p | w 0}
and D be the two-dimensional disc centered at 0 of radius diam(K), one has

K ⊆ H ∩ H ∩ (D × [-diam(K), diam(K)] d-2 ).
The intersection H ∩ H ∩ D is an angular sector of the disc D, whose angle is equal to θ := πarccos( v | w ) (see Figure 3). Therefore, we have

H d (K) H d (H ∩ H ∩ (D × [-diam(K), diam(K)] d-2 )) 2 d-2 diam(K) d tan(θ/2). (4.35)
Using the expression of cos(θ) in term of tan(θ/2) and recalling v | w 0,

tan(θ/2) = 1 + v | w 1 -v | w 1 + v | w (4.36)
The lemma follows directly from Equations (4.35)-(4.36).

Proof of Proposition 4.2. By definition of the bi-Lipschitz constant C exp ,

H d (L y (ψ)) ε/(2C d exp ρ ∞ ) and diam(L y (ψ)) C exp diam(X).
Applying the above lemma to the two outward normals v, w at p, we get

v | w + 1 H d (L y (ψ)) 2 4 d-2 diam(L y (ψ)) 2d ε 2 4 d-1 C 4d exp ρ 2 ∞ diam(X) 2d
.

We also record the following lemma for later use.

Lemma 4.4. Let y in Y and let p in be a point of L y (ψ) such that for some z = y, f z,y (p) = ψ(z)-ψ(y). Then, the point p := (exp c z ) -1 (exp c y (p)) belongs to L z (ψ) and the vector ∇f y,z (p) lies in the normal cone N p L z (ψ).

Proof. We introduce the point x = exp c y (p). The hypothesis is equivalent to c(x, y) + ψ(y) = c(x, z) + ψ(z). Since p belongs to L y (ψ), the point x belongs to Lag y (ψ). Then, for any z ∈ Y , c(x, z) + ψ(z) = c(x, y) + ψ(y) c(x, z ) + ψ(z ) thus establishing that x belongs to Lag z (ψ) or equivalently p ∈ L z (ψ). 4.2. Upper transversality estimates. We now turn to the proof of the quantitative transversality estimates. We begin with a bound which involves the condition number of differential of exponential maps, see Remark 4.1. The advantage of this first bound is that we do not have to assume that the points in Y are sampled from a continuous domain. A second transversality estimate is presented in §4.4.

Notation. We introduce notation that will be used in throughout this section. We fix a point y 0 in and also fix an arbitrary ordering of the remaining points, so that Y = {y 0 , y 1 , . . . , y N }. We define X := X y 0 and for every index i ∈ {1, . . . , N } we put

f i := f y i ,y 0 : p ∈ X → c(exp c
y 0 (p), y 0 )c(exp c y 0 (p), y i ). By the (Twist) condition, these functions f 1 , . . . , f N satisfy the non-degeneracy condition (ND), and we have the following inequalities: We also consider the same family of convex set as in Section 3:

ε nd := min i,j =0 min p∈Xy 0 ∇f i (p) -∇f j (p) C -1 exp ε tw > 0, (4.37)
K(λ) = {p ∈ X | ∀1 i N, f i (p) λ i }, so that K(λ ψ ) = (exp c y 0 ) -1 (Lag y 0 (ψ)).
Proposition 4.5. Assume that λ := λ ψ where ψ belongs to K ε/2 and let p be a point in K(λ). Then,

Case I: If f i (p) = λ i and f j (p) = λ j for i = j in {1, . . . , N }, then ∇f i (p) | ∇f j (p) ∇f i (p) ∇f j (p) 2 1 -δ 2 1 . (4.40) Case II: If p ∈ ∂ X and if f i (p) = λ i for some i in {1, . . . , N }, then ∀w ∈ N p X, ∇f i (p) | w ∇f i (p) w 2 1 -δ 2 1 . (4.41)
In the above inequalities,

δ 1 := ε nd δ 0 2C exp C ∇ C 2 cond .
By assumption, all the Laguerre cells associated to ψ contain a mass of at least ε/2. This allows us to apply Proposition 4.2, ensuring that normal vectors cannot be near-opposite, to all the Laguerre cells in their exponential charts. We denote L i := L y i (ψ) = (exp c y i ) -1 (Lag y i (ψ)) for brevity. The proposition also relies on two simple lemmas. The first lemma shows the effect of a diffeomorphism on the normal cone to a convex set, when its image is also convex. Lemma 4.6. Let K ⊂ R d be a compact, convex set, F be a C 1 diffeomorphism from an open neighborhood of K to an open subset of R d , and assume that F (K) is also a convex set. Then, for any point x in ∂K one has

N F (x) (F (K)) = [DF -1 F (x) ] * (N x K)
, where A * denotes the adjoint of A.

Proof. Consider x in ∂K, v ∈ N x K, and define ϕ(z) := F -1 (z)x | v . Since v is an outer normal to K at x, the restriction of ϕ to the set F (K) is non-positive. Since F (K) is convex, for any point y ∈ K, F (K) contains the segment [F (x), F (y)]. We therefore have

0 ϕ((1 -t)F (x) + tF (y)) ϕ(F (x)) + t ∇ϕ(F (x)) | F (y) -F (x) -o(t) = t [DF -1 F (x) ] * (v) | F (y) -F (x) -o(t).
where we have used ϕ(F (x)) = 0 and ∇ϕ(F (x)) = [DF -1 F (x) ] * (v) to obtain the equality at the end. Dividing by t and taking the limit as t goes to zero, we see that

∀y ∈ K, [DF -1 F (x) ] * (v) | F (y) -F (x) 0,
thus showing that [DF -1 F (x) ] * (v) belongs to the normal cone to F (K) at F (x). The converse inclusion follows from the symmetry of the problem.

The second lemma compares the angle between two vectors and the angle between their image under a linear map, using the generalized Wiedlandt inequality (see [START_REF] Householder | The theory of matrices in numerical analysis[END_REF]Section 3.4]). We identify R d with its tangent and cotangent spaces through the Euclidean structure. We denote the adjoint of the derivative of the exponential map exp c y at a point p in X y by

(D exp c y i ) * | p : T * exp c y i (p) Ω → T * p R d ∼ = R d , Lemma 4.7. Let y k = y ∈ Y , let x be a point in X and set p k := (exp c y k ) -1 (x) and p := (exp c y ) -1 (x) and A = (D exp c y k | p k ) * • [(D exp c y l | p l ) * ] -1 : T * p R d → T * p k R d .
Then, the following inequalities hold for all v, w in R d :

C -4 cond 1 + v | w v w 1 + Av | Aw Av Aw C 4 cond 1 + v | w v w .
Proof. Indeed, let θ be the angle between v and w and θ be the angle between Av and Aw, both in the interval (0, π). Let t := tan(θ/2), t := tan(θ /2). The generalized Wiedlandt inequality in [13, Section 3.4] asserts (1/ cond(A))t t cond(A)t. Expressing cos(θ) in term of t = tan(θ/2),

1 + cos(θ ) = 1 + 1 -t 2 1 + t 2 = 2 1 + t 2 cond(A) 2 (1 + cos(θ))
We conclude the second inequality above by using cond(A * 2 [A * 1 ] -1 ) cond(A 1 ) cond(A 2 ) and the definition of the constant C cond . For the first inequality, simply note that cond(A -1 ) = cond(A).

Proof of Proposition 4.5, case I. We let

V : = ∇f i (p) = ∇f y i ,y 0 (p), W : = ∇f j (p) = ∇f y j ,y 0 (p), v : = V V , w := W W .
Switching the indices i and j if necessary, we assume that V W . The proof depends on the sign of

W -V | V . Assume first W -V | V 0, and let α v := 1/ V , α w := 1/ W . Then, 1 -v | w = 1 2 v -w 2 = 1 2 α w (W -V ) -(α v -α w )V 2 = 1 2 α 2 w W -V 2 + 1 2 (α v -α w ) 2 V 2 -α w (α v -α w ) W -V | V
Using α w α v , and W -V ε nd we end up with

1 -v | w 2 1 -v | w 1 2 α 2 w W -V 2 1 2 ε 2 nd C 2 exp C 2 ∇ ε 2 nd δ 2 0 4C 2 exp C 2 ∇ C 4 cond = δ 2 1
where we have used (4.37) and (4.38), δ 0 1 and C cond 1. This establishes the desired bound when v | w ∈ [0, 1]. In the case v | w ∈ [-1, 0], we can apply Proposition 4.2 to show that 1

-v | w 2 1 + v | w δ 2 0 δ 2 1 , thus establishing the desired bound. Now suppose W -V | V 0. A slightly tedious computation gives v | w 2 = 1 - W -V 2 W 2 + W -V | v 2 W 2 = 1 - W -V 2 W 2 1 - W -V | v 2 W -V 2 1 - ε 2 nd C 2 exp C 2 ∇ 1 - W -V W -V | v , (4.42)
where we have used W -V | V 0 with (4.37) and (4.38) to get the last inequality. We will now apply Proposition 4.2 to the Laguerre cell L i . By Lemma 4.4, the point p i := (exp c y i ) -1 (exp c y 0 (p)) ∈ X y i belong to L i and the vectors V i := ∇f y 0 ,y i (p i ) and W i := ∇f y j ,y i (p i ) are both normals to L i at p i . Proposition 4.2 then shows that the vectors V i and W i satisfy

-1 + δ 2 0 V i | W i V i W i . ( 4

.43)

We transfer this inequality to the exponential chart of the original point y 0 using the linear map

A := ( D exp c y 0 p ) * • [( D exp c y i p i ) * ] -1 .
First, note that W -V = AW i and V = -AV i . Applying the generalized Wiedlandt inequality (Lemma 4.7) and ( 4.43) we have

1 - W -V | v W -V = 1 + AW i | AV i AW i AV i C -4 cond 1 + V i | W i V i W i C -4 cond δ 2 0 δ 2 1 . (4.44)
Combining this inequality with (4.42) we obtain (4.40) in this case as well.

Proof of Proposition 4.5, case II. Consider V := ∇f i (p) and let W be any vector in the normal cone N p X. When V | W 0, the inequality directly follows from Proposition 4.2, ensuring that normal vectors cannot be nearopposite. We now assume V | W 0 and we will apply Proposition 4.2 to the Laguerre cell of y i and transfer the result to the exponential chart of the point y 0 . Let p i = (exp c y i ) -1 (exp c y (p)). Then, by Lemma 4.4, p i belongs to L i and V i := ∇f y 0 ,y i (p i ) is a normal vector to L i at p i . We define a second normal vector by considering

A := ( D exp c y 0 p 0 ) * • [( D exp c y i p i ) * ] -1
and by setting

W i := A -1 W ∈ T * p i R d
. By Lemma 4.6, the vector W i belongs to the normal cone to X y i at p i . Moreover, since L i is contained in X y i and both sets contain p i , we have N p i X y i ⊆ N p i L i , thus ensuring that W i also belongs to the normal cone to L i at p i . Then, by Proposition 4.2 again,

V i | W i V i W i -1 + δ 2 0 .
As before, we transfer this inequality to the exponential chart of the original point y using the linear map A. We have V = ∇f i (p) = -AV i , and by construction W = AW i . We get the desired inequality by applying Lemma 4.7: where we have used the change-of-variable formula with x = exp c y 0 (p), so that ρ is the density of the pushforward measure (exp c y 0 ) -1 # (ρH d g ) with respect to the Lebesgue measure. We recall that K(λ ψ ) = (exp c y 0 ) -1 (Lag y 0 (ψ)), so that G y 0 (ψ) = Ĝ(λ ψ ) (as defined in (3.12)). The differentiability of Ĝ will be proven using Theorem 3.1 from the previous section.

1 - V | W V W = 1 + AV i | AW i AV i AW i C -4 cond 1 + V i | W i V i W i C -4 cond δ 2 0 δ 2
Let us fix a function ψ 0 in K ε and recall that λ 0 := λ ψ 0 . By Proposition 2.3 there exists a positive constant T tr such that every function ψ on Y satisfying ψψ 0 ∞ T tr belongs to K ε/2 . Then, by Proposition 4.5, we see that the functions f i satisfy the transversality condition (T) on the cube λ 0 + [-T tr , T tr ] N with constant

ε tr = δ 1 = ε nd δ 0 2C exp C ∇ C 2 cond , (4.45)
where we recall that

δ 0 = ε/(2 d ρ ∞ C exp 2d diam (X) d ).
Note also that since ρ is α-Hölder and since the exponential map is C 1,1 , the probability density ρ is also α-Hölder with constant

ρ C 0,α ( X) const( ρ C 0,α , C det ). ( 4

.46)

We can now apply Theorem 3.1. This ensures that the function Ĝ is of class C 1,α on the cube λ 0 +[-T tr , T tr ] N , so that ∂Φ/∂1 y 0 is C 1,α on a neighborhood of ψ 0 . Since this holds for any point y 0 ∈ Y and any function ψ 0 in K ε , we have established the C 2,α -regularity of Φ on K ε . The claimed dependency of Φ C 2,α (K ε ) follows from equations (4.45)-( 4.46) and from Theorem 3.1.

Our goal is now to deduce the formula for the gradient of G given in Theorem 4.1 (Equation (1.8)), from the formula for the gradient of Ĝ given in Theorem 3.1 (Equation (3.13)). This is done by looking more closely at the change of variable induced by the exponential map F := exp c y 0 : Ω → R d . For ease of notation we let h := c y 0c y i = f i • F -1 . By definition of the push-forward, we have for any bounded measurable function χ on Ω,

Ω χ(F (p))ρ(p)dH d (p) = Ω χ(x)ρ(x)dH d g (x).
Multiplying χ by the characteristic function of h -1 ([t, s]), this gives

f -1 i ([t,s]) χ(F (p))ρ(p)dH d (p) = h -1 ([t,s]) χ(x)ρ(x)dH d g (x).
Applying the coarea formula on both sides, we get

s t f -1 i (r) χ(F (p))ρ(p) ∇f i (p) dH d-1 (p)dr = s t h -1 (r) χ(x)ρ(x) ∇h(x) g dH d-1 g (x)dr. ( 4
.47) Using the C 1,1 smoothness of the functions f i and the (Twist) condition, we can see that for any χ in C 0 c (Ω), the two inner integrals

r → f -1 i (r) χ(F (p))ρ(p) ∇f i (p) dH d-1 (p) and r → h -1 (r) χ(x)ρ(x) ∇h(x) g dH d-1 g (x)
depend continuously on r. Using the continuity of these functions in r, equation (4.47) and the Fundamental Theorem of Calculus, we get that for any function χ in C 0 c (Ω) and any r in R,

f -1 i (r) χ(F (p))ρ(p) ∇f i (p) dH d-1 (p) = h -1 (r) χ(x)ρ(x) ∇h(x) g dH d-1 g (x)
.

By Tietze's extension theorem, every function in C 0 c (S) can be extended to a function in C 0 c (Ω). The previous equality therefore holds for any χ in C 0 c (S), and by density, it also holds for any function χ in L 1 (S). Applying this with χ equal to the indicator function of the interface between the Laguerre cell of y 0 and the cell of y i , we get the desired formula for the partial derivatives:

∂G i ∂ψ i (ψ) = ∂ Ĝ ∂λ i (λ ψ ) = Ly 0 (ψ)∩f -1 i (ψ(y i )-ψ(y 0 )) ρ(p) ∇f i (p) dH d-1 (p) = Lag y 0 (ψ)∩Lag y i (ψ) ρ(x) Dc y 0 (x) -Dc y i (x) g dH d-1 g (x).
4.4. Alternative upper transversality estimates. Finally, we state an alternate upper transversality estimate, under the assumption that the points in Y are sampled from some target domain Λ, along with some convexity conditions. Specifically, let Λ be a bounded, open subset in some Riemannian manifold, with Y ⊂ Λ. We then assume that for any x ∈ Ω cl , the mapping

y → -D x c(x , y)
is a diffeomorphism onto its range, and we denote the inverse by exp c x . We will also assume that (exp c x ) -1 (Λ) is convex for all x ∈ Ω, and finally that the following inequality holds: for any x, x ∈ Ω, p 0 , p 1 ∈ (exp c

x ) -1 (Λ), and

t ∈ [0, 1], -c(x, exp c x ((1 -t)p 0 + tp 1 )) + c(x , exp c x ((1 -t)p 0 + tp 1 ))) max{-c(x, exp c x (p 0 )) + c(x , exp c x (p 0 )), -c(x, exp c x (p 1 )) + c(x , exp c x (p 1 ))} (4.48)
For more on these conditions, see Remark 4.3 below. Proposition 4.8 can be applied to provide an alternative bound in the transversality condition (T) when the point p 0 ∈ ∂K(λ) is in the interior of X (so in particular, when dealing with Laguerre cells that do not intersect ∂X). The advantage of this bound is that it does not require knowledge of the condition number C cond .

We also slightly re-define the constants C ∇ and C exp so that in their definitions, the maximum of y ranges over the domain Λ instead of just Y .

Proposition 4.8. Suppose λ < T tr < ε 8C 2d-2 exp ρ ∞ H d-1 g (∂X)
, and p 0 ∈ K(λ)

with f i (p 0 ) = λ i and f j (p 0 ) = λ j for i = j in {1, . . . , N }. Then we have ∇f i (p 0 ) | ∇f j (p 0 ) ∇f i (p 0 ) ∇f j (p 0 ) 2 1 -δ 2 2 (4.49)
where

δ 2 := εε nd 4 √ 2C 2 ∇ C 2d+4 exp ρ ∞ H d-1 g (∂X)
.

The idea behind this upper bound is the following. First suppose ∇f i (p 0 ) is shorter than ∇f j (p 0 ) and the two vectors are collinear, and that λ = 0. Inspired by calculations from [12, Remark 2.5, Proof of Lemma 4.7], this would allow us to use (4.48) to obtain the inequality f i max{0, f j } everywhere, which in turn would mean that -c(•, y i ) + ψ(y i ) is less than one of -c(•, y 0 ) + ψ(y 0 ) or -c(•, y j ) + ψ(y j ) everywhere. However, combined with (Twist), this would imply that the Laguerre cell Lag y i (ψ) has measure zero, which is a contradiction. In the general case, we rotate ∇f i (p 0 ) so that it is collinear with ∇f j (p 0 ), and use the lower bound on Lag y i (ψ) to obtain a quantitative version of the above argument. We relegate the proof of the proposition to the Appendix in Section C. Remark 4.3. Under a set of standard conditions, we can obtain both (QC) and (4.48).

Let Ω and Λ be bounded and smooth domains in d dimensional Riemannian manifolds and take a cost c ∈ C 4 (Ω × Λ). Also assume

• c satisfies the (Twist) condition: for every x ∈ Ω, the map y ∈ Λ → -D x c(x, y) is a diffeomorphism onto its image Λ x := -D x c(x, Λ) and we define the c-exponential map exp c

x : Λ x → Λ by exp c x = (-D x c(x, •)) -1 .
• the cost c * (x, y) := c(y, x) satisfies the (Twist) condition: then for every y ∈ Λ, we can define the c * -exponential map exp c y :

Ω y → Ω by exp c y = (-D y c(•, y)) -1 . • (exp c x ) -1 (Λ) is convex for each x ∈ Ω. • (exp c y ) -1 (Ω) is convex for each y ∈ Λ. • det D 2 xy c(x, y) = 0 for all (x, y) ∈ Ω × Λ. • For any (x, y) ∈ Ω × Λ and η ∈ T * x Ω, V ∈ T x Ω with η(V ) = 0, -(c ij,pq -c ij,r c r,s c s,pq )c p,k c q,l V i V j η k η l 0, (A3w)
here indices before a comma are derivatives on Ω and after a comma on Λ, for fixed coordinate systems, and a pair of raised indices denotes the inverse of a matrix. This last condition (A3w) originates (in a stronger version) in [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF] related to regularity of optimal transport. [19, Theorem 3.2] in the Euclidean case and [16, Theorem 4.10] in the general manifold case show the above conditions imply (QC) and (4.48). In fact, they are equivalent as seen in [START_REF] Loeper | On the regularity of solutions of optimal transportation problems[END_REF]. This geometric interpretation is a key ingredient in showing regularity in the optimal transport problem in the vein of Caffarelli's classical work [START_REF] Caffarelli | The regularity of mappings with a convex potential[END_REF], see [START_REF] Figalli | Hölder continuity and injectivity of optimal maps[END_REF][START_REF] Guillen | On the local geometry of maps with c-convex potentials[END_REF].

Strong concavity of Kantorovich's functional

We establish in this section the strong concavity of Kantorovich's functional Φ over some suitable domain of R Y . As explained in the introduction, Φ is invariant under addition of a constant, so that we must restrict ourselves to the orthogonal complement E Y of the space of constant functions. Moreover, we will consider the set K ε defined by (4.32), which can be thought of as the space of strictly c-concave functions.

Theorem 5.1. Assume (Reg), (Twist), (QC). Let X be a compact, c-convex subset of Ω, and ρ be a continuous probability density on X satisfying (PW). Then, there exists a positive constant C, such that

∀ψ ∈ K ε , ∀v ∈ E Y , D 2 Φ(ψ)v | v -C • ε 3 v 2 ,
where C is defined in (5.58), and depends on ρ ∞ , H d-1 g (∂X), and C exp , C ∇ , and ε tw from Remark 4.1.

Remark 5.1. Note that unlike the domain X, the support of the density ρ does not need be c-convex. We provide in Appendix A an example of a radial measure on R d whose support is an annulus (hence is not simply connected) but whose Poincaré-Wirtinger constant (PW) is nonetheless positive.

The end of the section is devoted to the proof of Theorem 5.1. It relies on the fact that -D 2 Φ(ψ) can be regarded as the Laplacian matrix of a weighted graph on Y , whose first nonzero eigenvalue can be controlled from below using the Cheeger constant of the weighted graph. In turn, this weighted Cheeger constant can be controlled using the Poincaré-Wirtinger inequality.

5.1. Poincaré inequality and continuous Cheeger constant. We start by proving that the finiteness of the Poincaré-Wirtinger constant of the weighted domain (X, ρ) implies the positivity of the weighted Cheeger constant, defined in (5.50). In the following, a Lipschitz domain denotes the closure of an open set with Lipschitz boundary. Lemma 5.2. Assume (QC) and that X is compact and c-convex. Then (i) X is a Lipschitz domain, (ii) for any ψ ∈ K + and y in Y , Lag y (ψ) ∩ X is a Lipschitz domain.

Proof. By assumption, for any y ∈ Y one can write X = exp c y (X y ) where X y is a bounded convex subset of R d which must have nonempty interior since it supports an absolutely continuous probability measure. Moreover, the map exp c y is a diffeomorphism, hence bi-Lipschitz. This implies (i), while (ii) follows from the exact same arguments, remembering that ρ(Lag y (ψ)) > 0.

Given a Lipschitz domain A of X we denote

|∂A| ρ := ∂A∩int(X) ρ(x)dH d-1 g (x) and |A| ρ := A∩int(X) ρ(x)dH d g (x).
Lemma 5.3. Let X be a compact domain of Ω and ρ in C 0 (X) be a probability density with finite Poincaré-Wirtinger constant (PW). Then the weighted Cheeger constant of (X, ρ) is positive, that is

h(ρ) := inf A⊆X |∂A| ρ min(|A| ρ , |X \ A| ρ ) 2 C pw , (5.50)
where the infimum is taken over Lipschitz domains A ⊆ int(X) whose boundary has finite H d-1 g -measure. The proof is based on properties of functions with bounded variation. For more details on this topic, we refer the reader to [START_REF] Attouch | Variational analysis in sobolev and bv spaces: applications to pdes and optimization[END_REF]. Although the discussion in the reference is on Euclidean spaces, the relevant results easily extend to the Riemannian case, as exp c y serves as a global coordinate system on all of Ω.

Proof. Let A be a Lipschitz domain A of int(X). Since A has a Lipschitz boundary with finite area, its indicator function χ A has bounded variation in int(X). By the density theorem [2, Theorem 10.1.2], there exists a sequence of C 1 -functions f n on int(X) that converges to χ A in the sense of intermediate convergence (whose definition is not important here). By (PW),

f n -E ρ (f n ) L 1 (ρ) C pw ∇f n L 1 (ρ) . Since intermediate convergence is stronger than L 1 convergence, the continu- ity of ρ implies lim n→∞ f n -E ρ (f n ) L 1 (ρ) = χ A -E ρ (χ A ) L 1 (ρ) = 2 |A| ρ |X \ A| ρ .
Note that we used the fact that ρ is a probability measure, i.e. ρ(X) = 1. Proposition 10.1.2 of [START_REF] Attouch | Variational analysis in sobolev and bv spaces: applications to pdes and optimization[END_REF] implies that the total variation measure |Df n | narrowly converges to |Dχ A |, which with the continuity of ρ implies that

Ω |Df n |ρdH d g converges to Ω |Dχ A |ρdH d g = |A| ρ . The relation |Df n | = ∇f n g dH d g then gives lim n→∞ ∇f n L 1 (ρ) |∂A| ρ .
Combining the previous equations together leads to the desired inequality.

5.2.

Cheeger constant of a graph. The goal of this section is to give a lower bound of the second eigenvalue of -D 2 Φ(ψ) in terms of the Cheeger constant of the weighted graph induced by this matrix. An unoriented weighted graph can always be represented by its adjacency matrix (w yz ) (y,z)∈Y 2 , a symmetric matrix with zero diagonal entries. We introduce a few definitions from graph theory, following the conventions of [START_REF] Friedland | On cheeger-type inequalities for weighted graphs[END_REF]. The (weighted) Cheeger inequality bounds from below the first nonzero eigenvalue of the Laplacian of a weighted graph, denoted λ(w), from its Cheeger constant and its minimal degree. The formulation we use can be deduced from Corollary 2.2 of [START_REF] Friedland | On cheeger-type inequalities for weighted graphs[END_REF] and from the inequality 1-

√ 1 -x 2 x 2 /2.
Theorem 5.4 (Cheeger inequality). λ(w) 1 2 h 2 (w) • min y∈Y d y . We now proceed to the proof of the main theorem of this section.

w yz := - ∂ 2 Φ ∂1 y ∂1 z (ψ) = Lag y,z (ψ) ρ(x) Dc y (x) -Dc z (x) g dH d-1 g (x)
for y = z in Y , and with zero diagonal entries (w yy = 0). In the formula above, we used the notation Lag y,z (ψ) = Lag y (ψ) ∩ Lag z (ψ) for the facet between two Laguerre cells. By construction, the Laplacian matrix of this weighted graph is the Hessian matrix -D 2 Φ(ψ), so that Theorem 5.4 directly gives us a lower bound on the first nonzero eigenvalue of -D 2 Φ(ψ). To complete the proof, we need to bound the Cheeger constant and the minimum degree of the graph w from below.

Step where the second inequality comes from the fact that the facets Lag y,z (ψ) form a partition of the boundary ∂ Lag y (ψ) ∩ int(X) up to a H d-1 g -negligible set. To see that this is the case, it suffices to remark that in the exponential chart of y, the intersection of two distinct facets adjacent to y has a finite H d-2 g -measure, as implied by Lemma 3.4.

In order to apply the (continuous) Cheeger inequality, we need to replace the weighted area of the boundaries of Laguerre cells in (5.54) by the weighted volume of the cells. We have

H d-1 g (∂ Lag y (ψ)) C d-1 exp H d-1 ((exp c y ) -1 ∂ Lag y (ψ)) C d-1 exp H d-1 (∂X y ) C 2(d-1) exp H d-1
g (∂X). The first and third inequalities use the definition of the bi-Lipschitz constant C exp of the exponential map, while the second inequality uses the monotonicity of the H d-1 -measure of the boundary of a convex set with respect to inclusion (see [27, p.211]). Using the assumption Lag y (ψ) ρ ε, this gives us a (rather crude) reverse isoperimetric inequality

∂ Lag y (ψ) ρ ρ ∞ H d-1 g (∂ Lag y (ψ)) 1 ε C 2(d-
|S| w = y∈S d y 1 ε ρ ∞ C 2(d-1) exp H d-1 g (∂X) ε tw |A| ρ .
The same inequality holds for the complement |X \ S|. We combine the previous inequality with Equation (5.53) and with Lemma 5.3 to get a lower bound on the Cheeger constant h(w)

ε tw ε C 2(d-1) exp C ∇ H d-1 g (∂X) ρ ∞ C pw .
(5.56)

Note that, in order to apply Lemma 5.3 we implicitly used the fact that A is a Lipschitz domain (as a finite union of Lipschitz domains, see Lemma 5.2) whose boundary has finite H d-1 g -measure (by Equation (5.55)).

Step 2. In order to apply the Cheeger inequality, we still need to bound from below the weighted degree d y . By (5.51) one has, using the crucial fact that ∂ Lag y (ψ) ρ is the measure of ∂ Lag y (ψ) ∩ int(X),

d y = z =y w yz 1 2C ∇ z =y Lag y,z (ψ) ρ 1 2C ∇ ∂ Lag y (ψ) ρ .
Taking A = Lag y (ψ) in the definition of the Cheeger constant h(ρ), one gets ∂ Lag y (ψ) ρ h(ρ) min( Lag y (ψ) ρ , X \ Lag y (ψ) ρ ) h(ρ)ε.

The last inequality comes from the assumption that each Laguerre cell has a mass greater than ε and that X \ Lag y (ψ) also contains a Laguerre cell (except for the trivial case where Y is a singleton). We deduce

d y ε C ∇ C pw .
(5.57)

Step 3. Combining the Cheeger inequality with Equation (5.56) and (5.57) we have λ(w) Cε 3 where

C := ε 2 tw 2C 4(d-1) exp C 3 ∇ H d-1 g (∂X) 2 ρ 2 ∞ C 3 pw .
(5.58)

Since the graph induced by the Hessian is connected, the kernel of -D 2 Φ(ψ) is equal to the space of constant functions over Y , implying that Ker(-D 2 Φ(ψ)) ⊥ = E Y . Then, using the variational characterization of the first nonzero eigenvalue of the Laplacian matrix we get:

Cε 3 λ(w) = min v∈E Y -D 2 Φ(ψ) | v v 2 .

Convergence of the damped Newton algorithm

The goal of this section is to show the convergence of the Damped Newton algorithm for semi-discrete optimal transport. This follows in fact from a more general result. We establish in Section 6.1 the convergence of the damped Newton algorithm (Algorithm 1) under general assumptions on the functional. We finally apply this algorithm to the semi-discrete optimal transport problem, using the intermediate results (regularity and strict concavity of the Kantorovich functional) proven in Section 4 and 5.

6.1. General damped Newton algorithm. Let Y be a finite set and denote Y R the space of functions on Y . We consider P(Y ), the space of probability measures on Y , as a subset of Y R . Finally, we denote by E Y the space of functions on Y who sum to zero. In this section, we show that Algorithm 1 can be used to solve non-linear equations G(ψ) = µ where µ ∈ P(Y ) and the map G : R Y → P(Y ) satisfies some regularity and monotonicity assumptions. Proposition 6.1. Let G be a functional from R Y to P(Y ) which is invariant under addition of a constant. Let G(ψ) = y∈Y G y (ψ)1 y and

K ε = {ψ ∈ R Y | ∀y ∈ Y, G y (ψ) ε},
and assume that G has the following properties:

(i) (Regularity) For every positive ε, G is C 1,α on K ε . We let L ε be the smallest constant such that

∀ϕ = ψ ∈ K ε , G(ϕ) -G(ψ) ϕ -ψ + DG(ϕ) -DG(ψ) ϕ -ψ α L ε .
(ii) (Uniform monotonicity) For every ε > 0, there exists a positive constant κ ε such that G is κ ε -uniformly monotone on

K ε ∩ E Y : ∀ψ ∈ K ε , ∀v ∈ E Y , v | DG(ψ)v κ ε v 2 .
Now, let µ ∈ P(Y ) and let ψ 0 be a function on Y such that the constant ε 0 defined in (1.6) is positive. Set κ := min(κ 1 2 ε 0 , 1) and L := max(L 1 2 ε 0 , 1). Then, the iterates (ψ k ) of Algorithm 1 satisfy

G(ψ k+1 ) -µ (1 -τ k /2) G(ψ k ) -µ where τ k := min κ 1+ 1 α ε d 1 α L 1 α G(ψ k ) -µ , 1 .
(6.59)

In addition, as soon as τ k = 1 one has

G(ψ k+1 ) -µ L G(ψ k ) -µ 1+α κ 1+α .
In particular, the damped Newton's algorithm converges globally with linear speed and locally with superlinear speed (quadratic speed if α = 1).

Proof. We set ε := ε 0 , L := max(L ε/2 , 1) and κ := min(κ ε/2 , 1). First, we remark that for every ψ ∈ K ε 2 , the pseudo-inverse DG + (ψ) maps the subspace E Y to itself. The uniform monotonicity of G therefore implies that DG + (ψ) 1/κ, where • is the operator norm on R Y .

We start by the analysis of a single iteration of the algorithm. We let ψ := ψ k ∈ K ε , define v := DG(ψ) + (G(ψ)µ) and ψ τ := ψτ v. Since the pseudo-inverse DG + (ψ) is 1/κ-Lipschitz, one has v G(ψ)µ /κ. Now let τ 1 be the largest time before the curve ψ σ leaves K ε/2 . In particular, ψ τ 1 lies at the boundary of K ε/2 , meaning that there must exist a point y in Y such that G y (ψ τ 1 ) = ε/2. This implies that G(ψ τ 1 ) -G(ψ) ε/2, and using the Lipschitz bound on G we obtain a lower bound on τ

1 ε 2 G(ψ τ 1 ) -G(ψ) Lτ 1 v Lτ 1 κ G(ψ) -µ .
This implies that τ 1 is necessarily larger than κε/(2L G(ψ)µ ). We now established that the curve τ → ψ τ remains in K ε/2 before time τ 1 , implying that the function τ ∈ [0, τ 1 ] → G(ψ τ ) is uniformly C 1,α . Applying Taylor's formula we get

G(ψ τ ) = G(ψ -τ DG(ψ) + (G(ψ) -µ)) = (1 -τ )G(ψ) + τ µ + R(τ ) (6.60)
where, using v = DG(ψ) + (G(ψ)µ), and the α-Hölder property for DG

R(τ ) = τ 0 (DG(ψ σ ) -DG(ψ))vdσ L α + 1 τ α+1 v α+1 L G(ψ) -µ (1+α) κ (1+α) τ (1+α) (6.61)
For every y ∈ Y , using that µ y 2ε (by (1.6)) and G y (ψ) ε, one gets

G y (ψ τ ) (1 -τ )G y (ψ) + τ µ y + R y (τ ) (1 + τ )ε -R(τ ) .
If τ is chosen such that such that R(τ ) τ ε we will have G y (ψ τ ) ε for all points y in Y and therefore ψ τ will belong to K ε . Thanks to our estimate on R(τ ) this will be true provided that

τ τ 2 := min τ 1 , κ 1+ 1 α ε 1 α L 1 α G(ψ) -µ 1+ 1 α .
Finally we establish the second inequality required by Step 2 of the Algorithm.

To do that, we subtract µ from both sides in (6.60) to obtain

G(ψ τ ) -µ = (1 -τ )(G(ψ) -µ) + R(τ ). (6.62) In order to get G(ψ τ ) -µ (1 -τ 2 ) G(ψ) -µ , it is sufficient to establish that R(τ ) τ 2 G(ψ) -µ .
Using the estimation on R(τ ) again, we see that it suffices to take

τ τ 3 := min τ 2 , κ 1+ 1 α L 1 α G(ψ) -µ 2 1 α , 1 .
Finally, using L 1, κ 1 and G(ψ)µ d (since G(ψ) and µ are probability measures), we can establish that τ 3 τ k where the value of τ k is defined in (6.59). This ensures the first estimate on the improvement of the error between two successive steps.

By this first estimate, we see that there exists k 0 such that τ k = 1 for k k 0 . When this happens, one can use (6.62) to get G(ψ k+1 )µ R(τ ) . We obtain the second estimate of the theorem by plugging in (6.61). 6.2. Proof of Theorem 1.5. Proposition 6.1 can be directly applied to the gradient of the Kantorovich functional, or more precisely to

G(ψ) := y∈Y ρ(Lag y (ψ))1 y = ∇Φ(ψ) + µ
In that case, the set K ε is given by

K ε = {ψ ∈ R Y | ∀y ∈ Y, ρ(Lag y (ψ)) ε}.
We have assumed that the probability density ρ is C α (X) where X is a cconvex, compact subset of Ω. Then, by Theorem 4.1, for any ε > 0, the map G is uniformly C1,α over K ε . This ensures that the (Regularity) condition of Proposition 6.1 is satisfied. Furthermore, since we also assumed that ρ satisfies a weighted Poincaré-Wirtinger inequality, we can apply Theorem 5.1 to see that the (Uniform monotonicity) hypothesis of Proposition 6.1 is also satisfied. Applying Proposition 6.1, we deduce the desired convergence rates for Algorithm 1. 6.3. Numerical results. We conclude the article with a numerical illustration of this algorithm, for the cost c(x, y) = xy 2 and for a piecewise-linear density. The source density is piecewise-linear over a triangulation of [0, 3] with 18 triangles (displayed in Figure 6.3). It takes value 1 on the boundary ∂[0, 3] 2 and vanishes on the square [1, 2] 2 . In particular, the support of this density is not simply connected and not convex. The target measure is uniform over a uniform grid 1 n-1 {0, . . . , n -1} 2 . Figure 6.3 displays the iterates of the Newton algorithm, which in this case takes 25 iterations to solve the optimal transport problem with an error equal to the numerical precision of the machine. The source code of this algorithm is publicly available 1 .

We finally note that recent progress in computational geometry would allow one to implement Algorithm 1 for the quadratic cost on R 3 , refining [START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3d[END_REF] or [START_REF] Fernando De Goes | Power particles: an incompressible fluid solver based on power diagrams[END_REF]. It should also be possible to deal with optimal transport problems arising from geometric optics, such as the far-field reflector problem [START_REF] Machado | Intersection of paraboloids and application to minkowski-type problems[END_REF], whose associated cost satisfies the Ma-Trudinger-Wang condition [START_REF]Regularity of optimal maps on the sphere: the quadratic cost and the reflector antenna[END_REF].

:= (1/ω d-1 ) S d-1 f (z)dz, S d-1 |f (z) -F d-1 | dH d-1 (z) c d S d-1 ∇f (z) dH d-1 (z) (A.63)
for some positive constant c d . The second inequality is a Poincaré-Wirtinger inequality on the segment [0, R] weighted by ρ. Given a function in C 1 ([0, R]), and letting

F 1 := R 0 f (r)ρ(r)dr/ R 0 ρ(r)dr, R 0 f (r) -F 1 ρ(r)dr c ρ R 0 f (r) ρ(r)dr (A.64)
for some positive constant c ρ depending only on ρ. The inequality (A.64) can be deduced from Theorem 2.1 in [START_REF] Acosta | An optimal poincaré inequality in L 1 for convex domains[END_REF] and from the concavity of ρ on [r, R].

Proof. We now proceed to the proof of the Poincaré-Wirtinger inequality for (X, ρ). Let f : B(0, R) → R be a function of class C 1 . By polar coordinates and the definition of ρ, one has

F := B(0,R) f (x)ρ(x)dH d (x) = R 0 1 ω d-1 r d-1 S d-1 (r) f (z)ρ(r)dH d-1 (z)dr = R 0 f (r)ρ(r)dr,
where the function f (r) is the mean value of f over the sphere S d-1 (r),

f (r) = 1 ω d-1 r d-1 S d-1 (r) f (z)dH d-1 (z) = 1 ω d-1 S d-1 f (rz)dH d-1 (z).
Using the triangle inequality and the relation between ρ and ρ we get

B(0,R) |f (x) -F | ρ(x)dH d (x) = R 0 S d-1 (r) |f (z) -F | ρ(z)dH d-1 (z)dr R 0 ρ(r) f (r) -F dr + R 0 ρ(r) r d-1 ω d-1 S d-1 (r) f (z) -f (r) dH d-1 (z)dr (A.65)
We first deal with the second term of the sum. Using the Poincaré-Wirtinger inequality on the sphere (A.63), we have (A.66) By the calculation of F above, we see F is also the mean value of f weighted by ρ. We can therefore control the first term of the upper bound of (A.65) using the Poincaré-Wirtinger inequality on the segment (A. The next proposition gives an upper bound on the (d -2)-Hausdorff measure of Sing(K, τ ). Note that a more general version of this result can be found in [START_REF] Hug | Generalized curvature measures and singularities of sets with positive reach[END_REF], with optimal constant that depends on Minkowski's quermassintegrals. We provide below a straightforward and easy proof based on the notions of packing and covering numbers. where the constant depends on d and the diameter of K.

Before giving the proof of this proposition, we recall that the covering number Cov(K, η) of a subset K ⊆ R d is the minimum number of Euclidean balls of radius η required to cover K. The packing number of a subset K is given by Pack(K, η) := max{Card(X) | X ⊆ K and ∀x = y ∈ X, xy η}.

We will use the following comparisons between covering and packing numbers: Step 1. We first calculate an upper bound on the covering number of the unit bundle UK := {(x, n) ∈ ∂K × S d-1 | n ∈ N x K}. Given a positive radius r, we denote by K r the set of points that are within distance r of K. By convexity, the projection map p K : R d → K, mapping a point to its orthogonal projection on K, is well defined and 1-Lipschitz. We consider π : ∂K r → U(K)

x → p K (x),

xp K (x)

xp K (x)

The map π is surjective and has Lipschitz constant L := 1 + 4/r 2 . We deduce an upper bound on covering number of UK from the covering number of the level set ∂K r :

Cov(U(K), ε) Cov ∂K r , ε L .
Now, consider a sphere S with diameter 2 diam(K) that encloses the tubular neighborhood K r with r := diam(K). The projection map p K r is 1-Lipschitz, and p K r (S) = ∂K r . Using the same argument as above, we have: Using the comparison between packing and covering numbers, this means that we can cover Sing(K, τ ) with N ε balls of radius ε, such that N ε const(d, diam(K))/(τ ε d-2 ). By definition of the Hausdorff measure, we have where we have used that if y = exp c x 0 (q), then ( D exp c x 0 q ) * = D exp c y -Dxc(x 0 ,y) . As a result we obtain We now make the following observation. Let us write X i := (exp c y i ) -1 (X). Then for any t, s > 0, we can estimate the volume of X i ∩ {f y 0 ,y i -t} ∩ {f y j ,y i -s} by H d X i ∩ {f y 0 ,y i -t} ∩ {f y j ,y i -s} H d X i ∩ {f y 0 ,y i 0} ∩ {f y j ,y i 0} -H d (X i ∩ {-t < f y 0 ,y i 0}) -H d X i ∩ {-s < f y j ,y i 0} .

Using that L i ⊂ {f y 0 ,y i 0} ∩ {f y j ,y i 0}, we can bound the first term from below as

H d X i ∩ {f y 0 ,y i 0} ∩ {f y j ,y i 0} ε C d exp ρ ∞ .
For the second term, by the coarea formula, we can write H d (X i ∩ {-t < f y 0 ,y i 0}) 0 -t X i ∩{fy 0 ,y i =z} 1 ∇f y 0 ,y i (p) dH d-1 (p)dz

tH d-1 (∂X i ) ε nd tC d-1 exp H d-1 g (∂X) ε nd ,
where to obtain the second line we have again used the fact that for every z ∈ R, the set X i ∩ {f y 0 ,y i = z} is contained in the boundary of a convex subset of X i in conjunction with [START_REF] Kitagawa | An iterative scheme for solving the optimal transportation problem[END_REF]Remark 5.2]. By a similar bound on the third term, we see that as long as max{t, s} < εε nd 2C 2d-2 exp ρ ∞ H d-1 g (∂X) we have H d X i ∩ {f y 0 ,y i -t} ∩ {f y j ,y i -s} > 0, thus in particular, (by continuity of f y 0 ,y i and f y j ,y i ) there must exist a point p c ∈ X i for which max{f y 0 ,y i (p c ), f y j ,y i (p c )} -

εε nd 2C 2d-2 exp ρ ∞ H d-1 g (∂X)
.

Figure 1 .

 1 Figure 1. (Left) The domain X (with boundary in blue) is endowed with a probability density pictured in grayscale representing the density of population in a city. The set Y (in red) represents the location of bakeries. Here, X, Y ⊆ R 2 and c(x, y) = |x -y| 2 (Middle) The Voronoi tessellation induced by the bakeries (Right) The Laguerre tessellation: the price of bread the bakery near the center of X is higher than at the other bakeries, effectively shrinking its Laguerre cell.

  y) + ψ(y))ρ(x)dH d g (x) -y∈Y ψ(y)ν y = y∈Y Lag y (ψ) (c(x, y) + ψ(y))ρ(x)dH d g (x) -y∈Y ψ(y)ν y (1.3)

Definition 1 . 1 (

 11 Loeper's condition). The cost c satisfies Loeper's condition if for every y in Y there exists a convex open subset Ω y of R d and a C 1,1 diffeomorphism exp c y : Ω y → Ω such that the functions p ∈ Ω y → c(exp c y p, y)c(exp c y p, z) (QC)

. 29 )

 29 Combining (3.25), (3.28) and (3.29) we obtain H d (B t ) const • |t -s|, which implies (3.21) using the boundedness of σ.

CRemark 4 . 2 .

 42 det := max y∈Y det(D exp c y ) Lip(Xy) ,where cond(A) is the condition number of a linear transform A between finite dimensional normed spaces and det(A) is the determinant of A with respect to orthornormal bases. The quantitative transversality estimates involve all the above constants in an explicit way, see(4.45). Even in the Euclidean case, one needs a lower bound on the volume of Laguerre cells in order to establish the second-differentiability of the functional Φ. Indeed, let y ± = ±1, let y 0 = 0, Y = {y -, y 0 , y + } ⊆ R.Consider the cost c(x, y) = -xy, and the density ρ = 1 on X = [-1 2 ,1 2 ]. Let ϕ τ ∈ R Y be defined by ϕ(y ± ) = 1/2 and ϕ(y 0 ) = τ . A simple calculation gives, for τ 0, ∂Φ ∂1 y 0 (ψ τ ) = max(1 -2τ, 0), which is not differentiable at τ = 1/2, even though (Reg), (Twist), and (QC) are all satisfied.

Figure 3 .

 3 Figure 3. Bound on the volume of a convex set K as a function of the angle between two normal vectors v, w at the same point and its diameter (see Lemma 4.3).

  ψ : Y → R we associate the vector λ ψ : = (ψ(y 1 )ψ(y 0 ), . . . , ψ(y N )ψ(y 0 )) ∈ R N .(4.39)

1 , 4 . 3 .

 143 and by recalling that V | W 0. Proof of Theorem 4.1. By Theorem 1.1, the second-differentiability of Kantorovich's functional Φ will follow from the differentiability of the function G y 0 (ψ) := Lag y 0 (ψ)

Definition 5 . 1 .Definition 5 . 2 .

 5152 Let (w yz ) (y,z)∈Y 2 be a weighted graph over Y . The (weighted) degree of a vertex y is d y := z =y w yz . The (weighted) Laplacian is the matrix (L yz ) (y,z)∈Y 2 whose entries are L yz = -w yz for y = z and L yy = d y . The Cheeger constant of a weighted graph (w yz ) (y,z)∈Y 2 over a point set Y is given by h(w) := min S⊆Y |∂S| w min(|S| w , |Y \ S| w ) , where |∂S| w := y∈S,z ∈S w yz and |S| w := y∈S d y .

5. 3 .

 3 Proof of Theorem 5.1. Let ψ be a function in ψ ∈ K ε and consider the weighted graph (w yz ) (y,z)∈Y 2 given by

Figure 4 .

 4 Figure 4. Evolution of the Laguerre cells during the execution of the damped Newton algorithm for semi-discrete optimal transport. (Top) The source density ρ is piecewise linear over the domain X = [0, 3] 3 over the displayed triangulation: it takes value 1 on the boundary of the square [0, 3] 2 and 0 on the boundary of [1, 2] 2 . The target measure is uniform over a 30 2 uniform grid in [0, 1] 2 . (Bottom) Laguerre cells at steps k = 0, 2, 6, 9, 12, 15, 18, 21 and 25.

S d- 1

 1 (r) f (z)f (r) dH d-1 (z) c d S d-1 (r) ∇f (z)| z ⊥ dH d-1 (z),where ∇f (z)| z ⊥ is the orthogonal projection of the gradient on the tangent plane {z} ⊥ , so thatR 0 ρ(r) r d-1 ω d-1 S d-1 (r) f (z)f (r) dH d-1 (z) c d B(0,r) ∇f (x)| x ⊥ ρ(x)dH d (x).

  ) f (r) -F dr c ρ R 0 f (r) ρ(r)dr. Now, notice that:f (r) = lim h→0 f (r + h)f (r) h = lim h→0 1 ω d-1 S d-1 f ((r + h)z)f (rz) h dH d-1 (z), dH d (x). (A.67)From the simple inequality (a + b) 2 2(a 2 + b 2 ), we get∇f (x) | x x + ∇f (x)| x ⊥ √ 2 ∇f (x) .Using the bounds (A.66) and (A.67) in Equation (A.65), we get the desired inequality:B(0,R) |f (x) -F | ρ(x)dH d (x) √ 2(c d + c ρ ) B(0,R) ∇f (x) ρ(x)dH d (x).Appendix B. Measure of singular points on convex setsLet K be a convex compact set of R d . The normal cone to K at a point x in K is the setN x K = {v ∈ R d | ∀y ∈ K, yx | v 0}, (B.68)and its elements are said to be normal to K at x. Let τ be a positive parameter. A point x in the boundary of K is τ -singular if there exist two unit vectors u, v in its normal cone N x K such that u | v 2 1τ 2 . Note in particular that if x is τ -singular, the linear space spanned by its normal cone has dimension two or more. Given a parameter τ > 0, we consider the set of τ -singular points Sing(K, τ ) := {x ∈ ∂K | ∃u, v ∈ N x (K) ∩ S d-1 , u | v 2 1τ 2 }. (B.69)

Proposition B. 1 .

 1 Let K be a convex, compact set of R d . Then H d-2 (Sing(K, τ )) const(d, diam(K)) 1 τ ,

  Cov(K, η) Pack(K, η) Cov(K, The proof consists in comparing a lower bound and an upper bound of the packing number of the setU := {(x, n) ∈ R d × S d-1 | x ∈ Sing(K, τ ) and n ∈ N x (K)}.

1 .Step 2 .Step 3 .

 123 Cov (∂K r , η) Cov(S, η) const(d) • (diam(K)/η) d-Combining these bounds with the inclusion U ⊆ U(K) gives usCov(U, ε) const(d, diam(K)) ε d-1 .(B.71) We now establish a lower bound for Pack(U, 2ε). Let x be a τsingular point and u, v be two unit vectors such that u | v 2 1τ 2 . This implies that N x K ∩ S d-1 contains a spherical geodesic segment of length at least const • τ , giving us a lower bound on the packing number ofN x K ∩ S d-1 , namely Pack(N x K ∩ S d-1 , η) const • τ /η. Now,let X be a maximal set in the definition of the packing number Pack(Sing(K, τ ), 2ε) and for every x ∈ X, let Y x be a maximal set in the definition of the packing number Pack(N x (K) ∩ S d-1 , 2ε), so that Card(Y x ) const • τ /ε. Then, the set Z := {(x, y) | x ∈ X, y ∈ Y x } is a 2ε packing of U ,and the cardinality of this set is bounded from below by const • Card(X) • τ /ε. This gives Pack(U, 2ε) const • Pack(Sing(K, τ ), 2ε) • τ /ε. (B.72) Combining Equations (B.71), (B.72) and the comparison between packing and covering numbers (B.70), we get Pack(Sing(K, τ ), 2ε) const(d, diam(K)) τ ε d-2 .

H d- 2 (

 2 Sing(K, τ )) lim inf ε→0 N ε ε d-2 const(d, diam(K)) 1 τ .

2 .

 2 |f i (p)f i (p)| 2 = |-c(exp c y 0 (p), y i ) + c(exp c y 0 (p), y i )| 2 (C ∇ C 2 exp ) 2 V -V w 2 = 2(C ∇ C 2 exp ) 2 V 2 (1v | w ) 2(C ∇ C 2 exp ) 2 (C exp C ∇ ) 2 (1v | w ).Combining with (C.73) we then have for any p ∈ (exp c y 0 ) -1 (Ω),f i (p)λ i max{0, f j (p)λ j } + √ 2C 3 exp C 2 ∇ 1v | w or re-arranging and using that λ < T tr , 1v | w sup p∈(exp c y 0 ) -1 (Ω) (f i (p) -max{0, f j (p)} -2T tr ) 2 2(C 3 exp C 2 ∇ ) (C.74)

  Proving the continuity of G amounts to proving the continuity of the functions G y (ψ) := ρ(Lag y (ψ)) for any y in Y . Fix y in Y and remark that by definition, Lag y (ψ) = z =y∈Y H z (ψ) where

	Moreover, using Lebesgue's monotone con-
	vergence theorem one easily sees that g (resp. h) is right-continuous (resp.
	left-continuous). This concludes the proof
	Proof of Proposition 2.3.

  [START_REF] Acosta | An optimal poincaré inequality in L 1 for convex domains[END_REF]. The goal here is to bound from below the discrete Cheeger constant h(w) in terms of the continuous weighted Cheeger constant h(ρ) and the constants introduced in (4.33). By definition of the constants ε tw and C ∇ , for every y = z in Y , one hasε tw w yz | Lag y,z (ψ)| ρ 2C ∇ w yz .(5.51)Consider a subset S of Y , and let A = ∪ y∈S Lag y (ψ). Then, the intersection of the boundary of A with X is contained an union of facets of Laguerre cells, namely ∂A ∩ int(X) ⊆

							Lag y,z (ψ).	(5.52)
						y∈S,z ∈S
	The two inequalities (5.51) and (5.52) imply a lower bound on the numerator
	of the Cheeger constant:
			|∂A| ρ	y∈S,z ∈S	Lag y,z (ψ) ρ 2C ∇ |∂S| w .	(5.53)
	We now need to bound the denominator of the Cheeger constant from above,
	which requires us to control the weighted degrees d y . Note that
	d y =	z =y	w yz	1 ε tw z =y	| Lag y,z (ψ)| ρ	1 ε tw	|∂ Lag y (ψ)| ρ ,	(5.54)

  |A| ρ = y∈S Lag y (ψ) ρ we obtain

	Combining (5.54), (5.55) and
	exp	1)	g (∂X) Lag y (ψ) ρ H d-1	(5.55)

https://github.com/mrgt/PyMongeAmpere

In this section, we provide an (almost) explicit example of a probability density on R d whose support is an annulus, therefore not simply connected, but which still satisfies a weighted Poincaré-Wirtinger inequality.

where ω d-1 is the volume of the unit sphere S d-1 . Then, ρ satisfies the weighted Poincaré-Wirtinger inequality (PW) for some positive constant.

Here we provide the proof of the alternative upper transversality estimate in Section 4.4.

Proof of Proposition 4.8. Let us again write

and assume ε nd < V W and v | w > 0. Let us also define

x 0 : = exp c y 0 (p 0 ), q 0 : = -D x c(x 0 , y 0 ), q 1 := -D x c(x 0 , y j ).

A quick calculation yields

W ) + q 0 . Now we define the point q := [(D exp c y 0 | p 0 ) * ] -1 ( V w + q 0 ), since V W , the above calculation yields that q lies on the line segment between q 0 and q 1 ; since (exp c

x 0 ) -1 (Λ) is convex we have that q ∈ (exp c

x 0 ) -1 (Λ) as well. Thus we can define

x 0 (q ), f i (p) : = -c(exp c y 0 (p), y i ) + c(exp c y 0 (p), y 0 ) + c(exp c y 0 (p 0 ), y i )c(exp c y 0 (p 0 ), y 0 ) + λ i , and by (4.48) we will obtain for all p ∈ (exp c y 0 ) -1 (Ω),

while another quick calculation yields

Translating this back into coordinates in (exp c 0 ) -1 (X) and in terms of f i , f j , we see there exists a point p c ∈ (exp c 0 ) -1 (X) for which

.

Thus if we have T tr

, combining with (C.74) we will obtain the bound (4.49) as desired.