Entropic estimation of optimal transport maps

Aram-Alexandre Pooladian∗, Jonathan Niles-Weed∗†
---
∗Center for Data Science, New York University
---
†Courant Institute of Mathematical Sciences, New York University
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

We develop a computationally tractable method for estimating the optimal map
between two distributions over $\mathbb{R}^{d}$ with rigorous finite-sample
guarantees. Leveraging an entropic version of Brenier’s theorem, we show that
our estimator—the _barycentric projection_ of the optimal entropic plan—is
easy to compute using Sinkhorn’s algorithm. As a result, unlike current
approaches for map estimation, which are slow to evaluate when the dimension
or number of samples is large, our approach is parallelizable and extremely
efficient even for massive data sets. Under smoothness assumptions on the
optimal map, we show that our estimator enjoys comparable statistical
performance to other estimators in the literature, but with much lower
computational cost. We showcase the efficacy of our proposed estimator through
numerical examples. Our proofs are based on a modified duality principle for
entropic optimal transport and on a method for approximating optimal entropic
plans due to Pal (2019).

##  1 Introduction

The goal of optimal transport is to find a map between two probability
distributions that minimizes the squared Euclidean transportation cost. This
formulation leads to what is known as the Monge problem (Monge, 1781):

$\displaystyle\min_{T\in\mathcal{T}(P,Q)}\int\frac{1}{2}\|x-T(x)\|^{2}_{2}\,\mathrm{d}P(x),$
(1)

where $P$ and $Q$ are two probability measures on
$\Omega\subseteq\mathbb{R}^{d}$, and $\mathcal{T}(P,Q)$ is the set of
admissible maps111$\mathcal{T}(P,Q):=\\{T:\Omega\to\Omega\ |\
T_{\sharp}P:=P\circ T^{-1}=Q\\}$. A solution to the Monge problem is
guaranteed to exist if $P$ and $Q$ have finite second moments and $P$ is
absolutely continuous; moreover, the optimal map enjoys certain regularity
properties under stricter assumptions on $P$ and $Q$ (see Section 2.1 for more
information). Due to their versatility and mathematical simplicity, optimal
transport maps have found a wide range of uses in statistics and machine
learning, (Carlier et al., 2016; Chernozhukov et al., 2017; Arjovsky et al.,
2017; Huang et al., 2021; Finlay et al., 2020; Courty et al., 2014, 2017; Wang
et al., 2010; Onken et al., 2021; Makkuva et al., 2020), computer graphics
(Solomon et al., 2016, 2015; Feydy et al., 2017), and computational biology
(Schiebinger et al., 2019; Yang et al., 2020), among other fields.

In many applications, we are not given access to the full probability measures
$P$ and $Q$, but independent samples from them, denoted
$X_{1},\ldots,X_{n}\sim P$ and $Y_{1},\ldots,Y_{n}\sim Q$. When an optimal map
$T_{0}\in\mathcal{T}(P,Q)$ minimizing Eq. 1 exists, it is natural to ask
whether it is possible to estimate $T_{0}$ on the basis of these samples.
Hütter and Rigollet (2021) investigated this question and proposed an
estimator $\hat{T}_{n}$ which achieves

$\displaystyle\mathbb{E}\|\hat{T}_{n}-T_{0}\|^{2}_{L^{2}(P)}\lesssim
n^{-\frac{2\alpha}{2\alpha-2+d}}\log^{3}(n)\,,$ (2)

if $T_{0}\in\mathcal{C}^{\alpha}$ and $P$ and $Q$ satisfy additional technical
assumptions. Moreover, they showed that the rate in Eq. 2 is minimax optimal
up to logarithmic factors. Though statistically optimal, their estimator is
impractical to compute if $d>3$, since it relies on a gridding scheme whose
computational cost scales exponentially in the dimension. Recently, Deb et al.
(2021) and Manole et al. (2021) proposed plugin estimators that also achieve
the minimax estimation rate. Though simpler to compute than the estimator of
Hütter and Rigollet (2021), these estimators require at least $O(n^{3})$ time
to compute and cannot easily be parallelized, making them an unfavorable
choice when the number of samples is large.

We adopt a different approach by leveraging recent advances in computational
optimal transport based on entropic regularization (Peyré and Cuturi, 2019),
which replaces Eq. 1 by

$\inf_{\pi\in\Pi(P,Q)}\iint\frac{1}{2}\|x-y\|^{2}\,\mathrm{d}\pi(x,y)+\varepsilon
D_{\text{KL}}(\pi\|P\otimes Q)\,,$ (3)

where $\Pi(P,Q)$ denotes the set of couplings between $P$ and $Q$ and
$D_{\text{KL}}(\cdot\|\cdot)$ denotes the Kullback–Leibler divergence. This
approach, which was popularized by Cuturi (2013), has been instrumental in the
adoption of optimal transport methods in the machine learning community
because it leads to a problem that can be solved by Sinkhorn’s algorithm
(Sinkhorn, 1967), whose time complexity scales _quadratically_ in the number
of samples (Altschuler et al., 2017). Moreover, Sinkhorn’s algorithm is
amenable to parallel implementation on GPUs, making it very attractive for
large-scale problems (Feydy et al., 2019, 2020; Genevay et al., 2016, 2018;
Altschuler et al., 2019).

The efficiency and popularity of Sinkhorn’s algorithm raise the tantalizing
question of whether it is possible to use this practical technique to develop
estimators of optimal transport maps with convergence guarantees. In this
work, we develop such a procedure. Under suitable technical assumptions on $P$
and $Q$, we show that our estimator $\hat{T}$ enjoys the rate

$\mathbb{E}\|\hat{T}-T_{0}\|_{L_{2}(P)}^{2}\lesssim
n^{-\frac{(\alpha+1)}{2(d^{\prime}+\alpha+1)}}\log n$

for $d^{\prime}=2\lceil d/2\rceil$ if the inverse map $T_{0}^{-1}$ is
$\mathcal{C}^{\alpha}$ and $\alpha\in(1,3]$. This rate is worse than Eq. 2,
but our empirical results show that our estimator nevertheless outperforms all
other estimators proposed in the literature in terms of _both_ computational
and statistical performance. The estimator we analyze was originally suggested
by Seguy et al. (2018), who also showed consistency of the entropic plan in
the large-$n$ limit if the regularization parameter is taken to zero
sufficiently fast. However, to our knowledge, our work offers the first
finite-sample convergence guarantees for this proposal.

Our estimator is defined as the barycentric projection (Ambrosio et al., 2008)
of the entropically optimal coupling between the empirical measures arising
from the samples. This barycentric projection has been leveraged in other
works on map estimation as a straightforward way of obtaining a function from
a coupling between two probability measures (Deb et al., 2021). However, in
the context of entropic optimal transport, this operation has a more canonical
interpretation in light of Brenier’s theorem (Brenier, 1991). Brenier’s result
says that the solution to Eq. 1 can be realized as the gradient of the
function which solves the dual problem to Eq. 1; we show in Proposition 2 that
the barycentric projection of the entropically optimal coupling is the
gradient of the function which solves the dual problem to Eq. 3. In addition
to providing a connection to the classical theory of optimal transport, this
observation provides a canonical extension $\hat{T}$ to out-of-sample points.
Moreover, since Sinkhorn’s algorithm computes solutions to the dual of Eq. 3,
this interpretation shows that computing $\hat{T}$ is no more costly than
solving Eq. 3.

We analyze $\hat{T}$ by employing a strategy pioneered by Pal (2019) for
understanding the structure of the optimal entropic coupling. This technique
compares the solution to Eq. 3 to a coupling whose conditional laws are
Gaussian, with mean and covariance given by the solution to Eq. 1. To leverage
this comparison, we employ a duality principle in conjuction with short-time
expansions of the value of Eq. 3 developed in Conforti and Tamanini (2021) and
Chizat et al. (2020).

#### Paper outline

The paper is organized as follows: Section 2 reviews the relevant background
on optimal transport theory and entropic regularization for the quadratic
cost. We define our estimator and preview our main results in Section 3. Our
main statistical bounds appear in Section 4 and Section 5. A discussion of
computational considerations and numerical experiments are provided in Section
6.

#### Notation

The support of a probability distribution is given by
$\operatorname{\mathrm{supp}}(\cdot)$. For a convex function $\varphi$, we
denote its convex dual by $\varphi^{*}(y)=\sup_{x}x^{\top}y-\varphi(x).$ The
Kullback–Liebler divergence between two measures is denoted by
$D_{\text{KL}}(\mu\|\nu)=\int\log(\frac{\,\mathrm{d}\mu}{\,\mathrm{d}\nu})\,\mathrm{d}\mu.$
If $P$ possesses a density $p$ with respect to the Lebesgue measure, we denote
its differential entropy by $\text{Ent}(P)=\int p(x)\log(p(x))\,\mathrm{d}x$.
For a joint probability density $p(x,y)$, we denote the conditional density of
$y$ as $p^{x}(y)$. The square-root of the determinant of a matrix is
$J(\cdot):=\sqrt{\det(\cdot)}$. For $\alpha\geq 0$ and a closed set $\Omega$,
we write $h\in\mathcal{C}^{\alpha}(\Omega)$ if there exists an open set
$U\supseteq\Omega$ and a function $g:U\to\mathbb{R}$ such that $g|_{\Omega}=h$
and such that $g$ possesses $\lfloor\alpha\rfloor$ continuous derivatives and
whose $\lfloor\alpha\rfloor$th derivative is $(\alpha-\lfloor\alpha\rfloor)$
Hölder smooth. The total variation distance between two probability measures
$\mu$ and $\nu$ is $\mathrm{d_{TV}}(\mu,\nu)=\sup_{f:\|f\|_{\infty}\leq 1}\int
f(\mathrm{d}\mu-\mathrm{d}\nu)$. For $a\in\mathbb{R}^{d}$ and $r>0$, we write
$B_{r}(a)$ for the Euclidean ball of radius $r$ centered at $a$. A constant is
a quantity whose value may depend on the smoothness parameters appearing in
assumptions (A1) to (A3), the set $\Omega$, and the dimension, but on no other
quantities. We denote the maximum and minimum of $a$ and $b$ by $a\vee b$ and
$a\wedge b$, respectively. We use the symbols $c$ and $C$ to denote positive
constants whose value may change from line to line, and write $a\lesssim b$
and $a\asymp b$ if there exists constants $c,C>0$ such that $a\leq Cb$ and
$cb\leq a\leq Cb$, respectively.

##  2 Background on optimal transport theory

###  2.1 Optimal transport under the quadratic cost

Let $\mathcal{P}(\Omega)$ be the space of (Borel) probability measures with
support contained in $\Omega$, and $\mathcal{P}_{\text{ac}}(\Omega)$ be those
with densities. We first present Brenier’s Theorem, which guarantees the
existence of an optimal map between two distributions when the first measure
is absolutely continuous.

######  Theorem 1 (Brenier’s Theorem).

(Brenier, 1991) Let $P\in\mathcal{P}_{\text{ac}}(\Omega)$ and
$Q\in\mathcal{P}(\Omega)$. Then

  1. 1.

there exists a solution $T_{0}$ to Eq. 1, with $T_{0}=\nabla\varphi_{0}$, for
a convex function $\varphi_{0}$ solving

$\displaystyle\inf_{\varphi\in
L^{1}(P)}\int\varphi\,\mathrm{d}P+\int\varphi^{*}\,\mathrm{d}Q,$ (4)

where $\varphi^{*}$ is the convex conjugate to $\varphi$.

  2. 2.

If in addition $Q\in\mathcal{P}_{\text{ac}}(\Omega)$, then
$\nabla\varphi^{*}_{0}$ is the optimal transport map from $Q$ to $P$.

If $P$ does not have a density, then Eq. 1 may not have a solution, but this
problem can be ameliorated by passing to a convex relaxation of Eq. 1 due to
Kantorovitch (1942), which leads to the definition of the 2-Wasserstein
distance:

$\displaystyle\frac{1}{2}W_{2}^{2}(P,Q):=\min_{\pi\in\Pi(P,Q)}\int\frac{1}{2}\|x-y\|_{2}^{2}\,\mathrm{d}\pi(x,y),$
(5)

where

$\Pi(P,Q):=\\{\pi\in\mathcal{P}(\Omega\times\Omega)\ |\
\pi(A\times\Omega)=P(A),\pi(\Omega\times A)=Q(A)\\}$

is the set of couplings between $P$ and $Q$. Unlike Eq. 1, Eq. 5 always admits
a minimizer when $P$ and $Q$ have finite second moments. We call such a
minimizer an optimal plan, denoted $\pi_{0}$.

Equation 5 also possesses a dual formulation (see Villani, 2008; Santambrogio,
2015):

$\displaystyle\frac{1}{2}W^{2}_{2}(P,Q)=\sup_{(f,g)\in\Phi}\int
f\,\mathrm{d}P+\int g\,\mathrm{d}Q,$ (6)

where

$\Phi=\left\\{(f,g)\in L^{1}(P)\times L^{1}(Q)\ |\
f(x)+g(y)\leq\frac{1}{2}\|x-y\|_{2}^{2}\ \text{for all
}x,y\in\Omega\right\\}.$

Once again, if $P$ and $Q$ have finite second moments, then the supremum in
Eq. 6 is always achieved by a pair $(f_{0},g_{0})\in\Phi$, called optimal
potentials.

######  Remark 1.

It is straightforward to see that Eq. 8 and Eq. 4 are in fact equivalent: if
$(f_{0},g_{0})$ are solutions to Eq. 6, then
$\varphi_{0}=\frac{1}{2}\|\cdot\|_{2}^{2}-f_{0}$ solves Eq. 4 and, if
$P\in\mathcal{P}_{ac}$, $T_{0}=(\text{Id}-\nabla f_{0})$ solves Eq. 1.

###  2.2 Entropic optimal transport under the quadratic cost

Entropic optimal transport is defined by adding an entropic regularization
term to Eq. 5. Letting $P,Q\in\mathcal{P}(\Omega)$ and a regularization
parameter $\varepsilon>0$, the entropically regularized problem is

$\displaystyle S_{\varepsilon}(P,Q)$
$\displaystyle:=\inf_{\pi\in\Pi(P,Q)}\iint\frac{1}{2}\|x-y\|_{2}^{2}\,\mathrm{d}\pi(x,y)+\varepsilon
D_{\text{KL}}(\pi\|P\otimes Q)\,,$ (7) which, when $P$ and $Q$ have densities
$p$ and $q$ respectively, can also be written $\displaystyle
S_{\varepsilon}(P,Q)$
$\displaystyle=\inf_{\pi\in\Pi(P,Q)}\iint\frac{1}{2}\|x-y\|_{2}^{2}\,\mathrm{d}\pi(x,y)+\varepsilon\iint\log(\pi)\,\mathrm{d}\pi-\varepsilon\iint\log(p(x)q(y))\,\mathrm{d}\pi(x,y)$
$\displaystyle=\inf_{\pi\in\Pi(P,Q)}\iint\frac{1}{2}\|x-y\|_{2}^{2}\,\mathrm{d}\pi(x,y)+\varepsilon\iint\log(\pi)\,\mathrm{d}\pi-\varepsilon(\operatorname{Ent}{P}+\operatorname{Ent}{Q})\,.$

This problem also admits a dual (see Genevay, 2019), which is a relaxed
version of Eq. 6:

