Rates of Estimation of Optimal Transport Maps using Plug-in Estimators via Barycentric Projections
Year: 2021

Abstract:
Optimal transport maps between two probability distributions μ and ν on R have found extensive applications in both machine learning and statistics. In practice, these maps need to be estimated from data sampled according to μ and ν. Plugin estimators are perhaps most popular in estimating transport maps in the field of computational optimal transport. In this paper, we provide a comprehensive analysis of the rates of convergences for general plug-in estimators defined via barycentric projections. Our main contribution is a new stability estimate for barycentric projections which proceeds under minimal smoothness assumptions and can be used to analyze general plug-in estimators. We illustrate the usefulness of this stability estimate by first providing rates of convergence for the natural discretediscrete and semi-discrete estimators of optimal transport maps. We then use the same stability estimate to show that, under additional smoothness assumptions of Sobolev type or Besov type, kernel smoothed or wavelet based plug-in estimators respectively speed up the rates of convergence and significantly mitigate the curse of dimensionality suffered by the natural discrete-discrete/semi-discrete estimators. As a by-product of our analysis, we also obtain faster rates of convergence for plug-in estimators of W2(μ, ν), the Wasserstein distance between μ and ν, under the aforementioned smoothness assumptions, thereby complementing recent results in Chizat et al. (2020). Finally, we illustrate the applicability of our results in obtaining rates of convergence for Wasserstein barycenter between two probability distributions and obtaining asymptotic detection thresholds for some recent optimaltransport based tests of independence.

heading: 1 Introduction
text: Given two random variables X ∼ µ and Y ∼ ν, where µ, ν are probability measures on Rd, d ≥ 1, the problem of finding a “nice" map T0(·) such that T0(X) ∼ ν has numerous applications in machine learning such as domain adaptation and data integration [34, 35, 38, 48, 61, 112], dimension reduction [12, 66, 90], generative models [60, 81, 88, 110], to name a few. Of particular interest is the case when T0(·) is obtained by minimizing a cost function, a line of work initiated by Gaspard Monge [97] in 1781 (see (1.1) below), in which case T0(·) is termed an optimal transport (OT) map and has applications in shape matching/transfer problems [29, 47, 107, 121], Bayesian statistics [46, 75, 80, 108], econometrics [15, 28, 45, 50, 54], nonparametric statistical inference [39– 41, 113, 114]; also see [111, 128, 129] for book-length treatments on the subject. In this paper, we will focus on the OT map obtained using the standard squared Euclidean cost function, i.e.,
T0 := argmin T :T#µ=ν
E‖X − T (X)‖2, (1.1)
where T#µ = ν means T (X) ∼ ν for X ∼ µ. The estimation of T0 has attracted a lot of interest in recent years due to its myriad applications (as stated above) and interesting geometrical properties (see [19, 56, 91] and Definition 1.1 below). In practice, the main hurdle in constructing estimators for T0 is that the explicit forms of the measures µ, ν are unknown; instead only random samples
X1, . . . , Xm ∼ µ and Y1, . . . , Yn ∼ ν
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
are available. A natural strategy in this scenario is to estimate T0 using T̃m,n, where T̃m,n is computed as in (1.1) with µ and ν replaced by µ̃m and ν̃n which are empirical approximations of µ and ν based on X1, . . . , Xm and Y1, . . . , Yn respectively (see Definition 1.2). Such estimators are often called plug-in estimators and have been used extensively; see [7, 30, 67, 93, 94, 102, 116].
The main goal of this paper is to study the rates of convergence of general plug-in estimators of T0 under a unified framework. We show that when µ̃m and ν̃n are chosen as µ̂m and ν̂n respectively, where µ̂m and ν̂n are the standard empirical distributions supported on m and n atoms, i.e.,
µ̂m := 1
m m∑ i=1 δXi and ν̂n := 1 n n∑ j=1 δYj , (1.2)
T̃m,n (appropriately defined using Definition 1.2) converges at a rate of m−2/d + n−2/d for d ≥ 4 in the sense of (1.8). This rate happens to be minimax optimal under minimal smoothness assumptions (see [72, Theorem 6]) but suffers from the curse of dimensionality. We next show that, if µ and ν are known to admit sufficiently smooth densities, it is possible to apply kernel or wavelet based smoothing techniques on µ̂m and ν̂n to obtain plug-in estimators that mitigate the aforementioned curse of dimensionality.
Our next contribution pertains to the estimation of W 22 (µ, ν) (the squared Wasserstein distance), see (1.3) below, a quantity of independent interest in statistics and machine learning with applications in structured prediction [51, 89], image analysis [18, 59], nonparametric testing [16, 106], generative modeling [10, 96], etc. In this paper, we also obtain rates of convergence for plug-in estimators W 22 (µ̃m, ν̃n) of W 2 2 (µ, ν). We show that kernel smoothing µ̂m and ν̂n can be used to obtain plug-in estimators of W 22 (µ, ν) that mitigate the curse of dimensionality as opposed to a direct plug-in approach using µ̂m and ν̂n (as used in [30, Theorem 2]). This provides an answer to the open question of estimating W 22 (µ, ν) when µ, ν admit smooth densities laid out in [30].

