---
abstract: 'We propose a new statistical model, the *spiked transport model*, which formalizes the assumption that two probability distributions differ only on a low-dimensional subspace. We study the minimax rate of estimation for the Wasserstein distance under this model and show that this low-dimensional structure can be exploited to avoid the curse of dimensionality. As a byproduct of our minimax analysis, we establish a lower bound showing that, in the absence of such structure, the plug-in estimator is nearly rate-optimal for estimating the Wasserstein distance in high dimension. We also give evidence for a statistical-computational gap and conjecture that any computationally efficient estimator is bound to suffer from the curse of dimensionality.'
address:
- |
    [Jonathan Niles-Weed]{}\
    [Courant Institute of Mathematical Sciences]{}\
    [New York University]{}\
    [251 Mercer Street]{}\
    [New York, NY 10012-1185, USA]{}\
- |
    [Philippe Rigollet]{}\
    [Department of Mathematics]{}\
    [Massachusetts Institute of Technology]{}\
    [77 Massachusetts Avenue]{}\
    [Cambridge, MA 02139-4307, USA]{}\
author:
-  and 
bibliography:
- 'RigWee19ext.bib'
title: Estimation of Wasserstein distances in the Spiked Transport Model
---

Wasserstein distance; optimal transport; high-dimensional statistics

Introduction
============

Optimal transport is an increasingly useful toolbox in various data-driven disciplines such as machine learning [@AlaGraCut19; @PeyCut19; @SchHeiBon18; @FlaCutCou18; @ArjChiBot17; @GenPeyCut18; @JanCutGra19; @MonMulCut16; @RolCutPey16; @GraJouBer19; @StaClaSol17; @AlvJaaJeg18; @DelGamGor19; @CanRos12], computer graphics [@LavClaChi18; @SolGoePey15; @SolPeyKim16; @FeyChaVia17], statistics [@SegCut15; @AhiGouPar18; @RigWee18; @RigWee19; @WeeBer19; @ZemPan19; @PanZem19; @CazSegBig18; @RamTriCut17; @ClaSol18; @TamMun18; @KlaTamMun18; @BigCazPap17; @BarGorLes19; @KroSpoSuv19; @HutRig19; @Le-ParRig19] and the sciences [@SchShuTab19; @YanDamVen18; @delInoLou19; @LimLeeCha19]. A core primitive of this toolbox is the computation of Wasserstein distances between probability measures, and a natural statistical question is the estimation of Wasserstein distances from data.

A key object in this endeavor is the empirical measure $\mu_n$ associated to $\mu$. It is the empirical measure defined by $$\mu_n=\frac1n\sum_{i=1}^n\delta_{X_i}, \quad X_i \sim \mu \quad \text{i.i.d}\,.$$ Owing to their flexibility, Wasserstein distances are notoriously hard to estimate in high dimension since in such cases, the empirical distribution is a poor proxy for the underlying distribution of interest. Indeed, for $d$-dimensional distributions, Wasserstein distances between empirical measures generally converge at the slow rate $n^{-1/d}$ [@BoiLe-14; @FouGui15; @WeeBac18; @DerSchSch13; @Dud69] and thus suffer from the curse of dimensionality. For example, the following behavior is typical.