$\displaystyle S_{\varepsilon}(P,Q)=\sup_{\begin{subarray}{c}f\in L^{1}(P)\\\
g\in L^{1}(Q)\end{subarray}}\int f\,\mathrm{d}P+\int
g\,\mathrm{d}Q-\varepsilon\iint
e^{(f(x)+g(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}P(x)\,\mathrm{d}Q(y)+\varepsilon\,.$
(8)

Both Eq. 7 and Eq. 8 possess solutions if $P$ and $Q$ have finite second
moments; moreover, if we denote by $\pi_{\varepsilon}$ the solution to Eq. 7,
which we call the _optimal entropic plan_ , and
$(f_{\varepsilon},g_{\varepsilon})$ the solution to Eq. 8, which we call then
optimal entropic potentials, then we obtain the optimality relation (Csiszár,
1975):

$\,\mathrm{d}\pi_{\varepsilon}(x,y):=\tilde{\pi}_{\varepsilon}(x,y)\,\mathrm{d}P(x)\,\mathrm{d}Q(y):=\exp((f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon)\,\mathrm{d}P(x)\,\mathrm{d}Q(y)\,.$

A consequence of this relation is that we may choose optimal entropic
potentials satisfying

$\begin{split}\int
e^{\frac{1}{\varepsilon}(f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}P(x)&=1\quad\forall
y\in\mathbb{R}^{d}\\\ \int
e^{\frac{1}{\varepsilon}(f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q(y)&=1\quad\forall
x\in\mathbb{R}^{d}\,.\end{split}$ (9)

We will therefore always assume that Eq. 9 holds. Conversely, if a pair of
functions $(f_{\varepsilon},g_{\varepsilon})$ satisfies the duality conditions
Eq. 9 $P\otimes Q$ almost everywhere, then
$\exp((f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon)$
is the $P\otimes Q$ density of the optimal entropic plan.

Our proofs rely on a modified version of the duality relation given in Eq. 8,
in which the supremum is taken over a larger set of functions. Though it is a
straightforward consequence of Fenchel’s inequality, we have not encountered
this statement explicitly in the literature, so we highlight it here.

######  Proposition 1.

Assume $P$ and $Q$ possess finite second moments, and let $\pi_{\varepsilon}$
be the optimal entropic plan for $P$ and $Q$. Then

$S_{\varepsilon}(P,Q)=\sup_{\eta\in
L^{1}(\pi_{\varepsilon})}\int\eta\,\mathrm{d}\pi_{\varepsilon}-\varepsilon\iint
e^{(\eta(x,y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}P(x)\,\mathrm{d}Q(y)+\varepsilon\,.$
(10)

Comparing this proposition with Eq. 8, we see that we can always take
$\eta(x,y)=f(x)+g(y)$, in which case Eq. 10 reduces to Eq. 8. The novelty in
Proposition 1 therefore arises in showing that the quantity on the right side
of Eq. 10 is still bounded above by $S_{\varepsilon}(P,Q)$. We give the short
proof of Proposition 1 in Appendix C.

Several recent works have bridged the regularized and unregularized optimal
transport regimes, with particular interest in the setting where
$\varepsilon\to 0$. Convergence of $\pi_{\varepsilon}$ to $\pi_{0}$ was
studied by Carlier et al. (2017) and Léonard (2012), and recent work has
quantified the convergence of the plans (Ghosal et al., 2021; Bernton et al.,
2021) and the potentials (Nutz and Wiesel, 2021; Altschuler et al., 2021) in
certain settings. Convergence of $S_{\varepsilon}(P,Q)$ to
$\frac{1}{2}W_{2}^{2}(P,Q)$ has attracted significant research interest: under
mild conditions, Pal (2019) proves a first-order convergence result for
general convex costs (replacing $\frac{1}{2}\|\cdot\|^{2}$), and a second
order expansion was subsequently obtained by Chizat et al. (2020) and Conforti
and Tamanini (2021). We rely on the following bound, derived from Chizat et
al. (2020).

######  Theorem 2.

Suppose $P$ and $Q$ have bounded densities with compact support. Then

$\displaystyle
S_{\varepsilon}(P,Q)-\frac{1}{2}W^{2}_{2}(P,Q)+\varepsilon\log((2\pi\varepsilon)^{d/2})\leq-\frac{\varepsilon}{2}\left(\operatorname{Ent}(P)+\operatorname{Ent}(Q)\right)+\frac{\varepsilon^{2}}{8}I_{0}(P,Q),$
(11)

where $I_{0}(P,Q)$ is the integrated Fisher information along the Wasserstein
geodesic between $P$ and $Q$.

For completeness, a brief proof appears in Appendix A.

##  3 Estimator and main results

Given the optimal entropic plan $\pi_{\varepsilon}$ between $P$ and $Q$, we
define its barycentric projection to be

$T_{\varepsilon}(x):=\int
y\,\mathrm{d}\pi_{\varepsilon}^{x}(y)=\mathbb{E}_{\pi_{\varepsilon}}[Y\mid
X=x]\,.$ (12)

A priori, this map is only defined $P$-almost everywhere, making it unsuitable
for evaluation outside the support of $P$. In particular, since we will study
the barycentric projection obtained from the optimal entropic plan between
empirical measures, this definition does not extend outside the sample points.
However, the duality relation Eq. 9 implies that we may define a version of
the conditional density of $Y$ given $X=x$ for _all_ $x\in\mathbb{R}^{d}$ by

$\,\mathrm{d}\pi_{\varepsilon}^{x}(y)=e^{\frac{1}{\varepsilon}(f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q(y)=\frac{e^{\frac{1}{\varepsilon}(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q(y)}{\int
e^{\frac{1}{\varepsilon}(g_{\varepsilon}(y^{\prime})-\frac{1}{2}\|x-y^{\prime}\|^{2})}\,\mathrm{d}Q(y^{\prime})}\,,$

where $(f_{\varepsilon},g_{\varepsilon})$ are the entropic potentials. This
furnishes an extension of $T_{\varepsilon}$ to all of $\mathbb{R}^{d}$ by

$T_{\varepsilon}(x)\coloneqq\frac{\int
ye^{\frac{1}{\varepsilon}(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q(y)}{\int
e^{\frac{1}{\varepsilon}(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q(y)}\,.$
(13)

We call $T_{\varepsilon}$ the _entropic map_ between $P$ and $Q$, though we
stress that $(T_{\varepsilon})_{\sharp}P\neq Q$ in general. This natural
definition is motivated by the following observation, which shows that the
entropic map can also be defined as the map obtained by replacing the optimal
potential in Brenier’s theorem by its entropic counterpart.

######  Proposition 2.

Let $(f_{\varepsilon},g_{\varepsilon})$ be optimal entropic potentials
satisfying Eq. 9, and let $T_{\varepsilon}$ be the entropic map. Then
$T_{\varepsilon}=(Id-\nabla f_{\varepsilon})$.

###### Proof.

Eq. 9 implies

$f_{\varepsilon}(x)=-\varepsilon\log\int
e^{(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}Q(y)\,.$

Taking the gradient of this expression yields

$\displaystyle\nabla f_{\varepsilon}(x)$
$\displaystyle=-\varepsilon\frac{\int(-(x-y)/\varepsilon)e^{(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}Q(y)}{\int
e^{(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}Q(y)}$
$\displaystyle=x-\frac{\int
ye^{(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}Q(y)}{\int
e^{(g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}Q(y)}=x-T_{\varepsilon}(x)\,.$

∎

We write $P_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ and
$Q_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{Y_{i}}$ for the empirical
distributions corresponding to the samples from $P$ and $Q$, respectively. Our
proposed estimator is $T_{\varepsilon,(n,n)}$, the entropic map between
$P_{n}$ and $Q_{n}$, which can be written explicitly as

$T_{\varepsilon,(n,n)}(x)=\frac{\frac{1}{n}\sum_{i=1}^{n}Y_{i}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}{\frac{1}{n}\sum_{i=1}^{n}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}\,,$
(14)

where $g_{\varepsilon,(n,n)}$ is the entropic potential corresponding to
$Q_{n}$ in the optimal entropic plan between $P_{n}$ and $Q_{n}$, which can be
obtained as part of the output of Sinkhorn’s algorithm (see Peyré and Cuturi,
2019). In other words, once the optimal entropic plan is found, the map
$T_{\varepsilon,(n,n)}(x)$ can therefore be evaluated in linear time; this is
discussed more in Section 6. As in standard nonparametric estimation
(Tsybakov, 2009), the optimal choice of $\varepsilon$ will be dictated by the
smoothness of the target function.

To prove quantitative rates of convergence for $T_{\varepsilon,(n,n)}$, we
make the following regularity assumptions on $P$ and $Q$:

(A1) $P,Q\in\mathcal{P}_{\text{ac}}(\Omega)$ for a compact set $\Omega$, with
densities satisfying $p(x),q(x)\leq M$ and $q(x)\geq m>0$ for all
$x\in\Omega$,

(A2) $\varphi_{0}\in\mathcal{C}^{2}(\Omega)$ and
$\varphi_{0}^{*}\in\mathcal{C}^{\alpha+1}(\Omega)$ for $\alpha>1$,

(A3) $T_{0}=\nabla\varphi_{0}$, with $\mu
I\preceq\nabla^{2}\varphi_{0}(x)\preceq LI$ for $\mu,L>0$ for all
$x\in\Omega$,

In what follows, all constants may depend on the dimension, the set $\Omega$,
$M$, $m$, $\mu$, $L$, and $\|\varphi_{0}^{*}\|_{\mathcal{C}^{\alpha+1}}$.

######  Remark 2.

Assumptions (A1) to (A3) are similar to, but in general neither stronger nor
weaker than, those employed by Hütter and Rigollet (2021). On the one hand,
Hütter and Rigollet (2021) require that the set $\Omega$ be a bounded
Lipschitz domain, whereas we assume only compactness. On the other hand,
Hütter and Rigollet (2021) do not explicitly impose a smoothness requirement
on $\varphi_{0}^{*}$, though their assumptions do imply that it lies in the
class $\mathcal{C}^{2}(\Omega)$.

Our main result is the following.

######  Theorem 3.

Under assumptions (A1) to (A3), the entropic map
$\hat{T}=T_{\varepsilon,(n,n)}$ from $P_{n}$ to $Q_{n}$ with regularization
parameter $\varepsilon\asymp n^{-\frac{1}{d^{\prime}+\bar{\alpha}+1}}$
satisfies

$\mathbb{E}\|\hat{T}-T_{0}\|_{L^{2}(P)}^{2}\lesssim(1+I_{0}(P,Q))n^{-\frac{(\bar{\alpha}+1)}{2(d^{\prime}+\bar{\alpha}+1)}}\log
n\,,$

where $d^{\prime}=2\lceil d/2\rceil$ and $\bar{\alpha}=\alpha\wedge 3$.

When $d\to\infty$ and $\alpha\to 1$, we formally obtain the rate
$n^{-(1+o(1))/d}$. By contrast, Hütter and Rigollet (2021) show that, up to
logarithmic factors, the rate $n^{-2/d}$ is minimax optimal in this setting.
Theorem 3 therefore falls short of the minimax rate by a factor of
approximately 2 in the exponent; however, our numerical experiments in Section
6 show that $\hat{T}$ is competitive with minimax-optimal estimators in
practice.

To analyze $T_{\varepsilon,(n,n)}$, we adopt a two-step approach. We first
consider the one-sample setting and show that the entropic map
$T_{\varepsilon,n}$ between $P$ and $Q_{n}$ is close to $T_{0}$ in
expectation. We prove the following.

######  Theorem 4.

Under assumptions (A1) to (A3) there exists a constant $\varepsilon_{0}>0$
such that for $\varepsilon\leq\varepsilon_{0}$, the entropic map
$T_{\varepsilon,n}$ between $P$ and $Q_{n}$ satisfies

$\mathbb{E}\|T_{\varepsilon,n}-T_{0}\|_{L^{2}(P)}^{2}\lesssim\varepsilon^{1-d^{\prime}/2}\log(n)n^{-1/2}+\varepsilon^{(\bar{\alpha}+1)/2}+\varepsilon^{2}I_{0}(P,Q)\,,$
(15)

with $d^{\prime}=2\lceil d/2\rceil$ and $\bar{\alpha}=\alpha\wedge 3$.
Choosing $\varepsilon\asymp n^{-\frac{1}{d^{\prime}+\bar{\alpha}-1}}$, we get
the one-sample estimation rate

$\mathbb{E}\|T_{\varepsilon,n}-T_{0}\|^{2}_{L^{2}(P)}\lesssim(1+I_{0}(P,Q))n^{-\frac{\bar{\alpha}+1}{2(d^{\prime}+\bar{\alpha}-1)}}\,.$
(16)

######  Remark 3.

It can happen that $I_{0}(P,Q)$ is infinite, so the bounds of Theorems 3 and 4
are sometimes vacuous. However, Chizat et al. (2020) prove that
$I_{0}(P,Q)\leq C$ for a positive constant $C$ when $P$ and $Q$ satisfy (A1)
to (A3) for some $\alpha\geq 2$. Therefore, in the range $\alpha\in[2,3]$, we
obtain the rates in the theorems above without additional restrictions.

The proof of Theorem 4 is technical, and our approach is closely inspired by
Pal (2019) and empirical process theory arguments developed by Genevay et al.
(2019) and Mena and Niles-Weed (2019). We give a summary of our argument here,
and carry out the details in the following section.

Following Pal (2019), we define the _divergence_
$D[y|x^{*}]:=-x^{\top}y+\varphi_{0}(x)+\varphi_{0}^{*}(y)$, where
$\varphi_{0}$ solves Eq. 4. Though this quantity is a function of $x$ and $y$,
it is notationally convenient to write it in a way that highlights its
dependence on $x^{*}:=T_{0}(x)$. Indeed, we rely throughout on the following
fact

######  Lemma 1.

Under assumptions (A2) and (A3), for any
$x\in\operatorname{\mathrm{supp}}(P)$, we have

$D[y|x^{*}]=\frac{1}{2}(y-x^{*})^{\top}\nabla^{2}\varphi_{0}^{*}(x^{*})(y-x^{*})+o(\|y-x^{*}\|^{2})\quad\text{as
$y\to x^{*}$}\,,$ (17)

as well as the non-asymptotic bound

$\frac{1}{2L}\|y-x^{*}\|^{2}\leq
D[y|x^{*}]\leq\frac{1}{2\mu}\|y-x^{*}\|^{2}\,.$ (18)

###### Proof.

This follows from Taylor’s theorem and the fact that
$\nabla\varphi_{0}^{*}(x^{*})=T_{0}^{-1}(x^{*})=x$. ∎

We then define a conditional probability density in terms of this divergence:

$\displaystyle\quad
q_{\varepsilon}^{x}(y)=\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}e^{\frac{-1}{\varepsilon}D[y|x^{*}]}\,,\quad
Z_{\varepsilon}(x):=\frac{1}{\Lambda_{\varepsilon}}\int\exp\left(-\frac{1}{\varepsilon}D[y|x^{*}]\right)dy,$
(19)

for $\Lambda_{\varepsilon}=(2\pi\varepsilon)^{d/2}$. By virtue of Eq. 17, if
$\varphi_{0}^{*}$ is sufficiently smooth, then $q_{\varepsilon}^{x}$ will be
approximately Gaussian with mean $x^{*}$ and covariance
$\varepsilon\nabla^{2}\varphi_{0}^{*}(x^{*})^{-1}=\varepsilon\nabla^{2}\varphi_{0}(x)$.
We quantify this approximation via Laplace’s method, details about which
appear in Appendix B. Using variational arguments, reminiscent of those
employed by Bobkov and Götze (1999) in the study of transportation
inequalities, we then compare the measure $\pi_{\varepsilon,n}$ to the measure
$q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)$ and compute accurate
estimates of $T_{\varepsilon,n}$ via Laplace’s method.

A similar but much simpler argument establishes the following bound in the
two-sample case.

######  Theorem 5.

Let $T_{\varepsilon,(n,n)}$ be the entropic map from $P_{n}$ to $Q_{n}$, and
let $T_{\varepsilon,n}$ be as in Theorem 4. Under assumptions (A1) to (A3),
for $\varepsilon\leq 1$, $T_{\varepsilon,(n,n)}$ satisfies

$\mathbb{E}\|T_{\varepsilon,(n,n)}-T_{\varepsilon,n}\|_{L^{2}(P)}^{2}\lesssim\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,,$

where $d^{\prime}=2\lceil d/2\rceil$.

Combining Theorems 4 and 5 yields our main result.

###### Proof of Theorem 3.

We have

$\displaystyle\mathbb{E}\|T_{\varepsilon,(n,n)}-T_{0}\|^{2}_{L^{2}(P)}$
$\displaystyle\lesssim\mathbb{E}\|T_{\varepsilon,(n,n)}-T_{\varepsilon,n}\|^{2}_{L^{2}(P)}+\mathbb{E}\|T_{\varepsilon,n}-T_{0}\|^{2}_{L^{2}(P)}$
$\displaystyle\lesssim\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}+\varepsilon^{(\bar{\alpha}+1)/2}+\varepsilon^{2}I_{0}(P,Q)\,.$

Choosing $\varepsilon\asymp n^{-\frac{1}{d^{\prime}+\bar{\alpha}+1}}$ yields
the bound. ∎

##  4 One-sample estimates

In this section, we prove Theorem 4, which relates $T_{0}$ to the entropic map
between $P$ and $Q_{n}$:

$T_{\varepsilon,n}(x)=\frac{\int
ye^{\frac{1}{\varepsilon}(g_{\varepsilon,n}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q_{n}(y)}{\int
e^{\frac{1}{\varepsilon}(g_{\varepsilon,n}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q_{n}(y)}=\int
y\,\mathrm{d}\pi_{\varepsilon,n}^{x}(y)\,,$

where $\pi_{\varepsilon,n}$ is the optimal entropic plan for $P$ and $Q_{n}$.
We stress that since $T_{\varepsilon,n}$ is based on the entropic map from $P$
to $Q_{n}$, the second equality holds for $P$-almost every $x$.

Our main tool is the following inequality, which allows us to compare
$\pi_{\varepsilon,n}$ to the measure constructed from the conditional
densities $q_{\varepsilon}^{x}$. The proof relies crucially on Proposition 1
and on the second order-expansion provided in Theorem 2.

######  Proposition 3.

Assume (A1) to (A3), and let $a\in[L\varepsilon,1]$ for $\varepsilon\leq 1$.
Then

$\mathbb{E}\Big{\\{}\sup_{h:\Omega\to\mathbb{R}^{d}}\iint\left(h(x)^{\top}(y-T_{0}(x))-a\|h(x)\|^{2}\right)\,\mathrm{d}\pi_{\varepsilon,n}(x,y)\\\
-\iint(e^{h(x)^{\top}(y-T_{0}(x))-a\|h(x)\|^{2}}-1)q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)\Big{\\}}\\\
\lesssim\varepsilon
I_{0}+\varepsilon^{(\bar{\alpha}-1)/2}+\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,,$
(20)

where the supremum is taken over all bounded $h$.

###### Proof.

Given $h:\Omega\to\mathbb{R}^{d}$, write

$j_{h}(x,y)=h(x)^{\top}(y-T_{0}(x))-a\|h(x)\|^{2}\,.$

Choosing
$\eta(x,y)=\varepsilon(j_{h}(x,y)+\log(q_{\varepsilon}^{x}(y)/q(y)))+\|x-y\|^{2}/2$
and applying Proposition 1 with the measures $P$ and $Q_{n}$, we obtain

$\displaystyle\sup_{h:\Omega\to\mathbb{R}^{d}}\int
j_{h}\,\mathrm{d}\pi_{\varepsilon,n}$
$\displaystyle+\int\log\frac{q_{\varepsilon}^{x}(y)e^{\frac{1}{2\varepsilon}\|x-y\|^{2}}}{q(y)}\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle-\iint
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\,\mathrm{d}Q_{n}(y)\,\mathrm{d}P(x)+1\leq\varepsilon^{-1}S_{\varepsilon}(P,Q_{n})\,.$

We first expand
$\iint\log\frac{q_{\varepsilon}^{x}(y)e^{\frac{1}{2\varepsilon}\|x-y\|^{2}}}{q(y)}\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$,
where we use the fact that $\pi_{\varepsilon,n}$ has marginals $P$ and
$Q_{n}$:

$\displaystyle\iint\log\frac{q_{\varepsilon}^{x}(y)e^{\frac{1}{2\varepsilon}\|x-y\|^{2}}}{q(y)}\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle\qquad=\frac{1}{\varepsilon}\iint\bigg{[}f_{0}(x)+g_{0}(y)+\varepsilon\log\left(\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\right)-\varepsilon\log(q(y))\bigg{]}\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle\qquad=\frac{1}{\varepsilon}\Big{(}\int
f_{0}(x)\,\mathrm{d}P(x)+\int
g_{0}(y)\,\mathrm{d}Q_{n}(y)\Big{)}-\log(\Lambda_{\varepsilon})$
$\displaystyle\qquad\phantom{=}-\int\log(Z_{\varepsilon}(x))\,\mathrm{d}P(x)-\int\log(q(y))\,\mathrm{d}Q_{n}(y)\,,$

where $(f_{0},g_{0})$ solve (6). Replacing $Q_{n}$ by $Q$ yields

$\displaystyle\iint\log\frac{q_{\varepsilon}^{x}(y)e^{\frac{1}{2\varepsilon}\|x-y\|^{2}}}{q(y)}\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle=\frac{1}{2\varepsilon}W_{2}^{2}(P,Q)-\log(\Lambda_{\varepsilon})-\int\log(Z_{\varepsilon}(x))\,\mathrm{d}P(x)$
$\displaystyle\phantom{=}-\operatorname{Ent}(Q)+\int(g_{0}/\varepsilon-\log(q))(\mathrm{d}Q_{n}-\mathrm{d}Q).$

A change of variables (see Pal, 2019, Lemma 3(iv)) implies

$\frac{\operatorname{Ent}(Q)-\operatorname{Ent}(P)}{2}=\int\log
J(\nabla^{2}\varphi_{0}^{*}(x^{*}))\,\mathrm{d}P(x)\,,$

where we recall that $x^{*}=T_{0}(x)$. Substituting this identity into the
preceding expression yields

$\displaystyle\iint\log\frac{q_{\varepsilon}^{x}(y)e^{\frac{1}{2\varepsilon}\|x-y\|^{2}}}{q(y)}\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle=\frac{1}{2\varepsilon}W_{2}^{2}(P,Q)-\log(\Lambda_{\varepsilon})-\frac{1}{2}(\operatorname{Ent}(Q)+\operatorname{Ent}(P))$
$\displaystyle\phantom{=}+\int(g_{0}/\varepsilon-\log(q))(\mathrm{d}Q_{n}-\mathrm{d}Q)$
$\displaystyle\phantom{=}-\int\log(Z_{\varepsilon}(x)J(\nabla^{2}\varphi_{0}^{*}(x^{*})))\,\mathrm{d}P(x)\,.$

We therefore obtain

$\sup_{h:\Omega\to\mathbb{R}^{d}}\int
j_{h}\,\mathrm{d}\pi_{\varepsilon,n}-\iint
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\,\mathrm{d}Q_{n}(y)\,\mathrm{d}P(x)+1\\\
\leq\varepsilon^{-1}\Big{(}S_{\varepsilon}(P,Q_{n})-\frac{1}{2}W_{2}^{2}(P,Q)+\varepsilon\log(\Lambda_{\varepsilon})+\frac{\varepsilon}{2}(\operatorname{Ent}(Q)+\operatorname{Ent}(P))\Big{)}+\Delta_{1}\,,$

where
$\Delta_{1}:=\int(g_{0}/\varepsilon-\log(q))(\mathrm{d}Q-\mathrm{d}Q_{n})+\int\log(Z_{\varepsilon}(x)J(\nabla^{2}\varphi_{0}^{*}(x^{*})))\,\mathrm{d}P(x)$.
Applying Theorem 2, we may further bound

$\sup_{h:\Omega\to\mathbb{R}^{d}}\int
j_{h}\,\mathrm{d}\pi_{\varepsilon,n}-\iint
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\,\mathrm{d}Q_{n}(y)\,\mathrm{d}P(x)+1\leq\frac{\varepsilon}{8}I_{0}+\Delta_{1}+\Delta_{2}\,,$

where
$\Delta_{2}:=\varepsilon^{-1}(S_{\varepsilon}(P,Q_{n})-S_{\varepsilon}(P,Q))$.
Now we turn our attention to the second term on the left side. Since

$\displaystyle\iint
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\,\mathrm{d}Q(y)\,\mathrm{d}P(x)$
$\displaystyle=\iint_{\operatorname{\mathrm{supp}}(Q)}e^{j_{h}(x,y)}q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)$
$\displaystyle\leq\iint
e^{j_{h}(x,y)}q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)\,,$

we have

$\sup_{h:\Omega\to\mathbb{R}^{d}}\int
j_{h}\,\mathrm{d}\pi_{\varepsilon,n}-\iint(e^{j_{h}(x,y)}-1)q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)\leq\frac{\varepsilon}{8}I_{0}+\Delta_{1}+\Delta_{2}+\Delta_{3}\,,$

where

$\Delta_{3}:=\sup_{h:\Omega\to\mathbb{R}^{d}}\iint
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\,\mathrm{d}P(x)(\mathrm{d}Q_{n}-\mathrm{d}Q)(y)$

and where we have used the fact that $q_{\varepsilon}^{x}(y)$ is a probability
density.

It therefore remains only to show that

$\mathbb{E}[\Delta_{1}+\Delta_{2}+\Delta_{3}]\lesssim\varepsilon^{(\bar{\alpha}-1)/2}+\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,.$

First, a Laplace expansion (Corollary 2) implies

$\mathbb{E}\Delta_{1}=\int\log(Z_{\varepsilon}(x)J(\nabla^{2}\varphi_{0}^{*}(x^{*})))\,\mathrm{d}P(x)\lesssim\varepsilon^{(\bar{\alpha}-1)/2}$

Second, known results on the finite-sample convergence of the Sinkhorn
divergence (Corollary 3) yield

$\mathbb{E}\Delta_{2}\lesssim(\varepsilon^{-1}+\varepsilon^{-d^{\prime}/2})\log(n)n^{-1/2}\,,$

It therefore remains to bound $\Delta_{3}$, which an empirical process theory
argument (Proposition 6) shows

$\mathbb{E}\Delta_{3}\lesssim\varepsilon^{-d/2}n^{-1/2}$

as long as $a\in[L\varepsilon,1]$.

We obtain that

$\mathbb{E}[\Delta_{1}+\Delta_{2}+\Delta_{3}]\lesssim\varepsilon^{(\bar{\alpha}-1)/2}+(\varepsilon^{-1}+\varepsilon^{-d^{\prime}/2})\log(n)n^{-1/2}+\varepsilon^{-d/2}n^{-1/2}\,,$

and since $d^{\prime}\geq d$ and $\varepsilon\leq 1$, we obtain the bound

$\mathbb{E}[\Delta_{1}+\Delta_{2}+\Delta_{3}]\lesssim\varepsilon^{(\bar{\alpha}-1)/2}+\varepsilon^{-d^{\prime}/2}n^{-1/2}\log(n)\,,$

as desired. ∎

To exploit Proposition 3, we show that we can choose a function $h$ for which
the left side of the above expression scales like
$\|T_{\varepsilon,n}-T_{0}\|_{L^{2}(P)}^{2}$.

We first establish three lemmas, whose proofs are deferred.

######  Lemma 2.

Fix $x\in\operatorname{\mathrm{supp}}(P)$, and write $\bar{y}^{x}=\int
yq_{\varepsilon}^{x}(y)\,\mathrm{d}y$. There exists a positive constant $C$,
independent of $x$, such that for all $\varepsilon\in(0,1)$ and $\|v\|_{2}\leq
1$,

$\int
e^{(v^{\top}(y-\bar{y}^{x}))^{2}/(C\varepsilon)}q_{\varepsilon}^{x}(y)\,\mathrm{d}y\leq
2\,.$

In probabilistic language, Lemma 2 implies that if $Y^{x}$ is a random
variable with density $q_{\varepsilon}^{x}$, then
$\varepsilon^{-1/2}(Y^{x}-\mathbb{E}Y^{x})$ is subgaussian (Vershynin, 2018).
By applying standard moment bounds for subgaussian random variables, we then
arrive at the following result.

######  Lemma 3.

There exists a positive constant $C$ such that if $a\geq C\varepsilon$, then
for any $h:\mathbb{R}^{d}\to\mathbb{R}^{d}$ we have

$\iint
e^{h(x)^{\top}(y-T_{0}(x))-a\|h(x)\|^{2}}q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)\leq\int
e^{\frac{1}{4\varepsilon}\|\bar{y}^{x}-T_{0}(x)\|^{2}}\,\mathrm{d}P(x)\,.$

Finally, we show by an application of Laplace’s method that $\bar{y}^{x}$ is
close to $T_{0}(x)$.

######  Lemma 4.

Assume (A1) to (A3). For all $x\in\operatorname{\mathrm{supp}}(P)$,

$\|\bar{y}^{x}-T_{0}(x)\|^{2}\lesssim\varepsilon^{\alpha\wedge 2}\,.$

With these lemmas in hand, we can complete the proof.

###### Proof of Theorem 4.

We may assume $\varepsilon_{0}\leq 1$. Since $e^{t}-1\leq 2t$ for $t\in[0,1]$,
Lemma 4 implies that as long as $\varepsilon_{0}$ is sufficiently small, for
$\varepsilon\leq\varepsilon_{0}$,

$e^{\frac{1}{4\varepsilon}\|\bar{y}^{x}-T_{0}(x)\|^{2}}-1\lesssim\varepsilon^{(\alpha-1)\wedge
1}\leq\varepsilon^{(\bar{\alpha}-1)/2}\,,$

where the last inequality holds for $\alpha\geq 1$ and $\varepsilon\leq 1$.
Combining this fact with Lemma 3, we obtain that for any
$h:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and $a\geq C\varepsilon$,

$\iint(e^{h(x)^{\top}(y-T_{0}(x))-a\|h(x)\|^{2}}-1)q_{\varepsilon}^{x}(y)\,\mathrm{d}y\,\mathrm{d}P(x)\lesssim\varepsilon^{(\bar{\alpha}-1)/2}\,.$

For a sufficiently small constant $\varepsilon_{0}$, the interval
$[C\varepsilon,1]$ is non-empty for $\varepsilon\leq\varepsilon_{0}$, so
combining this fact with Proposition 3 yields that for $a\in[C\varepsilon,1]$
and $\varepsilon\leq\varepsilon_{0}$,

$\mathbb{E}\sup_{h:\Omega\to\mathbb{R}^{d}}\iint\left(h(x)^{\top}(y-T_{0}(x))-a\|h(x)\|^{2}\right)\,\mathrm{d}\pi_{\varepsilon,n}(x,y)\lesssim\varepsilon
I_{0}\\\
+\varepsilon^{(\bar{\alpha}-1)/2}+\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,.$
(21)

If we pick $h(x)=\frac{1}{2a}(T_{\varepsilon,n}(x)-T_{0}(x))$, the integral on
the left side equals

$\frac{1}{2a}\mathbb{E}\iint\left((T_{\varepsilon,n}(x)-T_{0}(x))^{\top}(y-T_{0}(x))-\frac{1}{2}\|T_{\varepsilon,n}(x)-T_{0}(x)\|^{2}\right)\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
(22)

By definition, $T_{\varepsilon,n}(x)=\int
y\,\mathrm{d}\pi_{\varepsilon,n}^{x}(y)$, so disintegrating
$\pi_{\varepsilon,n}(x,y)$ and recalling that the first marginal of
$\pi_{\varepsilon,n}$ is $P$ yields

$\iint\left(T_{\varepsilon,n}(x)-T_{0}(x))^{\top}(y-T_{0}(x))-\frac{1}{2}\|T_{\varepsilon,n}(x)-T_{0}(x)\|^{2}\right)\,\mathrm{d}\pi_{\varepsilon,n}(x,y)\\\
=\int\frac{1}{2}\|T_{\varepsilon,n}(x)-T_{0}(x)\|^{2}\,\mathrm{d}P(x)=\frac{1}{2}\|T_{\varepsilon,n}-T_{0}\|_{L^{2}(P)}^{2}\,.$

Combining this with Eq. 21 and Eq. 22 and picking $a=C\varepsilon$ yields

$\mathbb{E}\|T_{\varepsilon,n}-T_{0}\|_{L^{2}(P)}^{2}\lesssim\varepsilon^{2}I_{0}+\varepsilon^{(\bar{\alpha}+1)/2}+\varepsilon^{1-d^{\prime}/2}\log(n)n^{-1/2}\,,$

as desired. ∎

As a corollary to Theorem 4, we have the following population-level estimate
between $T_{\varepsilon}$ and $T_{0}$, which is potentially of independent
interest.

######  Corollary 1.

Assume (A1) to (A3), then

$\displaystyle\|T_{\varepsilon}-T_{0}\|^{2}_{L^{2}(P)}=\|\nabla
f_{\varepsilon}-\nabla
f_{0}\|^{2}_{L^{2}(P)}\lesssim\varepsilon^{2}I_{0}(P,Q)+\varepsilon^{(\bar{\alpha}+1)/2},$
(23)

where $\bar{\alpha}=3\wedge\alpha$.

##  5 Two-sample estimates

We now turn our attention to the two-sample case. Let
$\pi_{\varepsilon,(n,n)}$ be the optimal entropic plan between $P_{n}$ and
$Q_{n}$ and $(f_{\varepsilon,(n,n)},g_{\varepsilon,(n,n)})$ the corresponding
entropic potentials. We aim to show that

$\mathbb{E}\|T_{\varepsilon,(n,n)}-T_{\varepsilon,n}\|_{L^{2}(P)}^{2}\lesssim(\varepsilon^{-1}+\varepsilon^{-d^{\prime}/2})\log(n)n^{-1/2}\,,$

where $d^{\prime}=2\lceil d/2\rceil$.

As in Section 4, we proceed via duality arguments, but our task is
considerably simplified by the fact that the measure $Q_{n}$ remains fixed in
passing from $T_{\varepsilon,(n,n)}$ to $T_{\varepsilon,n}$. Let us write

$\gamma(x,y)=e^{\frac{1}{\varepsilon}(f_{\varepsilon,(n,n)}(x)+g_{\varepsilon,(n,n)}(y)-\frac{1}{2}\|x-y\|^{2})}=\frac{e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(y)-\frac{1}{2}\|x-y\|^{2})}}{\frac{1}{n}\sum_{i=1}^{n}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}$

for the $P_{n}\otimes Q_{n}$ density of $\pi_{\varepsilon,(n,n)}$, where the
second equality holds $P_{n}\otimes Q_{n}$ almost everywhere and furnishes an
extension of $\gamma$ to all $x\in\mathbb{R}^{d}$.

We employ the following analogue of Proposition 3, which does not require the
full force of assumptions (A1) to (A3).

######  Proposition 4.

The support of $P$ and $Q$ lies in $\Omega$, then

$\displaystyle\mathbb{E}\Big{\\{}\sup_{\chi:\Omega\times\Omega\to\mathbb{R}}\iint\chi(x,y)\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle-\iint(e^{\chi(x,y)}-1)\gamma(x,y)\,\mathrm{d}P(x)\,\mathrm{d}Q_{n}(y)\Big{\\}}$
$\displaystyle\lesssim(\varepsilon^{-1}+\varepsilon^{-d^{\prime}/2})\log(n)n^{-1/2}.$

The proof of Theorem 5 is now straightforward.

###### Proof.

As in the proof of Theorem 4, consider

$\chi(x,y)=h(x)^{\top}(y-T_{\varepsilon,(n,n)}(x))-a\|h(x)\|^{2}$

for $h$ and $a$ to be specified. By definition of $T_{\varepsilon,(n,n)}$, we
have

$\int
h(x)^{\top}(y-T_{\varepsilon,(n,n)}(x))\gamma(x,y)\,\mathrm{d}Q_{n}(y)=h(x)^{\top}\left(\int
y\gamma(x,y)\,\mathrm{d}Q_{n}(y)-T_{\varepsilon,(n,n)(x)}\right)\\\
=h(x)^{\top}\left(\frac{\frac{1}{n}\sum_{i=1}^{n}Y_{i}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-c(x,Y_{i}))}}{\frac{1}{n}\sum_{i=1}^{n}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-c(x,Y_{i}))}}-T_{\varepsilon,(n,n)}(x)\right)=0$

for all $x\in\mathbb{R}^{d}$. Moreover, since $\Omega$ is compact, there
exists a constant $C$ such that

$|h(x)^{\top}(y-T_{\varepsilon,(n,n)}(x))|\leq C\|h(x)\|\quad\forall
y\in\Omega\,.$

Hoeffding’s inequality therefore implies that if $a\geq C^{2}/2$, then this
choice of $\chi$ satisfies

$\iint(e^{\chi(x,y)}-1)\gamma(x,y)\,\mathrm{d}Q_{n}(y)\,\mathrm{d}P(x)\leq
0\,.$

Choosing $h(x)=\frac{1}{2a}(T_{\varepsilon,n}(x)-T_{\varepsilon,(n,n)}(x))$,
we conclude as in the proof of Theorem 4 that for $\varepsilon\leq 1$,

$\frac{1}{4a}\|T_{\varepsilon,n}-T_{\varepsilon,(n,n)}\|_{L^{2}(P)}^{2}\lesssim(\varepsilon^{-1}+\varepsilon^{-d^{\prime}/2})\log(n)n^{-1/2}\lesssim\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,,$

and picking $a$ to be a sufficiently large constant yields the claim. ∎

##  6 Computational aspects

Our reason for studying the entropic map as an optimal transport map estimator
arises from its strong computational benefits, which are a consequence of the
efficiency of Sinkhorn’s algorithm for entropic optimal transport (see Peyré
and Cuturi, 2019). In this section, we compare the computational complexity of
the entropic map to the estimators of Hütter and Rigollet (2021), Deb et al.
(2021), and Manole et al. (2021) in the two-sample setting. Finally, we
perform several experiments that demonstrate the computational advantages of
our procedure. Throughout this section, we use $\tilde{O}$ to hide poly-
logarithmic factors in the sample size $n$.

###  6.1 Estimator complexities from prior work

We first describe the wavelet-based estimator proposed by Hütter and Rigollet
(2021). Recall that this estimator is minimax optimal for all $\alpha>1$. The
implementation of this estimator requires various discretization/approximation
schemes. The authors of that work use a numerical implementation of the
Daubechies wavelets to approximate the optimal Brenier potential, and then
compute its convex conjugate by means of a discrete Legendre transform on a
discrete grid. The gradient of the resulting potential is then obtained using
finite differences, and this is extended to data outside the grid by linear
interpolation. Though computing this estimator takes time that scales only
linearly in the sample size $n$, the main bottleneck of this approach from a
computational standpoint is the computation of the Legendre transform on the
grid, which requires at least $cN^{d}$ operations, where $N$ denotes the
resolution of the grid. Since this resolution needs to be chosen fine enough
to be negligible, the exponential dependence in $d$ makes this approach
prohibitively expensive in most applications.

Another estimator recently analyzed in the literature is the “1-Nearest
Neighbor” estimator, which we denote by $\hat{T}^{\text{1NN}}_{(n,n)}$ (Manole
et al., 2021), which achieves the minimax rate when $T_{0}$ is Lipschitz
(i.e., $\alpha=1$). The estimator takes the form

$\displaystyle\hat{T}^{\text{1NN}}_{(n,n)}(x)=\sum_{i,j=1}^{n}(n\hat{\pi}_{ij})\bm{1}_{V_{i}}(x)Y_{j},$
(24)

where $\bm{1}$ is the indicator function for a set, and $(V_{i})_{i=1}^{n}$
are the Voronoi regions generated by $(X_{i})_{i=1}^{n}$, i.e.

$V_{i}=\\{x\in\Omega\ :\ \|x-X_{i}\|\leq\|x-X_{j}\|,\ \forall j\neq i\\},$

and $\hat{\pi}$ is the optimal coupling that solves Eq. 5 when the measures
are the empirical measures $P_{n}$ and $Q_{n}$. Solving for $\hat{\pi}$ can be
done through the Hungarian algorithm, and has time complexity ${O}(n^{3})$.
However, unlike the wavelet estimator described above, computing this
estimator does not require constructing a grid whose size scales exponentially
with dimension.

For the $\alpha>1$ case, both Manole et al. (2021) and Deb et al. (2021)
propose estimators based on barycentric projections of density estimators. For
these approaches, the idea is to construct nonparametric density estimates of
the measures $P$ and $Q$, resample points from these densities, and finally
perform the appropriate matching using the Hungarian algorithm once again.
Though tractable in low dimensions, this approach is limited by the difficulty
of sampling from nonparametric density estimates, which typically requires
time scaling exponentially in $d$.

In short, prior estimators proposed in the literature either have runtime
scaling exponentially in $d$ (in the case of the wavelet estimator the
estimators based on nonparametric density estimation) or cubicly in $n$ (in
the case of the 1NN estimator). By contrast, in the following section, we show
that our estimator can be computed in nearly $O(n^{2})$ time.

###  6.2 Computational complexity of the Entropic Map

We now turn to the computational analysis of our estimator, which has the
closed-form representation

$\hat{T}_{\varepsilon,(n,n)}(x)=\frac{\sum_{i=1}^{n}Y_{i}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}{\sum_{i=1}^{n}e^{\frac{1}{\varepsilon}(g_{\varepsilon,(n,n)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}\,.$
(25)

The computational burden of our estimator falls on computing the optimal
entropic potential evaluated at the data $g_{\varepsilon,(n,n)}(Y_{i})$.
Indeed, once we have this potential, it is clear that the remainder of
Equation 25 can be computed in $O(n)$ time.

The leading approach to compute optimal entropic potentials in practice is
_Sinkhorn’s algorithm_ (Peyré and Cuturi, 2019), an alternating maximization
algorithm that computes approximations of the entropic potentials by
iteratively updating $f$ and $g$ so that they satisfy one of the two dual
optimality conditions given in Eq. 9. Explicitly, defining $f^{(0)}=0$,
Sinkhorn’s algorithm performs the updates

$\displaystyle g^{(k)}(y)$
$\displaystyle=-\varepsilon\log\frac{1}{n}\sum_{i=1}^{n}e^{\frac{1}{\varepsilon}(f^{(k)}(X_{i})-\frac{1}{2}\|X_{i}-y\|^{2})}$
$\displaystyle f^{(k+1)}(y)$
$\displaystyle=-\varepsilon\log\frac{1}{n}\sum_{j=1}^{n}e^{\frac{1}{\varepsilon}(g^{(k)}(Y_{j})-\frac{1}{2}\|x-Y_{j}\|^{2})}\,.$

until termination. Since it is only necessary to compute $f^{(k)}$ and
$g^{(k)}$ on the support of $P_{n}$ and $Q_{n}$, respectively, each iteration
can be implemented in $O(n^{2})$ time.

Note that this update rule guarantees that

$\int
e^{\frac{1}{\varepsilon}(f^{(k)}(x)+g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}P_{n}(x)=1$

for all $y$ at each iteration. By contrast, the second optimality condition in
Eq. 9 is _not_ satisfied at each iteration, though Sinkhorn (1967) showed that

$\int
e^{\frac{1}{\varepsilon}(f^{(k)}(x)+g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q_{n}(y)\to
1$

as $k\to\infty$, and therefore that the iterates of Sinkhorn’s algorithm
converge to optimal entropic potentials.

To analyze the running time of our estimator, we will leverage recent analyses
of the convergence rate of Sinkhorn’s algorithm (Cuturi, 2013; Altschuler et
al., 2017; Dvurechensky et al., 2018) to explicitly quantify the error
incurred by terminating after a finite number of steps. For $k\geq 0$, we
consider the entropic map estimator obtained after $k$ iterates of Sinkhorn’s
algorithm:

$\displaystyle
T^{(k)}(x)=\frac{\sum_{i=1}^{n}Y_{i}e^{\frac{1}{\varepsilon}(g^{(k)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}{\sum_{i=1}^{n}e^{\frac{1}{\varepsilon}(g^{(k)}(Y_{i})-\frac{1}{2}\|x-Y_{i}\|^{2})}}\,.$
(26)

Despite the fact that $g^{(k)}$ is _not_ an entropic potential for the
original problem, the following theorem shows that $T^{(k)}$ is nevertheless
an acceptable estimator if $k$ is sufficiently large.

######  Theorem 6.

Suppose assumptions (A1) to (A3) hold, and we choose $\varepsilon$ as in
Theorem 3. Then for any $k\gtrsim n^{7/(d^{\prime}+\bar{\alpha}+1)}\log n$,

$\mathbb{E}\|T^{(k)}-T_{0}\|^{2}_{L^{2}(P)}\lesssim(1+I_{0}(P,Q))n^{-\frac{(\bar{\alpha}+1)}{2(d^{\prime}+\bar{\alpha}+1)}}\log
n\,,$

where $d^{\prime}=2\lceil d/2\rceil$ and $\bar{\alpha}=3\wedge\alpha$. In
particular, an estimator achieving the same rate as the estimator in Theorem 3
can be computed in
$\tilde{O}(n^{2+7/(d^{\prime}+\bar{\alpha}+1)})=n^{2+o_{d}(1)}$ time.

###### Proof.

We begin by decomposing the error and applying Theorem 4:

$\displaystyle\mathbb{E}\|T^{(k)}-T_{0}\|^{2}_{L^{2}(P)}$
$\displaystyle\lesssim\mathbb{E}\|T^{(k)}-T_{\varepsilon,n}\|^{2}_{L^{2}(P)}+\mathbb{E}\|T_{\varepsilon,n}-T_{0}\|^{2}_{L^{2}(P)}$
$\displaystyle\lesssim\mathbb{E}\|T^{(k)}-T_{\varepsilon,n}\|^{2}_{L^{2}(P)}+\varepsilon^{1-d^{\prime}/2}\log(n)n^{-1/2}+\varepsilon^{(\bar{\alpha}+1)/2}+\varepsilon^{2}I_{0}(P,Q)\,.$

We proceed almost exactly as in Theorem 5, and consider

$\chi(x,y)=h(x)^{\top}\left(y-{T}^{(k)}\right)-a\|h(x)\|^{2},$

for $h:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and $a$ to be specified. For
$x\in\mathbb{R}^{d}$, $y\in\operatorname{\mathrm{supp}}(Q_{n})$, define

$\tilde{\gamma}(x,y)=\frac{\exp\left(\frac{1}{\varepsilon}(g^{(k)}(y)-\tfrac{1}{2}\|x-y\|^{2})\right)}{\frac{1}{n}\sum_{i=1}^{n}\exp\left(\frac{1}{\varepsilon}(g^{(k)}(Y_{i})-\tfrac{1}{2}\|x-Y_{i}\|^{2})\right)}\,.$
(27)

By construction, $\int\tilde{\gamma}(x,y)\,\mathrm{d}Q_{n}(y)=1$ for all
$x\in\mathbb{R}^{d}$, and $T^{(k)}(x)=\int
y\tilde{\gamma}(x,y)\,\mathrm{d}Q_{n}(y)$. Therefore, for any
$h:\mathbb{R}^{d}\to\mathbb{R}^{d}$,

$\int
h(x)^{\top}\left(y-{T}^{(k)}(x)\right)\tilde{\gamma}(x,y)\,\mathrm{d}Q_{n}(y)=0$

for all $x\in\mathbb{R}^{d}$. Moreover, since $\Omega$ is compact, there
exists a constant $C$ such that

$\displaystyle|h(x)^{\top}(y-T^{(k)}(x))|\leq C\|h(x)\|\quad\forall
x,y\in\Omega\,.$

Hoeffding’s inequality therefore implies that for $a$ sufficiently large, this
choice of $\chi$ satisfies

$\displaystyle\iint(e^{\chi(x,y)}-1)\tilde{\gamma}(x,y)\,\mathrm{d}Q_{n}\,\mathrm{d}P(x)\leq
0\,.$

Now, define a probability measure $\tilde{P}$ with the same support as $P_{n}$
by setting

$\frac{\mathrm{d}\tilde{P}(x)}{\mathrm{d}P_{n}(x)}=\int
e^{\frac{1}{\varepsilon}(f^{(k)}(x)+g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q_{n}(y)\,,$
(28)

and let

$f^{(k+1)}(x)=-\varepsilon\log\frac{1}{n}\sum_{i=1}^{n}\exp\left(\varepsilon^{-1}(g^{(k)}(Y_{i})-\tfrac{1}{2}\|x-Y_{i}\|^{2})\right)\,.$
(29)

We claim that
$\tilde{\gamma}(x,y)=\exp(\frac{1}{\varepsilon}(f^{(k+1)}(x)+g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})$
is the $\tilde{P}\otimes Q_{n}$ density of the optimal entropic plan between
$\tilde{P}$ and $Q_{n}$. We have already observed that
$\int\tilde{\gamma}(x,y)\,\mathrm{d}Q_{n}(y)=1$ for all $x\in\mathbb{R}^{d}$
by construction, so it suffices to note that for all
$y\in\operatorname{\mathrm{supp}}(Q_{n})$,

$\displaystyle\int\tilde{\gamma}(x,y)\,\mathrm{d}\tilde{P}(x)$
$\displaystyle=\int\frac{e^{\varepsilon^{-1}(g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})}}{\int
e^{\varepsilon^{-1}(g^{(k)}(y^{\prime})-\frac{1}{2}\|x-y^{\prime}\|^{2})}\,\mathrm{d}Q_{n}(y^{\prime})}\,\mathrm{d}\tilde{P}(x)$
$\displaystyle=\int\frac{e^{\varepsilon^{-1}(f^{(k)}(x)+g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})}}{\int
e^{\varepsilon^{-1}(f^{(k)}(x)+g^{(k)}(y^{\prime})-\frac{1}{2}\|x-y^{\prime}\|^{2})}\,\mathrm{d}Q_{n}(y^{\prime})}\,\mathrm{d}\tilde{P}(x)$
$\displaystyle=\int
e^{\varepsilon^{-1}(f^{(k)}(x)+g^{(k)}(y)-c(x,y))}\,\mathrm{d}P_{n}(x)=1\,.$

Therefore $(f^{(k+1)},g^{(k)})$ satisfy Eq. 9, so $\tilde{\gamma}$ is indeed
the $\tilde{P}\otimes Q_{n}$ density of the optimal entropic plan between the
two measures.

Applying Proposition 5, we obtain for any $\varepsilon\leq 1$

$\mathbb{E}\sup_{h:\mathbb{R}^{d}\to\mathbb{R}^{d}}\iint
h(x)^{\top}\left(y-{T}^{(k)}\right)-a\|h(x)\|^{2}\,\mathrm{d}\pi_{\varepsilon,n}\lesssim\varepsilon^{-1}\delta+\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,,$
(30)

where $\delta:=\mathrm{d_{TV}}(\tilde{P},P_{n})$. Choosing
$h(x)=\frac{1}{2a}\left(T_{\varepsilon,n}(x)-{T}^{(k)}(x)\right)$, we conclude
as in Theorem 5, resulting in

$\displaystyle\mathbb{E}\|T^{(k)}-T_{\varepsilon,n}\|^{2}_{L^{2}(P)}\lesssim\varepsilon^{-1}\delta+\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}\,.$

All together, we have

$\displaystyle\mathbb{E}\|T^{(k)}-T_{0}\|^{2}_{L^{2}(P)}$
$\displaystyle\lesssim\varepsilon^{-1}\delta+\varepsilon^{-d^{\prime}/2}\log(n)n^{-1/2}+\varepsilon^{(\bar{\alpha}+1)/2}+\varepsilon^{2}I_{0}(P,Q)\,.$

The first term will be negligible if $\delta\lesssim\varepsilon^{3}$.

By definition, $\tilde{P}$ is the first marginal of the joint distribution
with $P_{n}\otimes Q_{n}$ density
$e^{\frac{1}{\varepsilon}(f^{(k)}(x)+g^{(k)}(y)-\frac{1}{2}\|x-y\|^{2})}$. By
Altschuler et al. (2017, Theorem 2), if $k$ satisfies

$k\lesssim\delta^{-2}\log(n\cdot\max_{i,j}e^{\frac{1}{2\varepsilon}\|x_{i}-y_{j}\|^{2}})\lesssim\delta^{-2}\varepsilon^{-1}\log
n\,,$

then $\mathrm{d_{TV}}(\tilde{P},P_{n})\leq\delta$. Choosing
$\delta=\varepsilon^{3}\asymp n^{-3/(d^{\prime}+\bar{\alpha}+1)}$ yields the
claim. ∎

###  6.3 Empirical performance

We test two implementations of Sinkhorn’s algorithm, one from the Python
Optimal Transport (POT) library (Flamary et al., 2021), and an implementation
that uses the KeOps library optimized for GPUs. Both implementations employ
log-domain stabilization to avoid numerical overflow issues arising from the
small choice of $\varepsilon$.

For simplicity, we employ the same experimental setup as Hütter and Rigollet
(2021). We generate i.i.d. samples from a source distribution $P$, which we
always take to be $[-1,1]^{d}$, and from a target distribution
$Q=(T_{0})_{\sharp}P$, where we define $T_{0}:\mathbb{R}^{d}\to\mathbb{R}^{d}$
to be an optimal transport map obtained by applying a monotone scalar function
coordinate-wise.

In Fig. 1, we visualize the output of our estimator in $d=2$. The figures
depict the effect of evaluating the estimator $\hat{T}_{\varepsilon}$ and the
true map $T_{0}$ on additional test points
$X_{1}^{\prime},\dots,X^{\prime}_{m}$ drawn i.i.d. from $P$.

(a) $T_{0}(x)=\exp(x)$ coordinate-wise

(b) $T_{0}(x)=3x^{2}\text{sign}(x)$ coordinate-wise

Figure 1: Visualization of $\hat{T}_{\varepsilon}$ and $T_{0}(x)$ in 2
dimensions.

####  6.3.1 Comparison to a tractable minimax estimator

Among the previously discussed estimators, the 1-Nearest Neighbor estimator
analyzed in Manole et al. (2021) is the most tractable, and the only one
remotely comparable to our method. As discussed in Section 6.1, this approach
uses the Hungarian algorithm which has a runtime of $O(n^{3})$. However, since
it is not parallelizable, we compare its performance to the non-parallel CPU
implementation of Sinkhorn’s algorithm from the POT library.

We perform a simple experiment comparing our approach to theirs: let
$P=[-1,1]^{d}$ and let $T_{0}(x)=\exp(x)$, acting coordinate-wise. We vary $d$
and $n$, and track runtime performance of both estimators, as well as the Mean
Squared Error (MSE) of the map estimate222We calculate MSE by performing Monte
Carlo integration over the space $[-1,1]^{d}$., averaged over 20 runs. For our
estimator, we choose $\varepsilon$ as suggested by Theorem 3.

(a) MSE comparison

(b) Runtime comparison

Figure 2: Dashed lines are our estimator, solid lines are
$\hat{T}^{\text{1NN}}$, and $T_{0}(x)=\exp(x)$

We observe that in $d=2$, the MSE of the two estimators are comparable, though
our error deteriorates for large $n$, which reflects our slightly sub-optimal
estimation rate. However, as $d$ increases to moderate dimensions, our
estimator consistently outperforms $\hat{T}^{\text{1NN}}$ in both MSE and
runtime with the choice of $\varepsilon$ in Theorem 3. For both estimators,
the CPU runtime begins to become significant (on the order of seconds) when
$n$ exceeds 1500.

####  6.3.2 Parallel estimation on massive data sets

Figure 2 makes clear that computation of both estimators slows for $n\gg
10^{3}$ when implemented on a CPU. However, Sinkhorn’s algorithm can be easily
parallelized. Unlike the 1-Nearest Neighbor estimator—and all other transport
map estimators of which we are aware—our proposal therefore runs extremely
efficiently on GPUs. We again average performance over 20 runs, and choose
$\varepsilon$ as in the previous example, with $T_{0}$ again as the
exponential map (coordinate-wise). We see that even when $n=10^{4}$ and
$d=10$, it takes roughly a third of a second to perform the optimization.

(a) MSE comparison

(b) Runtime comparison

Figure 3: Performance of a parallel implementation of our estimator on large
data sets.

##  7 Conclusion and future directions

We have presented the first finite-sample analysis of an entropic estimator
for optimal transport maps. The resulting estimator is easily parallelizable
and fast to compute, even on massive data sets. Though the theoretical rates
we obtain fall short of minimax optimal, we demonstrate that our estimator
empirically outperforms several other statistically optimal proposals in the
literature.

Based on the empirical success of our estimator, we conjecture that our
analysis is loose, and that it may be possible to show that a properly tuned
version of an entropic estimator achieves the minimax optimal rate, at least
when $\alpha$ is small. We also conjecture that our assumptions may be
significantly loosened, and that similar results hold without stringent
conditions on the densities or their support.

## Acknowledgements

AAP was supported in part by the Natural Sciences and Engineering Research
Council of Canada, and the National Science Foundation under NSF Award 1922658
and grant DMS-2015291. JNW gratefully acknowledges the support of National
Science Foundation grant DMS-2015291.

##  Appendix A Second-order error estimate

The proof of Theorem 2 hinges on the dynamic formulations of $W_{2}^{2}(P,Q)$
and $S_{\varepsilon}(P,Q)$ (Benamou and Brenier, 2000; Chizat et al., 2020;
Conforti and Tamanini, 2021). We begin with the former:

$\displaystyle\frac{1}{2}W_{2}^{2}(P,Q)=\inf_{\rho,v}\int_{0}^{1}\int_{\mathbb{R}^{d}}\frac{1}{2}\|v(t,x)\|^{2}_{2}\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t,$
(31)

subject to $\partial_{t}\rho+\nabla\cdot(\rho v)=0$, called the continuity
equation, with $\rho(0,\cdot)=p(\cdot)$ and $\rho(1,\cdot)=q(\cdot)$. We let
$(\rho_{0},v_{0})$ denote the joint minimizers to Eq. 31 satisfying these
conditions.

Similarly, there exists a dynamic formulation for $S_{\varepsilon}$ (see
(Chizat et al., 2020; Conforti and Tamanini, 2021) for more information): for
two measures with bounded densities and compact support,

$\displaystyle S_{\varepsilon}(P,Q)+\varepsilon\log(\Lambda_{\varepsilon})$
$\displaystyle=\inf_{\rho,v}\int_{0}^{1}\int_{\mathbb{R}^{d}}\left(\frac{1}{2}\|v(t,x)\|^{2}_{2}+\frac{\varepsilon^{2}}{8}\|\nabla_{x}\log(\rho(t,x)\|^{2}_{2}\right)\rho(t,x)\,\mathrm{d}x\,\mathrm{d}t$
(32)
$\displaystyle-\frac{\varepsilon}{2}(\operatorname{Ent}(P)+\operatorname{Ent}(Q),$

subject to the same conditions as Eq. 31, where
$\Lambda_{\varepsilon}=(2\pi\varepsilon)^{d/2}$.

If we plug in the minimizers from Eq. 31 into Eq. 32, we get exactly the
result of Eq. 11 by optimality

$\displaystyle S_{\varepsilon}(P,Q)+\varepsilon\log(\Lambda_{\varepsilon})$
$\displaystyle\leq\int_{0}^{1}\int_{\mathbb{R}^{d}}\frac{1}{2}\|v_{0}(t,x)\|^{2}_{2}\rho_{0}(t,x)\,\mathrm{d}x\,\mathrm{d}t+\frac{\varepsilon^{2}}{8}I_{0}(P,Q)-\frac{\varepsilon}{2}(\operatorname{Ent}(P)+\operatorname{Ent}(Q)),$
$\displaystyle=\frac{1}{2}W_{2}^{2}(P,Q)+\frac{\varepsilon^{2}}{8}I_{0}(P,Q)-\frac{\varepsilon}{2}(\operatorname{Ent}(P)+\operatorname{Ent}(Q)),$

where we identify
$I_{0}(P,Q)=\int_{0}^{1}\int_{\mathbb{R}^{d}}\|\nabla_{x}\log\rho_{0}(t,x)\|_{2}^{2}\rho_{0}(t,x)\,\mathrm{d}x\,\mathrm{d}t$.

##  Appendix B Laplace’s method

Laplace’s method is a technique used to approximate integrals of the following
form

$\displaystyle
I(\varepsilon):=\frac{1}{\Lambda_{\varepsilon}}\int\exp\left(-\frac{1}{\varepsilon}f(x)\right)\,\mathrm{d}x\,,$
(33)

when $\varepsilon\to 0$, with $f$ convex and sufficiently regular and where
$\Lambda_{\varepsilon}=(2\pi\varepsilon)^{d/2}$. This approximation relies on
expanding $f$ around its global minimum; assuming that $f$ is twice-
differentiable, the behavior of $f$ near its minimum will be quadratic, so
that Eq. 33 will resemble a Gaussian integral for $\varepsilon$ sufficiently
small.

In what follows, we write $\mathrm{d}^{2}f(0,y),\mathrm{d}^{3}f(0,y)$ for the
second and third total derivative of $f$ at $x$, respectively. That is, for
$y\in\mathbb{R}^{d}$

$\,\mathrm{d}^{2}f(x,y):=y^{\top}\nabla^{2}f(x)y,\
\,\mathrm{d}^{3}f(x,y):=\sum_{i,j,k=1}^{d}\frac{\partial^{3}f(x)}{\partial
y_{i}\partial y_{j}\partial y_{k}}y_{i}y_{j}y_{k}\,.$

We also define the set $B_{r}(a):=\\{y\in\mathbb{R}^{d}\ |\ \|y-a\|\leq r\\}$,
for some $r>0$ and $a\in\mathbb{R}^{d}$. Recall that for a positive definite
matrix $S$, we define $J(S):=\sqrt{\det(S)}$.

######  Theorem 7.

Let $I(\varepsilon)$ be as in Eq. 33, with $f\in\mathcal{C}^{\alpha+1}$,
$m$-strongly convex, $M$-smooth, and $\alpha>1$. Assume $f$ has a global
minimum at $x^{*}$. Then there exist positive constants $c$ and $C$ depending
on $m,M,\alpha,d,\text{and }\|f\|_{\mathcal{C}^{\alpha+1}}$ such that for all
$\varepsilon\in(0,1)$,

$\displaystyle c\leq J(\nabla^{2}f(x^{*}))I(\varepsilon)\leq
1+C(\varepsilon^{(\alpha-1)/2\wedge 1})\,.$ (34)

###### Proof.

Without loss of generality, we may assume that $x^{*}=0$. For the remainder of
the proof, we let $A:=\nabla^{2}f(0)$. Let
$\tau=C_{m,M,d,\alpha}\sqrt{\log(2\varepsilon^{-1})}$, where the constant is
to be decided later. We split the desired integral into two parts:

$\displaystyle I(\varepsilon)$
$\displaystyle=\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}e^{\frac{-1}{\varepsilon}f(y)}\,\mathrm{d}y+\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)^{c}}e^{\frac{-1}{\varepsilon}f(y)}\,\mathrm{d}y=:I_{1}(\varepsilon)+I_{2}(\varepsilon)\,.$

##### Lower bounds

Note that $I_{2}(\varepsilon)\geq 0$, so it suffices to only prove
$I_{1}(\varepsilon)\geq\frac{c}{\sqrt{\det(A)}}$ for some constant $c>0$.

Since $f\in\mathcal{C}^{\alpha+1}$, we have the following Taylor expansion

$-f(y)\geq-\frac{1}{2}y^{\top}Ay-C\|y\|^{(\alpha+1)\wedge
3}\geq-\frac{M}{2}\|y\|^{2}-C\|y\|^{(\alpha+1)\wedge 3}$

for some constant $C>0$. Using this expansion, we arrive at

$\displaystyle I_{1}(\varepsilon)$
$\displaystyle=\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}\exp\left[-\frac{M}{2\varepsilon}\|y\|^{2}-\frac{C}{\varepsilon}\|y\|^{(\alpha+1)\wedge
3}\right]\,\mathrm{d}y$
$\displaystyle\geq\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}\exp\left[-\frac{M}{2\varepsilon}\|y\|^{2}-\frac{C}{\varepsilon}(\tau\sqrt{\varepsilon})^{(\alpha+1)\wedge
3}\right]\,\mathrm{d}y\,.$

Performing a change of measure and rearranging, we get

$\displaystyle J(A)I_{1}(\varepsilon)$ $\displaystyle\geq
e^{-C(\tau\sqrt{\varepsilon})^{(\alpha+1)\wedge
3}/\varepsilon}\frac{J(A)}{(2M\pi)^{d/2}}\int_{B_{\tau\sqrt{M}(0)}}e^{-\frac{1}{2}\|y\|^{2}}\,\mathrm{d}y$
$\displaystyle\gtrsim e^{-C(\tau\sqrt{\varepsilon})^{(\alpha+1)\wedge
3}/\varepsilon}J(A)\mathbb{P}(\|Y\|\leq\tau\sqrt{M})\,,$

where $Y\sim N(0,I_{d})$. Since $\alpha>1$, the quantity
$C(\tau\sqrt{\varepsilon})^{(\alpha+1)\wedge 3}/\varepsilon$ is bounded as
$\varepsilon\to 0$, so we may bound
$e^{-C(\tau\sqrt{\varepsilon})^{(\alpha+1)\wedge 3}/\varepsilon}$ from below
by a constant. Since $J(A)$ and $\mathbb{P}(\|Y\|\leq\tau\sqrt{M})$ are both
also bounded from below, we obtain that $J(A)I_{1}(\varepsilon)\geq c>0$, as
desired.

##### Upper bounds

We first show that the contribution from $I_{2}(\varepsilon)$ is negligible.
The strong convexity of $f$ implies

$f\geq\frac{m}{2}\|y\|^{2},$

leading us to the upper bound

$\displaystyle I_{2}(\varepsilon)$
$\displaystyle\leq\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)^{c}}e^{-\frac{m}{2\varepsilon}\|y\|^{2}}\,\mathrm{d}y$
$\displaystyle=\frac{1}{(2m\pi)^{d/2}}\int_{B_{\tau}(0)^{c}}e^{-\frac{1}{2}\|y\|^{2}}\,\mathrm{d}y$
$\displaystyle\leq\frac{1}{(2m\pi)^{d/2}}e^{-\frac{1}{4}\tau^{2}}\int
e^{-\frac{1}{4}\|y\|^{2}}\,\mathrm{d}y$ $\displaystyle\lesssim
e^{-\frac{1}{4}\tau^{2}}\,,$

where in the penultimate inequality we have used the fact that
$e^{-\frac{1}{4}\|y\|^{2}}\leq e^{-\frac{1}{4}\tau^{2}}$ on $B_{\tau}(0)^{c}$.
Taking $C_{m,M,d,\alpha}$ sufficiently large in the definition of $\tau$, we
can make this term smaller than $\varepsilon$.

For upper bounds on $I_{1}(\varepsilon)$, we proceed in a similar fashion. If
$f\in\mathcal{C}^{\alpha+1}$ for $\alpha\in(1,2]$, then we employ the bound

$-f(y)\leq-\frac{1}{2}y^{\top}Ay+C\|y\|^{\alpha+1}\,,$

yielding

$\displaystyle I_{1}(\varepsilon)$
$\displaystyle=\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}e^{-\frac{1}{\varepsilon}f(y)}\,\mathrm{d}y\leq\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}e^{-\frac{1}{2\varepsilon}y^{\top}Ay+\frac{C}{\varepsilon}\|y\|^{\alpha+1}}\,\mathrm{d}y\,.$

Performing the change of variables $u=\sqrt{1/\varepsilon}y$, we arrive at

$\displaystyle
I_{1}(\varepsilon)\leq\frac{1}{(2\pi)^{d/2}}\int_{B_{\tau}(0)}e^{\frac{-1}{2}u^{\top}Au}e^{C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}}\,\mathrm{d}u$

Since $\alpha>1$, the term $C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}$ is
bounded above on $B_{\tau}(0)$, so that there exists a positive constant
$C^{\prime}$ such that

$e^{C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}}\leq
1+C^{\prime}\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}\quad\forall u\in
B_{\tau}(0)\,.$

We obtain

$\displaystyle I_{1}(\varepsilon)$
$\displaystyle\leq\frac{1}{(2\pi)^{d/2}}\int_{B_{\tau}(0)}e^{\frac{-1}{2}u^{\top}Au}(1+C^{\prime}\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1})\,\mathrm{d}u$
$\displaystyle\leq\frac{1}{(2\pi)^{d/2}}\int
e^{\frac{-1}{2}u^{\top}Au}(1+C^{\prime}\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1})\,\mathrm{d}u.$

Performing another change of variables yields

$I_{1}(\varepsilon)\leq\frac{1}{(2\pi)^{d/2}J(A)}\int(1+C^{\prime}\varepsilon^{(\alpha-1)/2}\|A^{-1/2}u\|^{\alpha+1})e^{-\frac{1}{2}\|u\|^{2}}\,\mathrm{d}u$

We obtain

$J(A)I_{1}(\varepsilon)\leq 1+C^{\prime\prime}\varepsilon^{(\alpha-1)/2}\,.$

Combining this with the bound on $J(A)I_{2}(\varepsilon)$ yields the bound for
$\alpha\leq 2$.

When $\alpha>2$, we use the same technique but expand to the third order,
yielding

$\displaystyle I_{1}(\varepsilon)$
$\displaystyle=\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}e^{-\frac{1}{\varepsilon}f(y)}\,\mathrm{d}y$
$\displaystyle\leq\frac{1}{\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(0)}e^{-\frac{1}{2\varepsilon}y^{\top}Ay-\frac{1}{6\varepsilon}\,\mathrm{d}^{3}f(0,y)+\frac{C}{\varepsilon}\|y\|^{\alpha+1}}\,\mathrm{d}y$
$\displaystyle=\frac{1}{(2\pi)^{d/2}}\int_{B_{\tau}(0)}e^{-\frac{1}{2}u^{\top}Au-\frac{\varepsilon^{1/2}}{6}\,\mathrm{d}^{3}f(0,u)+C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}}\,\mathrm{d}u$

Since
$-\frac{\varepsilon^{1/2}}{6}\,\mathrm{d}^{3}f(0,u)+C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}$
is bounded on $B_{\tau}(0)$, we have

$e^{-\frac{\varepsilon^{1/2}}{6}\,\mathrm{d}^{3}f(0,u)+C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}}\leq
1-\frac{\varepsilon^{1/2}}{6}\,\mathrm{d}^{3}f(0,u)+C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}+R(u)\,,$

where $R$ is a positive remainder term satisfying
$R(u)\lesssim\varepsilon(\,\mathrm{d}^{3}f(0,u))^{2}+\varepsilon^{\alpha-1}\|u\|^{2(\alpha+1)}$.
We obtain

$I_{1}(\varepsilon)\leq\frac{1}{(2\pi)^{d/2}}\int_{B_{\tau}(0)}\left(1-\frac{\varepsilon^{1/2}}{6}\,\mathrm{d}^{3}f(0,u)+C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}+R(u)\right)e^{-\frac{1}{2}u^{\top}Au}\,\mathrm{d}u\,.$

