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

Alessandro Rudi
email: alessandro.rudi@inria.fr

Francis Bach
email: francis.bach@inria.fr

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

Mikhail Belkin

Samory Kpotufe

B Muzellec
email: boris.muzellec@inria.fr

F.-X Vialard

Muzellec Vacher

Bach Rudi

Vialard

A Dimension-free Computational Upper-bound for Smooth Optimal Transport Estimation

It is well-known that plug-in statistical estimation of optimal transport suffers from the curse of dimensionality. Despite recent efforts to improve the rate of estimation with the smoothness of the problem, the computational complexity of these recently proposed methods still degrade exponentially with the dimension. In this paper, thanks to an infinitedimensional sum-of-squares representation, we derive a statistical estimator of smooth optimal transport which achieves a precision ε from Õ(ε -2 ) independent and identically distributed samples from the distributions, for a computational cost of Õ(ε -4 ) when the smoothness increases, hence yielding dimension-free statistical and computational rates, with potentially exponentially dimension-dependent constants.

Introduction

The comparison between probability distributions is a fundamental task and has been extensively used in machine learning. For this purpose, optimal transport (OT) has recently gained traction in different subfields of machine learning (ML), such as natural language processing (NLP) [START_REF] Xu | Distilled Wasserstein learning for word embedding and topic modeling[END_REF][START_REF] Chen | Adversarial text generation via feature-mover's distance[END_REF], generative modeling (Arjovsky et al., 2017;[START_REF] Tolstikhin | Wasserstein auto-encoders[END_REF][START_REF] Salimans | Improving GANs using optimal transport[END_REF], multi-label classification [START_REF] Frogner | Learning with a Wasserstein loss[END_REF], domain adaptation [START_REF] Redko | Optimal transport for multisource domain adaptation under target shift[END_REF], clustering [START_REF] Ho | Multilevel clustering via Wasserstein means[END_REF], and has had an impact in other areas such as imaging sciences [START_REF] Feydy | Optimal Transport for Diffeomorphic Registration[END_REF][START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF]. Indeed, OT is a tool to compare data distributions which has arguably many more geometric properties than other available divergences [START_REF] Peyré | Computational optimal transport[END_REF].

