Limit distribution theory for smooth $p$-Wasserstein distances

The Wasserstein distance is an instance of the Kantorovich optimal transport problem [Kan42], which defines a metric on a space of probability measures. Specifically, for $1\leq p<\infty$, the $p$-Wasserstein distance between two Borel probability measures $\mu$ and $\nu$ on $\mathbb{R}^{d}$ with finite $p$th moments is defined by

where $\Pi(\mu,\nu)$ is the set of couplings (or transportation plans) of $\mu$ and $\nu$. The Wasserstein distance has seen a surge of applications in statistics, machine learning, and applied mathematics, ranging from generative modeling [ACB17, GAA⁺17, TBGS18], image recognition [RTG00, SL11], and domain adaptation [CFT14, CFTR16] to robust optimization [GK16, MEK18, BMZ21] and partial differential equations [JKO98, San17]. The widespread applicability of the Wasserstein distance is driven by an array of desirable properties, including its metric structure ($\mathsf{W}_{p}$ metrizes weak convergence plus convergence of $p$th moments), a convenient dual form, robustness to support mismatch, and a rich geometry it induces on a space of probability measures. We refer to [Vil03, Vil08, AGS08, San15] as standard references on optimal transport theory.

However, statistical aspects of Wasserstein distances are bottlenecked by the curse of dimensionality, whereby the number of data points needed to accurately estimate them grows exponentially with dimension. Specifically, for the empirical distribution $\hat{\mu}_{n}$ of $n$ independent observations from a distribution $\mu$ on $\mathbb{R}^{d}$, it is known that $\operatorname{\mathbb{E}}\big{[}\mathsf{W}_{p}(\hat{\mu}_{n},\mu)\big{]}$ scales as $n^{-1/d}$ for $d>2p$ under moment conditions [DSS13, BLG14, FG15, WB19a, Lei20]. This slow rate renders performance guarantees in terms of $\mathsf{W}_{p}$ all but vacuous when $d$ is large. It is also a roadblock towards a more delicate statistical analysis concerning limit distributions, bootstrap, and valid inference.

Gaussian smoothing was recently introduced as a means to alleviate the curse of dimensionality of empirical $\mathsf{W}_{p}$ [GGNP20, GG20, GGK20, SGK21, NGK21]. For $\sigma>0$, the smooth $p$-Wasserstein distance is defined as $\mathsf{W}_{p}^{(\sigma)}(\mu,\nu):=\mathsf{W}_{p}(\mu*\gamma_{\sigma},\nu*\gamma_{\sigma})$, where $*$ denotes convolution and $\gamma_{\sigma}=N(0,\sigma^{2}I_{d})$ is the isotropic Gaussian distribution with variance parameter $\sigma^{2}$. For sufficiently sub-Gaussian $\mu$, [GGNP20] showed that the expected smooth distance between $\hat{\mu}_{n}$ and $\mu$ exhibits the parametric convergence rate, i.e., $\operatorname{\mathbb{E}}\big{[}\mathsf{W}^{(\sigma)}_{1}(\hat{\mu}_{n},\mu)\big{]}=O(n^{-1/2})$ in any dimension. This is a significant departure from the $n^{-1/d}$ rate in the unsmoothed case. [GG20] further showed that $\mathsf{W}^{(\sigma)}_{1}$ maintains the metric and topological structure of $\mathsf{W}_{1}$ and is able to approximate it within a $\sigma\sqrt{d}$ gap. The structural properties and fast empirical convergence rates were later extended to $p>1$ in [NGK21]. Other follow-up works explored relations between $\mathsf{W}_{p}^{(\sigma)}$ and maximum mean discrepancies [ZCR21], analyzed its rate of decay as $\sigma\to\infty$ [CNW21], and adopted it as a performance metric for nonparametric mixture model estimation [HMS21].

A limit distribution theory for $\mathsf{W}^{(\sigma)}_{1}$ was developed in [GGK20, SGK21], where the scaled empirical distance $\sqrt{n}\,\mathsf{W}^{(\sigma)}_{1}(\hat{\mu}_{n},\mu)$ was shown to converge in distribution to the supremum of a tight Gaussian process in every dimension $d$ under mild moment conditions. This result relies on the dual formulation of $\mathsf{W}_{1}$ as an integral probability metric (IPM) over the class of 1-Lipschitz functions. Gaussian smoothing shrinks the function class to that of 1-Lipschitz functions convolved with a Gaussian density, which is shown to be $\mu$-Donsker in every dimension, thereby yielding the limit distribution. Extending these results to empirical $\mathsf{W}_{p}^{(\sigma)}$ with $p>1$, however, requires substantially new ideas due to the lack of an IPM structure. Consequently, works exploring $\mathsf{W}_{p}^{(\sigma)}$ with $p>1$, such as [ZCR21, NGK21], did not contain limit distribution results for it and this question remained largely open.