The symmetry of $B_{\tau}(0)$ and the fact that
$\,\mathrm{d}^{3}f(0,u)e^{-\frac{1}{2}u^{\top}Au}$ is an odd function of $u$
imply

$\int_{B_{\tau}(0)}\,\mathrm{d}^{3}f(0,u)e^{-\frac{1}{2}u^{\top}Au}\,\mathrm{d}u=0\,,$

so

$\displaystyle I_{1}(\varepsilon)$
$\displaystyle\leq\frac{1}{(2\pi)^{d/2}}\int(1+C\varepsilon^{(\alpha-1)/2}\|u\|^{\alpha+1}+R(u))e^{-\frac{1}{2}u^{\top}Au}\,\mathrm{d}u$
$\displaystyle=\frac{1}{(2\pi)^{d/2}J(A)}\int(1+C\varepsilon^{(\alpha-1)/2}\|A^{-1/2}u\|^{\alpha+1}+R(A^{-1/2}u))e^{-\frac{1}{2}\|u\|^{2}}\,\mathrm{d}u$
$\displaystyle\leq
1+C^{\prime\prime}\varepsilon^{(\alpha-1)/2}+C^{\prime\prime}\varepsilon\,,$

which is the desired bound. ∎

######  Corollary 2.

Assume (A1) to (A3). For all $\alpha\in(1,3]$, there exist positive constants
$c$ and $C$ such that