State of the art. An approach that has first been advocated in the machine learning community as a way to efficiently approximate empirical OT and to make it differentiable consists in adding an entropic regularization term to the OT problem [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. Rates on the sample complexity of entropic optimal transport have then been studied by [START_REF] Genevay | Sample complexity of Sinkhorn divergences[END_REF] and [START_REF] Mena | Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem[END_REF], and were proven to be of the order O(

1 ε d/2 √ n ),
for small values of ε. Although the dependency in the number of samples is in 1/ √ n, the constant degrades exponentially with respect to the dimension. Entropic OT and Sinkhorn divergences [START_REF] Genevay | Learning generative models with Sinkhorn divergences[END_REF] were then leveraged as a tool to study the sample complexity of the (unregularized) OT problem itself: the most advanced results in this direction were derived by [START_REF] Chizat | Faster Wasserstein distance estimation with the Sinkhorn divergence[END_REF], who show that with few assumptions on the regularity of the Kantorovich potentials, the squared Wasserstein distance can be estimated using O(ε -d/2+2 ) samples and Õ(ε -(d +5.5) ) operations (with d = 2 d 2 ) with high probability. Our work can be related to a current research direction which consists in developing estimators of the Wasserstein distance for classes of smooth distributions, with smoothness parameter m, that have better performances than in the general case. Related to this trend, [START_REF] Liang | On how well generative adversarial networks learn densities: Nonparametric and parametric results[END_REF][START_REF] Liang | Estimating certain integral probability metric (IPM) is as hard as estimating under the IPM[END_REF] showed minimax rates for a class of integral probability metrics (IPM) that includes the Wasserstein-1 distance, as a function of the smoothness of the distributions. However, (i) for p > 1, the Wasserstein-p distance W p is not an IPM and (ii) no estimators with matching rates are proposed in those two works. So far, two main contributions leveraging smoothness that are applicable to the W 2 distance can be found. [START_REF] Hütter | Minimax rates of estimation for smooth optimal transport maps[END_REF] derive minimax rates for the estimation of the OT maps and propose an estimator which necessitates, for an L 2 error on the maps of order ε, O(ε -2m-1+d/2 2m ) samples. While statistically almost optimal, this estimator is not computationally feasible as it requires to project the potentials on a space of smooth, strongly-convex functions. Instead, [START_REF] Weed | Estimation of smooth densities in Wasserstein distance[END_REF] derive estimators for the densities requiring O(ε -d+2m 1+m ) samples and, under the assumption that an efficient resampler is available, derive an estimator of the OT distance that can be calculated in Õ(ε -(2d+d/2) ) time.

While the contributions above do succeed in taking advantage of the smoothness from a statistical point of view, they do not manage to take advantage of the smoothness from a computational point of view. Actually, statistical-computational gaps are known to exist for some instances of high-dimensional OT, such as the spiked transport model of [START_REF] Niles | Estimation of Wasserstein distances in the spiked transport model[END_REF].

Contributions. In this paper, we bridge the statistical-computational gap of smooth OT estimation and we provide a positive answer to the question whether smoothness of the optimal potentials can be computationally beneficial to an efficient statistical estimator. More precisely, we propose an algorithm which, for a given accuracy ε, needs O(ε -2 ) samples and has a computational complexity of Õ(ε -max(4, 7d m-d ) ). Note that the computational complexity improves with the regularity of the distributions and, when m ≥ 3d, it is Õ(ε -4 ), i.e., independent of the dimension d in the exponent (but not in the constants). We thus show that smoothness can be leveraged in the computational estimation of optimal transport. Moreover, we consider different scenarios beyond i.i.d. sampling, such as the case where we are able to compute exact integrals or where we can evaluate the densities in given points, by representing the problem in terms of kernel mean embeddings [START_REF] Muandet | Kernel Mean Embedding of Distributions: A Review and Beyond[END_REF]. This allows to make a unified analysis for all the cases. The total error is then the sum of the error induced by approximating via the kernel mean embedding plus the error induced by subsampling the constraints. Interestingly, in the other scenarios the computational cost to achieve an error ε can be smaller than ε -4 , as reported below. This is particularly interesting in the case we can evaluate the densities in given points and avoids using expensive Monte-Carlo sampling techniques to obtain i.i.d. samples (see Section 6 for more details). Our results are summarized below.

Theorem 1 Let ε > 0. Let µ, ν satisfy Assumption 1 for some m > d. Let OT be the proposed estimator defined in Eq. ( 6) and computed as in Appendix F with the same parameters as in Corollary 3. The cost to achieve | OT -OT(µ, ν)| ≤ ε for the three scenarios is:

1. (Exact integral) Time: Õ(ε -7d m-d ). Space: Õ(ε -4d m-d ).
2. (Evaluation) Time:

Õ(ε -7d m-d ). Space: Õ(ε -4d m-d ). #evaluations of µ, ν: Õ(ε -d m+1 ). 3. (Sampling) Time: Õ(ε -max(4, 7d m-d ) ). Space: Õ(ε -4d m-d ). #samples of µ, ν: Õ(ε -2 ).
The second key contribution of this paper is to provide a new representation theorem for solutions of smooth optimal transport. The inequality constraint in the dual OT problem can be replaced with an equality constraint involving a finite sum-of-squares in a Sobolev space.

In comparison with [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], it is a non-trivial extension of their representation result to the case of a continuous set of global minimizers instead of a finite set.

Sketch of the result and derivation of the algorithm

In this paper, we consider the optimal transport problem for the quadratic cost on bounded subsets X, Y of the Euclidean space R d . The set of probability measures on X is denoted by P(X). The optimal transport problem with quadratic cost c(x, y) = 1 2 x -y 2 can be stated in its dual formulation as

OT(µ, ν) = sup u,v∈C(R d ) u(x)dµ(x) + v(y)dν(y) subject to c(x, y) ≥ u(x) + v(y), ∀(x, y) ∈ X × Y, (1)
As a standard result in optimal transport theory, the supremum is attained and the functions u , v are referred to as the Kantorovich potentials (see [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling[END_REF].

The proposed approach to approximate OT(µ, ν) is the result of two main ingredients: (1) a suitable way to represent smooth functions and to approximate their integral in µ, ν, (2) a way to enforce efficiently the dense set of constraints on u, v.

Preliminary step: Representing smooth functions and integrals. We represent smooth functions via a reproducing Kernel Hilbert space (RKHS) [START_REF] Aronszajn | Theory of reproducing kernels[END_REF][START_REF] Steinwart | Support Vector Machines[END_REF], for which functions can be represented as linear forms. In Section 4 we show that under smoothness assumptions on µ and ν (Assumption 1) we have u ∈ H X and v ∈ H Y where H X and H Y are two suitable RKHSs on X and Y , associated with two bounded continuous feature maps φ X : X → H X and φ Y : Y → H Y . Note that RKHSs offer several advantages. First, leveraging the reproducing property, we can represent the integrals in the functional of Eq. (1) as inner products in terms of the kernel mean embeddings w µ ∈ H X and w ν ∈ H Y where w µ = X φ X (x)dµ(x) and w ν = Y φ Y (y)dν(y). Indeed, by the reproducing property, for all u ∈ H X , we have:

X u(x)dµ(x) = X u, φ X (x) H X dµ(x) = u, X φ X (x)dµ(x) H X = u, w µ H X ,
and the same reasoning holds for the integral on ν, i.e., Y v(y)dν(y) = v, w ν H Y , for all v ∈ H Y . This construction is known as kernel mean embedding [START_REF] Muandet | Kernel Mean Embedding of Distributions: A Review and Beyond[END_REF]. Moreover, RKHSs allow the so-called kernel trick [START_REF] Steinwart | Support Vector Machines[END_REF], i.e., to express the resulting algorithm in terms of kernel functions that in our case correspond to The main step: Dealing with a dense set of inequalities. Even assuming that we are able to compute integrals in closed form and restricting to m-times differentiable u, v, the main challenge is to deal with the dense set of inequalities c(x, y) ≥ u(x) + v(y) that u, v must satisfy, for any (x, y) ∈ X × Y . Indeed, an intuitive approach would be to subsample the set, i.e., to take points (x 1 , ỹ1 ), . . . , (x , ỹ ) in X × Y and consider only the constraints c(x j , ỹj ) ≥ u(x j ) + v(ỹ j ) for j = 1, . . . , . This approach, however, is only able to leverage the Lipschitzianity of u, v (Aubin-Frankowski and Szabo, 2020) and leads to an error in the order of -1/d that does not allow to break the dependence in d in the exponent, and yields rates that are equivalent to the plugin estimator.

k X (x, x ) = φ X (x), φ X (x ) H X and k Y (y, y ) = φ Y (y), φ Y (y ) H Y ,
In this paper, we leverage a more refined technique to approximate the dense set of inequalities, introduced by [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] for the problem of non-convex optimization, and that allows to break the curse of dimensionality for smooth problems. The idea behind this technique is the consideration that, while a dense set of inequalities is poorly approximated by subsampling, the situation is different in the case of a dense set of equality constraints, for which an optimal rate of O( -m/d ) is achievable for m-times differentiable constraints [START_REF] Wendland | Approximate interpolation with applications to selecting smoothing parameters[END_REF]. The construction works in two steps: first, substitute the inequality constraints with equality constraints that are equivalent, and then subsample. In the next two paragraphs we explain how to apply this approach to the problem of OT.

Removing the inequalities: positive definite operator characterization. To apply the construction recalled above to our scenario, we first consider the following problem. Let H XY be a Hilbert space on X × Y and φ : X × Y → H XY . Denote by k XY the kernel k XY ((x, y), (x , y )) = φ(x, y), φ(x , y ) H XY for any (x, y), (x , y ) ∈ X × Y and by S + (H XY ) the space of positive operators on H XY . We define

max u∈H X ,v∈H Y , A∈S + (H XY ) u, w µ H X + v, w ν H Y subject to ∀(x, y) ∈ X × Y, c(x, y) -u(x) -v(y) = φ(x, y), Aφ(x, y) H XY , (2)
where the inequality in (1) is substituted with an equality w.r.t. a positive definite operator A on H XY . Note that Problem (1) is a relaxation of Problem (2): indeed, if for a given pair u ∈ H X , v ∈ H Y there exists a positive definite A satisfying the equality above, then c(x, y) -u(x) -v(y) = φ(x, y), Aφ(x, y) H XY ≥ 0, ∀(x, y) ∈ X × Y, so the couple (u, v) is admissible for (1). However, even for an admissible couple in (1) satisfying u ∈ H X , v ∈ H Y , a positive operator A may not exist. Indeed, note that the technique of representing a positivity constraint in terms of a positive matrix has a long history in the community of polynomial optimization [START_REF] Bernard | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF][START_REF] Bernard | Moments, Positive Polynomials and Their Applications[END_REF], which shows that in general the resulting problem is not equivalent to the original one, for any chosen degree of polynomial approximation. This fact leads to the so-called sum of squares hierarchies, also used for optimal transport [START_REF] Henrion | Graph recovery from incomplete moment information[END_REF]. Instead, using kernels, [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] showed that there exists a positive operator with finite rank that matches the constraints and makes the two problems equivalent, when the constraint is attained on a finite set of points. However, such existence results cannot be used for the problem in (2), since in the case of optimal transport the set of zeros corresponds to the graph of the optimal transport map and is a smooth manifold, when µ, ν are smooth [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF].

A crucial point of our contribution is then to prove that, with a quadratic cost c(x, y) = maximizing (1), making the two problems equivalent. The result is reported in Theorem 5. The proof is derived using the Fenchel dual characterization of u , v and gives a sharp control of the rank of A.

Subsampling the constraints and approximating the integrals. We restrict the constraint of (2) to (x 1 , ỹ1 ), . . . , (x , ỹ ) ⊂ X × Y for ∈ N. However, we need to add a penalization for u, v and A to avoid overfitting, since the error induced by subsampling the constraints is proportional to the trace of A [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] and, in our case, also to the norms of u, v, as derived in Theorem 9 in Section 5. Finally, in two of the three scenarios of interest for the paper, i.e., (i) when we can only evaluate µ, ν pointwise, or (ii) when we have only i.i.d. samples from µ, ν, we do not have access to the kernel mean embeddings w µ ∈ H X , w ν ∈ H Y . Therefore, we need to use some estimators ŵµ ∈ H X , ŵν ∈ H Y that are derived in Section 6. The resulting problem is the following, for some regularization parameters λ 1 , λ 2 > 0:

max u∈H X ,v∈H Y , A∈S + (H XY ) u, ŵµ H X + v, ŵν H Y -λ 1 Tr(A) -λ 2 ( u 2 H X + v 2 H Y ) subject to ∀j ∈ [ ], c(x j , ỹj ) -u(x j ) -v(ỹ j ) = φ(x j , ỹj ), Aφ(x j , ỹj ) H XY .
(3)

Let û, v be the maximizers of the problem above (unique since the problem is strongly concave in u, v). The estimator for OT we consider corresponds to

OT = û, ŵµ H X + v, ŵν H Y . (4)
Finite-dimensional characterization. In Appendix F, following Marteau-Ferey et al.

(2020), we derive the dual problem of Eq. (3

). Define Q ∈ R × as Q i,j = k X (x i , xj )+k Y (ỹ i , ỹj )
and z j = ŵµ (x j ) + ŵν (ỹ j ) -2λ 2 c(x j , ỹj ) for i, j ∈ [ ] and q 2 = ŵµ 2

H X + ŵν 2
H Y , and let I ∈ R × be the identity matrix. Let K ∈ R × be defined as K i,j = k XY ((x i , ỹi ), (x j , ỹj )) and define Φ i ∈ R as the i-th column of R, the upper triangular matrix corresponding to the Cholesky decomposition of K (i.e., R satisfies K = R R). The dual problem writes:

min γ∈R 1 4λ 2 γ Qγ -1 2λ 2 j=1 γ j z j + q 2 4λ 2 such that j=1 γ j Φ j Φ j + λ 1 I 0. (5)
In the same section in Appendix F, we derive an explicit characterization of û, ŵ, Â in terms of γ, the solution of the problem above and we characterize OT as follows:

OT = q 2 2λ 2 -1 2λ 2 j=1 γj ( ŵµ (x j ) + ŵν (ỹ j )). ( 6
)
As it is possible to observe from the problem above and the characterization of OT, the only quantities necessary to compute γ and OT are the kernels k X , k Y , k XY and the evaluation of the functions ŵµ ∈ H X , ŵν ∈ H Y at the points xj and ỹj respectively for j ∈ [ ]. In Appendix F, we consider a Newton method with self-concordant barrier to solve the problem above [START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF]. To illustrate that this algorithm can indeed be implemented in practice, we run simulations on toy data in Appendix G. The total cost of the procedure to achieve error ε for the computation of OT is the following (see Theorem 17, Page 27 in the appendix):

O C + E + 3.5 log ε time, O( 2 ) memory, ( 7
)
where C is the cost for computing q 2 and E is the cost to compute one z j . Depending on the operations that we are able to perform on µ, ν and k X , k Y , we have three scenarios.

In Section 6 we specify how to compute the vectors ŵµ , ŵν , in Corollary 3 we report only the conditions of applicability and the resulting cost. In the next section and then in Appendix F we quantify instead how to choose , λ 1 , λ 2 to achieve | OT -OT(µ, ν)| ≤ ε with high probability and we provide a complete computational complexity in ε.

Theoretical Guarantees

Here we quantify the convergence rate of OT to OT. To simplify the exposition, in this section we will make a classical assumption on the smoothness of the densities [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF]. Note however that the results of the paper hold under a more general assumption on the smoothness of the potentials (see Theorem 5). Assumption 1 is particularly adapted to our context since it guarantees that the Kantorovich potentials have a similar order of differentiability [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF]. The main result, Theorem 9, is expressed for a general set of couples

Assumption 1 (m-times differentiable densities) Let m, d ∈ N. Let µ, ν ∈ P(R d ).
(x j , ỹj ), j ∈ [ ].
Here, we specify it for the case where the couples are sampled independently and uniformly at random.

λ 1 ≥ C 1 -m/2d+1/2 log δ , λ 2 ≥ w µ -ŵµ H X + w ν -ŵν H Y + λ 1 ,
then, with probability 1 -δ, we have

| OT -OT(µ, ν)| ≤ C 2 λ 2 .
Note that while the rate does not depend exponentially in d as we will see in the rest of the section, the constants 0 , C 1 , C 2 depend exponentially in d in the worst case, as [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] for the case of global optimization. From the theorem above it is clear that the approximation error of OT is the sum of the error induced by the kernel mean embeddings plus the error induced by the subsampling of the inequality. Note here that the result of the theorem above holds also if the couples are i.i.d. from ρ = µ ⊗ ν, as discussed in Remark 10. This can be beneficial if we do not know X, Y or we do not know how to sample from them.

In the next corollary we will specialize the result depending on the considered scenarios.

Corollary 3 Under the same assumptions as Theorem 2, let k

X = k m+1 , k Y = k m+1 , k XY = k m
where k s for s > 0 is defined in Eq. ( 8) and λ 1 ≥ C 1 -(m-d)/2d log δ . Compute OT with ŵµ , ŵν chosen according to one of the three scenarios below, as in Section 6. There exist C, C , C 2 , C 2 s.t. with probability at least 1 -δ, 1. (Exact integral) When we are able to compute exactly k X (x, x )dµ(x )dµ(x) and also X k X (x, x )dµ(x) for any x ∈ X (and analogously for ν). Choose λ 2 = λ 1 . Then,

| OT -OT(µ, ν)| ≤ C 2 -(m-d)/2d log δ .
2. (Evaluation) When we are only able to evaluate µ, ν on given points and to compute

X k X (x, z)k X (x , z)dz, X X k X (x, z)k X (z, z )k X (x , z )dzdz . Evaluate µ in n µ points sampled uniformly from X (and n ν for ν). Let λ 2 = λ 1 + C(n µ + n ν ) -(m+1)/d log nµ+nν δ , | OT -OT(µ, ν)| ≤ C 2 (n -(m+1)/d µ log nµ δ + n -(m+1)/d ν log nν δ + -(m-d)/2d log δ ).
3. (Sampling) When we are only able to sample from µ, ν, by using n µ i.i.d. samples from µ and n ν from ν. Choose

λ 2 = λ 1 + C (n µ + n ν ) -1/2 log nµ+nν δ
. Then,

| OT -OT(µ, ν)| ≤ C 2 (n -1/2 µ log nµ δ + n -1/2 ν log nν δ + -(m-d)/2d log δ ).

Notations and background

Let n ∈ N, we denote by [n] the set {1, . . . , n}. For a set Z, and a positive definite kernel k : Z × Z → R (i.e., so that all matrices of pairwise kernel evaluations are positive semidefinite), we can define a reproducing kernel Hilbert spaces [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]) H of real functions on Z, endowed with an inner product [START_REF] Aronszajn | Theory of reproducing kernels[END_REF].

: Z → H corresponding to z → k(z, •), so that k(z, z ) = φ(z), φ(z ) H
In this paper we use Sobolev spaces, defined on Z ⊆ R d , with d ∈ N, an open set. For s ∈ N, denote by H s (Z) the Sobolev space of functions whose weak derivatives up to order s are square-integrable, i.e., H s (Z) : [START_REF] Adams | Sobolev Spaces[END_REF]. A remarkable property of H s (Z) that we are going to use in the proofs is that

= {f ∈ L 2 (Z) | f H s (Z) < ∞} and f H s (Z) := |α|≤s D α f L 2 (Z)
H s (Z) ⊂ C k (Z) for any s > d/2 + k and k ∈ N. Moreover H m+1 (Z) ⊂ H m (Z), ∀m ∈ N.
Proposition 4 (Sobolev kernel, [START_REF] Wendland | Scattered Data Approximation[END_REF]

) Let Z ⊂ R d , d ∈ N, be an open bounded set. Let s > d/2. Define k s (z, z ) = c s z -z s-d/2 K s-d/2 ( z -z ), ∀z, z ∈ Z, (8)
where K : R → R is the Bessel function of the second kind (see, e.g., Eq. 5.10 of the same book) and c s = 2 1+d/2-s Γ(s-d/2) . Then the function k s is a kernel. Denoting by H Z the associated RKHS, when Z has Lipschitz boundary, then H Z = H s (Z) and the norms are equivalent.

In the particular case of s = d/2 + 1/2, we have k s (z, z ) = exp(-z -z ). Note that the constant c s is chosen such that k s (z, z) = 1 for any z ∈ Z.

Positive operator representation for Kantorovich potentials

We start with the following representation result on the structure of the optimal potentials, which is one of our main contributions in this paper.

Theorem 5 Let µ ∈ P(X), ν ∈ P(Y ) satisfying Assumption 1a and let (u , v ) be Kantorovich potentials such that u ∈ H s+2 (X) and v ∈ H s+2 (Y ) for s > d + 1. There exist functions w 1 , ...,

w d ∈ H s (X × Y ) such that 1 2 x -y 2 -u (x) -v (y) = d i=1 w i (x, y) 2 , ∀(x, y) ∈ X × Y. Proof Denote by h the function h(x, y) = c(x, y) -u (x) -v (y) for all (x, y) ∈ X × Y . Let f (x) = 1 2 x 2 -u (x),
x ∈ X. By Brenier's theoremfor quadratic optimal transport [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF]Santambrogio, 2015, Theorem 1.22), (i) T = ∇f , where T is the optimal transport map from µ to ν,

(ii) f is convex on X, (iii) h is characterized by h(x, y) = f (x) + f (y) -x y, where f : y ∈ Y → sup x∈X x y - f (x) is the Fenchel-Legendre conjugate of f . Moreover, f (y) = 1 2 y 2 -v (y).
Further, from the properties of Fenchel-Legendre conjugacy [START_REF] Rockafellar | Convex Analysis[END_REF], Section 26), we have

