Minimax estimation of smooth optimal transport maps

TL;DR

This paper establishes the first minimax estimation rates for smooth optimal transport maps in general dimensions, using wavelet-based estimators and stability analysis, supported by numerical experiments.

Contribution

It introduces a new estimator based on semi-dual optimal transport with wavelet truncation, achieving near minimax optimality under smoothness assumptions.

Findings

Estimator achieves near minimax optimality.

Numerical experiments support theoretical results.

First minimax rates for transport maps in general dimension.

Abstract

Brenier's theorem is a cornerstone of optimal transport that guarantees the existence of an optimal transport map $T$ between two probability distributions $P$ and $Q$ over $\mathbb{R}^d$ under certain regularity conditions. The main goal of this work is to establish the minimax estimation rates for such a transport map from data sampled from $P$ and $Q$ under additional smoothness assumptions on $T$. To achieve this goal, we develop an estimator based on the minimization of an empirical version of the semi-dual optimal transport problem, restricted to truncated wavelet expansions. This estimator is shown to achieve near minimax optimality using new stability arguments for the semi-dual and a complementary minimax lower bound. Furthermore, we provide numerical experiments on synthetic data supporting our theoretical findings and highlighting the practical benefits of smoothness regularization. These are the first minimax estimation rates for transport maps in general dimension.

Equations