The present paper closes this gap and provides a comprehensive limit distribution theory for empirical $\mathsf{W}_{p}^{(\sigma)}$ with $p>1$. Our main limit distribution results are summarized in the following theorem, where the ‘null’ refers to when $\mu=\nu$, while ‘alternative’ corresponds to $\mu\neq\nu$. In all what follows, the dimension $d\geq 1$ is arbitrary.

Parts (i) and (ii) show that the null limit distributions are non-Gaussian. On the other hand, Parts (iii) and (iv) establish asymptotic normality of empirical $\mathsf{W}_{p}^{(\sigma)}$ under the alternative. Notably, these result have the correct centering, $\mathsf{W}_{p}^{(\sigma)}(\mu,\nu)$, which enables us to construct confidence intervals for $\mathsf{W}_{p}^{(\sigma)}(\mu,\nu)$.

The proof strategy for Theorem 1.1 differs from existing approaches to limit distribution theory for empirical $\mathsf{W}_{p}$ for general distributions. In fact, an analog of Theorem 1.1 is not known to hold for classic $\mathsf{W}_{p}$ in this generality, except for the special case where $\mu,\nu$ are discrete (see a literature review below for details). The key insight is to regard $\mathsf{W}_{p}^{(\sigma)}$ as a functional defined on a subset of a certain dual Sobolev space. We show that the smooth empirical process converges weakly in the dual Sobolev space and that $\mathsf{W}_{p}^{(\sigma)}$ is Hadamard (directionally) differentiable w.r.t. the dual Sobolev norm. We then employ the extended functional delta method [Sha90, R0̈4] to obtain the limit distribution of one- and two-sample empirical $\mathsf{W}_{p}^{(\sigma)}$ under both the null and the alternative.

The limit distributions in Theorem 1.1 are non-pivotal in the sense that they depend on the population distributions $\mu$ and $\nu$, which are unknown in practice. To facilitate statistical inference using $\mathsf{W}_{p}^{(\sigma)}$, we employ the bootstrap to estimate the limit distributions and prove its consistency for each case of Theorem 1.1. Under the alternative, the consistency follows from the linearity of the Hadamard derivative. Under the null, where the Hadamard (directional) derivative is nonlinear, the bootstrap consistency is not obvious but still holds. This is somewhat surprising in light of [Düm93, FS19], where it is demonstrated that the bootstrap, in general, fails to be consistent for functionals whose Hadamard directional derivatives are nonlinear (cf. Proposition 1 in [Düm93] or Corollary 3.1 in [FS19]). Nevertheless, our application of the bootstrap differs from [Düm93, FS19] so there is no contradiction, and the specific structure of the Hadamard derivative of $\mathsf{W}_{p}^{(\sigma)}$ allows to establish consistency under the null (see the discussion after Proposition 3.8 for more details). These bootstrap consistency results enable constructing confidence intervals for $\mathsf{W}_{p}^{(\sigma)}$ and using it to test the equality of distributions.

As an application of the limit theory, we study implicit generative modeling under the minimum distance estimation (MDE) framework [Wol57, Pol80, PS80]. MDE extends the maximum-likelihood principle beyond the KL divergence and applies to models supported on low-dimensional manifolds [ACB17] (whence the KL divergence is not well-defined), as well as to cases when the likelihood function is intractable [GMR93]. For MDE with $\mathsf{W}_{p}^{(\sigma)}$, we establish limit distribution results for the optimal solution and the smooth $p$-Wasserstein error. Our results hold for arbitrary dimension, again contrasting the classic case where analogous distributional limits for MDE with $\mathsf{W}_{p}$ are known only for $p=d=1$ [BJGR19]. Remarkably, when $p=2$, the Hilbertian structure of the underlying dual Sobolev space allows showing asymptotic normality of the MDE solution.

