Boris Muzellec
email: boris.muzellec@inria.fr

Adrien Vacher
email: adrien.vacher@u-pem.fr

Francis Bach
email: francis.bach@inria.fr

François-Xavier Vialard
email: francois-xavier.vialard@u-pem.fr

Alessandro Rudi
email: alessandro.rudi@inria.fr

Near-optimal estimation of smooth transport maps with kernel sums-of-squares

It was recently shown that under smoothness conditions, the squared Wasserstein distance between two distributions could be efficiently computed with appealing statistical error upper bounds. However, rather than the distance itself, the object of interest for applications such as generative modeling is the underlying optimal transport map. Hence, computational and statistical guarantees need to be obtained for the estimated maps themselves. In this paper, we propose the first tractable algorithm for which the statistical L 2 error on the maps nearly matches the existing minimax lower-bounds for smooth map estimation. Our method is based on solving the semi-dual formulation of optimal transport with an infinite-dimensional sum-of-squares reformulation, and leads to an algorithm which has dimension-free polynomial rates in the number of samples, with potentially exponentially dimension-dependent constants.

Introduction

Optimal transport (OT) provides a principled method to compare probability distributions, by finding the optimal way of coupling one to another based on a cost function on their supports. This optimization problem yields two useful quantities: the transport cost itself, which corresponds to the Wasserstein distance when the ground cost is a distance, and the minimizer, which is a map that pushes forward the first distribution onto the second, known as the transport map. While OT has gained attention lately in the statistics and machine learning communities due to the properties of the Wasserstein distance, transportation maps are playing an increasingly important role in data sciences. Indeed, many applications such as generative modeling [START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF][START_REF] Salimans | Improving GANs using optimal transport[END_REF][START_REF] Bernton | Inference in generative models using the wasserstein distance[END_REF][START_REF] Makkuva | Optimal transport mapping via input convex neural networks[END_REF][START_REF] Onken | Ot-flow: Fast and accurate continuous normalizing flows via optimal transport[END_REF], domain adaptation [START_REF] Courty | Optimal transport for domain adaptation[END_REF][START_REF] Courty | Joint distribution optimal transportation for domain adaptation[END_REF], shape matching [START_REF] Su | Optimal mass transport for shape matching and comparison[END_REF][START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF] or predicting cell trajectories [START_REF] Schiebinger | Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming[END_REF][START_REF] Dai | Predicting cell lineages using autoencoders and optimal transport[END_REF] among others can be formulated as the problem of finding a map from a reference distribution to a target distribution.

Over the past decades, efforts have particularly concentrated on the problem of computing OT distances. Two cases must be distinguished: the case of discrete measures (i.e., measures supported on a finite number of points), and the general case, including e.g. measures that admit a density with respect to the Lebesgue measure. For discrete measures supported on N points, OT distances and plans may be computed by solving a linear program (LP) in O(N 3 log(N )) time using the network simplex algorithm (see, e.g., [START_REF] Ravindra K Ahuja | Network flows: theory, algorithms, and applications[END_REF]. Using entropic optimal transport [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] as a proxy, this computational cost can be further reduced to O(N 2 ). However, efficient methods for general measures remain to be found. The naive approach, which consists in using the OT distance between samples of the distributions (the so-called plugin estimator), fails as the dimension grows. Indeed, the sharpest statistical bounds known for the plugin estimator require ε -2 d samples to achieve precision ε [START_REF] Chizat | Faster Wasserstein distance estimation with the Sinkhorn divergence[END_REF]). Yet, theoretical estimators were recently derived in the case where the problem is smooth [START_REF] Weed | Estimation of smooth densities in Wasserstein distance[END_REF][START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF], yielding rates of estimation of the OT distance that improve as the smoothness grows. However, the estimators presented in those works are not algorithmically tractable. Very recently, [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] closed this statistical-computational gap by designing an estimator of the squared Wasserstein distance relying on kernel sums-of-squares [START_REF] Marteau-Ferey | Non-parametric models for non-negative functions[END_REF], and in particular [START_REF] Rudi | Finding global minima via kernel approximations[END_REF]) that may be computed with polynomial dimension-free rates in the number of samples, but with a constant term that is potentially exponential in the dimension.

On the other hand, fewer results on computationally efficient statistical estimators of optimal maps are available in the literature. Compared to the problem of estimating OT distances, the estimation of an optimal transport map based on samples is particularly difficult since the estimated maps need to be evaluated on out-of-sample points. When smoothness assumptions are made, one can design theoretical statistical estimators of the transport maps that are statistically minimax optimal with respect to the smoothness for the L 2 error [START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF], but such estimators cannot be computed in practice as they require projecting on the space of smooth and strongly convex functions, which is NP-hard. In a few empirical studies, the potentials were explicitly parameterized either by a neural network [START_REF] Seguy | Large-scale optimal transport and mapping estimation[END_REF] or by a Gaussian kernel [START_REF] Genevay | Stochastic optimization for large-scale optimal transport[END_REF]. However, neither computational nor statistical guarantees are provided in those works. Likewise, [START_REF] Paty | Regularity as regularization: Smooth and strongly convex brenier potentials in optimal transport[END_REF] recently managed to derive an explicit algorithm when the potentials are assumed to have smooth gradients and to be strongly convex, but the corresponding rate of approximation is not known. A recent line of research studies estimators of optimal maps that rely on barycentric projection, which can be developed in regularized and non-regularized settings [START_REF] Gunsilius | On the convergence rate of potentials of Brenier maps[END_REF][START_REF] Deb | Rates of estimation of optimal transport maps using plug-in estimators via barycentric projections[END_REF][START_REF] Pooladian | Entropic estimation of optimal transport maps[END_REF]. These estimators can be used without any further assumption on the source and target measures. Here again, either the rates do match minimax rates of estimation of [START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF] without computationally feasible estimators, or the rates are not optimal. In particular, the estimator of [START_REF] Pooladian | Entropic estimation of optimal transport maps[END_REF] is based on entropic regularization and is computationally friendly, but only works for low regularity of the maps, and is not minimax optimal.

Contributions. The contributions of this paper are twofold. First, relying on the sumof-squares (SoS) tight reformulation of OT proposed in [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF], we bridge the statistical-computational gap on the estimation of the optimal potentials. Second, we propose improvements on the algorithms proposed by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] to enhance the practicality of SoS for optimal transport.