$\displaystyle c\leq J(\nabla^{2}\varphi_{0}^{*}(x^{*}))Z_{\varepsilon}(x)\leq
1+C\varepsilon^{(\alpha-1)/2}\,,$ (35)

for all $\varepsilon\in(0,1)$ and $x\in\operatorname{\mathrm{supp}}(P)$.

###### Proof.

Take $f(\cdot)=D[\cdot|x^{*}]$ which is $1/L$-strongly convex, and
$1/\mu$-smooth, with minimizer $x^{*}$ (see Eq. 18). The claim now follows
from Theorem 7. ∎

##  Appendix C Omitted proofs

###  C.1 Proof of Proposition 1

It suffices to show that

$S_{\varepsilon}(P,Q)\geq\sup_{\eta\in
L^{1}(\pi_{\varepsilon})}\int\eta\,\mathrm{d}\pi_{\varepsilon}-\varepsilon\iint
e^{(\eta(x,y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}P(x)\,\mathrm{d}Q(y)+\varepsilon\,,$

since the other direction follows from choosing $\eta(x,y)=f(x)+g(y)$ and
using Eq. 8.

Write

$\gamma(x,y)=e^{\frac{1}{\varepsilon}(f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})}$

for the $P\otimes Q$ density of $\pi_{\varepsilon}$. The inequality

$a\log a\geq ab-e^{b}+a$

holds for all $a\geq 0$ and $b\in\mathbb{R}$, as can be seen by noting that
the right side is a concave function of $b$ which achieves its maximum at
$b=\log a$. Applying this inequality with $a=\gamma(x,y)$ and
$b=\frac{1}{\varepsilon}(\eta(x,y)-\frac{1}{2}\|x-y\|^{2})$ and integrating
with respect to $P\otimes Q$ yields

$\displaystyle\int\log\gamma\,\mathrm{d}\pi_{\varepsilon}\geq\frac{1}{\varepsilon}\bigg{(}\int\eta\,\mathrm{d}\pi_{\varepsilon}$
$\displaystyle-\int\frac{1}{2}\|x-y\|^{2}\,\mathrm{d}\pi_{\varepsilon}(x,y)\bigg{)}$
$\displaystyle-\iint
e^{(\eta(x,y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}P(x)\,\mathrm{d}Q(y)+1$

Multiplying by $\varepsilon$ and using the fact that

$\displaystyle\int\varepsilon\log\gamma\,\mathrm{d}\pi_{\varepsilon}$
$\displaystyle=\int(f_{\varepsilon}(x)+g_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})\,\mathrm{d}\pi_{\varepsilon}=S_{\varepsilon}(P,Q)-\int\frac{1}{2}\|x-y\|^{2}\,\mathrm{d}\pi_{\varepsilon}(x,y)$

yields the claim. ∎

###  C.2 Proof of Proposition 4

Proposition 4 follows from the following more general result, which recovers
Proposition 4 by choosing $\tilde{P}=P_{n}$.

######  Proposition 5.

Let $P$ and $Q$ be probability measures with support contained in $\Omega$,
and denote by $P_{n}$ and $Q_{n}$ corresponding empirical measures. If
$\tilde{P}$ is a probability measure with support in $\Omega$ such that
$\mathrm{d_{TV}}(\tilde{P},P_{n})\leq\delta$ for some $\delta\geq 0$, then

$\displaystyle\mathbb{E}\Big{\\{}\sup_{\chi:\Omega\times\Omega\to\mathbb{R}}\iint\chi(x,y)\,\mathrm{d}\pi_{\varepsilon,n}(x,y)$
$\displaystyle-\iint(e^{\chi(x,y)}-1)\tilde{\gamma}(x,y)\,\mathrm{d}P(x)\,\mathrm{d}Q_{n}(y)\Big{\\}}$
$\displaystyle\lesssim\varepsilon^{-1}\delta+(\varepsilon^{-1}+\varepsilon^{-d^{\prime}/2})\log(n)n^{-1/2}\,,$

where $\pi_{\varepsilon,n}$ is the optimal entropic plan for $P$ and $Q_{n}$
and $\tilde{\gamma}$ is the $\tilde{P}\otimes Q_{n}$ density of the optimal
entropic plan for $\tilde{P}$ and $Q_{n}$.

###### Proof.

Write $\tilde{f}_{\varepsilon}$ and $\tilde{g}_{\varepsilon}$ for the optimal
entropic potentials for $\tilde{P}$ and $Q_{n}$, so that
$\tilde{\gamma}(x,y)=\exp\left(\varepsilon^{-1}(\tilde{f}_{\varepsilon}(x)+\tilde{g}_{\varepsilon}(y)-\frac{1}{2}\|x-y\|^{2})\right)$.
Plugging in
$\eta(x,y)=\varepsilon\chi(x,y)+\tilde{f}_{\varepsilon}(x)+\tilde{g}_{\varepsilon}(y)$
into Proposition 1 gives

$\displaystyle\sup_{\chi:\Omega\times\Omega\to\mathbb{R}}\iint\chi\,\mathrm{d}\pi_{\varepsilon,n}-\iint(e^{\chi(x,y)}-1)\tilde{\gamma}(x,y)\,\mathrm{d}P(x)\,\mathrm{d}Q_{n}(y)$
$\displaystyle\leq\varepsilon^{-1}\bigg{(}S_{\varepsilon}(P,Q_{n})$
$\displaystyle-\int\tilde{f}_{\varepsilon}\,\mathrm{d}P-\int\tilde{g}_{\varepsilon}\,\mathrm{d}Q_{n}\bigg{)}\,,$

where we have used that $\tilde{\gamma}$ is a probability density with respect
to $P\otimes Q_{n}$ by Eq. 9.

Let $f_{\varepsilon,n}$ and $g_{\varepsilon,n}$ be the optimal entropic dual
potentials for $P$ and $Q_{n}$. As in the proof of Lemma 7, the optimality of
$\tilde{f}_{\varepsilon}$ and $\tilde{g}_{\varepsilon}$ for the pair
$(\tilde{P},Q_{n})$ implies

$\displaystyle\int\tilde{f}_{\varepsilon}\,\mathrm{d}\tilde{P}+\int\tilde{g}_{\varepsilon}\,\mathrm{d}Q_{n}$
$\displaystyle\geq\int f_{\varepsilon,n}\,\mathrm{d}\tilde{P}+\int
g_{\varepsilon,n}\,\mathrm{d}Q_{n}$ $\displaystyle-\varepsilon\iint
e^{\frac{1}{\varepsilon}(f_{\varepsilon,n}(x)+g_{\varepsilon,n}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}\tilde{P}(x)\,\mathrm{d}Q_{n}(y)+\varepsilon$
$\displaystyle=\int f_{\varepsilon,n}\,\mathrm{d}\tilde{P}+\int
g_{\varepsilon,n}\,\mathrm{d}Q_{n}\,,$

since $\int
e^{\frac{1}{\varepsilon}(f_{\varepsilon,n}(x)+g_{\varepsilon,n}(y)-\frac{1}{2}\|x-y\|^{2})}\,\mathrm{d}Q_{n}(y)=1$
by the dual optimality condition in Eq. 9. Therefore

$\displaystyle
S_{\varepsilon}(P,Q_{n})-\int\tilde{f}_{\varepsilon}\,\mathrm{d}P-\int\tilde{g}_{\varepsilon}\,\mathrm{d}Q_{n}$
$\displaystyle\leq\int(f_{\varepsilon,n}-\tilde{f}_{\varepsilon})(\mathrm{d}P-\mathrm{d}\tilde{P})(y)\,$
$\displaystyle=\int(f_{\varepsilon,n}-\tilde{f}_{\varepsilon})(\mathrm{d}P-\mathrm{d}{P}_{n})(y)$
$\displaystyle\qquad+\int(f_{\varepsilon,n}-\tilde{f}_{\varepsilon})(\mathrm{d}P_{n}-\mathrm{d}\tilde{P})(y)$

By Genevay et al. (2019, Proposition 1), we may choose $f_{\varepsilon,n}$ and
$\tilde{f}_{\varepsilon}$ to satisfy
$\left\|f_{\varepsilon,n}\right\|_{\infty},\|\tilde{f}_{\varepsilon}\|_{\infty}\lesssim
1$, so we may bound the second term as

$\int(f_{\varepsilon,n}-\tilde{f}_{\varepsilon})(\mathrm{d}P_{n}-\mathrm{d}\tilde{P}_{n})(y)\lesssim\mathrm{d_{TV}}(\tilde{P}_{n},P_{n})\leq\delta\,.$

Also, since $f_{\varepsilon,n}$ is independent of $P_{n}$,

$\mathbb{E}f_{\varepsilon,n}(\mathrm{d}P-\mathrm{d}P_{n})(y)=0\,.$

Altogether, we obtain

$\mathbb{E}\sup_{\chi:\Omega\times\Omega\to\mathbb{R}}\iint\chi\,\mathrm{d}\pi_{\varepsilon,n}-\iint(e^{\chi(x,y)}-1)\gamma(x,y)\,\mathrm{d}P(x)\,\mathrm{d}Q_{n}(y)\\\
\lesssim\varepsilon^{-1}\bigg{(}\delta+\mathbb{E}\int\tilde{f}_{\varepsilon}(\mathrm{d}P_{n}-\mathrm{d}P)(y)\bigg{)}\,.$

We conclude by again appealing to Genevay et al. (2019, Proposition 1): since
$\tilde{f}_{\varepsilon}$ is an optimal entropic potential for the pair of
compactly distributed probability measures $(\tilde{P},Q_{n})$, its
derivatives up to order $s$ are bounded by $C_{s,d,K}(1+\varepsilon^{1-s})$ on
any compact set $K$ for any $s\geq 0$. Taking $K$ to be a suitably large ball
containing $\Omega$ and applying Lemma 7 with $s=2\lceil d/2\rceil$ yields the
claim. ∎

###  C.3 Proofs from Section 4

###### Proof of Lemma 2.

Fix $x\in\operatorname{\mathrm{supp}}(P)$ and let $x^{*}:=T_{0}(x)$, and for
notational convenience, write $Y$ for the random variable with density
$q_{\varepsilon}^{x}$, and denote its mean by $\bar{Y}$. It suffices to show
the existence of a constant $K$ such that for any unit vector $v$,

$\mathbb{E}e^{(v^{\top}(Y-x^{*}))^{2}/4L\varepsilon}\leq K\,.$ (36)

Indeed, by Young’s and Jensen’s inequalities, this implies

$\mathbb{E}e^{(v^{\top}(Y-\bar{Y}))^{2}/8L\varepsilon}\leq
e^{(v^{\top}(\bar{Y}-x^{*}))^{2}/4L\varepsilon}\mathbb{E}e^{(v^{\top}(Y-x^{*}))^{2}/4L\varepsilon}\leq
K^{2}\,,$

and hence by another application of Jensen’s inequality that

$\mathbb{E}e^{(v^{\top}(Y-\bar{Y}))^{2}/C\varepsilon}\leq 2$

for $C=8LK^{2}$.

We prove Eq. 36 using the strong convexity of $D[y|x^{*}]$. By Eq. 18,

$\displaystyle\mathbb{E}e^{(v^{\top}(Y-x^{*}))^{2}/4L\varepsilon}$
$\displaystyle=\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int
e^{-\frac{1}{\varepsilon}D[y|x^{*}]+\frac{1}{4L\varepsilon}\|y-x^{*}\|^{2}}\,\mathrm{d}y$
$\displaystyle\leq\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int
e^{-\frac{1}{4L\varepsilon}\|y-x^{*}\|^{2}}\,\mathrm{d}y$
$\displaystyle=\frac{(2L)^{d/2}}{Z_{\varepsilon}(x)}$ $\displaystyle\lesssim
1\,,$

where the final inequality uses Corollary 2.

∎

###### Proof of Lemma 3.

Let us first fix an $x\in\operatorname{\mathrm{supp}}(P)$, and write $Y$ for
the random variable with density $q_{\varepsilon}^{x}$ and $\bar{Y}$ for its
mean, and write $x^{*}:=T_{0}(x)$. Lemma 2 implies (see Vershynin, 2018,
Proposition 2.5.2) that there exists a positive constant $C$, independent of
$x$, such for any $v\in\mathbb{R}^{d}$,

$\mathbb{E}e^{(v^{\top}(Y-x^{*}))}=e^{v^{\top}(\bar{Y}-x^{*})}\mathbb{E}e^{(v^{\top}(Y-\bar{Y}))}\leq
e^{v^{\top}(\bar{Y}-x^{*})+C\varepsilon\|v\|^{2}}\leq
e^{\frac{1}{4\varepsilon}\|\bar{Y}-x^{*}\|^{2}+(C+1)\varepsilon\|v\|^{2}}\,,$

where the last step uses Young’s inequality. Equivalently, for
$a>(C+1)\varepsilon$, we have for all $x\in\operatorname{\mathrm{supp}}(P)$
and $v\in\mathbb{R}^{d}$

$\int
e^{(v^{\top}(y-x^{*}))-a\|v\|^{2}}q_{\varepsilon}^{x}(y)\,\mathrm{d}y\leq
e^{\frac{1}{4\varepsilon}\|\bar{y}^{x}-x^{*}\|^{2}}$

Applying this inequality with $v=h(x)$ and integrating with respect to $P$
yields the claim. ∎

###### Proof of Lemma 4.

It suffices to prove the claim for $\alpha\in(1,2]$. Let us fix an
$x\in\operatorname{\mathrm{supp}}(P)$. Since
$\varphi_{0}^{*}\in\mathcal{C}^{\alpha+1}(\Omega)$, Taylor’s theorem implies

$D[y|x^{*}]=-x^{\top}y+\varphi_{0}(x)+\varphi_{0}^{*}(y)=\frac{1}{2}(y-x^{*})^{\top}\nabla^{2}\varphi_{0}^{*}(x^{*})(y-x^{*})+R(y|x^{*})\,,$

where the remainder satisfies

$|R(y|x^{*})|\lesssim\|y-x^{*}\|^{1+\alpha}\,.$ (37)

We aim to bound

$\displaystyle\left\|\bar{y}^{x}-x^{*}\right\|$
$\displaystyle=\left\|\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int(y-x^{*})e^{-\frac{1}{\varepsilon}D[y|x^{*}]}\,\mathrm{d}y\right\|$

Let $\tau=C\sqrt{\log(2\varepsilon^{-1})}$ for a sufficiently large constant
$C$. As in the proof of Theorem 7, the contribution to the integral from the
set $B_{\tau\sqrt{\varepsilon}}(x^{*})^{c}$ is negligible; indeed, Eq. 18
implies

$\displaystyle\left\|\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(x^{*})^{c}}(y-x^{*})e^{-\frac{1}{\varepsilon}D[y|x^{*}]}\,\mathrm{d}y\right\|$
$\displaystyle\qquad\leq\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(x^{*})^{c}}\|y-x^{*}\|e^{-\frac{1}{2L\varepsilon}\|y-x^{*}\|^{2}}\,\mathrm{d}y$
$\displaystyle\qquad=\frac{\varepsilon^{(d+1)/2}}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int_{B_{\tau}(0)^{c}}\|y\|e^{-\frac{1}{2L}\|y\|^{2}}\,\mathrm{d}y$
$\displaystyle\qquad\leq\frac{\varepsilon^{(d+1)/2}}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\left(\int\|y\|^{2}e^{-\frac{1}{2L}\|y\|^{2}}\,\mathrm{d}y\right)^{1/2}\left(\int_{B_{\tau}(0)^{c}}e^{-\frac{1}{2L}\|y\|^{2}}\,\mathrm{d}y\right)^{1/2}$
$\displaystyle\qquad\lesssim\varepsilon^{-1/2}\mathbb{P}[\|Y\|\geq\tau]\,,\quad\quad
Y\sim\mathcal{N}(0,I_{d})\,,$

and this quantity can be made smaller than $\varepsilon$ by choosing the
constant in the definition of $\tau$ sufficiently large.

It remains to bound

$\left\|\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(x^{*})}(y-x^{*})e^{-\frac{1}{\varepsilon}R(y|x^{*})}e^{-\frac{1}{2\varepsilon}(y-x^{*})^{\top}\nabla^{2}\varphi_{0}^{*}(x^{*})(y-x^{*})}\,\mathrm{d}y\right\|=\\\
\left\|\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}}\int_{B_{\tau\sqrt{\varepsilon}}(x^{*})}(y-x^{*})\left(e^{-\frac{1}{\varepsilon}R(y|x^{*})}-1\right)e^{-\frac{1}{2\varepsilon}(y-x^{*})^{\top}\nabla^{2}\varphi_{0}^{*}(x^{*})(y-x^{*})}\,\mathrm{d}y\right\|\,,$
(38)