Analysis of empirical Wasserstein distances, or more generally empirical optimal transport distances, has been an active research area in the statistics and probability theory literature. In particular, significant attention was devoted to rates of convergence and exact asymptotics [Dud69, AKT84, Tal92, Tal94, DY95, BdMM02, BGV07, DSS13, BB13, BLG14, FG15, WB19a, AST19, Led19, BL21, Lei20, CRL⁺20, MNW21, DGS21, MBNWW21]. As noted before, the empirical Wasserstein distance suffers from the curse of dimensionality, namely, $\operatorname{\mathbb{E}}[\mathsf{W}_{p}(\hat{\mu}_{n},\mu)]=O(n^{-1/d})$ whenever $d>2p$. This rate is known to be sharp in general [Dud69]. The recent work by [CRL⁺20, MNW21] discovered that the rate can be improved under the alternative, namely, $\operatorname{\mathbb{E}}\big{[}\big{|}\mathsf{W}_{p}(\hat{\mu}_{n},\nu)-\mathsf{W}_{p}(\mu,\nu)\big{|}\big{]}=O(n^{-\alpha/d})$ for $d\geq 5$ if $\nu\neq\mu$, where $\alpha=p$ for $1\leq p<2$ and $\alpha=2$ for $2\leq p<\infty$. Their insight is to use the duality formula for $\mathsf{W}_{p}^{p}$ and exploit regularity of optimal transport potentials. [MNW21] also derive matching minimax lower bounds up to log factors under some technical conditions.

Another central problem that has seen a rapid development is limit distribution theory for empirical Wasserstein distances. However, except for the two special cases discussed next, to the best of our knowledge, there is no proven analog of our Theorem 1.1 for classic Wasserstein distances, i.e., a comprehensive limit distribution theory for empirical $\mathsf{W}_{p}$ that holds for general $d$ and $p$. The first case for which the limit distribution is well understood is when $d=1$. Then, $\mathsf{W}_{p}$ reduces to the $L^{p}$ distance between quantile functions for $1\leq p<\infty$, and further simplifies to the $L^{1}$ distance between distribution functions when $p=1$. Building on such explicit expressions, [dBGM99] and [dBGU05] derived null limit distributions in $d=1$ for $p=1$ and $p=2$, respectively. More recently, under the alternative ($\mu\neq\nu$), [DBGL19] derived a central limit theorem (CLT) when $d=1$ and $p\geq 2$. The second case where a limit distribution theory for empirical $\mathsf{W}_{p}$ is available is when $\mu,\nu$ are discrete. If the distributions are finitely discrete, i.e., $\mu=\sum_{j=1}^{m}r_{j}\delta_{x_{j}}$ and $\nu=\sum_{j=1}^{k}s_{j}\delta_{y_{j}}$ for two simplex vectors $r=(r_{1},\dots,r_{m})$ and $s=(s_{1},\dots,s_{k})$, then $\mathsf{W}_{p}(\mu,\nu)$ can be seen as a function of those simplex vectors $r$ and $s$. Leveraging this, [SM18] applied the delta method to obtain limit distributions for empirical $\mathsf{W}_{p}$ in the finitely discrete case. An extension to countably infinite supports was provided in [TSM19], while [dBGSL21a] treated the semidiscrete case where $\mu$ is finitely discrete but $\nu$ is general.

Except for these two special cases, limit distributions for Wasserstein distances are less understood. To avoid repetitions, we focus here our discussion on the one sample case. In [dBL19], a CLT for $\sqrt{n}\big{(}\mathsf{W}_{2}^{2}(\hat{\mu}_{n},\nu)-\operatorname{\mathbb{E}}\big{[}\mathsf{W}_{2}^{2}(\hat{\mu}_{n},\nu)\big{]}\big{)}$ is derived in any dimension, but the limit Gaussian distribution degenerates to $0$ when $\mu=\nu$; see also [dBGSL21b] for an extension to general $1<p<\infty$. Notably, the centering constant there is the expected empirical Wasserstein distance $\operatorname{\mathbb{E}}\big{[}\mathsf{W}_{2}^{2}(\hat{\mu}_{n},\nu)\big{]}$, which in general can not be replaced with the (more natural) population distance $\mathsf{W}_{2}^{2}(\mu,\nu)$. The recent preprint [MBNWW21] addressed this gap and established a CLT for $\sqrt{n}\big{(}\mathsf{W}_{2}^{2}(\tilde{\mu}_{n},\nu)-\mathsf{W}_{2}^{2}(\mu,\nu)\big{)}$ for a wavelet-based estimator $\tilde{\mu}_{n}$ of $\mu$, assuming that the ambient space is $[0,1]^{d}$ and that $\mu,\nu$ are absolutely continuous w.r.t. the Lebesgue measure with smooth and strictly positive densities. Following arguments similar to [dBL19], they first derive a CLT for $\sqrt{n}\big{(}\mathsf{W}_{2}^{2}(\tilde{\mu}_{n},\nu)-\operatorname{\mathbb{E}}\big{[}\mathsf{W}_{2}^{2}(\tilde{\mu}_{n},\nu)\big{]}\big{)}$ and then use the strict positivity of the densities and higher order regularity of optimal transport potentials to control the bias term as $\operatorname{\mathbb{E}}\big{[}\mathsf{W}_{2}^{2}(\tilde{\mu}_{n},\nu)\big{]}-\mathsf{W}_{2}^{2}(\mu,\nu)=o(n^{-1/2})$.