1. An estimator for smooth OT potentials. Relying on an SoS representation of Brenier dual constraints, we design an algorithm to compute an estimator of the OT potentials based on samples. Remarkably, the complexity of this algorithm has polynomial rates in the number of samples (with potentially exponential constants w.r.t. the underlying dimension). For this estimator, we derive a statistical upper bound on the L 2 error on the resulting maps that nearly matches the lower bound of [START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF] when the smoothness is high. To achieve these nearly optimal rates, we refine the regularization that was proposed by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] and carry out a careful analysis to bound the solutions of the empirical problem. The central argument consists in using the local strong convexity of the so-called semi-dual formulation of optimal transport. It enables us to relate the L 2 error both to the quality of the (scalar) approximation of the original problem and to the norm of the empirical minimizers.

In turn, this allows us to obtain a statistical rate for the L 2 error with well-chosen regularizers.

2. Algorithmic improvements. We propose several improvements on the algorithms introduced by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF]. First, as our objective is to compute an estimator of the OT potentials themselves, our algorithm should provide convergence guarantees for the minimizing solution, and not only the minimum. To do so, we reformulate the objective as a strongly convex optimization problem by using a square norm regularization instead of the trace, and by replacing hard constraints with quadratic penalties. A second advantage of this strongly convex relaxation is that it allows using first order optimization methods such as accelerated gradient descent, which are more scalable than the Newton methods employed by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] whose cost quickly blows up with the number of samples. Next, we further reduce computational costs by introducing a Nyström approximation of the features in the dual constraints. Finally, we propose a criterion to choose the hyperparameters on which our estimator relies, which was lacking in the work of Vacher et al.

Notations and background

Let (µ, ν) be two probability distributions on two bounded domains X, Y ⊂ R d . We study the squared Wasserstein distance, whose dual Kantorovitch formulation is given by

W 2 2 (µ, ν) = sup u,v∈C(R d ) u, µ + v, ν subject to x -y 2 2 ≥ u(x) + v(y), ∀(x, y) ∈ X × Y , (1)
where C(R d ) is the space of continuous functions on R d . For convenience, we shall use the equivalent Brenier formulation, where the Kantorovitch potentials u, v are replaced by the

Brenier potentials f (•) := • 2 2 -u(•), g(•) := • 2 2 -v(•) OT(µ, ν) = inf f,g∈C(R d ) f, µ + g, ν subject to f (x) + g(y) ≥ x, y , ∀(x, y) ∈ X × Y .
(2)

These two quantities are related by W 2 2 (µ, ν) = . 2 2 , µ + ν -OT(µ, ν). Problem (2) (as well as problem (1)) is delicate to solve numerically due to the nonnegativity constraint which has to be satisfied on a continuous set. For instance, a feasible idea consists in sampling the inequality constraint on a finite set and trying to extend it to the continuous set. Unfortunately, this strategy may only leverage Lipschiztness, even if the functions involved are smoother than Lipschitz, yielding an approximation rate of order n -1 2d , 2d being the dimension of X × Y . The method proposed by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] overcomes this difficulty and is able to leverage the smoothness of the potentials. Following their work, we require assumptions on the support and the smoothness of the densities themselves. Using Cafarelli's regularity theory [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF], Assumption 1 ensures that the Brenier potentials have a similar order of differentiability. In particular, defining the Sobolev space of order m over

Z ⊂ R d H m (Z) := f ∈ L 2 (Z) | f H m := |α|≤m D α f L 2 (Z) < ∞ , (3)
the optimal Brenier potentials (f * , g * ) belong to the Sobolev space H m+2 (X) and H m+2 (Y ) respectively [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF]. When m > d/2 -2, the Sobolev embedding theorem gives that these spaces are continuously embedded in the space of continuous functions and thus are reproducing kernel Hilbert spaces (RKHS, see e.g. [START_REF] Paulsen | An Introduction to the Theory of Reproducing Kernel Hilbert Spaces[END_REF][START_REF] Steinwart | Support Vector Machines[END_REF]. The approach proposed in Vacher et al. to estimate the Wasserstein distance is based on the kernel Sum-Of-Squares (SoS) approach and the tools introduced by [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] to deal with optimization problems subject to a dense set of inequalities. The procedure consists in two steps: (1) show that the optimization problem is equivalent to a problem where the inequality constraint is substituted by an equality constraint with respect to an SoS term; (2) consider an empirical version of the resulting problem with only a finite number of equality constraints. [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] shows, for the case of non-convex optimization, that this procedure is adaptive to the degree of differentiability of the objective function leading to rates that overcome the curse of dimensionality for very smooth objectives. [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] shows that for the case of the problem in (2), the two steps correspond to the following Eq. (4) and Eq. ( 5) as reported in the following theorem and below.

Theorem 1 [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF]). Let µ ∈ P(X), ν ∈ P(Y ) satisfy Assumption 1 and assume that m > d + 1. Then, problem (2) can be reformulated as

OT(µ, ν) = inf f,g∈C(R d ) A∈S + (H m (X×Y )) f, µ + g, ν subject to f (x) + g(y) = x, y + φ(x, y), Aφ(x, y) H m (X×Y )) , (4)
where

S + (H m (X × Y )) is the set of linear, positive, self-adjoint operators on H m (X × Y ).
The key contribution of this representation theorem is to replace the inequality constraint by an equality constraint which is easier to deal with, as proposed and analyzed in [START_REF] Rudi | Finding global minima via kernel approximations[END_REF]. Given access to samples

x 1 , • • • , x n ∼ µ and y 1 , • • • , y n ∼ ν with associated empirical distributions μ and ν, we solve for λ 1 , λ 2 ≥ 0 the empirical problem min f ∈H m+2 (X),g∈H m+2 (Y ), A∈S + (H m (X×Y )) f, μ + g, ν + λ 1 Tr(A) + λ 2 ( f 2 H m+2 (X) + g 2 H m+2 (Y ) ) subject to ∀j ∈ [n], f (x j ) + g(y j ) -x j , y j = φ(x j , y j ), Aφ(x j , y j ) H m (X×Y ) , (5)
where φ :

X × Y → H m (X × Y ) is the feature map of H m (X × Y ).
Using the techniques of Marteau-Ferey et al. ( 2020) and [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], the authors show that problem (5) that in the case where m > 2d, the energy of the empirical potentials can be controlled with high probability by setting

λ 1 = λ 2 ∼ 1 √ n
, leading to a dimension-free approximation of OT(µ, ν) at a 1 √ n rate with high probability, with a computational complexity of O(n 3.5 log( 1 ε )), for a precision of ε.

As problem ( 5) is strongly convex with respect to the potentials f, g, the uniqueness of the empirical potentials fn , ĝn is ensured. Further, their existence may also be proven. However, the question remains of recovering a rate of convergence of the empirical potentials toward the optimal Brenier potentials (f * , g * ) with respect to some norm. [START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF] prove that when d > 3, any estimator of the transport map Tn can achieve at best an L 2 error that scales as