where we have used that

$\int_{B_{\tau\sqrt{\varepsilon}}(x^{*})}(y-x^{*})e^{-\frac{1}{2\varepsilon}(y-x^{*})^{\top}\nabla^{2}\varphi_{0}^{*}(x^{*})(y-x^{*})}\,\mathrm{d}y=0\,.$

By (37),

$\frac{1}{\varepsilon}|R(y|x^{*})|\lesssim\frac{1}{\varepsilon}\|y-x^{*}\|^{1+\alpha}\lesssim
1\quad\quad\forall y\in B_{\tau\sqrt{\varepsilon}}(x^{*})\,,$

and since $|e^{t}-1|\lesssim|t|$ for $|t|\lesssim 1$, we obtain that

$\left|e^{-\frac{1}{\varepsilon}R(y|x^{*})}-1\right|\lesssim\frac{1}{\varepsilon}\|y-x^{*}\|^{1+\alpha}\quad\quad\forall
y\in B_{\tau\sqrt{\varepsilon}}(x^{*})\,.$

Therefore (38) is bounded above by

$\displaystyle\frac{C}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}\varepsilon}\int_{\mathbb{R}^{d}}\|y-x^{*}\|^{2+\alpha}e^{-\frac{1}{2\varepsilon}(y-x^{*})^{\top}\nabla^{2}\varphi_{0}^{*}(x^{*})(y-x^{*})}\,\mathrm{d}y\lesssim\varepsilon^{\alpha/2}\,,$