heading: 1.1 Background on optimal transport
text: In this section, we present some basic concepts and results associated with the OT problem that will play a crucial role in the sequel. Let Pac(Rd) denote the set of all Lebesgue absolutely continuous probability measures on Rd andP2(Rd) be the set of probability measures with finite second moments. Then the 2-Wasserstein distance (squared) between µ, ν ∈ P2(Rd) is defined as:
W 22 (µ, ν) := min π∈Π(µ,ν)
∫ ‖x− y‖2 dπ(x, y), (1.3)
where Π(µ, ν) is the set of probability measures on Rd×Rd with marginals µ and ν. The optimization problem in (1.3) is often called the Kantorovich relaxation (see [76, 77]) of the optimization problem in (1.1). The existence of a minimizer in (1.3) follows from [129, Theorem 4.1].
Proposition 1.1 (Brenier-McCann polar factorization theorem, see [91, 128]). Suppose µ ∈ Pac(Rd). Then there exists a µ-a.e. (almost everywhere) unique function T0(·) : Rd → Rd, which is the gradient of a real-valued d-variate convex function, say ϕ0(·) : Rd → R, such that T0#µ = ν. Further, the distribution defined as π(A×B) = µ(A ∩ (T0)−1(B)) for all Borel sets A,B ⊆ Rd is the unique minimizer in (1.3) provided µ, ν ∈ P2(Rd).
Definition 1.1 (OT map and potential function). The function T0 : Rd → Rd in Proposition 1.1 which satisfies T0#µ = ν will be called the OT map from µ to ν. A convex function ϕ0(·) in Proposition 1.1 satisfying∇ϕ0 = T0 will be termed an OT potential.
The next and final important ingredient is the alternate dual representation of (1.3) which gives: 1
2 W 22 (µ, ν) =
1
2
∫ ‖x‖2 dµ(x) + 1
2
∫ ‖y‖2 dν(y)−min
f∈F Sµ,ν(f), where (1.4) Sµ,ν(f) = ∫ f dµ+ ∫ f∗ dν. (1.5)
Here F denotes the space of convex functions on Rd which are also elements of L1(µ) and f∗(·) is the standard Legendre-Fenchel dual defined as:
f∗(x) := sup y∈Rd [y>x− f(y)], for x ∈ dom(f). (1.6)

heading: 1.2 Estimating OT map via barycentric projection
text: Recall the setting from the Introduction. Let µ̃m, ν̃n ∈ P2(Rd). Here µ̃m, ν̃n need not be absolutely continuous and can be very general. Intuitively, µ̃m and ν̃n can be viewed as some empirical approximation of µ and ν respectively.
Example 1.2 (Simple choices of µ̃m and ν̃n). Let X1, . . . , Xm i.i.d.∼ µ and Y1, . . . , Yn i.i.d.∼ ν; in which case a natural choice would be to set µ̃m = µ̂m and ν̃n = ν̂n where µ̂m and ν̂n are the empirical distributions of X1, . . . , Xm and Y1, . . . , Yn respectively, as defined in (1.2). This is the standard choice adopted in the discrete-discrete Kantorovich relaxation; see [104, Section 2.3]. Another popular choice is µ̃m = µ̂m, ν̃n = ν or µ̃m = µ, ν̃n = ν̂n. This is the semi-discrete Kantorovich problem and is popular when one of the measures is fully specified; see [26, 55].
A natural way to estimate T0(·), as defined in (1.1), would be to approximate it using the OT map from µ̃m to ν̃n. However as µ̃m and ν̃n may not be elements of Pac(Rd), Proposition 1.1 does not apply and an OT map may not exist from µ̃m to ν̃n. Such is the case in Example 1.2 in the discrete-discrete case when m 6= n. To circumvent this issue, we leverage the notion of barycentric projections (see [3, Definition 5.4.2]) defined below:
Definition 1.2 (Barycentric projection). Define the set
Γ̃min := argmin π∈Π(µ̃m,ν̃n)
∫ ‖x− y‖2 dπ(x, y).
The optimization problem above is the plug-in analog of the optimization problem on the right hand side of (1.3). Given any γ ∈ Γ̃min, define the barycentric projection of γ as the conditional mean of y given x under γ, i.e.,
T̃m,n(x) ≡ T̃ γm,n(x) := ∫ y y dγ(x, y)∫ y dγ(x, y) , for x ∈ supp (µ̃m) . (1.7)
In general, Γ̃min need not be a singleton which is why we index the barycentric projection T̃ γm,n(·) by γ ∈ Γ̃min. Note that T̃ γm,n(·) need not be a transport map; however, if an OT map exists then it must be equal to T̃ γm,n(·) (µ̃m-a.e.). Our goal is to obtain stochastic upper bounds for
sup γ∈Γ̃min ∫ ∥∥T̃ γm,n(x)− T0(x)∥∥2 dµ̃m(x). (1.8) In addition, our proof techniques also yield rates of convergence for∣∣W 22 (µ̃m, ν̃n)−W 22 (µ, ν)∣∣. (1.9) In this paper, we will focus on d ≥ 2. Due to the canonical ordering of R, the case d = 1 can be handled easily using the classical Hungarian embedding theorem [82].