E( ∇f * -Tn 2 L 2 (µ) ) ∼ n -m+1 m+d/2 . ( 6
)
We prove in the next section that the solutions ( fn , ĝn ) of the problem (5) can nearly match this rate up to an exponent ε, where ε can be chosen arbitrary close to 0, for well-chosen regularization parameters λ 1 , λ 2 .

3 Nearly optimal minimax rates

The key to derive the nearly minimax optimal rates is to have estimates on λ 1 , λ 2 that are sharper than 1 √ n . Indeed, in the first step of the proof we show an upper-bound of the form

(λ 1 +λ 2 )( fn H m+2 (X) + ĝn H m+2 (Y ) + Ân ) 1 √ n ( ∇ fn -∇f * 2 L 2 (µ) + ∇ĝ n -∇g * 2 L 2 (ν) ) α ,
where α is positive, decreases with the dimension d to zero and increases with the smoothness m and ( fn , ĝn , Ân ) are minimizers of problem (5). In order to control in turn the errors

∇f -∇f * 2 L 2 (µ) , ∇g -∇g * 2 L 2 (ν)
with the regularizers λ 1 , λ 2 , we need extra convexity. To this end, we introduce the semi-dual formulation of optimal transport [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF] which replaces the potential pair (f, g) by the couple

(f, f * ) where f * is the Fenchel-Legendre transform f * (y) = sup x∈X x y -f (x). It reads OT(µ, ν) = inf f ∈C(R d ) J µ,ν (f ), (7)
where J µ,ν (f ) = f, µ + f * , ν . We simply denote it by J(f ) when no confusion is possible.

As shown in the next lemma, the functional J(f ) gains stronger convexity in comparison with the linear objective of Brenier's formulation.

Lemma 2. The semi-dual objective

J(f ) is convex in f ∈ C(X).
Assuming that there exists an optimal potential f * such that ∇f * pushes µ onto ν and that f is a convex C 1 function with M -Lipschitz gradient, we have

∇f -T * 2 L 2 (µ) ≤ 2M (J(f ) -J(f * )) . ( 8
)
Note that there is no assumption on the measures µ and ν themselves. However, we assume the Lipschitz smoothness of the gradient of f and the existence of an optimal map, which is ensured by Brenier's theorem if µ has density w.r.t. the Lebesgue measure. A similar result is proven by [START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF] and we include for completeness a simple proof in Apppendix A. Thanks to this gain of convexity, we obtain an upper bound of the errors by the regularizers λ 1 , λ 2 of the form

∇ fn -∇f * 2 L 2 (µ) + ∇ĝ n -∇g * 2 L 2 (ν) (λ 1 + λ 2 ) 1 √ n . (9)
Thanks to these two connected bounds, we can calibrate the regularizers λ 1 , λ 2 and eventually obtain rates sharper than 1 √ n .

Theorem 3. Let δ, ε ∈]0, 1[ 2 and assume that the regularizers are given by

λ 1 n = λ 2 n = λ n = log( 2 δ ) 2 n m+1 m+d/2+ε + C 1 log( n δ ) n m-d 2d , (10)
where C 1 is a constant that does not depend on n and δ. Denoting fn (resp. ĝn ) the transport potential associated to µ (resp. ν) in problem (5), we have with probability at least

1 -δ that for n ≥ n 0 (X, Y, d, m), ∇ fn -∇f * 2 L 2 (µ) + ∇ĝ n -∇g * 2 L 2 (ν) ≤ C 2 λ n , (11)
where C 2 is a positive constant that is independent from n and δ but grows to infinity when ε goes to 0. In particular, when m is sufficiently large, the minimax rate is nearly attained:

∇ fn -∇f * 2 L 2 (µ) + ∇ĝ n -∇g * 2 L 2 (ν) ≤ C 2 log( 2 δ ) n m+1 m+d/2+ε . ( 12
)
Proof Let us denote by ( fn , ĝn , Ân ) the solutions of the empirical problem (5) and by (f * , g * , A * ) the solutions of the deterministic problem (4). We define the energy of the potentials

R 2 = f * 2 H m+2 (X) + g * 2 H m+2 (Y ) and its empirical counterpart R2 n = fn 2 H m+2 (X) + ĝn 2 H m+2 (Y )
. For probability measures (α, β), we shall denote throughout the proof the linear objective of the dual as E α,β (f, g) = f, α + g, β . The strategy of the proof is as follows 1. We use the optimality conditions of the empirical and deterministic OT problems to upper bound λ n and the norms of the empirical objects ( fn , ĝn , Ân ) by the errors

∇ fn -f * L 2 (µ) , ∇ĝ n -g * L 2 (ν) .
2. We use the strong convexity of the semi-dual to upper bound in turn the errors by the regularizer λ n and the norms of ( fn , ĝn , Ân ).

3. We obtain two connected bounds and we optimize over λ n to obtain sharp rates while maintaining the empirical objects ( fn , ĝn , Ân ) bounded.

Upper bound on the regularizer. The optimality conditions in the empirical problem (5) yield

λ n (Tr( Ân ) -Tr(A * ) + R2 n -R 2 ) ≤ E μ,ν (f * , g * ) -E μ,ν ( fn , ĝn ) . (13)
Conversely, we wish to test the empirical potentials against the measures µ, ν and use the optimality condition of the deterministic problem.

Proposition 4. There exist constants C, G such that for all 0 < δ < 1, defining κn,δ = C( Rn + Tr( Ân )+G)n -m-d 2d log( n δ )

m-d 2d we have with probability at least 1-δ that for n ≥ n 0 (X, Y, d, m), the empirical potentials ( fn , ĝn + κn,δ ) are admissible candidates:

∀(x, y) ∈ X × Y, fn (x) + ĝn (y) + κn,δ ≥ x y . ( 14
)
The proof is left in Appendix A. We deduce that with probability at least 1 -δ

0 ≤ E µ,ν ( fn , ĝn ) -E µ,ν (f * , g * ) + κn,δ . (15)
Adding the equation above with the empirical optimality conditions (13) yields

λ n (Tr( Ân ) -Tr(A * ) + R2 n -R 2 ) ≤ E μ,ν (f * , g * ) -E μ,ν ( fn , ĝn ) + E µ,ν ( fn , ĝn ) -E µ,ν (f * , g * ) + κn,δ = f * -fn , μ + g * -ĝn , ν + fn -f * , µ + ĝn -g * , ν + κn,δ = fn -f * , µ -μ + ĝn -g * , ν -ν + κn,δ
Now we use the following lemma upper-bounding the linear forms of the r.h.s. with respect to the errors

∇ fn -f * L 2 (µ) , ∇ĝ n -g * L 2 (ν) . Lemma 5. Let (u, v) ∈ H m+2 (X) × H m+2 (X), (δ, ε) ∈]0, 1[×]0, 1[
and µ be a probability measure over X. If µ has a density on X bounded away from zero, then there exists a constant C independent of n, δ, u, v such that with probability at least 1 -δ, we have

u -v, µ -μ ≤ C µ log( 2 δ ) √ n ∆ m+2-d/2-ε m+1 u,v,µ R(u, v) d/2+ε-1 m+1 + ∆ u,v,µ , (16)
where

∆ u,v,µ = ∇u -∇v L 2 (µ) and R(u, v) = u H m+2 (X) + v H m+2 (X) .
The proof of this lemma is left in the Appendix A and is mainly an application of the Gagliardo-Nirenberg inequality and the Poincaré inequality. Applying the lemma to (f * , fn , µ) and (g * , ĝn , ν), and using the fact that for (α, x, y) ∈ (R + ) 3 , we have x α + y α ≤ 2(x + y) α , we recover that for n ≥ n 0 (X, Y, d, m) we have with probability at least 1 -3δ

λ n (Tr( Ân ) -Tr(A * ) + R2 n -R 2 ) ≤ 2C log( 2 δ ) √ n ( Rn + R) d/2+ε-1 m+1 a m+2-d/2-ε m+1 n + 2C log( 2 δ ) √ n a n + κn,δ , (17)
where

C = C µ + C ν and we denoted a n = ∇ fn -∇f * L 2 (µ) + ∇ĝ n -∇g * L 2 (ν) .
Bounding the error a n . First we bound the errors from above via the strong convexity of the semi-dual. Then, using the optimality conditions, we recover the same linear upper bound as in the previous paragraph.

a) Semi-dual upper bound. A straightforward way to upper-bound the errors would be to apply Lemma 2 to the empirical potentials fn , ĝn . Unfortunately, there is no guarantee that the empirical potentials are convex. Instead, we define for t ∈ [0, 1] the following interpolating potential fn (t)(.) = f * (.) + t( fn (.) -f * (.)) .