where in the last step we have applied Corollary 2. We therefore obtain

$\|\mathbb{E}_{q_{\varepsilon}^{x}}(Y)-x^{*}\|\
\lesssim\varepsilon^{\alpha/2}+\varepsilon\,.$

Taking squares, we get the desired result. ∎

##  Appendix D Supplementary results

######  Proposition 6.

For any $x\in\operatorname{\mathrm{supp}}(P)$, if $a\in[L\varepsilon,1]$, then

$\mathbb{E}\sup_{h:\Omega\to\mathbb{R}^{d}}\int
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}(\mathrm{d}Q_{n}-\mathrm{d}Q)(y)\lesssim(1+\varepsilon^{-d/2})n^{-1/2}\,,$

where the implicit constant is uniform in $x$.

###### Proof.

To bound this process, we employ the following two lemmas:

######  Lemma 5.

If $a\geq L\varepsilon$, then for any $v\in\mathbb{R}^{d}$,

$v^{\top}(y-x^{*})-a\|v\|^{2}-\frac{1}{\varepsilon}D[y|x^{*}]\leq-\frac{\varepsilon
L}{2}\|v\|^{2}\,.$

###### Proof.

By Eq. 18, $D[y|x^{*}]\geq\frac{1}{2L}\|y-x^{*}\|^{2}$. Combining this fact
with Young’s inequality yields