T -1 = ∇f . Hence, since u ∈ H s+2 (X) and v ∈ H s+2 (Y ), we have (iv) T = ∇f (resp. T -1 = ∇f ) is a H s+1 -diffeomorphism from X to Y (resp. Y to X).
Since f ∈ H s+2 (X) and s > d/2 and X is a bounded open set with locally Lipschitz boundary (see Lemma 11), we have H s+2 (X) ⊂ C 2 (X) [START_REF] Adams | Sobolev Spaces[END_REF] and the Hessian H f (x) is well defined for any x ∈ X. Since, by (iv) , T = ∇f is a diffeomorphism, then, by Fenchel-Legendre conjugacy, f is strictly convex [START_REF] Rockafellar | Convex Analysis[END_REF]. Hence by compactness of X, f has a Hessian H f (x) which is bounded away from 0. This implies:

(v) There exists ρ > 0 such that H f (x) ρ Id for all x ∈ X.

We will now use the decomposition (iii) along with a reparameterization of h to obtain a decomposition as a sum of squares. Let

h(x, z) def = h(x, T (z)), ∀(x, z) ∈ X × X.
The effect of this change of coordinates is illustrated in Figure 1. Since f is differentiable, by the properties of the Fenchel-Legendre conjugate it holds that f (∇f (z)) = ∇f (z) z -f (z) for any z ∈ X [START_REF] Rockafellar | Convex Analysis[END_REF]. Therefore, from (i) we have that

h(x, z) = f (x) -f (z) -∇f (z) (x -z).
Now, since X is convex, we can characterize f in terms of its Taylor expansion: where R is the integral reminder R(x, z)

f (x) = f (z) + ∇f (z) (x -z) + (x -z) R(x, z)(x -z), ∀x, z ∈ X,
def = 1 0 (1 -t)H f ((1 -t)x + tz)dt. This implies h(x, z) = (x -z) R(x, z)(x -z), ∀x, z ∈ X.
From (v), we get ∀x, z, R(x, z) ρ 2 Id. In particular, for all x, z ∈ X, the matrix R(x, z) admits a positive square root R(x, z). Further, since