Since f * and g * are mutual Legendre transforms and that they both have Lipschitz gradients, (f * , g * ) are strongly convex and we denote by γ a strong convexity constant of the optimal Brenier potentials. Hence for t ≤

γ 2 fn-f W 2,∞ (X)
, the interpolate fn (t)(.) is γ 2 strongly convex. Hence we choose tf = min 1, γ

2 fn -f W 2,∞ (X) , (19)
and we apply Lemma 2 to fn ( tf )(.), which yields

∇ fn ( tf ) -∇f * 2 L 2 (µ) ≤ 2 Lf (J( fn ( tf )) -J(f * )) , (20)
where Lf = fn ( tf ) W 2,∞ (X) . By convexity of the semi-dual J, the r.h.s. is upper bounded by 2 Lf tf (J( fn ) -J(f * )) and the l.h.s. is equal to t2 f ∇ fn -∇f 2 L 2 (µ) . Hence, we have

∇ fn -∇f 2 L 2 (µ) ≤ 2 Lf tf (J( fn ) -J(f * )) . (21)
Finally, we upper-bound the quantity Lf tf . Proposition 6. Denoting K X the embedding constant such that ∀u

∈ W 2,∞ (X), u W 2,∞ (X) ≤ K X u H m+2 (X)
, the quantity Lf tf can be upper bounded as

Lf tf ≤ 2K X ( fn H m+2 (X) + f * H m+2 (X) ) 1 + f * H m+2 (X) γ . ( 22
)
The proof is left in Appendix A. Conversely, applying this result to ĝn , we obtain

Lg tg ≤ 2K Y ( ĝn H m+2 (Y ) + g * H m+2 (Y ) ) 1 + g * H m+2 (Y ) γ , ( 23
)
where

K Y is the embedding constant from H m+2 (Y ) to W 2,∞ (Y ). In particular Lf tf + Lg tg ≤ 2(K X + K Y )( fn H m+2 (X) + ĝn H m+2 (Y ) + f * H m+2 (X) + g * H m+2 (Y ) ) × 1 + f * H m+2 (X) ) + g * H m+2 (Y ) γ (24) ≤ 2 √ 2(K X + K Y )( Rn + R) 1 + √ 2R γ . (25)
b) Linear upper bound. Recall that the couple ( fn , ĝn + κn,δ ) are admissible candidates and in particular J( fn ) ≤ E µ,ν ( fn , ĝn + κn,δ ), which implies

J( fn ) -J(f * ) ≤ E µ,ν ( fn , ĝn ) -E µ,ν (f * , g * ) + κn,δ . (26)
Combining the equation above with the empirical optimality condition 0 ≤ E μ,ν (f * , g * ) -E μ,ν ( fn , ĝn ) + λ n (Tr(A * ) + R 2 ) (13), we obtain

J( fn ) -J(f * ) ≤ f * -fn , µ -μ + g * -ĝn , ν -ν + λ n (Tr(A * ) + R 2 ) + κn,δ . (27)
As in the previous paragraph, we use Lemma 5 to upper-bound the linear forms of the r.h.s., and denoting a n = ∇ fn -∇f * L 2 (µ) + ∇ĝ n -∇g * L 2 (ν) , we obtain

J( fn ) -J(f * ) ≤2C log( 2 δ ) √ n ( Rn + R) d/2+ε-1 m+1 a m+2-d/2-ε m+1 n + a n + λ n (Tr(A * ) + R 2 ) + κn,δ , (28)
where C is a constant independent of n, δ.

Combining the bounds on the error a n . From the previous paragraphs, we obtain the following bounds

λ n (Tr( Ân ) + R2 n ) ≤ 2C log( 2 δ ) √ n [( Rn + R) d/2+ε-1 m+1 a m+2-d/2-ε m+1 n + a n ] + λ n (Tr(A * ) + R 2 ) + κn,δ , (29)
a 2 n ≤ 2( Rn + R)KC log( 2 δ ) √ n ( Rn + R) d/2+ε-1 m+1 a m+2-d/2-ε m+1 n + a n + λ n (Tr(A * ) + R 2 ) + κn,δ , (30)
where we defined

K = √ 2(K X + K Y ) 1 + 2 √ 2R γ
. The inequalities above mutually relate the empirical norms Rn , Tr( Ân ) and the error a n . The next result shows that for a well-chosen regularizer λ n , the empirical norms are bounded independently on δ and that for this choice of regularizer, the error a 2 n is at most of the order λ n .