heading: 1.3 Contributions
text: 1. We provide a new and flexible stability estimate Theorem 2.1 which yields a unified approach
to obtaining rates of convergence for general plug-in estimators of the OT map T0(·). Unlike existing stability estimates, Theorem 2.1 holds for the barycentric projection (which is the same as the OT map when it exists) and does not require any smoothness assumptions on µ̃m, ν̃n or T̃ γm,n(·); also see Remark 2.1 for a comparison with the existing literature.
2. In Sections 2.1 and 2.2, we use Theorem 2.1 to bound (1.8) and (1.9):
• In Section 2.1, we show that in both the discrete-discrete and semi-discrete Kantorovich relaxation problems (see Example 1.2), the rate of convergence of (1.8) is m−2/d + n−2/d for d ≥ 4 when T0 is assumed to be Lipschitz (see Theorem 2.2), which is the minimax rate (see [72, Theorem 6]). To the best of our knowledge, rates of convergence for these natural estimators weren’t previously established in the literature.
• In Section 2.2 and Appendix A, we show that the curse of dimensionality in the above rates can be mitigated provided µ and ν admit (uniform) Sobolev smooth densities (see Section 2.2) or Besov smooth densities (see Appendix A). In Section 2.2, our plug-in estimator is obtained by choosing µ̃m (and ν̃n) as the convolution of µ̂m (and ν̂n) and a smooth kernel with an appropriate bandwidth. Under this choice, the rate of convergence in (1.8) is m−( s+2 d ∧ 1 2 ) +
n−( s+2 d ∧ 1 2 ), where s denotes the degree of Sobolev smoothness (see Theorem 2.5). Clearly, if 2(s + 2) ≥ d, the rate of convergence becomes dimension-free and mitigates the curse of dimensionality. We also show the same rates of convergence mentioned above hold for (1.9) (see e.g., Proposition 2.6) which makes a strong case in favor of incorporating smoothness in the construction of plug-in estimators as was conjectured in [30]. In Appendix A, our plug-in estimator is obtained using natural wavelet based density estimators. The rate of convergence in (1.8) turns out to be n− 1+s d+2s where s denotes the degree of Besov smoothness (see Theorem A.1). Note that by choosing s large enough, the exponent in the rate can be made arbitrarily close to 1/2, thereby reducing the curse of dimensionality.
3. In Section 2.3, we use a discretization technique from [131] to construct discrete approximations to the smoothed µ̃m and ν̃n from the previous paragraph that in turn yield computable plug-in estimators for T0 (provided one can sample from µ̃m and ν̃n) that also achieve the same statistical guarantees as the smoothed plug-in estimator from Section 2.2 (see Theorem 2.7). However the number of atoms required in the discretizations and correspondingly the computational complexity increases with the degree of smoothness; this highlights a statistical and computational trade-off.
4. We provide implications of our results in popular applications of OT such as estimating the barycenter of two multivariate probability distributions (see Theorem B.1 in Appendix B.1) and in nonparametric independence testing (see Theorem B.3 in Appendix B.2).

heading: 1.4 Related work
text: Many recent works have focused on obtaining consistent estimators of T0 using the plug-in principle, see [26, 55] (in the semi-discrete problem) and [41, 68, 132] (in the discrete-discrete problem). In [55], the authors studied the rate of convergence of the semi-discrete optimal transport map from ν (absolutely continuous) to µ̂m. This paper complements the aforementioned papers by studying the rates of convergence for general plug-in estimators in a unified fashion. In two other papers [9, Theorem 1.1] and [87, Section 4], the authors use a “Voronoi tessellation" approach to estimate T0, however the rates obtained in this paper, even in the absence of smoothness, are strictly better than those in [9, 87]. Perhaps the most closely related paper to ours would be [67]. In [67], the author uses variational techniques to arrive at stability estimates while we exploit the Lipschitz nature of the OT map (see Definition 1.1). Further the rates in this paper have exponents s+2d ∧ 1 2 which are strictly better than the exponents s+22(s+2)+d obtained in [67, Proposition 1] under the same smoothness assumptions (Sobolev type of order s, see Definition 2.4). In another line of work [72], the authors use theoretical wavelet based estimators (not of the plug-in type) of T0 to obtain nearly minimax optimal rates of convergence. However these estimators, by themselves, are not transport maps between two probability measures, which makes them harder to interpret. In contrast, our focus is on obtaining rates of convergence for plug-in estimators, which are transport maps between natural approximations of µ and ν. Such plug-in type strategies are a lot more popular in computational OT [7, 30, 67, 93, 94, 102, 116].
In terms of obtaining rates of convergence for (1.9), some attempts include [109, 116] where parametric rates are obtained when µ, ν are known to be finitely supported or are both Gaussian. In a related problem, bounds for W 22 (µ̂m, µ) were obtained in [6, 42, 49, 100, 123, 131]. Using these bounds, for m = n, it is easy to get a n−1/d rate of convergence for (1.9). This rate was recently improved to n−2/d in [30] under no smoothness assumptions. Our rates coincide with the n−2/d rate from [30] under no smoothness assumptions. But further, we show in this paper that the curse of dimensionality in the above rate can be mitigated by incorporating smoothness into the plug-in procedure.