√ • is C ∞ on the closed set {A ∈ S + (R d ) : A ρ 2 Id} and ∂ 2 ∂x i ∂x j f ∈ H s (X) for all i, j ∈ [d]
, the functions r i,j : (x, z) → e i R(x, z)e j are in H s (X × X) for all i, j ∈ [d] (see Proposition 1 and Assumption 2b of [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], where (e 1 , ..., e d ) is the canonical ONB of R d . Define now the functions wi (x, z)

def = d j=1 r i,j (x, z)(e j (x -z)), ∀x, z ∈ X, i ∈ [d].
From the above arguments, it holds that wi ∈ H s (X×X), i ∈ [d], and h(x, z)

= d i=1 w2 i (x, z). Now, since T is a H s+1 -diffeomorphism from X to Y , changing parameterization again we have h(x, y) = d i=1 w 2 i (x, y), ∀(x, y) ∈ X × Y , with w i (x, y) = wi (x, T -1 (y)), ∀(x, y) ∈ X × Y . We conclude by proving that w i ∈ H s (X × Y ) for all i ∈ [d]. Indeed, from (iv) T -1 is a H s+1 -diffeomorphism from Ȳ to X
and it is bi-Lipschitz (since T and T -1 are Lipschitz due to the continuity of their Hessian and the boundedness of X, Y ). Define the map Q as (x, y) → (x, T -1 (y)) and note that Q ∈ H s+1 (X × Y, R 2d ), by construction. Note that wi ∈ H s (X × X) and has bounded weak derivatives of order 1, since s > d + 1 and H s (X × X) is bounded [START_REF] Adams | Sobolev Spaces[END_REF]. The conditions above on wi , Q guarantee that w i = wi • Q belongs to H s (X × Y ) (see, e.g., Theorem 1.2 of [START_REF] Campbell | The weak inverse mapping theorem[END_REF].

Theorem 5 implies the existence of A ∈ S + (H XY ) by effect of the reproducing property, when we consider a RKHS H XY containing H s (X × Y ). The proof is in Appendix A, Page 17.

Corollary 6 Let H XY be a RKHS such that H s (X × Y ) ⊆ H XY . Under the hypothesis of Theorem 5, there exists a positive operator A ∈ S + (H XY ) with rank at most d, such that c(x, y) -u (x) -v (y) = φ(x, y), A * φ(x, y) H XY .

(9)

Finally, the following corollary shows the effect of Assumption 1 on the existence of A .

Note indeed that such an assumption implies smoothness of the Kantorovich potentials [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF]. If m > d they are smooth enough to apply Theorem 5 and the corollary above. The proof of the following corollary is in Appendix A, Page 17.

Corollary 7 (Effects of Asm. 1) Let µ, ν satisfy Assumption 1 for m > d. Let (u , v ) be Kantorovich potentials for µ, ν. Let

H X = H m+3 (X), H Y = H m+3 (Y ), H XY = H m (X × Y ).
Then u ∈ H X , v ∈ H Y and there exists A ∈ S + (H XY ) satifying Eq. ( 9) and rank A ≤ d.

Subsampling the constraints

Given ∈ N and a set of points Z = {(x 1 , ỹ1 ), . . . , (x , ỹ )}, define the fill distance [START_REF] Wendland | Scattered Data Approximation[END_REF],

h = sup x∈X,y∈Y min j∈[ ] (x, y) -(x j , ỹj ) . ( 10
)
The following theorem, that is an adaptation of Theorem 4 from [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], quantifies the error of subsampling the constraints in terms of the fill distance.

Theorem 8 Let X, Y satisfy Assumption 1a. Let s ≥ 3 and s

> d. Let H X ⊆ H s (X), H Y ⊆ H s (Y ), H XY = H s (X × Y ). Let Z ⊂ X × Y be a set of points of cardinality and fill distance h and let u ∈ H X , v ∈ H Y , A ∈ S + (H XY ) satisfy c(x j , ỹj ) -u(x j ) -v(ỹ j ) = φ(x j , ỹj ), Aφ(x j , ỹj ) H XY , ∀j ∈ [ ]. (11)
There exist h 0 , C 0 depending only on s, d, X, Y such that, when h ≤ h 0 , then

c(x, y) ≥ u(x) + v(y) -ε, ∀x, y ∈ X × Y, where ε = Qh s-d , where Q = C 0 ( u H X + v H Y + Tr(A)). Note that h 0 , C 0 depend only on d, m, X, Y .
The proof of the theorem above is reported in Appendix B, Page 17. Using the theorem above we are able to show that, given a maximizer (û, v, Â) of Problem 3 the couple (û-ε/2, v -ε/2) is admissible for Problem 1. This is a crucial step to bound | OT -OT| in terms of the fill distance h and it is stated next.

Theorem 9 Let µ, ν and X, Y ⊂ R d satisfy Assumption 1a. Let s > d and let H X , H Y , H XY be RKHS on X, Y, X ×Y such that H X ⊆ H s (X), H Y ⊆ H s (Y ), H XY = H s (X ×Y ).
Assume that there exist two Kantorovich potentials u , v such that u ∈ H X , v ∈ H Y and there exists A ∈ S + (H XY ) that satisfies Eq. (9). Let OT be computed according to Eq. (6) using a set of ∈ N points Z ⊂ X × Y with fill distance h . Let h 0 , C 0 as in Theorem 8. Let

η = C 0 h s-d , γ = w µ -ŵµ H X + w ν -ŵν H Y . When h ≤ h 0 , λ 2 > 0 and λ 1 is chosen such that λ 1 ≥ 2η, then | OT -OT| ≤ 6λ 1 Tr(A ) + 6 (γ+η) 2 λ 2 + 6 λ 2 u 2 H X + v 2 H Y ).
The proof of the theorem above is in Appendix B, Page 17. Theorem 9 together with Theorem 5 (in particular Corollary 7) allow to prove Theorem 2. To bound h in terms of , we used a result recalled in Lemma 12.

Proof of Theorem 2. First note that k

X = k m+1 , k Y = k m+1 , k XY = k m imply that H X = H m+1 (X), H Y = H m+1 (X), H XY = H m (X × Y ).
Then, by Corollary 7 we have that under Assumption 1 for any Kantorovich potentials u , v there exists a finite rank A ∈ S + (H XY ) such that c(x, y) -u (x) -v (y) = φ(x, y), A * φ(x, y) H XY , and that u ∈ H m+3 (X) ⊆ H m+1 (X) = H X and analogously v ∈ H Y . Among them, we select the triplet minimizing u 2 H X + v 2 H Y + Tr(A ) and we denote by Q µ,ν the resulting minimum. The proof is obtained by using this triplet in Theorem 9 and bounding h as follows. First note that X × Y is a convex bounded set, since X, Y have the same property. As recalled in Lemma 11, Page 20 of the appendix, convex bounded sets have the so-called uniform interior cone condition. This guarantees that i.i.d. points sampled from a distribution ρ, that has a density bounded away from zero, achieve the following bound on the fill distance: 2d) , with probability at least 1 -δ, when ≥ 0 . Here 0 , C, C are constants that depend only on d, X × Y and the constant c 0 for which the density of ρ is bounded away from zero. The final constants C 1 , C 2 depend also on Q µ,ν .