Proposition 7. Assume that Equations ( 29) and (30) hold for n ≥ n 0 . If we set the regularizer as

λ n = log( 2 δ ) 2 n m+1 m+d/2+ε + C 1 log( n δ ) n m-d 2d , ( 31
)
where C 1 is the constant independent of n, δ, then Rn , Tr( Ân ) are bounded independently of δ and there exists a constant C 2 independent of n, δ, such that

a 2 n ≤ C 2 λ n . ( 32
)
The proof is presented in Appendix A. Hence setting λ n to the indicated value yields that with probability at least 1 -3δ, for n ≥ n 0 (X, Y, d, m) we have

a 2 n ≤ Cλ n , (33)
with C a constant independent of n, δ.

Discussion and remarks. Let us make a few remarks regarding the comparison of our results to those of [START_REF] Hütter | Minimax estimation of smooth optimal transport maps[END_REF] and [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF].

Remark 8. Our bound does not exactly match the lower-bound n -m+1 m+d/2 of Hütter and Rigollet ( 2021) because of the extra exponent ε. Indeed, we cannot take ε = 0 since the constant C 2 would then go to infinity. However, it would be interesting to recover the growth rate of C 2 when ε goes to 0 ; if it does not increase too fast, we could hope for a trade-off in n on ε and maybe recover the exact lower-bound.

Remark 9. Note that the minimax rate is nearly attained only for the highly smooth case m > 2d, and that when m goes toward d we cannot even guarantee that empirical potentials converges to the original potentials. This gap compared with Hütter and Rigollet (2021) is explained by the fact that in their case, the authors directly solve the semi-dual problem while in our case, we approximate the cost constraint through the SoS reformulation which can only be achieved when the smoothness is sufficiently high.

Remark 10. Comparing ours results to [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF], one may notice that we require lower values for the regularizers. In particular, one may wonder if this could lead to gains in terms of the estimation of the squared Wasserstein distance itself. This is not the case, as this effect is due to the fact that we are here considering the Brenier version of optimal transport. Indeed, it holds W 2 2 (µ, ν) = . 2 2 , µ + ν -OT(µ, ν). Hence, to estimate W 2 2 (µ, ν) one must estimate the extra moment term, that fluctuates as1 √ n .

Algorithms and Numerical Experiments

In this section, we provide and test algorithms to solve Eq. ( 5) and compute the estimators f and ĝ. Following [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] and [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF], we may rewrite Eq. ( 5) using a finite-dimensional representation of A that relies on finite-dimensional features

Φ i ∈ R , i ∈ [n], and a PSD matrix B ∈ S + (R ): min f ∈H m+2 (X),g∈H m+2 (Y ), B∈S + (R ) f, μ + g, ν + λ 1 TrB + λ 2 ( f 2 H m+2 (X) + g 2 H m+2 (Y ) ) subject to ∀j ∈ [ ], f (x j ) + g(ỹ j ) -xj , ỹj = Φ T j BΦ j . ( 34
)
We refer to [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] for the definition of Φ j , j ∈ [ ] and the derivation of Eq. ( 34). Note that compared to Eq. ( 5) we have two sets of samples that do not necessarily coincide:

x 1 . . . x n ∈ X and y 1 , . . . , y n ∈ Y which are the support points of the empirical distributions μ and ν, and (x 1 , ỹi ) . . . (x , ỹ ) ∈ X × Y which are used to subsample the dual constraints.

One can then apply the method proposed by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] and solve ( 5) by using barrier methods on its dual formulation:

min γ∈R 1 4λ 2 γ T Qγ - 1 2λ 2 j=1 γ j z j + q 2 2λ 2 - δ log det   j=1 γ j Φ j Φ T j + λ 1 I   subject to j=1 γ j Φ j Φ T j + λ 1 I 0, (35)
where

Q ∈ R × is defined as Q i,j = k(x i , ỹi )+k(x i , ỹi ), i, j ∈ [ ], z j = xj • ỹj +w μ(x j )+w ν (ỹ j ), q 2 = w μ(x j ) 2 H + w ν (ỹ j ) 2
H , and δ > 0 is the barrier parameter. 1 The authors showed that solving this problem to precision ε has a O(n 3.5 log 1 ε ) computational cost. This method has two drawbacks: first, its cost is prohibitive when the number of samples becomes large, which is necessary to ensure better statistical approximation, and second, there are no guarantees on the convergence to the minimizers f and ĝ as only the objective value is guaranteed to converge to the optimal value. Both issues can be alleviated by introducing a strongly convex and unconstrained relaxation of problem (5).

Let us now consider a variant of problem (5) with a squared norm penalty (instead of a trace penalty) on B:

inf f ∈H X ,g∈H Y , B∈S + (R ) f, ŵµ H X + g, ŵν H Y + λ 1 2 B 2 F + λ 2 ( f 2 H X + g 2 H Y ) subject to ∀j ∈ [ ], f (x j ) + g(ỹ j ) -xj , ỹj = Φ T j BΦ j .
(36)

Lemma 11. Equation (36) admits the following dual problem:

min γ∈R 1 4λ 2 γ T Qγ - 1 2λ 2 j=1 γ j z j + 1 2λ 1 (- i=1 γ i Φ i Φ T i ) + 2 F + q 2 4λ 2 , ( 37
)
where (M) + denotes the positive part of M ∈ R , i.e. its projection on S + (R ).

Proof in Appendix B, Page 26. Although the bounds described in Section 3 and by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] rely on a trace control of B, we conjecture that similar results could be derived using a square Hilbert-Schmidt norm control.

Faster algorithms with strong convexity. Replacing the trace penalty with a square norm penalty led to the unconstrained convex dual problem (37), without adding a barrier term. However, this problem is not necessarily strongly convex. Indeed, Q may have arbitrarily small eigenvalues, and the positive part term vanishes when its argument is a negative matrix. We now propose a strongly convex relaxation of (37), obtained by replacing the constraints in (36) with a quadratic penalization. That is, for δ > 0, we consider:

inf f ∈H X ,g∈H Y , B∈S + (R d ) f, ŵµ H X + g, ŵν H Y + λ 1 2 B 2 F + λ 2 ( f 2 H X + g 2 H Y ) + δ 2 j=1 (f (x j ) + g(ỹ j ) -xj , ỹj -Φ T j BΦ j ) 2 . ( 38
)
Let us derive the corresponding dual problem.

Proposition 12. Strong duality holds, and problem (38) admits the following dual formulation:

inf γ∈R 1 4λ 2 γ T Qγ - 1 2λ 2 j=1 γ j z j + 1 2λ 1 (- i=1 γ i Φ i Φ T i ) + 2 F + 2δ γ 2 + q 2 4λ 2 , ( 39
)
with the following primal-dual relations

f = 1 2λ 2 ( i=1 γ i φ X (x i ) -ŵµ ) g = 1 2λ 2 ( i=1 γ i φ Y (ỹ i ) -ŵν ). ( 40
)
This problem is strongly convex, with constant α = δ + 1 2λ 2 λ min (Q), and smooth, with constant

L = δ + 1 2λ 2 λ max (Q) + 1 λ 1 λ max (K • K) (where K = Φ T Φ).
(39) is a tight relaxation of (37), as the dual of the constrained problem may be recovered by taking δ → ∞. As mentioned above, for most kernels the spectrum of Q decays rapidly to zero. Therefore, problem (37) may not be considered strongly convex for practical purposes. On the other hand, problem (39) has a strong convexity constant that is bounded from under by δ , thanks to the quadratic regularization term on γ that appears in Eq. ( 39). Hence, we may solve (39) using accelerated gradient descent algorithms. The bottleneck of the computation is the positive part, which is obtained by computing the SVD of a × matrix in O( 3 ) time. The size of this matrix can be reduced by considering a low-rank approximation of K.

Note however that relaxing the constraints in (36) with a quadratic penalization introduces a gap with the theory developed in Section 3. However, we expect that bounds could be obtained in the case of quadratic penalties, much like soft constraints compared with hard constraints in kernel ridge regression. We leave this for future work.

Nyström approximation. The computational bottleneck of solving ( 39) is forming the matrix i=1 γ i Φ i Φ T i and computing the SVD at each each iteration, both in O( 3) time. This computational cost can be reduced using a Nyström approximation (see e.g. [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF] of the kernel matrix K. Indeed, such an approximation is of the form

K nys = VW -1 V T , (41)
where W ∈ R r×r is the submatrix corresponding to r < randomly sampled columns (and corresponding rows) of K, and V ∈ R ×r . Hence, writing W = L T L with L ∈ R r×r , we may consider features

Φ nys i ∈ R r , i ∈ [ ] corresponding to the columns of R nys def = L ¯T V T , (42)
verifying (R nys ) T R nys = K nys . Replacing K with K nys in the derivation of the dual leads to 39)), and computing its SVD costs O(r 3 ) (compared to O( 3)). This brings down the cost of a gradient step to O(r 3 + r 2 ), compared to O( 3 ) without using approximations.