heading: 2 Main results
text: Recall ϕ0(·) from Definition 1.1. The following is our main result.
Theorem 2.1 (Stability estimate). Suppose that µ, ν ∈ Pac(Rd) ∩ P2(Rd) and µ̃m, ν̃n ∈ P2(Rd). Assume that T0(·) (as defined in (1.1)) is L-Lipschitz (L > 0). Then,
sup γ∈Γ̃min
∫ ‖T̃ γm,n(x)− T0(x)‖2 dµ̃m(x) ≤ Lmax {∣∣∣∣ ∫ Ψ∗µ̃m,ν̃n d(ν̃n − ν†m)∣∣∣∣, ∣∣∣∣ ∫ Ψ∗µ̃m,ν†m d(ν̃n − ν†m) ∣∣∣∣}
+ 2L ∫ ϕ∗0(y) d(ν̃n − ν†m)(y), (2.1)
where ν†m := T0#µ̃m, ϕ ∗ 0(·) is defined as in (1.6), and with S·,·(·) defined as in (1.5), Ψµ̃m,ν̃n(·) := argminf∈F Sµ̃m,ν̃n(f), Ψµ̃m,ν†m(·) := argminf∈F Sµ̃m,ν†m(f), and D denotes the space of realvalued convex functions on Rd.
The proof of Theorem 2.1 (see Appendix C.1) starts along the same lines as the proof of the curvature estimate in [56, Proposition 3.3]. This is followed by some careful manipulations of W 22 (·, ·) (as in (1.3)) and an application of the conditional version of Jensen’s inequality, see (C.3). The final step of the proof uses the dual representation in (1.4) with techniques similar to some intermediate steps in the proof of [92, Proposition 2] and [30, Lemma 3].
Remark 2.1 (Comparison with other stability estimates). Theorem 2.1 provides some important advantages to existing stability estimates in the literature. One of the earliest results in this direction can be found in [56, Proposition 3.3] but their bound involves a push-forward constraint which makes it hard to use for rate of convergence analysis. A bound similar to Theorem 2.1 is presented in [55, Lemma 5.1] but there the authors assume the existence of an OT map from µ̃m to ν̃n. Therefore, it does not apply to the discrete-discrete problem where µ̃m = µ̂m and ν̃n = ν̂n with m 6= n. Overcoming all these limitations is an important contribution of Theorem 2.1 and allows us to deal with popular plug-in estimators all in one go. The stability estimate in [72, Proposition 10] on the other hand requires µ̃m, ν̃n to be sufficiently smooth and hence it does not hold for discrete-discrete or semi-discrete plug-in estimators (see Example 1.2). Further their result requires all the measures involved to be compactly supported unlike the much milder requirements of Theorem 2.1. However, a shortcoming of Theorem 2.1 is that it is hard to obtain rates faster than n−1/2 using it directly, whereas [72] can obtain rates arbitrarily close to n−1. This is a price we pay for analyzing natural and popular plug-in estimators as opposed to the (more intractable) wavelet based estimators in [72].
Remark 2.2 (How to use Theorem 2.1 to obtain rates of convergence?). Note that the second term on the right hand side of (2.1), under appropriate moment assumptions, is Op(m−1/2 +n−1/2) (free of dimension) by a direct application of Markov’s inequality. We therefore focus on the first term. By (1.5), Ψ∗µ̃m,ν̃n(·), Ψ ∗ µ̃m,ν † m
(·) ∈ F . Further, by Caffarelli’s regularity theory [20–22], depending on the “smoothness" of µ̃m, ν̃n, it can be shown that there exists a further class of functions Fs (see Remarks 2.3 and 2.6) such that Ψ∗µ̃m,ν̃n(·), Ψ ∗ µ̃m,ν † m
(·) ∈ F ∩ Fs. Thus, we can bound the first term on the right hand side of (2.1) as:
max {∣∣∣∣ ∫ Ψ∗µ̃m,ν̃n d(ν̃n − ν†m)∣∣∣∣, ∣∣∣∣ ∫ Ψ∗µ̃m,ν†m d(ν̃n − ν†m) ∣∣∣∣} ≤ sup
f∈F∩Fs ∣∣∣∣ ∫ f d(ν̃n − ν†m)∣∣∣∣. (2.2) The right hand side of (2.2) can now be bounded using the corresponding Dudley’s entropy integral bounds using empirical process techniques; see [126, Lemmas 19.35-19.37].
To conclude, the two main steps in our strategy are identifying the family of functions Fs and computing Dudley’s entropy integral. Further, the more the smoothness of µ̃m, ν̃n, the smaller is the class of functions Fs and smaller the supremum on the right hand side of (2.2). This shows why better rates can be expected under smoothness assumptions.

