Annealed quantitative estimates for the quadratic 2D-discrete random matching problem

TL;DR

This paper analyzes the quadratic 2D-discrete random matching problem on Riemannian manifolds, providing quantitative approximations of optimal transport plans for correlated and Markov chain samples.

Contribution

It introduces a novel approximation method for optimal transport plans using PDE linearization, applicable to correlated samples and Markov chains.

Findings

Optimal transport plans are well-approximated by PDE-based maps.

The method applies to correlated random points with exponential decay of mixing coefficients.

Results extend to discrete-time Markov chains with invariant measures.

Abstract

We study a random matching problem on closed compact $2$-dimensional Riemannian manifolds (with respect to the squared Riemannian distance), with samples of random points whose common law is absolutely continuous with respect to the volume measure with strictly positive and bounded density. We show that given two sequences of numbers $n$ and $m=m(n)$ of points, asymptotically equivalent as $n$ goes to infinity, the optimal transport plan between the two empirical measures $\mu^n$ and $\nu^{m}$ is quantitatively well-approximated by $\big(\mathrm{Id},\exp(\nabla h^{n})\big)_\#\mu^n$ where $h^{n}$ solves a linear elliptic PDE obtained by a regularized first-order linearization of the Monge-Amp\`ere equation. This is obtained in the case of samples of correlated random points for which a stretched exponential decay of the $\alpha$-mixing coefficient holds and for a class of discrete-time Markov chains having a unique absolutely continuous invariant measure with respect to the volume measure.

Equations

$$\underset{\pi \in {\Pi }_{nm}}{\min }\sum _{i=1}^n\sum _{j=1}^m{\pi }_{ij}{\mathrm{d}}^2(X_i,Y_j)$$

$$\Pi_{nm}:=\bigg{\{}\pi\in[0,1]^{n\times m}\ \Big{|}\ \sum_{i=1}^{n}\pi_{ij}=\tfrac{1}{m}\quad\text{and}\quad\sum_{j=1}^{m}\pi_{ij}=\tfrac{1}{n}\bigg{\}}.$$

$${\mu }^n:=\frac{1}{n}\sum _{i=1}^n{\delta }_{X_i}\text{and}{\nu }^m:=\frac{1}{m}\sum _{j=1}^m{\delta }_{Y_j},$$

$$\underset{\sigma \in {\mathcal{S}}_n}{\min }\sum _{i=1}^n{\mathrm{d}}^2(X_i,Y_{\sigma (i)}),$$

$$\mathbb{E}\big{[}W_{2}^{2}(\mu^{n},\nu^{n})\big{]}\sim\frac{\log(n)}{n},$$

$$\lim_{n\rightarrow\infty}\frac{n}{\log(n)}\mathbb{E}\big{[}W_{2}^{2}(\mu^{n},\nu^{n})\big{]}=\frac{|\mathcal{M}|}{2\pi},$$

$$\lambda \le \rho \le \Lambda .$$

$$({\mu }^n\circ T^n)\det (\nabla T^n)={\nu }^n.$$

$${\mu }^n\circ T^n\approx {\mu }^n+\nabla {\mu }^n\cdot \nabla h^n\text{and}\det (\nabla T^n)\approx 1+\Delta h^n.$$

$$-\nabla \cdot \rho \nabla h^n={\mu }^n-{\nu }^n.$$

$$\lim_{n\uparrow\infty}\bigg{|}\mathbb{E}\big{[}W_{2}^{2}(\mu^{n},\nu^{n})\big{]}-\mathbb{E}\bigg{[}\int_{\mathbb{T}^{2}}|\nabla h^{n}|^{2}\bigg{]}\bigg{|}=0\quad\text{and}\quad\lim_{n\uparrow\infty}\int_{\mathbb{T}^{2}}|T^{n}-(\text{Id}+\nabla h^{n})|^{2}=0\text{ a.s.}$$

$${\text{P}}_t\mu :={\int }_{\mathcal{M}}{p}_t(\cdot ,y)\mathrm{d}\mu (y)\text{and}{\text{P}}_{t}f:={\int }_{\mathcal{M}}{p}_t(\cdot ,y)f(y)\mathrm{dm}(y).$$

$$\alpha_{\ell}:=\sup_{k\geq 1}\sup\bigg{\{}\frac{\text{cov}(f,g)}{\|f\|_{\mathrm{L}^{\infty}}\|g\|_{\mathrm{L}^{\infty}}},\quad f\in\mathrm{L}^{\infty}(\sigma\{X_{j},j\leq k\})\quad\text{and}\quad g\in\mathrm{L}^{\infty}(\sigma\{X_{j},j\geq\ell+k\})\bigg{\}},$$

$$\beta_{ij}:=\sup_{|F|\leq 1}\bigg{|}\int_{\mathcal{M}\times\mathcal{M}}F\,\mathrm{d}(\mathbb{P}_{(X_{i},X_{j})}-\mathbb{P}_{X_{i}}\otimes\mathbb{P}_{X_{j}})\bigg{|}\quad\text{for any $i,j\geq 1$}.$$

$$\underset{n\ge 1}{\sup }\frac{1}{n}\sum _{1\le i<j\le n}{\beta }_{ij}<\infty ,$$

$$\alpha_{\ell}\leq a\exp\big{(}-b\ell^{\eta}\big{)}\quad\text{for any $\ell\geq 1$.}$$

$$\delta^{2}\big{(}(x,y),(z,w)\big{)}:=\mathrm{d}^{2}(x,z)+\mathrm{d}^{2}(y,w).$$

$$\mathrm{H}^{\varepsilon}:=\Big{\{}f:\mathcal{M}\rightarrow\mathbb{R}\ \Big{|}\|f\|^{2}_{\mathrm{H}^{\varepsilon}}:=\sum_{k\geq 1}\lambda^{2\varepsilon}_{k}|\hat{f}(k)|^{2}<\infty\Big{\}},$$

$$\dot{\mathrm{H}}^{1}:=\Big{\{}f:\mathcal{M}\rightarrow\mathbb{R}\ \Big{|}\ \int_{\mathcal{M}}|\nabla f|^{2}<\infty\quad\text{with }\int_{\mathcal{M}}f\mathrm{d}\mathrm{m}=0\Big{\}}.$$

$$-\nabla \cdot \rho \nabla h^{n,t}={\mu }^{n,t}-{\nu }^{m,t},$$

$$\gamma^{n,t}:=\big{(}\mathrm{Id},\exp(\nabla h^{n,t})\big{)}_{\#}\mu^{n,t},$$

$$\underset{\pi }{\inf }{W}_2^2(\pi ,{\gamma }^{n,t})\le {\mathcal{C}}_n\frac{\log (n)}{n}\sqrt{\frac{\log \log (n)}{\log (n)}}\text{with supn≥1E[C1Cn]≤1,}$$

$$\underset{\pi }{\inf }{W}_2^2(\pi ,{\gamma }^{n,t})\le {\mathcal{C}}_n\frac{\log (n)}{n}\sqrt{\frac{\log \log (n)}{{\log }^{1-\frac{2}{\eta }}(n)}}\text{with }\underset{n\ge 1}{\sup }\mathbb{E}[\exp (\frac{1}{C}{\mathcal{C}}_n^{\frac{1}{2}}))]\le 2,$$

$$\pi^{n}=\big{(}\text{Id},T^{n}\big{)}_{\#}\mu^{n}.$$

$$-\nabla \cdot \rho \nabla f^{n,t}={\mu }^{n,t}-{\rho }_t,$$

$$\int_{\mathcal{M}}\mathrm{d}^{2}\big{(}T^{n},\mathrm{exp}(\nabla f^{n,t})\big{)}\mathrm{d}\mathrm{m}\leq\mathcal{C}_{n}\frac{\log(n)}{n}\sqrt{\frac{\log\log(n)}{\log(n)}}\quad\text{with \,$\sup_{n\geq 1}\mathbb{E}[\tfrac{1}{C}\mathcal{C}_{n}]\leq 1$.}$$

$$\int_{\mathcal{M}}\mathrm{d}^{2}\big{(}T^{n},\mathrm{exp}(\nabla f^{n,t})\big{)}\mathrm{d}\mathrm{m}\leq\mathcal{C}_{n}\frac{\log(n)}{n}\sqrt{\frac{\log\log(n)}{\log^{1-\frac{2}{\eta}}(n)}}\quad\text{with }\sup_{n\geq 1}\mathbb{E}[\exp(\tfrac{1}{C}\mathcal{C}_{n}^{\frac{1}{2}}))]\leq 2.$$

$$\underset{k\ge 1}{\sup }\Vert {\varphi }_k{\Vert }_{{\mathrm{L}}^{\infty }}<\infty .$$

$$(K\mu )(A)={\int }_{\mathcal{M}}K(\cdot ,A)\mathrm{d}\mu \text{for every Borel set A⊂M and μ∈P(M),}$$

$$(K\psi )(x)={\int }_{\mathcal{M}}\psi (y)K(x,\mathrm{d}y)\text{for any x∈M.}$$

Full text

NC has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 948819)  

FM is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the SPP 2265 Random Geometric Systems. FM has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2044 -390685587, Mathematics Münster: Dynamics–Geometry–Structure. FM has been funded by the Max Planck Institute for Mathematics in the Sciences.

**Keywords: **Optimal transport $\cdot$ Matching problem $\cdot$ Quantitative estimates

1. 1 Introduction and statement of the main results

2. 1.1 The random matching problem and its asymptotic

3. 1.2 Linearization ansatz

4. 1.3 Formulation of the main results

5. 1.4 Extension to a class of Markov chains

6. 1.5 Open problems

7. 2 Structure of the proof

8. 3 Proofs

9. 3.1 Notations and preliminary results

10. 3.2 $\mathrm{L}^{q}$-type estimates

11. 3.3 Fluctuation estimates

12. 3.4 Contractivity estimates

13. 3.5 Proof of Theorem 1.2: Approximation of the transport plan

14. A Probabilistic and PDE tools

15. B Matching cost for point clouds

16. C Proof for the class of Markov chains

1. Introduction and statement of the main results

1.1. The random matching problem and its asymptotic

The random matching problem is a popular optimization problem at the interface between analysis and probability with applications in many different fields such as statistical physics [16, 44], computer science [40] and economics [20, 24]. Within the mathematical literature, it has been subject of intense studies due to its interactions with many areas, including for instance graph theory [38] and geometric probability [50]. In this paper we focus on one of its simple versions. Let $\{X_{k}\}_{1\leq k\leq n}$ and $\{Y_{k}\}_{1\leq k\leq m}$ (with possibly $m>n$) be two families of random points on a compact Riemannian manifold $\mathcal{M}$ (endowed with the Riemannian distance $\mathrm{d}$). We are interested in the quadratic matching problem

[TABLE]

where

[TABLE]

Classically, (1.1) can be phrased in terms of a transport problem. Indeed, letting

[TABLE]

be the empirical measures associated to the two point clouds, the linear programming problem (1.1) amounts to determine the quadratic Wasserstein distance $\ W_{2}^{2}(\mu^{n},\nu^{m})$.