inf γ∈R 1 4λ 2 γ T Qγ - 1 2λ 2 j=1 γ j z j + 1 2λ 1 (- j=1 γ j Φ (nys) j (Φ nys j ) T ) + 2 F + 2δ γ 2 + q 2 4λ 2 . ( 43
) Forming the matrix i=1 γ i Φ (nys) i (Φ nys i ) T costs O(r 2 ) (compared to O( 3 ) in (
The empirical results in Figures 1 to 3 tend to show that one can pick quite small ranks while retaining good performance: indeed, in those experiments the maximum employed rank is 100, whereas we increase the number of samples up to 1000. This is coherent with results on the Nyström approximation applied to statistical learning [START_REF] Rudi | Less is more: Nyström computational regularization[END_REF]. However, the effect of Nyström approximation is currently not measured in our statistical analysis.

Selecting the hyperparameters. Contrary to the optimization setting [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], there is no direct way of assessing the quality of the estimator ÔT (e.g. of the "smaller is better" type) for a given choice of hyperparameters, as the algorithm may output values that are larger or smaller than the true value of OT. Alternatively, we may use the output potentials to estimate transportation maps T1 , T2 from µ to ν and ν to µ, and compare the mapped points to the target distributions. That is, we compute

T1 : x → ∇f (x) = 1 2λ 2   j=1 γ j ∇ x k(x j , x) -∇ x ŵµ (x)   T2 : x → ∇f (x) = 1 2λ 2   j=1 γ j ∇ x k(ỹ j , x) -∇ x ŵν (x)   , (44)
and select the hyperparameters for which

MMD def = MMD 1 n µ i=1 δ T1 (x i ) , ν + MMD 1 n ν i=1 δ T2 (y i ) , μ (45)
is the smallest, where MMD is the maximum mean discrepancy [START_REF] Gretton | A kernel two-sample test[END_REF] for the kernel k, defined as

MMD(µ, ν) def = (k(x, x) -2k(x, y) + k(y, y))d(µ ⊗ ν)(x, y). ( 46
)
The choice of MMD as a criterion has two main motivations: first, it is a divergence that can be efficiently computed and whose plugin estimator converges at rate O(n -1/2 ), and second, as observed by [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] in Appendix F, when the filing samples xj , ỹj , j ∈ [ ] correspond to all pairs of samples in μ and ν, the dual problem (35) may be reformulated as a regularized optimal transport problem with MMD marginal penalties.

Numerical results. We now report empirical evaluation of our estimators. In Figures 1 to 3, we sample data from two Gaussian distributions (whose covariances follow a Wishart distribution), and solve problem (39) for hyperparameters pairs (λ 1 , λ 2 ) on the {10 -7 , 10 -6 , 10 -5 , 10 -4 , 10 -3 , 10 -2 } 2 grid. As explained above, we then compute the forward and backward transportation maps (44) corresponding to each (λ 1 , λ 2 ) pair, and report the performance of the pair minimizing (45), varying the number of samples from 25 to 1000, averaging over 20 random draws. For each number of samples, we report the error of the OT distance estimator (like Vacher et al. ( 2021)) compared to the plugin estimator, the mean squared error of T1 and T2 compared to the ground truths T 1 and T 2 (which admit an analytical expression (see [START_REF] Takatsu | Wasserstein geometry of Gaussian measures[END_REF]) evaluated on the samples from µ ⊗ ν:

MSE = 1 n n i=1 | T1 (x i ) -T 1 (x i )| 2 + 1 n n i=1 | T2 (y i ) -T 2 (y i )| 2 , ( 47
)
the best hyperparameter pair and the MMD metric (45). We may obtain several takeaways from Figures 1 to 3. First, we see that gains of using the kernel SoS approach to OT increase with the dimension: while SoS-OT matches the performance of the sampled estimator in 2D (Fig. 1), it becomes increasingly more performant than the plugin estimator in 4D (Fig. 2) and 8D (Fig. 3). Second, we observe that the map estimator converges to the true map, and that the decrease of the L 2 error ( 47) is well correlated with the MMD criterion (45). Finally, as expected, the optimal regularization parameters (selected from a gridsearch) decrease as the number of sample increases. Filling pairs (x i , ỹi ) are drawn from µ ⊗ ν. The number of filling pairs is equal to the number of µ and ν samples, reported on the x-axis. Note that the std in the top left is actually decreasing, but that the size of the shaded area visually increases due to logarithmic scaling.

A Additional Proofs

A.1 Proof of lemma 2

Proof First note that the Legendre transform is pointwise convex. For (f, g) and λ ∈ [0, 1],

we have for all y

(λf + (1 -λ)g) * (y) = sup x x y -(λf (x) + (1 -λ)g(x)) (48) ≤ λf * (y) + (1 -λ)g * (y) . (49)
We can apply the Pinelis inequality (see Proposition 2 of [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF] to control the second term. We have with probability at least 1 -δ

w µ -w μ H d/2+ε (X) ≤ k d,ε log( 2 δ ) √ n , (64)
where k d,ε is a constant that goes to infinity as 1/ε when ε → 0. Now let us deal with the first term. There exists a constant C such that

u -v -r H d/2+ε (X) ≤ C( ∇u -∇v H d/2+ε-1 (X) + u -v -r L 2 (X) ) . (65)
For the second term, Poincaré-Wirtinger inequality [START_REF] Meyers | Integral inequalities of Poincaré and Wirtinger type for BV functions[END_REF] states that there exists a constant C such that u -v -r L 2 (X) ≤ C ∇u -∇v L 2 (X) . For the first term, we apply the Gagliardo-Nirenberg inequality which yields

∇u -∇v H d/2+ε-1 (X) ≤ C 1 ∇u -∇v d/2+ε-1 m+1 H m+1 (X) ∇u -∇v 1- d/2+ε-1 m+1 L 2 (X) + C 2 ∇u -∇v L 2 (X) , (66)
for C 1 , C 2 constants independent of u, v and δ. Conversely, there exists C such that

∇u -∇v H m+1 (X) ≤ C u -v H m+2 (X) . (67)
Finally, denoting

µ 0 = inf x∈X dµ(x) dλ(x) , we have . L 2 (X) ≤ . L 2 (µ) µ 0
and we obtain with probability

1 -δ u -v, µ -μ ≤ K 0 log( 2 δ ) √ n K 1 u -v d/2+ε-1 m+1 H m+2 (X) ∇u -∇v m+2-d/2-ε m+1 L 2 (µ) + K 2 ∇u -∇v L 2 (µ) ,
where

K 0 = k d,ε C, K 1 = (C ) d/2+ε-1 m+1 µ m+2-d/2-ε m+1 0 and K 2 = C(C +C 2 ) µ 0 . Denoting C µ = K 0 (K 1 + K 2 ), we get u -v, µ -μ ≤ C µ log( 2 δ ) √ n u -v d/2+ε-1 m+1 H m+2 (X) ∇u -∇v m+2-d/2-ε m+1 L 2 (µ) + ∇u -∇v L 2 (µ) . (68)
A.4 Proof of Proposition 6

Proof Recall that Lf = fn ( tf ) W 2,∞ (X) and that fn ( tf ) = f * + tf ( fn -f * ). Hence Lf can be upper-bounded as

Lf ≤ f * W 2,∞ (X) + tf fn -f * W 2,∞ (X) , (69)
and as a consequence Lf tf ≤ f * W 2,∞ (X) tf + fn W 2,∞ (X) + f * W 2,∞ (X) . Now recall that tf is given by tf = min 1, γ

2 fn -f * W 2,∞ (X) , (70)
and in particular

1 tf ≤ 2 fn -f * W 2,∞ (X) γ + 1 ≤ 2 fn W 2,∞ (X) + f * W 2,∞ (X) γ + 1 , (71)
which yields

Lf tf ≤ f * W 2,∞ (X) 2 f * W 2,∞ (X) γ + 2 + fn W 2,∞ (X) 2 f * W 2,∞ (X) γ + 1 (72) ≤ 2( f * W 2,∞ (X) + fn W 2,∞ (X) ) f * W 2,∞ (X) γ + 1 (73) ≤ 2K X ( f * H m+2 (X) + fn H m+2 (X) ) f * W 2,∞ (X) γ + 1 , (74)
where K X is the embedding constant of W 2,∞ (X) into H m+2 (X) [START_REF] Adams | Sobolev Spaces[END_REF].

A.5 Proof of Proposition 7

Proof The objective of the proof is to find a parameter λ n as small as possible that keeps the quantities Rn and 93) and Eq. ( 94) into Eq. ( 91), we obtain Eq. (37).

B.2 Proof of Proposition 12

Proof Dual formulation. It holds

inf f ∈H X ,g∈H Y , B∈S + (R d ) f, ŵµ H X + g, ŵν H Y + λ 1 2 B 2 F + λ 2 ( f 2 H X + g 2 H Y ) + δ 2 j=1 (f (x j ) + g(ỹ j ) -xj • ỹj -Φ T j BΦ j ) 2 , = sup γ∈R inf f ∈H X ,g∈H Y , B∈S + (R d ) f, ŵµ H X + g, ŵν H Y + λ 1 2 B 2 F + λ 2 ( f 2 H X + g 2 H Y ) + j=1 γ j (Φ T j BΦ j -f (x j ) -g(ỹ j ) + xj • ỹj ) - 2δ γ 2 .
(95) Indeed, the second problem above is a convex-concave min-max problem, and inverting min and max and solving for γ directly yields the original problem. Let us rewrite the inner inf as a function of γ. Canceling the gradients in f and g, we get

f = 1 2λ 2 ( i=1 γ i φ X (x i ) -ŵµ ) g = 1 2λ 2 ( i=1 γ i φ Y (ỹ i ) -ŵν ), (96)
i.e. the primal-dual relations Eq. ( 40). Let us now derive the optimality condition on B: as in the proof of of Lemma 11, we have

inf B∈S + (R ) i=1 γ i Φ T i BΦ i + λ 1 2 B 2 F = - 1 2λ 1 (- i=1 γ i Φ i Φ T i ) + 2 F . (97)
Plugging Eq. ( 40) and Eq. ( 97) into Eq. ( 95), we obtain Eq. ( 39). Let us now derive the smoothness and strong convexity constants of

H(γ) def = 1 4λ 2 γ T Qγ - 1 2λ 2 j=1 γ j z j + 1 2λ 1 (- i=1 γ i Φ i Φ T i ) + 2 F + 2δ γ 2 + q 2 4λ 2 . ( 98
)
Strong convexity. Let us recall that H is α-strongly convex i.f.f. ∀γ, γ , α γ -γ 2 ≤ ∇ H(γ) -∇H(γ ), γ -γ . We have

H(γ) -∇H(γ ), γ -γ = l δ γ -γ 2 + 1 2λ 2 (γ -γ ) T Q(γ -γ ) + 1 λ 2 diag(Φ T ([ i=1 γ i Φ i Φ T i ] --[ i=1 γ i Φ i Φ T i ] -)Φ) T (γ -γ ).

  Assumption 1 (m-times differentiable densities). Let m, d ∈ N. Let µ, ν ∈ P(R d ). a) µ, ν have densities. Their supports, resp. X, Y ⊂ R d are convex, bounded and open with a Lipschitz boundary; b) the densities are finite, bounded away from zero, with Lipschitz derivatives up to order m.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: 2D Gaussian data. (Top left) OT error | OT -OT |. (Top right) Best hyperparameters (selected via gridsearch), log scale. (Bottom left) Transportation map mean square error, log scale. (Bottom right) MMD between transported samples, log scale. Shaded areas correspond to ± std. Algorithm: accelerated gradient descent (δ = 10 3 ), Sobolev kernel (s = 20, bandwidth = 1), with Nystrom approximation (rank = 100). Filling pairs (x i , ỹi ) are drawn from µ ⊗ ν. The number of filling pairs is equal to the number of µ and ν samples, reported on the x-axis.

+

  Tr( Ân ) bounded ; the rate in λ n immediately follows the boundedness. Denoting b n = Rn + R and c n = λ n (Tr(A * ) + R 2 ) + κn,δ , we have the following system of inequalities split the analysis into two parts: the indexes for which a n ) ≤ c n (Case 2) .

  the first case, we sub-split the analysis in two parts: the indexes for which On these indexes, we can re-write Equation (75) as   a 2 n ≤ 4KC b n log( 2 δ ) √ n a n λ n ( R2 n + Tr( Ân )) ≤ 4(2C + 1)