heading: 2.1 Natural non-smooth plug-in estimators
text: In this case, we discuss the rates of convergence for the discrete-discrete problem and the semi-discrete problem, where no smoothness is available on µ̃m and ν̃n.
Theorem 2.2. Suppose that T0(·) is L-Lipschitz, ν is compactly supported and E exp(t‖X1‖α) <∞ for some t > 0, α > 0.
(Discrete-discrete): Set µ̃m = µ̂m and ν̃n = ν̂n. Then the following holds:
sup γ∈Γ̃min
∫ ‖T̃ γm,n(x)− T0(x)‖2 dµ̃m(x) = Op ( r (m,n) d × (log (1 + max{m,n})) td,α ) , (2.3)
where r(m,n)d :=  m−1/2 + n−1/2 for d = 2, 3, m−1/2 log (1 +m) + n−1/2 log (1 + n) for d = 4, m−2/d + n−2/d for d ≥ 5,
(2.4)
and
td,α :=  (4α)−1(4 + ((2α+ 2dα− d) ∨ 0)) for d < 4, (α−1 ∨ 7/2)− 1 for d = 4, 2(1 + d−1) for d > 4.
The same bound holds for |W 22 (µ̃m, ν̃n)−W 22 (µ, ν)| without assuming T0(·) is Lipschitz. (Semi-discrete): Set µ̃m = µ, ν̃n = ν̂n or µ̃m = µ̂m, ν̃n = ν. Then the left hand side of (2.3) is Op(r (n,n) d × (log (1 + n))td,α) or Op(r (m,m) d × (log (1 +m))td,α) respectively.
A stronger result can be proved if both µ and ν are compactly supported.
Corollary 2.3. Consider the setting from Theorem 2.2 and assume further that µ is compactly supported. Then, with r(m,n)d defined as in (2.4), we have:
E [ sup
γ∈Γ̃min
∫ ‖T̃ γm,n(x)− T0(x)‖2 dµ̃m(x) ] ≤ Cr(m,n)d ,
for some constant C > 0, in both the discrete-discrete and semi-discrete settings from Theorem 2.2.
A brief description of the proof technique of Theorem 2.2 using Theorem 2.1 is provided in Remark 2.3 below, and the actual proof is presented in Appendix C.1.
Remark 2.3 (Proof technique). The proof of Theorem 2.2 proceeds via the strategy outlined in Remark 2.2. We first show that Fs (see Remark 2.2) can be chosen as a certain sub-class of convex functions which are in L2(ν). We then use Dudley’s entropy integral type bounds which in turn requires the bracketing entropy [126, Page 270] of Fs, recently proved in [83, Equation 26]. This strategy is slightly different from that used in the proof of [30, Theorem 2], where the authors assume that µ is compactly supported whereas we only assume the finiteness of E exp(t‖X1‖α) for some t > 0, α > 0. The compactness assumption on µ allows one to further restrict Fs to the class of Lipschitz functions. This additional restriction does not seem to be immediate without the compactness assumption.
As discussed in Section 1.3, the exponents obtained in Theorem 2.2 are minimax optimal, up to multiplicative logarithmic factors, under bare minimal smoothness assumptions (see [72, Theorem 6]). To the best of our knowledge, rates for the discrete-discrete case for m 6= n and those for the semi-discrete case were not known previously in the literature. Our rates are also strictly better than those (for different estimators, based on space tessellations) obtained in [9, 87] and require less stringent assumptions than those in [30]. In the next section, we show how smoothness assumptions can be leveraged to mitigate the curse of dimensionality in Theorem 2.2.