$$\underset{T}{\min }{\int }_{\text{mathbbm}{R}^d}\Vert T(x)-x{\Vert }_2^2\mathrm{d}P(x),\text{s.t. }{T}_{\#}P=Q,$$

$$T_{\#}P(A)=P(T^{-1}(A)),\text{for all Borel sets }A.$$

$$\underset{\Gamma }{\min }\int \Vert x-y{\Vert }_2^2\mathrm{d}\Gamma (x,y)\text{s.t. }\Gamma \in \Pi (P,Q),$$

$$\Vert \hat{T}-T_0{\Vert }_{L^2(P)}^2:=\int \Vert \hat{T}(x)-T_0(x){\Vert }_2^2\mathrm{d}P(x).$$

$${\int }_{\text{mathbbm}{R}^d}\Vert \hat{T}(x)-x{\Vert }_2^2\mathrm{d}P(x)-{\int }_{\text{mathbbm}{R}^d}\Vert T_0(x)-x{\Vert }_2^2\mathrm{d}P(x).$$

$$\Vert T-T_0{\Vert }_{L^2(P)}^2\lesssim {\int }_{\text{mathbbm}{R}^d}\Vert T(x)-x{\Vert }_2^2\mathrm{d}P(x)-{\int }_{\text{mathbbm}{R}^d}\Vert T_0(x)-x{\Vert }_2^2\mathrm{d}P(x),$$

$${\partial }^{♭}f=\frac{\partial }{{\partial }_{x_1}^{{♭}_1}}\dots \frac{\partial }{{\partial }_{x_d}^{{♭}_d}}f,\text{ and }|♭|=\sum _{i=1}^d{♭}_i.$$

$$\underset{T}{\min }$$

$${\mathcal{P}}_{\mathrm{M}}(T):=\frac{1}{2}{\int }_{\text{mathbbm}{R}^d}\Vert T(x)-x{\Vert }_2^2\mathrm{d}P(x),$$

$$\underset{\Gamma }{\min }$$

$${\tilde{\mathcal{P}}}_{\mathrm{K}}(\Gamma ):=\frac{1}{2}\int \Vert x-y{\Vert }_2^2\mathrm{d}\Gamma (x,y),$$

$$W_2^2(P,Q):=\underset{\Gamma \in \Pi (P,Q)}{\min }{\tilde{\mathcal{P}}}_{\mathrm{K}}(\Gamma ).$$

$${\tilde{\mathcal{P}}}_{\mathrm{K}}(\Gamma )=$$

$$=$$

$$\underset{\Gamma }{\max }{\mathcal{P}}_{\mathrm{K}}(\Gamma )\text{s.t. }\Gamma \in \Pi (P,Q),\text{where }{\mathcal{P}}_{\mathrm{K}}(\Gamma ):=\int ⟨x,y⟩\mathrm{d}\Gamma (x,y).$$

$$\underset{f,g}{\min }$$

$$\displaystyle\begin{array}[]{ll}&f(x)+g(y)\geq\langle x,y\rangle\,,\ P\otimes Q\text{-a.e},\\ &f\in L^{1}(P),g\in L^{1}(Q).\end{array}$$

$$g_f(y)=\underset{x\in \text{mathbbm}{R}^d}{\sup }⟨x,y⟩-f(x)=f^{\ast }(y),$$

$$\underset{f}{\min }$$

$${\mathcal{T}}_{\alpha }(M,R)=\{T\in \mathcal{T}(M):T\text{ is }\lfloor \alpha \rfloor \text{-times differentiable and }\Vert T{\Vert }_{C^{\alpha }(\tilde{\Omega })}\le R\}.$$

$$\underset{\hat{T}}{\inf }\underset{P\in \mathcal{M},T_0\in {\mathcal{T}}_{\alpha }}{\sup }\text{mathbbm}E\left[\int \Vert \hat{T}(x)-T_0(x){\Vert }_2^2\mathrm{d}P(x)\right]\gtrsim n^{-\frac{2\alpha }{2\alpha -2+d}}\vee \frac{1}{n}.$$

$$\underset{P\in \mathcal{M},T_0\in {\mathcal{T}}_{\alpha }}{\sup }\text{mathbbm}E\left[\int \Vert \hat{T}(x)-T_0(x){\Vert }_2^2\mathrm{d}P(x)\right]\lesssim n^{-\frac{2\alpha }{2\alpha -2+d}}(\log (n){)}^2\vee \frac{1}{n}.$$

$$\underset{\hat{T}}{\inf }\underset{P\in \mathcal{M},T_0\in {\mathcal{T}}_{\alpha }}{\sup }\text{mathbbm}E\left[\int \Vert \hat{T}(x)-T_0(x){\Vert }_2^2\mathrm{d}P(x)\right]\gtrsim n^{-\frac{2\alpha }{2\alpha -2+d}}\vee \frac{1}{n}.$$

$$g(x)=\prod _{i\in [d]}\xi (x_i),x=(x_1,\dots ,x_d),$$

$$g_j(x)=\frac{\kappa }{m^{\alpha +1}}g(m(x-x^{(j)})),$$

$${\varphi }_k(x)=\frac{1}{2}\Vert x{\Vert }^2+\sum _{j\in [m{]}^d}{\tau }_j^{(k)}{g}_j(x),k=0,\dots ,K.$$

$${\lambda }_{\min }(D^2{\varphi }_k(x))\ge 1-\sum _{j\in [m{]}^d}{\lambda }_{\max }(D^2{g}_j(x))\ge \frac{1}{2}.$$

$$d(T_k,T_{k^{\prime }}{)}^2={\int }_{[0,1{]}^d}\Vert T_k-T_{k^{\prime }}{\Vert }^2\mathrm{d}x,,0\le k,k^{\prime }\le K.$$

$${\int }_{[0,1{]}^d}\Vert \nabla {\varphi }_k(x)-\nabla {\varphi }_{k^{\prime }}(x){\Vert }^2\mathrm{d}x$$

$$=$$

$${\int }_{[0,1{]}^d}\Vert \nabla {\varphi }_k(x)-\nabla {\varphi }_{k^{\prime }}(x){\Vert }^2\mathrm{d}x\gtrsim n^{-\frac{2\alpha }{2\alpha -2+d}},$$

$$D(Q\|P)=\mathbbm{E}\log\big{(}\frac{\mathrm{d}Q}{\mathrm{d}P}(W)\big{)}\,,\quad W\sim Q.$$

$$\frac{\mathrm{d}{Q}_k}{\mathrm{d}P}(y):=\frac{1}{\det D^2{\varphi }_k((\nabla {\varphi }_k{)}^{-1}(y))}\text{mathbbm}1((\nabla {\varphi }_k{)}^{-1}(y)\in [0,1{]}^d).$$

Full text

62G05,

\startlocaldefs\endlocaldefs\arxiv

1905.05828

1 Introduction

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 $\mathbbm{R}^{d}$, how can we transport $P$ to $Q$ while minimizing the total distance traveled by this transport. A classical instantiation of this problem over $\mathbbm{R}^{d}$ is to find a map $T_{0}\colon\mathbbm{R}^{d}\to\mathbbm{R}^{d}$ that minimizes the objective

[TABLE]

which is known as the Monge problem, where $T_{\#}P$ denotes the push-forward of $P$ under $T$, that is,

[TABLE]

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 $\Gamma_{0}$ in the set of all possible probablity measures on $\mathbbm{R}^{d}\times\mathbbm{R}^{d}$ whose marginals coincide with $P$ and $Q$, which we denote by $\Pi(P,Q)$. This leads to the optimization problem

[TABLE]

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 $\Gamma_{0}$ to (1.3) is concentrated on the graph of a map $T_{0}$. That is, using a suggestive informal notation, $\Gamma_{0}(x,y)=P(x)\,\delta_{y=T_{0}(x)}$, where $\delta$ denotes a point mass. Moreover, $T_{0}$ is the gradient of a convex function $f_{0}$. While cost functions other than $\|x-y\|_{2}^{2}$ could be of interest, such as $\|x-y\|_{2}^{p}$ for $p\geq 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},\dots,X_{n}$ from $P$ and $Y_{1},\dots,Y_{n}$ from $Q$, based on which we would like to find an estimator $\hat{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 efficiently due to recent algorithmic advances [Cut13, AWR17, PC19, DGK18]. However, even if the resulting optimizer $\hat{\Gamma}$ 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\notin\{X_{1},\ldots,X_{n}\}$. In contrast, we would like to obtain an estimator $\hat{T}$ with guarantees in $L^{2}(P)$, that is, with convergence of

[TABLE]

Note that such an estimator of $T_{0}$ could be obtained by post-processing the above optimizer $\hat{\Gamma}$, 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}(\hat{P},\hat{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 suffers from the so-called curse of dimensionality. For example, if both $P$ and $Q$ are absolutely continuous with respect to the Lebesgue measure in $\mathbbm{R}^{d},d\geq 4$, it is known that up to logarithmic factors, $\mathbbm{E}[|W_{2}(\hat{P},\hat{Q})-W_{2}(P,Q)|]\asymp 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 $\|\hat{T}-T_{0}\|_{L^{2}(P)}$ to the sub-optimality gap

[TABLE]

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

[TABLE]

provided $T_{0}=\nabla f_{0}$ is Lipschitz continuous on $\mathbbm{R}^{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.14) 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], $\ell^{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},\ldots,Y_{n}$ is subject to measurement noise. There, the main statistical difficulty 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 difficult 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 differences: 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\alpha/(2\alpha-2+d)}$ (up to log factors) under assumptions quantitatively comparable to Gunsilius’s, while his results imply a sub-optimal rate of $n^{-(2\alpha-2)/(2\alpha-2+d)}$, see Section E in the appendix.

We note that both in the application of Caffarelli’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 effect, 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).

Notation. For any positive integer $m$, define $[m]:=\{1,\ldots,m\}$. We write $|A|$ for the cardinality of a set $A$. The relation $a\lesssim b$ is used to indicate that two quantities are the same up to a constant $C$, $a\leq Cb$. The relation $\gtrsim$ is defined analogously, and we write $a\asymp b$ if $a\lesssim b$ and $a\gtrsim b$. We denote by $c$ and $C$ constants that might change from line to line and that may depend on all parameters of the statistical problem except $n$. We abbreviate with $a\vee b$, $a\wedge b$ the maximum and minimum of $a\in\mathbbm{R}$ and $b\in\mathbbm{R}$, respectively. For $a\in\mathbbm{R}$, the floor and ceiling functions are denoted by $\lfloor a\rfloor$ and $\lceil a\rceil$, indicating rounding $a$ to the next smaller and larger integer, respectively. We use $\operatorname{supp}f$ to denote the support of a function or measure $f$, and $\operatorname{diam}\Omega$ for the diameter of a set $\Omega\subseteq\mathbbm{R}^{d}$. We denote by $B_{1}$ the unit-ball with respect to the Euclidean distance in $\mathbbm{R}^{d}$, where $d$ should be clear from the context. For a real symmetric matrix $A$ and $\lambda\in\mathbbm{R}$, we write $\lambda\preceq A$ if all eigenvalues of $A$ are bounded below $\lambda$, and similarly for $A\preceq\lambda$. Moreover, we denote the smallest and largest eigenvalues of $A$ by $\lambda_{\min}(A)$, $\lambda_{\max}(A)$, respectively.

For $p\in[1,\infty]$, we denote by $\ell^{p}$ either the space $\mathbbm{R}^{d}$ endowed with the usual $\ell^{p}$ norms $\|\,.\,\|_{p}$, or, by abuse of notation, the spaces of multi-dimensional sequences $\gamma:\mathbb{Z}^{m}\to\mathbbm{R}$ with $\|\gamma\|_{\ell^{p}}^{p}=\sum_{k\in\mathbb{Z}^{m}}|\gamma_{k}|^{p}$ for $m\in\mathbb{N}$. Further, for $p\in[1,\infty]$, we denote by $L^{p}$ the Lebesgue spaces of equivalence classes of functions on $\mathbbm{R}^{d}$ or subsets $\Omega\subseteq\mathbbm{R}^{d}$ with respect to the Lebesgue measure $\lambda$ on $\mathbbm{R}^{d}$, whose norms we denote by $\|\,.\,\|_{L^{p}(\mathbbm{R}^{d})}$ and $\|\,.\,\|_{L^{p}(\Omega)}$, respectively. By abuse of notation, for a different measure $P$, we denote the associated Lebesgue norms by $\|\,.\,\|_{L^{p}(P)}$. We abbreviate with “a.e.” any statement that holds “almost everywhere” with respect to the Lebesgue measure.

For a differentiable one-dimensional function $f\colon\mathbbm{R}\supseteq\Omega\to\mathbbm{R}$, we denote its derivative by $\frac{d}{dx}f$. For a function $f\colon\mathbbm{R}^{d}\supseteq\Omega\to\mathbbm{R}$, we denote by $\partial_{i}=\partial/(\partial x_{i})$ its weak derivative in the sense of distributions in direction $x_{i}$, which coincides with the usual (point-wise) derivative if $f$ is differentiable in $\Omega$. For a multi-index $\flat\in\mathbb{N}^{d}$, we set

[TABLE]

The symbol $\partial f$ is also used to denote the sub-differential of a convex function $f$, while we use the symbols $\nabla f$ for the gradient of a function $f$ and $Dg$ for the derivative of a vector-valued function $g\colon\mathbbm{R}^{d_{1}}\supseteq\Omega\to\mathbbm{R}^{d_{2}}$, $\nabla f=(\partial_{1}f,\dots,\partial_{d}f)^{\top}$ and $Dg=(\nabla g_{1},\dots,\nabla g_{d_{2}})^{\top}$ respectively, and $D^{2}f=D\nabla f$ denotes the Hessian of $f$.

If $\Omega\subseteq\mathbbm{R}^{d}$ is a closed set with non-empty interior and $\alpha>0$, the Hölder spaces on $\Omega$ as defined in Appendix B are denoted by $C^{\alpha}(\Omega)$ and their associated norms by $\|\,.\,\|_{C^{\alpha}(\Omega)}$. Similarly, the $p$-Sobolev spaces of order $\alpha$ for $p\in[1,\infty]$ are denoted by $W^{\alpha,p}(\Omega)$ with norms $\|\,.\,\|_{W^{\alpha,p}(\Omega)}$, as defined in Appendix B.

We say that $\Omega\subseteq\mathbbm{R}^{d}$ is a Lipschitz domain if its boundary can be locally expressed as the sublevel set of Lipschitz functions [Tri06, Definition 1.103].

2 Brenier’s theorem and the semi-dual problem

We begin by recalling the Monge and Kantorovich problems given in Section 1. Let $P,\,Q$ be two Borel probability measures on $\mathbbm{R}^{d}$ with finite second moments.

The Monge (primal) problem is defined as

[TABLE]

and the push-forward $T_{\#}P$ is defined as $T_{\#}P(A)=P(T^{-1}(A))$ for all Borel sets $A$.

Its relaxation, the Kantorovich (primal) problem, is given by

[TABLE]

and $\Pi(P,Q)$ denotes the set of couplings between $P$ and $Q$, that is, the set of probability measures $\Gamma$ on $\mathbbm{R}^{d}\times\mathbbm{R}^{d}$ such that $\Gamma(A\times\mathbbm{R}^{d})=P(A)$ and $\Gamma(\mathbbm{R}^{d}\times A)=Q(A)$ for all Borel set $A\subset\mathbbm{R}^{d}$.

The value of problem (2.3) is the square of the 2-Wasserstein distance, denoted by

[TABLE]

Note that we can expand the objective in (2.3) as

[TABLE]

Since the first two terms above do not depend on $\Gamma$, we obtain the equivalent optimization problem

[TABLE]

We focus on this equivalent formulation for the rest of the paper because it is more convenient to work with.

Problem (2.8) is a linear optimization problem, albeit an infinite-dimensional one. Hence, it is natural to consider its dual problem:

[TABLE]

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.9) can be further simplified: Assume we are given a candidate function $f$ in (2.9) 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$:

[TABLE]

Plugging solution (2.13) back into the optimization problem leads to the so-called semi-dual problem,

[TABLE]

where the supremum in (2.13) is interpreted as an essential supremum with respect to $P$. By transitioning to the semi-dual, we effectively solved for all constraints in (2.9), 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.8) exists and is of the form $\Gamma_{0}=(\mathrm{id},T_{0})_{\#}P$, where $T_{0}=\nabla f_{0}$ is the gradient of a convex function $f_{0}\colon\mathbbm{R}^{d}\to\mathbbm{R}$. In fact, $f_{0}$ can be chosen to be a minimizer of the semi-dual objective in (2.14).

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.

3 Main results

Let $X_{1:n}=(X_{1},\ldots,X_{n})$ and $Y_{1:n}=(Y_{1},\ldots,Y_{n})$ be $n$ independent copies of $X\sim P$ and $Y\sim 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\geq 2$.

A1****Source distribution.

Let $\mathcal{M}=\mathcal{M}(M)$ be the set of all probability measures $P$ with support $\Omega=[0,1]^{d}$ that admit a density $\rho_{P}$ with respect to the Lebesgue measure such that $\displaystyle M^{-1}\leq\rho_{P}(x)\leq M$ for almost all $x\in\Omega$. Assume that the source distribution $P$ is in $\mathcal{M}$.

A2****Transport map.

Let $\tilde{\Omega}=[-1,2]^{d}$ denote the enlargement of $\Omega$ by $1$ in every direction. Let $\mathcal{T}=\mathcal{T}(M)$ be the set of all differentiable functions $T\colon\tilde{\Omega}\to\mathbbm{R}^{d}$ such that $T=\nabla f$ for some differentiable convex function $f\colon\tilde{\Omega}\to\mathbbm{R}^{d}$ and

1. (i)

$\displaystyle|T(x)|\leq M$* for all $x\in\tilde{\Omega}$,* 2. (ii)

$\displaystyle M^{-1}\preceq DT(x)\preceq M$* for all $x\in\tilde{\Omega}$,* 3. (iii)

$\displaystyle\operatorname{supp}T_{\#}P=\Omega=[0,1]^{d}$* .*

For $R>1$ and $\alpha>1$, assume that $T_{0}\in\mathcal{T}_{\alpha}=\mathcal{T}_{\alpha}(M,R)$, where

[TABLE]

Our main result is the following theorem. It characterizes, up to logarithmic factors, the minimax rate of estimation of an $\alpha$-smooth transport map $T_{0}\in\mathcal{T}_{\alpha}$ in the setup described above.

Theorem 2**.**

Fix $\alpha\geq 1$, then

[TABLE]

where the infimum is taken over all measurable functions $\hat{T}$ of the data $X_{1:n}=(X_{1},\dots,X_{n})$, $Y_{1:n}=(Y_{1},\dots,Y_{n})$. Moreover, if $P\in\mathcal{M}$ and $T_{0}\in\mathcal{T}_{\alpha}$, there exists an estimator $\hat{T}$, given in Section 5, that is near minimax optimal. More specifically, there exists an integer $n_{0}=n_{0}(d,\alpha,M,R)$ such that for any $n\geq n_{0}$, it holds,

[TABLE]

Remark 3**.**

Assumption A1 can be significantly relaxed with respect to the geometry of $\Omega$ 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 $\operatorname{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 $\tilde{\Omega}$ 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 $\operatorname{supp}(P),\operatorname{supp}(Q)$ and $\tilde{\Omega}$—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}\preceq DT(x)$, in particular on the support of $P$, $x\in\Omega$. This yields convergence estimates for the optimal transport map as shown in [Gig11], see (1.6), 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-off. All remaining assumptions in A1 and A2, including the existence of extensions of $T_{0}$ to a superset $\tilde{\Omega}$, 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^{-\frac{2\alpha}{2\alpha-2+d}}$. Note first that it appears from this rate that estimation of transport maps, like the estimation of smooth functions suffers from the curse of dimensionality. However, as $\alpha\to\infty$, 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 $\alpha\to 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 Caffarelli’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 $\tilde{Q}=\hat{T}_{\#}P$, estimation rates for $\hat{T}$ immediately translate into rates for estimating $Q$ in the $2$-Wasserstein distance [WB17, SP18, WB19]. In fact, $W_{2}^{2}(\tilde{Q},Q)$ can be bounded by $\mathbbm{E}_{P}[\|\hat{T}-T_{0}\|_{2}^{2}]$, since $(\hat{T},T_{0})_{\#}P$ is a candidate transport plan between $Q$ and $\tilde{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.

4 Lower bound

In this section we begin by proving the lower bound (3.1) as it sheds light on the source of the non-standard exponent $\frac{2\alpha}{2\alpha-2+d}$ in the minimax rate. We prove the following theorem.

Theorem 6**.**

Fix $\alpha\geq 1$. It holds that

[TABLE]

where the infimum is taken over all measurable functions of $X_{1},\dots,X_{n}$, $Y_{1},\dots,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=\mathsf{Unif}([0,1]^{d})\in\mathcal{M}$, the uniform distribution on the hypercube. For $\alpha>1$, let $\xi\colon\mathbbm{R}\to\mathbbm{R}$ be a non-zero function in $C^{\infty}(\mathbbm{R})$ with support contained in $[0,1]$ such that there exists $x_{0}\in[0,1]$ with $\xi(x_{0})\neq 0$, $\frac{d}{dx}\xi(x_{0})\neq 0$, for example a bump-function. Define $g\colon\mathbbm{R}^{d}\to\mathbbm{R}$ by

[TABLE]

and note that $\nabla g(x_{0},\dots,x_{0})\neq 0$ and $\operatorname{supp}(g)=[0,1]^{d}$ by the above assumptions on $\xi$.

Let $m=\lceil\theta n^{\frac{1}{2\alpha-2+d}}\rceil$ be a positive integer where $\theta$ 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^{(j)}\,:\,j\in[m]^{d}\}\subset[0,1]^{d}$ to have coordinates $x^{(j)}_{i}=(j_{i}-1)/m$, $i=1,\ldots,d$ and let

[TABLE]

for a constant $\kappa>0$ to be chosen later. Note that $\operatorname{supp}(g_{j})\subseteq x^{(j)}+[0,1/m]^{d}$ and hence that the supports of the functions $\{g_{j}\}_{j\in[m]^{d}}$ are pairwise disjoint.

Next, let $\flat\in\mathbb{N}^{d}$ be a multi-index and observe that the differential operator $\partial^{\flat}$ applied to $g_{j}$ yields $\partial^{\flat}g_{j}(\cdot)=m^{|\flat|-\alpha-1}\partial^{\flat}g(m(\cdot-x_{j}))$. Since $\xi\in C^{\alpha+1}$, if $\alpha>1$, a second-order Taylor expansion yields that $g_{j}$ has uniformly vanishing Hessian: $\|D^{2}g_{j}\|_{L^{\infty}(\mathbbm{R}^{d})}\to 0$ as $m\to\infty$. In particular, in that case, there exists $m_{0}$ such that $\|D^{2}g_{j}(x)\|_{\mathrm{op}}\leq 1/2$ for all $x\in\mathbbm{R}^{d}$, $m\geq m_{0}$, $j\in[m]^{d}$. If $\alpha=1$, the same can be obtained by choosing $\kappa$ small enough. By the same reasoning, we can also guarantee $\|\nabla g_{j}\|_{L^{\infty}(\mathbbm{R}^{d})}\leq 1$.

For $m^{d}\geq 8$, the Varshamov-Gilbert lemma, Lemma 31, guarantees the existence of binary vectors $\tau^{(0)},\tau^{(1)},\ldots,\tau^{(K)}\in\{0,1\}^{[m]^{d}}$, $\tau^{(0)}=(0,\dots,0)$, $K\geq 2^{m^{d}/8}$ such that $\|\tau^{(k)}-\tau^{(k^{\prime})}\|_{2}^{2}\geq m^{d}/8$ for $0\leq k\neq k^{\prime}\leq K$. With this, we define the following collection of Kantorovich potentials:

[TABLE]

It is easy to see (Lemma 38) that for any $k=0,\ldots,K$ and $m\geq m_{0}$, $\nabla\phi_{k}$ is a bijection from $[0,1]^{d}$ to $[0,1]^{d}$. Moreover, by Weyl’s inequality and the above bound $\|D^{2}g_{j}(x)\|_{\mathrm{op}}\leq 1/2$, for all $k$,

[TABLE]

Similarly, we obtain $\|\nabla\phi_{k}\|_{L^{\infty}([0,1]^{d})}\leq 2$ and $\lambda_{\mathrm{max}}(D^{2}\phi_{k}(x))\leq 2$ for all $x\in\mathbbm{R}^{d}$. Hence, $T_{k}:=\nabla\phi_{k}\in\mathcal{T}_{\alpha}(M,R)$ for $M>2$ and $\kappa$ small enough. We now check the conditions of Theorem 30, where we consider the distance measure

[TABLE]

First, observe that for $0\leq k\neq k^{\prime}\leq K$, it holds that

[TABLE]

This yields

[TABLE]

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

[TABLE]

In view of Lemma 38, for any $k=0,\ldots,K$, the measure $Q_{k}=(\nabla\phi_{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

[TABLE]

Moreover, let $X\sim P$ and $Y\sim Q_{k}$ be two random variables. It holds

[TABLE]

Recall that $D^{2}\phi_{k}=I_{d}+\sum_{j\in[m]^{d}}\tau^{(k)}_{j}D^{2}g_{j}$ where $I_{d}$ denotes the identity matrix in $\mathbbm{R}^{d}$. Therefore, since the functions $g_{j}$ have disjoint support, we have for all $x\in[0,1]^{d}$ that

[TABLE]

where $\lambda_{l}(A)$ denotes the $l$th eigenvalue of a matrix $A$. Since $\log(1+z)\geq z-z^{2}/2$ for all $z>0$,

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm. Thus,

[TABLE]

On the one hand, by the divergence theorem and the fact that $g_{j}$ has bounded support,

[TABLE]

where $\partial[0,1]^{d}$ denotes the boundary of the unit hypercube and $v(x)$ its outward-pointing unit normal vector. On the other hand

[TABLE]

The above three displays yield

[TABLE]

for $\theta$ 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 $\phi_{0}(x)=\|x\|_{2}^{2}/2$ and $\phi_{1}(x)=\phi_{0}(x)+(\tilde{\theta}/\sqrt{n})g(x)$ for $\tilde{\theta}$ chosen to ensure $\phi_{0},\phi_{1}\in\mathcal{T}_{\alpha}$, applying Theorem 32 in Appendix D. The separation condition is given by

[TABLE]

and the KL divergence between the associated probability distributions can be estimated by

[TABLE]

for $\tilde{\theta}$ large enough. ∎

Looking back at this proof, we get a better understanding of the exponent in the minimax rate $n^{-\frac{2\alpha}{2\alpha-2+d}}$. Given that $n^{-\frac{2(\beta-k)}{2\beta+d}}$ is the minimax rate of estimation of the $k$th derivative of a $\beta$-smooth density in $L^{2}$ [MG79], the rate that we obtain is formally that of an “antiderivative” ($k=-1)$ of a $\beta=\alpha-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 $\alpha+1$-smooth Kantorovich potential, see (4.8), which is $\alpha-1$-smooth. This follows directly from the Monge-Ampère equation (4.7). On the other hand, we measure the performance of the estimator in terms of the $L^{2}(P)$ distance between $\hat{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\colon\mathbbm{R}^{d}\to\mathbbm{R}^{d}$, the existence of an antiderivative needs to be guaranteed a priori, as in the case of transport maps by assuming $T=\nabla f$ for $f\colon\mathbbm{R}^{d}\to\mathbbm{R}$.

Similar rates arise in the estimation of the invariant measure of a diffusion 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 different 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].

5 Upper bounds

In this section, we give an estimator $\hat{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 tradeoff.

5.1 Assumptions

Throughout, we fix two constants $M\geq 2,\beta>1$.

B1****Source distribution.

Let $\mathcal{M}=\mathcal{M}(M)$ be the set of all probability measures $P$ whose support $\Omega_{P}\subseteq MB_{1}$ is a bounded and connected Lipschitz domain, and that admit a density $\rho_{P}$ with respect to the Lebesgue measure such that $\displaystyle M^{-1}\leq\rho_{P}(x)\leq M$ for almost all $x\in\Omega_{P}$. Assume that the measure $P\in\mathcal{M}$.

B2****Transport map.

For any $P\in\mathcal{M}$ with support $\Omega_{P}$, let $\tilde{\Omega}_{P}$ denote a convex set with Lipschitz boundary such that $\tilde{\Omega}_{P}\subseteq MB_{1}$, and $\Omega_{P}+M^{-1}B_{1}\subseteq\tilde{\Omega}_{P}$. Let $\tilde{\mathcal{T}}=\tilde{\mathcal{T}}(M)$ be the set of all differentiable functions $T\colon\tilde{\Omega}_{P}\to\mathbbm{R}^{d}$ such that $T=\nabla f$ for some differentiable convex function $f\colon\tilde{\Omega}_{P}\to\mathbbm{R}^{d}$ and

1. (i)

$\displaystyle|T(x)|\leq M$* for all $x\in\tilde{\Omega}_{P}$,* 2. (ii)

$\displaystyle M^{-1}\preceq DT(x)\preceq M$* for all $x\in\tilde{\Omega}_{P}$.*

For $R>1$ and $\alpha>1$, assume that

[TABLE]

These new conditions have two implications. First, they imply regularity of the Kantorovich potential $f_{0}$, where $T_{0}=\nabla f_{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 $\mathcal{X}=\mathcal{X}(M)$ be the set of all twice continuously differentiable functions $f\colon\tilde{\Omega}_{P}\to\mathbbm{R}$ such that

1. (i)

$\displaystyle|f(x)|\leq 2M^{2}$* and $\displaystyle|\nabla f(x)|\leq M$ for all $x\in\tilde{\Omega}_{P}$,* 2. (ii)

$\displaystyle M^{-1}\preceq D^{2}f(x)\preceq M$* for all $x\in\tilde{\Omega}_{P}$.*

Then there exists a Kantorovich potential $f_{0}\in\mathcal{X}(M)$ such that $T_{0}=\nabla f_{0}$,

[TABLE]

Further, $Q=\nabla(f_{0})_{\#}P$ has a connected and bounded Lipschitz support $\Omega_{Q}\subseteq MB_{1}$ and admits a density $\rho_{Q}$ with respect to the Lebesgue measure that satisfies $M^{-(d+1)}\leq\rho_{Q}(y)\leq M^{d+1}$ for all $y\in\Omega_{Q}$.

Note that the simplified Assumptions A1 and A2 from Section 3 follow from B1 and B2 in the case $\Omega_{P}=\Omega=[0,1]^{d}$ and $\tilde{\Omega}_{P}=\tilde{\Omega}=[-1,2]^{d}$. Additionally, the simplified assumptions restrict the class of transport maps to those such that $\Omega_{Q}=[0,1]^{d}$ and for which $\beta=\alpha$. Indeed, noting that $\|T_{0}\|_{W^{\alpha,2}(\tilde{\Omega})}\lesssim\|T_{0}\|_{C^{\alpha}(\tilde{\Omega})}$, we can fold the two smoothness conditions into one.

5.2 Estimator

To construct an estimator for $T_{0}$, we observe that if we had access to a Kantorovich potential $f_{0}$, then $T_{0}=\nabla f_{0}$ by Brenier’s Theorem, Theorem 1. In turn, $f_{0}$ is the minimum of the semi-dual objective (2.14). 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

[TABLE]

where, as in Section 3, $X_{1:n}=(X_{1},\dots,X_{n})$ and $Y_{1:n}=(Y_{1},\dots,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 $\mathcal{S}(f)=Pf+Qf^{\ast}$ for $f\in L^{1}(P)$, where $f^{\ast}$ denotes the convex conjugate of $f$. Replacing both $P$ and $Q$ by their empirical counterparts, we obtain the empirical semi-dual,

[TABLE]

In order to incorporate smoothness regularization into the minimization of (5.5), 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}(\mathbbm{R}^{d})$, leading to nested subspaces

[TABLE]

that roughly correspond to increasing frequency ranges of the continuous Fourier transform of a function $f\in L^{2}(\mathbbm{R}^{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}(\tilde{\Omega}_{P})$ of those basis functions of the wavelet expansion whose support has non-trivial intersection with $\tilde{\Omega}_{P}$. This is a finite-dimensional vector space as long as the elements of the wavelet basis have compact support. The cut-off 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 sufficient 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 $\tilde{\Omega}_{P}$ by demanding $f\in\mathcal{X}(2M)$.

Combined, both wavelet regularization and strong convexity lead to the set

[TABLE]

of candidate potentials, based on which we define the estimators

[TABLE]

for the Kantorovich potential and transport map, respectively.

Note that since we consider candidate potentials only on the compact set $\tilde{\Omega}_{P}$, $f^{\ast}$ above is defined as

[TABLE]

where $\iota_{\tilde{\Omega}_{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}\in\mathbb{N}$ and $J$ such that for $n\geq n_{0}$,

[TABLE]

where $n_{0},\,C$, and $J$ may depend on $d,M,R,\Omega_{P},\Omega_{Q},\tilde{\Omega}_{P}$, $n_{0}$ may additionally depend on $\beta$, $C$ may additionally depend on $\alpha$, and $J$ depends on $n$.

The cutoff $J$ depends on $\alpha$ if $d\geq 3$, but in the cases $d=1$ and $d=2$, $J$ can be chosen independently from $\alpha$.

Remark 9**.**

A few remarks are in order.

1. (i)

Similar upper bounds hold with high probability and can be inferred from the proof. 2. (ii)

As written, the estimator $\hat{f}_{J}$ is not directly implementable since the calculation of $f^{\ast}$ involves computing a maximum over a continuous subset of $\mathbbm{R}^{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. 3. (iii)

The numerical experiments in Section 6 suggest that restricting the optimization in (5.8) to $\mathcal{X}(2M)$, while necessary for our proofs, might not be necessary in practice. 4. (iv)

The estimator employed in Theorem 8 can be made adaptive to the unknown smoothness parameter $\alpha$ 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 **[H1̈9]**. 5. (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}\geq 1$, i.e., there is an exponential dependence on the dimension $d$. Further, the dependence on $R$ is captured in (5.47) below and amounts to $R^{\frac{2(d-2)}{2\alpha-2+d}}\log(R)$ in the case $d\geq 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 tradeoff that achieves the minimax lower bound of Theorem 6 up to log factors.

5.3 Stability of optimal transport maps

In this section, we leverage the assumed regularity of the optimal transport map to relate the suboptimality gap in the semi-dual objective function $\mathcal{S}$ and the $L^{2}$-distance of interest.

Proposition 10**.**

Under assumptions B1 – B2, for all $f\in\mathcal{X}(2M)$ as defined in Proposition-Definition 7, we have

[TABLE]

and

[TABLE]

Proof.

It follows from Proposition-Definition 7(ii) and a second-order Taylor expansion that $f$ is of quadratic type [Kol11] around every $x\in\Omega_{P}$:

[TABLE]

for all $L\geq 2M$. It turns out that these conditions are sufficient to obtain the desired result.

The upper bound in (5.13) is of the form

[TABLE]

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

[TABLE]

The squared distance term vanishes for $L=4M$: by the triangle inequality $\|\nabla f_{0}(x)-\nabla f(x)\|_{2}\leq 4M$ and since $x\in\Omega_{P}$, it holds that

[TABLE]

Together with the fact that $Q=(\nabla f_{0})_{\#}P$, this yields

[TABLE]

Moreover, by strong duality, we have

[TABLE]

The above two displays yield $\mathcal{S}(f)-\mathcal{S}(f_{0})\geq(8M)^{-1}\|\nabla f-\nabla f_{0}\|_{L^{2}(P)}^{2}$. In the same way, using the lower bound in (5.13), we get that $\mathcal{S}(f)-\mathcal{S}(f_{0})\leq 2M\|\nabla f-\nabla f_{0}\|_{L^{2}(P)}^{2}$, which concludes the proof of (5.11).

It turns out that (5.12) is even easier to prove. Indeed, by Proposition-Definition 7(ii) and Lemma 17, we get that the upper bound in (5.13) is also true for $f^{*}$ on all of $\mathbbm{R}^{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].

5.4 Control of the stochastic error via empirical processes

In light of Proposition 10, the performance of our estimator $\hat{T}_{J}=\nabla f_{J}$ defined in (5.8) requires the control of $\mathcal{S}(\hat{f}_{J})-\mathcal{S}(f_{0})$, which can be achieved using tools from empirical process theory. To that end, for any $f$, define

[TABLE]

and let $\bar{f}_{J}\in\mathcal{F}_{J}$. We observe that by optimality of $\hat{f}_{J}$ for $\hat{\mathcal{S}}$,

[TABLE]

To proceed, we control the localized empirical process

[TABLE]

for some fixed $\tau^{2}>0$. More precisely, we prove the following result in Appendix F.2.

Proposition 11**.**

Let assumptions B1 – B2 be fulfilled and define $\mathcal{F}_{J}$ as in (5.7). For any $\tau>0$, define

[TABLE]

Then, there exists $C_{1}=C_{1}(d,\Omega_{P},\tilde{\Omega}_{P},\Omega_{Q},M)>0$ such that with probability at least $1-\exp(-t)$,

[TABLE]

where

[TABLE]

Equipped with this result, we can apply van de Geer’s localization technique. To simplify the presentation, assume that $\bar{f}=\bar{f}_{J}\in\mathop{\mathrm{argmin}}_{f\in\mathcal{F}_{J}}\mathcal{S}(f)$ exists. If not, we may repeat the proof with an $\varepsilon$-approximate minimizer and let $\varepsilon\to 0$. Throughout the proof, we write $\hat{f}=\hat{f}_{J}$ and $\|\cdot\|=\|\cdot\|_{L^{2}(P)}$.

Fix $\sigma>0$ to be defined later and set

[TABLE]

Note that since $s\in[0,1]$ and $\mathcal{F}_{J}$ is convex, we have $\hat{f}_{s}\in\mathcal{F}_{J}$.

On the one hand, $\hat{f}_{s}$ is localized in the sense that

[TABLE]

By Proposition 10 and the triangle inequality respectively, this yields

[TABLE]

Therefore, $\hat{f}_{s}\in\mathcal{F}_{J}(\tau^{2})$. For the same reason, we also have that $\bar{f}\in\mathcal{F}_{J}(\tau^{2})$.

On the other hand, $\hat{f}_{s}$, akin to $\hat{f}$, has empirical risk smaller than $\bar{f}$. Indeed, by convexity of $\hat{\mathcal{S}}$ and the fact that $\hat{f}$ minimizes $\hat{\mathcal{S}}$ over $\mathcal{G}$, we obtain

[TABLE]

which yields

[TABLE]

Together with Jensen’s inequality and Proposition 10 respectively, the above display yields

[TABLE]

Next, note for $s$ as in (5.22), we have that $\|\nabla\hat{f}_{s}-\nabla\bar{f}\|\geq\sigma/2$ iff $\|\nabla\hat{f}-\nabla\bar{f}\|\geq\sigma$. Hence

[TABLE]

Recalling Proposition 11, we take $\sigma$ such that

[TABLE]

so that we get

[TABLE]

In particular, we can check that (5.26) is fulfilled if we choose $\sigma$ such that

[TABLE]

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

[TABLE]

Moreover, integrating the tail with respect to $t$ readily yields by Fubini’s theorem that

[TABLE]

We have proved the following result.

Proposition 12**.**

Let B1 – B2 hold and define $\mathcal{F}_{J}$ as in (5.7). Then, writing

[TABLE]

the estimator $\hat{T}_{J}$ defined in (5.8) satisfies

[TABLE]

Moreover, with probability at least $1-\exp(-t)$,

[TABLE]

5.5 Control of the approximation error

Next, we control the approximation error $\inf_{f\in\mathcal{F}_{J}}\|\nabla f-\nabla f_{0}\|_{L^{2}(P)}$ that appears in Proposition 12. In fact, it is sufficient to control $\|\nabla\bar{f}-\nabla f_{0}\|_{L^{2}(P)}$ where $\bar{f}=\Pi_{J}\operatorname{ext}f_{0}$ is the truncation of $f_{0}$ to its first $J$ wavelet scales after extending $f_{0}$ to all $\mathbbm{R}^{d}$. In light of Theorem 23, we may assume that $\operatorname{ext}\bar{f}$ has the same $C^{\beta}$- and $W^{\alpha,2}$-norm as $\bar{f}$ up to a constant depending on $\tilde{\Omega}_{P}$.

To control the approximation associated with truncating a wavelet decomposition, we rely on the following lemma for Besov functions.

Lemma 13**.**

Let $f\in B_{p,q}^{s}(\mathbbm{R}^{d})$ and denote by $\Pi_{J}$ its projection onto the first scale $J$ wavelet coefficients. That is, if

[TABLE]

where $\Psi^{j,g}_{k}$ are multi-dimensional Daubechies wavelets and $G^{j}$ the associated index sets as in Section B. Then, for all $1\leq p,q\leq\infty$, $s\geq 0$,

[TABLE]

Moreover, for every $q^{\prime}\in[1,\infty]$, $s^{\prime}>0$, and $1\leq p\leq p^{\prime}\leq\infty$ such that $s-d/p>s^{\prime}-d/p^{\prime}$,

[TABLE]

In particular: If $f\in C^{\alpha+1}$ for $\alpha>1$, then $\|f-\Pi_{J}f\|_{C^{2}}\lesssim 2^{-J(\alpha-1)}\|f\|_{C^{\alpha+1}}$ and if $f\in W^{\alpha+1,2}$, for $\alpha>0$, then $\|f-\Pi_{J}f\|_{W^{1,2}}\lesssim 2^{-J\alpha}\|f\|_{W^{\alpha+1,2}}$.

Proof.

Write $\gamma$ for the wavelet coefficients of $f$. The statements (5.31) and (5.32) follow immediately from the wavelet characterization of the Besov norms, (B.11).

To prove the remaining statements, note that for every $j$, because $\|\,.\,\|_{\ell^{p^{\prime}}}\leq\|\,.\,\|_{\ell^{p}}$ and $\|\,.\,\|_{q}$ and $\|\,.\,\|_{q^{\prime}}$ are comparable up to a constant due to $|G^{j}|\leq 2^{d}$ being finite,

[TABLE]

Then, plugging this into the wavelet expansion of $\Pi_{J}f-f$, we obtain

[TABLE]

Finally, to obtain the special cases, note that $\|\,.\,\|_{B^{s}_{\infty,\infty}}\lesssim\|\,.\,\|_{C^{s}}\lesssim\|\,.\,\|_{B^{s}_{1,\infty}}$ and $\|\,.\,\|_{W^{s,2}}=\|\,.\,\|_{B^{s}_{2,2}}$ by Theorem 22. ∎

The above lemma together with Proposition-Definition 7 allows us to check that $\bar{f}\in\mathcal{F}_{J}$. Indeed, by Weyl’s inequality, we have for any $x\in\tilde{\Omega}_{P}$ that

[TABLE]

It follows from Lemma 13 that $\|\Pi_{J}\operatorname{ext}f_{0}-f_{0}\|_{C^{2}(\tilde{\Omega}_{P})}\lesssim 2^{-(\beta-1)J}(M^{3}+R)\leq 1/(2M)$, if

[TABLE]

and $C_{3}=C_{3}(d,\beta,\tilde{\Omega}_{P})$ is large enough. This yields $\lambda_{\min}(D^{2}\Pi_{J}\operatorname{ext}f_{0}(x))\geq 1/(2M)$ and $\bar{f}$ is strongly convex. Similarly, we get $\lambda_{\max}(D^{2}\Pi_{J}\operatorname{ext}f_{0}(x))\leq 2M$ and hence that $\bar{f}\in\mathcal{F}_{J}$. Thus,

[TABLE]

where we used Assumption B2 and Lemma 13.

We have thus proved that

[TABLE]

5.6 Bias-variance tradeoff

We are now in a position to complete the proof of Theorem 8. Combining the bounds (5.29) and (5.45), we get

[TABLE]

We conclude the proof by optimizing with respect to $J$. It yields

[TABLE]

We note that since $\alpha>1$, in the first two cases, $d\in\{1,2\}$, the cut-off $J$ can be picked independently from $\alpha$. Finally, high-probability bounds can be obtained in a similar manner.

6 Numerical experiments

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, $\hat{T}_{\mathrm{wav}}$ below, is an approximation to the estimator $\hat{T}_{J}$ in (5.8), illustrating that (5.8) can be implemented in low dimensions and that this approximation achieves favorable rates in $d=3$. The second, $\hat{T}_{\mathrm{ker}}$ 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.

6.1 Estimators

6.1.1 Baseline estimator

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

[TABLE]

respectively, and calculate the optimal transport matrix $\Gamma\in\mathbbm{R}^{n\times n}$

[TABLE]

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

[TABLE]

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 $\hat{\Gamma}$ is in fact a (rescaled) permutation matrix, leading to a matching $\hat{\pi}:[n]\to[n]$, and hence to $\hat{T}_{\mathrm{emp}}(X_{i})=Y_{\pi(i)}$.

Because $\hat{T}_{\mathrm{emp}}(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:

[TABLE]

6.1.2 Wavelet estimator

Next, we turn to an approximation of (5.8). Assume that in addition to a superset $\tilde{\Omega}_{P}\supseteq\Omega_{P}$, we are also given a superset $\tilde{\Omega}_{Q}\supseteq\Omega_{Q}$, and that both $\tilde{\Omega}_{P}$ and $\tilde{\Omega}_{Q}$ are boxes (hypercubes). We consider all functions originally defined over $\tilde{\Omega}_{P}$ and $\tilde{\Omega}_{Q}$ as given by their samples on grids with resolution $N\in\mathbb{N}$, $\bm{x}=(x_{i})_{i\in[N]^{d}}\in(\mathbbm{R}^{d})^{N^{d}}$ and $\bm{y}=(y_{i})_{i\in[N]^{d}}\in(\mathbbm{R}^{d})^{N^{d}}$, respectively. In particular, we write $\bm{f}=(f_{i})_{i\in[N]^{d}}\in\mathbbm{R}^{N^{d}}$ for the discretization of the potential $f$ and $\bm{T}=(T_{i})_{i\in[N]^{d}}\in(\mathbbm{R}^{d})^{[N]^{d}}$ for the discretization of the transport map $T$ on the grid $\bm{x}$. Here, we pick $N=65$.

We employ the following discretization/approximation schemes:

1. (i)

The restrictions to functions up to wavelet scale $J\in\mathbb{N}$ can be obtained by parametrizing $\bm{f}$ by the inverse discrete wavelet transform up to order $J$, which we write as $\bm{f}=W^{\top}_{J}\gamma_{J}$ for wavelet coefficients $\gamma_{J}\in\mathbbm{R}^{m_{J}}$ with $m_{J}\in\mathbb{N}$. For these experiments, we use db4 Daubechies wavelets, i.e., Daubechies wavelets with four vanishing moments. 2. (ii)

An approximation to the (continuous) Legendre transform is given by the discrete Legendre transform,

[TABLE]

for $\bm{f}\in(\mathbbm{R}^{d})^{N^{d}}$ and $j\in[N]^{d}$. 3. (iii)

The gradient of $\bm{f}$ on grid-points can be calculated by a finite-difference scheme, which we write as $\nabla_{\bm{x}}\bm{f}$. 4. (iv)

To obtain values of these approximations at non-grid points, we appeal to linear interpolation, which we write as $\mathcal{P}_{\bm{x}\to\{X_{1},\dots,X_{n}\}}\bm{f}$ for the interpolated value of $\bm{f}$ defined on the grid $\bm{x}$ at the point $\{X_{1},\dots,X_{n}\}\in\mathbbm{R}^{d}$, collected in one vector.

Given *i.i.d. *observations $X_{1:n}=\{X_{1},\dots,X_{n}\}$ and $Y_{1:n}=\{Y_{1},\dots,Y_{n}\}$ from $P$ and $Q$, respectively, we arrive at the following optimization problem as approximation to (5.8):

[TABLE]

Note that we dropped all boundedness and convexity constraints that were given by $\mathcal{X}(2M)$ in (5.8). In practice, this can lead to degraded estimation quality of the gradient of $W^{\top}_{J}\gamma_{J}$ near the boundary of $\Omega_{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

[TABLE]

With this, we denote by $\hat{T}_{\mathrm{wav}}^{(J)}$ the function obtained by linearly interpolating $\hat{\bm{T}}_{J}$,

[TABLE]

For the purpose of these experiments, we select the wavelet scale $\hat{J}$ by an oracle choice, i.e., as the minimizer of an approximation to the population semi-dual (2.14), see Section H.1, while in practice, one would resort to cross-validation methods for this purpose. Finally, we set

[TABLE]

To obtain an error measure for $\hat{T}_{\mathrm{wav}}$ that is easily comparable to $\mathrm{MSE}_{n}(\hat{T}_{\mathrm{emp}})$, we consider the empirical $L^{2}(\hat{P})$ norm on the sample points,

[TABLE]

Note that the objective function in (6.7) can be calculated in linear time with respect to the underlying grid, that is, in $O(N^{d})$, thanks to efficient 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 $\hat{\gamma}_{J}$, even though the objective is not smooth for every $\gamma_{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 $\hat{T}_{\mathrm{wav}}$ 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 $\hat{T}_{\mathrm{wav}}$ to the estimator $\hat{T}_{J}$ considered in Section 5. We remark, however, that our approximations closely follow the definition of $\hat{T}_{J}$ and that the omitted constraints defining the set $\mathcal{X}(2M)$ could be incorporated into the approximation by means of a finite-difference discretization as well, albeit at an additional computational cost.

6.1.3 Kernel smoothing estimator

While the estimator $\hat{T}_{\mathrm{wav}}$ 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 $\hat{T}_{\mathrm{emp}}$. 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 $\hat{T}_{\mathrm{emp}}$ by smoothing it via kernel-ridge regression [Mur12].

Let $\mathcal{H}$ denote a reproducing kernel Hilbert space (RKHS) with associated kernel $k(x,y)$, $x,y\in\mathbbm{R}^{d}$, and norm $\|T\|_{\mathcal{H}}$ for $T\in\mathcal{H}$. Here, we consider the RKHS given by Gaussian radial basis functions,

[TABLE]

We fit an RKHS function $T$ to the pairs $(X_{i},\tilde{Y}_{i}=\hat{T}_{\mathrm{emp}}(X_{i}))$ by solving the regularized kernel regression objective

[TABLE]

for $\nu_{\mathrm{ridge}}>0$. By the representer theorem [Mur12], (6.13) has a solution

[TABLE]

where

[TABLE]

We measure its performance by

[TABLE]

Similar to $\hat{T}_{\mathrm{wav}}$, we select the tuning parameters $\nu_{\mathrm{kernel}}$ and $\nu_{\mathrm{ridge}}$ by an oracle procedure, picking those parameters from a finite grid that minimize (6.16) on an independent hold-out sample $\tilde{X}_{1:n}$, and denote the resulting estimator by $\hat{T}_{\mathrm{ker}}$.

6.2 Setup

For $d\in\mathbb{N}$, we consider the following examples of smooth potentials and transport maps:

[TABLE]

where for (id), $P^{(1)}=Q^{(1)}=\mathsf{Unif}([0,1]^{d})$. For (exp), $P^{(2)}=\mathsf{Unif}([0,1]^{d})$, and the target measure is defined as $Q^{(2)}=(T_{0}^{(2)})_{\#}P^{(2)}$. Note that these potentials and transport maps are $C^{\infty}$ and strongly convex on any compact convex subset of $\mathbbm{R}^{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\in\{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.13) with off-the-shelf methods scales with $O(n^{3})$. For the sake of this comparison, we simply restrict our experiments on $\hat{T}_{\mathrm{emp}}$ and $\hat{T}_{\mathrm{ker}}$ to sample sizes $\leq 10^{4}$. Likewise, we do not compute $\hat{T}_{\mathrm{wav}}$ in $d=10$ due to the large computational cost.

6.3 Results

6.3.1 Qualitative comparison in 2D

To give a qualitative idea of the considered estimators compared to the baseline $\hat{T}_{\mathrm{emp}}$, we visualize the first coordinate of the transport map estimators for case (id) with $d=2$ and $n=100$ observations from $P=Q=\mathsf{Unif}([0,1]^{2})$ in Figure 1, together with the ground truth transport map $(T_{0}^{(1)})_{1}$. To depict the coupling $\hat{T}_{\mathrm{emp}}$ 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 1(c) produces a visibly smoother map compared to the unregularized $\hat{T}_{\mathrm{emp}}$ in Figure 1(b). The kernel estimator in Figure 1(d) is even smoother and visually very similar to the ground truth in Figure 1(a), due to the possibility of employing a large amount of regularization.

6.3.2 Quantitative comparison in 3D and 10D

To obtain a quantitative comparison, for both test cases and $d\in\{3,10\}$, we compute $\hat{T}_{\mathrm{emp}}$, $\hat{T}_{\mathrm{wav}}$ (only $d=3$), and $\hat{T}_{\mathrm{ker}}$ 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 $\hat{T}_{\mathrm{emp}}$ roughly follow a $n^{-2/3}$ rate. This corresponds to the decay of the average $\ell^{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}(\hat{P},P)$. We are, however, not aware of a generalization of this rate to the error measure $\mathrm{MSE}$ in the case of a transport map different from the identity.

The error curves for the wavelet estimator $\hat{T}_{\mathrm{wav}}$ 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 $\hat{T}_{\mathrm{wav}}$ all lie below the baseline estimator. The kernel estimator $\hat{T}_{\mathrm{ker}}$ performs even better, attaining rates close to $n^{-1}$ and yielding consistently better rates than $\hat{T}_{\mathrm{wav}}$.

We observe that the favorable behavior of $\hat{T}_{\mathrm{wav}}$ suggests that the restriction of candidate potentials to $\mathcal{X}(2M)$ in (5.8) might not be necessary and could possibly be omitted.

In 10D, for both test cases, $\hat{T}_{\mathrm{emp}}$ 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 $\hat{T}_{\mathrm{wav}}$ and $\hat{T}_{\mathrm{ker}}$.

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 $\hat{T}_{\mathrm{ker}}$ under regularity assumptions on $T_{0}$, for which the minimax rates established in this work can serve as a benchmark.

Acknowledgments. The authors would like to thank Richard Nickl for pointing out relevant references on the estimation of invariant measures in the multivariate Gaussian white noise model, and Jonathan Niles-Weed for suggesting references on the rate between random matchings. The authors are also indebted to the Associate Editor and three referees for providing insightful comments that greatly improved the quality of the paper.

Philippe Rigollet was supported by NSF awards IIS-1838071, DMS-1712596, DMS-TRIPODS- 1740751, and ONR grant N00014-17- 1-2147.

Appendix A Convex analysis

In this section, we recall some useful facts from convex analysis. We refer the reader to [HL01] for a comprehensive treatment.

Recall that a set $U\subseteq\mathbbm{R}^{d}$ is convex if for all $x,y\in U$, $t\in[0,1]$, $tx+(1-t)y\in U$. A function $f\colon U\to\mathbbm{R}\cup\{+\infty\}$ is convex if for all $x,y\in U$, $t\in[0,1]$, it holds that

[TABLE]

Moreover, we call a function $\mu$-strongly convex if for all $x,y\in U$, $t\in[0,1]$, we have

[TABLE]

For twice differentiable functions, it is often more convenient to employ the following analytic criterion for strong convexity.

Lemma 14** ([HL01, Theorem B.4.3.1]).**

Let $U\subseteq\mathbbm{R}^{d}$ be an open convex set and $f\colon U\to\mathbbm{R}$ twice differentiable. Then, $f$ is $\mu$-strongly convex if and only if $\lambda_{\min}(D^{2}f(x))\geq\mu$ for all $x\in U$.

Note that even if $U\subsetneq\mathbbm{R}^{d}$ is a proper convex subset of $\mathbbm{R}^{d}$, we can always consider a function $f\colon U\to\mathbbm{R}\cup\{+\infty\}$ to be defined on all of $\mathbbm{R}^{d}$ by setting it to $+\infty$ outside of $U$. To that end, let $\iota_{U}$ be the indicator function defined by

[TABLE]

We define the extension of $f$ outside $U$ by $f+\iota_{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 $\operatorname{dom}(f)=\{x\in\mathbbm{R}^{d}:f(x)<+\infty\}$ the domain of a convex function.

We now recall two important notions associated with convex functions $f$.

First, the subdifferential of $f$ at $x\in\operatorname{dom}(f)$ is defined as

[TABLE]

As indicated by the following lemma, the subdifferential reduces to the gradient for differentiable functions.

Lemma 15** ([HL01, Corollary D.2.1.4]).**

Let $f\colon\mathbbm{R}^{d}\to\mathbbm{R}\cup\{+\infty\}$ be a convex function. If $f$ is differentiable at $x\in\mathbbm{R}^{d}$ with gradient $\nabla f(x)$, then $\partial f(x)=\{\nabla f(x)\}$.

Conversely, if $\partial f(x)=\{a\}$ consists of only a single element, then $f$ is differentiable at $x$ with gradient $\nabla f(x)=a$.

Second, the convex conjugate, or Legendre-Fenchel conjugate, is defined for any function $f\colon\mathbbm{R}^{d}\to\mathbbm{R}\cup\{+\infty\}$ as

[TABLE]

By considering $f+\iota_{U}$, this definition extends to functions $f\colon U\to\mathbbm{R}\cup\{+\infty\}$.

We recall the following standard facts about the convex conjugate, stated here without proof (see [HL01, Part E] for details).

Lemma 16**.**

If $f\colon\mathbbm{R}^{d}\to\mathbbm{R}\cup\{+\infty\}$ is convex and lower semi-continuous, then $f^{\ast\ast}=f$.

Lemma 17**.**

Let $U$ be a closed, convex set. If $f\colon U\to\mathbbm{R}$ is $\mu$-strongly convex, then $\operatorname{dom}(f^{\ast})=\mathbbm{R}^{d}$, and $\nabla f^{\ast}$ is $\mu^{-1}$-Lipschitz:

[TABLE]

Lemma 18**.**

Denote by $f\colon\mathbbm{R}^{d}\to\mathbbm{R}\cup\{+\infty\}$ a lower semi-continuous and convex function. Then, for $x,y\in\mathbbm{R}^{d}$,

[TABLE]

Lemma 19**.**

If $f,g\colon U\to\mathbbm{R}$, $f(x)\leq g(x)$ for all $x\in U$, then $f^{\ast}(y)\geq g^{\ast}(y)$ for all $y\in\mathbbm{R}$.

Lemma 20**.**

Given two functions $f,g:U\to\mathbbm{R}$ for a compact set $U$. Abusing notation, denote the convex conjugate of the function $f+\iota_{U}$ by $f^{\ast}$. Then, $\|f^{\ast}-g^{\ast}\|_{L^{\infty}(\mathbbm{R}^{d})}\leq\|f-g\|_{L^{\infty}(U)}$.

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\in\mathbbm{R}^{d}$, $c\in\mathbbm{R}$ and let $U\subseteq\mathbbm{R}^{d}$ be a closed, convex set. Define

[TABLE]

Then,

[TABLE]

where $d^{2}$ denotes the squared distance $d^{2}(x,U)=\inf_{y\in U}\|x-y\|_{2}^{2}$.

Proof.

Note first that for any $y\in\mathbbm{R}^{d}$,

[TABLE]

where

[TABLE]

Moreover,

[TABLE]

Writing

[TABLE]

we see that the infimum in (A.9) is achieved by the projection $\bar{x}$ of $(y-b)/a$ onto the closed convex set $U+t$. Moreover, the value of the objective at $\bar{x}$ is given by

[TABLE]

Together with (A.8) and (A.9), this completes the proof of the lemma. ∎

Appendix B Wavelets and function spaces

In this section, we give a brief overview of the different function spaces used in the paper and how their norms can be related to their wavelet coefficients.

First, we recall the definition of Hölder and Sobolev spaces. Let $\Omega\subseteq\mathbbm{R}^{d}$ be a closed set with non-empty interior and denote by $C_{u}(\Omega)$ the set of uniformly continuous functions on $\Omega$. The Hölder norms for any $f\in C_{u}(\Omega)$ are defined as follows. For any integer $k\geq 0$ and $f$ that admits continuous derivatives up to order $k$, define

[TABLE]

and for any real number $\alpha>0$, define

[TABLE]

The space $C^{\alpha}(\Omega)$ is then defined as the set of functions for which this norm is finite. For a vector-valued function $T\colon\Omega\to\mathbbm{R}^{d}$, $T=(T_{1},\dots,T_{d})^{\top}$, we similarly define the norms as the sum over the individual norms,

[TABLE]

Similarly, for an integer $k\geq 0$ and $p\in[1,\infty]$, the Sobolev norms are defined as

[TABLE]

where the derivative $\partial^{\flat}$ is to be understood in the sense of distributions and the Sobolev space $W^{m,2}(\Omega)$ is the space of all functions for which this norm is finite. This definition can be extended to $\alpha>0$, for example by defining $W^{\alpha,2}(\Omega)$ as the Besov space $B^{\alpha}_{2,2}(\Omega)$, 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 $\psi_{\mathfrak{M}}\in C^{r}(\mathbbm{R})$ and $\psi_{\mathfrak{F}}\in C^{r}(\mathbbm{R})$ a compactly supported wavelet and scaling function, respectively, for example Daubechies wavelets. This implies that

[TABLE]

is an orthonormal basis of $L^{2}(\mathbbm{R})$. To obtain a basis of $L^{2}(\mathbbm{R}^{d})$, for $j\in\mathbb{N}$, set

[TABLE]

and for $g\in G^{j}\cup G^{0}$,

[TABLE]

This gives the orthonormal basis

[TABLE]

Wavelet coefficients, defined as the expansion coefficients with respect to the above basis for $L^{2}(\mathbbm{R}^{d})$ functions, can be used to characterize the so-called Besov spaces: Let $1\leq p,q\leq\infty$, $s\geq 0$, and let the regularity of the above wavelets satisfy

[TABLE]

With this, if $f$ admits the wavelet representation

[TABLE]

we define the family of norms $\|\,.\,\|_{B^{s}_{p,q}}$ by

[TABLE]

In particular, by the orthonormality of the wavelets $\{\Psi_{k}^{j,g}\}_{j,g,k}$,

[TABLE]

For a bounded Lipschitz domain $\Omega\subseteq\mathbbm{R}^{d}$, we can define the Besov spaces $B^{s}_{p,q}(\Omega)$ by restrictions of functions on $\mathbbm{R}^{d}$ with norm

[TABLE]

Note that since we work with compactly supported wavelets, for a function $f$ with compact support, only a finite number of wavelet coefficients are non-vanishing in (B.10). In particular, for non-zero wavelet coefficients are contained in a set $\Lambda(j)$ with $|\Lambda(j)|\lesssim 2^{jd}$.

The following theorem collects some basic properties of Besov spaces and their relationship to Hölder and Sobolev spaces.

Theorem 22**.**

Let $\Omega\subseteq\mathbbm{R}^{d}$ be a bounded Lipschitz domain.

1. (i)

[GN16, Proposition 4.3.6]** Let $s,s^{\prime}>0$, $1\leq p,p^{\prime},q,q^{\prime}\leq\infty$. Then, the following inclusions hold in the sense of continuous embeddings:

1. (a)

$B^{s}_{p,q}\subseteq B^{s}_{p,q^{\prime}}$, if $q\leq q^{\prime}$, 2. (b)

$B^{s}_{p,q}\subseteq B^{s^{\prime}}_{p,q^{\prime}}$, if $s>s^{\prime}$, 2. (ii)

[Tri06, Theorem 1.122]** Let $k\in\mathbb{N}$. Then, $W^{k,2}(\Omega)=B_{2,2}^{k}(\Omega)$ and we define $W^{\alpha,2}(\Omega):=B_{2,2}^{\alpha}(\Omega)$ for $\alpha>0$ not an integer. 3. (iii)

[Tri06, Theorem 1.122]**, **[GN16, Proposition 4.3.20]** Let $\alpha>0$. If $\alpha$ is not integer, then

[TABLE]

If $\alpha$ is integer, then

[TABLE]

Note that the proof in [GN16] of (B.16) is given in 1D, but extends naturally to arbitrary dimensions by summability of the wavelet coefficients.

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 coefficients of their extensions.

Theorem 23** (Extension operator, [Tri06, Theorem 1.105], [Ryc99]).**

Let $\Omega\subseteq\mathbbm{R}^{d}$ be a bounded Lipschitz domain. Then, there exists a linear extension operator $\operatorname{ext}$ that preserves $L^{2}$-, $C^{\beta}$-, and $W^{\alpha,2}$-norms. That is, there exists an extension operator $\operatorname{ext}=\operatorname{ext}(\Omega)$ such that for $A\in\{L^{2},C^{\beta},W^{\alpha,2}:\beta>0,\alpha>0\}$, there exist constants $C=C(\Omega,A)$ with

[TABLE]

We conclude this section by a lemma that provides uniform control of a function by its wavelet coefficients. It is useful to control bracketing entropy numbers.

Lemma 24**.**

Let $f\in V_{J}(\mathbbm{R}^{d})$ with wavelet coefficients $\gamma^{j,g}_{k}$ for compactly supported mother and father wavelets. Then, $\|f\|_{\infty}\lesssim 2^{Jd/2}\|\gamma\|_{\infty}\leq 2^{Jd/2}\|\gamma\|_{2}.$

Proof.

Let $x\in\mathbbm{R}^{d}$ and write

[TABLE]

where we used Hölder’s inequality and the fact that only a finite number of $k$ enters the summation for each fixed $x\in\mathbbm{R}^{d}$. The last inequality is trivial. ∎

Appendix C Metric entropy and suprema of stochastic processes

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 $\mathcal{F}$ with respect to a norm $\|.\|$ by $N_{[\,]}(\delta,\mathcal{F},\|\,.\,\|)$ and define the Dudley integral as

[TABLE]

Theorem 25** (Bernstein chaining, [vW07, Lemma 3.4.2] ).**

Let $\mathcal{F}$ be a class of measurable functions such that $\mathbbm{E}[f^{2}]<\sigma^{2}$ and $\|f\|_{\infty}\leq M$ for all $f\in\mathcal{F}$. Then,

[TABLE]

Theorem 26** (Concentration, see [Mas07, Equation (5.50)], going back to [Bou02]).**

Under the same assumptions as Theorem 25,

[TABLE]

Lemma 27** (Bracketing numbers from $L^{\infty}$ covering numbers).**

Let $P$ be a probability measure and $A\subseteq\mathcal{F}$ be a set of functions with $N(A,L^{\infty}(P),\delta/2)\leq\phi(\delta)$. Then, $N_{[]}(A,L^{2}(P),\delta)\leq\phi(\delta)$.

Moreover, if for every function in the original class $f\in\mathcal{F}$, $\mathbbm{E}[f^{2}]<\sigma^{2}$ and $\|f\|_{L^{\infty}(P)}<M$, then every bracket $[f_{1},f_{2}]$ in a $\delta$-cover above (note that $f_{1}$ and $f_{2}$ need not be members of $\mathcal{F}$) satisfies

[TABLE]

Proof.

Denote by $\{f_{1},\dots,f_{N}\}$ the centers of a minimal $\delta/2$-covering of $\mathcal{F}$ in $L^{\infty}(P)$. Let $i\in\{1,\dots,N\}$. Then, each $\delta/2$ ball around $f_{i}$ is contained in the bracket $[f_{i}-\delta/2,f_{i}+\delta/2]$. Moreover, the $L^{2}(P)$ diameter of the above bracket is bounded by its $L^{\infty}(P)$ diameter, which is $\delta$, so the collection of those brackets yields the desired covering with brackets.

The rest of the lemma follows from

[TABLE]

∎

The following lemma can be shown by directly specifying a grid or by a volume argument such as [Ver18, Proposition 4.2.12].

Lemma 28** (Covering numbers for norm balls).**

Fix $p\in\mathbb{N}$ and denote by $B_{\infty}(A)$ the $\ell^{\infty}$ ball of $\mathbbm{R}^{d}$ with radius $A$. Then, $N(B_{\infty}(A),\|\,.\,\|_{\infty},\delta)\leq\left({3A}/{\delta}\right)^{d}$.

Lemma 29** ([Rau10, Lemma 10.3]).**

For $\alpha>0$,

[TABLE]

Appendix D Tools for lower bounds

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\geq 2$, $\Theta=\{T_{0},\dots,T_{K}\}$ a collection of hypotheses, and let $d$ be a pseudometric, i.e., a bi-variate function on $\Theta$ such that

1. (a)

$d(T_{j},T_{k})\geq 0$; 2. (b)

$d(T_{j},T_{j})=0$; 3. (c)

$d(T_{j},T_{k})=d(T_{k},T_{j})$; 4. (d)

$d(T_{j},T_{\ell})\leq d(T_{j},T_{k})+d(T_{k},T_{\ell})$;

for all $0\leq j,k,\ell\in\{0,\dots,K\}$.

Suppose $\Theta$ fulfills:

1. (i)

$d(T_{j},T_{k})\geq 2s>0$, for all $0\leq j<k\leq K$; 2. (ii)

$P_{j}\ll P_{0}$, for all $j\in[K]$, and

[TABLE]

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,\dots,K$.

Then,

[TABLE]

Lemma 31** (Varshamov-Gilbert lemma).**

[Tsy09, Lemma 2.9]** Let $D\geq 8$. There exists a subset $\tau^{(0)},\dots,\tau^{(K)}$ of $\{0,1\}^{D}$ such that $\tau^{(0)}=(0,\dots,0)$,

[TABLE]

and

[TABLE]

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\in\{0,1\}$. Denoting $s=d(T_{0},T_{1})/2$, if

[TABLE]

then

[TABLE]

Appendix E Alternative assumptions via Caffarelli’s regularity theory

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 Caffarelli’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 $\Omega$ with $C^{2}$ boundary as uniformly convex if it can be written as the sublevel set, $\Omega=\{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 $\Omega$ on its boundary $\partial\Omega$ [Vil09].

C1****Supports.

Assume $\Omega_{P},\Omega_{Q}$ are two uniformly convex, bounded, open subsets of $\mathbbm{R}^{d}$ with $C^{\lfloor\alpha-1\rfloor+2}$ boundaries.

C2****Densities.

With $\Omega_{P},\Omega_{Q}$ as in Assumption C1, assume that $\rho_{P}\in C^{\alpha-1}(\overline{\Omega}_{P})$ and $\rho_{Q}\in C^{\alpha-1}(\overline{\Omega}_{Q})$ are two probability densities on $\mathbbm{R}^{d}$ with $\operatorname{supp}(\rho_{P})=\overline{\Omega}_{P}$, $\operatorname{supp}(\rho_{Q})=\overline{\Omega}_{Q}$ that are bounded from above and below on their support. Further, denote the associated probability distributions by $P$ and $Q$, respectively.

Theorem 33** (Caffarelli’s global regularity theory).**

[Vil09, Theorem 12.50(iii)]**, [Caf96]. Let $\alpha>1$ and assume that Assumptions C1 and C2 hold. Then, for the optimal transport potential $f_{0}$, i.e., the solution to (2.14) for $P$ and $Q$, which is unique up to an additive constant, it holds that $f_{0}\in C^{\alpha+1}(\overline{\Omega_{P}})$.

Moreover, we need the following extension lemma to smoothly extend the potential (and associated transport map) to a larger set as required by Assumption B2. Similar statements have been shown in [Gho02, Yan14, AM19].

Lemma 34**.**

Let $\Omega$ be a convex compact subset of $\mathbbm{R}^{d}$, $f:\Omega\to\mathbbm{R}$ a convex $C^{\alpha}(\Omega)$ function for $\alpha>2$ with

[TABLE]

Then, there exists an $\epsilon>0$ and an extension $\tilde{f}$ of $f$ to $\Omega_{\epsilon}=\Omega+\epsilon B(0,1)$ such that $f\in C^{\alpha}(\Omega_{\epsilon})$ and

[TABLE]

Proof.

The claim follows from extending $f$ via Theorem 23 to all of $\mathbbm{R}^{d}$ while preserving the Hölder class $C^{\alpha}$, and observing that the strong convexity condition (E.2) on an enlargement of $\Omega$ can be ensured by uniform continuity of the extension. ∎

Corollary 35**.**

Let $\alpha>1$. Under Assumptions C1 and C2, writing $f_{0}$ and $T_{0}=\nabla f_{0}$ for the Kantorovich potential and optimal transport map between $P$ and $Q$, respectively, there exists an estimator $\hat{T}$ such that

[TABLE]

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 $\overline{\Omega}_{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 $\rho_{P}$ and $\rho_{Q}$. By Theorem 33, $f_{0}\in C^{\alpha+1}(\overline{\Omega}_{P})$ and thus $T_{0}=\nabla f_{0}\in C^{\alpha}(\overline{\Omega}_{P},\mathbbm{R}^{d})$. Denoting by $\tilde{f}_{0}$ the extension of $f_{0}$ to $\tilde{\Omega}_{P}=(\overline{\Omega}_{P})_{\epsilon}$ from Lemma 34, we observe that $\tilde{\Omega}_{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 $\tilde{T}_{0}=\tilde{f}_{0}$ and the compactness of $\tilde{\Omega}_{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.7) 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^{\alpha-1}$-norms of $\rho_{P}$ and $\rho_{Q}$ alone. Although the statement of Theorem 33 strongly suggests a bound on $\|f_{0}\|_{C^{\alpha+1}(\overline{\Omega}_{P})}$ in terms of $\|\rho_{P}\|_{C^{\alpha-1}(\overline{\Omega}_{P})}$ and $\|\rho_{Q}\|_{C^{\alpha-1}(\overline{\Omega}_{Q})}$, to the best of our knowledge, such bounds are not available in the literature for the “global” version of Caffarelli’s regularity theory as stated Theorem 33. On the other hand, such bounds hold “locally”, that is, $\|f_{0}\|_{C^{\alpha+1}(U)}$ can be bounded for strictly open subsets of $\overline{\Omega}_{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.

Appendix F Omitted proofs

F.1 Proof of Proposition-Definition 7

We proceed in order of statement of the results.

(i) It follows from B2 that there exists $f_{0}$ such that $\nabla f_{0}=T_{0}$ and by B2(i) that $|\nabla f_{0}(x)|\leq M$ for all $x\in\tilde{\Omega}_{P}$. Moreover, since $f_{0}$ is defined up an additive constant, assume that $f_{0}(x_{0})=0$ for some $x_{0}\in\tilde{\Omega}_{P}$. A first-order Taylor expansion yields that for any $x\in\tilde{\Omega}_{P}$,

[TABLE]

where in the second inequality, we used B2(i), and in the third one, we used that $\operatorname{diam}(\tilde{\Omega}_{P})\leq 2M$ according to B2.

(ii) follows immediately from B2(ii).

Next, as in (i), (5.2) follows from a first-order expansion: For any $x\in\tilde{\Omega}_{P}$, note that on the one hand

[TABLE]

so that $\|f_{0}\|_{L^{\infty}}\leq M^{2}$ and $\|f_{0}\|_{L^{2}}\leq 2M^{3}$. Together with the fact that $T_{0}=\nabla f_{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 $\operatorname{int}(\tilde{\Omega}_{P})$ and can be extended to $+\infty$ outside of $\tilde{\Omega}_{P}$ and thus be also be considered a strongly convex function on $\mathbbm{R}^{d}$. By Lemmas 17 and 18, we conclude that $T_{0}=\nabla f_{0}$ is a bijection from $\Omega_{P}$ onto its image $\Omega_{Q}=\operatorname{supp}((\nabla f_{0})_{\#}P)$, with both $\nabla f_{0}$ and $(\nabla f_{0})^{-1}$ being continuously differentiable; in other words, $\nabla f_{0}$ is a $C^{1}$-diffeomorphism between $\Omega_{P}$ and $\Omega_{Q}$. Since $C^{1}$-diffeomorphisms preserve Lipschitz domains [HMT07, Theorem 4.1] and connectedness, we can conclude that $\Omega_{Q}$ is a connected and bounded Lipschitz domain, and the fact that $\Omega_{Q}\subseteq MB_{1}$ follows from $T(x)\leq M$ for all $x\in\tilde{\Omega}_{P}$.

Finally, we turn to check the condition on the density $\rho_{Q}$. To that end, note that by the change of variables formula, $Q$ has the density

[TABLE]

This readily yields the desired bound in light of B1 which assumes boundedness of $\rho_{P}$ and (ii) which gives boundedness of the Hessian $D^{2}f_{0}$.

F.2 Proof of Proposition 11

We first bound the expectation of the supremum of interest and then obtain the high-probability bound via concentration. To begin, note that

[TABLE]

which yields

[TABLE]

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 $\Omega_{Q}$, which allows us to use wavelet expansions for harmonic analysis.

First, because $(\hat{Q}-Q)(g+c)=(\hat{Q}-Q)g$ for every function $g$ and constant $c\in\mathbbm{R}$, we may assume without loss of generality that

[TABLE]

Note that $f^{*}$ and $f_{0}^{*}$ are defined over all of $\mathbbm{R}^{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 $\operatorname{ext}f^{*}$ (resp. $\operatorname{ext}f^{*}_{0}$) the Lipschitz extension of the restriction of $f^{*}$ (resp. $f^{*}_{0}$) to $\Omega_{Q}$. The existence of a linear extension operator is guaranteed by Theorem 23. In particular, we control the norms of $\operatorname{ext}f^{*}$ and $\operatorname{ext}f^{*}_{0}$.

For a function $g$ with wavelet expansion

[TABLE]

as defined in Section B, we define the two $L^{2}$-projections

[TABLE]

With this, we write $T_{2}=T_{2,1}+T_{2,2}$, where

[TABLE]

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}_{\infty,\infty}(\Omega_{Q})$-norm. In turn, their projections also have bounded $B^{2}_{\infty,\infty}(\Omega_{Q})$-norm, and by Theorem 23, this implies that both $\Pi_{J}\operatorname{ext}(f^{\ast}-f_{0}^{\ast})$ and $\Pi_{>J}\operatorname{ext}(f^{\ast}-f_{0}^{\ast})$ are in $C^{s}(\mathbbm{R}^{d})$, for $s<2$, and hence they are continuous.

Bound on $T_{2,1}$-term.

Recall that it follows from Proposition-Definition 7 that $\Omega_{Q}$ is a connected Lipschitz domain, so we can apply the Poincaré-Wirtinger inequality, Lemma 39, together with (5.12) from Proposition 10, to get

[TABLE]

where we used that we assumed $f^{\ast}-f_{0}^{\ast}$ to be centered.

Hence, $f\in\mathcal{F}_{J}(\tau^{2})$ implies $\|f^{\ast}-f_{0}^{\ast}\|_{W^{1,2}(\Omega_{Q})}\lesssim\tau$, and therefore due to the properties of the extension operator $\operatorname{ext}$,

[TABLE]

for some constant $C_{4}=C_{4}(\Omega_{Q},M)$. Since $\Pi_{J}$ is a non-expansive operator on Besov spaces, it follows from the above display that

[TABLE]

Bounding the empirical process over this standard function class can now be performed as follows. Observe first that for any function $h\in V_{J}(\mathbbm{R}^{d})$ with wavelet decomposition

[TABLE]

the condition $\|h\|_{W^{1,2}(\mathbbm{R}^{d})}\leq C_{4}\,\tau$ is equivalent to

[TABLE]

Next, by symmetrization, for independent copies of Rademacher random variables $\varepsilon_{i}$, $i=1,\dots,n$,

[TABLE]

where both suprema are taken over the set $\mathfrak{G}_{J}=\{h\in V_{J}(\mathbbm{R}^{d}),\;\|h\|_{W^{1,2}(\mathbbm{R}^{d})}\leq C_{4}\tau\}$. To control the Rademacher process, fix $h\in\mathfrak{G}_{J}$ with wavelet decomposition (F.12). By the Cauchy-Schwarz inequality,

[TABLE]

where here and below, all sums without indices are over $\{0\leq j\leq J,\,g\in G^{j},\,k\in\mathbb{Z}^{d}\}$. Since the right-hand side in the display above does not depend on $h$, we get by Jensen’s inequality that

[TABLE]

Since $Y_{i}\in\Omega_{Q}$, a compact set by assumption, for given $j$ and $g$, $\Psi^{j,g}_{k}(Y_{i})$ is non-zero only for $k\in\Lambda(j)$ where $\Lambda(j)$ depends on the diameter of $\Omega_{Q}$, and $|\Lambda(j)|\lesssim 2^{jd}$. Together with the independence of the $\varepsilon_{i}$, this yields

[TABLE]

By Proposition-Definition 7 and the fact that the $\Psi^{j,g}_{k}$ form an orthonormal basis in $L^{2}(\mathbbm{R}^{d})$,

[TABLE]

Thus,

[TABLE]

We have proved that

[TABLE]

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^{\infty}$-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 $\mathcal{F}_{J}(\tau^{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

[TABLE]

and let $f_{1},f_{2}\in\mathcal{F}_{J}(\tau^{2})$ have wavelet coefficients given by the sequences $\gamma_{1}$ and $\gamma_{2}$. Observe that by linearity of the projection and extension operators, we have for any cutoff $J^{\prime}\geq 0$,

[TABLE]

To control the first term in the right-hand side above, in the following lemma, we establish bounds on the potentials in $\tilde{F}_{J}(\tau^{2})$. Its proof is deferred to the next section.

Lemma 37**.**

There exists a constant $C_{5}=C_{5}(\Omega_{P},\tilde{\Omega}_{P},\Omega_{Q},M)$ such that for all $f\in\mathcal{X}(2M)$,

[TABLE]

Moreover, under assumptions B1 – B2, $\|\operatorname{ext}f_{0}^{\ast}\|_{B^{2}_{\infty,\infty}(\mathbbm{R}^{d})}\leq C_{5}/2$ as well.

It follows from Lemma 37 that there exists $C_{5}$ such that for any $f\in\mathcal{F}(\tau^{2})$, we have

[TABLE]

Therefore, by Lemma 13, we get

[TABLE]

if we choose $J^{\prime}=\left\lceil\log\left(C_{5}/\varepsilon\right)/2\right\rceil$ so that $C_{5}\,2^{-2J^{\prime}}\leq\varepsilon$.

To control the second term, we get from Lemma 24 that

[TABLE]

Moreover, using respectively the fact that $\operatorname{ext}$ is a Lipschitz operator on Besov spaces, $\Omega_{Q}$ is bounded, and the convex conjugate is non-expansive in $L^{\infty}$ by Lemma 20, we have

[TABLE]

where in the last inequality, we used Lemma 24. We have proved that

[TABLE]

for some constant $C_{6}=C_{6}(\tilde{\Omega}_{P},\Omega_{Q})$. This inequality allows us to control the $L^{\infty}$-bracketing numbers of $\tilde{\mathcal{F}}_{J}(\tau^{2})$ using $\ell^{\infty}$-covering numbers for the wavelet coefficients. To control the latter, note that for all $f\in\mathcal{F}_{J}(\tau^{2})\subset\mathcal{X}(2M)$, it holds $\|f\|_{B^{2}_{\infty,\infty}(\mathbbm{R}^{d})}\lesssim M^{2}$ so that $\|\gamma\|_{\ell^{\infty}}\lesssim M^{2}$. Moreover, these wavelet coefficients are in a space of dimension at most $C\,2^{Jd}$, $C(\tilde{\Omega}_{P},M)>0$ because $\mathcal{F}_{J}(\tau^{2})\subseteq V_{J}(\tilde{\Omega}_{P})$. Hence, choosing $\varepsilon=\delta/4$, Lemmas 27, 28, and the previous display yield

[TABLE]

To apply the chaining bound of Theorem 25, note that by (F.10) and Lemma 37, respectively, combined with Lemma 13, we have

[TABLE]

Thus,

[TABLE]

where $\mathcal{D}_{[\,]}=\mathcal{D}_{[\,]}(C_{4}\,2^{-J}\tau,\tilde{\mathcal{F}}_{J}(\tau^{2}),L^{2}(Q))$ is the Dudley integral defined in (C.1). Moreover, by (F.27) and Lemma 29, we have

[TABLE]

and therefore,

[TABLE]

Together with (F.18), we have $T_{2}\lesssim\phi_{J}(\tau^{2})$ with $\phi_{J}$ defined as in (5.21) and we absorbed $\mathfrak{r}_{J}$ for $d\geq 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}\lesssim\phi_{J}(\tau^{2})$.

Final bound and concentration. Collecting the above bounds on $T_{1}$ and $T_{2}$, we get

[TABLE]

To obtain a bound that holds with high probability, we apply the concentration result of Theorem 26. For this, note that (F.2) can be written as

[TABLE]

where $P\otimes Q$ denotes the product measure and $(f\otimes g)(x,y)=f(x)+g(y)$. Following the same argument as in the proof of Lemma 37, we get

[TABLE]

Moreover, similarly to Proposition 10, we have for any $f\in\mathcal{F}_{J}(\tau^{2})$ that

[TABLE]

We can therefore apply Theorem 26 and conclude that with probability at least $1-e^{-t}$, it holds that

[TABLE]

which concludes our proof.

F.3 Proof of Lemma 37

By the boundedness conditions in the definition of $\mathcal{X}(2M)$ in Proposition-Definition 7, and the boundedness of $\tilde{\Omega}_{P}$ and $\Omega_{Q}$, we have for $y\in\Omega_{Q}$ that

[TABLE]

Moreover, by Lemma 17, $\nabla f^{\ast}$ is $2M$-Lipschitz. Therefore, since $\Omega_{Q}$ is bounded,

[TABLE]

We can similarly deduce the second claim by Proposition-Definition 7 since $f_{0}\in\mathcal{X}(M)$, possibly choosing a larger $C_{5}$.

Appendix G Additional lemmas

Lemma 38**.**

In the notation of the proof of Theorem 6, for any $k=0,\ldots,K$ and $m\geq m_{0}$, $\nabla\phi_{k}$ is a bijection from $[0,1]^{d}$ to $[0,1]^{d}$.

Proof of Lemma 38.

By construction, $\phi_{k}$ is strongly convex and has Lipschitz continuous derivatives, hence so does its convex conjugate $\phi_{k}^{\ast}$ by Lemma 17. In particular, $\phi_{k}^{\ast}$ is defined on all of $\mathbbm{R}^{d}$ and for each $y$,

[TABLE]

Hence, the equation $\nabla\phi_{k}(x)=y$ has a unique solution $x(y)$ for every $y\in\mathbbm{R}^{d}$ which implies that $\nabla\phi_{k}$ is injective and that for any $y\in[0,1]^{d}$, there exists $x=x(y)\in\mathbbm{R}^{d}$ such that $\nabla\phi_{k}(x)=y$. It remains to check that $x(y)\in[0,1]^{d}$ for all $y\in[0,1]^{d}$. To that end, note that $\phi_{k}(x)=\|x\|^{2}/2$ for $x\notin[0,1]^{d}$. Hence $x(y)=y$ whenever $y\notin[0,1]^{d}$. In particular, if $y\in[0,1]^{d}$, we must have $x(y)\in[0,1]^{d}$. This completes the proof. ∎

Lemma 39** (Poincaré inequality, [Eva10, Section 5.8.1], [Leo17, Theorem 13.27]).**

Let $\Omega\subseteq\mathbbm{R}^{d}$ be a bounded and connected Lipschitz domain. Then, there exists a constant $C=C(d,\Omega)$ such that for any function $f\in W^{1,2}(\Omega)$,

[TABLE]

Appendix H Numerical Experiments, continued

In this section, we give additional details on the numerical experiments in Section 6, using the same notation used there to define $\hat{T}_{\mathrm{emp}}$, $\hat{T}_{\mathrm{wav}}$ and $\hat{T}_{\mathrm{ker}}$.

H.1 Implementation details

All simulations are done with Python 3.8.0 and Numpy 1.17.3, where some calculations are accelerated by the just-in-time compiler of Numba 0.47.0, in particular the calculation of the Linear-Time Legendre transform [Luc97]. The discrete wavelet transform is calculated with PyWT 1.1.1 [LGW*+*19]. We use a second-order finite difference operator for the calculation of the numerical gradient provided by the findiff package, version 0.8.0. The optimization (6.7) is performed with the L-BFGS algorithm [LN89] as implement in Scipy 1.4.1, stopping at a relative decrease in objective function value of less than $10^{-9}$ and a maximum iteration number of $10000$. The baseline estimator is computed with the ot.emd function of the Python Optimal Transport package, version 0.6.0. The kernel regression problem (6.13) is solved via scikit-learn, version 0.22.2 [PVG*+*11]. Plots were made using Matplotlib 3.1.2 and Seaborn 0.9.0.

The boxes in the calculation of $\hat{\bm{f}}_{J}$ are picked to be $\tilde{\Omega}_{P}=\tilde{\Omega}_{Q}=[-0.5,1.5]^{d}$ in the case (id) and $\tilde{\Omega}_{P}=[-0.5,1.5]^{d}$, $\tilde{\Omega}_{Q}=[0,4]^{d}$ in the case (exp). To compute $\hat{\bm{f}}_{J}$ for different $J$, we initialize the optimization at all scales of $J$ with the all zeros vector. An approximation to the ground truth semi-dual objective (2.14) is obtained by integrating both $W^{\top}_{J}\gamma_{J}$ and the interpolation $\mathcal{P}_{\bm{y}\to T_{0}(\bm{x})}\mathcal{L}(W^{\top}_{J}\gamma_{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 $\nu_{\mathrm{kernel}}$ and $\nu_{\mathrm{ridge}}$ are chosen via an oracle procedure, picking the best $\nu_{\mathrm{kernel}}\in\{10^{-9},10^{-8.5},\dots,10^{-5}\}$, $\nu_{\mathrm{ridge}}\in\{10^{-5},10^{-4.5},\dots,10^{-1}\}$ as determined by evaluating $\mathrm{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® Xeon® E5-2670 v3 (2.30GHz), the calculation of all wavelet scales of $\hat{T}_{\mathrm{wav}}$ 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$ different parameter settings for $\hat{T}_{\mathrm{ker}}$ with $n=10000$ takes about five hours.

H.2 Numerical error dominates for large sample size

To investigate the observed flattening out of the $\hat{T}_{\mathrm{wav}}$ 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.

H.3 Improved gradient estimation by computing convex envelope

Write $\hat{\bm{f}}_{J}=W^{\top}_{J}\hat{\gamma}_{J}$. We observe that while $\hat{\bm{f}}_{J}$ by itself often qualitatively yields good results on $\Omega_{P}$ (the support of $P$), the lack of additional global regularity conditions in the optimization (6.7) often leads to poor approximation outside of the support, i.e., on $\tilde{\Omega}_{P}\setminus\Omega_{P}$. In particular, certain regions could simply not enter the optimization at all due to the form of the discrete Legendre transform $\mathcal{L}$ in (6.5). 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 $\hat{\bm{f}}_{J}$.

We give a qualitative visual example of this in Figure 5. Depicted for $d=2$ are first the ground truth potential $f_{0}(x)=\frac{1}{2}\|x\|_{2}^{2}$ corresponding to the identity transport map (Figure 5(a)) and its derivative in the first coordinate (Figure 5(b)), where we set $P=\mathsf{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=\mathsf{Unif}([0,1]^{2})$ are shown in Figures 5(c) (potential) and 5(d) (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 $\Omega_{P}$, the outside appears very ragged. This is remedied by applying $\mathcal{L}$ twice, see Figures 5(e) (potential) and 5(f) (numerical gradient).

We also quantitatively compare the estimation accuracy of a linear interpolation of $\nabla_{\bm{x}}\hat{\bm{f}}_{\hat{J}}$, denoted by $\hat{T}_{\mathrm{wav,\,direct}}$ to that of $\hat{T}_{\mathrm{wav}}$ 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.

H.4 Performance for varying regularization strength

In order to further investigate the performance of the proposed regularized estimators $\hat{T}_{\mathrm{wav}}$ and $\hat{T}_{\mathrm{ker}}$, we show error plots for a selection of fixed regularization parameters in Figure 6. In the case of $\hat{T}_{\mathrm{wav}}$, we plot $\hat{T}_{\mathrm{wav}}^{(J)}$ for the whole range of $J\in\{0,\dots,3\}$ that can be calculated by means of the discrete wavelet transform when $N=65$. In the case of $\hat{T}_{\mathrm{ker}}$, we plot a subset of all considered combinations $\nu_{\mathrm{kernel}}\in\{10^{-9},10^{-8.5},\dots,10^{-5}\}$, $\nu_{\mathrm{ridge}}\in\{10^{-5},10^{-4.5},\dots,10^{-1}\}$.

First, for the wavelet estimators $\hat{T}_{\mathrm{wav}}^{(J)}$, we observe that the best bias-variance trade-off is achieved by $J=0$ for smaller values of $n$ and by $J=1$ for $n\geq 10^{3.5}$. Moreover, the curve for $J=3$, corresponding to no regularization, roughly follows the shape of the error curve of $\hat{T}_{\mathrm{emp}}$, until it flattens out due to the numerical error. This is most likely due to the grid approximations involved in the definition of $\hat{T}_{\mathrm{wav}}$.

Next, the kernel estimators $\hat{T}_{\mathrm{ker}}^{(\nu_{\mathrm{ridge}},\nu_{\mathrm{kernel}})}$ in 3D show mostly similar error curves for the range of parameters considered, although a trade-off that depends on the sample size can be observed between $(\nu_{\mathrm{ridge}},\nu_{\mathrm{kernel}})\in\{(10^{-7},10^{-4}),(10^{-7},10^{-3})\}$. Moreover, in the case of very strong regularization, $(\nu_{\mathrm{ridge}},\nu_{\mathrm{kernel}})=(10^{-5},10^{-5})$, the error curve flattens out for values of $n$ as low as $10^{2.5}$. In general, for $\hat{T}_{\mathrm{ker}}$, we observe that its performance is more sensitive to changes in $\nu_{\mathrm{kernel}}$ than to those in $\nu_{\mathrm{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 $\hat{T}_{\mathrm{emp}}$ for $(\nu_{\mathrm{ridge}},\nu_{\mathrm{kernel}})=(10^{-9},10^{-1})$ to a curve that matches most of the error curve shown in Figure 2(d) for $(\nu_{\mathrm{ridge}},\nu_{\mathrm{kernel}})=(10^{-7},10^{-4})$.