Finally, we note that our proof techniques differ from the aforementioned arguments for classic $\mathsf{W}_{p}$. Specifically, as opposed to the two-step approach of [MBNWW21] described above, we directly prove asymptotic normality for $\sqrt{n}\big{(}\mathsf{W}_{p}^{(\sigma)}(\hat{\mu}_{n},\nu)-\mathsf{W}_{p}^{(\sigma)}(\mu,\nu)\big{)}$ under the alternative. Their derivation does not apply to our case even when $p=2$ since their bias bound requires that the densities of $\mu$ and $\nu$ be bounded away from zero on their (compact) supports, which fails to hold after the Gaussian convolution. Our argument also differs from that of [SM18, TSM19], even though they also rely on the functional delta method. Specifically, since we do not assume that $\mu,\nu$ are discrete, $\mathsf{W}_{p}^{(\sigma)}$ can not be parameterized by simplex vectors, and hence the application of the functional delta method is nontrivial. Very recently, an independent work [HKSM22] used the extended functional delta method for the supremum functional [CCR20] to derive limit distributions for classic $\mathsf{W}_{p}$, with $p\geq 2$, for compactly supported distributions under the alternative in dimensions $d\leq 3$.

The rest of the paper is organized as follows. In Section 2, we collect background material on Wasserstein distances, smooth Wasserstein distances, and dual Sobolev spaces. In Section 3, we prove Theorem 1.1 and explore the validity of the bootstrap for empirical $\mathsf{W}_{p}^{(\sigma)}$. Section 4 presents applications of our limit distribution theory to MDE with $\mathsf{W}_{p}^{(\sigma)}$. Proofs for Section 3 and 4 can be found in Section 5. Section 6 provides concluding remarks and discusses future research directions. Finally, the Appendix contains additional proofs.

Let $|\cdot|$ and $\langle\cdot,\cdot\rangle$ denote the Euclidean norm and inner product, respectively. Let $B(x,r)=\{y\in\mathbb{R}^{d}:|y-x|\leq r\}$ denote the closed ball with center $x$ and radius $r$. We use $dx$ for the Lebesgue measure on $\mathbb{R}^{d}$, while $\delta_{x}$ denotes the Dirac measure at $x\in\mathbb{R}^{d}$. Given a finite signed Borel measure $\ell$ on $\mathbb{R}^{d}$, we identify $\ell$ with the linear functional $f\mapsto\ell(f):=\int fd\ell$. Let $\lesssim$ denote inequalities up to some numerical constants. For any $a,b\in\mathbb{R}$, we use the shorthands $a\vee b=\max\{a,b\}$ and $a\land b=\min\{a,b\}$.