heading: 2.2 Smooth kernel based plug-in estimator: mitigating the curse of dimensionality
text: In this section, we focus on kernel based density estimators for the probability densities associated with µ and ν (see [57, 58, 99, 103, 115]). We will show, using Theorem 2.1, that the corresponding estimators of T0(·) achieve (near) dimension-free rates under sufficient smoothness assumptions. We first introduce the Sobolev class of functions which we will exploit in this subsection to construct estimators that achieve rates of convergence which mitigate the curse of dimensionality under sufficient smoothness.
Definition 2.4 (Uniform Sobolev class of functions). Let Ω ⊆ Rd and f(·) be uniformly continuous on Ω and admits uniformly continuous derivatives up to order s on Ω for some s ∈ N. For any m := (m1, . . . ,md) ∈ Nd, let
∂mf := ∂
∂m1x1 . . .
∂
∂mdxd f, |m| := d∑ i=1 mi.
For any k ≤ s, we further define, ‖f‖Ck(Ω) := ∑ |m|≤k ‖∂mf‖L∞(Ω).
The space Cs(Ω) is defined as the set of functions f(·) for which ‖f‖Ck(Ω) <∞ for all k ≤ s.
For this subsection, assume that µ and ν admit Sobolev smooth densities fµ(·) and fν(·) in the uniform norm (see Definition 2.4 above). Given Ω ⊆ Rd and s ∈ N, let Cs(Ω) denote the set of Sobolev smooth functions on Ω of order s.
Assumption (A1) (Regularity of the densities). Suppose that
1. fµ and fν are supported on compact and convex subsets of Rd, say X and Y respectively. 2. There exists s,M > 0 such that fµ(·) ∈ Cs(X ;M) and fν(·) ∈ Cs(Y;M) where Cs(X ;M) is
the space of real valued functions supported on X such that for all f(·) ∈ Cs(X ;M), we have M−1 ≤ f(x) ≤ M for all x ∈ X and ‖f‖Cs(X ) ≤ M . Here ‖·‖Cs(X ) is the standard uniform Sobolev norm as defined in Definition 2.4. The space Cs(Y;M) is defined analogously.
We now define our estimators for fµ(·) and fν(·) using the standard kernel density estimation technique (see [125, Section 1.2]). Set
f̂µ(x) := 1
mhdm m∑ i=1 Kd ( Xi − x hm ) , (2.5)
for some bandwidth parameter hm > 0 and d-variate kernel Kd(·). We assume that Kd(·) is the dfold product of univariate kernels, i.e., there exists a kernelK(·) such that for u = (u1, . . . , ud) ∈ Rd, Kd(u) = ∏d i=1K(ui). We define f̂ν(·) similarly with the same univariate kernel and bandwidth.
Assumption (A2) (Regularity of the kernel). Assume that K(·) is a symmetric, bounded, s+ 1 times differentiable kernel on Rd with all s+ 1 derivatives bounded and integrable. Further, suppose that K(·) is of order 2s+ 2, i.e.,∫ ujK(u) du = 1(j = 0), for j = {0, 1, 2, . . . , 2s+1}, and ∫ |u|2s+2|K(u)| du <∞.
The above assumptions on K(·) are standard for estimating smooth densities and their derivatives of different orders in the kernel density estimation literature; see e.g. [4, 57, 58, 69, 125]. There are several natural ways to construct kernels satisfying Assumption (A2), see [125, Section 1.2.2]; an example is also provided in Example 2.4 below.
Example 2.4 (Example of a kernel satisfying Assumption (A2)). Let ψm(·) be the m-th Hermite polynomial on R (see [84]). Then the kernel function defined as
K(u) := 2s+2∑ m=0 ψm(0)ψm(u) exp(−u2/2)
satisfies Assumption (A2).
It is evident from Assumption (A2) that K(·) may take some negative values, in which case, f̂µ(·) (respectively f̂ν(·)) may not be a probability density. Consequently the barycentric projection (see Definition 1.2) between f̂µ(·) and f̂ν(·) is not well-defined. We get around this by projecting
f̂µ(·) and f̂ν(·) on an appropriate space of “smooth" probability densities (see (2.6)), via an integral probability metric (see Definition 2.5 below; also see [98, 105, 117] for examples, computational procedures and applications of such metrics).
Definition 2.5 (Integral probability metric). Given a classH of bounded functions on Rd and two probability densities g1(·) and g2(·) on Rd, the integral probability metric/distance between g1(·) and g2(·) with respect toH is defined as
dIP(g1, g2;H) := sup ψ(·)∈H ∣∣∣∣ ∫ ψ(x)(g1(x)− g2(x)) dx∣∣∣∣. Sufficient conditions onH for dIP(·, ·;H) to be a metric on the space of probability measures (not on the space of probability densities as they can be altered on set of Lebesgue measure 0 without altering the underlying probability measures) on Rd have been discussed in [98]. Observe that the measure dIP(g1, g2;H) is well defined even when g1(·) and g2(·) are not probability densities.
In Theorem 2.5 below, we useH = Cs+2(X ,M ′). Note that any function in Cs+2(X ,M ′) can be extended to a function in Cs+2(Rd;M ′) (see [72, Theorem 23] and [124, Theorem 1.105]). The fact that this choice of F results in a metric follows from the argument in [98, Page 8].
We are now in a position to describe the projection estimators for fµ(·) and fν(·), and the rates achieved by the corresponding plug-in estimator.
Theorem 2.5. Assume that T0(·) is L-Lipschitz and fµ, fν are Lebesgue densities satisfying Assumption (A1). Also suppose that K(·) satisfies Assumption (A2). Define hm := m− 1 d+2s logm , hn := n − 1d+2s log n and T := ∫ |Kd(u)| du+ 1. Fix any M ′ > 0. Consider any probability density f̃M ′
µ (·) ∈ Cs(X ;TM) (where M is defined as in Assumption (A1)) which satisfies
dIP
( f̃M ′
µ , f̂µ;C s+2(X ;M ′) ) ≤ inf f(·)∈Cs(X ;TM) f≥0, ∫ f=1 dIP ( f̂µ, f ;C s+2(X ;M ′) ) + r (m,n) d,s (2.6)
where r(m,n)d,s is defined as in (2.7) and dIP(·, ·;Cs+2(X ;M ′)) is the integral probability metric defined in Definition 2.5. We define f̃M ′
ν (·) analogously as in (2.6) with X , f̂µ(·) replaced by Y , f̂ν(·). Then the following conclusions hold.
1. SetM ′ := 8(1+TM). If µ̃m and ν̃n are the probability measures corresponding to the probability densities f̃M ′ µ (·) and f̃M ′ ν (·), then the following holds for some constant C > 0:
E [ sup
γ∈Γ̃min
∫ ‖T̃ γm,n(x)− T0(x)‖2 dµ̃m(x) ] ≤ Cr(m,n)d,s ,
where r(m,n)d,s :=  m−1/2 + n−1/2 for d < 2(s+ 2), m−1/2 (log (1 +m)) d + n−1/2 (log (1 + n))
d for d = 2(s+ 2), m− s+2 d + n− s+2 d for d ≥ 2(s+ 2). (2.7)
The same bound also holds for E|W 22 (µ̃m, ν̃n)−W 22 (µ, ν)|.
2. f̂µ(·) satisfies
lim n→∞ max
{ P ( ‖f̂µ‖Cs(X̃ ) ≥ TM ) ,P ( sup x∈X̃ |f̂µ(x)− fµ(x)| ≥ ε )} = 0 (2.8)
for any ε > 0, where X̃ is any compact subset of X o. The same conclusion holds for f̂ν(·) with X replaced by Y .
In Theorem 2.5, we have shown that the plug-in estimator for T0(·) using f̃M ′ µ (·) and f̃M ′
ν (·) (with M ′ = 8(1 + TM)) achieves rates that mitigate the curse of dimensionality under sufficient smoothness. In fact, f̃M ′
µ (·) can be viewed as an approximate minimizer of dIP(f̂µ, ·;Cs+2(X ,M ′)) over an appropriate class of Sobolev smooth probability densities. This is carried out because f̂µ(·) by itself may not be a probability density.
Further note that µ̃m, ν̃n as specified in Theorem 2.5 are both smooth, and consequently Γ̃min is a singleton and the supremum in Theorem 2.5 can be dropped. A brief description of the proof technique for Theorem 2.5 is presented in Remark 2.6 below and the actual proof is given in Appendix C.1.
Remark 2.6 (Proof technique). The proof of Theorem 2.5 proceeds along the same lines as Remark 2.3. We first show that Fs (see Remark 2.2) can be chosen as a certain subset of Cs+2(Y◦). We then use Dudley’s entropy integral type bounds which in turn requires the bracketing entropy [126, Page 270] of the class of compactly supported Sobolev smooth functions which can be found in [127, Corollary 2.7.2].
We now explain the implications of both the parts of Theorem 2.5 in the following two remarks.
Remark 2.7 (Mitigating the curse of dimensionality). Theorem 2.5 shows that, under enough smoothness, i.e., when 2(s+ 2) > d, both the upper bounds for (1.8) and (1.9) are Op(n−1/2). This shows that, for large dimensions, provided µ and ν admit smooth enough densities, it is possible to construct plug-in estimators that mitigate the curse of dimensionality. Note that a similar estimator was analyzed in [67, Proposition 1] when m = n. However, the rates obtained in Theorem 2.5 are strictly better than those in [67, Proposition 1]. For m = n, when d < 2(s + 2), [67] obtained a rate of n− s+2 2(s+2)+d which is worse than n−1/2 obtained in Theorem 2.5. For the other regimes, [67] obtains rates (up to log factors) of n−1/4 and n− 1
(s+2)(d+2(s+2)) which are both worse than the respective rates of n−1/2 and n− s+2 d in Theorem 2.5.
Remark 2.8 (Computational aspects of Theorem 2.5). Note that f̃M ′ µ (·) (with M ′ = 8(1 + TM)) is hard to compute whereas f̂µ(·) is computable easily in linear time. Note that if f̂µ(·) itself were a probability density in Cs(X ;TM), then we would have f̂µ = f̃M ′
µ . While Theorem 2.5 does not establish that, it does come close in part 2, from which we can easily derive the following:
lim n→∞
P(f̂µ(·) /∈ Cs(X̃ ;TM)) = 0.
The above shows that f̂µ(·) is indeed bounded below by (TM)−1 on X̃ (any compact subset of the interior of X ), and additionally belongs to Cs(X̃ ;TM) with probability converging to 1. This leads us to conjecture that the natural density version of f̂µ(·), i.e.,
max{f̂µ(·), 0}∫ max{f̂µ(x), 0} dx
should serve as a good proxy for f̃M ′
µ (·) and lead to rates of convergence that mitigate the curse of dimensionality. From a computational perspective, the density specified above is easy to simulate from using an accept-reject algorithm without computing the integral in the denominator (see [101, Algorithm 4.3]). However, our current proof technique does not provide rates of convergence for the above density estimator based on f̂µ(·).
Another important implication of Theorem 2.5 is the bound obtained on |W2(µ̃m, ν̃n)−W2(µ, ν)| when µ 6= ν. We first present the result and then describe the implication.
Proposition 2.6. Consider the setting in Theorem 2.5. Then, provided µ 6= ν, the following holds:
|W2(µ̃m, ν̃n)−W2(µ, ν)| = Op(r(m,n)d,s ).
Proposition 2.6 (see Appendix C.1 for a proof) shows an interesting distinction between the µ 6= ν case and the µ = ν case. For µ = ν, the best possible exponent is n− 1+s 2s+d for d ≥ 3 (see [131, Theorem 3] where the result was established under more general Besov smoothness assumptions). On the contrary, when µ 6= ν, Proposition 2.6 establishes a rate of n− s+2d for the Wasserstein distance which is strictly better than the minimax achievable rate mentioned above when µ = ν. This observation complements [30, Corollary 1] where the authors make a similar remark for the special case of s = 0.