\[prop:empirical\] Let $\mu$ be a probability measure on $[-1, 1]^d$. If $\mu_n$ is an empirical measure comprising $n$ i.i.d. samples from $\mu$, then for any $p \in [1, \infty)$, $${\mathbb E}W_{p}\left(\mu_n, \mu\right) \leq {r_{p, d}(n)} {:=}c_{p} \sqrt d \left\{
\begin{array}{ll}
n^{-1/2p} & \text{if $d < 2p$} \\
n^{-1/2p} (\log n)^{1/p} & \text{if $d = 2p$} \\
n^{-1/d} & \text{if $d > 2p$}
\end{array}
\right.$$

In many settings, this bound is known to be tight up to logarithmic factors. In fact, the rate in \[prop:empirical\] has been shown to be essentially minimax optimal for the problem of estimating $\mu$ in Wasserstein distance [@SinPoc18], using any estimator, not necessarily the empirical distribution $\mu_n$.

Since $W_p$ satisfies the triangle inequality, \[prop:empirical\] readily yields that $W_p({\mu^{(1)}}, {\mu^{(2)}})$ for two probability measures ${\mu^{(1)}}, {\mu^{(2)}}$ on $[-1, 1]^d$ can be estimated at the rate $n^{-1/d}$ by the plug-in estimator $W_p({\mu^{(1)}}_n, {\mu^{(2)}}_n)$ when $d > 2p$, and a lower bound of the same order can be shown for the plug-in estimator when ${\mu^{(1)}}$ and ${\mu^{(2)}}$ are, for example, the uniform measure on $[-1, 1]^d$. However, while the above results give a strong indication that the Wasserstein distance $W_p({\mu^{(1)}}, {\mu^{(2)}})$ itself is also hard to estimate in high dimension, they do not preclude the existence of estimators that are better than $W_p({\mu^{(1)}}_n, {\mu^{(2)}}_n)$. Indeed, until recently, the best known lower bound for the problem of estimating the distance itself was of order $n^{-3/2d}$ [@Do-NguNgu11]. A concurrent and independent result [@Lia19] closes this gap for $p = 1$ and indicates that estimating the $W_1$ itself is essentially as hard as estimating the measure itself. Our results show that, in fact, estimating the distance $W_p({\mu^{(1)}}, {\mu^{(2)}})$ is essentially as hard as estimating a measure $\mu$ in $W_p$-distance, for any $p\ge 1$. As a result, any estimator of the distance itself must suffer from the curse of dimensionality.

One goal of *statistical optimal transport* is to develop new models and methods that overcome this curse of dimensionality by leveraging plausible structure in the problem. Early contributions in this direction include assuming smoothness [@HutRig19; @WeeBer19; @Lia19] or sparsity [@ForHutNit19]. In this work, we propose a new model, called the *spiked transport model*, to formalize the assumption that two distributions differ only on a low dimensional subspace of ${\mathbb{R}}^d$. Such an assumption forms the basis of several popular alternatives to Wasserstein distances such as the Sliced Wasserstein distance [@RabPeyDel11] or random one-dimensional projections [@PitKokDah07]. More recently, several numerical methods that exploit a form of low-dimensional structure were proposed together with illustration of their good numerical performance [@KolNadSim19; @DesHuSun19; @PatCut19].

To exploit the low-dimensional structure of the spiked transport model, we consider a standard method in statistics often called “projection pursuit" [@Kru69; @Kru72; @FriTuk74]. This general method aims to alleviate the challenges of high-dimensional data analysis by considering low-dimensional projections of the dataset that reveal “interesting” features of the data. We show that a suitable instantiation of this method to the present problem, which we call *Wasserstein Projection Pursuit* (WPP), leads to near optimal rates of estimation of the Wasserstein distance in the spiked transport model and permits to alleviate the curse of dimensionality from which the plug-in estimator suffers.

While our results establish a clear statistical picture, it is unclear how to implement WPP efficiently. An efficient relaxation of this estimator was recently proposed by @PatCut19, and a natural question is to analyze its performance in the spiked transport model. Instead of pursuing this direction we bring strong evidence that, in fact, *no* computationally efficient estimator is likely to be able to take advantage of the low-dimensional structure inherent to the spiked transport model. Our computational lower bounds come from the well-established statistical query framework [@Kea98]. In particular, they indicate a fundamental tradeoff between statistical and computational efficiency [@BerRig13; @BreBreHul18; @BanPerWei18; @MaWu15; @CaiLiaRak17]: computationally efficient methods to estimate Wasserstein distances, are bound to suffer the curse of dimensionality.

The rest of this paper is organized as follows. We introduce our model and main results in \[sec:model,sec:main\]. In \[sec:low-dim\], we define a low-dimensional version of the Wasserstein distance and establish its connection to the spiked transport model. proves the equivalence between transportation inequalities and subgaussian concentration properties of the Wasserstein distance. In \[sec:upper\] we propose and analyze an estimator for the Wasserstein distance under the spiked transport model. We establish a minimax lower bound in \[sec:lower\]. Finally, we prove a statistical query bound on the performance of efficient estimators in \[sec:comp\]. Supplementary proofs and lemmas appear in the appendices.

[Notation.]{}

We denote by $\|\cdot\|$ the Euclidean norm on ${\mathbb{R}}^d$. The symbols ${\left\|\cdot\right\|_{\mathrm{op}}}$ and ${\left\|\cdot\right\|_{\mathrm{F}}}$ denote the operator norm and Frobenius norm, respectively. If $X$ is a random variable on ${\mathbb{R}}$, we let $\|X\|_p {:=}({\mathbb E}|X|^p)^{1/p}$. Throughout, we use $c$ and $C$ to denote positive constants whose value may change from line to line, and we use subscripts to indicate when these constants depend on other parameters. We write $a \lesssim b$ if $a \leq C b$ holds for a universal positive constant $C$.

Model and methods {#sec:model}
=================

In this section, we describe the spiked transport model and Wasserstein projection pursuit.

Wasserstein distances
---------------------

Given two probability measures $\mu$ and $\nu$ on ${\mathbb{R}}^d$, let $\Gamma_{\mu,\nu}$ denote the set of *couplings* between $\mu$ and $\nu$ so that $\gamma \in \Gamma_{\mu,\nu}$ iff $\gamma(U\times {\mathbb{R}}^d)=\mu(U)$ and $\gamma({\mathbb{R}}^d\times V) = \nu(V)$.

For any $p \ge 1$, the *$p$-Wasserstein distance* $W_p$ between $\mu$ and $\nu$ is defined as $$\label{eq:wp}
W_p(\mu,\nu) :=\inf_{\gamma\in\Gamma_{\mu,\nu}} \left(\int_{{\mathbb{R}}^d\times {\mathbb{R}}^d} \|x-y\|^p\, {\mathrm{d}}\gamma(x,y)\right)^{1/p}\,.$$

The definition of Wasserstein distances may be extended to measures defined on general metric spaces but such extensions are beyond the scope of this paper. We refer the reader to @Vil09 for a comprehensive treatment.

Spiked transport model
----------------------

We introduce a new model that induces a low-dimensional structure on the optimal transport between two measures ${\mu^{(1)}}$ and ${\mu^{(2)}}$ over ${\mathbb{R}}^d$. To that end, fix a subspace ${\mathcal{U}}\subseteq {\mathbb{R}}^d$ of dimension $k \ll d$ and let $X^{(1)}, X^{(2)} \in {\mathcal{U}}$ be two random variables with arbitrary distributions. Next, let $Z$ be a third random variable, independent of $({X^{(1)}},{X^{(2)}})$ and such that $Z$ is supported on the orthogonal complement ${\mathcal{U}}^\perp$ of ${\mathcal{U}}$. Finally, let $$\label{eq:spiked_transport}
\begin{aligned}
{\mu^{(1)}}& {:=}{\mathrm{Law}({X^{(1)}}+ Z)} \\
{\mu^{(2)}}& {:=}{\mathrm{Law}({X^{(2)}}+ Z)}\,.
\end{aligned}$$ Though ${\mu^{(1)}}$ and ${\mu^{(2)}}$ are high-dimensional distributions, they differ only on the low-dimensional subspace ${\mathcal{U}}$. Borrowing terminology from principal component analysis [@Joh01], we say that the pair $({\mu^{(1)}},{\mu^{(2)}})$ satisfies the *spiked transport model* and we call ${\mathcal{U}}$ the *spike*.

We pose the following question: given $n$ independent observations from both ${\mu^{(1)}}$ and ${\mu^{(2)}}$, is it possible to estimate the Wasserstein distance between them at a rate faster than $n^{-1/d}$?

Concentration assumptions
-------------------------

In order to establish sharp statistical results for estimation of the Wasserstein distance, it is necessary to adopt smoothness and decay assumptions on the measures in question [see, e.g., @BobLed14; @FouGui15]. We focus on a family of such conditions known as *transport inequalities*, the study of which is a central object in the theory of concentration of measure [@Led01].

A probability measure $\mu$ on ${\mathbb{R}}^d$ is said to satisfy the ${{T_{p}(\sigma^2)}}$ transport inequality if $$\label{eq:tpineq}
W_p(\nu, \mu) \leq \sqrt{2 \sigma^2 D(\nu \| \mu)}$$ for all probability measures $\nu$ on ${\mathbb{R}}^d$.

These inequalities interpolate between several well known assumptions in high-dimensional probability. For example inequality ${{T_{1}(\sigma^2)}}$ is essentially equivalent to the assertion that $\mu$ is subgaussian, and $T_2(\sigma^2)$ is implied by (and often equivalent to) a stronger log-Sobolev inequality [@GozLeo10].

Our main results on the estimation of $W_p({\mu^{(1)}}, {\mu^{(2)}})$ are established under the assumption that both $\mu^{(1)}$ and $\mu^{(2)}$ satisfy a $T_p(\sigma^2)$ transport inequality. We show in Section \[sec:concentration\] that  is precisely equivalent to requiring that the random variable $W_p(\mu_n, \mu)$ is subgaussian.

Wasserstein projection pursuit
------------------------------

To take advantage of the spiked transport model, we employ a natural estimation method that we call Wasserstein Projection Pursuit (WPP).

Let $\mu$ and $\nu$ be two probability distributions on ${\mathbb{R}}^d$. Given a $k \times d$ matrix $U$ with orthonormal rows, let $\mu_U$ (resp. $\nu_U)$ denote the distribution of $U Y$ where $Y \sim \mu$ (resp. $Y\sim \nu$). We define $$\tilde W_{p,k}(\mu,\nu)  {:=}\max_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} W_p(\mu_U, \nu_U)\,,$$ where the maximization is taken over the *Stiefel manifold* ${{\mathcal{V}}_k({\mathbb{R}}^d)}$ of $k \times d$ matrices with orthonormal rows.

Given empirical measures ${\mu^{(1)}}_n$ and ${\mu^{(2)}}_n$ associated to ${\mu^{(1)}}$ and ${\mu^{(2)}}$ that satisfy the spiked transport model, we propose the following WPP estimator of $W_p({\mu^{(1)}}, {\mu^{(2)}})$: $${\hat W_{p,k}}=\tilde W_{p,k}({\mu^{(1)}}_n,{\mu^{(2)}}_n)\,,$$ In the next section, we show that this estimator is near-minimax-optimal.

Main results {#sec:main}
============

As a theoretical justification for Wasserstein projection pursuit, we prove that our procedure successfully avoids the curse of dimensionality under the spiked transport model. Our results primarily focus on the estimation of the the Wasserstein distance itself but we also obtain as a byproduct of Wasserstein projection pursuit an estimator for the spike ${\mathcal{U}}$ using standard perturbation results.

Estimation of the Wasserstein distance
--------------------------------------

The following theorem shows that Wasserstein projection pursuit takes advantage of the low-dimensional structure of the spiked transport model when estimating the Wasserstein distance.

\[thm:intro-upper\] Let $(\mu^{(1)},\mu^{(2)})$ satisfy the spiked transport model . For any $p \in [1,2]$, if $\mu^{(1)}$ and $\mu^{(2)}$ satisfy the ${{T_{p}(\sigma^2)}}$ transport inequality, then the WPP estimator ${\hat W_{p,k}}$ satisfies $${\mathbb E}\big|{\hat W_{p,k}}- W_p(\mu^{(1)}, \mu^{(2)})\big| \leq c_k \cdot \sigma \Big({r_{p, k}(n)} + \sqrt{\frac{d \log n}{n}}\Big)\,.$$

Strikingly, the rate ${r_{p, d}(n)}$ achieved by the naïve plug-in estimator (see \[prop:empirical\]) has been replaced by ${r_{p, k}(n)}$—in other words, this estimator enjoys the rate typical for $k$-dimensional rather than $d$-dimensional measures. The only dependence on the ambient dimension is in the second term, which is of lower order than the first whenever $p > 1$ or $k > 2$. A more general version of this theorem appears in \[sec:upper\].

Estimation of the spike
-----------------------

We show that if the distance between ${\mu^{(1)}}$ and ${\mu^{(2)}}$ is large enough, Wasserstein projection pursuit recovers the subspace ${\mathcal{U}}$. For simplicity, we state here the result when $k = 1$ and defer the full version to \[sec:upper\].

Let $(\mu^{(1)},\mu^{(2)})$ satisfy the spiked transport model with $k=1$ and let ${\mathcal{U}}$ be spanned by the unit vector $u \in {\mathbb{R}}^d$. Fix $p \in [1,2]$ and assume that $\mu^{(1)}$ and $\mu^{(2)}$ satisfy the ${{T_{p}(\sigma^2)}}$ transport inequality and that $W_p({\mu^{(1)}}, {\mu^{(2)}}) \gtrsim 1$. Then the estimator $$\hat u {:=}\operatorname*{argmax}_{v \in {\mathbb{R}}^d, \|v\| = 1}  W_p({\mu^{(1)}}_v, {\mu^{(2)}}_v)$$ satisfies $${\mathbb E}\sin^2\big(\measuredangle(\hat u, u)\big) \lesssim \sigma\cdot \Big(n^{-1/2p} + \sqrt{\frac{d \log n}{n}}\Big)\,.$$

Lower bounds
------------

To show that \[thm:intro-upper\] has the right dependence on $n$ and $d$, we exhibit two lower bounds, which imply that neither term in \[thm:intro-upper\] can be avoided.

To show the optimality of the first term, we define $${r'_{p, d}(n)} {:=}c_{p,d}  \left\{
\begin{array}{ll}
n^{-1/2p} & \text{if $d < 2p$} \\
n^{-1/2p}  & \text{if $d = 2p$} \\
(n \log n)^{-1/d} & \text{if $d > 2p$}
\end{array}
\right.$$

\[intro-lower\] Fix $p \geq 1$. For any estimator $\hat W$, there exists a pair of measures $\mu$ and $\nu$ supported on $[0, 1]^d$ such $${\mathbb E}|\hat W - W_p({\mu^{(1)}}, {\mu^{(2)}})| \geq {r'_{p, d}(n)}\,.$$

This lower bound readily implies that the plug-in estimator for the Wasserstein distance is optimal up to logarithmic factors. By embedding $[0, 1]^k$ into $[0, 1]^d$, this result likewise implies that the term ${r_{p, k}(n)}$ in \[thm:intro-upper\] is essentially optimal. A proof appears in \[sec:lower\]

Independently, @Lia19 recently obtained a similar result in the case $p=1$. More specifically, he proved that when $d \ge 2$, for any estimator $\hat W$, there exist probability measures ${\mu^{(1)}}$ and ${\mu^{(2)}}$ such that the following lower bound holds: $${\mathbb E}|\hat W - W_1({\mu^{(1)}}, {\mu^{(2)}})| \gtrsim \frac{\log \log n}{\log n} n^{-1/d}$$ In particular, while our lower bound is slightly stronger and holds for all $p \geq 1$, both our result and that of @Lia19 fail to match the naïve upper bound of order $n^{1/d}$ by logarithmic factors when $d$ is large. The presence of a logarithmic factor in our lower bound comes from a reduction to estimating the total variation distance. In that case, as in several other instances of functional estimation problems, the presence of this factor is, in fact, optimal, and has been dubbed *sample size enlargement* [@JiaHanWei18]. Closing this gap in the context of estimation of the Wasserstein distance is an interesting and fundamental question.

The only appearance of the ambient dimension $d$ is in the second term of \[thm:intro-upper\]. The following theorem shows that this dependence cannot be eliminated, even when $k = 1$.

\[deviation\_lb\] Let $p \in [1, 2]$ and $\sigma > 0$, and assume $k=1$. For all estimators $\hat W$, there exists a pair of measures ${\mu^{(1)}}$ and ${\mu^{(2)}}$ satisfying the spiked transport model and $T_p(\sigma^2)$ such that $${\mathbb E}|\hat W - W_p({\mu^{(1)}}, {\mu^{(2)}})| \gtrsim \sigma \sqrt{\frac{d} n}\,.$$

The proof of \[deviation\_lb\] is deferred to the appendix. For problems where $d$ is large, dependence on $d$ may be a crippling limitation. In that case, we conjecture that assuming a sparse spike, in the same spirit as sparse PCA, can mitigate this effect and bring interpretability to the estimated spike.

A computational-statistical gap
-------------------------------

The WPP estimator achieving the rate in \[thm:intro-upper\] is computationally expensive to implement, which raises the question of whether an efficient estimator exists achieving the same rate. We give evidence in the form of a statistical query lower bound that no such estimator exists. The statistical query model considers algorithms with access to an oracle ${\mathrm{VSTAT}}(t)$, where $t > 0$ is a parameter which plays the role of sample size. We show that any such algorithm for estimating the Wasserstein distance needs an *exponential* number of queries to an oracle with *exponential* sample size parameter, even under the spiked transport model. By contrast, \[thm:intro-upper\] implies that a non-efficient estimator needs a number of samples only polynomial in the dimension.

Let $p \in [1, 2]$, and consider probability measures ${\mu^{(1)}}$ and ${\mu^{(2)}}$ satisfying the spiked transport model. There exists a positive constant $c$ such that any statistical query algorithm which estimates $W_p({\mu^{(1)}}, {\mu^{(2)}})$ to accuracy $1/\mathrm{poly}(d)$ with probability at least $2/3$ requires at least $2^{cd}$ queries to ${\mathrm{VSTAT}}(2^{cd})$.

Low-dimensional Wasserstein distances {#sec:low-dim}
=====================================

Motivated by projection pursuit, we define the following version of the Wasserstein distance which measures the discrepancy between low-dimensional projections of the measures.

For $k \in [d]$, the *$k$-dimensional Wasserstein distance* between ${\mu^{(1)}}$ and ${\mu^{(2)}}$ is $${\tilde W_{p,k}}({\mu^{(1)}}, {\mu^{(2)}}) {:=}\sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} W_p({\mu^{(1)}}_U, {\mu^{(2)}}_U)\,.$$

This definition has been proposed independently and concurrently by a number of other recent works [@KolNadSim19; @DesHuSun19; @PatCut19].

We will use throughout the following basic fact about the $k$-dimensional Wasserstein distance [@PatCut19 Proposition 1].

\[prop:triangle\] ${\tilde W_{p,k}}$ is a metric on the set of probability measures over ${\mathbb{R}}^d$ with finite $p$th moment.

The definition of ${\tilde W_{p,k}}$ is chosen so that, under the spiked transport model , the $k$-dimensional Wasserstein distance agrees with the normal Wasserstein distance.

\[prop:wpk-is-wp\] Under the spiked transport model , $${\tilde W_{p,k}}({\mu^{(1)}}, {\mu^{(2)}}) = W_p({\mu^{(1)}}, {\mu^{(2)}})\,.$$

follows from the following statement, which pertains to distributions that are allowed to have a different component on the space orthogonal to the spike ${\mathcal{U}}$. Suppose ${\nu^{(1)}}$ and ${\nu^{(2)}}$ satisfy $$\label{eq:almost-spike}
\begin{aligned}
{\nu^{(1)}}& = {\mathrm{Law}(X^{(1)} + Z^{(1)})} \\
{\nu^{(2)}}& = {\mathrm{Law}(X^{(2)} + Z^{(2)})}\,,
\end{aligned}$$ where as before $X^{(1)}$ and $X^{(2)}$ are supported on a subspace ${\mathcal{U}}$ and $Z^{(1)}$ and $Z^{(2)}$ are supported on its orthogonal complement ${\mathcal{U}}^\perp$, and where we assume that $X^{(i)}$ and $Z^{(i)}$ are independent for $i \in \{1, 2\}$. Note that unlike in the spiked transport model, the components $Z^{(1)}$ and $Z^{(2)}$ on the orthogonal complement of ${\mathcal{U}}$ need not be identical. The following result shows that under this relaxed model, the $k$-dimensional Wasserstein distance between ${\nu^{(1)}}$ and ${\nu^{(2)}}$ still captures the true Wasserstein distance between the distributions as long as the distributions of $Z^{(1)}$ and $Z^{(2)}$ are sufficiently close.

\[prop:almost-spike\] Under the relaxed spiked transport model , $$|{\tilde W_{p,k}}({\nu^{(1)}}, {\nu^{(2)}}) - W_p({\nu^{(1)}}, {\nu^{(2)}})| \leq W_p({\mathrm{Law}(Z^{(1)})}, {\mathrm{Law}(Z^{(2)})})\,.$$

Concentration {#sec:concentration}
=============

A key step to establish the upper bound of \[sec:upper\] consists in establishing good concentration properties for the Wasserstein distance between a measure and its empirical counterpart. The main assumption we adopt is that the measures in question satisfy a *transport inequality*. Since the pioneering work of [@Mar96; @Mar96a] and [@Tal96], transport inequalities have played a central role in the analysis of the concentration properties of high-dimensional measures.

We require two definitions.

Given a Polish space  ${\mathcal{X}}$ equipped with a metric $\rho$, denote by ${\mathcal{P}}({\mathcal{X}})$ the space of all Borel probability measures ${\mathcal{X}}$. Let ${\mathcal{P}}_p({\mathcal{X}}) {:=}\{\mu \in {\mathcal{P}}({\mathcal{X}}) : \int \rho(x, \cdot)^p {\, \mathrm{d}}\mu(x) < \infty\}$.

A measure $\mu \in {\mathcal{P}}_p({\mathcal{X}})$ satisfies the *${{T_{p}(\sigma^2)}}$ inequality* for some $\sigma > 0$ if $$W_p(\nu, \mu) \leq \sqrt{2 \sigma^2 D(\nu \| \mu)} \quad \quad \forall \nu \in {\mathcal{P}}_p({\mathcal{X}})\,,$$ where $W_p$ is the Wasserstein-$p$ distance on $({\mathcal{X}},\rho)$ and $D$ is the Kullback-Leibler divergence.

\[def:subg\] A random variable $X$ on ${\mathbb{R}}$ is $\sigma^2$-subgaussian if $${\mathbb E}e^{\lambda(X - {\mathbb E}X)} \leq e^{\lambda^2 \sigma^2/2} \qquad \forall \lambda \in {\mathbb{R}}\,.$$

In this section, we present a surprisingly simple equivalence between transport inequalities and subgaussian concentration for the Wasserstein distance. The essence of this result is present in the works of @GozLeo10 [see @GozLeo10; @Goz09; @GozLeo07], and similar bounds have been obtained by @BolGuiVil07. Nevertheless, we could not find this simple fact stated in a form suitable for our purposes in the literature. For any measure $\mu$, recall that the random measure $\mu_n {:=}\frac 1n \sum_{i=1}^n \delta_{X_i}$, where $X_i \sim \mu$ i.i.d., denotes its associated empirical measure.

\[thm:concentration\_is\_tp\] Let $p \in [1,2]$. A measure $\mu \in {\mathcal{P}}_p({\mathcal{X}})$ satisfies ${{T_{p}(\sigma^2)}}$ if and only if the random variable $W_p(\mu_n, \mu)$ is $\sigma^2/n$-subgaussian for all $n$.

Because ${{T_{p'}(\sigma^2)}}$ implies ${{T_{p}(\sigma^2)}}$ when $p' \geq p$, \[thm:concentration\_is\_tp\] also implies that a measure satisfying ${{T_{p'}(\sigma^2)}}$ has good concentration for $W_p$ if $p \leq p'$. In the opposite direction, if $p > p'$, then satisfying ${{T_{p'}(\sigma^2)}}$ still yields a weaker concentration bound. A modification of the proof of \[thm:concentration\_is\_tp\] yields the following result.

\[thm:p-pprime-concentration\] Let $p' \in [1, 2]$ and $p \geq 1$. If $\mu \in {\mathcal{P}}_{p'}({\mathcal{X}})$ satisfies ${{T_{p'}(\sigma^2)}}$, then $W_p(\mu_n, \mu)$ is $\sigma^2/n^{1 - \left(\frac{2}{p'} - \frac 2 p\right)_+}$ subgaussian.

The conclusion of \[thm:p-pprime-concentration\] is interesting whenever $ \left(\frac{2}{p'} - \frac 2 p\right)_+ < 1$. For example, if we assume merely that $\mu$ satisfies ${{T_{1}(\sigma^2)}}$, \[thm:p-pprime-concentration\] only yields a nontrivial concentration result for $W_p(\mu_n, \mu)$ when $p < 2$; by contrast, if $\mu$ satisfies $T_2(\sigma^2)$, then \[thm:p-pprime-concentration\] implies a concentration result for $W_p(\mu_n, \mu)$ for all $p < \infty$.

In \[sec:upper\], we require concentration properties not of the Wasserstein distance itself but of the $k$-dimensional Wasserstein distance. The following result shows that low-dimensional projections inherit the concentration properties of the $d$-dimensional measure.

Let $U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}$. For any $p \in [1, 2]$ and $\sigma > 0$, if $\mu$ satisfies ${{T_{p}(\sigma^2)}}$, then so does $\mu_U$.

The projection $x \mapsto Ux$ is a contraction. The result then follows from @Goz07 [Corollary 20] [see also @Mau91].

We conclude this section by giving some simple conditions under which the ${{T_{1}(\sigma^2)}}$ inequality is satisfied. The following characterization is well known. Denote by ${\mathrm{Lip}}({\mathcal{X}})$ the space of all functions $f: {\mathcal{X}}\to {\mathbb{R}}$ satisfying $|f(x) - f(y)| \leq d(x, y)$ for all $x, y \in {\mathbb{R}}$.

\[prop:bobgot\] A measure $\mu \in {\mathcal{P}}_1({\mathcal{X}})$ satisfies ${{T_{1}(\sigma^2)}}$ if and only if $f(X)$ is $\sigma^2$-subgaussian for all $f \in {\mathrm{Lip}}({\mathcal{X}})$.

It is common to extend \[def:subg\] to random vectors as follows.

A random variable $X$ on ${\mathbb{R}}^d$ is *$\sigma^2$-subgaussian* if $u^\top X$ is $\sigma^2$-subgaussian for all $u \in {\mathbb{R}}^d$ satisfying $\|u\| = 1$.

Subgaussian random vectors yield a large collection of random variables satisfying a $T_1$ inequality.

\[lem:t1\_is\_subg\] If $\mu$ on ${\mathbb{R}}^k$ satisfies ${{T_{p}(\sigma^2)}}$ for any $p \geq 1$, then $X \sim \mu$ is $\sigma^2$-subgaussian. Conversely, if $X \sim \mu$ on ${\mathbb{R}}^k$ is $\sigma^2$-subgaussian, then $\mu$ satisfies ${{T_{1}(C k \sigma^2)}}$ for a universal constant $C > 0$.

If the entries of $X$ are independent, then the result holds with $C = 1$ by a result of @Mar96. The presence of the factor $k$ in the converse statement is unavoidable; unlike $T_2$ inequalities, $T_1$ inequalities do not exhibit dimension-free concentration [@GozLeo10].

Upper bounds {#sec:upper}
============

In this section, we establish that under the spiked transport model, Wasserstein projection pursuit produces a significantly more accurate estimate of the Wasserstein distance than the plug-in estimator.

Let ${\mu^{(1)}}$ and ${\mu^{(2)}}$ be two measures generated according to the spiked transport model . For $i \in \{1, 2\}$, we let ${\mu^{(i)}}_n {:=}\frac 1n \sum_{j=1}^n \delta_{X^{(i)}_j}$, where $X^{(i)}_j \sim {\mu^{(i)}}$ are i.i.d. We define $$\hat W_{p,k} {:=}{\tilde W_{p,k}}({\mu^{(1)}}_n, {\mu^{(2)}}_n)\,.$$

Our main upper bound shows that $\hat W_{p,k}$ converges to the true Wasserstein distance $W_p(\mu, \nu)$ at a rate much faster than $n^{-1/d}$.

\[thm:main\] Let $p' \in [1, 2]$ and $p \geq 1$. Under the spiked transport model, if ${\mu^{(1)}}$ and ${\mu^{(2)}}$ satisfy ${{T_{p'}(\sigma^2)}}$, then $${\mathbb E}|\hat W_{p,k} - W_p({\mu^{(1)}}, {\mu^{(2)}})| \lesssim \sigma \left({r_{p, k}(n)} + c_p \cdot n^{(\frac{1}{p'} - \frac{1}{p})_+} \sqrt\frac{dk \log n}{n}\right)\,.$$

can also be extended to the misspecified model proposed in 

\[thm:almost\] Let $p' \in [1, 2]$ and $p \geq 1$. Under the relaxed spiked transport model , if ${\nu^{(1)}}$ and ${\nu^{(2)}}$ satisfy ${{T_{p'}(\sigma^2)}}$, then $${\mathbb E}|\hat W_{p,k} - W_p({\nu^{(1)}}, {\nu^{(2)}})| \lesssim \sigma \left({r_{p, k}(n)} + c_p \cdot n^{(\frac{1}{p'} - \frac{1}{p})_+} \sqrt\frac{dk \log n}{n}\right) + {\varepsilon}\,,$$ where ${\varepsilon}= W_p({\mathrm{Law}(Z^{(1)})}, {\mathrm{Law}(Z^{(2)})})$.

, which follows almost immediately from the proof of \[thm:main\], establishes that Wasserstein projection pursuit brings statistical benefits even in the situation where the spiked transport model holds only approximately.

follows from the following two propositions. We first show that the quality of the proposed estimator $\hat W_{p,k}$ can be bounded by the sum of two terms depending only on ${\mu^{(1)}}_n$ and ${\mu^{(2)}}_n$ individually.

\[prop:split\_obj\] $${\mathbb E}| \hat W_p - W_p({\mu^{(1)}}, {\mu^{(2)}})| \leq {\mathbb E}{\tilde W_{p,k}}({\mu^{(1)}},  {\mu^{(1)}}_n) + {\mathbb E}{\tilde W_{p,k}}({\mu^{(2)}}, {\mu^{(2)}}_n)$$

Since $\hat W_p = {\tilde W_{p,k}}({\mu^{(1)}}_n, {\mu^{(2)}}_n)$ and $W_p({\mu^{(1)}}, {\mu^{(2)}}) = {\tilde W_{p,k}}({\mu^{(1)}}, {\mu^{(2)}})$, the claim is immediate from \[prop:triangle\].

The following proposition allows us to bound both terms of \[prop:split\_obj\] by the desired quantity.

\[prop:supremum\] Let $p' \in [1, 2]$ and $p \geq 1$. If $\mu$ satisfies ${{T_{p'}(\sigma^2)}}$, then $${\mathbb E}{\tilde W_{p,k}}(\mu, \mu_n) \lesssim \sigma \left(r_{p,k}(n) + c
_p \cdot n^{(\frac{1}{p'} - \frac{1}{p})_+}\sqrt{\frac{dk \log n}{n}}\right)\,.$$

Wasserstein distances are invariant under translating both measures by the same vector. Therefore, we can assume without loss of generality that $\mu$ has mean $0$. Likewise, by homogeneity, we assume $\sigma = 1$.

Let $Z_U {:=}W_p(\mu_U, (\mu_n)_U)$. We first show that the process $Z_U$ is Lipschitz.

There exists a random variable $L$ such that for all $U, V \in {\mathcal{V}}_k({\mathbb{R}}^d)$, $$|Z_U - Z_V| \leq L {\left\|U - V\right\|_{\mathrm{op}}}$$ and ${\mathbb E}L \lesssim \sqrt{dp}$.

Let $X \sim \mu$. Then $$\begin{aligned}
|Z_U - Z_V| & \leq W_{p}(\mu_U, \mu_V) + W_{p}((\mu_n)_U, (\mu_n)_V) \\
& \leq ({\mathbb E}\|(U-V) X\|^p)^{1/p} + \Big(\frac 1 n\sum_{i=1}^n \|(U-V) X_i\|^p\Big)^{1/p} \\
& \leq {\left\|U-V\right\|_{\mathrm{op}}} \Big(({\mathbb E}\|X\|^p)^{1/p} + \Big(\frac 1 n\sum_{i=1}^n \|X_i\|^p\Big)^{1/p}\Big)\,.\end{aligned}$$

We obtain that $$|Z_U - Z_V| \leq L {\left\|U-V\right\|_{\mathrm{op}}}$$ where $L = ({\mathbb E}\|X\|^p)^{1/p} + \left(\frac 1 n\sum_{i=1}^n \|X_i\|^p\right)^{1/p}$. By Jensen’s inequality, we have ${\mathbb E}L \leq 2 ({\mathbb E}\|X\|^p)^{1/p}$. Together with \[prop:subg\_vec\_norm\], it yields the claim.

By \[thm:p-pprime-concentration\], for all $U \in {\mathcal{V}}_k({\mathbb{R}}^d)$, the random variable $Z_U$ is $n^{-1 + \left(\frac{2}{p'} - \frac 2 p\right)_+}$ subgaussian. Therefore, by a standard ${\varepsilon}$-net argument, if we denote by ${\mathcal{N}}({\mathcal{V}}_k, {\varepsilon}, {\left\|\cdot\right\|_{\mathrm{op}}})$ the covering number of ${\mathcal{V}}_k$ with respect to the operator norm, we obtain $${\mathbb E}\sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} (Z_U - {\mathbb E}Z_U) \lesssim \inf_{{\varepsilon}> 0} \left\{{\varepsilon}{\mathbb E}L+  n^{(\frac 1p - \frac{1}{p'})_+}\sqrt{\frac{\log {\mathcal{N}}({\mathcal{V}}_k, {\varepsilon}, {\left\|\cdot\right\|_{\mathrm{op}}})}{n}}\right\}\,.$$ shows that there exists a universal constant $c$ such that $\log {\mathcal{N}}({\mathcal{V}}_k, {\varepsilon}, {\left\|\cdot\right\|_{\mathrm{op}}}) \leq dk \log \frac {c\sqrt k} {\varepsilon}$ for ${\varepsilon}\in (0, 1]$. Choosing ${\varepsilon}= \sqrt{k/n}$ yields $$\begin{aligned}
{\mathbb E}\sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} (Z_U - {\mathbb E}Z_U) & \lesssim \sqrt{\frac{dkp}{n}} + n^{(\frac{1}{p'} - \frac{1}{p})_+}\sqrt{\frac{dk \log n}{n}} \\
& \leq c_{p} \cdot  n^{(\frac{1}{p'} - \frac{1}{p})_+}\sqrt{\frac{dk \log n}{n}}\,.\end{aligned}$$

Applying \[subgaussian\_wp\_rate\] yields $$\begin{aligned}
{\mathbb E}\sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} W_p(\mu_U, (\mu_n)_U) & \leq  \sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} {\mathbb E}W_p(\mu_U, (\mu_n)_U)+{\mathbb E}\sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} (Z_U - {\mathbb E}Z_U) \\
& \lesssim r_{p,k}(n) + c_p\cdot n^{(\frac{1}{p'} - \frac{1}{p})_+} \sqrt{\frac{dk \log n}{n}}\,,\end{aligned}$$ as claimed.

We also obtain a Davis-Kahan-type theorem on subspace recovery. Given two subspaces ${\mathcal{U}}_1$ and ${\mathcal{U}}_2$, the *minimal angle* [@Deu95; @Dix49] between them is defined to be $$\measuredangle({\mathcal{U}}_1, {\mathcal{U}}_2) {:=}\arccos \left(\sup_{u_1 \in {\mathcal{U}}_1, u_2 \in {\mathcal{U}}_2} \frac{u_1^\top u_2}{\|u_1\| \|u_2\|}\right)\,.$$ If $\measuredangle({\mathcal{U}}_1, {\mathcal{U}}_2) = 0$, then ${\mathcal{U}}_1 \cap {\mathcal{U}}_2 \neq \{0\}$, so that ${\mathcal{U}}_1$ and ${\mathcal{U}}_2$ are at least partially aligned. In the important special case that ${\mathcal{U}}_1$ and ${\mathcal{U}}_2$ are each one dimensional, this definition reduces to the angle between the subspaces.

The following result indicates that as long as ${\mu^{(1)}}$ and ${\mu^{(2)}}$ are well separated, Wasserstein projection pursuit also yields a subspace with at least partial alignment to ${\mathcal{U}}$.

\[thm:sin\] Let $p' \in [1, 2]$ and $p \geq 1$. Assume that ${\mu^{(1)}}$ and ${\mu^{(2)}}$ satisfy the spiked transport model and ${{T_{p'}(\sigma^2)}}$ Let $\hat {\mathcal{U}}{:=}\mathrm{span}(\hat U)$, where $$\hat U {:=}\operatorname*{argmax}_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} W_p\big(({\mu}_n^{(1)})_U, ({\mu}_n^{(2)})_U\big)\,.$$ Then $${\mathbb E}\sin^2\big(\measuredangle(\hat {\mathcal{U}}, {\mathcal{U}})\big) \lesssim \frac{\sigma \left({r_{p, k}(n)} + c_p \cdot n^{(\frac{1}{p'} - \frac{1}{p})_+} \sqrt\frac{dk \log n}{n}\right)}{W_p({\mu^{(1)}}, {\mu^{(2)}})}\,.$$

A proof of \[thm:sin\] appears in the appendix. Note that the a bound on the minimal angle is a rather weak guarantee. Indeed, $\measuredangle(\hat {\mathcal{U}}, {\mathcal{V}})\to 0$ implies that the subspaces $\hat {\mathcal{U}}$ and ${\mathcal{U}}$ share at least a common line asymptotically but not more. When $k=1$, this ensures recovery of the subspace, but this no longer holds true for higher dimensional spikes. In retrospect such a guarantee is all that we can hope for under the mere assumption that $W_p({\mu^{(1)}}, {\mu^{(2)}})>0$. Indeed, it may be the case that these distributions differ only on a one dimensional space. Stronger guarantees may be achieved by assuming that $W_p({\mu^{(1)}}_{V}, {\mu^{(2)}}_{V})>0$ for a large family of $V$, but we leave them for future research.

A lower bound on estimating the Wasserstein distance {#sec:lower}
====================================================

In this section, we prove that the rate ${r_{p, d}(n)}$ is optimal for estimating the Wasserstein distance, up to logarithmic factors. The core idea of our lower bound is to relate estimating the Wasserstein distance to the problem of estimating total variation distance, sharp rates for which are known [@ValVal11; @JiaHanWei18]. To obtain sufficient control over the Wasserstein distance as a function of total variation, we prove a refined bound incorporating both total variation and the $\chi^2$ divergence (\[random\_injection\]). We then show a modified lower bound (\[tv\_lb\]) for a testing problem involving the total variation distance over the class of distributions on $[m]$ close to the uniform measure in $\chi^2$ divergence. In the interest of generality, we formulate our results for any compact metric space ${\mathcal{X}}$ whose covering numbers satisfy $$\label{covering-assumption}
c {\varepsilon}^{-d} \leq {\mathcal{N}}({\mathcal{X}}, {\varepsilon}) \leq C {\varepsilon}^{-d}$$ for all ${\varepsilon}\leq {\mathrm{diam}}({\mathcal{X}})$. This condition clearly holds for compact subsets of ${\mathbb{R}}^d$ and more generally for metric spaces with Minkowski dimension $d$. We adopt the assumption ${\mathrm{diam}}({\mathcal{X}}) = 1$ without loss of generality.

Let ${\mathcal{P}}$ be the set of distributions supported on ${\mathcal{X}}$ and let $R(n, {\mathcal{P}})$ denote the minimax risk over ${\mathcal{P}}$, $$R(n, {\mathcal{P}}) {:=}\inf_{\hat W} \sup_{\mu, \nu \in {\mathcal{P}}} {\mathbb E}_{\mu, \nu} |\hat W - W_p(\mu, \nu)|\,.$$

The bound $R(n, {\mathcal{P}}) \gtrsim n^{-1/2p}$ is an almost trivial consequence of the fact that the distribution $\frac 12 \delta_{-1} + \frac 12 \delta_1$ cannot be distinguished from $(\frac 12 + {\varepsilon}) \delta_{-1} + (\frac 12 - {\varepsilon}) \delta_1$ on the basis of $n$ samples when ${\varepsilon}\asymp n^{-1/2}$. The interesting part of \[intro-lower\] is the rate when $d > 2p$. We prove the following.

\[dist-est-lb-general\] Let $d > 2p \geq 2$ and assume ${\mathcal{X}}$ satisfies . Then $$R(n, {\mathcal{P}}) \geq C_{d, p} (n \log n)^{-1/d}\,.$$

Before proving \[dist-est-lb-general\], we establish the two propositions described above. allows us to reduce \[dist-est-lb-general\] to an estimation problem involving total variation distance, and \[tv\_lb\] is a lower bound on the minimax rate for that total variation estimation problem.

\[random\_injection\] Assume $d > 2p \geq 2$, and let $m$ be a positive integer. Let $u$ be the uniform distribution on $[m]:=\{1, \ldots, m\}$. There exists a random function $F: [m] \to X$ such that for any distribution $q$ on $[m]$, $$c m^{-1/d} {d_{\mathrm{TV}}(q, u)}^{\frac 1 p} \leq W_p(F_\sharp q, F_\sharp u) \leq C_{d, p} m^{-1/d} (\chi^2(q, u))^{1/d} {d_{\mathrm{TV}}(q, u)}^{\frac 1 p - \frac 2 d}$$ with probability at least $.9$.

shows that the condition ${\mathcal{N}}({\mathcal{X}}, {\varepsilon}) \geq c {\varepsilon}^{-d}$ implies the existence a set ${\mathcal{G}}_m {:=}\{x_1, \dots, x_m\} \subseteq {\mathcal{X}}$ such that $d(x_i, x_j) \gtrsim m^{-1/d}$ for all $i \neq j$. We select $F$ uniformly at random from the set of all bijections from $[m]$ to ${\mathcal{G}}_m$.

To show the lower bound, we note that any points $x, y \in {\mathcal{G}}_m$ satisfy $$d(x, y)^p \gtrsim m^{-p/d} {\mathds{1}}\{x \neq y\}\,,$$ which implies that for any coupling $\pi$ between $F_\sharp q$ and $F_\sharp u$ $$\begin{aligned}
\int d(x, y)^p {\, \mathrm{d}}\pi(x, y)& \gtrsim m^{-p/d} {\mathbb P}_\pi[X \neq Y]\\
& \geq m^{-p/d} {d_{\mathrm{TV}}(F_\sharp q, F_\sharp u)} \\
&= m^{-p/d}{d_{\mathrm{TV}}(q, u)}\,.\end{aligned}$$ The lower bound therefore holds with probability $1$.

We now turn to the upper bound. We employ a dyadic covering bound [@WeeBac18 Proposition 1]. For any $k^*$ there exists a dyadic partition $\{{\mathcal{Q}}^k\}_{1 \leq k \leq k^*}$ of $X$ with parameter $\delta = 1/3$ such that $|{\mathcal{Q}}^k| \leq {\mathcal{N}}({\mathcal{X}}, 3^{-(k+1)})$. We obtain that for any $k^* \geq 0$, $$W_p^p(F_\sharp q, F_\sharp u) \leq 3^{-k^* p} + \sum_{k=1}^{k^*} 3^{-(k-1)p} \sum_{Q_i^k \in {\mathcal{Q}}^k} |F_\sharp q(Q_i^k) - F_\sharp u(Q_i^k)|\,,$$

By \[random\_partition\], for any $k$, $${\mathbb E}\sum_{Q_i^k \in {\mathcal{Q}}^k} |F_\sharp q(Q_i^k) - F_\sharp u(Q_i^k)| \leq 2{d_{\mathrm{TV}}(q, u)} \wedge C_{d, p}\left(\frac{3^{kd}\chi^2(q, u)}{m}\right)^{1/2} \,.$$ Let $k_0$ be a positive integer to be fixed later.

By applying the first bound, we obtain $${\mathbb E}\sum_{k > k_0} 3^{-(k-1)p} \sum_{Q_i^k \in {\mathcal{Q}}^k} |F_\sharp q(Q_i^k) - F_\sharp u(Q_i^k)| \lesssim 3^{-k_0 p} {d_{\mathrm{TV}}(q, u)}\,.$$ Applying the second bound and recalling that $d/2 > p$ yields $$\begin{aligned}
{\mathbb E}\sum_{k \leq k_0} 3^{-(k-1)p} \sum_{Q_i^k \in {\mathcal{Q}}^k} |F_\sharp q(Q_i^k) - F_\sharp u(Q_i^k)| & \leq C_{d, p} \left(\frac{\chi^2(q, u)}{m}\right)^{1/2} \sum_{k \leq k_0} 3^{k (d/2 - p)} \\
&  \leq C_{d, p} \left(\frac{\chi^2(q, u)}{m}\right)^{1/2} 3^{k_0(d/2 - p)}\,,\end{aligned}$$

We obtain for any $k^* \geq 0$ that $${\mathbb E}W_p^p(F_\sharp q, F_\sharp u) \leq C_{d, p} \left(\frac{\chi^2(q, u)}{m}\right)^{1/2} 3^{k_0(d/2 - p)} + C \cdot 3^{-k_0 p} {d_{\mathrm{TV}}(q, u)} + 3^{-k^* p}\,,$$ and taking $k^* \to \infty$ it suffices to bound the first two terms.

Let $k_0$ to be the smallest positive integer such that $$3^{k_0 d} \geq m \frac{{d_{\mathrm{TV}}(q, u)}^{2}}{\chi^{2}(q, u)}\,.$$ Then $$3^{k_0 d/2} \left(\frac{\chi^2(q, u)}{m}\right)^{1/2} \leq C_d \cdot {d_{\mathrm{TV}}(q, u)}\,,$$ and hence $${\mathbb E}W_p^p(F_\sharp q, F_\sharp u)\leq C_{d, p} \cdot 3^{-k_0 p} {d_{\mathrm{TV}}(q, u)} \leq C_{d, p} m^{-p/d} (\chi^2(q, u))^{p/d} {d_{\mathrm{TV}}(q, u)}^{1 - \frac{2p}{d}}\,.$$ The claim follows from Markov’s inequality.

We now show that there are composite hypotheses that are well separated in total variation distance but nevertheless hard to distinguish on the basis of samples.

\[tv\_lb\] Fix a positive integer $n$ and a constant $\delta \in [0, 1/10]$. Given a positive integer $m$, let ${\mathcal{D}}_{m}$ be the set of probability distributions $q$ on $[m]$ satisfying $\chi^2(q, u) \leq 9$. Denote by ${\mathcal{D}}_{m,\delta}^-$ the subset of ${\mathcal{D}}_m$ of distributions satisfying ${d_{\mathrm{TV}}(q, u)} \leq \delta$ and by ${\mathcal{D}}_m^+$ the subset of ${\mathcal{D}}_m$ satisfying ${d_{\mathrm{TV}}(q, u)} \geq 1/4$. If $m = \lceil C \delta^{-1} n \log n \rceil$ for a sufficiently large universal constant $C$ and $n$ is sufficiently large, then $$\inf_\psi \big\{\sup_{q \in {\mathcal{D}}_{m}^+} {\mathbb P}_q[\psi =1] + \sup_{q \in {\mathcal{D}}_{m,\delta}^-} {\mathbb P}_q[\psi =0] \big\} \geq  .9\,,$$ where the infimum is taken over all (possibly randomized) tests based on $n$ samples.

The proof of \[tv\_lb\] follows a strategy due to @ValVal10 and @WuYan19, and our argument is a modification of theirs which permits simultaneous control of total variation and the $\chi^2$ divergence. We give the proof in \[tv\_lb\_proof\].

We now give a proof of the main theorem.

Let $\hat W$ be any estimator for the Wasserstein distance between distributions on ${\mathcal{X}}$ constructed on the basis of $n$ samples from each distribution.

Let $u$ be the uniform distribution on $[m]$, for some $m$ to be specified. Let $c^*$ be the constant appearing in the lower bound of \[random\_injection\] and define $\Delta_d = \frac 1 {16} c^* m^{-1/d}$. Given $n$ samples $X_1, \dots, X_n$ from an unknown distribution on $[m]$, define the randomized test $$\psi = \psi(X_1, \dots, X_n) {:=}{\mathds{1}}\{\hat W(F(X_1), \dots, F(X_n); F(Y_1), \dots, F(Y_n)) \leq 2 \Delta_d\}\,,$$ where $F$ is the random function constructed in \[random\_injection\] and where $Y_i$ are i.i.d. from $u$.

By \[random\_injection\], if $\delta \leq \delta_{d, p} {:=}\left(\frac{c^*}{176 C_{d, p}}\right)^{\frac{1}{1/p-2/d}}$, any $q \in {\mathcal{D}}_{m, \delta}^-$ satisfies the bound $W_p(F_\sharp q, F_\sharp u) \leq \Delta_d$ with probability at least $.9$. Likewise, for $q \in {\mathcal{D}}_{m}^+$, the bound $W_p(F_\sharp q, F_\sharp u) \geq 3  \Delta_d$ also holds with probability at least $.9$.

Define the event $A=\{|\hat W - W_p(F_\sharp q, F_\sharp u)|  \geq \Delta_d\}$. We obtain, for any $q \in {\mathcal{D}}_{m, \delta}^-$, $$\begin{aligned}
{\mathbb E}_F {\mathbb P}_{F_\sharp q, F_\sharp u} [A] & \geq {\mathbb E}_F {\mathbb P}_{F_\sharp q F_\sharp u} [\hat W > 2 \Delta_d \text{ and } W_p(F_\sharp q, F_\sharp u) \leq \Delta_d] \\
& \geq {\mathbb E}_F {\mathbb P}_{F_\sharp q F_\sharp u} [\hat W > 2\Delta_d] - {\mathbb P}[W_p(F_\sharp q, F_\sharp u) > \Delta_d] \\\
& \geq {\mathbb P}_{q} [\psi = 0] - .1\,,\end{aligned}$$ and analogously for $q \in {\mathcal{D}}_{m}^+$, $${\mathbb E}_F {\mathbb P}_{F_\sharp q, F_\sharp u} [A] \geq {\mathbb P}_q [\psi = 1] - .1\,.$$

For any estimator $\hat W$, we have $$\begin{aligned}
\sup_{\mu, \nu \in {\mathcal{P}}}
{\mathbb P}_{\mu, \nu}[
|\hat W - W_p(\mu, \nu)| \geq \Delta_d
]
& \geq
\frac 1 2
\big
(
\sup_{q \in {\mathcal{D}}_m^+} {\mathbb E}_F {\mathbb P}_{F_\sharp q, F_\sharp u} [A]+
 \sup_{q \in {\mathcal{D}}_{m,\delta}^-} {\mathbb E}_F {\mathbb P}_{F_\sharp q, F_\sharp u}[A]
\big
) \\
&  \geq \frac 12 \big(\sup_{q \in {\mathcal{D}}_m^+} {\mathbb P}_q[\psi = 1] + \sup_{q \in {\mathcal{D}}_{m,\delta}^-} {\mathbb P}_q[\psi = 0]\big) - .1\,.\end{aligned}$$

Choosing $m = \lceil C \delta^{-1} n \log n \rceil$ for a sufficiently large constant $C$ and applying \[tv\_lb\] yields that $\sup_{\mu, \nu \in {\mathcal{P}}}
{\mathbb P}_{\mu, \nu}[|\hat W - W_p(F_\sharp q, F_\sharp u)|  \geq \Delta_d] \geq .8$, and Markov’s inequality yields the claim.

Computational-statistical gaps for the spiked transport model {#sec:comp}
=============================================================

clarify the statistical price for estimating the Wasserstein distance for high-dimensional measures. shows that the curse of dimensionality can be avoided under the spiked transport model. The WPP estimator exploits the low-dimensional structure in the spiked transport model, thereby beating the worst-case rate presented in \[sec:lower\]. However, it is not clear how to make the estimator we propose computationally efficient. In this section, we give evidence that this obstruction is a fundamental obstacle, that is, that no computationally efficient estimator can beat the curse of dimensionality.

The statistical query model, first introduced in the context of PAC learning [@Kea98], is a well known computational framework for analyzing statistical algorithms. Instead of being given access to data points from a distribution, a statistical query (SQ) algorithm can approximately evaluate the expectation of arbitrary functions with respect to the distribution. This model naturally captures the power of noise-tolerant algorithms [@Kea98] and is strong enough to implement nearly all common machine learning procedures [see, e.g. @BluDwoMcS05].

We recall the following definition.

Given a distribution ${\mathcal{D}}$ on ${\mathbb{R}}^d$, for any sample size parameter $t > 0$ and function $f: {\mathbb{R}}^d \to [0, 1]$, the oracle ${\mathrm{VSTAT}}(t)$ returns a value $v \in [p - \tau, p + \tau]$, where $p = {\mathbb E}f(X)$ and $\tau = \frac 1 t \vee \sqrt{\frac{p(1-p)}{t}}$.

A query to a ${\mathrm{VSTAT}}(t)$ oracle can be simulated by using a data set of size approximately $t$. Our main result proves a lower bound against an oracle with sample size parameter $t = 2^{cd}$ for a positive constant $c$. Simulating such an oracle would require a number of samples exponential in the dimension. Nevertheless, we show that even under this strong assumption, at least $2^{cd}$ queries to the oracle are required. This result suggests that any *computationally efficient* procedure to estimate the Wasserstein distance under the spiked transport model requires an exponential number of samples. By contrast, \[sec:upper\] establishes that, information theoretically, only a polynomial number of samples are required.

We now state our main result.

\[sq\_lb\] There exists a positive universal constant $c$ such that, for any $d$, estimating $W_1({\mu^{(1)}}, {\mu^{(2)}})$ for distributions ${\mu^{(1)}}$ and ${\mu^{(2)}}$ on ${\mathbb{R}}^d$ satisfying the spiked transport assumption with $k = 1$ to accuracy $\Theta(1/\sqrt d)$ with probability at least $2/3$ requires at least $2^{cd}$ queries to ${\mathrm{VSTAT}}(2^{cd})$.

Our proof is based on a construction due to [@DiaKanSte17] [see also @BubLeePri19]. We defer the details to the appendix.

Proofs of Lower Bounds
======================

Proof of \[deviation\_lb\] {#deviation_lb_proof}
--------------------------

We reduce from the spiked covariance model. By homogeneity, we may assume that $\sigma = 1$. Let $\mu^{(1)}$ be the standard Gaussian measure on ${\mathbb{R}}^d$, and for $\mu^{(2)}$ we take either the standard Gaussian measure or the distribution of a centered Gaussian with covariance $I + \beta v v^\top$, where $\|v\| = 1$ and $\beta > 0$ is to be specified. As long as $\beta \lesssim 1$, the measure $\mu^{(2)}$ is a $O(1)$-Lipschitz pushforward of the Gaussian measure. Hence, it satisfies $T_p(O(1))$ [@Goz07 Corollary 20].

Note that if $\mu^{(2)}$ has covariance $I + \beta v v^\top$, then $$W_p({\mu^{(1)}}, {\mu^{(2)}}) \geq W_1({\mu^{(1)}}, {\mu^{(2)}}) = W_1({\mathcal{N}}(0, 1), {\mathcal{N}}(0, 1+\beta)) \gtrsim \beta\,.$$ However, @CaiMaWu15 [Proposition 2] establish that the minimax testing error for $$\mathrm{H}_0: {\mathcal{N}}(0, I) \text{ vs. } \mathrm{H}_1: {\mathcal{N}}(0, I + \beta v v^\top), \|v\| = 1$$ is bounded below by a constant when $\beta \lesssim \sqrt{d/n}$. A standard application of Le Cam’s two-point method [@Tsy09] yields the claim.

Proof of \[tv\_lb\] {#tv_lb_proof}
-------------------

We require the existence of two distributions on ${\mathbb{R}}_+$, which will serve as the building blocks of our construction.

\[priors\] For any integer $L \geq 0$ and ${\varepsilon}\in [0, 1/6]$, there exists a pair of random variables $U$ and $V$ with the following properties:

-   ${\mathbb E}U^j = {\mathbb E}V^j \quad \forall j \leq L$

-   $U, V \in [0, 16{\varepsilon}^{-1} L^2]$ almost surely

-   ${\mathbb E}U = {\mathbb E}V = 1$ and ${\mathbb E}U^2 = {\mathbb E}V^2 \leq 6$.

-   ${\mathbb E}|U - 1| \leq 12 {\varepsilon}$ but ${\mathbb E}|V - 1| \geq 1$.

The proof is deferred to \[sec:priors\_proof\].

The proof follows closely the approach of @WuYan19 [Proposition 1]. We first employ a standard argument showing that we can consider the *Poissonalized* setting. It is trivial to see that given samples $X_1, \dots, X_n$ from a distribution $q$ on $[m]$, the *counts* $N_i = N_i(X_1, \dots, X_n) {:=}|\{j \in [n]: X_j = i\}|$ are sufficient for $q$. We therefore consider tests $\psi$ based on count vectors. Note that, under $q$, the count vector $(N_1, \dots, N_m)$ has distribution ${\mathrm{Multinomial}}(n, q)$.

Define $$R_n {:=}\inf_\psi \big\{\sup_{q \in {\mathcal{D}}_{m}^+} {\mathbb P}_q[\psi =1 ] + \sup_{q \in {\mathcal{D}}_{m,\delta}^-} {\mathbb P}_q[\psi =0 ]\}\,.$$ We aim to prove a lower bound on $R_n$.

Let $\rho > 0$, and let for $n \geq 1$ let $\{\psi_n\}$ be a set of near optimal tests for a fixed sample size; i.e. $$\sup_{q \in {\mathcal{D}}_{m,\delta}^- \cup {\mathcal{D}}_m^+} {\mathbb P}_q[\psi_n \neq {\mathds{1}}\{q \in {\mathcal{D}}_{m, \delta}^-\}] \leq R_n + \rho\,.$$ Define set of approximate probability vectors $$\begin{aligned}
\tilde {\mathcal{D}}_{m,\delta}^- &{:=}\left\{q \in {\mathbb{R}}^m_+ : \left|\sum_{i=1}^m q_i - 1\right| \leq \delta , \frac{q}{\sum_{i=1}^m q_i} \in {\mathcal{D}}_{m,\delta}^-\right\} \\
\tilde {\mathcal{D}}_{m}^+ & {:=}\left\{q \in {\mathbb{R}}^m_+ : \left|\sum_{i=1}^m q_i - 1\right| \leq \delta , \frac{q}{\sum_{i=1}^m q_i} \in {\mathcal{D}}_m^+\right\}\,.\end{aligned}$$

We let $\tilde {\mathcal{D}}_{m, \delta} {:=}\tilde {\mathcal{D}}_{m,\delta}^- \cup \tilde {\mathcal{D}}_{m}^+$. Given $q \in \tilde {\mathcal{D}}_{m, \delta}$, define the renormalization $\bar q=\sum_{i=1}^m q_i/q$. We then define $$\tilde R_n {:=}\inf_\psi \big\{\sup_{q \in \tilde {\mathcal{D}}_{m}^+} {\mathbb P}_q[\psi =1] + \sup_{q \in \tilde {\mathcal{D}}_{m,\delta}^-} {\mathbb P}_q[\psi =0]\big\}\,,$$ where the infimum is taken over all estimators based on the counts $N_1, \dots, N_m$ and where ${\mathbb P}_q$ indicates the probability when $N_1, \dots, N_m$ are independent and $N_i \sim {\mathrm{Pois}}(n q_i)$ for all $i \in [m]$. We set $N = \sum_{i=1}^m N_i$, and note that, conditioned on $N = n'$, the count vector $(N_1, \dots, N_m)$ has distribution ${\mathrm{Multinomial}}(n', \bar q)$.

We define a test $\tilde \psi$ based on these Poissonalized counts by setting $$\tilde \psi(N_1, \dots, N_m) {:=}\psi_N(N_1, \dots, N_m)\,.$$ This definition along with the near optimality of $\psi_{n'}$ for $n' \geq 0$ implies $$\begin{aligned}
\sup_{q \in \tilde {\mathcal{D}}_{m}^+} {\mathbb P}_q[\tilde \psi =1] + \sup_{q \in \tilde {\mathcal{D}}_{m,\delta}^-} {\mathbb P}_q[\tilde \psi =0] & \leq \sum_{n' \geq 0} R_{n'} {\mathbb P}_q[N = n'] + \rho\\
& \leq R_{n/2} + {\mathbb P}_q[N < n/2] + \rho\,,\end{aligned}$$ where the last inequality follows from the fact that $R_{n'} \leq 1$ for all $n' \geq 0$ and $R_{n'}$ is non-increasing in $n'$. Since $N = {\mathrm{Pois}}(n \sum_{i=1}^m q_i)$ and $\sum_{i=1}^m q_i \geq 3/4$, a standard Chernoff bound implies ${\mathbb P}[N < n/2] \leq \exp(-Cn)$. Since $\rho$ was arbitrary, we obtain that $$\tilde R_n \leq R_{n/2} + \exp(-Cn)\,.$$

To prove a lower bound on $\tilde R_n$, we consider random vectors $$\begin{aligned}
Q & = \frac 1 m (U_1, \dots, U_m) \\
Q'& = \frac 1 m (V_1, \dots, V_m)\,,\end{aligned}$$ where $U_i$ and $V_i$ for $i \in [m]$ are independent copies of $U$ and $V$ constructed in \[priors\] with ${\varepsilon}= \frac{1}{24}\delta$. Conditioned on $Q$ and $Q'$, let $N$ and $N'$ be count vectors with independent entries generated by $N_i \sim {\mathrm{Pois}}(nQ_i)$ and $N'_i \sim {\mathrm{Pois}}(nQ'_i)$. Let us denote by $\mathrm{P}$ and $\mathrm{P'}$ the distributions of $N$ and $N'$ respectively. Under $\mathrm{P}$ and $\mathrm{P'}$, the entries of $N$ and $N'$ are i.i.d. Poisson mixtures, so applying @WuYan19 [Lemma 4] yields $${d_{\mathrm{TV}}(\mathrm{P}, \mathrm{P'})} \leq m \left(\frac{8 e n L}{{\varepsilon}m}\right)^L \,.$$ Let $E = \{Q \in \tilde {\mathcal{D}}_{m, \delta}^-\}$ and $E' =\{Q' \in \tilde {\mathcal{D}}_{m}^+\}$. By \[approximate-probabilities\], ${\mathbb P}[E^C]$ and ${\mathbb P}[{E'}^C]$ are each at most $C \frac{L^4}{\delta^2 m}$.

Let $\pi_E$ be the law of $Q$ conditioned on $E$, and define $\pi'_{E'}$ analogously, and let $\mathrm{P}_E$ and $\mathrm{P}'_{E'}$ be the laws of $N$ and $N'$ under these priors. We obtain for any estimator $\psi$ based on count vectors $$\begin{aligned}
\sup_{q \in \tilde {\mathcal{D}}_{m}^+} {\mathbb P}_q[\tilde \psi =1] + \sup_{q \in \tilde {\mathcal{D}}_{m,\delta}^-} {\mathbb P}_q[\tilde \psi =0] & \geq  \int {\mathbb P}_{q'}[\psi = 1] {\, \mathrm{d}}\pi'_{E'}(q') + \int {\mathbb P}_q[\psi =0] {\, \mathrm{d}}\pi_E(q) \\
& \geq 1 - {d_{\mathrm{TV}}(\mathrm{P}_E, \mathrm{P}'_{E'})} \\
& \geq 1 - {d_{\mathrm{TV}}(\mathrm{P}, \mathrm{P'})} - C \frac{L^4}{\delta^2 m}\,.\end{aligned}$$ Choosing $L = c\frac{\delta m}{n}$ for a sufficiently small constant $c$ yields that $$\tilde R_n \geq 1-  m \exp(- C \frac{\delta m}{n}) - C \frac{\delta^2 m^3}{n^4}\,.$$ Therefore $$R_n \geq 1 - m \exp(- C \frac{\delta m}{n}) - C \frac{\delta^2  m^3}{n^4} - \exp(-C n)\,,$$ and choosing $m = \lceil C \delta^{-1} n \log n \rceil$ for $C$ a sufficiently large constant and $n$ sufficiently large yields the claim.

Proof of Proposition \[priors\] {#sec:priors_proof}
-------------------------------

First, the reduction of @WuYan19 [Lemma 7] implies that it suffices to construct random variables $Y$ and $Y'$ such that

-   ${\mathbb E}Y^j = {\mathbb E}{Y'}^j \quad \forall 0 \leq j < L$

-   $Y, Y' \in [1, 16 {\varepsilon}^{-1} L^2]$ a.s.

-   ${\mathbb E}Y = {\mathbb E}Y' \leq 6$

-   ${\mathbb E}\frac 1 Y \geq 1 - 6 {\varepsilon}$ but ${\mathbb E}\frac{1}{Y'} \leq \frac 1 2$.

Indeed, applying their construction yields $U$ and $V$ satisfying the first three requirements of \[priors\] as well as ${\mathbb P}[U = 0] \leq 6 {\varepsilon}$ and ${\mathbb P}[V = 0] \geq \frac 1 4$. Since the supports of $U$ and $V$ lie in $\{0\} \cup [1, +\infty)$, we have ${\mathbb E}|U - 1| = 2 {\mathbb P}[U = 0] \leq 12 {\varepsilon}$ and ${\mathbb E}|V - 1| = 2 {\mathbb P}[V = 0] \geq 1$, as desired. We therefore focus on constructing such a $Y$ and $Y'$.

By @WuYan19 [Lemma 7] combined with @Tim94 [Section 2.11.1], there exist random variables $X$ and $X'$ supported on $[1, 16 L^2]$ such that ${\mathbb E}X^j = {\mathbb E}{X'}^j$ for $0 \leq j < L$ and ${\mathbb E}\frac 1 X - {\mathbb E}\frac{1}{X'} \geq \frac 1 2$. Let $\mathrm{P}_{\varepsilon}$ and $\mathrm{P}'_{\varepsilon}$ denote the distribution of ${\varepsilon}^{-1} X$ and ${\varepsilon}^{-1} X'$, respectively.

Let $$\begin{aligned}
\Delta_{\varepsilon}&{:=}\int \frac{1}{(y-1)(y-2)} {\, \mathrm{d}}\mathrm{P}_{\varepsilon}(y)\\
\Delta'_{\varepsilon}&{:=}\int \frac{1}{(y'-1)(y'-2)} {\, \mathrm{d}}\mathrm{P}'_{\varepsilon}(y')\\
{\mathcal{Z}}_{\varepsilon}& {:=}\int \frac{1}{y - 2} {\, \mathrm{d}}\mathrm{P}_{\varepsilon}(y) - \int \frac{1}{y' - 1} {\, \mathrm{d}}\mathrm{P}'_{\varepsilon}(y')\,.\end{aligned}$$ We define two new distributions $\mathrm{Q}$ and $\mathrm{Q'}$ by $$\begin{aligned}
\mathrm{Q}(\mathrm{d}y) & = \delta_1(\mathrm{d}y) + \frac{1}{{\mathcal{Z}}_{\varepsilon}}\left(\frac{1}{(y-1)(y - 2)} \mathrm{P}_{\varepsilon}(\mathrm{d}y) - \Delta_{\varepsilon}\delta_1(\mathrm{d}y)\right) \label{eq:py}\\
\mathrm{Q'}(\mathrm{d}y') & = \delta_2(\mathrm{d}y') + \frac{1}{{\mathcal{Z}}_{\varepsilon}}\left(\frac{1}{(y'-1)(y' - 2)} \mathrm{P}'_{\varepsilon}(\mathrm{d}y') - \Delta_{{\varepsilon}'} \delta_2(\mathrm{d}y')\right)\label{eq:py-prime}\end{aligned}$$

By \[rescaled-prior-bounds\], $$\begin{aligned}
\Delta_{\varepsilon}, \Delta_{{\varepsilon}'} & \in [0, \frac{9}{5} {\varepsilon}^2] \\
{\mathcal{Z}}_{\varepsilon}& \geq \frac{3}{10} {\varepsilon}\,,\end{aligned}$$ which implies in particular that both $Q$ and $Q'$ are probability distributions.

Let $Y \sim Q$ and $Y' \sim Q'$. We first check the last three conditions. Clearly $Y$ and $Y'$ are supported on $[1, 16 {\varepsilon}^{-1} L^2]$, and \[prior-first-moments-match\] implies that ${\mathbb E}Y = {\mathbb E}Y' \leq 6$. We have ${\mathbb E}\frac 1 Y \geq {\mathbb P}[Y = 1] = 1 - \frac{\Delta_{\varepsilon}}{{\mathcal{Z}}_{\varepsilon}} \geq 1 - 6 {\varepsilon}$, and since $Y' \geq 2$ almost surely the bound ${\mathbb E}\frac{1}{Y'} \leq \frac 1 2$ is immediate.

It remains to check the moment-matching condition. Any polynomial $p(y)$ of degree at most $L-1$ can be written $$p(y) = (y-1)(y-2)q(y) + \alpha y + \beta\,,$$ where $q(y)$ has degree less than $L-1$. Then $${\mathbb E}p(Y) - {\mathbb E}p(Y') = {\mathbb E}(Y-1)(Y-2)q(Y) - {\mathbb E}(Y'-1)(Y'-2)q(Y') + \alpha ({\mathbb E}Y - {\mathbb E}Y')\,.$$ The last term vanishes because ${\mathbb E}Y = {\mathbb E}Y'$, and $${\mathbb E}(Y-1)(Y-2)q(Y) - {\mathbb E}(Y'-1)(Y'-2)q(Y') = \frac{1}{{\mathcal{Z}}_{\varepsilon}} \left({\mathbb E}q({\varepsilon}^{-1} X) - {\mathbb E}q({\varepsilon}^{-1} X')\right) = 0\,,$$ since ${\mathbb E}X^j = {\mathbb E}{X'}^j$ for all $j < L - 1$.

Therefore ${\mathbb E}p(Y) = {\mathbb E}p(Y')$ for all polynomials of degree at most $L - 1$.

Proof of \[sq\_lb\]
-------------------

We first establish the existence of a probability distribution on ${\mathbb{R}}$ which agrees with ${\mathcal{N}}(0, 1)$ on many moments, but is far from ${\mathcal{N}}(0, 1)$ in Wasserstein distance.

\[one\_dimensional\_distribution\] There exists a $O(1)$-subgaussian distribution $A$ on ${\mathbb{R}}$ that satisfies the following requirements.

-   $A$ agrees with ${\mathcal{N}}(0, 1)$ on the first $2m - 1$ moments.

-   $W_1(A, {\mathcal{N}}(0, 1)) = \Omega(1/\sqrt m)$.

-   ${\chi^2(A, {\mathcal{N}}(0, 1))} = \exp(O(m))$.

By @DiaKanSte17 [Corollary 4.4], for any $\delta \in (0, 1)$ we can find an atomic measure $Q$ supported on $m$ points in $[-O(\sqrt m), O(\sqrt m)]$ such that $A {:=}Q * {\mathcal{N}}(0, \delta)$ matches ${\mathcal{N}}(0, 1)$ on the first $2m - 1$ moments. By @DiaKanSte17 [Lemma 4.6], this distribution satisfies the bound ${\chi^2(A, {\mathcal{N}}(0, 1))} = \exp(O(m))/\sqrt \delta$. Moreover, the distribution $Q$ is supported on the zeros of a rescaled version of the $m$th Hermite polynomial, and by @BubLeePri19 [Lemma A.1] these zeros are $\Omega(1/\sqrt m)$ apart.

Since the points in the support of $Q$ are $\Omega(1/\sqrt m)$ apart, there exists a constant $c$ such that the union of balls of radius $c/\sqrt m$ centered at the support of $Q$ covers at most half of the interval $[-1, 1]$, and since the Gaussian density is bounded below on this interval, a constant fraction of the mass of ${\mathcal{N}}(0, 1)$ is located at distance at least $c/\sqrt m$ from $Q$. Hence $W_1(Q, {\mathcal{N}}(0, 1)) = \Omega(1/\sqrt m)$. Clearly $W_1(A, Q) = O(\sqrt \delta)$.

Therefore, if we choose $\delta = O(1/m)$, then $$W_1(A, {\mathcal{N}}(0, 1)) \geq W_1(Q, {\mathcal{N}}(0, 1)) - W_1(A, Q) = \Omega(1/\sqrt m)$$ and ${\chi^2(A, {\mathcal{N}}(0, 1))} = O(\sqrt m) \exp(O(m)) = \exp(O(m))$.

Finally, we show that $A$ is $O(1)$-subgaussian, and therefore satisfies ${{T_{1}(C)}}$ for a positive constant $C$. Standard facts [see @Ver18] imply that it suffices to show that if $X \sim A$, then $\|X\|_k = O(\sqrt k)$ for all $k$. This clearly holds for $k \leq 2m -1$, since ${\mathcal{N}}(0, 1)$ is itself $1$-subgaussian and the first $2m -1$ moments of $A$ and ${\mathcal{N}}(0, 1)$ agree. On the other hand, for $k \geq 2m$, if $Y \sim Q$ and $Z \sim {\mathcal{N}}(0, 1)$, then $X = Y + \sqrt \delta Z \sim A$, and $$\|X + \sqrt \delta Z\|_k \leq \|X\|_\infty + \sqrt \delta \|Z\|_k \lesssim \sqrt m + \sqrt{k/m} \lesssim \sqrt k\,,$$ as desired.

The separation $W_1(A, {\mathcal{N}}(0, 1)) = \Omega(1/\sqrt m)$ in \[one\_dimensional\_distribution\] is easily seen to be tight. Indeed, @RigWee19 [Corollary 2] show that if $\mu$ and $\nu$ are $O(1)$-subgaussian and agree on their first $O(m)$ moments, then $W_1(\mu, \nu) = O(1/\sqrt m)$.

By planting the distribution constructed in \[one\_dimensional\_distribution\] in a random direction, we obtain two high-dimensional measures satisfying the spiked transport model.

\[A\_spike\] Let $v$ be a unit vector in ${\mathbb{R}}^d$, and denote by $P_v$ the distribution on ${\mathbb{R}}^d$ of the random variable $X v + Z$, where $X \sim A$ and $Z \sim {\mathcal{N}}(0, I_d - v v^\top)$ is independent of $X$. Then $\mu^{(1)}=P_v$ and $\mu{(2)}={\mathcal{N}}(0, I_d)$ satisfy the spiked transport model , and $W_1(P_v, {\mathcal{N}}(0, I_d)) = \Omega(1/\sqrt m)$.

If we let $\xi \sim {\mathcal{N}}(0, 1)$, then ${\mathcal{N}}(0, I_d)$ is the law of $\xi v + Z$, where $\xi$ and $Z$ are independent and $Z \sim {\mathcal{N}}(0, I_d - v v^\top)$. Denoting by ${\mathcal{U}}$ the span of $v$, we see that $P_v$ and ${\mathcal{N}}(0, I_d)$ satisfy  with $X^{(1)} = X v$ and $X^{(2)} = \xi v$. By \[prop:wpk-is-wp,one\_dimensional\_distribution\], $W_1(P_v, {\mathcal{N}}(0, I_d)) = W_1(A, {\mathcal{N}}(0, 1)) = \Omega(1/\sqrt m)$.

The proof of \[sq\_lb\] follows from a framework due to @FelGriRey17, from which the following result is extracted. Given distributions $P_1$, $P_2$, and $Q$, define $$\chi^2_{Q}(P_1, P_2) {:=}\int \left(\frac{dP_1}{dQ} - 1\right)\left(\frac{dP_2}{dQ} - 1\right) \, {\, \mathrm{d}}Q\,.$$ We call a set ${\mathcal{P}}$ of distributions $(\gamma, \beta)$ correlated with respect to $Q$ if for all $P_i, P_j \in {\mathcal{P}}$, $$\chi^2_Q(P_i, P_j) \leq \left\{\begin{array}{ll}
\beta & \text{if $i = j$,} \\
\gamma & \text{if $i \neq j$.}
\end{array}\right.$$ We then have the following.

\[decision\_lb\] Let ${\mathcal{Q}}$ be a set of distributions, and let $Q$ be a reference distribution. Suppose that there exists a set ${\mathcal{P}}\subseteq {\mathcal{Q}}$ such that ${\mathcal{P}}$ is $(\gamma, \beta)$ correlated with respect to $Q$. Then any SQ algorithm that distinguishes queries from $P = Q$ and $P \in {\mathcal{Q}}$ requires at least $|{\mathcal{P}}| \gamma/3\beta$ queries to ${\mathrm{VSTAT}}(1/2 \gamma)$

By choosing $\gamma' = \gamma$ in @FelGriRey17 [Lemma 3.10], we obtain that the set ${\mathcal{P}}$ satisfies $\mathrm{SDA}({\mathcal{P}}, Q, 2\gamma) \geq |{\mathcal{P}}|\gamma/(\beta - \gamma) \geq |{\mathcal{P}}| \gamma/\beta$. Then @FelGriRey17 [Theorem 3.7] implies that distinguishing $Q$ from ${\mathcal{Q}}$ requires at least $|{\mathcal{P}}| \gamma/3\beta$ queries to ${\mathrm{VSTAT}}(1/2 \gamma)$.

We can now prove the lower bound.

Let ${\mathcal{Q}}$ be the set $\{P_v : v \in {\mathbb{R}}^d, \|v\| = 1\}$. The Johnson-Lindenstrauss lemma [@JohLin84] implies that for any $\delta \in (0, 1)$, there exists a set of $2^{\Omega(\delta^2 d)}$ unit vectors in ${\mathbb{R}}^d$ with pairwise inner product at most $\delta$. Denote by ${\mathcal{S}}$ a set of such vectors, and set ${\mathcal{P}}{:=}\{P_v : v \in {\mathcal{S}}\} \subseteq {\mathcal{Q}}$.

Write $\Delta$ for ${\chi^2(A, {\mathcal{N}}(0, 1))}$, and recall that $\Delta = \exp(O(m))$. By @DiaKanSte17 [Lemma 3.4], $\chi^2_{{\mathcal{N}}(0, I)}(P_v, P_{v'}) \leq |v \cdot v'|^{2m} \Delta$, and the set ${\mathcal{P}}$ is therefore $(\delta^{2m} \Delta, \Delta)$ correlated. By \[decision\_lb\], any SQ algorithm that distinguishes between $\mu = {\mathcal{N}}(0, 1)$ and $\mu \in {\mathcal{P}}$ with probability at least $2/3$ requires $2^{\Omega(\delta^2 d)}\delta^{2m}$ queries to ${\mathrm{VSTAT}}(\delta^{-2m}/2 \Delta)$.

Let $\delta > 0$ be a constant small enough that $\delta^{-2m}/2\Delta = \exp(\Omega(m))$. If $m = c d$ for a sufficiently small positive constant $c$, then $2^{\Omega(\delta^2 d)} \delta^{2m} = \exp(\Omega(d))$. An SQ algorithm to distinguish queries from ${\mathcal{N}}(0, 1)$ from those from a distribution in ${\mathcal{Q}}$ with probability at least $2/3$ therefore requires $\exp(\Omega(d))$ queries to ${\mathrm{VSTAT}}(\exp(\Omega(d))$. Therefore, by \[A\_spike\], any SQ algorithm which estimates $W_1$ under the spiked transport model to accuracy $\Theta(1/\sqrt m) = \Theta(1/\sqrt d)$ requires $\exp(\Omega(d))$ queries to ${\mathrm{VSTAT}}(\exp(\Omega(d))$, as claimed.

Supplementary material {#sec:omitted}
======================

Additional lemmas {#subsec:lemmas}
-----------------

\[lem:vk\_op\_covering\] There exists a universal constant $c$ such that for all ${\varepsilon}\in (0, 1)$, the covering number of ${\mathcal{V}}_k({\mathbb{R}}^d)$ satisfies $$\log {\mathcal{N}}({\mathcal{V}}_k, {\varepsilon}, {\left\|\cdot\right\|_{\mathrm{op}}}) \leq dk \log \frac{c\sqrt k} {\varepsilon}\,.$$

If we denote by ${\mathcal{S}}_d$ the unit sphere in ${\mathbb{R}}^d$, then ${\mathcal{V}}_k({\mathbb{R}}^d) \subseteq {\mathcal{S}}_d^{\times k}$, and $${\left\|U - V\right\|_{\mathrm{op}}}^2 \leq \sum_{i = 1}^k \|U_i - V_i\|^2\,,$$ where $\{U_i\}$ (resp. $\{V_i\}$) are the rows of $U$ (resp. $V$). The claim follows immediately from the fact that $\log {\mathcal{N}}({\mathcal{S}}_d, {\varepsilon}, \|\cdot\|) \leq d \log \frac c {\varepsilon}$ for a universal constant $c$.

\[prop:subg\_vec\_norm\] Let $\mu$ be a centered distribution satisfying ${{T_{1}(\sigma^2)}}$. If $X \sim \mu$, then $({\mathbb E}\|X\|^p)^{1/p} \lesssim \sqrt{d p}$.

By the monotonicity of $L_p$ norms, we may assume that $p$ is a positive even integer. The claim then follows from \[prop:norm\_subG\] and standard bounds on the moments of subgaussian random variables [@Ver18].

The following proposition provides a sort of converse to \[prop:wpk-is-wp\]: if the $k$-dimensional Wasserstein distance agrees with the Wasserstein distance, then there is an optimal coupling between the measures which acts only on a $k$-dimensional subspace.

\[prop:low-dim-coupling\] If ${\tilde W_{p,k}}(\mu, \nu) = W_p(\mu, \nu)$, then there exists a coupling $\gamma$ optimal for $\mu$ and $\nu$ such that $X - Y$ is supported on a $k$-dimensional subspace.

Let $\gamma$ be an optimal coupling for $\mu$ and $\nu$, and write $$U = \operatorname*{argmax}_{V \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} W_p(\mu_V, \nu_V)\,.$$ Then $$\begin{aligned}
W_p^p(\mu, \nu)  =  \int \|x - y\|^p {\, \mathrm{d}}\gamma(x, y)  \geq \int \|U (x - y)\|^p {\, \mathrm{d}}\gamma(x, y)  \geq {\tilde W_{p,k}}^p(\mu, \nu)\,,\end{aligned}$$ where the first inequality uses the fact that ${\left\|U\right\|_{\mathrm{op}}} = 1$ and the second uses the fact that if $(X, Y) \sim \gamma$, then $U X$ and $U Y$ are distributed according to $\mu_U$ and $\nu_U$, respectively. If the first equality holds, then we must have $\|X - Y\| = \|U(X - Y)\|$ $\gamma$-almost surely, and if the second inequality holds then this coupling is optimal for $\mu_U$ and $\nu_U$

\[prop:second\_moment\_to\_subG\] If a real-valued random variable $X$ is $\sigma^2$-subgaussian, then $$\label{eq:second_moment}
{\mathbb E}e^{\lambda^2 X^2} \leq e^{4\lambda^2 \sigma^2} \quad \quad \forall |\lambda| \leq \frac{1}{2\sigma}\,.$$

Conversely, if  holds and ${\mathbb E}X = 0$, then $X$ is $8\sigma^2$-subgaussian.

For the first claim, let $Z$ be a standard Gaussian random variable independent of $X$. Then, if $\lambda^2 \sigma^2 < 1/4$, $${\mathbb E}e^{\lambda^2 X^2} = {\mathbb E}e^{\sqrt 2 \lambda Z X} \leq {\mathbb E}e^{\lambda^2 \sigma^2 Z^2} = \frac{1}{(1 - 2 \lambda^2 \sigma^2)^{1/2}} \leq e^{4 \lambda^2 \sigma^2}\,,$$ where the last step uses the inequality $(1-x)^{-1} \leq e^{2x}$ for $0 \leq x \leq \frac 12$.

Conversely, suppose that  holds and ${\mathbb E}X = 0$. Then, for $|\lambda| \leq \frac{1}{2 \sigma}$, we have $${\mathbb E}e^{\lambda X} \leq {\mathbb E}\lambda X + {\mathbb E}e^{\lambda^2 X^2} \leq e^{4 \lambda^2 \sigma^2}\,.$$ If $|\lambda| > \frac{1}{2 \sigma}$, then $${\mathbb E}e^{\lambda X} \leq {\mathbb E}e^{2 \lambda^2 \sigma^2 + \frac{X^2}{8 \sigma^2}} \leq e^{2 \lambda^2 \sigma^2 + 1} \leq e^{4 \lambda^2 \sigma^2}\,.$$

\[prop:norm\_subG\] If $X$ on ${\mathbb{R}}^k$ is $\sigma^2$-subgaussian, then ${\varepsilon}\|X\|$ is $(8 k \sigma^2)$-subgaussian, where ${\varepsilon}$ is a Rademacher random variable independent of $X$.

For $|\lambda| < \frac{1}{2 \sqrt k\sigma}$, $$\begin{aligned}
{\mathbb E}\exp(\lambda^2 ({\varepsilon}\|X\|)^2) & = {\mathbb E}\exp(\lambda^2 \|X\|^2)  \leq {\mathbb E}\exp(k \lambda^2 (e_i^\top X)^2)  \leq \exp(4 k \lambda^2 \sigma^2)\,.\end{aligned}$$ Applying \[prop:second\_moment\_to\_subG\] yields the claim.

\[subgaussian\_wp\_rate\] Let $p' \in [1, 2]$. If $\rho_n$ on ${\mathbb{R}}^k$ satisfies ${{T_{p'}(\sigma^2)}}$, then for any $p \in [1, \infty)$, $${\mathbb E}W_{p}\left(\rho_n, \rho\right) \leq {r_{p, k}(n)} {:=}c_{p} \sigma \sqrt k \left\{
\begin{array}{ll}
n^{-1/2p} & \text{if $k < 2p$} \\
n^{-1/2p} (\log n)^{1/p} & \text{if $k = 2p$} \\
n^{-1/k} & \text{if $k > 2p$}
\end{array}
\right.$$

implies that $X \sim \rho$ is $\sigma^2$ subgaussian. in particular, $X$ satisfies $({\mathbb E}\|X\|^q)^{1/q} \leq \sigma \sqrt{qk}$. Choosing $q = 3p$ and using @Lei18 [Theorem 3.1] implies the claim.

\[packing-covering-duality\] If ${\mathcal{N}}_{\varepsilon}({\mathcal{X}}) \geq c {\varepsilon}^{-d}$, then for any positive integer $m$ there exists a subset of ${\mathcal{X}}$ of cardinality $m$ such that each pair of points is separated by at least $\frac 1 2 \left(\frac c m\right)^{1/d}$.

Let $S$ be a maximal subset of ${\mathcal{X}}$ such that each pair of points in $S$ is separated by at least $\frac 1 2 \left(\frac c m\right)^{1/d}$. Then, for any point $x \in X$, there must be an $s \in S$ such that $d(x, s) < \frac 1 2 \left(\frac c m\right)^{1/d}$, since otherwise $x$ could be added to $S$. Therefore, balls of radius $\frac 1 2 \left(\frac c m\right)^{1/d}$ around the points $s \in S$ cover $X$, which implies $$|S| \geq {\mathcal{N}}(X, \left(\frac c m\right)^{1/d}) \geq m\,,$$ as claimed.

\[random\_partition\] Let $Q_1, \dots, Q_\ell$ be a partition of ${\mathcal{X}}$. If $F$ is a uniform random bijection from $[m]$ to ${\mathcal{G}}_m$, then $${\mathbb E}\sum_{i=1}^\ell \left| F_\sharp q(Q_i) - F_\sharp u (Q_i) \right| \leq \|q - u\|_1 \wedge \sqrt{\frac{\ell \cdot \chi^2(q, u)} m}$$

The first bound follows immediately from the triangle inequality: $$\sum_{i=1}^\ell \left| F_\sharp q(Q_i) - F_\sharp u (Q_i) \right| \leq \sum_{i=1}^\ell \sum_{j: F(j) \in Q_i} \left|q(j) - \frac 1 m \right| = \|q - u\|_1\,.$$

We now show the second bound. Fix $i \in [\ell]$. Under the random choice of $F$, the quantity $F_\sharp q(Q_i) - F_\sharp u (Q_i) = \sum_{j: F(j) \in Q_i} \left(q(j) - \frac{1}{m}\right)$ is a sum of $|Q_i|$ terms selected uniformly *without replacement* from the set $\Delta {:=}\{q(j) - \frac 1 m: j \in [m]\}$. By @Hoe63 [Theorem 4], since $x \mapsto |x|$ is continuous and convex, $${\mathbb E}\left| \sum_{j: F(j) \in Q_i} \left(q(j) - \frac{1}{m}\right)\right| \leq {\mathbb E}\left| \sum_{j=1}^{|Q_i|} \Delta_j \right|\,,$$ where $\Delta_j$ are selected uniformly *with replacement* from $\Delta$. Note that ${\mathbb E}\Delta_j = 0$, and ${\mathbb E}\Delta_j^2 = \frac 1 m \sum_{j=1}^m (q(j) - \frac 1 m )^2 = m^{-2} \chi^2(q, u)$. Jensen’s inequality then yields $${\mathbb E}\left| \sum_{j=1}^{|Q_i|} \Delta_j \right| \leq \frac{\sqrt{|Q_i|\chi^2(q, u)} }{m} \,.$$ The Cauchy-Schwarz inequality finally implies $${\mathbb E}\sum_{i=1}^\ell \left| \sum_{j: F(j) \in Q_i} \left(q(j) - \frac{1}{m}\right)\right| \leq \sum_{i=1}^\ell\frac{\sqrt{|Q_i|\chi^2(q, u)}}{m}  \leq \sqrt{\frac {\ell \cdot \chi^2(q, u)} m} \,.$$

\[rescaled-prior-bounds\] Let ${\varepsilon}\in [0, 1/6]$. If $X, X' \geq 1$ almost surely and ${\mathbb E}\frac 1 X - {\mathbb E}\frac{1}{X'} \geq \frac 12$, then $${\mathbb E}\frac{1}{({\varepsilon}^{-1} X - 1)({\varepsilon}^{-1} X - 2)}, {\mathbb E}\frac{1}{({\varepsilon}^{-1} X' - 1)({\varepsilon}^{-1} X' - 2)} \in [0, \frac 9 5 {\varepsilon}^2]$$ and $${\mathbb E}\frac{1}{{\varepsilon}^{-1} X - 2} - {\mathbb E}\frac{1}{{\varepsilon}^{-1} X' - 1} \geq \frac{3}{10}{\varepsilon}\,.$$

If $x \geq 1$ and ${\varepsilon}\leq 1/6$, then $(x - {\varepsilon})(x - 2{\varepsilon}) \geq \frac{5}{9}$, so $$\frac{1}{({\varepsilon}^{-1} x - 1)({\varepsilon}^{-1} x - 2)} = \frac{{\varepsilon}^2}{(x - {\varepsilon})(x - 2{\varepsilon})} \leq  \frac 9 5 {\varepsilon}^2\,,$$ and since $X, X' \geq 1$ almost surely the first claim is immediate.

Similarly, if $x' \geq 1$ and ${\varepsilon}\leq 1/6$, we have $\frac{{\varepsilon}}{x'(x'-{\varepsilon})} \leq \frac 1 5$; hence almost surely $$\frac{1}{{\varepsilon}^{-1} X - 2} - \frac{1}{{\varepsilon}^{-1} X' - 1} \geq \frac{{\varepsilon}}{X} - \frac{{\varepsilon}}{X'} - \frac{{\varepsilon}^2}{X'(X' - {\varepsilon})} \geq {\varepsilon}\left(\frac 1 X - \frac{1}{X'} - \frac 1 5\right)\,,$$ and taking expectations yields the claim.

\[prior-first-moments-match\] Let $\mathrm{Q}$ and $\mathrm{Q}'$ be defined as in and . If $Y \sim \mathrm{Q}$ and $Y' \sim \mathrm{Q'}$, then $${\mathbb E}Y = {\mathbb E}Y' \leq 6\,.$$

It follows directly from the definition that $$\begin{split}
{\mathbb E}Y & = 1 + \frac 1{{\mathcal{Z}}_{\varepsilon}} \left(\int \frac{y}{(y-1)(y-2)} \mathrm{P_{\varepsilon}}(\mathrm d y) - \Delta_{\varepsilon}\right)\\
& = 1 + \frac 1{{\mathcal{Z}}_{\varepsilon}} \int \frac{1}{y-2} \mathrm{P_{\varepsilon}}(\mathrm d y)
\end{split}$$ and analogously $$\begin{split}
{\mathbb E}Y' & = 2 + \frac 1{{\mathcal{Z}}_{\varepsilon}} \left(\int \frac{y'}{(y'-1)(y'-2)} \mathrm{P}'_{{\varepsilon}}(\mathrm d y') - 2\Delta_{{\varepsilon}'}\right)\\
& = 2 + \frac 1{{\mathcal{Z}}_{\varepsilon}} \int \frac{1}{y'-1} \mathrm{P}'_{{\varepsilon}}(\mathrm d y')\,,
\end{split}$$ which implies $$\begin{split}
{\mathbb E}Y' - {\mathbb E}Y & = 1 - \frac 1{{\mathcal{Z}}_{\varepsilon}} \left(\int \frac{1}{y - 2} {\, \mathrm{d}}\mathrm{P}_{\varepsilon}(y) - \int \frac{1}{y' - 1} {\, \mathrm{d}}\mathrm{P'}_{\varepsilon}(y')\right) \\
& = 1 - \frac 1{{\mathcal{Z}}_{\varepsilon}} ({\mathcal{Z}}_{\varepsilon}) = 0\,.
\end{split}$$

To verify the inequality, note that $$\int \frac{1}{y-2} \mathrm{P_{\varepsilon}}(\mathrm d y) = {\mathbb E}\frac{1}{{\varepsilon}^{-1}X - 2} = {\varepsilon}{\mathbb E}\frac{1}{X - 2{\varepsilon}}\,,$$ and the fact that $X \geq 1$ almost surely and ${\varepsilon}\leq 1/6$ implies that this quantity is bounded by $\frac 32 {\varepsilon}$. By \[rescaled-prior-bounds\], ${\mathcal{Z}}_{\varepsilon}\geq \frac{3}{10} {\varepsilon}$; therefore $${\mathbb E}Y = 1 + \frac 1{{\mathcal{Z}}_{\varepsilon}} \int \frac{1}{y-2} \mathrm{P_{\varepsilon}}(\mathrm d y) \leq 1 + \frac{\frac 32 {\varepsilon}}{\frac{3}{10} {\varepsilon}} = 6\,.$$

\[approximate-probabilities\] If $Q = \frac 1m (U_1, \dots, U_m)$ and $Q' = \frac 1m (V_1, \dots, V_m)$, where $U$ and $V$ satisfy the assumptions of \[priors\] with ${\varepsilon}= \frac 1 {24} \delta$, and recall that $E = \{Q \in \tilde {\mathcal{D}}_{m, \delta}^-\}$ and $E' = \{Q \in \tilde {\mathcal{D}}_{m}^+\}$. Then $${\mathbb P}[E^C]+ {\mathbb P}[{E'}^C] \lesssim \frac{L^4}{\delta^2 m}\,.$$

We first show that $E$ holds as long as $\frac 1 m \sum_{i=1}^m U_i^2 \leq 8$ and $\frac 1 m \sum_{i=1}^m |U_i - 1| \leq \delta$. Indeed, if we set $s {:=}\sum_{i=1}^m Q_i$, where $Q_i=U_i/m$, then the second condition together with Jensen’s inequality imply that $|s - 1| \leq \delta$. Therefore, recalling that $\bar Q=Q/s$, we get $$\chi^2(\bar Q, u) + 1 = m \sum_{i=1}^m {\bar Q_i}^2 = \frac{1}{s^2m} \sum_{i=1}^m U_i^2 \leq \frac{8}{s^2} \leq \frac{8}{1-2 \delta} \leq 10\,,$$ where the final inequality follows from the assumption that $\delta \leq 1/10$. Hence $\chi^2(\bar Q, u) \leq 9$.

Likewise, $${d_{\mathrm{TV}}(\bar Q, u)} = \frac 1 2 \sum_{i=1}^m \left|\bar{Q_i} - \frac 1 m\right| \leq \frac 12 \sum_{i=1}^m \left|\bar{Q_i} - Q_i\right| + \frac 1 2 \sum_{i=1}^m \left|Q_i - \frac 1 m \right| \leq \frac 12 (|s-1|+ \delta) \leq \delta\,.$$

By Chebyshev’s inequality combined with the assumption that ${\mathbb E}U^2 \leq 6$, $${\mathbb P}\left[\frac 1 m \sum_{i=1}^m U_i^2 > 8 \right] \leq \frac{\operatorname{Var}(U^2)}{4 m}\,,$$ and since $U \leq 16 {\varepsilon}^{-1} L^2$ almost surely, $$\operatorname{Var}(U^2) \leq {\mathbb E}U^4 \leq (16 {\varepsilon}^{-1} L^2)^2 {\mathbb E}U^2 \lesssim \delta^{-2} L^4\,.$$ Moreover, since ${\varepsilon}= \frac 1{24} \delta$ and ${\mathbb E}|U - 1| \leq 12 {\varepsilon}$, $${\mathbb P}\left[\frac 1m \sum_{i=1}^m |U_i - 1| > \delta\right] \leq \frac{4 \operatorname{Var}(|U-1|)}{\delta^2 m} \lesssim \frac{1}{\delta^2 m}\,.$$ Combining these bounds yields $${\mathbb P}[E^C] \lesssim \frac{L^4}{\delta^2 m}\,.$$

The argument for $E'$ is analogous. It suffices for $E'$ to hold that $\frac 1m \sum_{i=1}^m V_i^2 \leq 8$ along with the conditions $\frac 1 m \sum_{i=1}^m |V_i - 1| \geq \frac 3 4$ and $\left|(\frac 1m \sum_{i=1}^m V_i) - 1\right| \leq \delta$. As above, the first condition is violated with probability at most $C \frac{L^4}{\delta^2 m}$, and the latter two conditions are violated with probability at most $C \frac{1}{\delta^2 m}$. The claim follows.

Omitted proofs {#subsec:omitted}
--------------

Denote by $\rho^{(1)}$ and $\rho^{(2)}$ the distributions of $Z^{(1)}$ and $Z^{(2)}$, respectively. It suffices to show ${\tilde W_{p,k}}({\nu^{(1)}}, {\nu^{(2)}}) + W_p(\rho^{(1)}, \rho^{(2)}) \geq W_p({\nu^{(1)}}, {\nu^{(2)}})$. Let $U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}$ be such that the column span of $U$ agrees with the subspace ${\mathcal{U}}$ on which the random variables $X^{(1)}$ and $X^{(2)}$ are supported in the definition of ${\mu^{(1)}}$ and ${\mu^{(2)}}$. We let $(Y^{(1)}, Y^{(2)})$ be a pair of random variables such that marginally $Y^{(i)} \sim {\mu^{(i)}}_U$ for $i \in \{1, 2\}$ and $$W^p_p({\mu^{(1)}}_U, {\mu^{(2)}}_U) = {\mathbb E}\|Y^{(1)} - Y^{(2)}\|^p\,.$$ Likewise, let $(W^{(1)}, W^{(2)})$ be a coupling of $\rho^{(1)}$ and $\rho^{(2)}$ such that $$W^p_p(\rho^{(1)}, \rho^{(2)}) = {\mathbb E}\|W^{(1)} - W^{(2)}\|^p\,.$$ Then  implies that $(Y^{(1)} + W^{(1)}, Y^{(2)} + W^{(2)})$ is a coupling of ${\mu^{(1)}}$ and ${\mu^{(2)}}$. Therefore $$\begin{aligned}
W_p({\mu^{(1)}}, {\mu^{(2)}}) & \leq \left({\mathbb E}\|Y^{(1)} + W^{(1)} - (Y^{(2)} + W^{(2)})\|^p\right)^{1/p} \\
& \leq \left({\mathbb E}\|Y^{(1)} - Y^{(2)}\|^p\right)^{1/p} + \left({\mathbb E}\|W^{(1)} - W^{(2)}\|^p\right)^{1/p} \\
& = W_p({\mu^{(1)}}_U, {\mu^{(2)}}_U) + W_p(\rho^{(1)}, \rho^{(2)})\\
& \leq {\tilde W_{p,k}}({\mu^{(1)}}, {\mu^{(2)}}) + W_p(\rho^{(1)}, \rho^{(2)})\,.\end{aligned}$$

This can be deduced from a series of well known facts. By \[prop:bobgot\], it suffices to show that for $f$ Lipschitz, $${\mathbb E}e^{\lambda (f(X) - {\mathbb E}f(X))} \leq e^{C^2 k \lambda^2 \sigma^2/2}\,.$$ By symmetrization, we can bound $${\mathbb E}e^{\lambda (f(X) - {\mathbb E}f(X))} \leq {\mathbb E}e^{\lambda {\varepsilon}(f(X) - f(X'))} \leq {\mathbb E}e^{\lambda {\varepsilon}\|X - X'\|}\,,$$ where $X' \sim \mu$ is independent of $X$. Continuing, we obtain $${\mathbb E}e^{\lambda {\varepsilon}\|X - X'\|} \leq {\mathbb E}e^{2 \lambda {\varepsilon}\|X\|} \leq e^{16 k \lambda^2 \sigma^2}$$ by \[prop:norm\_subG\]. This proves the claim with $C = 32$.

First, assume that $\mu$ satisfies the ${{T_{p}(\sigma^2)}}$ inequality. Then, by @GozLeo10 [Proposition 1.9], the product measure $\mu^{\otimes n}$ satisfies the inequality ${{T_{p}(n^{2/p-1}\sigma^2 )}}$ on the space ${\mathcal{X}}^n$ equipped with the metric $\rho_p(x, y) {:=}\left(\sum_{i=1}^n \rho(x_i, y_i)^p \right)^{1/p}$. Therefore, $\mu^{\otimes n}$ also satisfies ${{T_{1}(n^{2/p - 1}\sigma^2 )}}$ with respect to this same metric.

We now note that the function $$(x_1, \dots x_n) \mapsto W_p\left(\frac 1n \sum_{i=1}^n \delta_{x_i}, \mu \right)$$ is $n^{-1/p}$-Lipschitz with respect to $\rho_p$. Indeed, we have $$\begin{aligned}
\Big|W_p\Big(\frac 1n \sum_{i=1}^n \delta_{x_i}, \mu \Big) - W_p\Big(\frac 1n \sum_{i=1}^n \delta_{y_i}, \mu \Big)\Big| & \leq W_p\Big(\frac 1n \sum_{i=1}^n \delta_{x_i},\frac 1n \sum_{i=1}^n \delta_{y_i}\Big) \\
& \leq \Big(\frac 1n \sum_{i=1}^n \rho(x_i, y_i)^p \Big)^{1/p}\\
& = n^{-1/p} \rho_p(x, y)\,.\end{aligned}$$

Combining these observations with \[prop:bobgot\], we obtain $${\mathbb E}e^{\lambda(W_p(\mu_n, \mu) - {\mathbb E}W_p(\mu_n, \mu))} \leq e^{\lambda^2 n^{-2/p} n^{2/p-1} \sigma^2/2} = e^{\lambda^2 n^{-1} \sigma^2/2} \qquad \forall \lambda \in {\mathbb{R}}\,.$$

In the other direction, suppose that $W_p(\mu_n, \mu)$ is $\sigma^2/n$-subgaussian. A Chernoff bound implies the concentration inequality $${\mathbb P}[W_p(\mu_n, \mu) - {\mathbb E}W_p(\mu_n, \mu) \geq t] \leq e^{nt^2/2\sigma^2}\,.$$ We follow the proof of @Goz09 [Theorem 3.4], and only sketch the argument here. It is easy to verify that the function $W_p(\cdot, \mu)$ is lower semi-continuous with respect to the weak topology [see @Vil09 proof of Theorem 4.1], which implies that, for any $t \geq 0$, the set ${\mathcal{O}}_t {:=}\{\nu \in {\mathcal{P}}({\mathcal{X}}) : W_p(\nu, \mu) > t\}$ is open. Moreover, as $\mu \in {\mathcal{P}}_p$, @Vil09 [Theorem 6.9] combined with Varadarajan’s theorem yields that ${\mathbb E}W_p(\mu_n, \mu) \to 0$ as $n \to \infty$ [see @Vil09 proof of Theorem 22.22].

Applying Sanov’s theorem then yields $$\begin{aligned}
- \inf_{\nu \in {\mathcal{O}}_t} D(\nu \| \mu) & \leq \liminf_{n \to \infty} \frac 1n \log {\mathbb P}[W_p(\mu_n, \mu) > t]  \\
& \leq \liminf_{n \to \infty} \frac 1n \cdot \frac{n \max\{(t- {\mathbb E}W_p(\mu_n, \mu))^2, 0\}}{2 \sigma^2} \\
& = \frac{t^2}{2 \sigma^2}\,.\end{aligned}$$ We obtain that, for any $t \geq 0$, $$D(\nu \| \mu) \geq - \frac{t^2}{2 \sigma^2} \qquad \forall \nu \in {\mathcal{P}}({\mathcal{X}}) \text{ s.t. } W_p(\nu, \mu) > t\,.$$ Therefore $\mu$ satisfies ${{T_{p}(\sigma^2)}}$.

As in the proof of \[thm:concentration\_is\_tp\], we have that $\mu^{\otimes n}$ satisfies ${{T_{1}(n^{2/p'-1} \sigma^2)}}$ with respect to $\rho_{p'}(x, y)$. We also that that $$(x_1, \dots x_n) \mapsto W_p\left(\frac 1n \sum_{i=1}^n \delta_{x_i}, \mu \right)$$ is $n^{-1/p}$-Lipschitz with respect to $\rho_{p'}$, since the fact that $p' \leq p$ implies $$\Big|W_p\Big(\frac 1n \sum_{i=1}^n \delta_{x_i}, \mu \Big) - W_p\Big(\frac 1n \sum_{i=1}^n \delta_{y_i}, \mu \Big)\Big| \leq n^{-1/p} \rho_p(x, y) \leq n^{-1/p} \rho_{p'}(x, y)\,.$$ We conclude as in the first part of the proof of \[thm:concentration\_is\_tp\].

We show that for any measures ${\nu^{(1)}}$, ${\nu^{(2)}}$, if ${\mathcal{V}}= \mathrm{span}(\hat V)$ where $V = \operatorname*{argmax}_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} W_p({\nu^{(1)}}_U, {\nu^{(2)}}_U)$, then $$\sin^2\big(\measuredangle({\mathcal{V}}, {\mathcal{U}})\big) \lesssim \frac{{\tilde W_{p,k}}({\mu^{(1)}}, {\nu^{(1)}}) + {\tilde W_{p,k}}({\mu^{(2)}}, {\nu^{(2)}})}{W_p({\mu^{(1)}}, {\mu^{(2)}})}\,.$$ Combined with \[prop:supremum\], this implies the claim.

guarantees the existence of an optimal coupling between ${\mu^{(1)}}$ and ${\mu^{(2)}}$ such that if $(X, Y) \sim \gamma$, then $X - Y$ lies in ${\mathcal{U}}$ almost surely. Let $U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}$ be such that the column span of $U$ is ${\mathcal{U}}$. We have $$\begin{aligned}
W_p({\mu^{(1)}}_{V}, {\mu^{(2)}}_{V}) & \leq ({\mathbb E}\|V (X - Y) \|^p)^{1/p} \\
& = ({\mathbb E}\|V U^\top U (X - Y) \|^p)^{1/p} \\
& \leq {\left\|V U^\top\right\|_{\mathrm{op}}} ({\mathbb E}\|X - Y\|^p)^{1/p} = {\left\|V U^\top\right\|_{\mathrm{op}}} W_p({\mu^{(1)}}, {\mu^{(2)}})\,.\end{aligned}$$

By the definition of $V$, $$\begin{aligned}
W_p({\mu^{(1)}}, {\mu^{(2)}}) - W_p({\mu^{(1)}}_V, {\mu^{(2)}}_V) &  = W_p({\mu^{(1)}}_U, {\mu^{(2)}}_U) - W_p({\mu^{(1)}}_V, {\mu^{(2)}}_V) \\
& \leq W_p({\mu^{(1)}}_U, {\mu^{(2)}}_U) - W_p({\nu^{(1)}}_U, {\nu^{(2)}}_U) \\
& \quad \quad \quad + W_p({\nu^{(1)}}_V, {\nu^{(2)}}_V) - W_p({\mu^{(1)}}_V, {\mu^{(2)}}_V) \\
& \le 2 \sup_{U \in {{\mathcal{V}}_k({\mathbb{R}}^d)}} |W_p({\mu^{(1)}}_U, {\mu^{(2)}}_U) - W_p({\nu^{(1)}}_U, {\nu^{(2)}}_U)| \\
& \leq 2{\tilde W_{p,k}}({\mu^{(1)}}, {\nu^{(1)}}) + 2{\tilde W_{p,k}}({\mu^{(2)}}, {\nu^{(2)}})\,.\end{aligned}$$ Therefore $$(1 - {\left\|V U^\top\right\|_{\mathrm{op}}}) W_p({\mu^{(1)}}, {\mu^{(2)}}) \le 2 {\tilde W_{p,k}}({\mu^{(1)}}, {\nu^{(1)}}) + 2{\tilde W_{p,k}}({\mu^{(2)}}, {\nu^{(2)}})\,.$$

Moreover, $$\begin{aligned}
\sin^2 \big(\measuredangle({\mathcal{V}}, {\mathcal{U}})\big) & = 1 - \sup_{\substack{v \in {\mathcal{V}}, u \in {\mathcal{U}}\\ \|u\| = \|v\| = 1}}(u^\top v)^2 \\
& = 1 - \sup_{\substack{x, y \in {\mathbb{R}}^d \\ \|x\| = \|y\| = 1}} ((Ux)^\top (Vy))^2 \\
& = 1 - {\left\|V U^\top\right\|_{\mathrm{op}}}^2 \\
& \leq 2(1 - {\left\|V U^\top\right\|_{\mathrm{op}}})\,,\end{aligned}$$ since ${\left\|V U^\top\right\|_{\mathrm{op}}} \leq 1$. The claim follows.

Acknowledgements {#acknowledgements .unnumbered}
================

The authors thank Oded Regev for helpful suggestions.