In the special case $n=m$ the Birkhoff-von Neumann’s Theorem provides a correspondence between (1.1) and the usual bipartite matching

[TABLE]

where $\mathcal{S}_{n}$ denotes the set of injective maps $\sigma:\{1,\dots,n\}\rightarrow\{1,\dots,n\}$. Indeed, since $\Pi_{nn}$ is a convex polytope, minimizers in (1.1) have to be searched among extremal points. By the Birkhoff-von Neumann Theorem [8, Lemma 2.1.3], the latter are nothing but permutation matrices (up to a factor $\tfrac{1}{n}$).

A first natural question is to understand the asymptotic of (1.1) as $n,m\uparrow\infty$. For the same number of samples $n=m$ and independently and identically distributed (i.i.d.) on the unit square $[0,1]^{d}$, the scaling of the cost (1.1) has been well understood in the mathematical and statistical physics literature. A simple heuristic argument, see for instance [44], suggests that given a point $X_{i}$, we can find a point $Y_{j}$ within a volume of order $O(n^{-1})$ with high probability. For this reason, the typical inter-point distance is of order $O(n^{-\frac{1}{d}})$ suggesting that the scaling of (1.1) is of order $O(n^{-\frac{2}{d}})$. Although attractive, this heuristic turns out to be unfortunately false in low dimension showing a critical behavior when $d=2$. This critical case is the one on which we focus in this paper. Ajtai, Komlós and Tusnády in [1] were the first to show that, for i.i.d. uniform samples, a logarithmic correction is needed, deriving111We use the notation $A\lesssim B$ if there exists a global constant $C>0$, which may only depend on $d$, such that $A\leq CB$. We write $A\sim B$ if both $A\lesssim B$ and $B\lesssim A$ hold.

[TABLE]

extended later by Talagrand in [53] for clouds of i.i.d. points which are distributed accordingly to more general common law. A recent breakthrough was obtained within the physics community by Caracciolo, Lucibello, Parisi and Sicuro in [16], and further developed by Caracciolo and Sicuro in [17] and by Sicuro in [49], where the asymptotics of the cost are formally derived thanks to a novel PDE approach and optimal transport theory rather than combinatorics. A couple of years later, in general $2$-dimensional compact Riemannian manifolds without boundary, the first-order asymptotic has been rigorously justified by Ambrosio, Stra and Trevisan in [6] for i.i.d. uniform samples and recently extended by Ambrosio, Goldman and Trevisan in [5] for samples distributed accordingly to more general laws which are absolutely continuous (with Hölder continuous density) w.r.t. the volume measure $\mathrm{d}\mathrm{m}$, leading to

[TABLE]

where $|\mathcal{M}|$ denotes the Lebesgue measure of $\mathcal{M}$. The case $n\neq m$ with $n,m\uparrow\infty$ with similar rates is also covered, see [5, Theorem 1.2].

The novel approach introduced in [16], later revised in [11] consists on a linearization of the Monge-Amperé equation that allows for an explicit description of the cost thanks to the linearized proxies (see Section 1.2 for more details). The aim of this work is to quantitatively justify the linearization ansatz in terms of convergence of the approximating minimizers of (1.1) towards the optimal ones. In particular, we are interested in the case where the points are identically distributed with a common law $\rho\,\mathrm{d}\mathrm{m}$ (where we recall that $\mathrm{d}\mathrm{m}$ denotes the volume measure) and $\rho$ satisfies for some $\lambda,\Lambda>0$

[TABLE]

To the best of our knowledge, there are only few results on the asymptotic behavior of the transport map and they are so far limited to the case of i.i.d. uniform samples in the study of the semi-discrete matching problem (that is couplings between $\mu^{n}$ and $\mathrm{d}\mathrm{m}$), see the work of Ambrosio, Glaudo and Trevisan in [4]. In connection with this work, quantitative estimates on the optimal map for the matching between the Lebesgue measure and Poisson clouds have been obtained by Goldman, Huesmann and Otto in [29] and Goldman and Huesmann in [28].

Our extension in this paper is fourfold: First, we look at more general distribution of points and we consider the case of general densities $\rho$ satisfying (1.6). Second, we do not assume independence and we consider samples which may possess correlations. Third, we do not restrict the analysis to the semi-discrete matching problem and we also investigate the ansatz for the full matching problem (1.1). Finally, we investigate the case where the points are not identically distributed and we extend our result to a class of Harris positive recurrent Markov chains.

We finally mention that the effectiveness of the linearization ansatz introduced in [16] is not only limited to the case of i.i.d. distributed points on bounded domains, but it can be employed in many different settings. See for instance [15] for an interesting application to the matching on unbounded domains, [35, 36, 37] for an application to Gaussian matching, [33] for an application in stochastic matrix theory, [58, 59, 60, 32] for an application to a continuous instance of the matching problem, i.e. when the empirical measure is replaced by the occupation measure of a stochastic process. It is worth to further mention that these techniques can also be employed when considering the matching problem with $p$-costs in higher dimension, see [30].

1.2. Linearization ansatz