h ≤ (C -1 log(C /δ)) -1/(
We conclude with the following remark that is useful when we do not know the supports X, Y or we are not able to sample i.i.d. points from the uniform distribution on them.

Remark 10 (Sampling from X × Y using µ ⊗ ν) Since, under Assumption 1, we have µ and ν bounded away from 0, we can compute OT by using i.i.d. samples from ρ = µ ⊗ ν, obtaining the same guarantees as Theorem 2. Indeed, by inspecting the proof above, it is only required that ρ has support X × Y , and has a density that is bounded away from 0. However, the constants 0 , C 1 , C 2 will depend also on how far the density of ρ is bounded away from 0.

Estimators for the kernel mean embeddings of µ, ν

In this section we consider three classes of estimators ( ŵµ , ŵν ) for the kernel mean embeddings w µ ∈ H X and w ν ∈ H Y defined as w µ = φ X (x)dµ(x) and w ν = φ Y (y)dν(y). As we observed in the introduction to Eq. ( 5), the operations we need to perform on ŵµ , ŵν to compute the algorithm are the evaluation of the norm ŵµ 2 H X and the evaluation of ŵµ (x i ) for i ∈ [ ] (and the same for ŵν ). For each class we will specify such operations. Here, we assume that φ X , φ Y are uniformly bounded maps (as for the Sobolev kernel, where

sup x∈X φ X (x) 2 H X = sup x∈X k(x, x) = 1, see Proposition 4).
Clearly a class of estimators must only be chosen if we are able to perform the required operations.

Exact integral. Here we take ŵµ := w µ and ŵν = w ν and we report only the operations for µ since the ones for ν are the same. This is the estimator that leads to the best rates as shown in Theorem 1 and Corollary 3. However it requires to perform the most difficult operations, i.e., (1) ŵµ (

x i ) = w µ , φ X (x i ) H X = X φ X (x i ), φ X (x) dµ(x) = X k X (x i , x)dµ(x) for all i ∈ [ ]; (2) ŵµ 2 H X = w µ , w µ H X = X φ X (x), φ X (x ) H X dµ(x)dµ(x ) = X k X (x, x )dµ(x)dµ(x ).
Moreover the costs C, E in Eq. ( 7) are C = O(1), E = O(1).

Evaluation estimator. Here we assume we are able to evaluate µ(x j ) in a given set of points x 1 , . . . , x nµ with n µ ∈ N (analogously for ν on y 1 , . . . , y nν with n ν ∈ N). We define the estimators as ŵµ = X φ X (x)ĝ µ (x)dx (analogously for ŵν ). Here ĝµ ∈ H X is the kernel least squares estimator [START_REF] Narcowich | Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting[END_REF] of the density µ defined as and(2) ŵµ 2

ĝµ (x) = j∈[nµ] α j k X (x j , x) where α = K -1 X c µ and K X ∈ R nµ×nµ , (K X ) i,j = k X (x i , x j ) for all i, j ∈ [n µ ], while c µ ∈ R nµ , c µ = (µ(x 1 ), . . . , µ(x n )). In this case the steps are: (1) ŵµ (x i ) = X φ X (x i ), φ X (x) ĝµ (x)dx = j∈[nµ] α j X k X (x i , x)k X (x j , x)dx, i ∈ [ ],
H X = i,t∈[nµ] α i α t X X k X (x i , x)k(x, x )k(x t , x )dxdx . Moreover the costs C, E in Eq. (7) are C = O(n 3 µ + n 3 ν ), E = O(n µ + n ν ). Sample estimator.
Here we assume we are able to sample from µ, ν. Let x 1 , . . . , x nµ , n µ ∈ N, be sampled i.i.d. from µ and y 1 , . . . , y nν , n ν ∈ N be sampled i.i.d. from ν. The estimators are ŵµ = 1 nµ j∈[nµ] φ X (x j ), and analogously for ŵν . In this case the operations are:

(1) ŵµ (

x i ) = 1 nµ j∈[nµ] k X (x j , xi ), and (2) ŵµ 2 H X = 1 n 2 µ i,j∈[nµ] k X (x i , x j ). Moreover the costs C, E in Eq. (7) are C = O(n 2 µ + n 2 ν ), E = O(n µ + n ν ).
The effects of the estimators above are studied in Theorem 1 and Corollary 3, reported in Section 1 and proven in Appendix E, Page 22, using standard tools from approximation theory and machine learning with kernel methods [START_REF] Wendland | Scattered Data Approximation[END_REF][START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF][START_REF] Muandet | Kernel Mean Embedding of Distributions: A Review and Beyond[END_REF].

where

|g| H s (X×Y ) := |α|=s D α g L 2 (X×Y ) ≤ g H s (X×Y ) for any g ∈ H s (X × Y ). Let r A (x, y) = φ(x j , ỹj ), Aφ(x j , ỹj ) H XY , we have |f | H s (X×Y ) ≤ |c| H s (X×Y ) + |u1 X | H s (X×Y ) + |v1 Y | H s (X×Y ) + |r A | H s (X×Y ) . Now |c| H s (X×Y ) = 0 since c is the quadratic cost and s ≥ 3 and |u1 X | H s (X×Y ) = |u| H s (X) ≤ u H s (X) , |v1 Y | H s (X×Y ) = |v| H s (Y ) ≤ v H s (Y )
. We recall now that for any two Banach spaces satisfying B ⊆ A, we have w A ≤ w B for any

w ∈ B. Then u H s (X) ≤ u H X , v H s (Y ) ≤ v H Y . Finally |r A | H s (X×Y ) ≤ r A H s (X×Y ) ≤ C 1 Tr(A)
via Lemma 9 page 41 from [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] and Proposition 1 of the same paper, where C 2 depends only on s, X, Y . Now since |f (x, y)| ≤ ε for all (x, y) ∈ X × Y , and since r A (x, y) ≥ 0 for any x, y ∈ X × Y since A is a positive operator, we have

-ε ≤ f (x, y) ≤ f (x, y) + r A (x, y) = c(x, y) -u(x) -v(y), ∀(x, y) ∈ X × Y.
The final result is obtained by defining

C 0 = C 1 (1 + C 2 ).
Proof of Theorem 9. Denote by V (u, v) the functional of Problem 1 and by Vλ 1 ,λ 2 (u, v, A) the functional of Problem 3. Then OT = V (u , v ). Denote by ∆(u, v) the quantity

∆(u, v) = | u, ŵµ -w µ H X + v, ŵν -w ν H Y |. Denote by R2 (u, v) the quantity R2 (u, v) = u 2 H X + v 2 H Y for any u ∈ H X , v ∈ H Y .
Step 0. Admissibility of u , v , A and existence of a maximizer. Note that (u , v , A ) is an admissible point for Problem 3, since the triple satisfies c (x, y) -u (x)v (y) = φ(x, y), A φ(x, y) ∀(x, y) ∈ X × Y , and Problem 3 applies the same constraints but on a subset of X × Y . Moreover λ 1 , λ 2 > 0, this is enough to guarantee the existence of a maximizer for Problem 3. Indeed, as we recall in Lemma 14, Page 21, in the appendix, a form of representer theorem holds for Problem 3 (see Lemma 13, Page 21 in the appendix), moreover the functional is coercive on the finite-dimensional space induced by the data, it has an upper bound and the problem has an admissible point. Now denote by (û, v, Â) a maximizer of Problem 3 and define OT := Vλ 1 ,λ 2 (û, v, Â). By construction OT is the maximum of Vλ 1 ,λ 2 . Then by definition of OT in Eq. ( 6),

OT = OT + λ 1 Tr( Â) + λ 2 R 2 (û, v).
Step 1. Subsampling the inequality. The assumption on h and the fact that the constraints of Eq. (3) satisfy Eq. ( 11) allow to apply Theorem 8, from which there exists a couple (û ε , vε ) that is admissible for Problem 1, where ûε = û -ε 2 and vε = v -ε 2 , for any

ε ≥ C 0 h s-d ( û H X + v H Y + Tr( Â)). In particular, since u H X + v H Y ≤ 2R(û, v), we choose ε = η(R(û, v) + Tr( Â)), with η = C 0 h s-d .
Step 2. Bounding OT -OT. Since OT is the maximum of Problem 1 then

OT = V (u , v ) ≥ V (û ε , vε ) = V (û, v) -ε ≥ OT -∆(û, v) -ε. (12)
Analogously, since (û, v, Â) maximize Problem 3, then

OT = Vλ 1 ,λ 2 (û, v, Â) ≥ Vλ 1 ,λ 2 (u , v , A ) = V (u , v ) -V (u , v ) -Vλ 1 ,λ 2 (u , v , A ) ≥ OT -∆(u , v ) -λ 1 Tr(A ) + λ 2 R 2 (u , v ).
By expressing OT in terms of OT in the inequality above and combining it with Eq. ( 12), we have

λ 1 Tr( Â -A ) + λ 2 (R 2 (û, v) -R 2 (u , v )) -∆(u , v ) ≤ OT -OT ≤ ∆(û, v) + ε. ( 13
)
Step 3. Bound on ∆. Note that for two RKHS H, K and any u, v ∈ H, w, z ∈ K, the following identity holds: u, v H + w, z K = (u, w), (v, z) H⊕K , where H ⊕ K is a RKHS [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. This implies

( u, v H + w, z K ) 2 = (u, w) 2 H⊕K (v, z) 2 H⊕K = ( u 2 H + v 2 K )( w 2 H + z 2 K ).
Applying this inequality to ∆ leads to the following result:

∆(a, b) ≤ R(a, b)γ. ( 14
)
Step 4. Bounding Tr( Â), R(û, v). The bound Eq. ( 13) implies

λ 1 Tr( Â) + λ 2 R 2 (û, v) ≤ λ 1 Tr(A ) + λ 2 R 2 (u , v ) + ∆(u , v ) + ∆(û, v) + ε.
By bounding ∆ via Eq. ( 14), expanding the definition of ε and reordering the terms in the inequality above, we obtain

αTr( Â) + λ 2 R 2 (û, v) -βR(û, v) ≤ λ 1 Tr(A ) + λ 2 R 2 (u , v ) + γR(u , v ),
with α = λ 1 -η and β = γ + η. By completing the square in R 2 (û, v), the inequality above is rewritten as

αTr( Â) + λ 2 (R(û, v) -β 2λ 2 ) 2 ≤ λ 1 Tr(A ) + λ 2 R 2 (u , v ) + γR(u , v ) + β 2 4λ 2 .
Since α ≥ λ 1 /2, by the assumption λ 1 ≥ 2η, the inequality above implies

λ 1 2 Tr( Â) ≤ λ 2 S, R(û, v) ≤ β 2λ 2 + √ S. (15)
with

S := λ 1 λ 2 Tr(A ) + R 2 (u , v ) + γ λ 2 R(u , v ) + β 2 4λ 2 2 .

Conclusion.

From the lower bound in Eq. ( 13) and the bound Eq. ( 14) on ∆, we have that OT

-OT ≥ -λ 1 Tr(A ) -λ 2 R 2 (u , v ) -γR(u , v ) ≥ -λ 2 S.
From the upper bound in Eq. ( 13), the bound for ∆ in Eq. ( 14), the bound for Tr( Â), R(û, v) in Eq. ( 15), the definition of ε, and the fact that λ 1 ≥ 2η, we have ) 2 , since γ ≤ β, and that (a + b + c) 2 ≤ 3a 2 + 3b 2 + 3c 2 for any a, b, c ≥ 0, we have

OT -OT ≤ βR(û, v) + λ 1 2 Tr( Â) ≤ β 2 2λ 2 + β √ S + λ 2 S.
β 2 2λ 2 + β √ S + λ 2 S ≤ 2λ 2 ( β 2λ + √ S) 2 ≤ 2λ 2 R(u , v ) + β λ 2 + λ 1 λ 2 Tr(A ) 2 ≤ 6λ 2 (R 2 (u , v ) + β 2 λ 2 2 + λ 1 λ 2 Tr(A )).

Appendix C. Additional results on convex sets and random points

We first recall that the following property about bounded sets in Euclidean spaces.

Lemma 11 [START_REF] Krieg | Random points are optimal for the approximation of sobolev functions[END_REF]) Let Z ⊂ R q , q ∈ N be a non-empty open set. The following holds:

1. If Z is a bounded Lipschitz domain [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF], then it satisfies the uniform interior cone condition [START_REF] Wendland | Scattered Data Approximation[END_REF].

2. If Z is a convex bounded set, then it is a bounded Lipschitz domain.

Proof For the first point, see Lemma 5 of [START_REF] Krieg | Random points are optimal for the approximation of sobolev functions[END_REF] or Theorem 1.2.2.2 of [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF]. For the second point, see Lemma 4 of [START_REF] Krieg | Random points are optimal for the approximation of sobolev functions[END_REF] or Lemma 7 in [START_REF] Dekel | Whitney estimates for convex domains with applications to multivariate piecewise polynomial approximation[END_REF]. Also, the fact that a convex bounded open set has the uniform interior cone condition is implied by Proposition 11.26 of [START_REF] Wendland | Scattered Data Approximation[END_REF], that proves it for the more general class of sets called star shaped sets w.r.t. a ball.

