Minimax estimation of smooth densities in Wasserstein distance

TL;DR

This paper establishes the first minimax-optimal rates for nonparametric density estimation under Wasserstein distance, revealing dependence on density bounds and providing computationally feasible measures.

Contribution

It introduces the first minimax rates for Wasserstein-based density estimation and constructs discretely supported measures achieving these rates.

Findings

Derived minimax-optimal rates for general Wasserstein distances.

Showed the dependence of rates on whether densities are bounded below.

Provided discretely supported measures suitable for computation.

Abstract

We study nonparametric density estimation problems where error is measured in the Wasserstein distance, a metric on probability distributions popular in many areas of statistics and machine learning. We give the first minimax-optimal rates for this problem for general Wasserstein distances, and show that, unlike classical nonparametric density estimation, these rates depend on whether the densities in question are bounded below. Motivated by variational problems involving the Wasserstein distance, we also show how to construct discretely supported measures, suitable for computational purposes, which achieve the minimax rates. Our main technical tool is an inequality giving a nearly tight dual characterization of the Wasserstein distances in terms of Besov norms.

Equations

$$W_p^p(\mu ,\nu )=\underset{\pi \in \mathcal{M}(\mu ,\nu )}{\inf }\int \Vert x-y{\Vert }^p\mathrm{d}\pi (x,y),$$

$$W_1(\mu ,\nu )=\underset{f\in \mathrm{Lip}}{\sup }\int f(\mathrm{d}\mu -\mathrm{d}\nu ),$$

$${\Delta }_h^{1}f(x)$$

$${\Delta }_h^{k}f(x)$$

$${\omega }_n(f,t{)}_p=\underset{\Vert h\Vert \le t}{\sup }\Vert {\Delta }_h^{n}f{\Vert }_{L_p({\Omega }_{h,n})}.$$

$$\Vert f{\Vert }_{{\mathcal{B}}_{p,q}^s}^{\prime }:=\Vert f{\Vert }_{L_p}+{\Vert (2^{sj}{\omega }_n(f,2^{-j}{)}_p{)}_{j\ge 0}\Vert }_{{\ell }_q}$$

$$\Vert f{\Vert }_{{\mathcal{B}}_{p,q}^s}:=\Vert \alpha {\Vert }_{{\ell }_p}+{\Vert 2^{js}{2}^{dj(\frac{1}{2}-\frac{1}{p})}\Vert {\beta }_j{\Vert }_{{\ell }_p}\Vert }_{{\ell }_q}.$$

$${\mathcal{B}}_{p,q}^s(L)$$

$${\mathcal{B}}_{p,q}^s(L;m)$$