For a topological space $S$, $\mathcal{B}(S)$ and $\mathcal{P}(S)$ denote, respectively, the Borel $\sigma$-field on $S$ and the class of Borel probability measures on $S$. We write $\mathcal{P}:=\mathcal{P}(\mathbb{R}^{d})$ and for $1\leq p<\infty$, use $\mathcal{P}_{p}$ to denote the subset of $\mu\in\mathcal{P}$ with finite $p$th moment $\int_{\mathbb{R}^{d}}|x|^{p}d\mu(x)<\infty$. For $\mu,\nu\in\mathcal{P}$, we write $\mu\ll\nu$ for the absolute continuity of $\mu$ w.r.t. $\nu$, and use $d\mu/d\nu$ for the corresponding Radon-Nikodym derivative. Let $\mu*\nu$ denote the convolution of $\mu,\nu\in\mathcal{P}$. Likewise, the convolution of two measurable functions $f,g:\mathbb{R}^{d}\to\mathbb{R}$ is denoted by $f*g$. We write $\gamma_{\sigma}=N(0,\sigma^{2}I_{d})$ for the centered Gaussian distribution on $\mathbb{R}^{d}$ with covariance matrix $\sigma^{2}I_{d}$, and use $\phi_{\sigma}(x)=(2\pi\sigma^{2})^{-d/2}e^{-|x|^{2}/(2\sigma^{2})}$ ($x\in\mathbb{R}^{d}$) for the corresponding density. We write $\mu\otimes\nu$ for the product measure of $\mu,\nu\in\mathcal{P}$. Let $\smash{\stackrel{{\scriptstyle w}}{{\to}},\stackrel{{\scriptstyle d}}{{\to}}}$, and $\smash{\stackrel{{\scriptstyle\operatorname{\mathbb{P}}}}{{\to}}}$ denote weak convergence of probability measures, convergence in distribution of random variables, and convergence in probability, respectively. When necessary, convergence in distribution is understood in the sense of Hoffmann-Jørgensen (cf. Chapter 2 in [vdVW96]).

Throughout the paper, we assume that $(X_{1},Y_{1}),(X_{2},Y_{2}),\dots$ are the coordinate projections of the product probability space $\prod_{i=1}^{\infty}\big{(}\mathbb{R}^{2d},\mathcal{B}(\mathbb{R}^{2d}),\mu\otimes\nu\big{)}$. To generate auxiliary random variables, we extend the probability space as $(\Omega,\mathcal{A},\operatorname{\mathbb{P}})=\left[\prod_{i=1}^{\infty}\big{(}\mathbb{R}^{2d},\mathcal{B}(\mathbb{R}^{2d}),\mu\otimes\nu\big{)}\right]\times\big{(}[0,1],\mathcal{B}([0,1]),\mathrm{Leb}\big{)}$, where $\mathrm{Leb}$ denotes the Lebesgue measure on $[0,1]$. For $\beta\in(0,2]$, let $\psi_{\beta}(t)=e^{t^{\beta}}-1$ for $t\geq 0$, and recall that the corresponding Orlicz (quasi-)norm of a real-valued random variable $\xi$ is defined as $\|\xi\|_{\psi_{\beta}}:=\inf\{C>0:\operatorname{\mathbb{E}}[\psi_{\beta}(|\xi|/C)]\leq 1\}$. A Borel probability measure $\mu\in\mathcal{P}$ is called $\beta$-sub-Weibull if $\||X|\|_{\psi_{\beta}}<\infty$ for $X\sim\mu$. We say that $\mu$ is sub-Weibull if it is $\beta$-sub-Weibull for some $\beta\in(0,2]$. Finally, $\mu$ is sub-Gaussian if it is $2$-sub-Weibull.

For an open set $\mathcal{O}$ in a Euclidean space, $C_{0}^{\infty}(\mathcal{O})$ denotes the space of compactly supported, infinitely differentiable, real functions on $\mathcal{O}$. We write $C_{0}^{\infty}=C_{0}^{\infty}(\mathbb{R}^{d})$ and define $\dot{C}_{0}^{\infty}=\{f+a:f\in C_{0}^{\infty},a\in\mathbb{R}\}$. For any $p\in[1,\infty)$ and $\mu\in\mathcal{P}(\mathbb{R}^{d})$, let $L^{p}(\mu;\mathbb{R}^{k})$ denote the space of measurable maps $f:\mathbb{R}^{d}\to\mathbb{R}^{k}$ such that $\|f\|_{L^{p}(\mu;\mathbb{R}^{k})}=(\int_{\mathbb{R}^{d}}|f|^{p}d\mu)^{1/p}<\infty$; when $d=1$ we use the shorthand $L^{p}(\mu)=L^{p}(\mu;\mathbb{R}^{1})$. Recall that $\big{(}L^{p}(\mu;\mathbb{R}^{k}),\|\cdot\|_{L^{p}(\mu;\mathbb{R}^{k})}\big{)}$ is a Banach space. Finally, for a subset $A$ of a topological space $S$, let $\overline{A}^{S}$ denote the closure of $A$; if the space $S$ is clear from the context, then we simply write $\overline{A}$ for the closure.