Lemma 12 (Fill distance of i.i.d points on a u.i.c. set [START_REF] Reznikov | The covering radius of randomly distributed points on a manifold[END_REF])) Let Z ⊂ R q , q ∈ N be a non-empty open set satisfying the uniform interior cone condition (see Theorem 11). Let Z be a collection of points sampled independently and uniformly at random from a probability ρ that admits density (denote it by p) and such that p(z) ≥ c 0 > 0 for any z ∈ Z. Let δ ∈ (0, 1]. Then there exist 0 , C 1 , C 2 depending on c 0 , Z, ρ, q, r 0 such that for ≥ 0 , the following holds with probability at least 1 -δ:

h ≤ (C 1 -1 log(C 2 /δ)) 1/q .
Proof The uniform interior cone condition guarantees that there exists a cone C such that for any z ∈ Z there exists a spherical cone of radius r such that C z ⊆ Z that is congruent to C and with vertex in z. Then, for any r ≤ r 0 there exists c 1 such that

ρ(B(z, r) ∩ Z) ≥ ρ(B(z, r) ∩ C z ) ≥ c 0 vol(B(z, r) ∩ C z ) vol(Z) ≥ c 0 c 1 vol(Z) r q .
Then we can apply Theorem 2.1 of [START_REF] Reznikov | The covering radius of randomly distributed points on a manifold[END_REF] with Φ(r) = c 0 c 1 vol(Z) r q obtaining that there exists 0 , C 1 , C 2 depending on c 0 , c 1 , q, r 0 , Z such that with probability at least

1 -δ, h ≤ (C 1 -1 log(C 2 /δ)) 1/q .
Deriving the dual. Next, let us observe that problem (3) is convex, and admits a feasible point by Corollary 7 and a maximizer by Lemma 14. The same applies to (16). Therefore, strong duality holds. The Lagragian of Eq. ( 16) is

L(u, v, B, γ) = u, ŵµ + v, ŵν -λ 1 TrB -λ 2 u 2 H X -λ 2 v 2 H Y + i=1 γ i (c(x i , ỹi ) -u, φ X (x i ) -v, φ Y (ỹ i ) -Φ i BΦ i ). (17)
At the optimum, we have

∇ u L(u, v, B, γ) = 0 and ∇ v L(u, v, B, γ) = 0 , which yields u = 1 2λ 2 ( ŵµ - i=1 γ i φ X (x i )) v = 1 2λ 2 ( ŵν - i=1 γ i φ Y (ỹ i )). (18)
Let us now derive the optimality condition on B: we have

sup B∈S + (R ) - i=1 γ i Φ i BΦ i -λ 1 TrB = sup B∈S + (R ) B, -( i=1 γ i Φ i Φ i + λ 1 I ) = 0 if i=1 γ i Φ i Φ i + λ 1 I 0 -∞ otherwise. (19)
Plugging Eq. ( 18) and Eq. ( 19) in the Eq. ( 17), we get (5). Finally, using Eq. ( 18), we have OT = û, ŵµ H X + û, ŵν H X = q 2 2λ 2 -1 2λ 2 j=1 γj ( ŵµ (x j ) + ŵν (ỹ j )).

To solve (5), it is possible to use standard software packages [START_REF] Boyd | Convex Optimization[END_REF]. Alternatively, it can be made more scalable by adding a self-concordant barrier term to (5) and using interior point methods with Newton steps. For a given barrier penalization δ > 0, we thus aim to solve

min γ∈R 1 4λ 2 γ Qγ - 1 2λ 2 j=1 γ j z j + q 2 4λ 2 - δ log det( i=1 γ i Φ i Φ i + λ 1 I )) such that j=1 γ j Φ j Φ j + λ 1 I 0. (20)
Starting from an initial value δ 0 , the barrier method [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF] consists in iteratively solving Eq. (20) (using Newton iterations), and progressively decreasing δ. In Theorems 16 and 17, we precisely analyze the complexity of the barrier method applied to Eq. ( 5), and bound the number of operations required to obtain an estimator of OT with a desired accuracy.

Theorem 16 Using a dual interior point method, a solution of problem (3) with value precision O(τ ) can be obtained in O(C + E + 3.5 log( τ )) operations and O( 2 ) memory, where E is the cost of querying ŵµ and ŵν , and C is the cost of computing q 2 . Proof of Theorem 16. Removing terms that are constant in γ, problem ( 20) is equivalent to minimizing the dual functional

J(γ) def = 1 4λ 2 γ Qγ - 1 2λ 2 j=1 γ j z j - δ log det( i=1 γ i Φ i Φ i + λ 1 I ).
Its gradient is

J (γ) i = 1 2λ 2 (Qγ) i - 1 2λ 2 z i - δ Φ i (Φ diag(γ)Φ + λ 1 I ) -1 Φ i , i ∈ [ ],
and its Hessian

J (γ) ij = 1 2λ 2 Q ij - δ [Φ i (Φ diag(γ)Φ + λ 1 I ) -1 Φ j ] 2 , i, j ∈ [ ].
From there, we may minimize J(γ) using damped Newton iterations

γ = γ - [J (γ)] -1 J (γ) 1 + δ λ(γ)
, where λ 2 (γ) = J (γ) [J (γ)] -1 J (γ) is the Newton decrement, or using backtracking linesearch Newton iterations [START_REF] Boyd | Convex Optimization[END_REF].

Number of iterations. J (γ) can be computed and inverted in O( 3 ) operations, and assuming ŵµ (x i ), i ∈ [l] and ŵν (ỹ i ), i ∈ [l] are precomputed, J (γ) can be computed in O( 3 ) operations, hence the complexity per iteration is O( 3 ).

Let F (γ) def = 1 4λ 2 γ Qγ -1 2λ 2
j=1 γ j z j + q 2 4λ 2 be the objective function of (5). Since H(γ) def = -log det( i=1 γ i Φ i Φ i + λ 1 I )) is a self-concordant barrier function of concordance parameter , standard results on barrier methods imply that δ controls the deviation (in value) to the optimum of F [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF]. Moreover, a solution γ to (5) of precision τ > 0, i.e. satisfying F (γ) -F (γ) ≤ τ where γ is the optimum of (5), can computed in O( √ log τ ) Newton iterations using an interior point method by progressively decreasing δ using a suitable scheme until δ ≤ τ : see [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF].

Hence, taking into account the O(E ) cost of computing z j , j = 1, ..., and the O(C) cost required to compute q 2 , to achieve a precision τ > 0 a total of O(C + E + 3.5 log( τ )) operations and O( 2 ) memory are required.

In Theorem 16, we may use any of the kernel mean estimators presented in Section 6 and apply the corresponding computational costs C and E. However, the given bounds only apply to the precision in value, i.e., on V λ 1 ,λ 2 (û, v, Â) -V λ 1 ,λ 2 (u, v, A), and not on the solutions (u, v, A) themselves. In Theorem 17, we derive bounds on the algorithmic approximation of the estimator in (6) obtained by minimizing Eq. ( 20) to a precision τ > 0, as a function of its computational complexity.

Theorem 17 Under the same notation and assumptions as Theorem 9, let γ be obtained by running a barrier method on (5) to precision τ > 0, i.e., by iteratively solving (20) and progressively decreasing δ until δ ≤ τ , as described by [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF]. Define OT as OT = q 2 2λ 2 -1 2λ 2 j=1 γj ( ŵµ (x j ) + ŵν (ỹ j )).

Then OT satisfies the following bound

| OT -OT| ≤ 6λ 2 ( u 2 H X + v 2 H Y ) + 6 β 2 λ 2 + 6λ 1 TrA + 6λ 2 τ.
Further, the considered algorithm to compute OT has a cost of O(C + E + 3.5 log( τ )) in time and O( 2 ) in memory.

Proof of Theorem 17. Let γ be obtained by running a barrier method on Eq. ( 5) to precision τ , i.e. by solving (20) and decreasing δ until δ ≤ τ . Then, from properties of barrier methods [START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF] we can associate to γ the following primal feasible points, obtained by nullifying the Lagragian of (20) at γ:

ũ = 1 2λ 2 ( ŵµ - i=1 γi φ X (x i )) ṽ = 1 2λ 2 ( ŵν - i=1 γi φ Y (ỹ i )) Ã = ij=1 B ij φ(x i , ỹi ) ⊗ φ(x j , ỹj ),
with B = δ ( i=1 γi Φ i Φ i + λ 1 I ) -1 . In particular, from properties of interior point methods [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF][START_REF] Nesterov | Interior-Point Polynomial Algorithms in Convex Programming[END_REF], we have (i) (ũ, ṽ, Ã) is a feasible point of (3), (ii) The duality gap between the objective Vλ 1 ,λ 2 of (3) evaluated at (ũ, ṽ, Ã) and the objective F of (5) at γ is equal to δ, i.e. F (γ) -Vλ 1 ,λ 2 (ũ, ṽ, Ã) = δ.

Further, note that we have OT def = q 2 2λ 2 -1 2λ 2 j=1 γj ( ŵµ (x j ) + ŵν (ỹ j ))

= ũ, ŵµ H X + ṽ, ŵν H Y .