$\displaystyle v^{\top}(y-x^{*})-a\|v\|^{2}-\frac{1}{\varepsilon}D[y|x^{*}]$
$\displaystyle\leq\frac{\varepsilon L}{2}\|v\|^{2}+\frac{1}{2\varepsilon
L}\|y-x^{*}\|^{2}-a\|v\|^{2}-\frac{1}{\varepsilon}D[y|x^{*}]$
$\displaystyle\leq-\frac{\varepsilon L}{2}\|v\|^{2}\,,$

as claimed. ∎

By slight abuse of notation, for any $v\in\mathbb{R}^{d}$, write
$j_{v}:\Omega\to\mathbb{R}$ for the function

$j_{v}(y)=v^{\top}(y-T_{0}(x))-a\|v\|^{2}\,.$

Let
$\mathcal{J}=\\{e^{j_{v}}\frac{q_{\varepsilon}^{x}(y)}{q(y)}:v\in\mathbb{R}^{d}\\}$.

######  Lemma 6.

If $a\in[L\varepsilon,1]$, then

$\log N(\tau,\mathcal{J},\|\cdot\|_{L^{\infty}(Q)})\lesssim d\log(K/\tau)\,,$

where $K\lesssim(1+\varepsilon^{-d/2})$.

###### Proof.

Fix $\delta\in(0,1)$. Let $\mathcal{N}_{\delta}$ be a $\delta^{3/2}$-net with
respect to the Euclidean metric of a ball of radius $\delta^{-1/2}$ in
$\mathbb{R}^{d}$, and consider the set

$\mathcal{G}_{\delta}\coloneqq\left\\{e^{j_{v}}\frac{q_{\varepsilon}^{x}(y)}{q(y)}:v\in\mathcal{N}_{\delta}\right\\}\cup\\{e^{j_{w}}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\\}\,,$

where $w\in\mathbb{R}^{d}$ is an arbitrary vector of norm $\delta^{-1/2}$. By
Lemma 5, if $a>L\varepsilon$ and $\|v\|\geq R$, then

$\sup_{y\in\operatorname{\mathrm{supp}}(Q)}e^{j_{v}}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\leq\sup_{y\in\operatorname{\mathrm{supp}}(Q)}\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}q(y)}e^{-\frac{1}{2\varepsilon
L}R^{2}}\leq\sup_{y\in\operatorname{\mathrm{supp}}(Q)}\frac{2L}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}q(y)R^{2}}\,.$

Therefore, if $v\in\mathbb{R}^{d}$ satisfies $\|v\|\geq\delta^{-1/2}$, then

$\sup_{y\in\operatorname{\mathrm{supp}}(Q)}\left|e^{j_{w}(y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}-e^{j_{v}(y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\right|\leq\sup_{y\in\operatorname{\mathrm{supp}}(Q)}\frac{4\delta
L}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}q(y)}\leq K\delta$

for
$K=\sup_{y\in\operatorname{\mathrm{supp}}(Q)}\frac{1+4L}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}q(y)}\lesssim\varepsilon^{-d/2}$.

On the other hand, if $v\in\mathbb{R}^{d}$ satisfies $\|v\|\leq\delta^{-1/2}$,
pick $u\in\mathcal{N}_{\delta}$ satisfying $\|u-v\|\leq\delta^{3/2}$. We then
have