Recall that, for $1\leq p<\infty$, the $p$-Wasserstein distance $\mathsf{W}_{p}(\mu,\nu)$ between $\mu,\nu\in\mathcal{P}_{p}$ is defined in (1). Some basic properties of $\mathsf{W}_{p}$ are (cf. e.g., [Vil03, AGS08, Vil08, San15]): (i) the $\inf$ is attained in the definition of $\mathsf{W}_{p}$, i.e., there exists a coupling $\pi^{\star}\in\Pi(\mu,\nu)$ such that $\mathsf{W}_{p}^{p}(\mu,\nu)=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}|x-y|^{p}d\pi^{\star}(x,y)$, and the optimal coupling $\pi^{\star}$ is unique when $p>1$ and $\mu\ll dx$; (ii) $\mathsf{W}_{p}$ is a metric on $\mathcal{P}_{p}$; and (iii) convergence in $\mathsf{W}_{p}$ is equivalent to weak convergence plus convergence of $p$th moments: $\mathsf{W}_{p}(\mu_{n},\mu)\to 0$ if and only if $\mu_{n}\stackrel{{\scriptstyle w}}{{\to}}\mu$ and $\int|x|^{p}d\mu_{n}(x)\to\int|x|^{p}d\mu(x)$.

The proof of the limit distribution for empirical $\mathsf{W}_{p}^{(\sigma)}$ under the alternative hinges on duality theory for $\mathsf{W}_{p}$, which we summarize below. For a function $g:\mathbb{R}^{d}\to[-\infty,\infty)$ and a cost function $c:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$, the $c$-transform of $g$ is defined by

A function $g:\mathbb{R}^{d}\to[-\infty,\infty)$, not identically $-\infty$, is called $c$-concave if $g=f^{c}$ for some function $f:\mathbb{R}^{d}\to[-\infty,\infty)$.

The smooth Wasserstein distance convolves the distributions with an isotropic Gaussian kernel. Gaussian convolution levels out local irregularities in the distributions, while largely preserving the structure of classic $\mathsf{W}_{p}$. Recalling that $\gamma_{\sigma}=N(0,\sigma^{2}I_{d})$, the smooth $p$-Wasserstein distance is defined as follows.

The smooth Wasserstein distance was studied in [GGNP20, GG20, GGK20, SGK21, NGK21] for structural properties and empirical convergence rates. We recall two basic properties: (i) $\mathsf{W}_{p}^{(\sigma)}$ is a metric on $\mathcal{P}_{p}$ that generates the same topology as classic $\mathsf{W}_{p}$; (ii) for $\mu,\nu\in\mathcal{P}_{p}$ and $0\leq\sigma_{1}\leq\sigma_{2}<\infty$, we have $\mathsf{W}_{p}^{(\sigma_{2})}(\mu,\nu)\leq\mathsf{W}_{p}^{(\sigma_{1})}(\mu,\nu)\leq\mathsf{W}_{p}^{(\sigma_{2})}(\mu,\nu)+C_{p,d}\sqrt{\sigma_{2}^{2}-\sigma_{1}^{2}}$ for a constant $C_{p,d}$ that depends only on $p,d$. In particular, $\mathsf{W}_{p}^{(\sigma)}(\mu,\nu)$ is continuous and monotonically non-increasing in $\sigma\in[0,+\infty)$ with $\lim_{\sigma\downarrow 0}\mathsf{W}_{p}^{(\sigma)}(\mu,\nu)=\mathsf{W}_{p}(\mu,\nu)$. See [NGK21] for additional structural results, including an explicit expression for $C_{p,d}$ and weak convergence of smooth optimal couplings. For empirical convergence, it was shown in [NGK21] that under appropriate moment assumptions $\mathbb{E}\big{[}\mathsf{W}_{p}^{(\sigma)}(\hat{\mu}_{n},\mu)\big{]}=O(n^{-1/2})$ for $p>1$ in any dimension $d$. Versions of this result for $p=1$ and $p=2$ were derived earlier in [GGNP20, GGK20, SGK21].

Our proof strategy for the limit distribution results is to regard $\mathsf{W}_{p}$ as a functional defined on a subset of a certain dual Sobolev space. We will show that the smooth empirical process converges weakly in the dual Sobolev space and that $\mathsf{W}_{p}$ is Hadamard (directionally) differentiable w.r.t. the dual Sobolev norm. Given these, the limit distributions in Theorem 1.1 follow via the functional delta method. Here we briefly discuss (homogeneous) Sobolev spaces and their duals.