heading: 2.3 Discretized plug-in estimator under smoothness assumptions
text: In Section 2.1, we discussed how smoothness can be incorporated into the plug-in procedure to get faster rates of convergence. Such plug-in estimators are popular in the computational OT literature (see [7, 8, 25, 36]). However, even after f̃µ(·) ≡ f̃M ′ µ (·), f̃ν(·) ≡ f̃M ′
ν (·) are calculated, T̃ γm,n as in Theorem 2.5 cannot be computed explicitly from data if f̃µ(·) and f̃ν(·) are continuous densities. This is in contrast to T̃ γm,n from Theorem 2.2 in the discrete-discrete case which is explicitly computable using a standard linear program, but achieves worse rates of convergence. This is not unexpected. Thanks to the no free lunch principle, better statistical accuracy is naturally accompanied by heavier computational challenges. Therefore, our goal here is to construct estimators, under smoothness assumptions as in Section 2.2, which are computable in polynomial time (with complexity increasing with smoothness) provided f̃µ(·) and f̃ν(·) can be sampled from, and also attain rates that mitigate the curse of dimensionality.
Construction: We will illustrate the discretized estimator using the kernel based estimator from Section 2.2. Similar results also hold for the wavelet based estimator from Appendix A. Recall the kernel density estimators f̃µ(·) and f̃ν(·) (see (2.6)). Sample M ≥ 1 random points from both f̃µ(·) and f̃ν(·). Let µ̂m,M and ν̂n,M denote the standard empirical measures on the M points sampled from f̃µ(·) and f̃ν(·) respectively. Finally construct T̃m,n ≡ T̃ γm,n as in Definition 1.2 with µ̃m = µ̂m,M and ν̃n = ν̂n,M . It should be pointed out that a similar construction was also used in [131, Section 6] for estimating probability densities under the Wasserstein loss. Based on this construction, the main result of this section is as follows:
Theorem 2.7. Consider the setting in Theorem 2.5 and the same construction of T̃ γm,n as above. For simplicity, let’s also assume m = n. Accordingly set M = n s+2 2 . Then Γ̃min is a singleton and consequently the following conclusion holds for some constant C > 0:
E [∫ ‖T̃m,n(x)− T0(x)‖2 dµ̃m(x) ] ≤ Cr(n,n)d,s .
The same rates also hold for E|W 22 (µ̃m, ν̃n)−W 22 (µ, ν)|.
The proof of Theorem 2.7 is given in Appendix C.1. Once the empirical measures µ̂m,M and ν̂n,M have been obtained, an explicit computation of T̃m,n as described above requires O(M3) = O(n 3(s+2) 2 ) steps using the Hungarian algorithm, see [73]. This highlights the statistical versus computational trade-off, i.e., in order to mitigate the curse of dimensionality in convergence rates by exploiting smoothness, the computational complexity gets progressively worse by polynomial factors in n. It should be mentioned that (approximate) algorithms faster than the Hungarian algorithm stated above, can be found in [1, 36, 53] to name a few. Due to space constraints, we avoid a detailed discussion on this.
In the above construction, sampling from the smoothed kernel densities f̃µ(·) and f̃ν(·) is crucial. If we would simply draw M bootstrap samples from the empirical distributions µ̂m and ν̂n, the rates of convergence wouldn’t improve from those observed in Theorem 2.2 no matter how large M is.
Acknowledgments and Disclosure of Funding We would like to thank the reviewers for their constructive suggestions that greatly helped improve the quality of the paper. The third author is supported by NSF Grant DMS-2015376.