Let us now bound OT-OT. We will follow similar arguments than in the proof of Theorem 9, with an additional τ precision term. As in the proof of Theorem 9, let (ũ ε , ṽε ) = (ũ-ε/2, ṽ-ε).

Since (ũ, ṽ, Ã) satisfy the constraints of (3), from the same arguments (ũ ε , ṽε ) defines a feasible point for (1). Hence, we have the equivalent of Eq. ( 12): OT = V (u , v ) ≥ V (ũ ε , ṽε ) = V (ũ, ṽ) -ε ≥ OT -∆(ũ, ṽ) -ε.

Next, let γ be the optimum of (5), and F the objective function of (5). By strong duality (see Theorem 15), we have F (γ) = Vλ 1 ,λ 2 (û, v, Â). Further, by optimality of γ, we have F (γ) ≤ F (γ). Hence, using (ii) and δ ≤ τ , we have

Vλ 1 ,λ 2 (ũ, ṽ, Ã) = F (γ) -δ ≥ F (γ) -δ = Vλ 1 ,λ 2 (û, v, Â) -δ ≥ Vλ 1 ,λ 2 (u , v , A ) -τ = V (u , v ) -V (u , v ) -Vλ 1 ,λ 2 (u , v , A ) -τ ≥ OT -∆(u , v ) -λ 1 Tr(A ) + λ 2 R 2 (u , v ) -τ.
From there, the rest of the proof of Theorem 9 follows identically: developing Vλ 1 ,λ 2 (ũ, ṽ, Ã) and combining with Eq. ( 21), we have λ 1 Tr( Ã -A ) + λ 2 (R 2 (ũ, ṽ) -R 2 (u , v )) -∆(u , v ) -τ ≤ OT -OT ≤ ∆(ũ, ṽ) + ε.

Hence λ 1 Tr( Ã) + λ 2 R 2 (ũ, ṽ) ≤ λ 1 Tr(A ) + λ 2 R 2 (u , v ) + ∆(u , v ) + ∆(ũ, ṽ) + ε + τ.

Therefore, replacing S from the proof of (9) with S := S + τ , and û, v, Â with ũ, ṽ, Ã, the rest of the proof follows and eventually yields

| OT -OT| ≤ 6λ 2 (R 2 (u , v ) + β 2 λ 2 2 + λ 1 λ 2 TrA + τ ).
To conclude, note that as a consequence of Theorem 16, OT can be computed in O(C + E + 3.5 log( τ )) operations and O( 2 ) memory.

Recovering unregularized optimal transport. In Eq. (3) and in the case of empirical estimators ŵµ and ŵν (see Section 6), we considered the case where the sample pairs (x i , ỹi ), i ∈ [l] covering X × Y are distinct from the samples x i ∼ µ, i ∈ [n µ ] and y j ∼ ν, j ∈ [n ν ]. However, covering X × Y with the n µ n ν pairs given by the µ and ν samples

  that are known explicitly and are easily computable in O(d).

  a) µ, ν have densities. Their supports, resp. X, Y ⊂ R d are convex, bounded and open; b) the densities are finite, bounded away from zero, with Lipschitz derivatives up to order m.

Theorem 2

 2 Let µ, ν satisfy Assumption 1 for m > d and m ≥ 3. Let ∈ N and δ ∈ (0, 1]. Let (x j , ỹj ) be independently sampled from the uniform distribution on X × Y . Let OT be computed with k X = k m+1 , k Y = k m+1 and k XY = k m where k s for s > 0 is the Sobolev kernel in Eq. (8). Then, there exists 0 depending only on d, m, X, Y and C 1 , C 2 depending only on u , v and d, such that when ≥ 0 and λ 1 , λ 2 are chosen to satisfy