The restriction $f\in\dot{C}_{0}^{\infty}$ can be replaced with $f\in C_{0}^{\infty}$ in the definition of the dual norm $\|\cdot\|_{\dot{H}^{-1,p}(\rho)}$ since $\ell(f+a)=\ell(f)$ for any $\ell\in\dot{H}^{-1,p}(\rho)$.

We have defined the homogeneous Sobolev space $\dot{H}^{1,p}(\rho)$ as the completion of $\dot{C}_{0}^{\infty}$ w.r.t. $\|\cdot\|_{\dot{H}^{1,p}(\rho)}$. It is not immediately clear that the so-constructed space is a function space over $\mathbb{R}^{d}$. Below we present an explicit construction of $\dot{H}^{1,p}(\rho)$ when $d\rho/d\kappa$ is bounded away from zero for some reference measure $\kappa\gg dx$ satisfying the $p$-Poincaré inequality. To that end, we first define the Poincaré inequality.

The standard Poincaré inequality refers to the $2$-Poincaré inequality. It is known that any log-concave distribution (i.e., a distribution $\kappa$ of the form $d\kappa=e^{-V}dx$ for some convex function $V:\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$; cf. [LV07, SW14]) satisfies the $p$-Poincaré inequality for any $1\leq p<\infty$ [Bob99, Mil09]. In particular, the Gaussian distribution $\gamma_{\sigma}$ satisfies every $p$-Poincaré inequality (see also [Bog98, Corollary 1.7.3]).

The next lemma summarizes some basic results about the space $\dot{H}^{-1,p}(\rho)$ and $\dot{H}^{-1,p}(\rho)$-valued random variables that we use in the sequel. The proof can be found in Section A.1.

Consider a stochastic process $Y=(Y(f))_{f\in\dot{H}^{1,q}(\rho)}$ indexed by $\dot{H}^{1,q}(\rho)$, i.e., $\omega\mapsto Y(f,\omega)$ is measurable for each $f\in\dot{H}^{1,q}(\rho)$. The process can be thought of as a map from $\Omega$ into $\dot{H}^{-1,p}(\rho)$ as long as $Y$ has paths in $\dot{H}^{-1,p}(\rho)$, i.e., for each fixed $\omega\in\Omega$, the map $f\mapsto Y(f,\omega)$ is continuous and linear. The fact that the Borel $\sigma$-field on $\dot{H}^{-1,p}(\rho)$ coincides with the cylinder $\sigma$-field guarantees that a stochastic process indexed by $\dot{H}^{1,q}(\rho)$ with paths in $\dot{H}^{-1,p}(\rho)$ is Borel measurable as a map from $\Omega$ into $\dot{H}^{-1,p}(\rho)$.

In Section 3, we will explore limit distributions for empirical $\mathsf{W}_{p}^{(\sigma)}$. One of the key technical ingredients there is a comparison of the Wasserstein distance with a certain dual Sobolev norm, which we present next.

Proposition 2.1 follows directly from Theorem 5.26 of [DNS09]. Similar comparison inequalities appear in [Pey18, Led19, WB19b]. We include a self-contained proof of Proposition 2.1 in Section A.2 as some elements of the proof are key to our derivation of the null limit distribution for empirical $\mathsf{W}_{p}^{(\sigma)}$. The proof builds on the Benamou-Brenier dynamic formulation of optimal transport [BB00], which shows that $\mathsf{W}_{p}(\mu_{0},\mu_{1})$ is bounded from above by the length of any absolutely continuous path from $\mu_{0}$ to $\mu_{1}$ in $(\mathcal{P}_{p},\mathsf{W}_{p})$. The dual Sobolev norm emerges as a bound on the length of the linear interpolation $t\mapsto t\mu_{1}+(1-t)\mu_{0}$.

Throughout this section, we fix $1<p<\infty$, take $q$ as the conjugate index of $p$, and let $\sigma>0$. For $\mu,\nu\in\mathcal{P}_{p}$, let $X_{1},\ldots,X_{n}\sim\mu$ and $Y_{1},\ldots,Y_{n}\sim\nu$ be independent observations and denote the associated empirical distributions by $\hat{\mu}_{n}:=n^{-1}\sum_{i=1}^{n}\delta_{X_{i}}$ and $\hat{\nu}_{n}:=n^{-1}\sum_{i=1}^{n}\delta_{Y_{i}}$, respectively.