Compared to Vacher et al. (2021), the definition of z above is slightly different. This is due to the fact that we consider the "Brenier" formulation of dual OT (an infimum with an inner product cost), instead of the Kantorovich formulation (a supremum with a quadratic cost).

Acknowledgements

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). We also acknowledge support from the European Research Council (grants SEQUOIA 724063 and REAL 947908), and Région Ile-de-France.

It follows that the semi-dual functional J(f ) = f, µ + f * , ν is convex.

Denoting T 0 the optimal transport map from µ to ν, the semi-dual functional can be rewritten as J(f ) = f, µ + f * • T 0 , µ . Now, the Fenchel-Young inequality on f gives for every couple (x, y) ∈ X × Y y x ≤ f (x) + f * (y) , Equality holds for y = ∇f (x), and in this case, one also have x ∈ ∂f * (∇f (x)). To simplify notations, let us denote T (x) = ∇f (x). We get f (x) + f * (T (x)) = T (x) x .

(50)

Denoting f 0 an optimal potential, the optimality condition applied to f 0 gives f 0 (x) + f * 0 (T 0 (x))-T 0 (x) x = 0 , ∀x ∈ Supp(µ) and by integration J(f 0 ) = T 0 (x) x dµ . Therefore, we have

Subtracting Eq. ( 50) applied to f pointwise from the integrand in Eq. ( 51), we get

Applying the case of equality in Fenchel-Young (see above), we have x ∈ ∂f * (T (x)) where ∂f * is the sub-gradient of f * . Hence, denoting z 1 = T 0 (x) and z 2 = T (x), the integrand has the form

Since f is a C 1 convex function with a M -Lipschitz gradient, then f * is 1 M -strongly convex and ( 53) is lower bounded as

and in particular, we recover

A.2 Proof of Proposition 4

Proof For E ⊂ R d with lipschitz boundary and points (e i ) 1≤i≤n , define the fill distance h(E, e) as h(E, e) = sup

If f has smoothness m, the sampling inequalities [START_REF] Narcowich | Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting[END_REF] state that if ∀i ∈ [n], f (e i ) = 0 then there exists a positive constant B depending on the smoothness m, the dimension d and the geometry of

Now, recall that at the optimum, the function

Using Lemma 12 from [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF], if n ≥ n 0 (X, Y ) there exists a universal constant C 0 such that we have with probability at least 1

There remains to upper bound αn H m (X×Y )

where c(x, y)

where

where C 2 is a constant depending on X, Y, m, d. Hence we obtain

where

In particular, since ∀(x, y) ∈ X × Y, φ(x, y), Ân φ(x, y) H m (X×Y ) ≥ 0, we obtain that for n ≥ n 0 (X, Y, m, d), we have with probability at least 1 -δ

where

A.3 Proof of Lemma 5

Proof In order to apply the Poincaré inequality, we need to re-normalize the function u -v such that it integrates to 0. Denoting the residual r = X u(x) -v(x)dx, we have that u -v, µ -μ = u -v -r, µ -μ since µ and μ have the same total mass. For a probability measure α and a RKHS H with kernel k, we define kernel mean embedding w α = E α (k(X, .)) ∈ H. With this notation, using the reproducing property, we can re write u -v -r, µ -μ as u -v -r, w µ -w μ H d/2+ε (X) . Using Cauchy-Schwarz, we have the upper-bound

Combining these two equations yields

n , we recover

which implies that Rn and Tr( Ân ) are bounded independently on δ and that

where b * = sup n b n .

Case 1b. On these indexes, the inequalities we obtain are

The first equation implies that

Hence the upper bound we obtain on λ n ( R2 n + Tr( Ân )) is

= 4(2C + 1)(4KC )

gives that Rn and Tr( Ân ) are bounded (independently on δ) and yields the rate

Case 2. On these indexes, the second equation of (75) becomes

Now recall that c n is given by

where C 1 , G are constants that do not depend of n, δ. Hence, if we choose λ n such that

, the quantities Rn and Tr( Ân ) are bounded independently on δ and we recover a n of the form

where K is a constant independent on n and δ.

, we have that in any case, the energies are bounded and we get

where C is a constant independent on n and δ.

B Proofs from Section 4

B.1 Proof of Lemma 11

Proof Not including the PSD constraints on B, problem (34) admits the following Lagrangian:

Canceling the gradients in f and g, we get

Let us now derive the optimality condition on B: completing the square, we have

Given that the term with negative parts is non-negative but vanishes when the corresponding matrices are positive, we get the lower bound α ≥ δ + 1 2λ 2 λ min (Q). Smoothness. Finally, we have

Bounded those three terms independently yields