Review summaries:
This article proposes a study on the rates of estimation of Brenier maps using the plug in "barycentric mapping" estimator. As a consequence of this study authors obtain rates of convergence for $W_2^{2}$ in various context (depending on the regularity of the measures and on the setting: discrete-discrete of semi-discrete). They apply their results for Wasserstein barycenter estimation and Nonparametric independence testing. 

This paper analyze the statistical behavior of the general plug-in estimators of OT defined by barycentric projections. The authors provide a thorough analysis of the rate of convergence for the transport cost and the transport map. They also consider kernel smoothed plug-in estimators and relate its rate of convergence to the smoothness of the densities, which alleviates the curse of dimensionality suffered by the plug-in estimator.

Summary =======  Context: --------  This paper tackles the problem of estimating the rate of convergence of empirical Monge maps (through its Kantorovich barycentric map approximator) to the true Monge map for  absolutely continuous Lebesgue measures.    Contribution ------------  - The authors provide new rates of convergence of the barycentric map using standard  empirical plug-in estimators (O(n^{-2/d})) - If the measures are regular enough (Holder smooth), the dependence on the dimension completely vanishes (which was established for Gaussian distributions https://arxiv.org/pdf/1905.10155.pdf) - The authors propose to take advantage of this smoothness result to sample from smoothened distributions (using a Kernel density estimator) and provide the sufficient conditions to obtain a desired rate. Higher smoothness however requires more samples: a tradeoff between statistical accuracy and computational cost cannot be avoided.  

In standard OT problem, marginal distributions of the data are unknown. The current work considers plug-in estimation using Barycentric projection and derived rate of estimation. The idea of Barycentric projection is to plug the estimated marginal distributions of the data in the OT problem. Conditional mean of given us estimated and rate of estimation is derived. Various plug-in estimators of density of are considered such as empirical CDF and kernel density estimator.

This paper studies the problem of estimating an optimal transport map between two probability distributions on Euclidean space, given IID samples from each of the distributions. The main result (Theorem 2.1) gives a general upper bound, in terms of the dual representation of the 2-Wasserstein distance. This is used to give upper bounds, in probability, on the rate of decay of the estimation risk both for the plugging in empirical distributions in the absence of smoothness assumptions and for kernel estimates under smoothness assumptions. Finally, implications for two applications (Wasserstein barycenter estimation and nonparametric independence testing) are discussed.