$\left|e^{j_{u}(y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}-e^{j_{v}(y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}\right|\leq\frac{|j_{u}(y)-j_{v}(y)|}{q(y)}\,,$

where we have used Lemma 5 combined with the inequality

$|e^{a}-e^{b}|\leq|a-b|\quad\quad\forall a,b\leq 0\,.$

Since $\|u-v\|\leq\delta^{3/2}$ and $\|u\|+\|v\|\leq 2\delta^{-1/2}$, we have
for any $y\in\Omega$,

$|j_{u}(y)-j_{v}(y)|=|(u-v)^{\top}(y-T_{0}(x))-a(\|u\|^{2}-\|v\|^{2})|\lesssim\delta^{3/2}+\delta
a\,,$

where we have used the fact that $y$ and $T_{0}(x)$ lie in the compact set
$\Omega$. Therefore, as long as $a\leq 1$, this quantity is bounded by
$C\delta$ for a positive constant $C$.

All told, we obtain that for any $v\in\mathbb{R}^{d}$, there exists a
$g\in\mathcal{G}_{\delta}$ such that

$\left\|e^{j_{v}}\frac{q_{\varepsilon}^{x}(y)}{q(y)}-g\right\|_{L^{\infty}(Q)}\lesssim
K\delta\,,$

where $K\lesssim 1+\varepsilon^{-d/2}$. Moreover, Lemma 5 implies that, for
any $g\in\mathcal{G}_{\delta}$,

$\|g\|_{L^{\infty}(Q)}\leq\sup_{y\in\operatorname{\mathrm{supp}}(Q)}\frac{1}{Z_{\varepsilon}(x)\Lambda_{\varepsilon}q(y)}\leq
K.$

By a volume argument, we may choose $\mathcal{N}_{\delta}$ such that it
satisfies

$\log|\mathcal{N}_{\delta}|\lesssim\log(1/(\delta)).$

We therefore obtain for any $\tau\leq K$,

$\log
N(\tau,\mathcal{J},\|\cdot\|_{L^{\infty}(Q)})\leq\log|\mathcal{G}_{\tau/K}|\lesssim\log(K/\tau)\,,$

as claimed. ∎

Returning to the empirical process, we obtain by a chaining bound (Giné and
Nickl, 2016, Theorem 3.5.1)

$\displaystyle\mathbb{E}\sup_{h:\Omega\to\mathbb{R}^{d}}\int
e^{j_{h}(x,y)}\frac{q_{\varepsilon}^{x}(y)}{q(y)}(\mathrm{d}Q_{n}-\mathrm{d}Q)(y)$
$\displaystyle=\mathbb{E}\sup_{j\in\mathcal{J}}\int
j(y)(\mathrm{d}Q_{n}-\mathrm{d}Q)(y)$ $\displaystyle\lesssim
n^{-1/2}\int_{0}^{K}\sqrt{\log(K/\tau)}\,\mathrm{d}\tau$
$\displaystyle\lesssim Kn^{-1/2}\,.$

Recalling that $K\lesssim(1+\varepsilon^{-d/2})$ completes the proof. ∎

######  Lemma 7.

For a convex, compact $K\subseteq\mathbb{R}^{d}$, integer $s\geq d/2$, and
$M>0$, let $\mathcal{C}^{s}(K;M)$ be the set of $s$-Hölder smooth functions
with Hölder norm bounded by $M$. For any probability measure $\nu$ with
support contained in $K$ and corresponding empirical measure $\nu_{n}$, we
have that

$\mathbb{E}\sup_{g\in\mathcal{C}^{s}(K;M)}\int
g(y)(\mathrm{d}\nu_{n}(y)-\mathrm{d}\nu(y))\lesssim C_{K}M\log(n)n^{-1/2}\,.$

###### Proof.

We write $\mathcal{F}$ to be the set of functions in $\mathcal{C}^{s}(K;1)$. A
version of Dudley’s chaining bound (see, e.g., von Luxburg and Bousquet,
2003/04, Theorem 16) therefore implies for any $\delta\geq 0$,

$\displaystyle\mathbb{E}\sup_{g\in\mathcal{C}^{s}(K;M)}\int
g(y)(\mathrm{d}\nu_{n}(y)-\mathrm{d}\nu(y))$ $\displaystyle\lesssim
M\left(\delta+n^{-1/2}\int_{\delta}^{1}\sqrt{\log
N(\tau,\mathcal{F},\|\cdot\|_{\infty})}\,\mathrm{d}\tau\right)\,.$

Letting $s\geq d/2$ and applying standard covering number bounds for Hölder
spaces (van der Vaart and Wellner, 1996, Theorem 2.7.1) implies

$\mathbb{E}\sup_{g\in\mathcal{C}^{s}(K;M)}\int
g(y)(\mathrm{d}\nu_{n}(y)-\mathrm{d}\nu(y))\lesssim C_{K}\inf_{\delta\geq
0}M\left(\delta+n^{-1/2}\int_{\delta}^{1}\tau^{-1}\,\mathrm{d}\tau\right)\,.$

Taking $\delta=n^{-1/2}$ yields

$\mathbb{E}\sup_{g\in\mathcal{C}^{s}(K;M)}\int
g(y)(\mathrm{d}\nu_{n}(y)-\mathrm{d}\nu(y))\lesssim
C_{K}Mn^{-1/2}(1-\log(n^{-1/2}))\lesssim C_{K}Mn^{-1/2}\log n\,,$

as claimed. ∎

######  Corollary 3.

If $P$ and $Q$ are compactly supported, then

$\displaystyle\mathbb{E}S_{\varepsilon}(P,Q_{n})-S_{\varepsilon}(P,Q)\lesssim(1+\varepsilon^{1-d^{\prime}/2})\log(n)n^{-1/2}\,.$

###### Proof.

Let $(f_{\varepsilon,n},g_{\varepsilon_{n}})$ be the optimal dual potentials
for $P$ and $Q_{n}$. Following (Mena and Niles-Weed, 2019, Proposition 2),
observe that

$\displaystyle S_{\varepsilon}(\mu,\nu_{n})-S_{\varepsilon}(\mu,\nu)$
$\displaystyle=\int{f}_{(\varepsilon,n)}\,\mathrm{d}\mu+\int{g}_{(\varepsilon,n)}\,\mathrm{d}\nu_{n}-\sup_{f,g}\Big{\\{}\int
f\,\mathrm{d}\mu+\int g\,\mathrm{d}\nu$
$\displaystyle\quad\quad-\varepsilon\iint
e^{(f(x)+g(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}\mu(x)\,\mathrm{d}\nu(y)+\varepsilon\Big{\\}}$
$\displaystyle\leq\int{g}_{(\varepsilon,n)}(\mathrm{d}\nu_{n}(y)-\mathrm{d}\nu(y))\,,$

where the bound follows from choosing
$({f}_{(\varepsilon,n)},{g}_{(\varepsilon,n)})$ in the supremum and using

$\int
e^{({f}_{(\varepsilon,n)}(x)+{g}_{(\varepsilon,n)}(y)-\frac{1}{2}\|x-y\|^{2})/\varepsilon}\,\mathrm{d}\mu(x)=1\quad\forall
y\in\mathbb{R}^{d}$

by the dual optimality condition Eq. 9.

We conclude by applying Genevay et al. (2019, Proposition 1): the derivatives
of $g_{\varepsilon,n}$ up to order $s$ are bounded by
$C_{s,d,K}(1+\varepsilon^{1-s})$ on any compact set $K$ for any $s\geq 0$, so
we may take $K$ to be a suitably large ball containing the support of $P$ and
$Q$ and apply Lemma 7 with $s=2\lceil d/2\rceil=:d^{\prime}$. ∎

## References

  * Altschuler et al. (2017) Altschuler, J., Weed, J., and Rigollet, P.  Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.  In Guyon, I., von Luxburg, U., Bengio, S., Wallach, H. M., Fergus, R., Vishwanathan, S. V. N., and Garnett, R., editors, _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, 4-9 December 2017, Long Beach, CA, USA_ , pages 1961–1971, 2017.
  * Altschuler et al. (2019) Altschuler, J., Bach, F., Rudi, A., and Niles-Weed, J.  Massively scalable sinkhorn distances via the nyström method.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Altschuler et al. (2021) Altschuler, J. M., Niles-Weed, J., and Stromme, A. J.  Asymptotics for semi-discrete entropic optimal transport.  _arXiv preprint arXiv:2106.11862_ , 2021.
  * Ambrosio et al. (2008) Ambrosio, L., Gigli, N., and Savaré, G.  _Gradient flows in metric spaces and in the space of probability measures_.  Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition, 2008.
  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S., and Bottou, L.  Wasserstein GAN.  _arXiv preprint arXiv:1701.07875_ , 2017.
  * Benamou and Brenier (2000) Benamou, J.-D. and Brenier, Y.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  _Numer. Math._ , 84(3):375–393, 2000.
  * Bernton et al. (2021) Bernton, E., Ghosal, P., and Nutz, M.  Entropic optimal transport: geometry and large deviations.  _arXiv preprint arXiv:2102.04397_ , 2021.
  * Bobkov and Götze (1999) Bobkov, S. G. and Götze, F.  Exponential integrability and transportation cost related to logarithmic Sobolev inequalities.  _J. Funct. Anal._ , 163(1):1–28, 1999.
  * Brenier (1991) Brenier, Y.  Polar factorization and monotone rearrangement of vector-valued functions.  _Comm. Pure Appl. Math._ , 44(4):375–417, 1991\.
  * Carlier et al. (2016) Carlier, G., Chernozhukov, V., and Galichon, A.  Vector quantile regression: an optimal transport approach.  _The Annals of Statistics_ , 44(3):1165–1192, 2016.
  * Carlier et al. (2017) Carlier, G., Duval, V., Peyré, G., and Schmitzer, B.  Convergence of entropic schemes for optimal transport and gradient flows.  _SIAM Journal on Mathematical Analysis_ , 49(2):1385–1418, 2017.
  * Chernozhukov et al. (2017) Chernozhukov, V., Galichon, A., Hallin, M., and Henry, M.  Monge–kantorovich depth, quantiles, ranks and signs.  _The Annals of Statistics_ , 45(1):223–256, 2017\.
  * Chizat et al. (2020) Chizat, L., Roussillon, P., Léger, F., Vialard, F.-X., and Peyré, G.  Faster wasserstein distance estimation with the sinkhorn divergence.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Conforti and Tamanini (2021) Conforti, G. and Tamanini, L.  A formula for the time derivative of the entropic cost and applications.  _Journal of Functional Analysis_ , 280(11):108964, 2021.
  * Courty et al. (2014) Courty, N., Flamary, R., and Tuia, D.  Domain adaptation with regularized optimal transport.  In _ECML PKDD_ , pages 274–289, 2014.
  * Courty et al. (2017) Courty, N., Flamary, R., Tuia, D., and Rakotomamonjy, A.  Optimal transport for domain adaptation.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 39(9):1853–1865, 2017.
  * Csiszár (1975) Csiszár, I.  $I$-divergence geometry of probability distributions and minimization problems.  _Ann. Probability_ , 3:146–158, 1975.
  * Cuturi (2013) Cuturi, M.  Sinkhorn distances: Lightspeed computation of optimal transport.  In Burges, C. J. C., Bottou, L., Ghahramani, Z., and Weinberger, K. Q., editors, _Advances in Neural Information Processing Systems 26: 27th Annual Conference on Neural Information Processing Systems 2013. Proceedings of a meeting held December 5-8, 2013, Lake Tahoe, Nevada, United States._ , pages 2292–2300, 2013.
  * Deb et al. (2021) Deb, N., Ghosal, P., and Sen, B.  Rates of estimation of optimal transport maps using plug-in estimators via barycentric projections.  _arXiv preprint arXiv:2107.01718_ , 2021.
  * Dvurechensky et al. (2018) Dvurechensky, P., Gasnikov, A., and Kroshnin, A.  Computational optimal transport: Complexity by accelerated gradient descent is better than by sinkhorn’s algorithm.  _arXiv preprint arXiv:1802.04367_ , 2018.
  * Feydy et al. (2017) Feydy, J., Charlier, B., Vialard, F.-X., and Peyré, G.  Optimal transport for diffeomorphic registration.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 291–299. Springer, 2017.
  * Feydy et al. (2019) Feydy, J., Séjourné, T., Vialard, F.-X., Amari, S.-i., Trouvé, A., and Peyré, G.  Interpolating between optimal transport and mmd using sinkhorn divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 2681–2690. PMLR, 2019.
  * Feydy et al. (2020) Feydy, J., Glaunès, J., Charlier, B., and Bronstein, M.  Fast geometric learning with symbolic matrices.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Finlay et al. (2020) Finlay, C., Gerolin, A., Oberman, A. M., and Pooladian, A.-A.  Learning normalizing flows from entropy-kantorovich potentials.  _arXiv preprint arXiv:2006.06033_ , 2020.
  * Flamary et al. (2021) Flamary, R., Courty, N., Gramfort, A., Alaya, M. Z., Boisbunon, A., Chambon, S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., Gautheron, L., Gayraud, N. T., Janati, H., Rakotomamonjy, A., Redko, I., Rolet, A., Schutz, A., Seguy, V., Sutherland, D. J., Tavenard, R., Tong, A., and Vayer, T.  Pot: Python optimal transport.  _Journal of Machine Learning Research_ , 22(78):1–8, 2021.
  * Genevay (2019) Genevay, A.  _Entropy-regularized optimal transport for machine learning_.  PhD thesis, Paris Sciences et Lettres (ComUE), 2019.
  * Genevay et al. (2016) Genevay, A., Cuturi, M., Peyré, G., and Bach, F. R.  Stochastic optimization for large-scale optimal transport.  In Lee, D. D., Sugiyama, M., von Luxburg, U., Guyon, I., and Garnett, R., editors, _Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain_ , pages 3432–3440, 2016.
  * Genevay et al. (2018) Genevay, A., Peyré, G., and Cuturi, M.  Learning generative models with sinkhorn divergences.  In Storkey, A. J. and Pérez-Cruz, F., editors, _International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain_ , volume 84 of _Proceedings of Machine Learning Research_ , pages 1608–1617. PMLR, 2018.
  * Genevay et al. (2019) Genevay, A., Chizat, L., Bach, F., Cuturi, M., and Peyré, G.  Sample complexity of sinkhorn divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1574–1583. PMLR, 2019.
  * Ghosal et al. (2021) Ghosal, P., Nutz, M., and Bernton, E.  Stability of entropic optimal transport and schrödinger bridges.  _arXiv preprint arXiv:2106.03670_ , 2021.
  * Giné and Nickl (2016) Giné, E. and Nickl, R.  _Mathematical foundations of infinite-dimensional statistical models_.  Cambridge Series in Statistical and Probabilistic Mathematics, [40]. Cambridge University Press, New York, 2016.
  * Huang et al. (2021) Huang, C.-W., Chen, R. T. Q., Tsirigotis, C., and Courville, A.  Convex potential flows: Universal probability distributions with optimal transport and convex optimization.  In _International Conference on Learning Representations_ , 2021.
  * Hütter and Rigollet (2021) Hütter, J.-C. and Rigollet, P.  Minimax estimation of smooth optimal transport maps.  _The Annals of Statistics_ , 49(2):1166–1194, 2021.
  * Kantorovitch (1942) Kantorovitch, L.  On the translocation of masses.  _C. R. (Doklady) Acad. Sci. URSS (N.S.)_ , 37:199–201, 1942\.
  * Léonard (2012) Léonard, C.  From the Schrödinger problem to the Monge-Kantorovich problem.  _J. Funct. Anal._ , 262(4):1879–1920, 2012.
  * Makkuva et al. (2020) Makkuva, A., Taghvaei, A., Oh, S., and Lee, J.  Optimal transport mapping via input convex neural networks.  In _International Conference on Machine Learning_ , pages 6672–6681. PMLR, 2020.
  * Manole et al. (2021) Manole, T., Balakrishnan, S., Niles-Weed, J., and Wasserman, L.  Plugin estimation of smooth optimal transport maps.  _arXiv preprint arXiv:2107.12364_ , 2021.
  * Mena and Niles-Weed (2019) Mena, G. and Niles-Weed, J.  Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem.  In Wallach, H. M., Larochelle, H., Beygelzimer, A., d’Alché-Buc, F., Fox, E. B., and Garnett, R., editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 4543–4553, 2019.
  * Monge (1781) Monge, G.  Mémoire sur la théorie des déblais et des remblais.  _Histoire de l’Académie royale des sciences_ , 1:666–704, 1781.
  * Nutz and Wiesel (2021) Nutz, M. and Wiesel, J.  Entropic optimal transport: Convergence of potentials.  _arXiv preprint arXiv:2104.11720_ , 2021.
  * Onken et al. (2021) Onken, D., Fung, S. W., Li, X., and Ruthotto, L.  Ot-flow: Fast and accurate continuous normalizing flows via optimal transport.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 9223–9232, 2021.
  * Pal (2019) Pal, S.  On the difference between entropic cost and the optimal transport cost.  _arXiv preprint arXiv:1905.12206_ , 2019.
  * Peyré and Cuturi (2019) Peyré, G. and Cuturi, M.  Computational optimal transport.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Santambrogio (2015) Santambrogio, F.  _Optimal transport for applied mathematicians_ , volume 87 of _Progress in Nonlinear Differential Equations and their Applications_.  Birkhäuser/Springer, Cham, 2015.  Calculus of variations, PDEs, and modeling.
  * Schiebinger et al. (2019) Schiebinger, G., Shu, J., Tabaka, M., Cleary, B., Subramanian, V., Solomon, A., Gould, J., Liu, S., Lin, S., Berube, P., et al.  Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming.  _Cell_ , 176(4):928–943, 2019.
  * Seguy et al. (2018) Seguy, V., Damodaran, B. B., Flamary, R., Courty, N., Rolet, A., and Blondel, M.  Large-scale optimal transport and mapping estimation.  In _International Conference on Learning Representations_ , 2018.
  * Sinkhorn (1967) Sinkhorn, R.  Diagonal equivalence to matrices with prescribed row and column sums.  _The American Mathematical Monthly_ , 74(4):402–405, 1967.
  * Solomon et al. (2015) Solomon, J., De Goes, F., Peyré, G., Cuturi, M., Butscher, A., Nguyen, A., Du, T., and Guibas, L.  Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains.  _ACM Transactions on Graphics (TOG)_ , 34(4):66, 2015.
  * Solomon et al. (2016) Solomon, J., Peyré, G., Kim, V. G., and Sra, S.  Entropic metric alignment for correspondence problems.  _ACM Trans. Graph._ , 35(4):72:1–72:13, 2016\.
  * Tsybakov (2009) Tsybakov, A. B.  _Introduction to nonparametric estimation_.  Springer Series in Statistics. Springer, New York, 2009.  Revised and extended from the 2004 French original, Translated by Vladimir Zaiats.
  * van der Vaart and Wellner (1996) van der Vaart, A. W. and Wellner, J. A.  _Weak convergence and empirical processes_.  Springer Series in Statistics. Springer-Verlag, New York, 1996.  With applications to statistics.
  * Vershynin (2018) Vershynin, R.  _High-dimensional probability_ , volume 47 of _Cambridge Series in Statistical and Probabilistic Mathematics_.  Cambridge University Press, Cambridge, 2018.  An introduction with applications in data science, With a foreword by Sara van de Geer.
  * Villani (2008) Villani, C.  _Optimal transport: old and new_ , volume 338.  Springer Science & Business Media, 2008.
  * von Luxburg and Bousquet (2003/04) von Luxburg, U. and Bousquet, O.  Distance-based classification with Lipschitz functions.  _J. Mach. Learn. Res._ , 5:669–695, 2003/04.
  * Wang et al. (2010) Wang, W., Ozolek, J. A., Slepčev, D., Lee, A. B., Chen, C., and Rohde, G. K.  An optimal transportation approach for nuclear structure-based pathology.  _IEEE transactions on medical imaging_ , 30(3):621–631, 2010.
  * Yang et al. (2020) Yang, K. D., Damodaran, K., Venkatachalapathy, S., Soylemezoglu, A. C., Shivashankar, G., and Uhler, C.  Predicting cell lineages using autoencoders and optimal transport.  _PLoS computational biology_ , 16(4):e1007828, 2020.