$$\sup_{f\in\mathcal{B}^{s}_{p^{\prime},q}(L;m)}{\rm I}\kern-1.79993pt{\rm E}W_{p}(f,\hat{f})\lesssim\left\{\begin{array}[]{ll}n^{-\frac{1+s}{d+2s}}&d\geq 3\\ n^{-1/2}\log n&d=2\\ n^{-1/2}&d=1\,.\end{array}\right.$$

$$\sup_{f\in\mathcal{B}^{s}_{p^{\prime},q}(L;m)}{\rm I}\kern-1.79993pt{\rm E}W_{p}(f,\hat{f}^{\circ})\lesssim\left\{\begin{array}[]{ll}n^{-\frac{1+s}{d+2s}}\log n&d\geq 3\\ n^{-1/2}(\log n)^{2}&d=2\\ n^{-1/2}&d=1\,.\end{array}\right.$$

$$\adjustlimits{\inf}_{\tilde{\mu}}{\sup}_{f\in\mathcal{B}^{s}_{p^{\prime},q}(L;m)}{\rm I}\kern-1.79993pt{\rm E}W_{p}(f,\tilde{\mu})\gtrsim\left\{\begin{array}[]{ll}n^{-\frac{1+s}{d+2s}}&d\geq 2\\ n^{-1/2}&d=1\,,\end{array}\right.$$

$$M^{-1/p^{\prime }}\Vert f-g{\Vert }_{{\mathcal{B}}_{p,\infty }^{-1}}\lesssim W_p(f,g)\lesssim m^{-1/p^{\prime }}\Vert f-g{\Vert }_{{\mathcal{B}}_{p,1}^{-1}},$$

$$\adjustlimits{\inf}_{\tilde{\mu}}{\sup}_{f\in\mathcal{B}^{s}_{p^{\prime},q}(L)}W_{p}(f,\tilde{\mu})\gtrsim\left\{\begin{array}[]{ll}n^{-\frac{1+s/p}{d+s}}&d-s\geq 2p\\ n^{-1/2p}&d-s<2p\end{array}\right.$$

$$\sup_{f\in\mathcal{C}^{s}(L)}{\rm I}\kern-1.79993pt{\rm E}W_{p}(f,\hat{f})\lesssim\left\{\begin{array}[]{ll}n^{-\frac{1+s/p}{d+s}}&d-s>2p\\ n^{-\frac{1}{2p}}\log n&d-s=2p\\ n^{-\frac{1}{2p}}&d-s<2p\,.\end{array}\right.$$

$$f$$

$$g$$

$$W_{p}(f,g)\lesssim m^{-1/p^{\prime}}\Big{(}\|\alpha-\alpha^{\prime}\|_{\ell_{p}}+\sum_{j\geq 0}2^{-j}2^{dj(\frac{1}{2}-\frac{1}{p})}\|\beta_{j}-\beta_{j}^{\prime}\|_{\ell_{p}}\Big{)}\,,$$

$$\nabla \cdot V$$

$$\Vert V{\Vert }_{L^p}$$

$$-{\int }_{\Omega }\nabla h\cdot V\mathrm{d}x={\int }_{\Omega }h(f-g)\mathrm{d}x.$$

$$W_p^p(\mu ,\nu )=\underset{(\rho ,E)\in K_{\Omega }}{\inf }\left\{{\mathcal{B}}_p(\rho ,E):\rho (\cdot ,0)=\mu ,\rho (\cdot ,1)=\nu ,{\partial }_t\rho +{\nabla }_x\cdot E=0\right\},$$

$$\mathcal{B}_{p}(\rho,E):=\left\{\begin{array}[]{ll}\int_{\Omega\times[0,1]}\left\|\frac{dE}{d\rho}(x,t)\right\|^{p}\,\mathrm{d}\rho(x,t)&\text{ if $E\ll\rho$,}\\ +\infty&\text{ otherwise.}\end{array}\right.$$

$$\lambda(x):=\left\{\begin{array}[]{ll}\frac{1}{2}(2t)^{p}&\text{ if $t\leq 1/2$,}\\ 1-\frac{1}{2}(2-2t)^{p}&\text{ if $t>1/2$.}\end{array}\right.$$

$$\rho(x,t)\geq\left\{\begin{array}[]{ll}\frac{1}{2}(2t)^{p}m&\text{ if $t\leq 1/2$,}\\ \frac{1}{2}(2-2t)^{p}m&\text{ if $t>1/2$}\end{array}\right.$$

$$E(x,t)={\lambda }^{\prime }(t)V(x)\text{ for t∈[0,1]∖{1/2}.}$$

$${\Vert \frac{dE}{d\rho }(x,t)\Vert }^p\rho (x,t)\le \Vert V(x){\Vert }^p\frac{{\lambda }^{\prime }(t{)}^p}{\rho (x,t{)}^{p-1}}\lesssim \Vert V(x){\Vert }^p{m}^{1-p}.$$

$$W_p(f,g)$$

$$\lesssim m^{1/p-1}\Vert V{\Vert }_{L_p}$$

$$\lesssim m^{-1/p^{\prime }}\left(\Vert \alpha -{\alpha }^{\prime }{\Vert }_{{\ell }_p}+\sum _{j\ge j_0}{2}^{-j}{2}^{dj(\frac{1}{2}-\frac{1}{p})}\Vert {\beta }_j-{\beta }_j^{\prime }{\Vert }_{{\ell }_p}\right),$$

$$V=\sum _{\varphi \in \Phi }({\alpha }_{\varphi }-{\alpha }_{\varphi }^{\prime })V_{\varphi }+\sum _{j\ge 0}\sum _{\psi \in {\Psi }_j}({\beta }_{\psi }-{\beta }_{\psi }^{\prime })V_{\psi }.$$

$${\varphi }^{(i)}(x^{(i)})={\int }_0^1\cdots {\int }_0^1\varphi (x^{(i)},t_{i+1},\dots ,t_d)\mathrm{d}{t}_{i+1}\dots \mathrm{d}{t}_d,$$

$$(V_{\varphi }{)}_i(x)={\int }_0^{x_i}{\varphi }^{(i)}(x^{(i-1)},t_i)\mathrm{d}{t}_i-x_i{\varphi }^{(i-1)}(x^{(i-1)})1\le i\le d.$$

Full text

62F99,

62H99,

and

t1Parts of this research were conducted while JNW was at the Institute for Advanced Study and at Massachusetts Institute of Technology, supported by the Josephine de Karman Fellowship. t2Most of this work was conducted while QB was at University of Cambridge and was supported in part by The Alan Turing Institute under the EPSRC grant EP/N510129/1.

1 Introduction

Optimal transport is a fundamental problem in geometry, optimization, and analysis (Monge, 1781; Kantorovitch, 1942; Villani, 2008), with increasing applications in statistics. A great deal of recent work in the machine learning community has shown that optimal transport can be used to develop successful empirical methodologies (see, e.g., Peyré and Cuturi, 2017, for a survey) at a relatively low computational cost (Cuturi, 2013; Altschuler et al., 2017). This raises the need for a rigorous theoretical understanding of the statistical properties of optimal transport.

The Wasserstein distances are a metric on probability distributions defined using optimal transport. A metric is defined by asking how one can transport mass with distribution $\mu$ to have another distribution $\nu$, with minimal global transport cost. This problem also has the probabilistic interpretation, known as the Monge–Kantorovich formulation, of finding a joint distribution $\pi$ minimizing a cost for variables $X$ and $Y$ with given marginals. The Wasserstein distance emerges as the minimum value of this problem, and creates a natural tool to compare distributions, with $W_{p}$ corresponding to the $\|\cdot\|^{p}$ transport cost:

[TABLE]

where the set $\mathcal{M}(\mu,\nu)$ denotes the set of joint measures with marginals $\mu$ and $\nu$, respectively. This definition naturally extends to positive measures $\mu$ and $\nu$ with the same total mass.

In many modern applications, a Wasserstein distance is used as a loss function in an optimization problem over measures. Solving such problems involves optimizing functionals of the form $\nu\mapsto W_{p}(\nu,\mu)$ where $\mu$ is unknown. Given $n$ i.i.d. samples from $\mu$, much of the statistics literature adopts the plug-in approach and focuses on using the empirical distribution $\hat{\mu}_{n}$ to obtain the estimated functional $\nu\mapsto W_{p}(\nu,\hat{\mu}_{n})$. In this case, the rates of convergence are of order $n^{-1/d}$, and the sample size required for a particular precision is exponential in the dimension, a phenomenon known as the curse of dimensionality. Moreover, it is known that this exponential dependence is tight, in the sense that no better estimate is available in general (Singh and Póczos, 2018).

Our work adopts a different approach to show that the plug-in estimator is suboptimal for measures possessing a smooth density. Estimating the density of a distribution, based on independent samples, is one of the fundamental problems of statistics. The usual goal in these problems is to produce an estimate $\tilde{f}$ which is as close as possible to the unknown density $f$, measured either at one point of the sample space, or in $L_{p}$ norm. In this line of work, $f$ is usually assumed to belong to a large, nonparametric class defined via smoothness or regularity conditions, and typically the rates obtained in this setting show that sufficient smoothness can substantially mitigate the curse of dimensionality. This is the subject of a wide literature on nonparametric density estimation (see, e.g Tsybakov, 2009). In this work, we follow the same philosophy and derive similar rates for $W_{p}$ distances, over Besov classes of densities $\mathcal{B}^{s}_{p,q}$. We likewise show that the smoothness parameter $s$ improves the optimal exponent of $n$ in the Wasserstein setting.

The key challenge in proving tight rates for nonparametric density estimation in the Wasserstein distances is that the distances $W_{p}$ are both non-linear and non-local. Unlike $L_{p}$ norms or Besov norms, commonly considered in the nonparametric estimation literature, the $W_{p}$ distance is not a norm on the space of probability measures when $p>1$. Moreover, it is not even “norm-like”: we show (theorem 9) that when $p>1$, the Wasserstein distance $W_{p}$ cannot be controlled by any reasonable norm on densities. Moreover, the optimization problem (1) that defines the Wasserstein distance is highly non-local: modifications to $\mu$ and $\nu$ in a small neighborhood can change the global structure of the optimal coupling $\pi$. This makes analyzing the behavior of an estimator challenging.

Algorithmic aspects are an important part of optimal transport problems. For practical applications, the proposed estimates must therefore also be computationally tractable. We describe a method inspired by the parametric bootstrap (Efron and Tibshirani, 1993) to produce computationally tractable atomic estimators from any estimator that outperforms the empirical distribution, under minimal assumptions. We study the computational cost of this method, compared to the cost of using the empirical distribution with $n$ atoms, and exhibit a trade-off between computational cost and statistical precision.

1.1 Prior work

The question of establishing minimax rates for estimation in Wasserstein distances has been examined in several recent works. Singh and Póczos (2018) established that, in the absence of smoothness assumptions, the empirical distribution $\hat{\mu}$ is rate optimal in a variety of examples. Their proof relies on a dyadic partitioning argument (see, e.g. Weed and Bach, 2018), and does not appear to extend to the smooth case. Closer to our setting, under a smoothness assumption on the density of $\mu$, Liang (2017) and Singh et al. (2018) showed minimax rates of convergence for the Wasserstein-1 distance. To obtain these rates, these works focus on the dual form of $W_{1}$:

[TABLE]

where the supremum is taken over all 1-Lipschitz functions. This dual formulation puts the Wasserstein-1 distance into the category of integral probability metrics (Müller, 1997), for which both Liang (2017) and Singh et al. (2018) obtain general results. It has been shown that choosing functions which are smoother than Lipschitz in this definition can result in improved rates of convergence for empirical measures (Kloeckner, 2018). Crucially, the metric $W_{p}$ for $p>1$ is not an integral probability metric. Establishing sharp rates for general Wasserstein distances therefore requires different techniques. A separate line of work has focused instead on modifying the definition of the Wasserstein distance to include a regularizing term based on the mutual information of the coupling. It has been shown that this definition enjoys improved convergence rates relative to the unregularized version (Genevay et al., 2018; Mena and Niles-Weed, 2019).

Our proofs rely on establishing control of Wassserstein distances by Besov norms of negative smoothness. Similar results have been obtained elsewhere under different conditions. Shirdhonkar and Jacobs (2008) showed that the optimal transportation distance with cost $\|\cdot\|^{p}$ for $0<p<1$ can be characterized explicitly via an expression involving wavelet coefficients, which implies that these distances agree with a particular Besov norm (see section 2.1.1); however, their proof technique does not extend to $p>1$. Loeper (2006) (see also Maury et al., 2010) showed that the Wasserstein-2 distance between measures with densities bounded above dominates a negative Sobolev norm, and Peyre (2018) extended this result to show that $W_{2}$ is in fact equivalent to such a norm when the densities are in addition bounded below. Similar connections have been developed by Leeb and Coifman (2015) and Ledoux (2017). To our knowledge, ours is the first result to establish a connection to Besov norms of negative smoothness and general Wasserstein distances.

The use of wavelet estimators for density estimation has a long history in nonparametric statistics (Kerkyacharian and Picard, 1992; Härdle et al., 1998; Donoho et al., 1996; Doukhan and León, 1990; Walter, 1992). However, while wavelets have been used for computational purposes in the optimal transport community (Chen et al., 2012; Shirdhonkar and Jacobs, 2008; Dominitz et al., 2008; Rabin et al., 2011), the statistical properties of wavelet estimators with respect to Wasserstein distances have remained largely unexplored.

2 Main results

2.1 Problem description and preliminaries

Our observation consists of an i.i.d. sample of size $n$ drawn from a probability measure on $\mathbb{R}^{d}$ with smooth density $f$. Our goal is to compute an estimator $\tilde{\mu}_{n}$ that is close to $\mu_{f}$ in expected Wasserstein distance. As noted above, such an estimator can serve as a proxy for $\mu_{f}$ in statistical and computational applications. While estimation of the density $f$ in norms such as $L_{p}$ is a well studied problem in nonparametric statistics (Tsybakov, 2009), such estimates do not readily lend themselves to guarantees in Wasserstein distance.

For technical reasons, we restrict ourselves to the case of compactly supported measures. We further assume that this support is contained in $\Omega:=[0,1]^{d}$. The general case can be reduced to this one by a rescaling argument, as long as the statistician can obtain a priori bounds on the size of the support of the measure. We do not address the question of obtaining such bounds in this work, but note that they are often available in practice.

2.1.1 Wavelets and Besov spaces

We direct the reader to Härdle et al. (1998) and Meyer (1990) for an introduction to the theory of wavelets. In brief, we assume the existence of sets $\Phi$ and $\Psi_{j}$ for $j\geq 0$ of functions in $L_{2}(\Omega)$ satisfying the standard requirements of a wavelet basis. (See Supplemental Appendix E for our precise assumptions.)

Wavelets can be used to characterize the Besov spaces $\mathcal{B}^{s}_{p,q}(\Omega)$. We follow the approach of Cohen (2003) for defining such spaces on bounded domains. Suppose $s>0$ and $p,q\geq 1$, and let $n>s$ be an integer. Given $h\in\mathbb{R}^{d}$, set

[TABLE]

where these functions are defined on $\Omega_{h,n}:=\{x\in\Omega:x+nh\in\Omega\}$. For $t>0$, we then define

[TABLE]

The function $\omega_{n}$ measures the order-$n$ smoothness of $f$ in $L_{p}$. Finally, we define the space $\mathcal{B}^{s}_{p,q}(\Omega)$ to be the set of functions for which the quantity

[TABLE]

is finite.

Assuming that the elements of $\Phi$ and $\Psi_{j}$ have $r$ continuous derivatives for $r>s$ and that polynomials of degree up to $\lfloor s\rfloor$ lie in the span of $\Phi$, the norm $\|\cdot\|^{\prime}_{\mathcal{B}^{s}_{p,q}}$ is equivalent to a sequence norm based on wavelet coefficients. Given $f\in L_{p}(\Omega)$, denote by $\alpha=\{\alpha_{\phi}\}_{\phi\in\Phi}$ the vector defined by $\alpha_{\phi}:=\int f\phi$ and by $\beta_{j}=\{\beta_{\psi}\}_{\psi\in\Psi_{j}}$ the vector whose entries are given by $\beta_{\psi}:=\int f\psi$. Then $\|\cdot\|^{\prime}_{\mathcal{B}^{s}_{p,q}}$ is equivalent to $\|\cdot\|_{\mathcal{B}^{s}_{p,q}}$ defined by

[TABLE]

This expression can then be used directly to define a norm when $s<0$ (see Cohen, 2003, Theorem 3.8.1), as long as the elements of $\Phi$ and $\Psi_{j}$ have $r$ continuous derivatives for $r>|s|$ and polynomials of degree up to $\lfloor|s|\rfloor$ lie in the span of $\Phi$. In what follows, we therefore adopt (2) as our primary definition and assume throughout that the wavelet system has sufficient regularity that the equivalence of $\|\cdot\|^{\prime}_{\mathcal{B}^{s}_{p,q}}$ and $\|\cdot\|_{\mathcal{B}^{s}_{p,q}}$ holds.

2.1.2 Notation

The quantities $C$ and $c$ will refer to constants whose value may change from line to line. All constants throughout may depend on the choice of wavelet system and the dimension. Since we are interested in establishing optimal rates of decay with respect to the exponent (i.e., finding $\gamma$ such that the rate $n^{-\gamma}$ holds), we leave finer control on dimension-dependent constants to future work. We freely use the notation $a\lesssim b$ to indicate that there exists a positive constant $C$ for which $a\leq Cb$ holds. Again, such constants may depend on the multiresolution and dimension. The notation $a\asymp b$ indicates that $a\lesssim b$ and $b\lesssim a$.

Given $a,b\in\mathbb{R}$, we denote by $a\wedge b$ and $a\vee b$ the minimum and maximum of $a$ and $b$, respectively. The symbol $\|\cdot\|$ denotes the Euclidean norm on $\mathbb{R}^{d}$.

We set denote by $\mathcal{D}(K)$ the set of probability density functions on $\Omega$. Since the measures we consider are absolutely continuous with respect to the Lebesgue measure, we do not distinguish between a probability measure and its density. In particular, we write $W_{p}(f,g)$ for the $p$-Wasserstein distance between the measures with densities $f$ and $g$.

2.2 Minimax estimation of smooth densities

In this section, we give our main statistical results on the problem of estimating densities in Wasserstein distance. These results reveal several phenomena: (i) the minimax rate of estimation can improve significantly for smooth densities, and (ii) the optimal rates depend strongly on whether the density in question is bounded away from [math]. Indeed, we show that the optimal rate for general densities is strictly worse than the corresponding rate for densities bounded below, no matter the smoothness. While the first phenomenon is well known in nonparametric statistics, the second phenomenon does not occur in classical density estimation problems. As we explore further below, this behavior is fundamental to the Wasserstein distances.

We define two classes of probability densities on $\Omega$. Given $m,L>0$, set

[TABLE]

We note that if $s$ is sufficiently large and $L$ is sufficiently small then in fact $\mathcal{B}^{s}_{p,q}(L)\subseteq\mathcal{B}^{s}_{p,q}(L;m)$ for $m$ a constant. We assume throughout that $m<1$, since when $m\geq 1$, the class $\mathcal{B}^{s}_{p,q}(L;m)$ is trivial.

2.2.1 Bounded densities

Our first result gives an upper bound on the rate of estimation for functions in $\mathcal{B}^{s}_{p,q}(L;m)$.

Theorem 1**.**

For any $p\geq 1$ and $s\geq 0$, there exists an estimator $\hat{f}$ such that for any $p\leq p^{\prime}<\infty$, $1\leq q\leq\infty$ and $m>0$, the estimator satisfies

[TABLE]

The upper bound in theorem 1 is achieved by a wavelet estimator. As $s$ ranges between [math] and $\infty$, the upper bound interpolates between the dimension-dependent rate $n^{-1/d}$ and the fully parametric rate $n^{-1/2}$.

The estimator proposed in theorem 1 is adaptive to the lower bound $m$, but not to the smoothness $s$. An adaptive version of this estimator can be obtained at the price of an extra logarithmic factor.

Theorem 2**.**

For any $p\geq 1$ and $L>0$, there exists an estimator $\hat{f}^{\circ}$ such that for any $p\leq p^{\prime}<\infty$, $1\leq q\leq\infty$, $m>0$, and $s\geq 0$, the estimator satisfies

[TABLE]

Though this estimator is adaptive to $s$, $p^{\prime}$, and $q$, using it requires an a priori estimate of $L$ (or, more specifically, the $L_{p}$ norm of $f$). A similar phenomenon is present when designing adaptive density estimators in classical nonparametric statistics (Donoho et al., 1996).

Our lower bounds nearly match the upper bounds proved in Theorem 1, up to a logarithmic factor in the $d=2$ case.

Theorem 3**.**

For any $p,p^{\prime},q\geq 1$, and $s\geq 0$,

[TABLE]

where the infimum is taken over all estimators $\tilde{\mu}$ based on $n$ observations.

Our core technical contribution is the following result, which establishes a connection between Wasserstein distances and Besov norms of negative smoothness. Theorem 4 allows us to bypass the difficulties of the Wasserstein distance by bounding instead a nearly equivalent norm.

Theorem 4**.**

Let $p\in[1,\infty)$. Let $f,g$ be two densities in $L_{p}([0,1]^{d})$, and assume $M\geq f(x)\vee g(x)\geq m>0$ for almost every $x\in[0,1]^{d}$. Then

[TABLE]

where $\frac{1}{p}+\frac{1}{p^{\prime}}=1$.

Theorem 4 can be viewed as a partial extension of the dual formulation of $W_{1}$ to $W_{p}$ for $p>1$. Indeed, the inclusions $\mathcal{B}^{1}_{\infty,1}\subseteq\mathrm{Lip}\subseteq\mathcal{B}^{1}_{\infty,\infty}$, where $\mathrm{Lip}$ is the space of bounded Lipschitz functions, imply $\|f-g\|_{\mathcal{B}^{-1}_{1,\infty}}\lesssim W_{1}(f,g)\lesssim\|f-g\|_{\mathcal{B}^{-1}_{1,1}}$. Theorem 4 establishes the analogous result when $p>1$, but only when the densities $f$ and $g$ are bounded. A proof of this theorem appears in section 3.

We prove theorems 1, 2, LABEL: and 3 in section 4.

2.2.2 Unbounded densities

Obtaining minimax rates when the densities are no longer bounded below is significantly more challenging since, as theorem 9 below makes clear, it is no longer possible to control $W_{p}$ by any function norm in the absence of a lower bound on the densities in question. Moreover, the statistical properties of the estimation problem also change markedly: surprisingly, the density estimation problem over the class $\mathcal{B}^{s}_{p,q}(L)$ is strictly harder than the corresponding problem over $\mathcal{B}^{s}_{p,q}(L;m)$, even under a smoothness assumption. We prove the following lower bound.

Theorem 5**.**

For any $p,p^{\prime},q\geq 1$, and $s\geq 0$, if $L$ is a sufficiently large constant, then

[TABLE]

where the infimum is taken over all estimators $\tilde{\mu}$ based on $n$ observations.

Note that, when $p\geq 2$, this rate is worse than the upper bound given in theorem 1 for all $s>0$ and $d\geq 1$. This establishes that the class of densities bounded from below is strictly easier to estimate than the class of all densities, for all nontrivial smoothness parameters.

When $s\in[0,1)$, and we consider the Holder class $\mathcal{C}^{s}(L)$, we can obtain a nearly matching upper bound, up to a logarithmic factor. Moreover, the estimator we construct is a histogram. This property enables the use of such an estimator in practical applications. We take up this point in section 6.

Theorem 6**.**

Assume $p\geq 2$. For any $s\in[0,1)$, there exists a histogram estimator $\hat{f}$ such that

[TABLE]

The proofs of both theorems 6 and 5 appear in section 5.

2.3 Computational aspects of smooth density estimation

In many computational applications, it is significantly simpler to work with discrete measures supported on a finite number of points, since in general there is no closed form expression for the Wasserstein distance between continuous measures. Unfortunately, the estimators presented in section 2.2 are not of this form, so it is unclear whether smoothness of the underlying measure can be exploited in applications. However, a simple argument shows that optimal rates can be achieved by resampling from the smooth estimator we construct to obtain a discrete distribution supported on $M\geq n$ points which achieves the optimal rate for $s\in[0,1)$. We extract one simple result in this direction.

Theorem 7**.**

For any $s\in[0,1)$ and any $1\leq p<\infty$, there exists an estimator $\bar{\mu}_{n,M}$, supported on $M=o(n^{2})$ points, enjoying the same rate as in theorem 6, up to logarithmic factors. Moreover, $\bar{\mu}_{n,M}$ can be computed in time $O(M)$.

Additional computational considerations along with a proof of theorem 7 appear in section 6.

3 Controlling the Wasserstein distance by Besov norms

The main goal of this section is a proof of theorem 4, which establishes that the Wasserstein distance between two measures on $\Omega=[0,1]^{d}$ can be controlled by a Besov norm of the difference in their densities as long as their densities are bounded above and below. We also establish that no analogous result can hold for arbitrary densities. While we give upper and lower bounds, the Besov norms appearing in the two bounds do not agree. We do not know whether under some conditions the $W_{p}$ distance is in fact equivalent to a particular Besov norm $\|\cdot\|_{\mathcal{B}^{-1}_{p,q}}$ for some $q\in[1,+\infty]$.

The results of this section are closely related to results of Shirdhonkar and Jacobs (2008) and Peyre (2018), who established similar relations for $p<1$ and $p=2$, respectively. In section 3.1, we show the upper bound of theorem 4, and in section 3.2 we show the lower bound. Section 3.3 establishes that there is no general relationship between Wasserstein distances and Besov norms once the assumption that the density is bounded away from [math] is relaxed.

3.1 Upper bound

Let $f$ and $g$ be probability densities in $L_{p}(\Omega)$ for $p\in[1,\infty)$ with the following wavelet expansions.

[TABLE]

where we assume (Assumption E.2) that constant functions lie in the span of $\Phi$. For the upper bound, we do not need to assume any additional regularity—in particular, 1 holds for the Haar wavelet basis (see Triebel, 2010). In principle, the expansions in (3) hold only in $L_{2}$, but in fact convergence also holds in $L_{p}$ assuming that $f,g\in L_{p}(\Omega)$ (Härdle et al., 1998, Remark 8.4).

Proposition 1**.**

Let $1\leq p<\infty$. If $f(x)\vee g(x)\geq m>0$ for almost every $x\in[0,1]^{d}$, then

[TABLE]

where $\frac{1}{p}+\frac{1}{p^{\prime}}=1$.

By rescaling $f$ and $g$ by a positive real number, this proposition implies that the analogous claim holds for any two nonnegative functions $f$ and $g$ satisfying $\int_{\Omega}f=\int_{\Omega}g<\infty$.

Proof.

We will follow a strategy originally developed by Moser (1965) for the purpose of showing that all volume forms on a smooth, compact manifold are equivalent up to automorphism. We define a vector field $V$ on $\Omega$ satisfying

[TABLE]

where the first condition is intended in the distributional sense that

[TABLE]

for all $h\in C^{1}(\Omega)$. In particular, we require the boundary condition $V\cdot\mathbf{n}=0$ on $\partial\Omega$, where $\mathbf{n}$ is an outward-pointing normal. We defer the construction of this vector field to Proposition 2, below.

To show the theorem, we appeal to the following characterization of the Wasserstein distance. Denote by $K_{\Omega}$ the set of pairs of measures $(\rho,E)$ on $\Omega\times[0,1]$ where $\rho$ is scalar valued and $E$ is vector valued.

Theorem 8** (Benamou and Brenier, 2000; Brenier, 2003).**

For any measures $\mu$ and $\nu$ on $\Omega$ and $p\in[1,\infty)$,

[TABLE]

where

[TABLE]

Let us show how to prove the theorem. We choose $\rho$ and $E$ absolutely continuous with respect to the Lebesgue measure on $\Omega\times[0,1]$, and consequently identify them with their density. First, set $\rho(x,t)=(1-\lambda(t))f(x)+\lambda(t)g(x)$, where $\lambda:[0,1]\to[0,1]$ is defined by

[TABLE]

Clearly $\rho(\cdot,0)=f(\cdot)$ and $\rho(\cdot,1)=g(\cdot)$. Moreover, we have the following lower bound.

Lemma 1**.**

Suppose $f(x)\vee g(x)\geq m$ for almost all $x\in\Omega$. Then

[TABLE]

for almost every $x\in\Omega$. In particular, for $t\neq 1/2$, we have that $\rho(x,t)\geq\frac{1}{2}m\cdot\left(\frac{\lambda^{\prime}(t)}{p}\right)^{\frac{p}{p-1}}$.

We then define $E$ by

[TABLE]

Since $\nabla_{x}\cdot E=\lambda^{\prime}(t)(f-g)$ for almost all $x\in\Omega$ and $t\in[0,1]$, the pair $(\rho,E)$ defined in this way satisfies $\partial_{t}\rho+\nabla_{x}\cdot E=0$ in the distributional sense.

For almost every $t\in[0,1]$, we have the bound

[TABLE]

We obtain

[TABLE]

as claimed. ∎

All that remains is to establish the existence of the promised vector field $V$.

Proposition 2**.**

There exists a vector field $V$ satisfying the requirements of (4) and (5).

Proof.

We will proceed by defining vector fields $V_{\phi}$ for each $\phi\in\Phi$ and $V_{\psi}$ for each $\psi\in\Psi_{j}$, $j\geq 0$ satisfying $\nabla\cdot V_{\phi}=\phi$ and $\nabla\cdot V_{\psi}=\psi$, along with appropriate boundary conditions. The desired vector field $V$ will then be obtained as

[TABLE]

An application of Fubini’s theorem immediately yields that this definition satisfies (4) in the distributional sense. We conclude by obtaining the desired estimate for $\|V\|_{L_{p}}$ to show (5)

Definition of $V_{\phi}$ for $\phi\in\Phi$.

Given $x\in\mathbb{R}^{d}$, we write $x^{(i)}$ for the vector consisting of the first $i$ coordinates of $x$. For each $1\leq i\leq d$, define $\phi^{(i)}:\mathbb{R}^{i}\to\mathbb{R}$ by

[TABLE]

and set $\phi^{(0)}=0$. We define $V_{\phi}$ componentwise as

[TABLE]

We now verify that this definition satisfies the desired identity. The proof appears in Supplemental Appendix A.

Lemma 2**.**

The field $V_{\phi}$ satisfies $\nabla\cdot V_{\phi}=\phi$. Moreover, $\left(\sum_{\phi\in\Phi}(\alpha_{\phi}-\alpha^{\prime}_{\phi})V_{\phi}\right)\cdot\mathbf{n}=0$ on the boundary of $[0,1]^{d}$, where $\mathbf{n}$ is an outward-pointing normal.

Definition of $V_{\psi}$ for $\psi\in\Psi_{j}$, $j\geq 0$.

We adopt essentially the same construction as above. First, by Assumption E.3, $\psi$ can be written as $\bigotimes_{i=1}^{d}\psi_{i}$, where each $\psi_{i}$ is a univariate function. Assumptions E.1 and E.2 imply that $\int_{[0,1]^{d}}\psi(x)\,\mathrm{d}x=0$, so there exists an index $k\in[d]$ such that $\int_{[0,1]}\psi_{k}(x_{k})\,\mathrm{d}x_{k}=0$. We set

[TABLE]

and $(V_{\psi})_{i}=0$ for $i\neq k$.

Lemma 3**.**

The field $V_{\psi}$ satisfies $\nabla\cdot V_{\psi}=\psi$ and $V_{\psi}\cdot\mathbf{n}=0$ on the boundary of $[0,1]^{d}$, where $\mathbf{n}$ is an outward-pointing normal.

A proof appears in Supplemental Appendix A.

Norm estimates.

We now obtain an estimate for $\|V\|_{L_{p}}$. We require two lemmas, the proofs of which appear in Supplemental Appendix A.

Lemma 4**.**

For any sequence $\{\alpha_{\phi}\}_{\phi\in\Phi}$,

[TABLE]

Lemma 5**.**

For any sequence $\{\beta_{\psi}\}_{\psi\in\Psi_{j}}$,

[TABLE]

We obtain, for $V$ defined as in (6),

[TABLE]

which gives (5). ∎

3.2 Lower bound

We can prove a similar lower bound when $f$ and $g$ are bounded above. Unlike the assumption that the densities are bounded below required for 1, this assumption is relatively benign, insofar as it holds automatically for continuous densities on $[0,1]^{d}$. For 3, we require the wavelets in (3) to possess at least one continuous derivative (see Assumption E.2).

Proposition 3**.**

Let $1\leq p<\infty$. If $f(x)\vee g(x)\leq M$ for almost every $x\in[0,1]^{d}$, then

[TABLE]

Proof.

We use the following fact, due to Maury et al. (2010, Lemma 3.4 and Remark):

Lemma 6** (Maury et al., 2010).**

For all $h\in C^{1}(\Omega)$,

[TABLE]

Fix an index $j\geq 0$. Let $h$ be a function of the form

[TABLE]

for some vectors $\kappa$ and $\lambda$ satisfying $\|\kappa\|_{\ell_{p^{\prime}}}\leq 1$ and $\|\lambda\|_{\ell_{p^{\prime}}}\leq 2^{-j+dj(\frac{1}{2}-\frac{1}{p})}$.

We require the following bound, whose proof appears in Supplemental Appendix A.

Lemma 7**.**

If $\|\kappa\|_{\ell_{p^{\prime}}}\leq 1$ and $\|\lambda\|_{\ell_{p^{\prime}}}\leq 2^{-j+dj(\frac{1}{2}-\frac{1}{p})}$, then $\|\nabla h\|_{L_{p^{\prime}}(\Omega)}\lesssim 1$.

Applying Lemmas 6 and 7, we obtain

[TABLE]

Taking the supremum over $\kappa$ and $\lambda$ subject to the constraints $\|\kappa\|_{\ell_{p^{\prime}}}\leq 1$ and $\|\lambda\|_{\ell_{p^{\prime}}}\leq 2^{-j+dj(\frac{1}{2}-\frac{1}{p})}$ implies

[TABLE]

and taking the supremum over $j\geq 0$ yields the claim. ∎

3.3 Densities not bounded below

We now show that no statement like 1 can hold if one of the densities is not bounded away from zero. Indeed, in this case, under mild assumptions, it is impossible to control $W_{p}(f,g)$ by any function norm when $p>1$. This stands in sharp contrast to the fact that, when $p=1$, the dual formulation of $W_{1}$ implies that the Wasserstein distance is such a norm.

Theorem 9**.**

Let ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ be any norm on functions on $\Omega$, and suppose that there exists a function $h$ in $L_{1}(\Omega)$, not identically zero, satisfying

• •

$\int_{\Omega}h=0$**

• •

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|h\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}<\infty$**

• •

The sets $\overline{\{h>0\}}$ and $\overline{\{h<0\}}$ are disjoint.

Then for any $p>1$,

[TABLE]

Proof.

Let $h_{+}:=h\vee 0$ and $h_{-}:=-(h\wedge 0)$, so that $h=h_{+}-h_{-}$, and note that by assumption $\int_{\Omega}h_{+}=\int_{\Omega}h_{-}=:\rho>0$. For any $\lambda\in[0,1]$, set

[TABLE]

Note that $f_{\lambda},g_{\lambda}\in\mathcal{D}(\Omega)$, and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|f-g\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\frac{\lambda}{\rho}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|h\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. On the other hand, since the compact sets $\overline{\{h>0\}}$ and $\overline{\{h<0\}}$ are disjoint, there exist two sets $S$ and $T$ and $c>0$ such that $\operatorname*{supp}(h_{+})\subseteq S$ and $\operatorname*{supp}(h_{-})\subseteq T$ and $\|x-y\|\geq c$ for any $x\in S,y\in T$. Lemma 9, below, therefore implies $W_{p}(f,g)\geq c|\int_{S}f-g|^{1/p}=c\lambda^{1/p}$.

We obtain

[TABLE]

∎

4 Wavelet estimation in for bounded densities

In this section, we employ the results of section 3 to prove theorems 1, 2, LABEL: and 3. We show that the minimax rate over $\mathcal{B}^{s}_{p,q}(L;m)$ can be achieved by a wavelet estimator.

4.1 Upper bound

To prove Theorem 1, we introduce the following estimator based on a wavelet expansion of regularity $r>\max\{s,1\}$ (see Assumption E.2) truncated to level $J$, for some $J\geq 0$ to be chosen. Set

[TABLE]

and let $\tilde{f}:=\sum_{\phi\in\Phi}\tilde{\alpha}_{\phi}\phi+\sum_{0\leq j\leq J}\sum_{\psi\in\Psi_{j}}\tilde{\beta}_{\psi}\psi$. While such an estimator can already yield optimal rates in $L_{p}$ (Kerkyacharian and Picard, 1992), $\tilde{f}$ may fail to be a probability density, in which case the quantity $W_{p}(f,\tilde{f})$ is undefined. We therefore focus on the estimator

[TABLE]

where $\mathcal{D}$ is the set of probability densities on $\Omega$. By construction, $\hat{f}$ is a density, so that $W_{p}(f,\hat{f})$ is meaningful.

The following lemma shows that the quality of this estimator can be controlled by the distance between $f$ and $\tilde{f}$ in Besov norm.

Lemma 8**.**

For any $f\in\mathcal{B}^{s}_{p,q}(L;m)$,

[TABLE]

Proof.

By assumption, $f$ is bounded below by $m$. Theorem 4 therefore implies

[TABLE]

where the final inequality uses the definition of $\hat{f}$ and the fact that $f\in\mathcal{D}$. ∎

The proof of Theorem 1 now follows from standard facts in wavelet density estimation. We require the following proposition, whose proof appears in Supplemental Appendix A.

Proposition 4**.**

Let $f$ have wavelet expansion as in (3), and let $0\leq j\leq J$. If $n\geq 2^{dJ}$ and $p\geq 1$, then

[TABLE]

We are now in a position to prove Theorem 1.

Proof of Theorem 1.

Denote by $f_{J}$ the projection of $f$ to $\mathrm{span}\left(\Phi\cup\left(\bigcup_{j\leq J}\Psi_{j}\right)\right)$, i.e..

[TABLE]

The assumption that $f\in\mathcal{B}^{s}_{p^{\prime},q}(L)$ implies by lemma E.1 that

[TABLE]

for all $j\geq 0$; hence

[TABLE]

Lemma 8 implies

[TABLE]

by 4, as long as $2^{J}\leq n^{1/d}$. Choose $J$ so that $2^{J}\asymp n^{\frac{1}{d+2s}}$. If $d\geq 3$, the last term in the sum dominates and is of the same order as $2^{-J(s+1)}=n^{-{\frac{1+s}{d+2s}}}$. If $d\leq 2$, the approximation term is negligible and the sum is of order $n^{-1/2}$ if $d=1$, or $n^{-1/2}\log n$ if $d=2$. ∎

4.2 An adaptive estimator

By applying a truncation technique, we can also obtain an adaptive procedure which achieves the optimal rate (up to logarithmic factors) for any $s\geq 0$.

For any $j\geq 0$, set

[TABLE]

where $K$ is a constant defined in lemma B.1 depending on the wavelet system. For each $j\geq 0$, define by $\tilde{\beta}_{j}$ the vector of empirical wavelet coefficients at level $j$, that is, the vector indexed by $\Psi_{j}$ whose coordinates are given by

[TABLE]

The threshold $\tau_{j}$ is chosen in such a way that $\|\tilde{\beta}_{j}-\beta_{j}\|_{\ell_{p}}\ll\tau_{j}$ for all $j\leq J$. (See lemma B.1.)

We let $J$ be such that $2^{J}\asymp n^{1/d}$, and set

[TABLE]

Finally, in order to obtain an estimator that is a probability density, we defined the projected estimator $\hat{f}^{\circ}$ by

[TABLE]

The analysis of $\hat{f}^{\circ}$ is standard, and follows ideas for adaptive estimation developed by Donoho et al. (1996). We defer the details to Supplemental Appendix B.

4.3 Lower bound

Our lower bound follows almost directly from the bound proved by Kerkyacharian and Picard (1992) to establish minimax rates for density estimation in $L_{p}$ over Besov spaces. By the monotonicity of the Wasserstein-$p$ distances in $p$, it suffices to prove the lower bound for $p=1$.

Proof of theorem 3.

Given an index $J$ to be specified and a vector $\varepsilon\in\{\pm 1\}^{|\Psi_{J}|}$, we write

[TABLE]

As long as $n^{-1/2}2^{J(s+d/2)}\lesssim 1$, the functions $f_{\varepsilon}$ all lie in $\mathcal{B}^{s}_{p^{\prime},q}(L;m)$.

Moreover, note that for any $\varepsilon,\varepsilon^{\prime}\in\{\pm 1\}^{|\Psi_{J}|}$, 3 implies

[TABLE]

where $\rho(\varepsilon,\varepsilon^{\prime})$ is the Hamming distance between $\varepsilon$ and $\varepsilon^{\prime}$.

Moreover, when $\rho(\varepsilon,\varepsilon^{\prime})=1$, the Hellinger distance satisfies

[TABLE]

Therefore, since $W_{1}$ is a metric, a standard application of Assouad’s lemma (Tsybakov, 2009, Theorem 2.12) implies that

[TABLE]

where the infimum is taken over all estimators $\tilde{\mu}$ constructed from $n$ samples and where the final inequality is a consequence of lemma E.1. Choosing $J$ such that $2^{J}\asymp n^{\frac{1}{d+2s}}$ when $d\geq 2$ and $J=0$ when $d=1$ yields the claim. ∎

5 General smooth densities

In this section, we give proofs for our theorems applying to general smooth densities (theorems 5 and 6). Our main results are matching upper and lower bounds showing showing that the rate of estimation over the class $\mathcal{B}^{s}_{p,q}(L)$ is strictly worse than the rate over the class $\mathcal{B}^{s}_{p,q}(L;m)$ when $L$ is large enough that $\mathcal{B}^{s}_{p,q}(L)\not\subseteq\mathcal{B}^{s}_{p,q}(L;m)$.

5.1 Lower bounds

We assume that $L$ is large enough that $\mathcal{B}^{s}_{p^{\prime},q}(L)$ contains a function $g_{0}$ whose support lies entirely inside $(0,1/3)^{d}$. It is easy to see that this goal is indeed achievable by choosing $g_{0}$ to be suitable compactly supported smooth bump functions, as long as $L$ is a large enough constant.

The lower bound is based on the following lemma, which gives a lower bound on the Wasserstein distances for a pair of measures with disconnected support.

Lemma 9**.**

Let $\mu$ and $\nu$ be measures on $\mathbb{R}^{d}$. Suppose there exist two compact sets $S$ and $T$ such that $d(S,T)\geq c$ and such that the supports of $\mu$ and $\nu$ lie in $S\cup T$. Then $W_{p}(\mu,\nu)\geq c|\mu(S)-\nu(S)|^{1/p}$.

Proof.

Assume without loss of generality that $\mu(S)\geq\nu(S)$. Then any coupling between $\mu$ and $\nu$ must assign mass at least $\mu(S)-\nu(S)$ to $S\times T$, so that $W_{p}^{P}(\mu,\nu)\geq c^{p}|\mu(S)-\nu(S)|$. ∎

The proof of theorem 5 boils down to applying lemma 9 to appropriately chosen measures.

Proof of theorem 5.

We first prove the $n^{-1/2p}$ bound. Let $g_{0}\in\mathcal{B}^{s}_{p^{\prime},q}$ be supported in $[0,1/3]^{d}$ and let $g_{1}$ be a translation of $g_{0}$ supported on $[2/3,1]^{d}$. For $\lambda\in[-1,1]$, define $f_{\lambda}:=\frac{1}{2}((1+\lambda)g_{0}+(1-\lambda)g_{1})$. Then for any $\lambda\in[-1,1]$, the densities $f_{-\lambda}$ and $f_{\lambda}$ satisfy

[TABLE]

On the other hand, by lemma 9, choosing $S=[0,1/3]^{d}$ yields $W_{p}(f,{f^{\prime}})\gtrsim\lambda^{1/p}$. Therefore, if we choose $\lambda\asymp n^{-1/2}$, then the claim follows from the method of LeCam (1973).

We now prove the $n^{-\frac{1+s/p}{d+s}}$ bound. We proceed via Assouad’s lemma.

Let $g_{0}$ be as above. For $M>1$ to be specified, the characterization of $\mathcal{B}^{s}_{p^{\prime},q}$ by $\|\cdot\|^{\prime}_{\mathcal{B}^{s}_{p^{\prime},q}}$ implies that there exists a universal constant $c>0$ such that $h(x):=cM^{-s}g_{0}(Mx)$ also lies in $\mathcal{B}^{s}_{p^{\prime},q}(L)$. We denote by $\Gamma$ a set of vectors in $\mathbb{R}^{d}$ such that for any $\gamma_{1},\gamma_{2}\in\Gamma$, the supports of the functions $h(x-\gamma_{1})$, $h(x-\gamma_{2})$, and $g_{0}$ are all separated by at least $cM^{-1}$ for $c$ a small constant. By a volume argument, we can choose $\Gamma$ such that $|\Gamma|\asymp M^{d}$. We assume that $|\Gamma|$ is even.

We divide the elements of $|\Gamma|$ into pairs and label them $\{(\gamma^{+}_{i},\gamma^{-}_{i})\}_{i=1}^{|\Gamma|/2}$. For any $\varepsilon\in\{\pm 1\}^{|\Gamma|/2}$, define

[TABLE]

where $\kappa:=1-\frac{|\Gamma|}{2}\int h(x)\,\mathrm{d}x$ is chosen to ensure that $f_{\varepsilon}$ integrates to $1$. Since $\int h(x)\,\mathrm{d}x\asymp M^{-s-d}$, the constant $\kappa$ is positive and $f_{\varepsilon}$ is a density.

Given $\varepsilon,\varepsilon^{\prime}\in\{\pm 1\}^{|\Gamma|/2}$, define

[TABLE]

This is the subset of $\Gamma$ present in the density $f_{\varepsilon}$ but not in $f_{\varepsilon^{\prime}}$. We set

[TABLE]

If we denote by $\rho$ the Hamming distance, then the density $f_{\varepsilon}$ assigns mass $|\Delta(\varepsilon,\varepsilon^{\prime})|cM^{-s-d}\gtrsim\rho(\varepsilon,\varepsilon^{\prime})M^{-s-d}$ to $S$, and $f_{\varepsilon}^{\prime}$ assigns zero mass to this set. By construction, the rest of the support of $f_{\varepsilon}$ and $f_{\varepsilon}^{\prime}$ lies at distance at least $cM^{-1}$ from $S$. Therefore, by lemma 9,

[TABLE]

where the last inequality follows from the fact that $\rho(\varepsilon,\varepsilon^{\prime})\leq|\Gamma|/2\lesssim M^{d}$. Moreover, if $\rho(\varepsilon,\varepsilon^{\prime})=1$, then

[TABLE]

Therefore, if we choose $M\asymp n^{\frac{1}{s+d}}$, then Assouad’s lemma (Tsybakov, 2009, Theorem 2.12) and the fact that $W_{p}$ satisfies the triangle inequality imply

[TABLE]

as claimed. ∎

5.1.1 Upper bounds

When $p=1$, there is no longer any dependence on the lower bound $m$ in theorem 4, which suggests that the optimal rate of estimation in $W_{1}$ for general densities on $[0,1]^{d}$ should correspond with the one obtained in theorem 1 for bounded densities. And indeed, this is true: the results of Uppal et al. (2019) imply that the minimax rates of estimation for densities in $W_{1}$ is the same irrespective of whether the density is bounded below or not. However, the duality argument that Uppal et al. (2019) employ breaks down when $p\neq 1$, and the gap between the bounds in theorem 5 and theorem 1 hints that the situation when $p\geq 2$ is significantly more subtle.

In this section, we give an estimation matching the rate of theorem 5 rate up to a logarithmic factor when $s\in[0,1)$ and $f\in\mathcal{C}^{s}(L)$. Our estimator can only achieves the optimal rate for smoothness $s<1$ because we rely on particular properties of the Haar wavelet basis. While this basis is well suited to characterizing $\mathcal{C}^{s}(L)$ when $s<1$, it cannot exploit further regularity of the unknown density. We suspect that a similar result holds for other, more regular choices of wavelet systems, but our techniques do not allow us to prove a more general result.

For $j\geq 0$, denote by $\Psi_{j}$ the elements of the $d$-dimensional Haar wavelet basis at scale $2^{-j}$ (see Triebel, 2010, Section 2.3). Note that in the case of the Haar wavelet basis over the cube, the set $\Phi$ of scaling functions contains only the function which is identically $1$ on $[0,1]^{d}$.

For $j\geq 0$, let $\mathcal{Q}:=\bigcup_{j\geq 0}\mathcal{Q}_{j}$ be the dyadic decomposition of $[0,1]^{d}$, where $\mathcal{Q}_{j}$ consists of a partition of $[0,1]^{d}$ into cubes with sides of length $2^{-j}$. Each element of $\Psi_{j}$ is supported on a single cube in $\mathcal{Q}_{j}$. We write $V_{j}$ for the span of $\left\{1\cup\left(\bigcup_{k<j}\Psi_{k}\right)\right\}$ and denote by $K_{j}$ the orthogonal projection onto $V_{j}$. The functions in $V_{j}$ are precisely those which are constant on the elements of $\mathcal{Q}_{j}$.

The Haar basis has been used implicitly in existing bounds on the Wasserstein distance in what is known as the “dyadic partitioning argument” (Weed and Bach, 2018; Fournier and Guillin, 2015; Boissard and Le Gouic, 2014). Suppose that $\mu$ and $\nu$ are two measures on $[0,1]^{d}$. The bound of Weed and Bach (2018, Proposition 1) reads

[TABLE]

When $\mu$ and $\nu$ possess densities $f$ and $g$, respectively, the expression on the right side of the above inequality is an expansion of $f-g$ with respect to the Haar wavelet basis, since $|\mu(Q)-\nu(Q)|=\|K_{j}(f-g)\|_{L_{1}(Q)}$. Written in this way, the bound (8) therefore implies

[TABLE]

To achieve the optimal rate, we sharpen this bound significantly by observing that, when $f$ is bounded above and below on $Q$, the $L_{1}$ bound can be replaced by an $L_{p}$ bound. We present the refined result as theorem 10.

Theorem 10**.**

Let $J\geq 0$ be an integer, and let $f$ and $g$ be probability densities satisfying $g/f\in V_{J}$ (with the convention $0/0=0$). Then

[TABLE]

where

[TABLE]

The quantities $m_{Q}$ and $\Delta_{Q}$ control the variation of $f$ on $Q$. The first is a weighted version of the quantity $m$ appearing in 1. The second is small as long as $f$ does not change too much on $Q$, which is the case as long as $f$ is sufficiently smooth. The proof of theorem 10 appears in Supplemental Appendix D.

In practice, to apply theorem 10 to bound the Wasserstein distance between arbitrary probability measures $\mu$ and $\nu$, we can find a measure $\bar{\nu}$ such that $\mu$ and $\bar{\nu}$ satisfy the assumptions of theorem 10 and $W_{p}(\bar{\nu},\nu)$ is small. We adopt this strategy in the proof of the theorem 6, which appears in Supplemental Appendix C.

6 Computational aspects

One of the motivations for this line of work is found in applications of optimal transport techniques for data analysis and machine learning, with unknown distributions and access to an independent sample of size $n$.

Many so-called variational Wasserstein problems involve the problem of minimizing a functional $F:\nu\mapsto W_{p}(\nu,\mu)$ with unknown $\mu$. These problems, such as minimum Kantorovich estimators (Bassetti et al., 2006) and Wasserstein barycenters (Agueh and Carlier, 2011), are increasingly common in practical applications (Peyré and Cuturi, 2017), especially when the minimization is taken over a parametric class, with $\nu=\nu_{\theta}$ for $\theta\in\Theta$.

Solving variational Wasserstein problems in practice requires first obtaining an empirical estimate of the functional $F$ on the basis of data drawn from $\mu$, and then writing the resulting optimization problem in a computationally tractable form. The first issue is typically addressed by obtaining an estimator $\tilde{\mu}_{n}$ of $\mu$ and estimating the function via the plug-in principle. Indeed, the triangle inequality implies that

[TABLE]

where the supremum is taken over all probability distributions and where equality is achieved at $\mu=\nu$. Following this approach, guarantees in Wasserstein distance between $\mu$ and the estimator $\tilde{\mu}_{n}$ therefore yield uniform deviation bounds for these functionals over the set of probability measures $\mathcal{P}$ on $\mathbb{R}^{d}$.

To solve the resulting optimization problem, finite discretizations are often taken for $\nu$ and $\tilde{\mu}_{n}$ to render the resulting problem amenable to discrete optimization techniques (Cuturi, 2013; Altschuler et al., 2017). For this reason, the estimator $\tilde{\mu}_{n}$ is often taken to be the empirical distribution $\hat{\mu}_{n}$, since this measure is a finitely supported measure and enjoys the rate

[TABLE]

as long as $d>2p$, which is minimax optimal over the class of compactly supported probability measures (Singh and Póczos, 2018).

However, in sections 4 and 5, we show that under natural regularity assumptions for $\mu$, estimators based on density estimation statistically outperform the empirical distribution. Indeed, focusing on the regime $d>2p$, theorems 1 and 6 yield guarantees of the form

[TABLE]

where $\gamma^{*}(s)\geq 1$ increases as the smoothness of $\mu$ increases.

These results are summarized in the following table, highlighting that the optimal exponent $\gamma^{*}(s)/d$ interpolates between $1/d$ (for $s=0$) a dimension-free rate (for $s$ going to $\infty$). The value $s<1$ is of special interest, as it corresponds to the smoothness which can be exploited by a histogram estimator, which is most relevant for computational aspects. We take this point up further below.

[TABLE]

Note that the exponent $\gamma^{*}(s)$ is greater than 1 for all $s>0$. In the light of these results, there is an apparent tension between two objectives: statistical precision and computational efficiency. On the one hand, using atomic mesures as estimators of unknown distributions allows to efficiently compute Wasserstein distances: the optimal transport problem reduces to a linear program in finite dimension. Their statistical performance can however be suboptimal, as shown in the table above: for smooth densities the lower bound in $n^{-1/d}$ applies to all $n$-atomic distributions (Dudley, 1969). On the other hand, estimators with smooth densities have minimax optimal statistical precisions, but are defined through a projection in Besov norm of a wavelet-based density estimator. Even in the simple case of histograms (piecewise-constant densities), there is no explicit or simple way to solve optimal transport problems.

We therefore propose a procedure to leverage the regularity of the distribution, and to handle our proposed estimators in a computationally efficient manner. The idea is to exploit the best of both worlds, by creating an atomic measure with $M>n$ atoms, based on a density estimator. They statistically outperform the empirical measure, and optimal transport problems can be explicitly solved on these measures. If it is possible to efficiently sample from one of these estimators, it is always possible to extract an atomic distribution out of it, as described in the following

Definition 1**.**

Let $\tilde{\mu}_{n}$ be a probability measure from which one can draw sample points, and let $Z_{1},\ldots,Z_{M}$ be an i.i.d. sample from $\tilde{\mu}_{n}$. We denote by $\bar{\mu}_{n,M}$ the estimator resample distribution, the empirical distribution of the $Z_{i}$s.

The distribution $\bar{\mu}_{n,M}$ is “simply” the empirical distribution of a sample of size $M$ from a distribution. However, we retain $n$ in the notation, to highlight that the $Z_{i}$s are themselves drawn from an estimator based on a sample of size $n$ from an unknown $\mu$. We recall the following result for compactly supported distributinos.

Proposition 5** (Fournier and Guillin, 2015).**

For $d>2p$, the estimator resample distribution $\bar{\mu}_{n,M}$ satisfies ${\rm I}\kern-1.79993pt{\rm E}W_{p}(\tilde{\mu}_{n},\bar{\mu}_{n,M})\lesssim M^{-1/d}.$

As a consequence of this result, resampling from the estimated distribution yields an atomic measure as close in Wasserstein distance to the original estimator as desired, since $M$ can be chosen by the statistician. As the following makes clear, by choosing $M>n$, it is possible to obtain atomic distributions with estimation error as good as the original estimator $\tilde{\mu}_{n}$.

Corollary 1**.**

Assume $d>2p$. Let $\mu$ be in a nonparametric class such that there exists an estimator $\tilde{\mu}_{n}$ from which one can efficiently sample, and such that

[TABLE]

For any $\gamma\in[1,\gamma^{*}]$, the estimator resample distribution $\bar{\mu}_{n,M}$ with $M=n^{\gamma}$ satisfies

[TABLE]

Proof.

We have by the triangle inequality that

[TABLE]

having taken $M$ of order $n^{\gamma}$. ∎

It is now straightforward to give a proof of theorem 7.

Proof of theorem 7.

We consider the estimator resample distribution constructed from the histogram estimator give by theorem 6. If $d\leq 2p$, then Fournier and Guillin (2015, Theorem 1) implies that the empirical distribution $\hat{\mu}_{n}$ already achieves the rate appearing in theorem 6 up to a logarithmic factor; therefore, choosing $M=n$ suffices. If $d>2p$, we apply 1. Theorem 6 implies that, up to logarithmic factors, the histogram estimator $\mu_{\hat{f}}$ achieves the rate $n^{-\gamma^{*}/d}$ with $\gamma^{*}=\frac{1+\frac{s}{p}}{1+\frac{s}{d}}\wedge\frac{d}{2p}$. Choosing $\gamma=\gamma^{*}$ and noting that $\gamma^{*}<2$ yields the claim.

Finally, to show that constructing $\bar{\mu}_{n,M}$ takes time $O(M)$, we note that the histogram estimator constructed in theorem 6 is piecewise constant with $O(n)$ pieces. Constructing the histogram estimator takes $O(n)$ time, and sampling $M$ points from a histogram distribution can be done in time $O(n+M)=O(M)$ by the alias method (Kronmal and Peterson, 1979). ∎

In some examples, the choice of $\gamma^{*}$ can be left to the practitioner. For our estimators, this corresponds to choosing the depth of the wavelet decomposition. Taking piecewise constant estimators (histograms) limits the exponent $\gamma$ to $\gamma^{*}=\gamma^{*}(1)$. In any case, it is also possible to chose $M=n^{\gamma}$ for $\gamma\in(1,\gamma^{*}(1)]$, and let the approximation error $n^{-\gamma/d}$ dominate the statistical error $n^{-\gamma^{*}(1)/d}$.

Using the estimator resample distribution $\bar{\mu}_{n,M}$ instead of $\hat{\mu}_{n}$ requires solving optimal transport problems of size $M=n^{\gamma}$ instead of $n$. This naturally increases the computational cost. This motivates the question of quantifying the statistical and computational tradeoffs of our proposal. The dependency of the algorithmic cost of solving optimal transportation problems on the size of the distribution is the subject of a large literature (see Peyré and Cuturi, 2017), from which we extract a simple bound.

Proposition 6** (Altschuler et al., 2017; Dvurechensky et al., 2018).**

Given two distributions $\alpha$ and $\beta$ supported on at most $M$ atoms on a set of diameter $1$, an additive $\varepsilon$ approximation to $W^{p}_{p}(\alpha,\beta)$ can be computed in time $O(M^{2}\log(M)/\varepsilon^{2})$.

The following describes the interplay between statistical precision and computational efficiency for

Theorem 11**.**

Let $\mu$ be in a nonparametric class such that there exists an estimator $\tilde{\mu}_{n}$ from which one can efficiently sample, and such that

[TABLE]

Given a sample of size $n$ from $\mu$ and known $\nu$, for any $\gamma\in[1,\gamma^{*}]$ an estimate $\tilde{W}_{p,n}$ of $W_{p}(\nu,\mu)$ satisfying

[TABLE]

can be computed in time O\big{(}n^{\gamma(2+2/d)}\log(n)\big{)}.

Proof.

Taking $\tilde{\mu}_{n}$ such that ${\rm I}\kern-1.79993pt{\rm E}W_{p}(\mu,\tilde{\mu}_{n})\leq n^{-\gamma^{*}/d}$, $\bar{\mu}_{n,M}$ the empirical resample distribution of $\tilde{\mu}_{n}$, and $\bar{\nu}_{M}$ an $M$-atomic version of $\nu$ (obtained by sampling from or discretizing $\nu$), we have

[TABLE]

by taking $M=n^{\gamma}$. Taking $\varepsilon=n^{-\gamma/d}$, computing an $\varepsilon$-approximation $\tilde{W}_{p,n}^{p}$ to $W^{p}_{p}(\bar{\nu}_{M},\bar{\mu}_{n,M})$ given $\bar{\mu}_{n,M}$ and $\bar{\nu}_{M}$ takes time O\big{(}n^{\gamma(2+2/d)}\log(n)\big{)} by 6. Finally, the inequality

[TABLE]

yields the claim. ∎

Taking $\gamma=1$, $\tilde{\mu}_{n}=\hat{\mu}_{n}$ and $M=n$ with $\bar{\mu}_{n,M}=\hat{\mu}_{n}$ (without resampling) is always possible. It yields an algorithm that outputs a $n^{-1/d}$ approximation in time O\big{(}n^{2+2/d}\log(n)\big{)}. However, whenever another estimator $\tilde{\mu}_{n}$ with precision $n^{-\gamma/d}$ exists for $\gamma\in(1,\gamma^{*}(s)]$, it is possible to obtain a better approximation with error $\varepsilon_{n}(\gamma)=n^{-\gamma/d}$ in time of order $T_{n}(\gamma)=n^{\gamma(2+2/d)}$, up to terms logarithmic in $n$. This quantifies the computational cost for added statistical precision. We summarize these results in the following table, for $\gamma^{*}=\gamma^{*}(1)$ for histogram estimators (from which it is easy to sample). The parameter $\gamma$ can be taken in the full range from $1$ to $\gamma^{*}(1)$.

[TABLE]

In the high-dimensional limit, we obtain that a histogram estimator can improve the exponent in the precision by a factor $\gamma$ of nearly 2 at the price of increasing the exponent in the running time by nearly the same factor. The choice of $M$, which can be left to the statistician, determines the value of $\gamma$.

Appendix A Omitted proofs

A.1 Proof of lemma 1

The bound is invariant under the reparametrization $t\mapsto 1-t$, so it suffices to prove the claim for $t\leq 1/2$. By definition, we have

[TABLE]

where we have used the fact that $1-\lambda(t)\geq\lambda(t)$ when $t\leq 1/2$. Since $f(x)+g(x)\geq f(x)\vee g(x)\geq m$, the claim follows. ∎

A.2 Proof of lemma 2

The partial derivatives of this vector field satisfy

[TABLE]

so that $\nabla\cdot V=\phi$ on the interior of $[0,1]^{d}$.

On the boundary $x_{i}=0$ we have $(V_{\phi})_{i}=0$ for all $i$, and for $i>1$ on the boundary $x_{i}=1$ we have

[TABLE]

Finally, on $x_{1}=1$, we have $(V_{\phi})_{1}=\int_{0}^{1}\phi^{(1)}(t_{1})\,\mathrm{d}t_{1}=\int_{[0,1]^{d}}\phi(t)\,\mathrm{d}t$. Since $f$ and $g$ are both probability densities and the constant function $1$ is in the span of $\Phi$, we have $\int_{[0,1]^{d}}f(x)\,\mathrm{d}x=\int_{[0,1]^{d}}\sum_{\phi\in\Phi}\alpha_{\phi}\phi(x)\,\mathrm{d}x=1$, and, analogously, $\int_{[0,1]^{d}}\sum_{\phi\in\Phi}\alpha^{\prime}_{\phi}\phi(x)\,\mathrm{d}x=1$. We therefore obtain that $\sum_{\phi\in\Phi}(\alpha_{\phi}-\alpha^{\prime}_{\phi})(V_{\phi})_{1}=0$ on the face $x_{1}=0$; thus, the boundary conditions hold as well. ∎

A.3 Proof of lemma 3

By construction, $\nabla\cdot V_{\psi}=\frac{d(V_{\psi})_{k}}{dx_{k}}=\psi$ on the interior of $[0,1]^{d}$. The equality $(V_{\phi})_{i}=0$ holds on the boundaries $x_{i}=0$ and $x_{i}=1$ for $i\neq k$, and by construction $(V_{\phi})_{k}=0$ on $x_{k}=0$. Finally, since $\int_{0}^{1}\psi_{k}(t)\,\mathrm{d}t=0$, we also have $(V_{\phi})_{k}=0$ on $x_{1}=1$, and thus the boundary conditions are satisfied.∎

A.4 Proof of lemma 4

We will show that $\left\|\sum_{\phi\in\Phi}\alpha_{\phi}V_{\phi}\right\|_{L_{p}([0,1]^{d})}\lesssim\left\|\sum_{\phi\in\Phi}\alpha_{\phi}\phi\right\|_{L_{p}([0,1]^{d})}$ and conclude by appealing to the stability property of the scaling functions.

The definition of $V_{\phi}$ implies that

[TABLE]

where in the second inequality we use Jensen’s inequality and the fact that $x_{i}\leq 1$.

We obtain

[TABLE]

The claim then follows from lemma E.2.∎

A.5 Proof of lemma 5

By Assumption E.4, there exists an interval $I$ with $|I|\lesssim 2^{-j}$ such that $\operatorname*{supp}(\psi_{i})\subseteq I$. Since $(V_{\psi})_{i}(x_{i})=0$ if $x\notin I$, Hölder’s inequality implies that

[TABLE]

We therefore obtain by Assumption E.5 that

[TABLE]

The construction of $V_{\psi}$ and Assumption E.4 imply that $V_{\psi}(x)=0$ if $x\notin I_{\psi}$, where $\left\|\sum_{\psi\in\Psi_{j}}{\rm 1}\kern-2.40005pt{\rm I}\{x\in I_{\psi}\}\right\|_{L_{\infty}}\lesssim 1$. Hölder’s inequality therefore yields

[TABLE]

∎

A.6 Proof of lemma 7

This follows directly from the assumptions on the multiresolution analysis:

[TABLE]

where the inequalities follow, respectively, from the triangle inequality, Assumption E.6, lemma E.2, and the assumption on $\kappa$ and $\lambda$. ∎

A.7 Proof of 4

We first show the claim for $\|\beta_{j}-\tilde{\beta}_{j}\|_{\ell_{p}}$, $j\geq 0$. The inequalities of Rosenthal (1972) imply that there exists a constant $c_{p}$ such that for any $\psi\in\Psi_{j}$,

[TABLE]

where $\sigma_{\psi}^{2}:={\rm I}\kern-1.79993pt{\rm E}\left|\psi(X)-{\rm I}\kern-1.79993pt{\rm E}\psi(X)\right|^{2}$ and $X\sim f$.

Assumption E.5 implies that $\|\psi\|_{L_{\infty}}\lesssim 2^{dj/2}$, so

[TABLE]

where the last inequality follows from the assumption that $n\geq 2^{dj}$.

In order to establish that ${\rm I}\kern-1.79993pt{\rm E}\|\beta_{j}-\hat{\beta}_{j}\|_{\ell_{p}}^{p}\leq\frac{2^{dj}}{n^{p/2}}$, it therefore suffices to show that

[TABLE]

for $k\in\{2,p\}$.

We have

[TABLE]

where in the last step we use Jensen’s inequality combined with the fact that $\int\psi^{2}=1$ for all $\psi\in\Psi_{j}$ (see Assumption E.1). Finally, we obtain

[TABLE]

The fact that $f\in\mathcal{B}^{s}_{p^{\prime},q}(L;m)$ where $1\leq p\leq p^{\prime}$ implies that $\|f\|^{k/2}_{L_{k/2}}\leq\|f\|^{p^{\prime}}_{L_{p^{\prime}}}\lesssim 1$, and Assumptions E.4 and E.5 imply $\left\|\sum_{\psi\in\Psi_{j}}\psi(x)^{2}\right\|_{L_{\infty}}\lesssim 2^{dj}$. This proves the claim.

To show the analogous bound on $\|\alpha-\tilde{\alpha}\|_{\ell_{p}}$, we again write

[TABLE]

for any $\phi\in\Phi$, where $\sigma^{2}:=\left|\phi(X)-{\rm I}\kern-1.79993pt{\rm E}\phi(X)\right|^{2}$. Since $\|\phi\|_{L_{\infty}}\lesssim 1$ for any $\phi\in\Phi$, this bound simplifies to

[TABLE]

and hence

[TABLE]

by lemma E.1. ∎

Appendix B Adaptivity bound

Let us first justify the definition of $\tau_{j}$, by showing that it gives an upper bound on $\|\tilde{\beta}_{j}-\beta_{j}\|_{\ell_{p}}$ with high probability.

Lemma B.1**.**

Suppose $f\in\mathcal{B}^{s}_{p^{\prime},q}(L)$ and $\tau_{j}$ is defined as in (7). For any $j\geq 0$, if $n\geq 2^{dj}$, it holds

[TABLE]

Moreover,

[TABLE]

Proof.

For any $\psi\in\Psi_{j}$, Bernstein’s inequality (Boucheron et al., 2013) yields

[TABLE]

where $\sigma_{\psi}^{2}$ denotes the variance of $\psi(X)$ when $X\sim f$. In particular, this yields that

[TABLE]

Taking a union bound over all $\psi\in\Psi_{j}$ yields

[TABLE]

As in the proof of 4, there exists a constant $K_{0}$ depending only on the wavelet basis such that $\sum_{\psi\in\Psi_{j}}\sigma_{\psi}^{p}\leq K_{0}^{p}L^{p}2^{dj}$. By Assumptions E.5 and E.1, there exists a constant $K_{1}$ depending on the wavelet basis such that $\left(\sum_{\psi\in\Psi_{j}}\|\psi\|_{L_{\infty}}^{p}\right)^{1/p}\leq K_{1}2^{dj(\frac{1}{2}+\frac{1}{p})}$. Likewise, lemma E.1 implies the existence of a constant $K_{2}$ such that $|\Psi_{j}|\leq K_{2}2^{dj}$.

Therefore, setting $K=\max\{K_{0},K_{1},K_{2}\}$ and using the assumption that $n\geq 2^{dj}$ implies

[TABLE]

Choosing $t=4cdj$ and applying the definition of $\tau_{j}$ yields (B.1).

The bound (B.2) follows from integrating the estimate (B.1). ∎

We now turn to the proof of theorem 2. By an argument analogous to the one given in lemma 8, it suffices to bound

[TABLE]

We require two simple lemmas bounding the expectations of each summand.

Lemma B.2**.**

For any $j\geq 0$,

[TABLE]

Proof.

We compute

[TABLE]

where the last inequality follows from lemma B.1, (B.2). Combining this estimate with the trivial bound ${\rm I}\kern-1.79993pt{\rm E}\|\beta_{j}\|_{\ell_{p}}\mathds{1}_{\|\tilde{\beta}_{j}\|_{\ell_{p}}<\tau_{j}}\leq\|\beta_{j}\|_{\ell_{p}}$ gives the claim.

∎

Lemma B.3**.**

For any $j\geq 0$,

[TABLE]

Proof.

If $\|\beta_{j}\|_{\ell_{p}}<\frac{1}{2}\tau_{j}$, we use the Cauchy-Schwarz inequality and lemma B.1 to obtain

[TABLE]

For the second bound, we simply use use ${\rm I}\kern-1.79993pt{\rm E}\|\tilde{\beta}_{j}-\beta_{j}\|_{\ell_{p}}\mathds{1}_{\|\tilde{\beta}_{j}\|_{\ell_{p}}<\tau_{j}}\leq{\rm I}\kern-1.79993pt{\rm E}\|\tilde{\beta}_{j}-\beta_{j}\|_{\ell_{p}}\lesssim\tau_{j}$. ∎

We are now ready to complete the proof.

Proof of theorem 2.

By lemmas B.2 and B.3, we have the bound

[TABLE]

First, the definition of $\tau_{j}$ implies that

[TABLE]

so $\sum_{j\geq 0}2^{-j}2^{dj(\frac{1}{2}-\frac{1}{p})}\tau_{j}p_{j}^{1/2}\lesssim n^{-1/2}$. Since the claimed rate is no faster than $n^{-1/2}$, this term is not the dominant one.

If $f\in\mathcal{B}^{s}_{p^{\prime},q}(L;m)$, then $2^{dj(\frac{1}{2}-\frac{1}{p})}\|\beta_{j}\|_{\ell_{p}}\lesssim 2^{-js}$, which implies

[TABLE]

Summing this expression yields

[TABLE]

∎

Appendix C An upper bound for densities not bounded below

We require two basic facts about Binomial random variables; the first is a direct consequence of the inequalities of Rosenthal (1972) and the second is a standard Chernoff bound (Hagerup and Rüb, 1990).

Lemma C.1**.**

Let $X\sim\mathrm{Bin}(n,p)$. If $np\geq 1$, then

[TABLE]

Lemma C.2**.**

If $X\sim\mathrm{Bin}(n,p)$, then

[TABLE]

Throughout the proof, we let $\Psi_{j}$ for $j\geq 0$ denote the Haar wavelet basis. Fix a threshold $J$, which we will assume satisfies $n\asymp 2^{J(d+s)}$. We set

[TABLE]

and define

[TABLE]

For convenience, throughout the proof we write $\mu$ for the probability measure with density $f$ and $\nu$ for the probability measure with density $\hat{f}$. Note that for any $Q\in\mathcal{Q}_{j}$, $j\leq J$, we have $n\cdot\nu(Q)\sim\mathrm{Bin}(n,\mu(Q))$.

Let us now define a density $g$ by

[TABLE]

where we interpret $0/0=0$. Since $\mu(Q)=0\implies\nu(Q)=0$ almost surely, this density is well defined with probability $1$. By construction, $g/f$ is constant on each element of $\mathcal{Q}_{J}$, so $g/f\in V_{J}$.

Let $0\leq j<J$, and fix $Q\in\mathcal{Q}_{j}$. Theorem 10 implies that the contribution to $W_{p}(f,g)$ from $Q\in\mathcal{Q}_{j}$ can be bounded by either the $L_{1}$ bound $\|K_{j+1}(f-g)\|_{L_{1}(Q)}$ or the $L_{p}$ bound $m_{Q}^{1-p}\Delta_{Q}^{p}\|K_{j+1}(f-g)\|^{p}_{L_{p}(Q)}$. We can always apply the first bound, but we will apply the second when the density $f$ is bounded below on $Q$. We denote by $\mathcal{E}_{Q}$ the high-probability event that $\nu(Q)\geq\frac{1}{2}\mu(Q)$.

We divide into cases. As we shall see, the dominant term arises from the first case; the latter two are asymptotically negligible.

Case 1: $\inf_{x\in Q}f(x)<2^{-js}$

Under the assumption that $f\in\mathcal{C}^{s}$, we have that $\sup_{x\in Q}\leq\inf_{x\in Q}+\operatorname*{diam}(Q)^{s}\lesssim 2^{-js}$. Therefore $\mu(Q)=\int_{Q}f(x)\,\mathrm{d}x\lesssim 2^{-j(d+s)}$.

To evaluate $K_{j+1}(f-g)$, we note that for each $R\subseteq Q,R\in\mathcal{Q}_{j+1}$, the densities $K_{j+1}f$ and $K_{j+1}g$ are constant on R, with value $2^{d(j+1)}\mu(R)$ and $2^{d(j+1)}\nu(R)$. This implies

[TABLE]

In particular, recalling that $n\cdot\nu(R)\sim\mathrm{Bin}(n,\mu(R))$, we have that that

[TABLE]

Jensen’s inequality along with the fact that $\sup_{x\in Q}f(x)\leq 2^{-js}$ therefore yields

[TABLE]

Case 2: $\inf_{x\in Q}f(x)\geq 2^{-js}$ and $\mathcal{E}_{Q}$ holds

Under the assumption that $f\in\mathcal{C}^{s}$, we have that $\sup_{x\in Q}f(x)-\inf_{x\in Q}f(x)\leq\operatorname*{diam}(Q)^{s}\lesssim 2^{-js}$. We therefore have $\Delta_{Q}\lesssim 1$. Moreover, on $\mathcal{E}_{Q}$, we have

[TABLE]

Arguing as above, we have

[TABLE]

Since $\sup_{x\in Q}f(x)\asymp\inf_{x\in Q}f(x)$, we have $\mu(R)\asymp 2^{-jd}\inf_{x\in Q}f(x)$. Since $n\gtrsim 2^{J(d+s)}\geq 2^{j(d+s)}\geq 1/\mu(R)$, we can apply lemma C.1 to obtain

[TABLE]

Combining these bounds, we obtain

[TABLE]

Case 3: $\mathcal{E}_{Q}$ does not hold

The Cauchy-Schwarz inequality and the bound (C.1) imply

[TABLE]

By lemma C.2,

[TABLE]

and combining these bounds yields

[TABLE]

We are now prepared to bound $W_{p}(f,g)$. Combining the three cases above, we observe

[TABLE]

Under the assumption that $n\geq 2^{j(d+s)}$, the first term dominates. We therefore have

[TABLE]

and theorem 10 combined with Jensen’s inequality implies

[TABLE]

On the other hand, a simple truncation argument (lemma C.3, below) shows that ${\rm I}\kern-1.79993pt{\rm E}W_{p}(g,\hat{f})\lesssim 2^{-J(1+s/p)}\asymp n^{-\frac{1+s/p}{d+s}}$.

Since $W_{p}(f,\hat{f})\leq W_{p}(f,g)+W_{p}(g,{\hat{f}})$, combining these two bounds yields the claim. ∎

All that remains is to give the proof of lemma C.3.

Lemma C.3**.**

For $g$ and $\hat{f}$ defined as in the proof of theorem 6, we have

[TABLE]

Proof.

By construction, $\hat{f}$ is an element of $V_{J}$. Moreover, $g$ is defined so that $K_{J}g=K_{J}\hat{f}$. Combining these two facts yields that $\hat{f}=K_{J}g$. For $j\leq J$, this implies that $K_{j}\hat{f}=K_{j}g$, and for $j>J$, it yields

[TABLE]

where in the last step we use the definition of $g$. Since ${\rm I}\kern-1.79993pt{\rm E}\nu(Q)=\mu(Q)$, we obtain

[TABLE]

For $j>J$, the operator $K_{j}-K_{J}$ is a projection onto the span of $\{\Psi_{k}\}_{J\leq k<j}$. Under the assumption that $f\in\mathcal{C}^{s}=\mathcal{B}_{\infty,\infty}^{s}$, we obtain

[TABLE]

where the final inequality uses lemma E.2.

Therefore, applying the bound of Weed and Bach (2018) in the equivalent form (9), we obtain

[TABLE]

as desired. ∎

Appendix D Proof of Theorem 10

The assumption that $g/f\in V_{J}$ implies that, for any $j\leq J$, if $\int_{Q}f=0$, then $\int_{Q}g=0$. For $0\leq j\leq J$, we therefore define $f_{j}$ by setting

[TABLE]

where we interpret $0/0=0$. With this definition, $f_{0}=f$ and $f_{J}=g$. Moreover, for $0\leq j\leq J$, we have $K_{j}f_{j}=K_{j}f_{j+1}$.

The triangle inequality therefore implies

[TABLE]

We focus on bounding the quantity $W_{p}({f_{j}},{f_{j+1}})$, from which the claimed bound will follow.

We require the following lemma, whose proof is deferred.

Lemma D.1**.**

Fix $j\geq 0$. Let $p\geq 1$, and let $f_{j}$ and $f_{j+1}$ be any two densities such that $K_{j}f_{j}=K_{j}f_{j+1}$. Then

[TABLE]

By the definition of $f_{j}$, we have $\inf_{x\in Q}f_{j}(x)=\frac{\int_{Q}g}{\int_{Q}f}\inf_{x\in Q}f(x)=m_{Q}$. The claim then follows from the following lemma.

Lemma D.2**.**

For any $Q\in\mathcal{Q}_{j}$, the densities $f_{j}$ and $f_{j+1}$ satisfy

[TABLE]

Combining (D.1), lemma D.1, and lemma D.2 yields theorem 10. ∎

It remains to give proofs of the two required lemmas.

D.1 Proof of lemma D.1

If $K_{j}f_{j}=K_{j}f_{j+1}$, then the measures ${f_{j}}$ and ${f_{j+1}}$ assign the same mass to each cube in $\mathcal{Q}_{j}$. We can therefore obtain a coupling between ${f_{j}}$ and ${f_{j+1}}$ by optimally coupling ${f_{j}}|_{Q}$ and ${f_{j+1}}|_{Q}$ on each $Q\in\mathcal{Q}_{j}$ and combining the resulting couplings (see Weed and Bach, 2018, Lemma A.2). Doing so yields the bound

[TABLE]

To show that $W_{p}^{p}({f_{j}}|_{Q},{f_{j+1}}|_{Q})\leq 2^{-kp}\|f_{j}-f_{j+1}\|_{L_{1}(Q)}$, we use the fact that for two measures with support in a set $K$, it holds $W_{p}^{p}(\mu,\nu)\leq\operatorname*{diam}(K)^{p}d_{\mathrm{TV}}(\mu,\nu)$, and when $\mu$ and $\nu$ have densities, the total variation $d_{\mathrm{TV}}$ reduces to the $L_{1}$ norm.

To obtain the second bound, we apply 1. By lemma E.2, if we denote by $P_{j}$ the $L_{2}$ projection onto the span of $\Psi_{j}$, then the bound of 1 reads

[TABLE]

where the second inequality follows from lemma E.3. ∎

D.2 Proof of lemma D.2

For convenience, let us write $\mu$ for the measure with density $f$ and $\nu$ for the measure with density $g$. We recall that on $R\in\mathcal{Q}_{j+1}$, the densities $K_{j+1}f$ and $K_{j+1}g$ are constant with values $2^{d(j+1)}\mu(R)$ and $2^{d(j+1)}\nu(R)$. For any cube $R\in\mathcal{Q}_{j+1}$ such that $R\subseteq Q$, we have that $(f_{j}-f_{j+1})=\left(\frac{\nu(Q)}{\mu(Q)}-\frac{\nu(R)}{\mu(R)}\right)f$ on $R$, which implies

[TABLE]

The first sum is bounded by $|\nu(Q)-\mu(Q)|+\sum_{R}|\nu(R)-\mu(R)|\leq 2\sum_{R}|\nu(R)-\mu(R)|=2\|K_{j+1}(f-g)\|_{L_{1}(Q)}$.

For the second, we write

[TABLE]

Since $R$ is a cube of side length $2^{-(j+1)}$, we have the simple bounds

[TABLE]

We obtain

[TABLE]

and applying the definition of $K_{j+1}f$ and $K_{j+1}g$ yields

[TABLE]

as desired. ∎

Appendix E Assumptions on wavelets

For completeness, we have extracted the properties we require of the sets $\Phi$ and $\Psi_{j}$.

Assumption E.1** (Basis).**

$\Phi\cup\left\{\bigcup_{j\geq 0}\Psi_{j}\right\}$* forms an orthonormal basis for $L_{2}([0,1]^{d})$.*

Assumption E.2** (Regularity).**

The functions in $\Phi$ and $\Psi_{j}$ for $j\geq 0$ all lie in $\mathcal{C}^{r}(\Omega)$, and polynomials of degree at most $r$ on $\Omega$ lie in $\mathrm{span}(\Phi)$.

Assumption E.3** (Tensor construction).**

Each $\psi\in\Psi_{j}$ can be expressed as $\psi(\mathbf{x})=\prod_{i=1}^{d}\psi_{i}(x_{i})$ for some univariate functions $\psi_{i}$.

Assumption E.4** (Locality).**

For each $\psi\in\Psi_{j}$ there exists a rectangle $I_{\psi}\subseteq[0,1]^{d}$ such that $\operatorname*{supp}(\psi)\subseteq I_{\psi}$, $\operatorname*{diam}(I_{\psi})\lesssim 2^{-j}$, and $\left\|\sum_{\psi\in\Psi_{j}}{\rm 1}\kern-2.40005pt{\rm I}\{x\in I_{\psi}\}\right\|_{L_{\infty}}\lesssim 1$.

Assumption E.5** (Norm).**

$\|\psi\|_{L_{p}(\Omega)}\asymp 2^{dj(\frac{1}{2}-\frac{1}{p})}$* for all $\psi\in\Psi_{j}$.*

Assumption E.6** (Bernstein estimate).**

$\|\nabla f\|_{L_{p}(\Omega)}\lesssim 2^{j}\|f\|_{L_{p}(\Omega)}$* for any $f\in\mathrm{span}\left(\Phi\cup\left\{\bigcup_{0\leq k<j}\Psi_{k}\right\}\right)$.*

These requirements are achievable for $\Omega=[0,1]^{d}$ due to a classic construction due to Cohen et al. (1993). (See also Cohen, 2003, Chapter 2.12, for further details.)

We state without proof some straightforward consequences of our assumptions.

Lemma E.1**.**

$|\Phi|\lesssim 1$* and $|\Psi_{j}|\lesssim 2^{dj}$ for $j\geq 0$.*

Lemma E.2**.**

For any vector $\{\alpha_{\phi}\}_{\phi\in\Phi}$, $\|\sum_{\phi\in\Phi}\alpha_{\phi}\phi\|_{L_{p}(\Omega)}\asymp\|\alpha\|_{\ell_{p}}$. Likewise, for any $\{\beta_{\psi}\}_{\psi\in\Psi_{j}}$, $\|\sum_{\psi\in\Psi_{j}}\beta_{\psi}\psi\|_{L_{p}(\Omega)}\asymp 2^{dj(\frac{1}{2}-\frac{1}{p})}\|\beta\|_{\ell_{p}}$.

Lemma E.3**.**

If $P_{j}$ denotes the orthogonal projection onto the span of $\Psi_{j}$, then $\|P_{j}f\|_{L_{p}(\Omega)}\lesssim\|f\|_{L_{p}(\Omega)}$.