We now briefly reproduce the linearization ansatz introduced in [16]. For simplicity, we consider the case $\mathcal{M}=\mathbb{T}^{2}$, $n=m$ and i.i.d. samples with common distribution $\rho\,\mathrm{d}\mathrm{m}$. Let $T^{n}$ be an optimal transport map (whose existence is ensured by Brenier’s Theorem [14]) between $\mu^{n}$ and $\nu^{n}$. Based on the transport relation $T^{n}_{\#}\mu^{n}=\nu^{n}$ and a change of variables, $T^{n}$ solves (formally) the Monge-Ampère equation

[TABLE]

Since the cost is quadratic, by [48, Theorem 1.25], there exists a function $h^{n}$ such that $T^{n}=\mathrm{Id}+\nabla h^{n}$. Applying the law of large numbers, we further have the weak convergences $\mu^{n},\nu^{n}{\rightharpoonup}\rho\,\mathrm{d}\mathrm{m}$ as $n\uparrow\infty$ so that we expect $T^{n}\approx\mathrm{Id}$ as $n\uparrow\infty$. Thus, this suggests that the correction $\nabla h^{n}$ is small as $n\uparrow\infty$, allowing to perform (formally) the Taylor expansions

[TABLE]

Plugging (1.8) into (1.7), neglecting the higher order terms and replacing $\mu^{n}$ by $\rho$ yields

[TABLE]

This formal linearization suggests the following two conjectures

[TABLE]

Unfortunately, (1.10) cannot hold as it is, since the solution of (1.9) does not belong to $\mathrm{H}^{1}$ due to the roughness of the source term. To overcome this, following the strategy in [6], a regularization using the heat-semigroup at time $\sim\frac{1}{n}$ is made. Doing so, the first item of (1.10) turns out to be true, leading to the result (1.5) (see for instance [3] for a convergence rate).

1.3. Formulation of the main results

For the remainder of the paper $\mathcal{M}$ denotes a $2$-dimensional connected and compact Riemannian manifold without boundary (or the square $[0,1]^{2}$) endowed with the Riemannian distance $\mathrm{d}$. For $t>0$ we denote by $p_{t}$ the fundamental solution of the heat operator $\partial_{t}-\Delta$ on $\mathcal{M}$, where $\Delta$ denotes the Beltrami-Laplace operator. We define the heat semigroup $(\text{P}_{t})_{t>0}$ via its action on probability measures $\mu\in\mathcal{P}(\mathcal{M})$ and square integrable functions $f\in\mathrm{L}^{2}(\mathcal{M})$

[TABLE]

We first introduce the class of correlated point clouds that we consider for studying the matching problem (1.1). This class concerns point clouds $\{X_{i}\}_{i}$ for which the correlations between points decay at an exponential rate, where the correlations are measured in terms of the $\alpha$-mixing coefficients given by, for any $\ell\geq 1$

[TABLE]

and the $\beta$-mixing coefficients given by222We denote by $\mathbb{P}_{X}$ the law of a random variable $X$.

[TABLE]

Assumption 1.1** (Correlated point clouds).**

We consider point clouds $\{X_{i}\}_{i}\subset\mathcal{M}$ which are identically distributed according to $\rho\,\mathrm{d}\mathrm{m}$ where $\rho$ satisfies (1.6). We further assume decay of the correlations in the form of

[TABLE]

and there exist $a,b>0$ and $\eta\in(0,\infty]$ such that

[TABLE]

Our first main result concerns the approximation of transport plans coupling $\{\mu^{n}\}_{n}$ and $\{\nu^{m}\}_{m}$ defined in (1.2). We justify the formal linearization of the Monge-Ampère equation achieved in Section 1.2 in a annealed quantitative way (i.e. in expectation): We show that, a suitable regularization of, the plan \big{(}\text{Id},\exp(\nabla h^{n})\big{)}_{\#}\mu^{n}, with $h^{n}$ defined in (1.9), provides a good approximation when measuring the error with respect to the $W_{2}$-Wasserstein distance in the product space $\mathcal{M}\times\mathcal{M}$ endowed with the metric

[TABLE]

The density $\rho$ will need further regularity in form of fractional Sobolev spaces defined as, for some $\varepsilon>0$

[TABLE]

where $\{\lambda_{k},\phi_{k}\}_{k}$ denote the eigenvalues and eigenvectors of $-\Delta$ on $\mathcal{M}$ and $\hat{f}(k)=\int_{\mathcal{M}}f\,\phi_{k}$ denotes the Fourier modes of $f$. Finally, we denote by $\dot{\mathrm{H}}^{1}$ the $\mathrm{L}^{2}$-based Sobolev space

[TABLE]

Theorem 1.2** (Approximation of the transport plan).**

Let $\rho\in\mathrm{H}^{\varepsilon}$ for some $\varepsilon>0$ satisfying (1.6) and $\{\mu^{n}\}_{n}$ and $\{\nu^{m}\}_{m}$ be defined in (1.2) (for $m=m(n)$ with some given increasing map $m:\mathbb{N}\rightarrow\mathbb{N}$) with point clouds satisfying Assumption 1.1 and such that there exists $q\in[1,\infty)$ for which $\frac{m(n)}{n}\underset{n\uparrow\infty}{\rightarrow}q$. We consider333where we impose additional Neumann boundary conditions in the case $\mathcal{M}=[0,1]^{2}$ $h^{n,t}\in\dot{\mathrm{H}}^{1}$ the weak solution of

[TABLE]

for any $t\in(0,1)$ with $\mu^{n,t}:=\mathrm{P}_{t}\mu^{n}$ and $\nu^{m,t}:=\mathrm{P}_{t}\nu^{m}$.

There exist an exponent $\kappa>0$, a deterministic constant $C$ and a random variable $\mathcal{C}_{n}$ both depending on $\lambda,\Lambda$ and $\mathcal{M}$ for which given $t=\frac{\log^{\kappa}(n)}{n}$ and

[TABLE]

it holds

[TABLE]

where the $\inf$ runs over all optimal transport plans $\pi$ between $\mu^{n}$ and $\nu^{m}$.

Furthermore, if (1.14) holds with $\eta>2$, the assumption (1.13) can be dropped and it holds

[TABLE]

where the $\inf$ runs over all optimal transport plans $\pi$ between $\mu^{n}$ and $\nu^{m}$.

Our second main result concerns the particular case of the semi-discrete matching problem, i.e. optimal coupling between the common law $\rho\,\mathrm{d}\mathrm{m}$ and $\{\mu^{n}\}_{n}$. We know from McCann’s theorem [39] that there exists a unique optimal transport map $T^{n}$, that is the optimal transport plan $\pi^{n}$ can be written as

[TABLE]

We show that $T^{n}$ can be approximated in a annealed quantitative way in $\mathrm{L}^{2}$ by (a suitable regularized version of) the solution of (1.9) with $\nu^{m}$ replaced by $\rho$.

Theorem 1.3**.**

Let $\rho\in\mathrm{H}^{\varepsilon}$ for some $\varepsilon>0$ satisfying (1.6) and $\{\mu^{n}\}_{n}$ be defined in (1.2) with a point cloud satisfying Assumption 1.1. We consider $f^{n,t}\in{\rm\dot{H}^{1}}$ the weak solution of

[TABLE]

where, for all $t\in(0,1)$, we recall that $\mu^{n,t}=\mathrm{P}_{t}\mu^{n}$ and $\rho_{t}=\mathrm{P}_{t}\rho$. Finally, we denote by $T^{n}$ the optimal transport map from $\rho\,\mathrm{d}\mathrm{m}$ to $\mu^{n}$.

There exist an exponent $\kappa>0$, a deterministic constant $C$ and a random variable $\mathcal{C}_{n}$ both depending on $\lambda,\Lambda$ and $\mathcal{M}$ for which given $t=\frac{\log^{\kappa}(n)}{n}$, it holds

[TABLE]

Furthermore, if (1.14) holds with $\eta>2$, the assumption (1.13) can be dropped and it holds

[TABLE]

We finally mention that in the case where the eigenfunctions $\{\phi_{k}\}_{k}$ admit a uniform bound, the conclusions (1.21) and (1.24) can be improved. We comment on the proof at the end of Sub-section 3.5.

Remark 1.4**.**

Let $\{\mu^{n}\}_{n}$ and $\{\nu^{m}\}_{m}$ be as in Theorem 1.2. We assume that the family of eigenfunctions $\{\phi_{k}\}_{k}$ satisfies the uniform bound

[TABLE]

Then, (1.21) and (1.24) hold true for $\eta>1$ with a convergence rate $\frac{\log(n)}{n}\sqrt{\frac{\log\log(n)}{\log^{1-\frac{1}{\eta}}(n)}}$ with the same stochastic integrability. Note that (1.25) typically holds when the geometry of $\mathcal{M}$ is flat, see [54].

Theorem 1.2 and Theorem 1.3 are not restricted to the case of identically distributed point clouds and we present in the next section and in Appendix C a possible extension, using the same techniques, to a class of regular discrete-time Markov chains.

1.4. Extension to a class of Markov chains

We first recall some basic facts on discrete-time Markov chains on $\mathcal{M}$. Such a Markov process is described by its initial distribution $\mu_{0}\in\mathcal{P}(\mathcal{M})$ and its transition kernel $K$, that is a measurable map from $\mathcal{M}$ to the space of probability measures $\mathcal{P}(\mathcal{M})$. We recall that $K$ acts on $\mathcal{P}(\mathcal{M})$ in form of

[TABLE]

and likewise on bounded measurable function $\psi$ in form of

[TABLE]

Given an initial distribution $\mu_{0}\in\mathcal{P}(\mathcal{M})$, we recall that the law of a Markov chain $\{X_{n}\}_{n\geq 0}$ can be expressed with the help of the transition kernel, namely

[TABLE]

where $K^{n}$ stands for the $n^{\mathrm{th}}$-iteration of the kernel $K$.

We now introduce the class of discrete-time Markov chains we consider.

Assumption 1.5**.**

Let $\mathcal{M}$ be a $2$-dimensional compact Riemannian manifold. Let $\mu_{0}\in\mathcal{P}(\mathcal{M})$ and $\{X_{n}\}_{n\geq 1}\subset\mathcal{M}$ be a Markov chain with initial law $\mu_{0}$. We assume that the chain satisfies the two following conditions:

• (i)

We assume that there exists a measurable function $k:\mathcal{M}\times\mathcal{M}\rightarrow[0,1]$ and $\lambda,\Lambda>0$ such that for any Borel set $A\subset\mathcal{M}$

[TABLE]

• (ii)

We assume that there exists an unique invariant measure $\mu_{\infty}$, i.e. $K\mu_{\infty}=\mu_{\infty}$, and $V:\mathcal{M}\rightarrow[0,\infty)$, such that $\int V\,\mathrm{d}\mu_{0}<\infty$, as well as $\sigma\in(0,1)$ for which

[TABLE]

with

[TABLE]

We now comment on the consequences of the above assumptions. First, the condition (1.29) ensures that the invariant measure $\mu_{\infty}$ is absolutely continuous with respect to the volume measure with bounded density, namely

[TABLE]

We briefly give the argument for (1.31). Using the second item of (1.29), we have that the operator

[TABLE]

Furthermore, we note that the closed convex set $\mathcal{C}:=\{f\in\mathrm{L}^{1}\,|\,\lambda\leq f\leq\Lambda\text{ and }\int_{\mathcal{M}}f=1\}$ is invariant under the action of $\mathcal{I}$. Therefore Schauder’s fixed point theorem implies that $\mathcal{I}$ admits a fixed point in $\mathcal{C}$. Given such a fixed point $\rho$, it is clear that $\mu_{\infty}$ defined in (1.31) is an invariant measure according to (1.29).

Second, we note that by Harris’ theorem [31, Theorem 3.7], the condition (1.30) is satisfied when the kernel satisfies the following geometric drift condition: There exist a function $V:\mathcal{M}\rightarrow[0,\infty)$ and constants $\tau\in(0,1)$ and $C\geq 0$ such that

[TABLE]

and the kernel $K$ is bounded on the sublevel sets of $V$, i.e. for all $R\geq 0$, there exists a constant $0<a<1$ for which

[TABLE]

Moreover, the assumption (1.30) implies the geometric decay of the $\beta$-mixing coefficient (1.12), namely there exists a constant $C$ depending on $\lambda$, $\Lambda$ and $\mathrm{m}(\mathcal{M})$ such that

[TABLE]

which ensures that (1.13) and (1.14) (with $\eta=1$ and $b=-\log(\sigma)$) hold and ensures good concentration property of the Markov chain, see Proposition A.1. The estimates (1.32) can be seen as a direct consequence of the combination of the representation for the $\beta$-mixing in [21, Proposition 1] and the assumption (1.29).

Finally, the condition (1.30) quantifies the weak convergence of the law of the Markov chain to its stationary distribution, namely there exists a constant $C>0$ such that for any $f\in\mathrm{L}^{\infty}(\mathcal{M})$

[TABLE]

We shortly give the argument. We first notice that a direct inductive argument together with the semigroup property $K^{n_{1}+n_{2}}=K^{n_{1}}K^{n_{2}}$ for every $n_{1},n_{2}>0$ and Fubini’s theorem gives

[TABLE]

The combination of (1.28), (1.34) and (1.30) gives

[TABLE]

which leads to (1.33) using that $\|f-\mu_{\infty}(f)\|_{V}\leq\|f\|_{\mathrm{L}^{\infty}}$ and by Assumption 1.5 $\int_{\mathcal{M}}V(x)\,\mathrm{d}\mu_{0}<\infty$.

A classical example of a Markov chain satisfying Assumption 1.5 is given by iterated function systems with additive noise. For simplicity, let $\mathcal{M}=\mathbb{T}^{2}$. Let $\{\theta_{n}\}_{n\geq 1}$ be i.i.d. random variables with common law $h\,\mathrm{d}\mathrm{m}$ for some $h$ satisfying $\lambda\leq h\leq\Lambda$. Let $F:\mathbb{T}^{2}\rightarrow\mathbb{T}^{2}$ be a contraction, i.e. there exists a constant $L<1$ such that

[TABLE]

We define the iterated function system $\{X_{n}\}_{n\geq 1}$ according to the induction

[TABLE]

The kernel is given by

[TABLE]

so that $K$ satisfies (1.29) with

[TABLE]

Moreover, the condition (1.35) ensures the validity of (1.30), see for instance [2, Theorem 3.2]. Thus the Markov process $\{X_{n}\}_{n\geq 1}$ satisfies Assumption 1.5.

We show in Appendix C that the conclusions of Theorem 1.3 and Theorem 1.2 hold true for Markov chains satisfying Assumption 1.5.

1.5. Open problems

We conclude this section with open questions that arise in view of our results:

1. (1)

*Sharpness of the rates. *The convergence rate in (1.20) and (1.23) match with the one obtained for the case of uniformly distributed samples in [4]. However, even in the latter case, it has not been shown whether the rate is optimal and we suspect the opposite. A possible way to track the optimal rate could be to perform a second-order linearization of the Monge-Ampère equation (1.7). Following the same type of computations leading to (1.9) in the case $\rho=1$, a second-order linearization $q^{n}$ should solve

[TABLE]

where we recall that $h^{n}$ solves (1.9), providing the conjecture

[TABLE]

This direction of research is left for future investigations. 2. (2)

*Extension to $p$-costs. *A natural question is to investigate if our results hold for different cost functions as $p$-cost functions for $p>1$. The behavior of the cost has been optimally quantified in [1, 13, 9, 23]. However, to the best of our knowledge, quantitative estimates on the transport plan in the setting of general $p$-costs are not known, even for uniformly distributed samples. A possible approach would be to revise the linearization ansatz of [16] for general $p$-costs. As observed in [34], if the transport cost between two points $x,y$ is given by $\frac{1}{p}|x-y|^{p}$, Gangbo-McCann’s theorem [25, Theorem 1.2] ensures that there exists a map $h^{n}$ such that the optimal transport map $T^{n}$ takes the form $T^{n}=\text{Id}+|\nabla h^{n}|^{p^{\prime}-2}\nabla h^{n}$, where $p^{\prime}$ denotes the conjugate exponent. Therefore, following the same type of computations leading to (1.9), a first-order linearization should solve the following degenerate $p^{\prime}$-Laplace equation

[TABLE]

and we may expect

[TABLE]

2. Structure of the proof

This section is devoted to describe the main ideas and how are organized the proofs of Theorem 1.2 and 1.3. We mainly focus on the proof of Theorem 1.2 since the proof of Theorem 1.3 follows by the same strategy.

General strategy.

The proof of Theorem 1.2 follows the strategy employed in [4] to deal with independent and uniformly distributed random points. The main idea is to use the quantitative stability result for transport maps in [4, Theorem 3.2], stating that two transports map are close in the $\mathrm{L}^{2}$-topology if the target measures are close in the $W_{2}$-topology. We restate below the result for reader’s convenience.

Theorem 2.1** (Stability of transport maps).**

Let $\nu,\mu_{1},\mu_{2}\in\mathcal{P}(\mathcal{M})$ such that $\nu\ll\mathrm{m}$ and let $S,T:\mathcal{M}\rightarrow\mathcal{M}$ be the optimal transport maps respectively for the pairs of measures $(\nu,\mu_{1})$ and $(\nu,\mu_{2})$. We assume that $S=\exp(\nabla f)$ for some $f:\mathcal{M}\rightarrow\mathbb{R}$ with $\mathrm{C}^{1,1}$-regularity.

There exists a constant $c>0$ depending on $\mathcal{M}$ such that, provided

[TABLE]

we have

[TABLE]

The first step consists of using Theorem 2.1 to deduce a stability estimate (in terms of the quadratic Wasserstein distance) of transport plans in the special case where $\mu_{1}=\nu^{m}$, $\mu_{2}=\exp(\nabla h^{n,t})_{\#}\mu^{n}$ and $\nu=\mu^{n}$. In this step, we immediately face the issue of the lack of regularity of $h^{n,t}$ necessary to ensure the condition (2.1): Indeed, recalling that $h^{n,t}$ solves (1.18), it does not have the $\mathrm{C}^{1,1}$-regularity condition for non-smooth densities $\rho$ that we consider here. We overcome this issue introducing an additional regularization step: We smooth the operator $-\nabla\cdot\rho\nabla$ and implicitly $\gamma^{n,t}$ in form of

[TABLE]

for a regularization parameter $\delta$ to be optimized and $\rho_{\delta}:=\text{P}_{\delta}\rho$. Classical Schauder’s theory ensures that $h^{n,t}_{\delta}$ owns $\text{C}^{\infty}$-regularity. Doing so, we can use Theorem 2.1 to deduce, provided that

[TABLE]

a stability error estimate which reads

[TABLE]

where we refer to (3.107) for further details. Our argument then differs from [4] by splitting the error in two parts

[TABLE]

Our proof then proceeds in two steps, controlling separately the two terms in (2.5).

Control of the regularization error

To deal with the regularization error, we look at the difference $e^{n,t}_{\delta}:=h^{n,t}_{\delta}-h^{n,t}$ which solves, according to (1.18) and (2.2),

[TABLE]

Using an energy estimate, we get

[TABLE]

On the one hand, since $\rho\in\mathrm{L}^{\infty}$, we have $\rho_{\delta}\rightarrow\rho$ as $\delta\downarrow 0$ in every $\mathrm{L}^{q}$ with $q<\infty$. On the other hand, we learn from Meyers’ estimate (recalled in Proposition A.3) that there exists $\bar{q}>2$ such that $\nabla h^{n,t}\in\mathrm{L}^{\bar{q}}$. Consequently, we can treat (2.7) using Hölder’s inequality which provides

[TABLE]

The next step is to control the averaged Meyers’ norm \mathbb{E}\big{[}\|\nabla h^{n,t}\|^{2}_{\mathrm{L}^{\bar{q}}}\big{]}, that we show in Proposition 3.3 to be of order of

[TABLE]

where we recall that $\eta$ denotes the correlation length, see Assumption 1.1. The combination of (2.8), (2.9) and local Lipschitzianity of the exponential map yields

[TABLE]

where we refer to (3.104) for further details.

Control of the stability error

For the stability error, we first need to ensure (2.3). Our strategy follows the idea in [3] which consists of showing that (2.3) is satisfied with very high probability. In our case, a new difficulty comes from our regularization of $\rho$ and the regularization parameter $\delta$ has to be carefully optimized. We show that if $\delta$ is taken as an inverse power of $\log(n)$, (2.3) becomes very likely as $n\uparrow\infty$. More precisely, we show in Proposition 3.4 that given two exponents $\kappa_{1}$ and $\upsilon\gg\kappa_{1}$, there exists $\kappa_{2}$ such that given the choices

[TABLE]

we have

[TABLE]

This result is in the spirit of [3, Theorem 3.3] but, in our setting, the proof relies on Schauder’s theory rather than an explicit formula for $h^{n,t}_{\delta}$ as well as concentration inequalities in form of Proposition A.1 to treat the correlations. For further details on the strategy, we refer to Section 3.3. With (2.12) in hands, we can restrict the analysis to realizations satisfying $\|\nabla h^{n,t}_{\delta}\|_{\mathrm{L}^{\infty}}+\|\nabla^{2}h^{n,t}_{\delta}\|_{\mathrm{L}^{\infty}}\leq\tfrac{1}{\log^{\upsilon}(n)}$ where, for $n\gg 1$, (2.3) is satisfied which puts us in the validity of (2.4). We then treat each terms appearing in (2.4) separately:

• •

The optimal control of the cost $W^{2}_{2}(\mu^{n},\nu^{m})$ has been already well studied and is optimally estimated by

[TABLE]

We refer to Appendix B for a detailed statement, references and extensions to the cases of Assumption 1.1 and Assumption 1.5.

• •

The smoothing errors $W^{2}_{2}(\mu^{n},\mu^{n,t})$ and $W^{2}_{2}(\nu^{m},\nu^{m,t})$. Classical contractivity estimates are known and are usually applied to deal with these errors, see for instance [22, Theorem 3], which bound the errors by $t$. However, due to the choice of $t$ in (2.11), this result is of no use since $t$ is much larger than the magnitude of the cost, namely $t\gg\frac{\log(n)}{n}$. Instead, we follow the approach in [3], where the authors showed that in the particular case of empirical measures in dimension $2$, we can improve the rate and obtain the bound $\frac{\log\log(n)}{n}\ll\frac{\log(n)}{n}$. We extend this result to our setting of non-constant densities and correlated points. In Proposition 3.5, we derive

[TABLE]

As opposed to [3], our approach uses Fourier analysis together with additional cares to handle the correlations and non-constant densities.

• •

The error in the Moser coupling W^{2}_{2}\big{(}\nu^{m,t},\exp(\nabla h^{n,t}_{\delta})_{\#}\mu^{n,t}\big{)}. We follow the strategy in [3]. This error can be related to a Moser coupling between $\mu^{n,t}$ and $\nu^{m,t}$ (see for instance [56, Appendix p. 16]): The equation (1.18) gives a natural coupling between $\mu^{n,t}$ and $\nu^{m,t}$ which can be formulated using Benamou-Brenier’s theorem [10],

[TABLE]

Then, using the transport plan \big{(}\phi(1,\cdot),\exp(\nabla h^{n,t}_{\delta})\big{)}_{\#}\mu^{n,t} as a competitor, we get

[TABLE]

that we combine with a quantitative stability result for flows of vector fields, proved in [3, Proposition A.1], leading to

[TABLE]

where we recall that $\bar{q}$ denotes the Meyers’ exponent introduced in (2.7). For further details, we refer to (3.115). It then remains to appeal to Meyers’ estimate, see Proposition A.3, to (2.6) together with (2.9) in form of

[TABLE]

which finally yields

[TABLE]

To conclude, we see that all bounds involve errors in terms of the approximation of $\rho$ using the heat-semigroup, that we need to quantify. This is where the assumption $\rho\in\mathrm{H}^{\varepsilon}$ plays a role, in form of the quantitative estimate

[TABLE]

see (3.94) for a proof. Combining (2.4), (2.5), (2.10), (2.13), (2.14), (2.15), (2.16) with the choices of $\delta$ and $t$ in (2.11), we obtain Theorem 1.2. The proof of Theorem 1.3 is obtained using the same strategy where the first step is simpler, since we apply directly Theorem 2.1 with $\mu_{1}=\rho$, $\nu=\mu^{n}$ and $\mu_{2}=\exp(\nabla f^{n,t}_{\delta})_{\#}\mu^{n}$ where $f^{n,t}$ solves $-\nabla\cdot\rho_{\delta}\nabla f^{n,t}=\mu^{n,t}-\rho_{t}$.

3. Proofs

3.1. Notations and preliminary results

We provide in this section some notations and preliminary results needed in the proofs of Theorem 1.2 and Theorem 1.3. We recall that throughout the paper, we denote by $\mathcal{M}$ a $2$-dimensional compact connected Riemannian manifold (or the square $[0,1]^{2}$) endowed with the Riemannian distance $\mathrm{d}$.

**Wasserstein distance. **Given $\mu,\nu\in\mathcal{P}(\mathcal{M})$, we define the quadratic Wasserstein distance as

[TABLE]

where $\Pi(\mu,\nu)$ is the set of couplings between $\mu$ and $\nu$, that is the set of $\pi\in\mathcal{P}(\mathcal{M}\times\mathcal{M})$ having $\mu$ and $\nu$ as first and second marginal, respectively. We refer the reader to the monographs [55, 56] for a detailed overview on the subject of optimal transport. We recall the following simple, but useful Lipschitz contraction property of the Wasserstein distance.

Lemma 3.1** (Lipschitz property of the Wasserstein metric).**

Let $(D,\mathrm{d}_{D})$ be a complete and separable metric space, let $\mu,\nu\in\mathcal{P}(\mathcal{M})$ and let $T:\mathcal{M}\rightarrow D$ be a $L$-Lipschitz map. It holds

[TABLE]

Proof.

For any coupling $\pi\in\Pi(\mu,\nu)$, the push-forward $(T,T)_{\#}\pi$ is a coupling between $T_{\#}\mu$ and $T_{\#}\nu$. Moreover, it holds

[TABLE]

Taking the infimum among all possible couplings $\pi\in\Pi(\mu,\nu)$ leads to (3.2). ∎

**Heat semigroup and heat kernel. **We recall some basic facts on the heat semigroup and its generator, we refer the reader to [18, Chapter VI] for a more detailed overview. For $t>0$, we denote by $p_{t}$ the fundamental solution of the heat operator $\partial_{t}-\Delta$ on $\mathcal{M}$, where $\Delta$ denotes the Beltrami-Laplace operator. Classical Schauder’s theory ensures that $p_{t}$ is smooth and it is known, see for instance [52] and [3, Appendix B], that $p_{t}$ and its derivatives satisfy for some $C>0$ depending on $\mathcal{M}$

[TABLE]

The kernel $p_{t}$ admits the spectral decomposition

[TABLE]

converging in $\mathrm{L}^{2}(\mathcal{M}\times\mathcal{M})$, where we recall that $\{\lambda_{k},\phi_{k}\}$ denotes the eigenvalues and eigenvectors of $-\Delta$ on $\mathcal{M}$. Specifying (3.4) on the diagonal and using $\|\phi_{k}\|_{\mathrm{L}^{2}(\mathcal{M})}=1$, we obtain the trace formula

[TABLE]

We recall that $\{\phi_{k}\}_{k\geq 1}$ can be estimated in terms of the eigenvalues,

[TABLE]

We briefly recall the argument. Applying the Gagliardo-Nirenberg’s interpolation inequality [7, Theorem 3.70], it holds

[TABLE]

In combination with $-\Delta\phi_{k}=\lambda_{k}\phi_{k}$ and elliptic regularity, in form of

[TABLE]

we obtain (3.6).

We recall that $(\text{P}_{t})_{t>0}$ admits the spectral gap property, that is there exists a constant $C_{\text{sg}}>0$ such that

[TABLE]

Note that equivalently, (3.7) can be formulated in terms of the eigenvalues $\{\lambda_{k}\}_{k\geq 1}$ in form of

[TABLE]

simply by specifying (3.7) on $\{\phi_{k}\}_{k\geq 1}$. Finally, we recall the Riesz-transform bound

[TABLE]

where $(-\Delta)^{\frac{1}{2}}$ can be defined via its spectral representation

[TABLE]

with $(\cdot,\cdot)_{\mathrm{L}^{2}}$ the inner product in $\mathrm{L}^{2}$. We refer to the monograph [57] for a discussion of the inequalities (3.7) and (3.9), see Chapter 1 for the case of a Riemannian manifold without boundary and Chapter 2 for the case of a Riemannian manifold with (convex) boundary. In connection with the Wasserstein metric, the heat semigroup satisfies the following classical contraction property.

Lemma 3.2** (Semigroup contraction for absolutely continuous measures).**

Let $\rho\in\mathrm{L}^{\infty}$ satisfy (1.6). Given $\rho_{t}:=\mathrm{P}_{t}\rho$, it holds

[TABLE]

Proof.

Using $g$ defined via $-\Delta g=\rho_{t}-\rho$ together with (1.6), [3, Corollary 4.4] yields

[TABLE]

Writing $g=\int_{0}^{\infty}\text{P}_{\tau}(\rho_{t}-\rho)\,\mathrm{d}\tau=-\int_{0}^{t}\rho_{\tau}\,\mathrm{d}\tau$ together with $|\int_{0}^{t}\rho_{\tau}\,\mathrm{d}\tau|\leq\|\rho\|_{\mathrm{L}^{\infty}}t$ gives (3.11). ∎

3.2. $\mathrm{L}^{q}$-type estimates

As we have seen in (2.8), we need a sharp control of the averaged Meyers’ norm \mathbb{E}\big{[}\|\nabla h^{n,t}\|^{2}_{\mathrm{L}^{\bar{q}}}\big{]}. This will be obtained as an immediate corollary of the following proposition, for more details see (3.101).

Proposition 3.3** ($\mathrm{L}^{q}$-estimates).**

Let $\{\mu^{n}\}_{n}$ be defined in (1.2) with point clouds satisfying Assumptions 1.1. Let $\bar{q}$ be the Meyers exponent given in Theorem A.3 for the operator $-\nabla\cdot\rho\nabla$. The solution444with Neumann boundary conditions in case $\mathcal{M}$ has a boundary $f^{n,t}\in{\rm\dot{H}^{1}}$ of

[TABLE]

satisfies:

[TABLE]

and a random variable $\mathcal{C}_{n,t}$ satisfying for some generic constant $C_{q}$ depending on $q$,

[TABLE]

Furthermore, if (1.14) holds with $\eta\geq 1$ then the assumption (1.13) can be dropped and the stochastic integrability can be improved up to losing a $\log(n)$ factor, namely

[TABLE]

and a random variable $\mathcal{D}_{n,t}$ satisfying for some generic constant $D_{q}$ depending on $q$,

[TABLE]

Proof.

We proceed in four steps. In the first step, we prove a representation formula for $(-\Delta)^{\frac{1}{2}}$ that we will use as the core tool in the next steps. In the second step, we compare the two operators $-\nabla\cdot\rho\nabla$ and $-\Delta$, with help of Meyers’ estimate recalled in Theorem A.3. Doing so, we then have to bound the $\mathrm{L}^{q}$-norms ($2\leq q\leq\bar{q}$) of the gradient of the solution555which belongs to any $\mathrm{L}^{q}$ for any $q<\infty$ from Calderón-Zygmund’ theory, see for instance [26] to the Poisson equation with r.h.s. $\mu^{n,t}-\rho_{t}$ and Neumann boundary conditions. We control all the norms by the $\mathrm{L}^{4}$-norm that, in turn, we estimate using the Riesz-transform bound (3.9) and following some ideas from [6, Lemma 3.17]. In the third and fourth steps, we control the bound previously obtained in expectation where our main tool is Assumption 1.1 and the concentration inequalities in Proposition A.1.

Step 1. A representation formula for $(-\Delta)^{\frac{1}{2}}$. We show that given $f\in\text{C}^{2}$ such that $n_{\mathcal{M}}\cdot\nabla f=0$ on $\partial\mathcal{M}$ we have

[TABLE]

Note that $(-\Delta)^{\frac{1}{2}}f$, defined in (3.10), is well defined in $\mathrm{L}^{2}$. Indeed, using the fact that $(\phi_{k},f)_{\mathrm{L}^{2}}=\frac{1}{\lambda_{k}}(\phi_{k},\Delta f)_{\mathrm{L}^{2}}$ and (3.8), we have for any $N\leq M<\infty$

[TABLE]

which vanishes as $N\uparrow\infty$ uniformly in $M$.

We now justify (3.16). Observe that since $n_{\mathcal{M}}\cdot\nabla f=0$ on $\partial\mathcal{M}$,

[TABLE]

which is a direct consequence of two integration by parts using the heat-kernel representation $\text{P}_{s}f=\int_{\mathcal{M}}p_{s}(\cdot,y)f(y)\,\mathrm{d}\mathrm{m}(y)$. Therefore

[TABLE]

where the last integral is well-defined in $\mathrm{L}^{2}$ since from (3.7)

[TABLE]

We then use the spectral decomposition of the heat semigroup (3.4) to get

[TABLE]

The combination of (3.18) and (3.19) yields for any $\eta\in\mathrm{L}^{2}$

[TABLE]

Using (3.8), we have

[TABLE]

so that we can exchange integration and summation in (3.20) to obtain

[TABLE]

which gives (3.16) by arbitrariness of $\eta$. Finally, note that the r. h. s. integral in (3.16) is absolutely convergent thanks to the integration by parts $\Delta\mathrm{P}_{\tau}f=\mathrm{P}_{\tau}\Delta f$ and the bounds on the heat kernel (3.3), so that it defines a function in $\text{C}^{0}$.

Step 2. Comparison with the solution of the Poisson equation. We claim that for any $q\in[2,\min\{\bar{q},4\}]$ and $p<\infty$

[TABLE]

Let $g^{n,t}\in{\rm\dot{H}^{1}}$ be the solution of the following Poisson equation

[TABLE]

Re-expressing the right-hand side of (3.12) as $\nabla\cdot\nabla g^{n,t}$, we apply Meyers’ estimate recalled in Theorem A.3 and Hölder’s inequality to obtain:

[TABLE]

We now introduce the Paley-Littlewood functional

[TABLE]

We recall that the inverse of $\mathcal{L}$ is continuous, see [51], namely

[TABLE]

Combining the Riesz transform bound (3.9) with (3.24) yields

[TABLE]

We now claim that

[TABLE]

which requires a special care when $\mathcal{M}$ has a boundary. We use the definition of $\text{P}_{s}$ in form of $\partial_{s}\text{P}_{s}=\Delta\text{P}_{s}$ to get

[TABLE]

Recalling that $n_{\mathcal{M}}\cdot\nabla g^{n,t}=0$, (3.17) implies that $\Delta\text{P}_{\tau}g^{n,t}=\text{P}_{\tau}\Delta g^{n,t}$ which, combined with (3.22) and (3.16), gives

[TABLE]

In particular, it implies that $n_{\mathcal{M}}\cdot\nabla(-\Delta)^{\frac{1}{2}}g^{n,t}=0$. Therefore, one can use once more (3.17) and, together with (3.28), (3.27) turns into

[TABLE]

Using a last time (3.17) in form of $\text{P}_{s}\Delta\text{P}_{\tau}(\mu^{n,t}-\rho_{t})=\Delta\text{P}_{s}\text{P}_{\tau}(\mu^{n,t}-\rho_{t})$ that we combine with the semigroup property $\text{P}_{t}\text{P}_{t^{\prime}}=\text{P}_{t+t^{\prime}}$ yields

[TABLE]

which, together with (3.29) and (3.16) leads to (3.26).

The combination of (3.23), (3.25) and (3.26) together with Minkowski’s inequality gives

[TABLE]

which is (3.21).

Step 3. Proof of (3.13). The estimate (3.13) is a consequence of (3.21) applied with $p=1$ and

[TABLE]

with

[TABLE]

Indeed, plugging (3.30) in (3.21) gives

[TABLE]

Using that

[TABLE]

and analogously

[TABLE]

(3.32) implies (3.13).

We now show (3.30). First, using the definition (1.2) of $\mu^{n}$, we have

[TABLE]

such that expanding the square gives

[TABLE]

The estimate (3.30) is then a consequence of

[TABLE]

Indeed, the first item of (3.34) immediately implies

[TABLE]

while, using the assumption (1.14) and that $t^{\frac{2}{\log(n)}}\lesssim 1$,

[TABLE]

It remains to prove (3.34) and we start with the first item. Here, we use the fact that

[TABLE]

This can be seen from (3.10): Using that $\mathrm{P}_{\frac{s}{2}}$ is an auto-adjoint operator in $\mathrm{L}^{2}$ and that from (3.19) we learn $\mathrm{P}_{\frac{s}{2}}\phi_{n}=e^{-\lambda_{n}\frac{s}{2}}\phi_{n}$, we have

[TABLE]

Hence, using the semi-group property $\text{P}_{s+t}=\text{P}_{\frac{s}{2}}\text{P}_{\frac{s}{2}+t}$, we deduce

[TABLE]

We now bound $\|\text{P}_{\frac{s}{2}+t}(\delta_{y}-\rho)\|^{2}_{\mathrm{L}^{2}}$ in two different ways. First, using the bounds (3.3) of the heat-kernel, we have by the triangle inequality

[TABLE]

yielding to the first alternative in the first item of (3.34). Second, applying Poincaré’s inequality yields

[TABLE]

yielding to the second alternative in the first item of (3.34).

We now turn to the second item of (3.34). Here, we make use of the representation formula (3.16) applied to $f=\mathrm{P}_{s+t}(\delta_{y}-\rho)$. Using the semi-group property $\text{P}_{\tau}\text{P}_{s+t}=\text{P}_{\tau+s+t}$, this takes the form

[TABLE]

A direct application of the heat-kernel bounds (3.3) leads to

[TABLE]

which is the first alternative in the second item of (3.34). For the second alternative, we write

[TABLE]

so that, using (3.3),

[TABLE]

Step 4. Proof of (3.15). Following the same computations done in the previous step, the estimate (3.15) is a consequence of (3.21) and the moment bounds

[TABLE]

together with the second item of (3.34) and Lemma A.2, where we recall that $\zeta$ is defined in (3.31).

We now show (3.37). Using (3.33) and the assumption $\eta\geq 1$, we apply the concentration inequality in Proposition A.1 to the effect of: for any $\lambda\geq 0$

[TABLE]

with

[TABLE]

We then obtain (3.37) combining (3.38) with (3.34) and (3.35), together with an application of the Layer-cake formula.

∎

3.3. Fluctuation estimates

This section is devoted to justify (2.12) needed to ensure the condition (2.3) with very high probability. Our result is in the spirit of [3, Theorem 3.3]. However, our strategy differs from [3, Theorem 3.3] and is based on Schauder’s theory, with an additional special care on the dependences on $\delta$. We briefly sketch the main ingredients of the fluctuation estimates (2.12). By linearity, it is enough to show (2.12) for $f^{n,t}_{\delta}\in\dot{\mathrm{H}}^{1}$ being the solution of

[TABLE]

We then make use of the chain rule to expand the equation in

[TABLE]

and we define the auxiliary problem

[TABLE]

such that the difference $v^{n,t}_{\delta}:=f^{n,t}_{\delta}-u^{n,t}_{\delta}$ solves $-\Delta v^{n,t}_{\delta}=\tfrac{1}{\rho_{\delta}}\nabla\rho_{\delta}\cdot\nabla f^{n,t}_{\delta}$. Doing so, on the one hand, we can control $u^{n,t}_{\delta}$ which can be handled using the explicit formula in terms of the heat-kernel, the explicit bounds on the latter, cf (3.3), and the regularity of $\frac{1}{\rho_{\delta}}$. On the other hand, we use Schauder’s estimates to control $v^{n,t}_{\delta}$ and $f^{n,t}_{\delta}$. Using the fact that those estimates depend polynomially on $\|\rho_{\delta}\|_{\text{C}^{0,\alpha}}$, we can keep track on the dependences on $\delta$ that we can optimize later on. We finally mention that in the case where $\mathcal{M}$ has a boundary, $u^{n,t}_{\delta}$ cannot be directly defined by (3.41) since the r. h. s. does not have zero mean. In order to also include this case, we add a zero order term in the equation, see (3.42).

Proposition 3.4** (Fluctuation estimates).**

Let $\{\mu^{n}\}_{n}$ be defined in (1.2) with point clouds satisfying Assumption 1.1. For any parameter $\delta\in(0,1)$ and $\upsilon>0$, we define666with Neumann boundary conditions in case $\mathcal{M}$ has a boundary $u^{n,t}_{\delta}\in\mathrm{H}^{1}$ weak solution of

[TABLE]

and the two events

[TABLE]

There exists $\kappa>0$ such that for any $\kappa_{1}>0$ and the choice

[TABLE]

the solution $f^{n,t}_{\delta}\in\dot{\mathrm{H}}^{1}$ of (3.39) satisfies

[TABLE]

Furthermore, there exists $\kappa_{2}>0$ depending on $\upsilon$ such that for the choice $t=t_{n}:=\frac{\log^{\kappa_{2}}(n)}{n}$, we have

[TABLE]

Proof.

The proof of Proposition 3.4 is split into two steps. In the first step, we prove (3.45) where our main tool is Schauder’s theory and an explicit formula for $u^{n,t}_{\delta}$ defined in (3.40). In the second step, we show (3.46), where our main tool is the concentration inequalities in Proposition A.1 and the explicit formula for $u^{n,t}_{\delta}$ used in the first step.

Step 1. Proof of (3.45). We define the difference v^{n,t}_{\delta}:=f^{n,t}_{\delta}-\Big{(}u^{n,t}_{\delta}-\int_{\mathcal{M}}u^{n,t}_{\delta}\Big{)}\in\dot{\text{H}}^{1} and note that from (3.40) and (3.42), it solves

[TABLE]

We claim that there exists $\kappa>0$ such that with the choices of $\delta$ in (3.44), we have

[TABLE]

The estimate (3.45) then follows from (3.48) and the triangle inequality.

We now prove (3.48). We apply Schauder’s estimate (see for instance [27] and [45, Chapter 10]) to (3.47) to obtain

[TABLE]

While the second r. h. s. term is directly of order of $\frac{1}{\log^{\upsilon}(n)}$ in $\mathcal{B}_{\delta,n}$, the first r. h. s. term requires additional attention. Using (1.6) and (3.3), $\rho_{\delta}$ satisfies

[TABLE]

and together with the algebraic property of $\|\cdot\|_{\text{C}^{0,\alpha}}$, we have

[TABLE]

The latter is bounded using Schauder’s estimate applied this time on (3.39) (knowing that the dependence on $\|\rho_{\delta}\|_{\text{C}^{0,\alpha}}$ is at most polynomial): there exists $\kappa>0$ such that

[TABLE]

which yields the following control of the first r. h. s. term of (3.49)

[TABLE]

which is of order of $\frac{1}{\log^{\upsilon-(\kappa+2)\kappa_{1}}(n)}$ in $\mathcal{A}_{n}$.

Step 2. Proof of (3.46). We provide the arguments for (3.46) in case that (1.14) holds for $\eta<1$ and the event $\mathcal{B}^{c}_{\delta,n}$ that we reduce to $\{\|\nabla^{2}u^{n,t}_{\delta}\|_{\mathrm{L}^{\infty}}\leq\tfrac{1}{\log^{\upsilon}(n)}\}$, the other cases as well as the case of $\nabla u^{n,t}_{\delta}$ follow by a straightforward adaptation. The estimate (3.46) is a consequence of

[TABLE]

for some constants $C_{1},C_{2},C_{3},C_{4}>0$. To see this, let $\eta_{n}$ be defined by

[TABLE]

By compactness of $\mathcal{M}$, we can find a $\eta_{n}$-net $\{x_{k}\}_{1\leq k\leq N}\subset\mathcal{M}$ with $N\lesssim\eta^{-2}_{n}$. We note that

[TABLE]

Indeed, if for any $k\in\{1,\cdots,N\}$, $|\partial^{2}_{ij}u^{n,t}_{\delta}(x_{k})|\leq\tfrac{1}{2\log^{\upsilon}(n)}$ then for any $x\in\mathcal{M}$ there exists $j\in\{1,\cdots,N\}$ such that

[TABLE]

Applying $\mathbb{P}$ on (3.54) yields

[TABLE]

Using (3.52) and

[TABLE]

which can be proven following the arguments leading to the second item of (3.59) below, this yields ${\mathbb{P}}(\{\|\nabla^{2}u^{n,t}_{\delta}\|_{\mathrm{L}^{\infty}}\leq\tfrac{1}{\log^{\upsilon}(n)}\})=o(\tfrac{1}{n^{m}})$ for any $m\in\mathbb{N}$ for the choice $\kappa_{2}>4\kappa_{1}+\frac{1}{\eta}+2\upsilon$.

We now prove (3.52). We first exploit the following explicit representation formula777a simple change of variable gives that the kernel $\tilde{p}_{s}$ associated to $1-\Delta$ is given by $\tilde{p}_{s}=e^{-s}p_{s}$ for $u^{n,t}_{\delta}$,

[TABLE]

that we expand, using the definition (1.2) of $\mu^{n}$, in form of

[TABLE]

Then applying the concentration inequalities in Proposition A.1, we get

[TABLE]

where

[TABLE]

The estimate (3.52) then follows from the three following estimates

[TABLE]

that we prove separately in the next three sub-steps.

Sub-step 2.1. Proof of the first item of (3.59). Splitting the time integral into $\int_{0}^{t}+\int_{t}^{\infty}$ and subtracting and adding back $\tfrac{1}{\rho_{\delta}(x)}$ in the first integral as well as using the semigroup property of $\{\mathrm{P}_{s}\}_{s}$ in form of

[TABLE]

we decompose $\omega$ into a regular-part $\mathcal{J}_{1}$ and a singular-part $\mathcal{J}_{2}$:

[TABLE]

For the regular-part $\mathcal{J}_{1}$, we apply directly the heat-kernel bounds (3.3) and the first item of (3.50), to obtain from the triangle inequality and Minkowski’s inequality

[TABLE]

The first r. h. s. term $\mathcal{R}_{1}(x)$ is dominated using directly the heat-kernel bounds (3.3) and (3.50)

[TABLE]

For the second r. h. s. side term $\mathcal{R}_{2}(x)$, we first simplify the $y$-integral. Using that

[TABLE]

we have by Jensen’s inequality and the heat-kernel bounds (3.3)

[TABLE]

Thus,

[TABLE]

The combination of (3.61), (3.62) and (3.63) yields

[TABLE]

We now turn to the singular-part $\mathcal{J}_{2}$. We first apply Minkowski’s inequality in form of

[TABLE]

We then simplify the $y$-integral. To this aim, we bound the integrand in $\mathrm{L}^{\infty}$ using the heat-kernel bounds (3.3), (3.50) and

[TABLE]

in form of

[TABLE]

This yields together with (3.66)

[TABLE]

To conclude, the combination of (3.60), (3.65) and (3.69) shows the first item of (3.59).

Sub-step 2.2. Proof of the second item of (3.59). We use the decomposition (3.60). For the regular-part $\mathcal{J}_{1}$, we argue as in (3.61) for the first term whereas the second-term is estimated using the heat-kernel bounds (3.3) in form of

[TABLE]

so that

[TABLE]

Hence,

[TABLE]

For the singular-part $\mathcal{J}_{2}$, we use the bound (3.68) which directly yields

[TABLE]

The combination of (3.60), (3.70) and (3.71) gives the second item of (3.59).

Sub-Step 2.3. Proof of the third item of (3.59). According to the first item of (3.59), it suffices to give the argument for the second term in the definition (3.58) of $v^{2}$. We use the assumption (1.14) together with the two first items of (3.59) in form of

[TABLE]

which concludes since $\limsup_{n\uparrow\infty}\delta^{\frac{4}{\log(n)}}t^{-\frac{1}{\log(n)}}\lesssim 1$.

∎

3.4. Contractivity estimates

This section is devoted to the control of the smoothing errors $W^{2}_{2}(\mu^{n,t},\mu^{n})$ and $W^{2}_{2}(\nu^{m,t},\nu^{m})$ for the particular choice of $t$ given in Proposition 3.4. The first result is in the spirit of [3, Theorem 5.2] that we extend in the case of non-uniformly distributed and correlated points. This extension requires a finer analysis of the error and the proof relies on Berry-Esseen type inequalities in the spirit of [12, Theorem 5].

Proposition 3.5** (Semigroup contraction for empirical measures).**

Let $\{\mu^{n}\}_{n}$ be defined in (1.2) with point clouds satisfying Assumption 1.1. Given $t$ such that Proposition 3.4 holds, we have

[TABLE]

for some random variable $\mathcal{C}_{n}$ satisfying for $C<\infty$

[TABLE]

Furthermore, if (1.14) holds with $\eta\geq 1$ then the assumption (1.13) can be dropped and the stochastic integrability can be improved up to losing a $\log(n)$ factor, namely

[TABLE]

for some random variable $\mathcal{D}_{n}$ satisfying for $D<\infty$

[TABLE]

Proof.

According to the fluctuation estimates in Proposition 3.4 together with $W^{2}_{2}(\mu^{n,t},\mu^{n})\leq(\text{diam}(\mathcal{M}))^{2}$, we can restrict the analysis in $\mathcal{A}_{n}$ defined in (3.43). Note that for $n$ large enough, (1.6) yields

[TABLE]

We split the proof into three steps. In the first step, we prove a Berry–Esseen type smoothing inequality for $W^{2}_{2}(\mu^{n,t},\mu^{n})$ which decomposes the error in a deterministic part involving $\rho$ and a random part involving the Fourier coefficients $\{\widehat{\mu^{n}}(k)\}_{k}$ of $\mu^{n}$. In the second step, we prove (3.72). In the third step, we control the fluctuations of $\{\widehat{\mu^{n}}(k)\}_{k}$ using the concentration inequalities in Proposition A.1 and deduce (3.73).

Step 1. Berry-Esseen type inequality. Recalling that we denote by $\{\lambda_{k},\phi_{k}\}_{k}$ the eigenvalues and eigenfunctions of $-\Delta$ respectively, we prove that

[TABLE]

where

[TABLE]

We first apply the triangle inequality and use the classical contractivity estimate in [22, Theorem 3] to get

[TABLE]

We then apply Peyre’s estimate [46] to the second r. h. s. , which takes the form

[TABLE]

Now, given an arbitrary $f$ such that

[TABLE]

we split

[TABLE]

For the first r. h. s. term of (3.79), we expand the integral using (3.4). Thus, together with the semigroup property of $\{\mathrm{P}_{t}\}_{t>0}$ and Cauchy-Schwarz’s inequality, we get

[TABLE]

Using that from (3.74) we have $\mu^{n,t+\frac{1}{n}}\geq\frac{\lambda}{2}$ and recalling (3.78), we get

[TABLE]

Furthermore, noticing that $\mathrm{P}_{t}\phi_{k}=e^{-t\lambda_{k}}\phi_{k}$ which implies that, since $\mathrm{P}_{t}$ is self-adjoint

[TABLE]

we obtain

[TABLE]

This leads to

[TABLE]

For the second r. h. s. of (3.79), we introduce $w\in\dot{{\rm H}}^{1}$ satisfying

[TABLE]

so that from an integration by parts, Cauchy-Schwarz’ inequality and the combination of (3.74) and (3.78), we obtain

[TABLE]

Using then the explicit formula $w=-\int_{\frac{1}{n}}^{t+\frac{1}{n}}\rho_{\tau}\,\mathrm{d}\tau$ together with $|\int_{\frac{1}{n}}^{t+\frac{1}{n}}\rho_{\tau}\,\mathrm{d}\tau|\leq\Lambda t$, we get

[TABLE]

so that

[TABLE]

The combination of (3.79), (3.80), (3.81) and (3.76) leads to (3.75).

Step 2. Proof of (3.72). According to (3.75), it remains to show that

[TABLE]

Writing

[TABLE]

we first expand the square in form of

[TABLE]

For the first r. h. s. term of (3.84), we use the normalisation $\|\phi_{k}\|_{\mathrm{L}^{2}}=1$ together with (1.6) to the effect of

[TABLE]

and get

[TABLE]

For the second r. h. s. term of (3.84), we use the definition of the $\beta$-mixing coefficient (1.12) together with the assumption (1.13) in form of

[TABLE]

The combination of (3.84), (3.86) and (3.87) yields

[TABLE]

It remains to show that

[TABLE]

We only treat the second l. h. s. term of (3.88), the first term is controlled the same way. For any $x,y\in\mathcal{M}$, we expand

[TABLE]

We then disintegrate using the spectral decomposition of the heat kernel (3.4) in form of

[TABLE]

so that we obtain

[TABLE]

Step 3. Proof of (3.73). It is a consequence of the following fluctuation estimates

[TABLE]

together with Lemma A.2. Indeed, applying Minkowski’s inequality followed by (3.89) and $\lambda^{\frac{1}{\log(n)}}_{k}e^{-\frac{2}{n}\lambda_{k}}\lesssim e^{-\frac{1}{n}\lambda_{k}}$ yield

[TABLE]

and (3.73) follows from

[TABLE]

which is obtained the same way as in (3.88) using additionally the trace formula (3.5).

We now prove (3.89). It follows from the estimate on the probability tails

[TABLE]

for some $C>0$, together with a simple application of the layer-cake representation.

To see (3.91), we use (3.83) together with Proposition A.1 to obtain

[TABLE]

with

[TABLE]

The estimate (3.91) is then a consequence of

[TABLE]

The first item of (3.93) has been treated in (3.6). For the second item of (3.93), we use (3.85) and, combined with (1.14), we obtain

[TABLE]

∎

3.5. Proof of Theorem 1.2: Approximation of the transport plan

We only give the arguments for (1.20), (1.21) is proved the same way using the corresponding results (3.15), (3.73) and (B.2) in the case $\eta>2$ (where some additional comments are given if necessary along the proof). We split the proof into four steps. In the first step, we display some preliminary estimates useful all along the proof. In the second step, we deal with the approximation error that occurs in the process of regularizing $\rho$ into $\rho_{\delta}$. In the third step, we estimate the $W_{2}$-distance for the regularized quantity using the quantitative stability result in [4, Theorem 3.2], splitting the estimates in small pieces that we control in the fourth step. We finally comment on the proof of Remark 1.4 and Theorem 1.3, which are obtained with similar techniques.

Step 1. Preliminary estimates.

**Heat kernel regularization. **The assumption $\rho\in\mathrm{H}^{\varepsilon}$ provides

[TABLE]

Indeed, using the definition of the heat-kernel, Minkowski’s inequality and the spectral decomposition (3.4), we have

[TABLE]

Noticing that

[TABLE]

we get that

[TABLE]

$\mathrm{L}^{\infty}$**-estimates. **Let $\kappa_{1}>0$ and $\upsilon>\max\{(\kappa+2)\kappa_{1},\frac{1}{\eta}\}$, where $\kappa$ is given in Proposition 3.4. For the given choices

[TABLE]

provided in Proposition 3.4, we define $h^{n,t}_{\delta}\in\dot{\mathrm{H}}^{1}$ the weak solution of

[TABLE]

Note that by linearity, one can decompose $h^{n,t}_{\delta}=h^{(1)n,t}_{\delta}-h^{(2)n,t}_{\delta}$ with

[TABLE]

so that, considering $u^{n,t}_{\delta}$ as in (3.42) and likewise $v^{n,t}_{\delta}$ with $\mu^{n,t}$ replaced by $\nu^{m,t}$ and defining

[TABLE]

as well as

[TABLE]

we deduce from Proposition 3.4 and $m=m(n)\underset{n\uparrow\infty}{\sim}qn$ that

[TABLE]

$\mathrm{L}^{q}$**-estimates. **A similar decomposition as in (3.97) of (1.18) together with (3.13) and the choice of $t$ in (3.95) yields

[TABLE]

where $\mathcal{C}_{n}$ denotes, all along the proof, a random variable which satisfies (3.14) and may change from line to line.

$\mathrm{L}^{2}$**-regularization error. **Note that from (3.96) and (1.18)

[TABLE]

so that from an energy estimate, Hölder’s inequality, (3.101) and (1.6), we obtain

[TABLE]

where $\bar{q}$ denotes the Meyers’ exponent of the operator $-\nabla\cdot\rho\nabla$, see Theorem A.3.

Finally, since $\inf_{\pi}W^{2}_{2}(\pi,\gamma^{n,t})\leq(\text{diam}(\mathcal{M}))^{2}$ and (3.46) holds, we can restrict our analysis in $\mathcal{A}_{n}\cap\mathcal{B}_{\delta,n}$ that we do for the rest of the proof.

Step 2. Regularization error. We show that (3.103) survives when measuring the $W_{2}$-distance, namely

[TABLE]

Using the coupling \big{(}(\mathrm{Id},\exp(\nabla h^{n,t}_{\delta})),(\mathrm{Id},\exp(\nabla h^{n,t})\big{)}_{\#}\mu^{n,t} as a competitor in (3.1) and the fact that $\|\mu^{n,t}\|_{\mathrm{L}^{\infty}}\lesssim 1$ in $\mathcal{A}_{n}$, we have

[TABLE]

We then claim that

[TABLE]

which, combined with (3.105) and (3.103) yields (3.104).

We now justify (3.106). The difficulty arises from the fact that $\exp$ is not globally Lipschitz. To overcome this, we define

[TABLE]

for a given $\varsigma$ fixed later, and we split

[TABLE]

For the first right-hand side integral, note that from the choice of $C_{n}$ and (3.103), we can choose $\varsigma\ll 1$ uniformly in $n$ such that in $\mathrm{E}_{n}$ the quantity $|\nabla h^{n,t}_{\delta}-\nabla h^{n,t}|$ can be made arbitrary small. Since $\exp$ is Lipschitz-continuous in a neighborhood of the null vector, we deduce

[TABLE]

For the second right-hand side term, we simply apply Markov’s inequality in form of

[TABLE]

The combination of the two previous estimates gives (3.106).

To prove (1.21), we need to control arbitrary $p$-moments, according to Lemma A.2. The argument above can be easily adapted in this case by considering

[TABLE]

We then follow the same argument, choosing $\varsigma^{-1}=O(\sqrt{p})$.

Step 3. Quantitative stability. We show that

[TABLE]

where we recall that $\gamma^{n,t}_{\delta}$ is defined in (3.104).

Let $\pi$ be a coupling between $\mu^{n}$ and $\nu^{m}$. We introduce a regularization parameter $s<1$ and, smoothing the measure $\mu^{n}$ into $\mu^{n,s}:=\text{P}_{s}\mu^{n}$, the optimal transport plan $\pi^{n,s}$ from $\mu^{n,s}$ to $\nu^{m}$ is represented by a transport map $T^{n,s}$, according to McCann’s theorem [39], that is

[TABLE]

We then apply the triangle inequality in form of

[TABLE]

First, using (3.100), $\nabla h^{n,t}_{\delta}$ is Lipschitz-continuous and $\|\nabla h^{n,t}_{\delta}\|_{\mathrm{L}^{\infty}}$ can be made as small as possible for $n$ large. Since $\exp$ is Lipschitz-continuous in a neighborhood of the null vector, we learn from Lemma 3.1 that

[TABLE]

Second, we build a competitor for the second right-hand side term of (3.108): Defining

[TABLE]

we have

[TABLE]

Using again (3.100) we can apply, for large $n$, the quantitative stability result of transport maps, Theorem 2.1, to $\mu_{1}=\nu^{m}$, $\mu_{2}=\exp(\nabla h^{n,t}_{\delta})_{\#}\mu^{n,s}$ and $\nu=\mu^{n,s}$ to the effect of

[TABLE]

which turns into, using the triangle inequality,

[TABLE]

The combination of (3.108), (3.109) and (3.110) yields

[TABLE]

Since $\mu^{n,s}\underset{s\downarrow 0}{\rightharpoonup}\mu^{n}$, and consequently (up to extracting a subsequence) $\pi^{n,s}\underset{s\downarrow 0}{\rightharpoonup}\pi$, for some optimal transport plan $\pi$, according to the qualitative stability result [56, Theorem 5.20], we can pass to the limit as $s\downarrow 0$ in (3.111) which leads to (3.107).

Step 4. Proof of (1.20). We now fix $\kappa_{1}=\frac{1}{\eta}(\tfrac{\bar{q}}{2})\frac{1}{2\varepsilon}+1$ such that, applying (3.94), the regularization error (3.104) turns into, recalling that $\delta$ is given by (3.95),

[TABLE]

It remains to show that

[TABLE]

which together with (3.94) and (3.112) leads to (1.20). To show (3.113), we control each terms of (3.107) separately.

The three last terms are controlled using the contractivity estimate (3.72) and (B.1) which gives

[TABLE]

For the first two terms, we argue that

[TABLE]

which combined with (B.1) leads to (3.113).

Let us define the curve $\eta:s\in[0,1]\mapsto\eta_{s}:=s\mu^{n,t}+(1-s)\nu^{m,t}$ and note that from (1.18) we have

[TABLE]

Applying Benamou-Brenier’ theorem [10], we learn that

[TABLE]

Next, using that

[TABLE]

and applying [3, Proposition A.1] together with Hölder’s inequality yields

[TABLE]

Using Meyers’ estimate of Proposition A.3 to (3.102) together with (1.6) and (3.101) provides

[TABLE]

which, combined with (3.115), yields (3.114).

We finally point out that, in the case $\eta>2$, we use (3.73) and (B.2) and the same computations lead to (1.21).

Step 5. Proof of Theorem 1.3 and Remark 1.4. The proof of Theorem 1.3 follows the same strategy with the main difference that Step 3 is now dropped and Theorem 2.1 is directly applied with $\mu_{1}=\rho$, $\nu=\mu^{n}$ and $\mu_{2}=\exp(\nabla f^{n,t}_{\delta})_{\#}\mu^{n}$ where $f^{n,t}_{\delta}$ solves

[TABLE]

The improvement of Remark 1.4 follows from the improved contractivity estimate (3.73): Under the assumption 1.25, we have (keeping the notations as in Proposition 3.5)

[TABLE]

i.e. we do not have the loss $\log^{\frac{1}{\eta}}(n)$ in (3.73). Inspecting the proof of (3.73), the loss $\log^{\frac{1}{\eta}}(n)$ comes from estimating $v^{2}$ defined in (3.92). We obtain (3.116) by simply using (1.25) and (1.11) to upgrade the second item of(3.93) into

[TABLE]

Appendix A Probabilistic and PDE tools

This section is devoted to recall some probabilistic and analytical tools needed in the proofs. We first recall some concentration inequalities for sequences of random variable satisfying Assumptions 1.1. Originally proved for i.i.d. samples, see for instance [19, Theorem 3.6 & 3.7], the proofs in the correlated case can be found in [42, Theorem 1] and [41, Theorem 2].

Proposition A.1**.**

Let $n\in\mathbb{N}$, $M>0$, $\{X_{i}\}_{i}$ be a family of centred random variables such that $\sup_{i\geq 1}|X_{i}|\leq M$ for which (1.14) holds.

For any $\lambda>0$, it holds for some constants $(C_{i})_{i\in\{1,\cdot,5\}}$ depending on $a$, $b$ :

• (i)

If $\eta<1$,

[TABLE]

with

[TABLE]

• (ii)

If $\eta=1$,

[TABLE]

We then recall the link between algebraic moments and exponential moments. The proof is a direct consequence of the Taylor expansion of the exponential function.

Lemma A.2**.**

Let $X$ be a non-negative random variable. The following two statements are equivalent:

• (i)

There exists $C_{1}>0$ such that

[TABLE]

• (ii)

There exists $C_{2}>0$ such that

[TABLE]

We conclude this section by recalling the standard Meyers’ estimate for elliptic equations in divergence form, see for instance the original paper [43].

Theorem A.3** (Meyers estimate).**

Let $a:\mathcal{M}\rightarrow\mathbb{R}^{2\times 2}$ be measurable and uniformly elliptic. Consider $u\in\mathrm{H}^{1}$ the solution of the Neumann boundary problem

[TABLE]

for some $g\in\mathrm{L}^{q}$ with $q>2$. There exists $2<\bar{q}<q$ such that

[TABLE]

Appendix B Matching cost for point clouds

This section is devoted to recall the upper bounds on the matching cost, results which can be found in [12, Theorem 2] under mild $\beta$-mixing conditions. The case of Markov chains have been studied in [47, 23] where sharp upper bounds are obtained. We include a short proof for convenience.

Proposition B.1** (Matching cost).**

Let $\rho$ satisfying (1.6) and $\{\mu^{n}\}_{n}$ be defined in (1.2) with points cloud satisfying the Assumptions 1.1 or in the class of Markov chains satisfying the Assumptions 1.5. There exists a constant $C>0$ such that

[TABLE]

Furthermore, if (1.14) holds with $\eta\geq 1$ then the assumption (1.13) can be dropped and the stochastic integrability can be improved up to losing a $\log(n)$ factor, namely

[TABLE]

Proof.

Note that the proof of (B.1) can be found in [12, Theorem 2] when the point cloud satisfies the Assumptions 1.1. We first show how (B.1) can be extended to point clouds which are sampled from a Markov chain satisfying the Assumptions 1.5. Second, we show how the stochastic integrability can be improved to (B.2) when (1.14) holds with $\eta\geq 1$.

Step 1. Markov chains case. Recall that a Markov chain satisfying the Assumptions 1.5 admits an absolutely continuous invariant measure of the form $\mu_{\infty}=\rho\,\mathrm{d}\mathrm{m}$ with $\rho$ satisfying (1.6), that is $\lambda\leq\rho\leq\Lambda$. Recalling that we denote by $\{\lambda_{n},\phi_{n}\}_{n}$ the set of eigenvalues and normalized eigenfunctions of the Laplace-Beltrami operator $-\Delta$ on $\mathcal{M}$, we have by definition (1.2) of $\mu^{n}$, for any $k\geq 1$

[TABLE]

where we use interchangeably the notation $\widehat{\rho}(k)=\int\phi_{k}\rho_{k}\,\mathrm{d}\mathrm{m}=\mu_{\infty}(\phi_{k})$. Using the Berry-Esseen smoothing inequality [12, Theorem 5] together with (B.3), we get

[TABLE]

We now estimate the last two terms of (B.4) separately and we start with the third one. Using (1.33) and (3.93), we have

[TABLE]

where we recall that $\sigma\in(0,1)$. Thus, using in addition (3.5), we get

[TABLE]

We now turn to the second term of (B.4). Expanding the square provides

[TABLE]

We now estimate the two terms on the right hand side of (B.6). For the first term, an easy induction argument combining (1.29) and (1.28) show that for any $n\geq 1$, $\mathbb{P}_{X_{n}}\ll\mathrm{m}$ with $\lambda\leq\frac{\mathrm{d}\mathbb{P}_{X_{n}}}{\mathrm{d}\mathrm{m}}\leq\Lambda$. Therefore, we have

[TABLE]

and we deduce

[TABLE]

For the second term, we use (1.32) to obtain

[TABLE]

Combining the latter with (B.4), (B.5), (B.6) and (B.7) yields

[TABLE]

We finally conclude similarly as for (3.88).

Step 2. Higher stochastic integrability. We now prove (B.2). We argue using the moment estimate (3.89) which, together with Minkowski’s inequality and $\lambda_{k}^{\frac{1}{\log(n)}}e^{-\frac{1}{n}\lambda_{k}}\stackrel{{\scriptstyle<}}{{\sim}}e^{-\frac{1}{2n}\lambda_{k}}$ implies

[TABLE]

Finally, combining the latter with the Berry-Esseen smoothing inequality [12, Theorem 5] and arguing similarly as for (3.88) yields (B.2) thanks to Proposition A.2. ∎

Appendix C Proof for the class of Markov chains

We provide in this Section the arguments for extending Theorem 1.2 and Theorem 1.3 to the class of Markov chains introduced in Section 1.4. The proof follows the lines of the proof of Theorem 1.2, where the main difference is that we drop the assumption that the point clouds is identically distributed. That affects the proofs of the main ingredients (we recall that the scaling of the cost has already be proven in Proposition B.1), namely the $\mathrm{L}^{q}$ estimates in Proposition 3.3, the fluctuation estimates in Proposition 3.4 and the contractivity estimates in Proposition 3.5. We show in the following how to adapt the proofs for a given Markov chain $\{X_{n}\}_{n\geq 1}$ satisfying Assumption 1.5. In the following, we recall that $\mu_{\infty}=\rho\,\mathrm{d}\mathrm{m}$ denotes the unique invariant measure of the chain. We split the proof in three steps

Step 1. $\mathrm{L}^{q}$ estimates. We have to understand the extra error term coming from the deviation of $\mathbb{E}[\mu^{n}]$ from $\mu_{\infty}$. In view of (3.21), it is

[TABLE]

Using the definition (1.2) of $\mu^{n}$, the convergence to equilibrium (1.33) applied to $f=(-s\Delta)^{\frac{1}{2}}p_{s+t}(x,\cdot)$ and the heat-kernel estimates (3.3), we have for any $s\geq 0$ and $x\in\mathcal{M}$

[TABLE]

so that, recalling $t=\frac{\log^{\kappa}(n)}{n}$,

[TABLE]

Step 2. Fluctuation estimates. Here, the distribution of the Markov chain affects the concentration estimate (3.52). We show that, defining

[TABLE]

we have

[TABLE]

where the r.h.s can be estimated following the lines of the proof of (3.52). As before, we investigate the extra term coming from the deviation of $\mathbb{E}[\mu^{n}]$ from $\mu_{\infty}$. The estimate (C.3) follows from

[TABLE]

We argue as in (3.60), decomposing $u^{n,t}_{\delta}-\bar{u}^{n,t}_{\delta}$ into a regular-part and a singular part: for any $x\in\mathcal{M}$

[TABLE]

To estimate the second r.h.s integral of (C.5), we use (3.68). For the first r.h.s integral, that we denote by $\mathcal{J}$, we use the definition (1.2) of $\mu^{n}$, the convergence to equilibrium (1.33) and the heat-kernel bounds (3.3) to obtain

[TABLE]

Step 3. Contractivity estimate. Here, the law of the Markov chain affects the estimate (3.82). The extra error term coming from the deviation of $\mathbb{E}[\mu^{n}]$ from $\mu_{\infty}$ reads

[TABLE]

Using the definition (1.2) of $\mu^{n}$ and the convergence to equilibrium (1.33) applied with $f=\phi_{k}$ and the bound on the eigenfunctions (3.6), we have for any $k\leq n$

[TABLE]

so that, using the trace formula (3.5) and the heat-kernel estimates (3.3), we deduce

[TABLE]

Acknowledgments

The authors warmly thank Lorenzo Dello Schiavo, Antonio Agresti and Martin Huesmann for useful discussions and fruitful comments.