Wasserstein distances and the associated problem of optimal transport date back to the work of Gaspard Monge [Mon81] and have since then become important tools in pure and applied mathematics [Vil03,Vil09,San15]. Tools from optimal transport have been successfully employed in machine learning [AGCJ18, PC19, SHB + 18, FCCR18, ACB17, GPC18, JCG18, MMC16, RCP16, SHB + 17, GJB19, SCSJ17, AJJ18, ACB17, DBGGL19, CR12] computer graphics [LCCS18, SdGP + 15, SPKS16, FCVP17], statistics [SC15, AGP18, RW18, WB19, ZP19, PZ19, BGK + 17, CSB + 18, RTC17, CS18, TM18, KTM18, BCP17, dBGLL19, KSS19], and, more recently, computational biology [SST + 19, YDV + 18]. Monge asked the following question: Given two probability measures P, Q in d , how can we transport P to Q while minimizing the total distance traveled by this transport. A classical instantiation of this problem over d is to find a map which is known as the Monge problem, where T # P denotes the push-forward of P under T , that is, T # P (A) = P (T 1 (A)), for all Borel sets A.
(1.2)
The highly non-linear constraint in (1.2) made the mathematical treatment of the Monge problem seem elusive for a long time, until the seminal work of Kantorovich [Kan42,Kan48], who considered the following relaxation. Instead of looking for a map T 0 , look for a transport plan 0 in the set of all possible probablity measures on d ⇥ d whose marginals coincide with P and Q, which we denote by ⇧(P, Q). This leads to the optimization problem min Z kx yk 2 2 d (x, y) s.t. 2 ⇧(P, Q), (1.3) which is known as the Kantorovich problem, and whose value is the square of the 2-Wasserstein distance, W 2 2 (P, Q), between the two probability measures P and Q. The two optimization problems are indeed linked: Brenier's Theorem (Theorem 1 below), guarantees that under regularity assumptions on P , a solution 0 to (1.3) is concentrated on the graph of a map T 0 . That is, using a suggestive informal notation, 0 (x, y) = P (x) y=T 0 (x) , where denotes a point mass. Moreover, T 0 is the gradient of a convex function f 0 . While cost functions other than kx yk 2 2 could be of interest, such as kx yk p 2 for p 1, this work entirely focuses on the quadratic cost, which allows leveraging the well-established theory of convex functions and formulating key assumptions in terms of strong convexity.
Statistical optimal transport describes a body of questions that arise when the measures P and Q are unknown but samples are available. While the question of estimation of various quantities such as W 2 (P, Q), for example, are of central importance, for applications such as domain adaptation and data integration [DKF + 18, CFT14, CFHR17, SDF + 18, PCFH16, CFTR17, FHN + 19, RW19], the main quantity of interest is the transport map T 0 itself since it can be used to push almost every point in the support of P to a point in the support of Q. Moreover, the optimal transport map plays an important role in characterizing the Riemannian geometry that arises from endowing probability measures that have finite second moments with the 2-Wasserstein-distance. In particular, it can be used to define the right-inverse to the exponential map in that space [Gig11], which in turn enables the generalization of PCA (principal component analysis) to spaces of probability measures [BGK + 17, MPZ19]. The goal of this paper is to study the rates of estimation of a smooth transport map T 0 from samples.
To fix a concrete setup assume that we have at our disposal 2n independent observations X 1 , . . . , X n from P and Y 1 , . . . , Y n from Q, based on which we would like to find an estimator T for T 0 . This statistical problem poses several challenges: (i) The most straightforward estimator is obtained by replacing P and Q by their empirical counterparts [MC98]. It leads to a finite-dimensional linear problem that can be approximated very e ciently due to recent algorithmic advances [Cut13, AWR17, PC19, DGK18]. However, even if the resulting optimizer ˆ is actually a map (matching), which it is not in general, it is not defined outside the sample points. In particular, it does not indicate how to transport a point
x / 2 {X 1 , . . . , X n }. In contrast, we would like to obtain an estimator T with guarantees in L 2 (P ), that is, with convergence of k T T 0 k 2 L 2 (P ) := Z k T (x) T 0 (x)k 2 2 dP (x) .
(1.4)
Note that such an estimator of T 0 could be obtained by post-processing the above optimizer ˆ , for example by interpolation techniques, see [ABM16]. We also employ related techniques in Section 6 to obtain practical estimators for T 0 . However, we are not aware of a statistical analysis of these procedures.
(ii) It is known that the estimator W 2 2 ( P , Q) can be a poor proxy for W 2 2 (P, Q) if the underlying dimensionality of the distributions P and Q is large, as it su↵ers from the so-called curse of dimensionality. For example, if both P and Q are absolutely continuous with respect to the Lebesgue measure in d , d 4, it is known that up to logarithmic factors, [|W 2 ( P , Q) W 2 (P, Q)|] ⇣ n 1/d (see [NWR19]). In fact, we show in Theorem 6 that without further assumptions on either P, Q, or T 0 , no estimator can have an expected squared loss (1.4) uniformly better than n 2/d . (iii) While many heuristic approaches have been brought forward to address the previous point, a thorough statistical analysis of the rate of convergence has so far been lacking. This can be partly attributed to the structure (or lack thereof) of problem (1.3). Being a linear optimization problem, it lacks simple stability estimates that are key to establish statistical guarantees by relating k T T 0 k L 2 (P ) to the sub-optimality gap
In this paper, we aim to address these problems by imposing additional assumptions on the transport map T 0 that lead to a rate faster than n 2/d . One assumption we impose on the transport map T 0 is smoothness, a standard way of alleviating the curse of dimensionality in non-parametric estimation. Another key assumption is based on an observation of Ambrosio published in an article by Gigli [Gig11]. They show that the optimization problem (1.1) has quadratic growth, in the sense of a stability estimate kT T 0 k 2 L 2 (P ) . provided T 0 = rf 0 is Lipschitz continuous on d and T # P = Q. While this observation does not immediately lend itself to the analysis of an estimator due to the presence of the push-forward constraint, we show in Proposition 10 that under similar assumptions, the so-called semi-dual problem (see (2.5) below) admits a stability estimate. Due to the rising interest in Optimal Transport as a tool for statistics and machine learning, many empirical regularization techniques have been proposed, ranging from the computationally successful entropic regularization [Cut13, GCPB16, AWR17], `2-regularization [BSR17], smoothness regularization [PCFH16], to regularization techniques specifically adapted to the application of domain adaptation [CFT14,CFTR17,CFHR17]. Notably, [GCPB16] also consider regularization based on the semi-dual objective and reproducing kernel Hilbert spaces. However, the statistical performance of these regularization techniques to estimate transport maps from sampled data has been largely unanswered, with the following exceptions.
The estimation of transport maps has been studied in the one-dimensional case under the name uncoupled regression [RW19] where the sample Y 1 , . . . , Y n is subject to measurement noise. There, the main statistical di culty arises from the presence of this additional noise and boils down to obtaining deconvolution guarantees in the Wasserstein distance. Such guarantees were recently obtained under smoothness assumptions on the underlying density [CCDM11, DM13, DFM15] but they do not translate directly into rates of estimation for the optimal transport map beyond the 1D case. Note that in the presence of Gaussian measurement noise, the rates of estimation are likely to become logarithmic rather than polynomial as deconvolution is a statistically di cult task.
Concurrently to this work, [FLF19] study the estimation of linear transport maps and establish a fast rate of convergence in this parametric setup. Moreover, after finishing the first version of this paper, the authors became aware of the parallel work [Gun18]. There, Gunsilius analyzes the asymptotic variance of an estimator for f 0 under smoothness assumptions on the densities of the marginals P and Q and states the problem of obtaining estimation rates for the transport map T 0 explicitly as an open problem, which we address in this paper. A similarity between [Gun18] and this work is that the rates for the variance are obtained by applying empirical process theory to the semi-dual objective function. However, we note the following key di↵erences: First, Gunsilius obtains curvature estimates for the semi-dual objective via variational techniques, while we use strong convexity of the candidate potentials. Note that under slightly stronger regularity assumptions, i.e., uniform convexity of the supports, his assumptions would imply strong convexity for the ground truth potential, as shown in [Gig11]. Second, by assuming smoothness of the transport map instead of the distributions P and Q, our results are more flexible and can be applied in cases where neither distribution possesses Hölder smooth densities. Third, by appealing to Cafarelli's global regularity theory, [Vil09, Theorem 12.50(iii)], [Caf92b,Caf92a,Caf96], we can obtain a variance rate of n 2↵/(2↵ 2+d) (up to log factors) under assumptions quantitatively comparable to Gunsilius's, while his results imply a sub-optimal rate of n (2↵ 2)/(2↵ 2+d) , see Section E in the appendix.
We note that both in the application of Ca↵arelli's theorem in Appendix E and in the statement of our main result, Theorem 2 below, currently available analytical tools limit the extent to which minimax results can be established. In Appendix E, the lack of uniformity in Cafarelli's global regularity theory prevents us from claiming (near) minimax optimality over Hölder smooth densities, see Remark 36. Similarly, in Theorem 2, (near) minimax optimality is attained by considering fixed marginal supports because of the need for uniformly bounded constants in some of the classical inequalities we employ (for example, Poincaré inequality), see Remark 3. In e↵ect, further strengthening these minimax results poses an interesting open problem involving deep analytical questions.
The rest of this paper is organized as follows. In Section 2, we review some important concepts of optimal transport, mainly duality and Brenier's theorem. These are instrumental in the definition of our estimator, which is postponed to Section 5. Indeed, since the main goal of this paper is to establish minimax rates of convergence for smooth transport maps, we present these rates in Section 3 and prove lower bounds in the following Section 4, since this proof illustrates well the source of the nonstandard exponent in the rates. We then proceed to Section 5 where we define a minimax optimal estimator constructed as follows. First, we define an estimator for the optimal Kantorovich potential as the solution to the empirical counterpart of the semi-dual problem restricted to a class of wavelet expansions. Then, our estimator is defined as the gradient of this potential. We prove that it achieves the near-optimal rate in the same section. In Section 6, we present numerical experiments on synthetic data, introducing two estimators that exploit smoothness of the transport map. The first illustrates that a version of the estimator considered in Section 5 can in fact be implemented, at least in low dimensions. The second is heuristically motivated and based on kernel-smoothing the transport plan between empirical distributions, showing that practical gains in higher dimensions can be achieved for smooth transport maps. Finally, some useful facts from convex analysis (Section A), approximation theory for wavelets (Section B), empirical process theory (Section C), and tools for proving lower bounds (Section D), are appendix. Moreover, the appendix also contains a version of our upper bounds based on smoothness assumptions on the densities instead of the transport map (Section E), the deferred proofs (Section F), additional lemmas (Section G), and more details on the numerical experiments (Section H). , respectively. By abuse of notation, for a di↵erent measure P , we denote the associated Lebesgue norms by k . k L p (P ) . We abbreviate with "a.e." any statement that holds "almost everywhere" with respect to the Lebesgue measure.
For a di↵erentiable one-dimensional function f : ◆ ⌦ ! , we denote its derivative by d dx f . For a function f : d ◆ ⌦ ! , we denote by @ i = @/(@x i ) its weak derivative in the sense of distributions in direction x i , which coincides with the usual (point-wise)
The symbol @f is also used to denote the sub-di↵erential of a convex function f , while we use the symbols rf for the gradient of a function f and Dg for the derivative of a vector-valued function g :
We begin by recalling the Monge and Kantorovich problems given in Section 1. Let P, Q be two Borel probability measures on d with finite second moments.
The Monge (primal) problem is defined as
Note that we can expand the objective in (2.1) as
Since the first two terms above do not depend on , we obtain the equivalent optimization problem
We focus on this equivalent formulation for the rest of the paper because it is more convenient to work with. Problem (2.2) is a linear optimization problem, albeit an infinite-dimensional one. Hence, it is natural to consider its dual problem:
(2.3) f (x) + g(y) hx, yi , P ⌦ Q-a.e, f 2 L 1 (P ), g 2 L 1 (Q). The dual variables f and g are called potentials, and for an optimal pair (f 0 , g 0 ), f 0 is called a Kantorovich potential.
The dual problem (2.3) can be further simplified: Assume we are given a candidate function f in (2.3) above. Then, we can formally solve for the corresponding g to get an optimal g given by the Legendre-Fenchel conjugate (see Section A) of f :
(2.4)
Plugging solution (2.4) back into the optimization problem leads to the so-called semi-dual problem,
where the supremum in (2.4) is interpreted as an essential supremum with respect to P . By transitioning to the semi-dual, we e↵ectively solved for all constraints in (2.3), leaving us with an unconstrained convex problem that is not linear anymore. Under regularity assumptions, a solution to the semi-dual provides a solution to the Monge problem as indicated by the following theorem, which is a cornerstone of modern optimal transport.
Theorem 1 (Brenier's theorem, [KS84, Bre91, RR90]). Assume P is absolutely continuous with respect to the Lebesgue measure and that both P and Q have finite second moments. Then, a unique optimal solution to (2.2) exists and is of the form 0 = (id, T 0 ) # P , where T 0 = rf 0 is the gradient of a convex function f 0 : d ! . In fact, f 0 can be chosen to be a minimizer of the semi-dual objective in (2.5).
Brenier's theorem implies that a solution to the semi-dual problem readily gives an optimal transport map. Our strategy is to minimize an approximation of the semi-dual and establish stability results as well as generalization bounds to conclude that the minimizer to the approximation is close to the minimizer of the original problem.
Let X 1:n = (X 1 , . . . , X n ) and Y 1:n = (Y 1 , . . . , Y n ) be n independent copies of X ⇠ P and Y ⇠ Q = (T 0 ) # P respectively. Furthermore, assume that X 1:n and Y 1:n are mutually independent. Our goal is to estimate T 0 . To that end, we consider the following set of assumptions on P, Q and T 0 . Throughout, we fix a constant M 2. A1 (Source distribution). Let M = M(M ) be the set of all probability measures P with support ⌦ = [0, 1] d that admit a density ⇢ P with respect to the Lebesgue measure such that M 1 ⇢ P (x) M for almost all x 2 ⌦. Assume that the source distribution P is in M. For R > 1 and ↵ > 1, assume that T 0 2 T ↵ = T ↵ (M, R), where
Our main result is the following theorem. It characterizes, up to logarithmic factors, the minimax rate of estimation of an ↵-smooth transport map T 0 2 T ↵ in the setup described above.
where the infimum is taken over all measurable functions T of the data X 1:n = (X 1 , . . . , X n ), Y 1:n = (Y 1 , . . . , Y n ). Moreover, if P 2 M and T 0 2 T ↵ , there exists an estimator T , given in Section 5, that is near minimax optimal. More specifically, there exists an integer n 0 = n 0 (d, ↵, M, R) such that for any n n 0 , it holds,
Remark 3. Assumption A1 can be significantly relaxed with respect to the geometry of ⌦ and the density of P . In fact, the upper bounds are given under more general assumptions in Section 5. Similarly, the assumption A2(iii) that supp(Q) = [0, 1] d can also be relaxed.
However, the constants in the resulting upper bounds exhibit a dependence on the geometry of the supports of both P and Q as well as on the enclosing set ⌦ through functional analytical results used in the proofs. While it may be possible to make this dependence explicit in terms of geometric features of the sets supp(P ), supp(Q) and ⌦-see for example [Jon81, ENT11, TMS + 13] for such estimates under restrictive assumptions-providing a uniform control on these quantities in terms of easily interpretable properties of the sets is beyond the scope of this article. Instead, we chose to present Theorem 2 under these simplified assumptions to make the results more readable.
To discuss the remaining assumptions, we note that the most essential ones to obtain upper bounds are the following: first, the lower bound in A2 (ii), M 1 DT (x), in particular on the support of P , x 2 ⌦. This yields convergence estimates for the optimal transport map as shown in [Gig11], see (1.5), and might be necessary to obtain fast rates for transport map estimation since it provides curvature estimates commonly needed to prove error bounds for M-estimators [VdV00, Chapter 5]. Second, the Sobolev regularity of T 0 is what governs the approximation rates of T 0 by wavelet expansions (see Section 5.5 below) and thus enables fast rates via a bias-variance trade-o↵. All remaining assumptions in A1 and A2, including the existence of extensions of T 0 to a superset ⌦, are of technical nature and serve to give explicit bounds as needed in the proof of the upper bound. While one might be able to relax these assumptions, especially in specific problem instances, we do not pursue this here beyond the more general versions given in B1 and B2 below.
We conjecture that the logarithmic terms appearing in the upper bound are superfluous and arise as an artifact of our proof techniques. We briefly make a qualitative comment on the rate n 2↵ 2↵ 2+d . Note first that it appears from this rate that estimation of transport maps, like the estimation of smooth functions su↵ers from the curse of dimensionality. However, as ↵ ! 1, this curse of dimensionality may be mitigated by extra smoothness with the parametric rate n 1 as a limiting case. Note also that we can formally take the limit ↵ ! 1, which corresponds to the case where no additional smoothness condition holds beyond having a strongly convex Kantorovich potential with Lipschitz gradient. This is essentially the minimal structural condition arising from Brenier's theorem with additional bounds on the derivative of T 0 . In this case, one formally recovers the rate n 2/d and we conjecture that this is the minimax rate of estimation in the context where T 0 is only assumed to be the gradient of a strongly convex function with Lipschitz gradient. If either of these two additional requirements is not fulfilled, our stability results no longer hold. Remark 4. Since the transport map T 0 is the main focus of our results, our assumptions impose smoothness directly on T 0 . In fact, smoothness of T 0 can also be seen as a consequence of smoothness of the source and target distribution using Ca↵arelli's regularity theory [Caf92b,Caf92a,Caf96]. For completeness, we also give a version of our upper bound results under smoothness assumptions on P and Q in Theorem 33, Section E of the Appendix.
Remark 5. By prescribing a known base measure P , such as the uniform distribution on [0, 1] d , and considering Q = T# P , estimation rates for T immediately translate into rates for estimating
, since ( T , T 0 ) # P is a candidate transport plan between Q and Q. Up to log factors, our rates obtained below match those obtained in [WB19] for the estimation of a smooth density on [0, 1] d in the 2-Wasserstein distance.
In this section we begin by proving the lower bound (3.1) as it sheds light on the source of the non-standard exponent 2↵ 2↵ 2+d in the minimax rate. We prove the following theorem.
Theorem 6. Fix ↵ 1. It holds that inf
where the infimum is taken over all measurable functions of X 1 , . . . , X n , Y 1 , . . . , Y n .
Proof. The proof uses standard tools to establish minimax lower bounds, [Tsy09, Theorem 2.5, Lemma 2.9, Theorem 2.2], that we restate in Appendix D for the convenience of the reader as Theorem 30, Lemma 31, and Theorem 32, respectively. It relies on the following construction.
Set P = Unif([0, 1] d ) 2 M, the uniform distribution on the hypercube. For ↵ > 1, let ⇠ : ! be a non-zero function in C 1 ( ) with support contained in [0, 1] such that there exists x 0 2 [0, 1] with ⇠(x 0 ) 6 = 0, d dx ⇠(x 0 ) 6 = 0, for example a bump-function. Define g : d ! by
and note that rg(x 0 , . . . , x 0 ) 6 = 0 and supp(g) = [0, 1] d by the above assumptions on ⇠.
Let m = d✓n 1 2↵ 2+d e be a positive integer where ✓ is a universal constant to be chosen later. We form a regular discretization of the space [0, 1] d by defining the collection of vectors {x
)), for a constant > 0 to be chosen later. Note that supp(g j ) ✓ x (j) + [0, 1/m] d and hence that the supports of the functions {g j } j2[m] d are pairwise disjoint.
Next, let [ 2 N d be a multi-index and observe that the di↵erential operator @ [ applied to
a second-order Taylor expansion yields that g j has uniformly vanishing Hessian: kD 2 g j k L 1 ( d ) ! 0 as m ! 1. In particular, in that case, there exists m 0 such that kD 2 g j (x)k op 1/2 for all x 2 d , m m 0 , j 2 [m] d . If ↵ = 1, the same can be obtained by choosing small enough. By the same reasoning, we can also guarantee
For m d 8, the Varshamov-Gilbert lemma, Lemma 31, guarantees the existence of binary vectors
With this, we define the following collection of Kantorovich potentials:
It is easy to see (Lemma 38) that for any k = 0, . . . , K and m m 0 , r k is a bijection from [0, 1] d to [0, 1] d . Moreover, by Weyl's inequality and the above bound kD 2 g j (x)k op 1/2, for all k,
Similarly, we obtain kr k k L 1 ([0,1] d ) 2 and max (D 2 k (x)) 2 for all x 2 d . Hence, T k := r k 2 T ↵ (M, R) for M > 2 and small enough. We now check the conditions of Theorem 30, where we consider the distance measure
First, observe that for 0
which completes checking the separation condition (i) of Theorem 30.
To check condition (ii) of Theorem 30, recall the Kullback-Leibler (KL) divergence between two measures Q, P such that Q is absolutely continuous with respect to P is defined by
In view of Lemma 38, for any k = 0, . . . , K, the measure Q k = (r k ) # P is supported on [0, 1] d and in particular, it is absolutely continuous with respect to P . By the change of variables formula, it admits the density
Moreover, let X ⇠ P and Y ⇠ Q k be two random variables. It holds
g j where I d denotes the identity matrix in d . Therefore, since the functions g j have disjoint support, we have for all
where l (A) denotes the lth eigenvalue of a matrix A. Since log(1 + z) z z 2 /2 for all z > 0,
where k • k F denotes the Frobenius norm. Thus,
On the one hand, by the divergence theorem and the fact that g j has bounded support,
where @[0, 1] d denotes the boundary of the unit hypercube and v(x) its outward-pointing unit normal vector. On the other hand
The above three displays yield
for ✓ large enough. This completes checking (ii) in Theorem 30 and hence the proof of the first part of the minimax lower bound.
To show the remaining lower bound of 1/n, repeat the same argument as above with the two potentials 0 (x) = kxk 2 2 /2 and 1 (x) = 0 (x) + ( ✓/ p n)g(x) for ✓ chosen to ensure 0 , 1 2 T ↵ , applying Theorem 32 in Appendix D. The separation condition is given by Z
and the KL divergence between the associated probability distributions can be estimated by
for ✓ large enough.
Looking back at this proof, we get a better understanding of the exponent in the minimax rate n
is the minimax rate of estimation of the kth derivative of a -smooth density in L 2 [MG79], the rate that we obtain is formally that of an "antiderivative" (k = 1) of a = ↵ 1-smooth signal in this model. This is due to the fact that, on the one hand, the information structure, measured in terms of Kullback-Leibler divergence, is governed by the derivative of the signal T 0 , i.e., the Hessian of the ↵ + 1-smooth Kantorovich potential, see (4.3), which is ↵ 1-smooth. This follows directly from the Monge-Ampère equation (4.2). On the other hand, we measure the performance of the estimator in terms of the L 2 (P ) distance between T and T 0 , corresponding to the antiderivative of the Hessian. Of course, in the absence of the classical fundamental theorem of calculus in dimension d > 1 for arbitrary maps T : d ! d , the existence of an antiderivative needs to be guaranteed a priori, as in the case of transport maps by assuming T = rf for f : d ! . Similar rates arise in the estimation of the invariant measure of a di↵usion process when smoothness is imposed on the drift [DR07,Str18]. This is not surprising as the drift is the gradient of the logarithm of the density of the invariant measure in an overdamped Langevin process.
Finally, note that the multivariate case is singularly di↵erent from the traditional univariate case where the rate of estimation of linear functionals such as anti-derivatives is known to be parametric regardless of the smoothness of the signal [IH81].
In this section, we give an estimator T that achieves the near-optimal rate (3.2). We present this estimator under the following more general assumptions on the distribution and the geometry of the support of both P and Q = (T 0 ) # P . We also need slightly weaker conditions on the regularity of the transport map (Sobolev instead of Hölder regularity). After stating these weaker assumptions, we present our estimator and restate the main upper bound. Its proof relies on a separate control of approximation error and stochastic error, similar to a standard bias-variance tradeo↵.
Throughout, we fix two constants M 2, > 1.
Let M = M(M ) be the set of all probability measures P whose support ⌦ P ✓ MB 1 is a bounded and connected Lipschitz domain, and that admit a density ⇢ P with respect to the Lebesgue measure such that M 1 ⇢ P (x) M for almost all x 2 ⌦ P . Assume that the measure P 2 M. B2 (Transport map). For any P 2 M with support ⌦ P , let ⌦P denote a convex set with Lipschitz boundary such that ⌦P ✓ MB 1 , and ⌦ P + M 1 B 1 ✓ ⌦P . Let T = T (M ) be the set of all di↵erentiable functions T : ⌦P ! d such that T = rf for some di↵erentiable convex function
These new conditions have two implications. First, they imply regularity of the Kantorovich potential f 0 , where T 0 = rf 0 , and second, they imply some conditions on the pushforward measure Q = (T 0 ) # P that subsume the generalization of Assumption A2(iii). These results are gathered in the following proposition (see Section F.1 for a proof).
Proposition-Definition 7. Assume that P satisfies B1, T 0 satisfies B2 and let X = X (M ) be the set of all twice continuously di↵erentiable functions f : ⌦P ! such that
M for all x 2 ⌦P . Then there exists a Kantorovich potential f 0 2 X (M ) such that T 0 = rf 0 ,
Further, Q = r(f 0 ) # P has a connected and bounded Lipschitz support ⌦ Q ✓ MB 1 and admits a density ⇢ Q with respect to the Lebesgue measure that satisfies
Note that the simplified Assumptions A1 and A2 from Section 3 follow from B1 and B2 in the case
Additionally, the simplified assumptions restrict the class of transport maps to those such that ⌦ Q = [0, 1] d and for which = ↵. Indeed, noting that
, we can fold the two smoothness conditions into one.
To construct an estimator for T 0 , we observe that if we had access to a Kantorovich potential f 0 , then T 0 = rf 0 by Brenier's Theorem, Theorem 1. In turn, f 0 is the minimum of the semidual objective (2.5). Hence, we replace population quantities with sample ones in its definition to obtain an empirical loss function. Moreover, to account for the assumed smoothness of the transport map and to ensure stability of the objective, we constrain our minimization problem to smooth and strongly convex Kantorovich potentials, restricted to a compact superset of the support of P . Then, our estimator is the gradient of the solution to this stochastic optimization problem.
More precisely, for a measurable function f , let us write
where, as in Section 3, X 1:n = (X 1 , . . . , X n ) and Y 1:n = (Y 1 , . . . , Y n ) are n i.i.d. samples from P and Q, respectively, that are mutually independent as well. Recall from Section 2 that the semi-dual objective is defined as S(f
, where f ⇤ denotes the convex conjugate of f . Replacing both P and Q by their empirical counterparts, we obtain the empirical semi-dual,
In order to incorporate smoothness regularization into the minimization of (5.2), we consider the restriction of potentials f to a wavelet expansion of finite degree, a strategy that is frequently used in non-parametric estimation [HKPT98,GN16]. For the purpose of this section, it is enough to think about wavelets as a graded orthogonal basis of L 2 ( d ), leading to nested subspaces
that roughly correspond to increasing frequency ranges of the continuous Fourier transform of a function f 2 L 2 ( d ). Truncated wavelet decompositions yield good approximations for smooth functions and we control their approximation error in Lemma 13. We only consider the span V J ( ⌦P ) of those basis functions of the wavelet expansion whose support has non-trivial intersection with ⌦P . This is a finite-dimensional vector space as long as the elements of the wavelet basis have compact support. The cut-o↵ parameter J is chosen according to the regularity of f 0 in assumption B2, see Section 5.6, or alternatively can be chosen adaptively by a straightforward but technical extension using a penalization scheme [Mas07] that we omit for readability. Alternatively, other selection methods such as Lepski's method [Lep91, Lep92, Lep93] could be used. In order to ensure the necessary regularity and the compact support of the elements of the wavelet basis, we assume throughout that the wavelet basis is given by Daubechies wavelets of su cient order. For a more detailed treatment of wavelets, we refer the reader to Section B.
To ensure stability of the minimizer of the semi-dual with respect to perturbations of the input distributions P and Q, we further restrict the potentials f to mimic the assumptions in Proposition-Definition 7, in particular, we enforce upper and lower bounds on the Hessian D 2 f on ⌦P by demanding f 2 X (2M ).
Combined, both wavelet regularization and strong convexity lead to the set
of candidate potentials, based on which we define the estimators fJ 2 argmin
for the Kantorovich potential and transport map, respectively. Note that since we consider candidate potentials only on the compact set ⌦P , f ⇤ above is defined as
where ◆ ⌦P is the usual indicator function in convex analysis (see Section A).
With this, we can restate the upper bound of Theorem 2.
Theorem 8. Under assumptions B1 and B2, there exists n 0 2 N and J such that for n n 0 ,
where n 0 , C, and J may depend on d, M, R, ⌦ P , ⌦ Q , ⌦P , n 0 may additionally depend on , C may additionally depend on ↵, and J depends on n. The cuto↵ J depends on ↵ if d 3, but in the cases d = 1 and d = 2, J can be chosen independently from ↵.
Remark 9. A few remarks are in order. (i) Similar upper bounds hold with high probability and can be inferred from the proof. (ii) As written, the estimator fJ is not directly implementable since the calculation of f ⇤ involves computing a maximum over a continuous subset of d . However, this limitation can be overcome by a discretization of the space, although this is not practical even in moderate dimensions. We provide such an example implementation in Section 6, along with a more practical estimator for which we give no theoretical upper bounds.
(iii) The numerical experiments in Section 6 suggest that restricting the optimization in (5.4) to X (2M ), while necessary for our proofs, might not be necessary in practice. (iv) The estimator employed in Theorem 8 can be made adaptive to the unknown smoothness parameter ↵ using a standard penalization scheme, see [Mas07]. We omit this straightforward extension and instead focus on establishing minimax rates of estimation. For a more detailed account, we refer the reader to [H 19]. (v) For the sake of readability, we do not explicitly track the dependence on the parameter M .
However, an inspection of the proof yields that the final rate scales like M c 1 d+c 2 for constants c 1 , c 2 1, i.e., there is an exponential dependence on the dimension d. Further, the dependence on R is captured in (5.15) below and amounts to R
2↵ 2+d log(R) in the case d 3. We do not claim an optimal dependence of our rates on these parameters.
In the rest of this section, we present the proof of Theorem 8. We begin by stating our key result, which relates the semi-dual objective to the square of our measure of performance. This result also allows us to employ a fixed-point argument when controlling the risk of our estimator using empirical process theory. Combined with approximation results for truncated wavelet expansions, these lead to a bias-variance tradeo↵ that achieves the minimax lower bound of Theorem 6 up to log factors.
In this section, we leverage the assumed regularity of the optimal transport map to relate the suboptimality gap in the semi-dual objective function S and the L 2 -distance of interest.
Proposition 10. Under assumptions B1 -B2, for all f 2 X (2M ) as defined in Proposition-Definition 7, we have
(5.5)
Proof. It follows from Proposition-Definition 7(ii) and a second-order Taylor expansion that f is of quadratic type [Kol11] around every x 2 ⌦ P :
for all L 2M . It turns out that these conditions are su cient to obtain the desired result. The upper bound in (5.7) is of the form
Since the convex conjugate is order reversing and because the convex conjugate q ⇤ x of the quadratic function q x can be computed explicitly (Lemma 21), we have
The squared distance term vanishes for L = 4M : by the triangle inequality krf 0 (x) rf (x)k 2 4M and since x 2 ⌦ P , it holds that
Together with the fact that Q = (rf 0 ) # P , this yields
Moreover, by strong duality, we have
The above two displays yield S(f ) S(f 0 ) (8M ) 1 krf rf 0 k 2 L 2 (P ) . In the same way, using the lower bound in (5.7), we get that S(f ) S(f 0 ) 2M krf rf 0 k 2 L 2 (P ) , which concludes the proof of (5.5).
It turns out that (5.6) is even easier to prove. Indeed, by Proposition-Definition 7(ii) and Lemma 17, we get that the upper bound in (5.7) is also true for f ⇤ on all of d . In this case, we can simply take L = 2M and get similar results.
There are many ways to leverage strong convexity in order to obtain faster rates of convergence, often known as fixed-point arguments [Mas07,Kol11,GN16]. In this work, we employ van de Geer's "one-shot" localization technique originally introduced in [vdG87] and stated in a form close to our needs in [vdG02].
In light of Proposition 10, the performance of our estimator TJ = rf J defined in (5.4) requires the control of S( fJ ) S(f 0 ), which can be achieved using tools from empirical process theory. To that end, for any f , define
and let fJ 2 F J . We observe that by optimality of fJ for Ŝ,
To proceed, we control the localized empirical process sup
for some fixed ⌧ 2 > 0. More precisely, we prove the following result in Appendix F.2.
Proposition 11. Let assumptions B1 -B2 be fulfilled and define F J as in (5.3). For any ⌧ > 0, define
Then, there exists C 1 = C 1 (d, ⌦ P , ⌦P , ⌦ Q , M) > 0 such that with probability at least 1 exp( t),
where
(5.8)
Equipped with this result, we can apply van de Geer's localization technique. To simplify the presentation, assume that f = fJ 2 argmin f 2F J S(f ) exists. If not, we may repeat the proof with an "-approximate minimizer and let " ! 0. Throughout the proof, we write f = fJ and k
Fix > 0 to be defined later and set
Note that since s 2 [0, 1] and F J is convex, we have fs 2 F J .
On the one hand, fs is localized in the sense that
By Proposition 10 and the triangle inequality respectively, this yields
Therefore, fs 2 F J (⌧ 2 ). For the same reason, we also have that f 2 F J (⌧ 2 ). On the other hand, fs , akin to f , has empirical risk smaller than f . Indeed, by convexity of Ŝ and the fact that f minimizes Ŝ over G, we obtain
Together with Jensen's inequality and Proposition 10 respectively, the above display yields
Next, note for s as in (5.9), we have that kr fs r f k /2 i↵ kr f r f k . Hence
Recalling Proposition 11, we take such that
(5.10) so that we get
In particular, we can check that (5.10) is fulfilled if we choose such that
for a suitable choice of C 2 > 0.
With this, and applying Theorem 10 again, we get that with probability at least 1 e t , it holds
Moreover, integrating the tail with respect to t readily yields by Fubini's theorem that
We have proved the following result.
Proposition 12. Let B1 -B2 hold and define F J as in (5.3). Then, writing
the estimator TJ defined in (5.4) satisfies
(5.11)
Moreover, with probability at least 1 exp( t),
Next, we control the approximation error inf f 2F J krf rf 0 k L 2 (P ) that appears in Proposition 12. In fact, it is su cient to control kr f rf 0 k L 2 (P ) where f = ⇧ J ext f 0 is the truncation of f 0 to its first J wavelet scales after extending f 0 to all d . In light of Theorem 23, we may assume that ext f has the same C -and W ↵,2 -norm as f up to a constant depending on ⌦P .
To control the approximation associated with truncating a wavelet decomposition, we rely on the following lemma for Besov functions.
Lemma 13. Let f 2 B s p,q ( d ) and denote by ⇧ J its projection onto the first scale J wavelet coe cients. That is, if
where j,g k are multi-dimensional Daubechies wavelets and G j the associated index sets as in Section B. Then, for all 1 p, q 1, s 0,
(5.12)
(5.13)
Moreover, for every q 0 2 [1, 1], s 0 > 0, and
In particular:
Proof. Write for the wavelet coe cients of f . The statements (5.12) and (5.13) follow immediately from the wavelet characterization of the Besov norms, (B.2).
To prove the remaining statements, note that for every j, because k . k `p0 k . k `p and k . k q and k . k q 0 are comparable up to a constant due to
kf k B s p,q . Then, plugging this into the wavelet expansion of ⇧ J f f , we obtain
Finally, to obtain the special cases, note that k .
by Theorem 22.
The above lemma together with Proposition-Definition 7 allows us to check that f 2 F J . Indeed, by Weyl's inequality, we have for any x 2 ⌦P that
It follows from Lemma 13 that k⇧
and C 3 = C 3 (d, , ⌦P ) is large enough. This yields min (D 2 ⇧ J ext f 0 (x)) 1/(2M ) and f is strongly convex. Similarly, we get max (D 2 ⇧ J ext f 0 (x)) 2M and hence that f 2 F J . Thus,
where we used Assumption B2 and Lemma 13.
We have thus proved that inf
(5.14)
We are now in a position to complete the proof of Theorem 8. Combining the bounds (5.11) and (5.14), we get
We conclude the proof by optimizing with respect to J. It yields k TJ T 0 k 2 L 2 (P ) .
:
(5.15)
We note that since ↵ > 1, in the first two cases, d 2 {1, 2}, the cut-o↵ J can be picked independently from ↵. Finally, high-probability bounds can be obtained in a similar manner.
In this section, we provide numerical experiments with synthetic data that illustrate how leveraging the smoothness of the underlying transport map can lead to dramatically improved rates. We give two estimators exploiting smoothness: the first, Twav below, is an approximation to the estimator TJ in (5.4), illustrating that (5.4) can be implemented in low dimensions and that this approximation achieves favorable rates in d = 3. The second, Tker below, is a more practical heuristic two-step procedure based on smoothing the optimal matching between the empirical distributions via radial basis functions. We compare these to a baseline estimator given by the optimal transport plan between the empirical distributions.
Additional implementation details and comments on these experiments are provided in Section H of the Appendix.
In order to highlight the benefit of regularization, we consider the following simple estimator based on the optimal transport matrix between the empirical distributions as a baseline. Denote the empirical distributions of P and Q by
respectively, and calculate the optimal transport matrix 2 n⇥n ˆ = argmin
which can be solved exactly by linear optimization toolkits or approximated by entropic regularization [PC19]. An estimated transport function on the observations X i is then obtained by
which corresponds to the conditional mean of the coupling given X i . Note that since we assume that the sample sizes from P and Q are both n, the optimal ˆ is in fact a (rescaled) permutation matrix, leading to a matching ⇡ : [n] ! [n], and hence to Temp (X i ) = Y ⇡(i) .
Because Temp (X i ) above is only defined on the sample points and we do not want to introduce additional bias against the estimator, we consider the following error measure, approximating the L 2 (P ) norm analyzed in Theorem 2:
6.1.2 Wavelet estimator Next, we turn to an approximation of (5.4). Assume that in addition to a superset ⌦P ◆ ⌦ P , we are also given a superset ⌦Q ◆ ⌦ Q , and that both ⌦P and ⌦Q are boxes (hypercubes). We consider all functions originally defined over ⌦P and ⌦Q as given by their samples on grids with resolution N 2 N, We employ the following discretization/approximation schemes:
(i) The restrictions to functions up to wavelet scale J 2 N can be obtained by parametrizing f by the inverse discrete wavelet transform up to order J, which we write as f = W > J J for wavelet coe cients J 2 m J with m J 2 N. For these experiments, we use db4 Daubechies wavelets, i.e., Daubechies wavelets with four vanishing moments. (ii) An approximation to the (continuous) Legendre transform is given by the discrete Legendre transform,
(iii) The gradient of f on grid-points can be calculated by a finite-di↵erence scheme, which we write as r x f . (iv) To obtain values of these approximations at non-grid points, we appeal to linear interpolation, which we write as P x!{X 1 ,...,Xn} f for the interpolated value of f defined on the grid x at the point {X 1 , . . . , X n } 2 d , collected in one vector.
Given i.i.d. observations X 1:n = {X 1 , . . . , X n } and Y 1:n = {Y 1 , . . . , Y n } from P and Q, respectively, we arrive at the following optimization problem as approximation to (5.4):
Note that we dropped all boundedness and convexity constraints that were given by X (2M ) in (5.4). In practice, this can lead to degraded estimation quality of the gradient of W > J J near the boundary of ⌦ P , see Section H.3 in the Appendix, which we remedy by computing the convex envelope of the resulting estimator, yielding an estimator of f 0 that is convex. This envelope can, for example, be calculated by applying the Legendre transform twice, and we set
With this, we denote by T (J) wav the function obtained by linearly interpolating TJ ,
For the purpose of these experiments, we select the wavelet scale Ĵ by an oracle choice, i.e., as the minimizer of an approximation to the population semi-dual (2.5), see Section H.1, while in practice, one would resort to cross-validation methods for this purpose. Finally, we set Twav = T ( Ĵ) wav .
To obtain an error measure for Twav that is easily comparable to MSE n ( Temp ), we consider the empirical L 2 ( P ) norm on the sample points,
Note that the objective function in (6.2) can be calculated in linear time with respect to the underlying grid, that is, in O(N d ), thanks to e cient algorithms for the discrete wavelet transform [Mal99] and the Linear-Time Legendre Transform algorithm [Luc97]. It can be checked that the objective is convex, rendering first-order methods provably convergent. We use the L-BFGS Quasi-Newton method to find ˆ J , even though the objective is not smooth for every J due to the form of the discrete Legendre transform. In practice, we observe that it converges faster than simple (sub-)gradient descent methods.
Since Twav is mainly used to illustrate the practical behavior of a wavelet-based regularization of the semi-dual objective, we do not explicitly analyze the convergence of Twav to the estimator TJ considered in Section 5. We remark, however, that our approximations closely follow the definition of TJ and that the omitted constraints defining the set X (2M ) could be incorporated into the approximation by means of a finite-di↵erence discretization as well, albeit at an additional computational cost.
While the estimator Twav closely follows our theoretical analysis, its applicability is limited to low dimensions by the fact that a discretization grid is needed. As a heuristic alternative, we consider smoothing the assignment obtained by Temp . The idea of smoothing the empirical transport matrix has been previously used in practice, see for example [SST + 19], where regularized linear regression was used as a post-processing step, and kernels have been applied for smoothing potential functions in [GCPB16]. Here, we obtain a smoother version of Temp by smoothing it via kernel-ridge regression [Mur12].
Let H denote a reproducing kernel Hilbert space (RKHS) with associated kernel k(x, y), x, y 2 d , and norm kT k H for T 2 H. Here, we consider the RKHS given by Gaussian radial basis functions, k(x, y) = exp( ⌫ kernel kx yk 2 ), x,y 2 d , ⌫ kernel > 0.
We fit an RKHS function T to the pairs (X i , Ỹi = Temp (X i )) by solving the regularized kernel regression objective We measure its performance by
Similar to Twav , we select the tuning parameters ⌫ kernel and ⌫ ridge by an oracle procedure, picking those parameters from a finite grid that minimize (6.4) on an independent hold-out sample X1:n , and denote the resulting estimator by Tker .
For d 2 N, we consider the following examples of smooth potentials and transport maps:
where for (id), P (1) = Q (1) = Unif([0, 1] d ). For (exp), P (2) = Unif([0, 1] d ), and the target measure is defined as convex on any compact convex subset of d . For the purpose of qualitative comparisons, we consider case (id) in d = 2. In order to determine the quantitative behavior of the estimators, we study both cases for d 2 {3, 10}. The runtime of computing the optimal transport matching via linear programming scales unfavorably with the sample size, a shortcoming that could be remedied by employing recent numerical approximation techniques [ABRNW19]. Similarly, computing the kernel regression (6.3) with o↵-the-shelf methods scales with O(n 3 ). For the sake of this comparison, we simply restrict our experiments on Temp and Tker to sample sizes 10 4 . Likewise, we do not compute Twav in d = 10 due to the large computational cost.
To give a qualitative idea of the considered estimators compared to the baseline Temp , we visualize the first coordinate of the transport map estimators for case (id) with d = 2 and n = 100 observations from P = Q = Unif([0, 1] 2 ) in Figure 1, together with the ground truth transport map (T (1) 0 ) 1 . To depict the coupling Temp induced by the empirical distributions, we employ 1-nearest-neighbor (1-NN) interpolation to obtain a map on [0, 1] 2 . We observe that the wavelet-based regularization in Figure 1c produces a visibly smoother map compared to the unregularized Temp in Figure 1b. The kernel estimator in Figure 1d is even smoother and visually very similar to the ground truth in Figure 1a, due to the possibility of employing a large amount of regularization. Figure 2: Log-log plot of MSE plotted against n, showing the median error over 32 replicates. See Section 6.2 for details.
To obtain a quantitative comparison, for both test cases and d 2 {3, 10}, we compute Temp , Twav (only d = 3), and Tker over 32 replicates and a logarithmically spaced selection of sample sizes n, reporting the median error over the replicates in Figure 2. Here, the dashed lines indicate the result of linear regression on the logarithmically transformed sample sizes and error results for a selected subset of n.
In 3D, for both test cases, the error curves for the standard empirical measure-based estimator Temp roughly follow a n 2/3 rate. This corresponds to the decay of the average `2 cost of optimal matchings between samples from the uniform distribution on the cube [Tal14,Yuk06] (and Gaussian distributions [Led19]), and also matches the n 2/d rate for the convergence of W 2 2 ( P , P ). We are, however, not aware of a generalization of this rate to the error measure MSE in the case of a transport map di↵erent from the identity.
The error curves for the wavelet estimator Twav all follow a similar trend. For low sample sizes, we obtain rates faster than n 0.85 , showing the large statistical benefit of fitting only functions that have wavelet expansions of low order. For large sample sizes, the error flattens out, which can be explained by the numerical approximation errors dominating the statistical ones. This can be readily seen from repeating the experiment with a smaller grid resolution (N = 33), which shows the same trend for smaller values of n, see Section H.2 in the Appendix. Moreover, for all sample sizes we considered, the error curves for Twav all lie below the baseline estimator. The kernel estimator Tker performs even better, attaining rates close to n 1 and yielding consistently better rates than Twav .
We observe that the favorable behavior of Twav suggests that the restriction of candidate potentials to X (2M ) in (5.4) might not be necessary and could possibly be omitted.
In 10D, for both test cases, Temp shows a convergence rate of about n 0.25 , which is slightly better than the expected n 2/d = n 0.2 rate. It is vastly outperformed by the kernel-based estimator that achieves rates better than n 0.6 in both examples.
Further plots showing individual error curves for varying values of the regularization parameters can be found in Section H.4 of the appendix, illustrating the sample size-dependent performance gain achieved by Twav and Tker .
To summarize, in cases where smoothness of the transport map can be assumed, its estimation greatly benefits from smoothness regularization. In particular, these experiments suggest further research on proving error bounds for the kernel estimator Tker under regularity assumptions on T 0 , for which the minimax rates established in this work can serve as a benchmark.
For twice di↵erentiable functions, it is often more convenient to employ the following analytic criterion for strong convexity.
We define the extension of f outside U by f + ◆ U , which by abuse of notation we also denote by f . Note that if f is (strongly) convex on U then its extension outside U is also (strongly) convex. We call dom(f ) = {x 2 d : f (x) < +1} the domain of a convex function.
We now recall two important notions associated with functions f . First, the subdi↵erential of f at x 2 dom(f ) is defined as
As indicated by the following lemma, the subdi↵erential reduces to the gradient for di↵erentiable functions.
Lemma 15 ([HL01, Corollary D.2.1.4]). Let f
Conversely, if @f (x) = {a} consists of only a single element, then f is di↵erentiable at x with gradient rf (x) = a.
Second, the convex conjugate, or Legendre-Fenchel conjugate, is defined for any function f
Lemma 19. If f, g : U ! , f (x) g(x) for all x 2 U , then f ⇤ (y) g ⇤ (y) for all y 2 . Lemma 20. Given two functions f, g : U ! for a compact set U . Abusing notation, denote the convex conjugate of the function f
The next lemma provides an explicit form for the convex conjugate of a quadratic function. It is instrumental in the stability proof for the semi-dual objective function (see Proposition 10).
Lemma 21. Let a > 0, b, t 2 d , c 2 and let U ✓ d be a closed, convex set. Define
Then,
where d 2 denotes the squared distance d 2 (x, U ) = inf y2U kx yk 2 2 .
Proof. Note first that for any y 2 d ,
where
Moreover,
we see that the infimum in (A.2) is achieved by the projection x of (y b)/a onto the closed convex set U + t. Moreover, the value of the objective at x is given by
Together with (A.1) and (A.2), this completes the proof of the lemma.
In this section, we give a brief overview of the di↵erent function spaces used in the paper and how their norms can be related to their wavelet coe cients.
First, we recall the definition of Hölder and Sobolev spaces. Let ⌦ ✓ d be a closed set with non-empty interior and denote by C u (⌦) the set of uniformly continuous functions on ⌦. The Hölder norms for any f 2 C u (⌦) are defined as follows. For any integer k 0 and f that admits continuous derivatives up to order k, define
and for any real number ↵ > 0, define
The space C ↵ (⌦) is then defined as the set of functions for which this norm is finite. For a vectorvalued function T : ⌦ ! d , T = (T 1 , . . . , T d ) > , we similarly define the norms as the sum over the individual norms,
Similarly, for an integer k 0 and p 2 [1, 1], the Sobolev norms are defined as
where the derivative @ [ is to be understood in the sense of distributions and the Sobolev space W m,2 (⌦) is the space of all functions for which this norm is finite. This definition can be extended to ↵ > 0, for example by defining W ↵,2 (⌦) as the Besov space B ↵ 2,2 (⌦), which we define shortly. Next, we define wavelet bases and Besov spaces, following the definitions given in [Tri06, Section 3], and we refer the reader to this reference for further details on wavelets. Denote by M 2 C r ( ) and F 2 C r ( ) a compactly supported wavelet and scaling function, respectively, for example Daubechies wavelets. This implies that
To obtain a basis of L 2 ( d ), for j 2 N, set
This gives the orthonormal basis
Wavelet coe cients, defined as the expansion coe cients with respect to the above basis for L 2 ( d ) functions, can be used to characterize the so-called Besov spaces: Let 1 p, q 1, s 0, and let the regularity of the above wavelets satisfy
With this, if f admits the wavelet representation
we define the family of norms k . k B s p,q by
In particular, by the orthonormality of the wavelets { j,g k } j,g,k ,
For a bounded Lipschitz domain ⌦ ✓ d , we can define the Besov spaces B s p,q (⌦) by restrictions of functions on d with norm
Note that since we work with compactly supported wavelets, for a function f with compact support, only a finite number of wavelet coe cients are non-vanishing in (B.1). In particular, for non-zero wavelet coe cients are contained in a set ⇤(j) with |⇤(j)| . 2 jd .
The following theorem collects some basic properties of Besov spaces and their relationship to Hölder and Sobolev spaces.
Theorem 22. Let ⌦ ✓ d be a bounded Lipschitz domain. (i) [GN16, Proposition 4.3.6] Let s, s 0 > 0, 1 p, p 0 , q, q 0 1. Then, the following inclusions hold in the sense of continuous embeddings:
Note that the proof in [GN16] of (B.3) is given in 1D, but extends naturally to arbitrary dimensions by summability of the wavelet coe cients.
A very useful tool to handle function spaces on domains is the availability of extension operators that preserve the norms. The following theorem guarantees the existence of such an extension operator that is the same among all Besov spaces. This allows us to characterize Besov functions on domains via the wavelet coe cients of their extensions.
Theorem 23 (Extension operator, [Tri06, Theorem 1.105], [Ryc99]). Let ⌦ ✓ d be a bounded Lipschitz domain. Then, there exists a linear extension operator ext that preserves L 2 -, C -, and W ↵,2 -norms. That is, there exists an extension operator ext = ext(⌦) such that for A 2
We conclude this section by a lemma that provides uniform control of a function by its wavelet coe cients. It is useful to control bracketing entropy numbers.
Lemma 24. Let f 2 V J ( d ) with wavelet coe cients j,g k for compactly supported mother and father wavelets. Then,
Proof. Let x 2 d and write
where we used Hölder's inequality and the fact that only a finite number of k enters the summation for each fixed x 2 d . The last inequality is trivial.
Here, we collect some basic results about empirical processes that are needed in the proofs. Note that because all suprema we deal with are over subsets of finite-dimensional vector spaces, we do not consider issues of measurability in the remainder.
Denote the bracketing number of a set F with respect to a norm k.k by N [ ] ( , F, k . k) and define the Dudley integral as
Theorem 25 (Bernstein chaining, [vW07, Lemma 3.4.2] ). Let F be a class of measurable functions such that [f 2 ] < 2 and kf k 1 M for all f 2 F. Then,
Theorem 26 (Concentration, see [Mas07, Equation (5.50)], going back to [Bou02]). Under the same assumptions as Theorem 25,
Lemma 27 (Bracketing numbers from L 1 covering numbers). Let P be a probability measure and A ✓ F be a set of functions with N (A, L 1 (P ), /2) ( ). Then, N [] (A, L 2 (P ), ) ( ).
Moreover, if for every function in the original class f 2 F, [f 2 ] < 2 and kf k L 1 (P ) < M, then every bracket [f 1 , f 2 ] in a -cover above (note that f 1 and f 2 need not be members of F) satisfies
Proof. Denote by {f 1 , . . . , f N } the centers of a minimal /2-covering of F in L 1 (P ). Let i 2 {1, . . . , N}. Then, each /2 ball around f i is contained in the bracket [f i /2, f i + /2]. Moreover, the L 2 (P ) diameter of the above bracket is bounded by its L 1 (P ) diameter, which is , so the collection of those brackets yields the desired covering with brackets.
The rest of the lemma follows from
The following lemma can be shown by directly specifying a grid or by a volume argument such as [Ver18, Proposition 4.2.12].
For the convenience of the reader, in this section we restate the standard tools we use in the proof of Theorem 6 to establish lower bounds based on Fano's inequality and the Varshamov-Gilbert Lemma, taken from [Tsy09, Theorem 2.9, Lemma 2.9, Theorem 2.2].
Theorem 30 (Lower bounds from multiple hypotheses). [Tsy09, Theorem 2.9] Let K 2, ⇥ = {T 0 , . . . , T K } a collection of hypotheses, and let d be a pseudometric, i.e., a bi-variate function on ⇥ such that
for all 0 j, k, `2 {0, . . . , K}.
Suppose ⇥ fulfills:
(i) d(T j , T k ) 2s > 0, for all 0 j < k K;
(ii) P j ⌧ P 0 , for all j 2 [K], and
with 0 < C < 1/8, where P j = P T j denotes a probability distribution associated with every T j for all j = 0, 1, . . . , K.
Then, inf T sup
Lemma 31 (Varshamov-Gilbert lemma). [Tsy09, Lemma 2.9] Let D 8. There exists a subset ⌧ (0) , . . . , ⌧ (K) of {0, 1} D such that ⌧ (0) = (0, . . . , 0),
, for all 0 j < k K, and
Theorem 32 (Lower bounds from two hypotheses). [Tsy09, Theorem 2.2] Let T 0 , T 1 be two hypotheses with associated probability measures P j = P T j for j 2 {0, 1}. Denoting s = d(T 0 , T 1 )/2, if
In this section, we show how to apply Theorem 8 under smoothness assumptions on the source and target distribution instead of the transport map. This is enabled by Ca↵arelli's regularity theory [Caf92b,Caf92a,Caf96], Theorem 33 below, which gives regularity estimates on the transport map T 0 under regularity assumptions on the source and target densities and their supports.
For simplicity, we denote an open, convex set ⌦ with C 2 boundary as uniformly convex if it can be written as the sublevel set, ⌦ = {f < 0}, of a strongly convex function f . Note that uniform convexity can also be characterized by the positivity of the second fundamental form of ⌦ on its boundary @⌦ [Vil09]. C2 (Densities). With ⌦ P , ⌦ Q as in Assumption C1, assume that ⇢ P 2 C ↵ 1 (⌦ P ) and ⇢ Q 2 C ↵ 1 (⌦ Q ) are two probability densities on d with supp(⇢ P ) = ⌦ P , supp(⇢ Q ) = ⌦ Q that are bounded from above and below on their support. Further, denote the associated probability distributions by P and Q, respectively.
Theorem 33 (Ca↵arelli's global regularity theory). [Vil09, Theorem 12.50(iii)], [Caf96]. Let ↵ > 1 and assume that Assumptions C1 and C2 hold. Then, for the optimal transport potential f 0 , i.e., the solution to (2.5) for P and Q, which is unique up to an additive constant, it holds that f 0 2 C ↵+1 (⌦ P ). Corollary 35. Let ↵ > 1. Under Assumptions C1 and C2, writing f 0 and T 0 = rf 0 for the Kantorovich potential and optimal transport map between P and Q, respectively, there exists an estimator T such that
where X 1:n , Y 1:n denote i.i.d. observations from P and Q, respectively.
Proof. The statement follows by verifying the assumptions B1 and B2 on the original transport map and concluding by Theorem 8, using the same estimator.
Assumption B1 is satisfied for ⌦ P by the (stronger) conditions imposed on it in Theorem 33.
To show Assumption B2, denote by f 0 the optimal transport potential and by T 0 the optimal transport map for the problem given by ⇢ P and ⇢ Q . By Theorem 33, f 0 2 C ↵+1 (⌦ P ) and thus T 0 = rf 0 2 C ↵ (⌦ P , d ). Denoting by f0 the extension of f 0 to ⌦P = (⌦ P ) ✏ from Lemma 34, we observe that ⌦P is connected and has a Lipschitz boundary. This follows from the general fact that the Minkowski sum of a convex set with C 1 boundary and another convex set has a C 1 boundary [KP91]. From the regularity of T0 = f0 and the compactness of ⌦P , all requirements in B2 can now be easily checked.
An inspection of the proof of our lower bound results, Theorem 6, readily yields a lower bound under smoothness assumptions on P and Q as well. In particular, after changing the domain from [0, 1] d to the unit ball in d dimensions, the explicit form of the density of Q k in (4.2) guarantees that the candidate densities fulfill assumptions C1 and C2. This indicates that the above rate is optimal up to log factors.
Remark 36. In Corollary 35, we do not claim any uniform bound over a regularity class, such as an explicit dependence of the constant on the C ↵ 1 -norms of ⇢ P and ⇢ Q alone. Although the statement of Theorem 33 strongly suggests a bound on kf 0 k C ↵+1 (⌦ P ) in terms of k⇢ P k C ↵ 1 (⌦ P ) and k⇢ Q k C ↵ 1 (⌦ Q ) , to the best of our knowledge, such bounds are not available in the literature for the "global" version of Ca↵arelli's regularity theory as stated Theorem 33. On the other hand, such bounds hold "locally", that is, kf 0 k C ↵+1 (U ) can be bounded for strictly open subsets of ⌦ P , see [Vil09,Theorem 12.56], [Caf92b]. This deficiency is due to the non-constructive nature of the available proofs of Theorem 33 and could possibly be remedied, but we consider it out of the scope of this paper. We proceed in order of statement of the results. (i) It follows from B2 that there exists f 0 such that rf 0 = T 0 and by B2(i) that |rf 0 (x)| M for all x 2 ⌦P . Moreover, since f 0 is defined up an additive constant, assume that f 0 (x 0 ) = 0 for some x 0 2 ⌦P . A first-order Taylor expansion yields that for any x 2 ⌦P ,
where in the second inequality, we used B2(i), and in the third one, we used that diam( ⌦P ) 2M according to B2.
(ii) follows immediately from B2(ii). Next, as in (i), (5.1) follows from a first-order expansion: For any x 2 ⌦P , note that on the one hand
Together with the fact that T 0 = rf 0 , this allows us to shift the smoothness index by one speaking about potentials.
To show the statements about Q, by Lemma 14 and (ii) above, f 0 is strongly convex on int( ⌦P ) and can be extended to +1 outside of ⌦P and thus be also be considered a strongly convex function on d . By Lemmas 17 and 18, we conclude that T 0 = rf 0 is a bijection from ⌦ P onto its image ⌦ Q = supp((rf 0 ) # P ), with both rf 0 and (rf 0 ) 1 being continuously di↵erentiable; in other words, rf 0 is a C 1 -di↵eomorphism between ⌦ P and ⌦ Q . Since C 1 -di↵eomorphisms preserve Lipschitz domains [HMT07, Theorem 4.1] and connectedness, we can conclude that ⌦ Q is a connected and bounded Lipschitz domain, and the fact that ⌦ Q ✓ MB 1 follows from T (x) M for all x 2 ⌦P .
Finally, we turn to check the condition on the density ⇢ Q . To that end, note that by the change of variables formula, Q has the density
This readily yields the desired bound in light of B1 which assumes boundedness of ⇢ P and (ii) which gives boundedness of the Hessian D 2 f 0 .
We first bound the expectation of the supremum of interest and then obtain the high-probability bound via concentration. To begin, note that
We first focus on the T 2 -term. Once understood, the T 1 -term can be bounded similarly.
Bound on T 2 -term. We estimate T 2 from above by two suprema, corresponding to low and high frequencies. To this end, we center all functions and consider the extension of a restriction of the functions to ⌦ Q , which allows us to use wavelet expansions for harmonic analysis.
First, because ( Q Q)(g + c) = ( Q Q)g for every function g and constant c 2 , we may assume without loss of generality that Z
Note that f ⇤ and f ⇤ 0 are defined over all of d but we do not control their norms over the whole space. To overcome this limitation, with a slight abuse of notation, we denote by ext f ⇤ (resp. ext f ⇤ 0 ) the Lipschitz extension of the restriction of f ⇤ (resp. f ⇤ 0 ) to ⌦ Q . The existence of a linear extension operator is guaranteed by Theorem 23. In particular, we control the norms of ext f ⇤ and ext f ⇤ 0 . For a function g with wavelet expansion g
as defined in Section B, we define the two L 2 -projections
With this, we write T 2 = T 2,1 + T 2,2 , where
Note that the projected functions are again well-defined continuous functions: By assumption, all f and f 0 are strongly convex, hence their conjugates have Lipschitz-continuous gradients and therefore bounded B 2 1,1 (⌦ Q )-norm. In turn, their projections also have bounded B 2 1,1 (⌦ Q )-norm, and by Theorem 23, this implies that both ⇧ J ext(f ⇤ f ⇤ 0 ) and ⇧ >J ext(f ⇤ f ⇤ 0 ) are in C s ( d ), for s < 2, and hence they are continuous.
Bound on T 2,1 -term.
Recall that it follows from Proposition-Definition 7 that ⌦ Q is a connected Lipschitz domain, so we can apply the Poincaré-Wirtinger inequality, Lemma 39, together with (5.6) from Proposition 10, to get
where we used that we assumed
⌧ , and therefore due to the properties of the extension operator ext,
for some constant C 4 = C 4 (⌦ Q , M). Since ⇧ J is a non-expansive operator on Besov spaces, it follows from the above display that
Bounding the empirical process over this standard function class can now be performed as follows. Observe first that for any function h 2 V J ( d ) with wavelet decomposition h
Next, by symmetrization, for independent copies of Rademacher random variables " i , i = 1, . . . , n,
where both suprema are taken over the set
To control the Rademacher process, fix h 2 G J with wavelet decomposition (F.3). By the Cauchy-Schwarz inequality,
where here and below, all sums without indices are over {0 j J, g 2 G j , k 2 Z d }. Since the right-hand side in the display above does not depend on h, we get by Jensen's inequality that sup 1
Since Y i 2 ⌦ Q , a compact set by assumption, for given j and g, j,g k (Y i ) is non-zero only for k 2 ⇤(j) where ⇤(j) depends on the diameter of ⌦ Q , and |⇤(j)| . 2 jd . Together with the independence of the " i , this yields
By Proposition-Definition 7 and the fact that the j,g k form an orthonormal basis in
We have proved that
(F.4)
Bound on T 2,2 -term.
To control the term T 2,2 , we use a chaining bound for bracketing entropy (Theorem 25). To that end, we exhibit bounds on the L 1 -covering numbers of the corresponding function space, which in turn implies control of L 2 -bracketing numbers. The idea is to bound the covering numbers in the original space F J (⌧ 2 ) and exploit continuity properties of the transformations that lead to the function class in the definition of T 2,2 , in particular the operation of taking the convex conjugate.
Define the function space
and let f 1 , f 2 2 F J (⌧ 2 ) have wavelet coe cients given by the sequences 1 and 2 . Observe that by linearity of the projection and extension operators, we have for any cuto↵ J 0 0,
To control the first term in the right-hand side above, in the following lemma, we establish bounds on the potentials in FJ (⌧ 2 ). Its proof is deferred to the next section.
Lemma 37. There exists a constant C 5 = C 5 (⌦ P , ⌦P , ⌦ Q , M) such that for all f 2 X (2M ),
Moreover, under assumptions B1 -B2, k ext f ⇤ 0 k B 2 1,1 ( d ) C 5 /2 as well.
It follows from Lemma 37 that there exists C 5 such that for any f 2 F(⌧ 2 ), we have
Therefore, by Lemma 13, we get
if we choose J 0 = dlog (C 5 /") /2e so that C 5 2 2J 0 ". To control the second term, we get from Lemma 24 that
Moreover, using respectively the fact that ext is a Lipschitz operator on Besov spaces, ⌦ Q is bounded, and the convex conjugate is non-expansive in L 1 by Lemma 20, we have
where in the last inequality, we used Lemma 24. We have proved that
for some constant C 6 = C 6 ( ⌦P , ⌦ Q ). This inequality allows us to control the L 1 -bracketing numbers of FJ (⌧ 2 ) using `1-covering numbers for the wavelet coe cients. To control the latter, note that for all f 2 F J (⌧ 2 ) ⇢ X (2M ), it holds kf k B 2 1,1 ( d ) . M 2 so that k k `1 . M 2 . Moreover, these wavelet coe cients are in a space of dimension at most C 2 Jd , C( ⌦P , M) > 0 because F J (⌧ 2 ) ✓ V J ( ⌦P ). Hence, choosing " = /4, Lemmas 27, 28, and the previous display yield log
To apply the chaining bound of Theorem 25, note that by (F.2) and Lemma 37, respectively, combined with Lemma 13, we have
is the Dudley integral defined in (C.1). Moreover, by (F.5) and Lemma 29, we have
and therefore,
Together with (F.4), we have T 2 . J (⌧ 2 ) with J defined as in (5.8) and we absorbed r J for d 2 into the first term on the right-hand side.
Bounding T 1 . T 1 can be bounded completely analogously to how we bounded T 2 , with the exception that Lemma 20 is not needed. Thus, we obtain T 1 . J (⌧ 2 ). To obtain a bound that holds with high probability, we apply the concentration result of Theorem 26. For this, note that (F.1) can be written as Ŝ0 (f ) S 0 (f ) = (( P ⌦ Q) (P ⌦ Q))((f f 0 ) ⌦ (f ⇤ f ⇤ 0 )),
where P ⌦ Q denotes the product measure and (f ⌦ g)(x, y) = f (x) + g(y). Following the same argument as in the proof of Lemma 37, we get
Moreover, similarly to Proposition 10, we have for any f 2 F J (⌧ 2 ) that
We can therefore apply Theorem 26 and conclude that with probability at least 1 e t , it holds that optimization at all scales of J with the all zeros vector. An approximation to the ground truth semidual objective (2.5) is obtained by integrating both W > J J and the interpolation P y!T 0 (x) L(W > J J ) over [0, 1] d with Simpson's rule, exploiting the fact that integration over Q is equivalent to integration over P under the push-forward T 0 .
The parameters ⌫ kernel and ⌫ ridge are chosen via an oracle procedure, picking the best ⌫ kernel 2 {10 9 , 10 8.5 , . . . , 10 5 }, ⌫ ridge 2 {10 5 , 10 4.5 , . . . , 10 1 } as determined by evaluating MSE for an independent draw from P .
The data for all quantitative plots are obtained by taking the median over 32 i.i.d. replicates. Error bars are not shown since 95-percentile bootstrapped confidence intervals were not visible for most estimators at the present scale of the plots.
Running on one core of server processors such as an Intel R Xeon R E5-2670 v3 (2.30GHz), the calculation of all wavelet scales of Twav for one replicate takes between 10 and 70 minutes, depending on the sample size and hence the conditioning of the problem. We note that the runtime and space complexity of the algorithm is determined both by the discretization size N and to a lesser extent by the sample size n, as opposed to computing the optimal transport plan between empirical distributions, whose complexity in the regimes considered here is governed entirely by the sample size n. As a comparison, computing all 100 di↵erent parameter settings for Tker with n = 10000 takes about five hours.
To investigate the observed flattening out of the Twav error curves in Figure 2, we repeat the experiment for (id) and d = 3 with a lower resolution discretization, N = 33. In Figure 3, we compare the resulting error to the N = 65 case considered in Section 6. The error bottoms out for much lower values of n, suggesting that numerical accuracy is indeed responsible for this behavior. 4: Estimation errors for (id), d = 3, direct output of (6.2), Twav,direct , compared to its convex envelope Twav . Median over 32 replicates. For low sample sizes, taking the convex envelope improves the estimation rate.
Write fJ = W > J ˆ J . We observe that while fJ by itself often qualitatively yields good results on ⌦ P (the support of P ), the lack of additional global regularity conditions in the optimization (6.2) often leads to poor approximation outside of the support, i.e., on ⌦P \ ⌦ P . In particular, certain regions could simply not enter the optimization at all due to the form of the discrete Legendre transform L in (6.1). The numerical gradient estimate is somewhat sensitive to this, especially near the boundary, which prompts us to regularize the result further by instead considering its convex envelope, i.e., the largest convex function that lies entirely below the graph of fJ .
We give a qualitative visual example of this in Figure 5. Depicted for d = 2 are first the ground truth potential f 0 (x) = 1 2 kxk 2 2 corresponding to the identity transport map (Figure 5a) and its derivative in the first coordinate (Figure 5b), where we set P = Unif([0, 1] 2 ). The parts of those functions corresponding to the support of P are depicted in yellow, while those outside of the support are colored blue. Results of the optimization for n = 100 i.i.d. samples from P and Q = Unif([0, 1] 2 ) are shown in Figures 5c (potential) and 5d (first coordinate of numerical gradient) for J = 1. While the estimator yields a good approximation to f 0 and T 0 on the interior of ⌦ P , the outside appears very ragged. This is remedied by applying L twice, see Figures 5e (potential) and 5f (numerical gradient).
We also quantitatively compare the estimation accuracy of a linear interpolation of r x f Ĵ , denoted by Twav, direct to that of Twav by repeating the experiment for (id) and d = 3, plotting the results in Figure 4. One can observe that for low sample sizes, considering the convex envelope indeed significantly improves the estimation accuracy.
In order to further investigate the performance of the proposed regularized estimators Twav and Tker , we show error plots for a selection of fixed regularization parameters in Figure 6. In the case of Twav , we plot T (J) wav for the whole range of J 2 {0, . . . , 3} that can be calculated by means of the discrete wavelet transform when N = 65. In the case of Tker , we plot a subset of all considered combinations ⌫ kernel 2 {10 9 , 10 8.5 , . . . , 10 5 }, ⌫ ridge 2 {10 5 , 10 4.5 , . . . , 10 1 }. First, for the wavelet estimators T (J) wav , we observe that the best bias-variance trade-o↵ is achieved by J = 0 for smaller values of n and by J = 1 for n 10 3.5 . Moreover, the curve for J = 3, corresponding to no regularization, roughly follows the shape of the error curve of Temp , until it flattens out due to the numerical error. This is most likely due to the grid approximations involved in the definition of Twav .
Next, the kernel estimators T (⌫ ridge ,⌫ kernel ) ker in 3D show mostly similar error curves for the range of parameters considered, although a trade-o↵ that depends on the sample size can be observed between (⌫ ridge , ⌫ kernel ) 2 {(10 7 , 10 4 ), (10 7 , 10 3 )}. Moreover, in the case of very strong regularization, (⌫ ridge , ⌫ kernel ) = (10 5 , 10 5 ), the error curve flattens out for values of n as low as 10 2.5 . In general, for Tker , we observe that its performance is more sensitive to changes in ⌫ kernel than to those in ⌫ ridge .
Last, In 10D, we observe a wider range of error curves according to the regularization strength, ranging from a curve that closely matches that of Temp for (⌫ ridge , ⌫ kernel ) = (10 9
, 10 1 ) to a curve that matches most of the error curve shown in Figure 2d for (⌫ ridge , ⌫ kernel ) = (10 7 , 10 4 ).