Figure 1 :

 1 Figure 1: Left: dual constraint function h(x, y) corresponding to the measures in Figure 2 in Appendix G. Note that h attains its minimum on the graph of transportation map T , and is elsewhere positive. Right: changing parameterization flattens h in the neighborhood of the graph of T . The red dotted line represent the zeros of both functions, and coincides with the graph of T in the original parameterization (left).

  Then | OT-OT| ≤ β 2 2λ 2 +β √ S +λ 2 S.To conclude, by noting that S ≤ λ 1 λ 2 Tr(A )+(R(u , v )+

	β
	2λ 2 2

Proceedings of Machine Learning Research vol 134:1-30, 2021

x -y 2 and under the same assumptions on the densities and their supports, or under smoothness assumptions on the Kantorovich potentials, there exists a positive operator on a suitable Hilbert space that represents the function c(x, y) -u(x) -v(y) for a pair u, v

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.

Appendix A. Proofs of Corollary 6 and Corollary 7

Proof of Corollary 6. Define h(x, y) := c(x, y) -u (x) -v (y), ∀(x, y) ∈ X × Y . By Theorem 5, there exist w 1 , ..., w d ∈ H s (X × Y ) such that h = d i=1 w 2 i . Since H s (X × Y ) ⊆ H XY , then w 1 , . . . , w d ∈ H XY . Hence A = d i=1 w i ⊗ w i ∈ S + (H XY ). Moreover, by the reproducing property of H XY , A satisfies φ(x, y), A φ(x, y) H XY = d i=1 w j , φ(x, y) 2 H XY = d i=1 w(x, y) 2 = h(x, y), for all (x, y) ∈ X × Y . Finally, note that rank A ≤ d by construction.

Proof of Corollary 7. Theorem 3.3 in De [START_REF] De | The Monge-Ampère equation and its link to optimal transportation[END_REF] implies that Kantorovich potentials satisfy u ∈ C m+2,1 (X), v ∈ C m+2,1 (Y ), where C m+2,1 (Z) for an open set Z is the space of real functions over Z that are m + 2-times differentiable, with all the derivatives of order m + 2 that are Lipschitz continuous (Adams and Fournier, 2003, 1.29, Page 10). Since X is convex and bounded, then Lipschitz continuity of all the derivatives of order m + 2 implies that all the weak derivatives up to order m + 3 are in L ∞ (X). Since L ∞ (X) ⊂ L 2 (X), by the boundedness of X, we have C m+2,1 (X) ⊂ H m+3 2 (X). Analogously C m+2,1 (Y ) ⊂ H m+3 2 (Y ). Then u ∈ H m+3 (X), v ∈ H m+3 (Y ), and we can apply Theorem 5 and Corollary 6 with s = m + 1, when m > d, which guarantee the existence of A ∈ S + (H XY ), since H XY = H s (X × Y ).

Appendix B. Proofs of Theorem 8 and Theorem 9

Proof of Theorem 8. Let f (x, y) = c(x, y) -u(x)1 X (y) -v(y)1 Y (x) -φ(x, y), Aφ(x, y) for any x, y ∈ X × Y , where the function 1 X (x) and 1 Y (y) are respectively the constant function 1 over X and over Y . By construction f ∈ H s (X × Y ). Since X, Y are open bounded convex sets, then X × Y has the same property which, in turn, implies the so-called uniform interior cone condition (see Lemma 11, Page 20 in the appendix). Then we can apply Theorem 1 of [START_REF] Narcowich | Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting[END_REF] for which there exist h 0 , C 0 depending only on s, d, X, Y such that for any h ≤ h 0 the following holds

Appendix D. Representer theorem and coercivity for Problem (3)

We first adapt the proof of the representer theorem from [START_REF] Marteau-Ferey | Non-parametric models for non-negative functions[END_REF]. We use it to prove coercivity of the functional. Given the RKHS H X , H Y , H XY and with the same notation of Problem (3) define ĤX = span{ ŵµ , φ X (x 1 ), . . . , φ X (x )}, ĤY = span{ ŵν , φ Y (ỹ 1 ), . . . , φ Y (ỹ )} and ĤXY = span{φ(x 1 , y 1 ), . . . , φ(x , y )}.

Lemma 13 (Representer theorem, Marteau-Ferey et al. (2020)

Note that the components u 2 , v 2 , A 2 , A 3 do not impact the constraints or the functional but are only penalized by the regularizer, indeed

The same reasoning holds for v(x i ) and for u, w µ H X , v, w ν H Y . For A analogously we have that t = A 2 φ(x i , ỹi ) ∈ Ĥ⊥ XY so φ(x i , ỹi ), A 2 φ(x i , ỹi ) H XY = φ(x i , ỹi ), t H XY = 0 and similarly for A 3 , we have φ(x i , ỹi ), A 3 φ(x i , ỹi ) H XY = 0. Let's see what happens to the penalization terms. For the quadratic term, we have -u 2

H X and analogously for v. For the trace term we have Tr(A) = Tr(A 1 ) + Tr(A 3 ). Moreover A ∈ S + (H XY ) implies that A 1 0 and by the Schur complement property

Then if A 3 0 and different from zero we have -Tr(A) < -Tr(A 1 ). If A 2 is different from zero this implies by the Schur complement property that A 3 is a positive definite operator different from zero and so -Tr(A) < -Tr(A 1 ).

Lemma 14 (Problem (3) has a maximizer) When λ 1 , λ 2 > 0 and when there exists an admissible point

From the lemma above, we have that for any admissible point in S \ Ŝ there exists another admissible point in Ŝ that has a strictly larger value. Note moreover that Ŝ is a finite dimensional Hilbert space (with dimension at most 2 ( +1) 2 ) and that (minus) the functional of Problem ( 3) is coercive on it. Note moreover that Eq. ( 3) is bounded from above since 3) has also an admissible point by assumption, then it admits at least one maximizer (see, e.g., Proposition 3.2.1, Page 119 of [START_REF] Bertsekas | Convex Optimization Theory[END_REF].

Appendix E. Proof of Corollary 3 and Theorem 1

Proof of Corollary 3. The result is obtained by plugging in Theorem 2 the bounds on γ = w µ -ŵµ H X + w ν -ŵν H Y . We deal now with the three scenarios.

(Exact integral) In this case, since ŵµ := w µ , ŵν := w ν , then γ = 0.

(Evaluation) We use here the construction used to analyze kernel least squares, e.g., from [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF]; [START_REF] Rosasco | On learning with integral operators[END_REF] (see also [START_REF] Steinwart | Support Vector Machines[END_REF]. We will do the construction for w µ , since the case of w ν is analogous. We recall that φ X : X → H X is uniformly bounded and continuous on X, since we are considering the Sobolev kernel (see Proposition 4). Denote by T ∈ S + (H X ) the operator defined as T = φ X (x) ⊗ φ X (x)dx where dx is the Lebesgue measure. Note that T is trace class, indeed

In particular the equation above implies that f 2 L 2 (X) = T 1/2 f 2 H X for any f ∈ H X . Now, note that by assumption in this Corollary, we have chosen the kernel k X = k m+1 so H X corresponds to the Sobolev space H X = H m+1 (X) (see Proposition 4). Now by Assumption 1, µ has a density that we denote g µ , that is differentiable up to order m and such that all the derivatives of order m are Lipschitz continuous. Since X is convex and bounded, then Lipschitz continuity of all the derivatives of order m implies that all the weak derivatives up to order m + 1 are in L ∞ (X). Since L ∞ (X) ⊂ L 2 (X), by the boundedness of X, we have that all the weak derivatives of g µ up to order m + 1 belong to L 2 (X), i.e. g µ H m+1 (X) < ∞ so g µ ∈ H m+1 (X) = H X . So, by the reproducing property

Now, the estimator ŵµ is defined as ŵµ = X φ X (x)ĝ µ (x), where ĝµ ∈ H X is the Kernel Least Squares estimator (see [START_REF] Narcowich | Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting[END_REF][START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF] of the density of µ that is g µ . With the same reasoning as above, we see that ŵµ = T ĝµ . Now note that

Now T 1/2 2 op = T 1/2 2 op ≤ Tr(T ) < ∞ as we have seen above. Moreover g µ -ĝµ L 2 (X) is controlled by classical results on approximation theory, e.g. Proposition 3.2 of [START_REF] Narcowich | Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting[END_REF] (applied with α = 0 and q = 2). It is possible to apply such results as the set X is convex bounded and so it satisfies the required uniform interior cone condition [START_REF] Wendland | Scattered Data Approximation[END_REF]) (see Lemma 11, Page 20). The result guarantees that there exists two constants C, h 0 depending only on X, such that g µ -ĝµ L 2 (X) ≤ Ch m+1 g µ H X , where h is the fill distance (see Eq. ( 10)) of the sampled n µ points, with respect to X. Now by Lemma 12, Page 20 we have that h ≤ (C n µ log(C n µ /δ)) 1/d with probability at least 1 -δ for some constants C , C depending on d, X. Then, finally for some constant C we have

(Sampling) In this case we have ŵµ = 1 nµ i∈[nµ] φ X (x i ) where x 1 , . . . , x n are independently and identically distributed according to µ. Then ξ i = φ X (x i ) for i ∈ [n ν ] are i.i.d. random vectors and ŵµ = Eξ 1 . Since φ X : X → H X is uniformly bounded on X (see Proposition 4) denote by c that bound. We have that ξ i H X ≤ c almost surely and

Then 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 ŵµ = 1 nµ nµ i=1 ξ i , for which the following holds with probability 1 -δ w

Proof of Theorem 1. This theorem is a direct consequence of Corollary 3. Let ε > 0, , n µ , n ν ∈ N. We denote by f (x) g(x) the fact that there exists two constants 0 < C 1 ≤ C 2 not depending on x such that C 1 g(x) ≤ f ≤ C 2 g(x) (in our particular case this means that the constants will not depend on ε, , n µ , n ν ). For the rest of the proof we will fix ε -2d/(m-d) . Indeed, this choice implies that -(m-d)/2d = O(ε). We recall, moreover, from Eq. ( 7) (that is proven in Appendix F) that the computational time to achieve the estimator is Õ(C + E + 3.5 ) and the required memory is O( 2 ), where E, C are specified for each scenario in Section 6. Now we will quantify the complexity for the three scenarios.

(Exact integral) In this case, by Corollary 3

Note that, since for this scenario C = O(1), E = O(1), the complexity in time is

In space, analogously O( 2 ) = O(ε -4d/(m-d) ).

(Evaluation) In this case we choose n ν n µ O(ε -d/(m+1) ). Indeed, with this choice n

Then the time complexity is

(Sampling) In this case we choose n ν n µ O(ε -2 ). Indeed with this choice n -1/2 µ = O(ε) (analogously for n ν ). Then, by Corollary 3

Then the time complexity is

Also in this case the space complexity is O( 2 ) = O(ε -4d/(m-d) ).

Appendix F. Dual Algorithm and Computational Bounds

In this section, we describe a dual algorithmic procedure to compute Eq. ( 4), and bound the computational complexity and memory footprint it requires to achieve a given precision. We start by deriving Eq. ( 5), the dual formulation of Eq. ( 3), and express Eq. ( 4) as a function of the dual solution γ.

Theorem 15 The dual problem of Eq. (3) is Eq. (5). Further, the estimator OT can be expressed as a function of the solution γ of Eq. (5):

Proof of Theorem 15.

Finite-dimensional representation. Let us start by formulating (3) using a finitedimensional representation for A ∈ S + (H XY ). Following [START_REF] Marteau-Ferey | Non-parametric models for non-negative functions[END_REF], observe that problem (3) needs only to be solved w.r.t. A in the finite-dimensional Hilbert space spanned by {φ(x j , ỹj ) : j ∈ [ ]}. This result is formally proven for our setting in Appendix D. Therefore, it is sufficient to consider positive operators of the form

With this parameterization, we have TrA = Tr(CK) and φ(x i , ỹi ), Aφ(x j , ỹj ) = (KCK) ij , where K = [k XY ((x i , ỹi ), (x j , ỹj )] ij=1 . As [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], we can thus consider the Cholesky decomposition K = RR (or alternatively take the square root R = K 1/2 ), and represent problem (3) in terms of the columns {Φ j : j ∈ [ ]} of R and solve directly for B = RCR :

we may rewrite (5) as a regularized optimal transport problem:

i=1,...,nµ j=1,...,nν

were we reindexed

and where

Hence, ( 22) can be interpreted as a regularized optimal transport problem, where Γ ∈ R nµ×nν plays the role of the transportation plan, and the marginal violations are penalized with maximum mean discrepancy (MMD) terms [START_REF] Gretton | A kernel two-sample test[END_REF]. When λ 1 goes to 0, the SDP constraint becomes a Γ ∈ R nµ×nν + positivity constraint, and when λ 2 goes to 0, the MMD penalization terms enforce the hard constraints Γ1 nν = 1n µ nµ and Γ 1 nµ = 1n ν nν . In particular, when (λ 1 , λ 2 ) → (0, 0), we recover the unregularized OT problem between the empirical measures μ = 1 nµ nµ i=1 δ x i and ν = 1 nν nν j=1 δ y j , which can be formally verified by deriving the dual of (3) with λ 1 and/or λ 2 equal to 0. 1D transportation maps and dual constraint functions. We illustrate the algorithm described in Theorem 16 in a 1D setting in Figures 2 to 5, by representing the inferred transportation map T obtained from û, defined as T = x -∇ x û(x), and the corresponding constraint function ĥ(x, y) = 1 2 x -y 2 -û(x) -v(y). We sample x 1 , ..., x nµ i.i.d. from µ and y 1 , ..., y nν i.i.d. from ν, and use quasi-random samples (x 1 , ỹ1 ), ...(x , ỹ ) from a 2D Sobol sequence [START_REF] Meerovich | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], and illustrate the effect of varying n, the number of µ and ν samples, and , the number of space-filling samples. For k X , k Y and k XY , we use Gaussian kernels k(x, y) = exp(-x-y 2 2σ 2 ) of fixed bandwidth σ 2 = 0.1, and scale the regularization parameters as λ 1 = 1 and λ 2 = 1 √ n .

Appendix G. Numerical Experiments

Convergence of OT to OT. In Figure 6, we compare OT to the sampled optimal transport estimator on two 4D truncated Gaussian distributions µ and ν s.t. the optimal transportation map from one to another is linear. We progressively increase the number of µ and ν samples, averaging on 20 random draws for each number of samples. The number of filling sample pairs (x i , ỹi ) is = 100 + n, where n = n µ = n ν is the number samples from µ and ν. We select the best estimator OT using a grid search on (λ 1 , λ 2 ). As such, this simulation does not provide a method for selecting those parameters, but rather illustrates that a good pair of parameters exists. (left). Full lines correspond to the average mean absolute error (MAE), shaded areas to 25% -75% and 10% -90% MAE quantiles. The parameters λ 1 , λ 2 are selected via a grid search.