Approximation of Optimal Transport problems with marginal moments constraints

TL;DR

This paper introduces a method to approximate optimal transport problems using moment constraints, resulting in finite discrete solutions that scale well with problem dimension and converge rapidly as the number of moments increases.

Contribution

It demonstrates that moment-constrained OT problems can be solved with finite discrete measures, providing convergence rates and practical algorithms for approximation.

Findings

Finite discrete measures solve MCOT problems.

Number of points scales linearly with marginal laws.

Convergence rates of O(1/n) and O(1/n^2) for moments.

Abstract

Optimal Transport (OT) problems arise in a wide range of applications, from physics to economics. Getting numerical approximate solution of these problems is a challenging issue of practical importance. In this work, we investigate the relaxation of the OT problem when the marginal constraints are replaced by some moment constraints. Using Tchakaloff's theorem, we show that the Moment Constrained Optimal Transport problem (MCOT) is achieved by a finite discrete measure. Interestingly, for multimarginal OT problems, the number of points weighted by this measure scales linearly with the number of marginal laws, which is encouraging to bypass the curse of dimension. This approximation method is also relevant for Martingale OT problems. We show the convergence of the MCOT problem toward the corresponding OT problem. In some fundamental cases, we obtain rates of convergence in $O(1/n)$ or $O(1/n^2)$ where $n$ is the number of moments, which illustrates the role of the moment functions. Last, we present algorithms exploiting the fact that the MCOT is reached by a finite discrete measure and provide numerical examples of approximations.

Equations

$$\inf {\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}c(x,y)\mathrm{d}\pi (x,y)$$

$$\inf \sum _{1\le i,j\le N}{\pi }_{ij}c(x^i,x^j)$$

$$\forall 1\le i,j\le N,\sum _{j=1}^N{\pi }_{ij}={\mu }_i\text{mbox}and\sum _{i=1}^N{\pi }_{ij}={\nu }_j,$$

$$\inf {\int }_{({\mathbb{R}}^d{)}^M}c(x_1,\cdots ,x_M)\mathrm{d}\pi (x_1,\cdots ,x_M)$$

$$\inf {\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}c(x,y)d\pi (x,y)$$

$${\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}{\varphi }_i(x)\mathrm{d}\pi (x,y)={\int }_{{\mathbb{R}}^d}{\varphi }_i(x)\mathrm{d}\mu (x)\text{mbox}and{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}{\varphi }_j(y)\mathrm{d}\pi (x,y)={\int }_{{\mathbb{R}}^d}{\varphi }_j(y)\mathrm{d}\nu (y).$$

$$\Pi (\mu ,\nu ):=\left\{\pi \in \mathcal{P}(\mathcal{X}\times \mathcal{Y}):{\int }_{\mathcal{X}}\mathrm{d}\pi (x,y)=\mathrm{d}\nu (y),{\int }_{\mathcal{Y}}\mathrm{d}\pi (x,y)=\mathrm{d}\mu (x)\right\},$$

$$I=\inf \left\{{\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}\pi (x,y):\pi \in \Pi (\mu ,\nu )\right\}.$$

$$\displaystyle\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq M},(\psi_{n})_{1\leq n\leq N}):=\Bigg{\{}\pi\in\mathcal{P}(\mathcal{X}\times\mathcal{Y}):$$

$$\forall 1\le m\le M,1\le n\le N,{\int }_{\mathcal{X}\times \mathcal{Y}}|{\varphi }_m(x)|+|{\psi }_n(y)|\mathrm{d}\pi (x,y)<\infty ,$$

$$\displaystyle\int_{\mathcal{X}\times\mathcal{Y}}\phi_{m}(x){\mathrm{d}}\pi(x,y)=\int_{\mathcal{X}}\phi_{m}(x){\mathrm{d}}\mu(x),\;\int_{\mathcal{X}\times\mathcal{Y}}\psi_{n}(y){\mathrm{d}}\pi(x,y)=\int_{\mathcal{X}}\psi_{n}(y){\mathrm{d}}\mu(x)\Bigg{\}},$$

$$I^{M,N}=\inf \left\{{\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}\pi (x,y):\pi \in \Pi (\mu ,\nu ;({\varphi }_m{)}_{1\le m\le M},({\psi }_n{)}_{1\le n\le N})\right\}.$$

$${\mu }_m:={\int }_{\mathcal{X}}{\varphi }_m\mathrm{d}\mu \text{and}{\nu }_n:={\int }_{\mathcal{Y}}{\psi }_n\mathrm{d}\nu .$$

$${\int }_{{\mathbb{R}}^{N_0}}\Vert u\Vert \mathrm{d}\Lambda \#\pi (u)={\int }_{{\mathbb{R}}^d}\Vert \Lambda (z)\Vert \mathrm{d}\pi (z)<\infty ,$$

$$\forall 1\le i\le N_0,{\int }_{{\mathbb{R}}^d}{\Lambda }_i(z)\mathrm{d}\pi (z)=\sum _{k=1}^K{p}_k{\Lambda }_i(z_k),$$

$$\exists \gamma \in \Pi (\mu ,\nu ;({\varphi }_m{)}_{1\le m\le M},({\psi }_n{)}_{1\le n\le N}),{\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}\gamma (x,y)<\infty .$$

$$\sum _{k=1}^{2N+1}{w}_k{\delta }_{(x_k,y_k)}\in \Pi (\mu ,\nu ;({\varphi }_m{)}_{1\le m\le N},({\psi }_n{)}_{1\le n\le N})\text{ and }\sum _{k=1}^{2N+1}{w}_{k}c(x_k,y_k)<\infty .$$

$$|{\varphi }_n(x)|\le C(1+{\theta }_{\mu }(|x|){)}^s\text{mbox}and|{\psi }_n(y)|\le C(1+{\theta }_{\nu }(|y|){)}^s.$$

$$I_A^N=\underset{\begin{array}{c}\pi \in \Pi (\mu ,\nu ;({\varphi }_m{)}_{1\le m\le N},({\psi }_n{)}_{1\le n\le N})\\ {\int }_{\mathcal{X}\times \mathcal{Y}}({\theta }_{\mu }(|x|)+{\theta }_{\nu }(|y|))\mathrm{d}\pi (x,y)\le A\end{array}}{\inf }{\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}\pi (x,y).$$

$$A_0^{\prime }={\int }_{\mathcal{X}}{\theta }_{\mu }(|x|)\mathrm{d}\mu (x)+{\int }_{\mathcal{Y}}{\theta }_{\nu }(|y|)\mathrm{d}\nu (y)<\infty ,$$

$$I_A^N=\underset{\begin{array}{c}\pi \in \Pi (\mu ,\nu ;({\varphi }_m{)}_{1\le m\le N},({\psi }_n{)}_{1\le n\le N})\end{array}}{\inf }{\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}\pi (x,y).$$

$$\Lambda:\left\{\begin{array}[]{ccc}\mathcal{X}\times\mathcal{Y}&\to&\mathbb{R}^{2N+2}\\ (x,y)&\mapsto&\begin{pmatrix}\phi(x)\\ \psi(y)\\ 1\\ c(x,y)\end{pmatrix}\\ \end{array}\right.$$

$${\int }_{\mathcal{X}\times \mathcal{Y}}\Lambda (x,y)\mathrm{d}\gamma (x,y)=\sum _{k=1}^K{w}_k\Lambda (x_k,y_k).$$

$${\int }_{\mathcal{X}\times \mathcal{Y}}\left({\theta }_{\mu }(|x|)+{\theta }_{\nu }(|y|)\right)\mathrm{d}\tilde{\gamma }(x,y)<\infty .$$

$${\mathcal{K}}_1=\left\{x\in \mathcal{X},\text{s.t.}|x|\le M_1\right\},{\mathcal{K}}_2=\left\{y\in \mathcal{Y},\text{s.t.}|y|\le M_2\right\}.$$

$$\begin{array}{rl}{\pi }_l(({\mathcal{K}}_1\times {\mathcal{K}}_2{)}^c)& ={\int }_{\mathcal{X}\times \mathcal{Y}}{1}_{(x,y)\notin {\mathcal{K}}_1\times {\mathcal{K}}_2}\mathrm{d}{\pi }_l(x,y)\le {\int }_{\mathcal{X}\times \mathcal{Y}}{1}_{x\notin {\mathcal{K}}_1}+1_{y\notin {\mathcal{K}}_2}\mathrm{d}{\pi }_l(x,y)\\ & \le {\int }_{\mathcal{X}\times \mathcal{Y}}\frac{{\theta }_{\mu }(|x|)}{{\theta }_{\mu }(M_1)}+\frac{{\theta }_{\nu }(|y|)}{{\theta }_{\nu }(M_2)}\mathrm{d}{\pi }_l(x,y)\le \frac{A}{\min ({\theta }_{\mu }(M_1),{\theta }_{\nu }(M_2))},\end{array}$$

$$I_A^N\le {\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}{\pi }_{\infty }(x,y).$$

$$\underset{l\to \infty }{\mathrm{liminf}}\mathbb{E}\left[c(X_l,Y_l)\right]\ge \mathbb{E}\left[c(X_{\infty },Y_{\infty })\right].$$

$$I_A^N\ge \mathbb{E}\left[c({\tilde{X}}_{\infty },{\tilde{Y}}_{\infty })\right]={\int }_{\mathcal{X}\times \mathcal{Y}}c(x,y)\mathrm{d}{\pi }_{\infty }(x,y).$$

$$c:\left\{\begin{array}[]{ccl}\mathbb{R}\times\mathbb{R}&\to&\mathbb{R}_{+}\\ (x,y)&\mapsto&(x-y)^{2}+\varphi(|x|)+\varphi(|y|),\end{array}\right.$$

$$I=\underset{\begin{array}{c}\pi \in \mathcal{P}(\mathbb{R}\times \mathbb{R})\\ {\int }_{\mathbb{R}}x\mathrm{d}\pi (x,y)=1\\ {\int }_{\mathbb{R}}y\mathrm{d}\pi (x,y)=1\end{array}}{\inf }\left\{{\int }_{\mathbb{R}\times \mathbb{R}}c(x,y)\mathrm{d}\pi (x,y)\right\}.$$

Full text

Approximation of Optimal Transport problems with marginal moments constraints

Aurélien Alfonsi

  

Rafaël Coyaud

  

Virginie Ehrlacher

  

Damiano Lombardi

**Abstract: ** Optimal Transport (OT) problems arise in a wide range of applications, from physics to economics. Getting numerical approximate solution of these problems is a challenging issue of practical importance. In this work, we investigate the relaxation of the OT problem when the marginal constraints are replaced by some moment constraints. Using Tchakaloff’s theorem, we show that the Moment Constrained Optimal Transport problem (MCOT) is achieved by a finite discrete measure. Interestingly, for multimarginal OT problems, the number of points weighted by this measure scales linearly with the number of marginal laws, which is encouraging to bypass the curse of dimension. This approximation method is also relevant for Martingale OT problems. We show the convergence of the MCOT problem toward the corresponding OT problem. In some fundamental cases, we obtain rates of convergence in $O(1/n)$ or $O(1/n^{2})$ where $n$ is the number of moments, which illustrates the role of the moment functions. Last, we present algorithms exploiting the fact that the MCOT is reached by a finite discrete measure and provide numerical examples of approximations.

1 Introduction

The aim of this paper is to investigate a new direction to approximate optimal transport problems [36, 31]. Such problems arise in a variety of application fields ranging from economy [21, 14] to quantum chemistry [16] or machine learning [29] for instance. The simplest prototypical example of optimal transport problem is the two-marginal Kantorovich problem, which reads as follows: for some $d\in\mathbb{N}^{*}$, let $\mu$ and $\nu$ be two probability measures on $\mathbb{R}^{d}$ and consider the optimization problem

[TABLE]

where $c$ is a non-negative lower semi-continuous cost function defined on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ and where the infimum runs on the set of probability measures $\pi$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginal laws $\mu$ and $\nu$.

The most straightforward approach for the resolution of problems of the form (1.0.1) consists in introducing discretizations of the state spaces, which are fixed a priori. More precisely, $N$ points $x^{1},\cdots,x^{N}\in\mathbb{R}^{d}$ are chosen a priori and fixed, marginal laws $\mu$ and $\nu$ are approximated by discrete measures of the form $\mu\approx\sum_{i=1}^{N}\mu_{i}\delta_{x^{i}}$ and $\nu\approx\sum_{i=1}^{N}\nu_{i}\delta_{x^{i}}$ with some non-negative coefficients $\mu_{i}$ and $\nu_{i}$ for $1\leq i\leq N$. An optimal measure $\pi$ minimizing (1.0.1) is then approximated by a discrete measure $\pi\approx\sum_{1\leq i,j\leq N}\pi_{ij}\delta_{x^{i},x^{j}}$ where the non-negative coefficients $(\pi_{ij})_{1\leq i,j\leq N}\in\mathbb{R}_{+}^{N^{2}}$ are solution to the optimization problem

[TABLE]

and satisfy the following discrete marginal constraints:

[TABLE]

which boils down to a classical linear programming problem, which becomes computationally prohibitive when $N$ is large.

Several numerical methods have already been proposed in the literature for the resolution of optimal transport problems at a lower computational cost. Most of them rely on an a priori discretization of the state spaces as presented above. One of the most successful approach consists in minimizing a regularized cost involving the Kullback-Leibler divergence (or relative entropy) via iterative Bregman projections: the so-called Sinkhorn algorithm [11, 28, 34]. Let us also mention other approaches such as the auction algorithm [13], numerical methods based on Laguerre cells [22], multiscale algorithms [27, 32] and augmented Lagrangian methods using the Benamou-Brenier dynamic formulation [9, 10].

In this work, we are also interested in studying multi-marginal and martingale-constrained optimal transport problems.

Multimarginal optimal transport problems arise in a wide variety of contexts [36, 31], like for instance the computation of Wasserstein barycenters [1] or the approximation of the correlation energy for strongly correlated systems in quantum chemistry [33, 15, 16]. Such problems read as follows: let $M\in\mathbb{N}^{*}$ and $\mu_{1},\cdots,\mu_{M}$ be $M$ probability measures on $\mathbb{R}^{d}$ and consider the optimization problem

[TABLE]

where $c$ is a lower semi-continuous cost function defined on $(\mathbb{R}^{d})^{M}$ and where the infimum runs on the set of probability measures $\pi$ on $(\mathbb{R}^{d})^{M}$ with marginal laws given by $\mu_{1},\cdots,\mu_{M}$. Approximations of such multi-marginal problems on discrete state spaces can be introduced in a similar way to (1.0.2), leading to a linear programming problem of size $N^{M}$. For large values of $M$, such discretized problems become intractable numerically. The most successful method up to now for solving such problems, which avoids this curse of dimensionality, is based on an entropic regularization together with the Sinkhorn algorithm [11, 28].

Constrained martingale transport arise in problems related to finance [7]. Few numerical methods have been proposed so far for the resolution of such problems. In [2, 3], algorithms using sampling techniques preserving the convex order is proposed, which enables then to use linear programming solvers. Algorithms making use of an entropy regularization and the Sinkhorn algorithm have been studied in [17, 23].

In this paper, we consider an alternative direction to approximate optimal transport problems, with a view to the design of numerical schemes. In this approach, the state space is *not *discretized, but the approximation consists in relaxing the marginal laws constraints (or the martingale constraint) of the original problem into a finite number of moment constraints against some well-chosen test functions. More precisely, in the case of Problem (1.0.1), let us introduce some real-valued bounded functions $\phi_{1},\cdots,\phi_{N}$ defined on $\mathbb{R}^{d}$, which are called hereafter test functions, and consider the approximate optimal transport problem, called hereafter the Moment Constrained Optimal Transport (MCOT) problem:

[TABLE]

where the infimum runs over the set of probability measures $\pi$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ satisfying for all $1\leq i,j\leq N$,

[TABLE]

The aim of this paper is to study the properties of this alternative approximation of optimal transport problems, and its generalization for multi-marginal and martingale-constrained optimal transport problems. A remarkable feature of this approximation is that it circumvents the curse of dimensionality with respect to the number of marginal laws in the case of multimarginal optimal transport problems. Note that in the martingale constrained case, our method enables to consider the original formulation of the financial problem that has moment constraints (see for instance Example 2.6 of [24]), which is not the case of the previous methods.

Our first contribution in this paper is to characterize some minimizers of the MCOT problem. Thanks to the Tchakaloff theorem, we prove that, under suitable assumptions, there exists at least one minimizer which is a discrete measure charging a finite number of points. Interestingly, in the multi-marginal case, the number of charged points scales at most linearly in the number of marginals. In the particular case of problems issued from quantum chemistry applications, due to the inherent symmetries of the system, the number of charged points is independent from the number of marginals, and only scales linearly with the number of imposed moments. This formulation of the multimarginal optimal transport problem thus does not suffer from the curse of dimensionality. The result obtained in the quantum chemistry case is close in spirit to the one of [20] where the authors studied a multimarginal Kantorovich problem with Coulomb cost on finite state spaces.

These considerations motivate us to consider a new family of algorithms for the resolution of multi-marginal and martingale constrained optimal transport problems, in which an optimal measure is approximated by a discrete measure charging a relatively low number of points. The points and weights of this discrete measure are then optimized in order to satisfy a finite number of moment constraints and to minimize the cost of the original optimal transport problem.

Of course, another interesting issue consists in determining how the choice of the particular test functions influences the quality of the approximation with respect to the exact optimal transport problem. In this paper, we prove interesting one-dimensional results in this direction. More precisely, for piecewise affine test functions defined on a compact interval, and for the two-marginal optimal transport problems involved in the computation of the $W_{2}$ or the $W_{1}$ distance between two sufficiently regular measures, the convergence of the approximate optimal cost with respect to the optimal cost scales like $\mathcal{O}\left(\frac{1}{N^{2}}\right)$ where $N$ is the number of test functions. These results indicate that the choice of appropriate test functions has an influence on the rate of convergence of the approximate problem to the exact problem. Besides, there is very few results, up to our knowledge, concerning the speed of convergence of approximations of optimal transport problems. The study of these rates of convergence for more general set of test functions and of optimal transport problems is an interesting issue which is left for future research.

The article is organized as follows. Some preliminaries, including the Tchakaloff theorem, are recalled in Section 2. In Section 3, we introduce the approximate MCOT problem and prove under suitable assumptions that one of its minimizers reads as a discrete measure whose number of charged points is estimated depending on the number of moment constraints and on the nature of the optimal transport problem considered. Under additional assumptions, we prove that the MCOT problem converges to the exact optimal transport problem as the number of test functions grows in Section 4. Rates of convergence of the approximate problem to the exact problem depending on the choice of test functions are proved in Section 5. Finally, algorithms which exploits the particular structure of the MCOT problem are proposed in Section 6 and tested on some numerical examples.

2 Preliminaries

2.1 Presentation of the problem and notation

We begin this section by recalling the classical form of the 2-marginal optimal transport (OT) problem, which will be the prototypical problem considered in this paper, and introduce the notation used in the sequel.

Let $d_{x},d_{y}\in\mathbb{N}^{*}$. We assume that $\mathcal{X}\subset\mathbb{R}^{d_{x}}$ (resp. $\mathcal{Y}\subset\mathbb{R}^{d_{y}}$) is a $G_{\delta}$-set, i.e. a countable intersection of open sets. This property ensures by Alexandroff’s lemma (see e.g. [4], p. 88) that $\mathcal{X}$ (resp. $\mathcal{Y}$) is a Polish space for a metric which is equivalent to the original one on $\mathbb{R}^{d_{x}}$ (resp. $\mathbb{R}^{d_{y}}$). In particular, $\mathcal{X}$ can either be a closed or an open set of $\mathbb{R}^{d_{x}}$.

Let $\mu\in\mathcal{P}(\mathcal{X})$ and $\nu\in\mathcal{P}(\mathcal{Y})$ and let us define

[TABLE]

the set of probability couplings between $\mu$ and $\nu$. We consider a non-negative cost function $c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}\cup\{+\infty\}$, which we assume to be lower semi-continuous (l.s.c.). The Kantorovich optimal transport (OT) problem with the two marginal laws $\mu$ and $\nu$ associated to the cost function $c$ is the following minimization problem:

[TABLE]

Under our assumptions, it is known (see e.g. Theorem 1.7 in [31]) that there exists $\pi^{*}\in\Pi(\mu,\nu)$ such that $I=\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\pi^{*}(x,y)$. Problem (2.1.1) will be referred hereafter as the *exact *OT problem, with respect to the *approximate *problem which we define hereafter.

The aim of this paper is to study a relaxation of Problem (2.1.1) with a view to the design of numerical schemes to approximate the exact OT problem. More precisely, the *approximate *problem considered in this paper consists in relaxing the marginal constraints into a *finite number of *moments constraints. Let $M,N\in\mathbb{N}^{*}$ and $(\phi_{m})_{1\leq m\leq M}\subset L^{1}(\mathcal{X},\mu;\mathbb{R})$ (respectively $(\psi_{n})_{1\leq n\leq N}\subset L^{1}(\mathcal{Y},\nu;\mathbb{R})$) measurable real-valued functions that are integrable with respect to $\mu$ (resp. $\nu$). The functions $(\phi_{m})_{1\leq m\leq M}$ and $(\psi_{n})_{1\leq n\leq N}$ will be called *test functions *hereafter. We define for such families of functions

[TABLE]

which is the set of probability measures on $\mathcal{X}\times\mathcal{Y}$ that have the same moments as $\mu$ and $\nu$ for the test functions. We are then interested in the moment constrained optimal transport (MCOT) problem, which we defined as the following minimization problem :

[TABLE]

Since $\Pi(\mu,\nu)\subset\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq M},(\psi_{n})_{1\leq n\leq N})$, we clearly have $I^{M,N}\leq I$. In this paper, we are interested in the following question.

• •

Is the MCOT problem attained by some probability measure $\pi^{*}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq M},(\psi_{n})_{1\leq n\leq N})$?

• •

Under which assumptions does it hold: $I^{M,N}\to_{M,N\to+\infty}I$? Can the speed of convergence be estimated?

For simplicity, we will assume that $M=N$ in the whole paper and we will denote for $1\leq m,n\leq N$:

[TABLE]

For all $x\in\mathcal{X}$ (respectively for all $y\in\mathcal{Y}$), we define $\phi(x):=(\phi_{1}(x),...,\phi_{N}(x))\in\mathbb{R}^{N}$ (respectively $\psi(y):=(\psi_{1}(y),...,\psi_{N}(y))\in\mathbb{R}^{N}$) and $\Phi(x):=(1,\phi(x))\in\mathbb{R}^{N+1}$ (respectively $\Psi(y):=(1,\psi(y))\in\mathbb{R}^{N+1}$).

2.2 Tchakaloff’s theorem

In this section, we present a corollary of the Tchakaloff theorem which is the backbone of our results concerning the existence of a minimizer to the MCOT problem. A general version of the Tchakaloff theorem has been proved by Bayer and Teichmann [5] and Bisgaard [12]. The next proposition is rather an immediate consequence of Tchakaloff’s theorem, see Corollary 2 in [5]. We recall first that

Proposition 2.1**.**

Let $\pi$ be a measure on $\mathbb{R}^{d}$ concentrated on a Borel set $A\in\mathcal{F}$, i.e. $\pi(\mathbb{R}^{d}\setminus A)=0$. Let $N_{0}\in\mathbb{N}^{*}$ and $\Lambda:\mathbb{R}^{d}\to\mathbb{R}^{N_{0}}$ a Borel measurable map. Assume that the first moments of $\Lambda\#\pi$ exist, i.e.

[TABLE]

where $\|\cdot\|$ denotes the Euclidean norm of $\mathbb{R}^{N_{0}}$. Then, there exist an integer $1\leq K\leq N_{0}$, points $z_{1},...,z_{K}\in A$ and weights $p_{1},...,p_{K}>0$ such that

[TABLE]

where $\Lambda_{i}$ denotes the $i$-th component of $\Lambda$.

We recall here that $\Lambda\#\pi$ is the push-forward of $\pi$ through $\Lambda$, and is defined as $\Lambda\#\pi(A)=\pi(\Lambda^{-1}(A))$ for any Borel set $A\subset\mathbb{R}^{N_{0}}$. Let us note here that even if $\pi$ is a probability measure, we may have $\sum_{k=1}^{K}p_{k}\not=1$. In the sequel, we will apply this proposition to functions $\Lambda$ such that the constant $1$ is a coordinate of $\Lambda$, which will ensure $\sum_{k=1}^{K}p_{k}=1$.

Last, let us remark that the number of points $K$ needed may be significantly different from $N_{0}$. Lemma A.1 gives, for any $N\in\mathbb{N}^{*}$, an example with $N_{0}=2N+1$ and $K=N+1$.

2.3 An admissibility property

A natural requirement for the MCOT Problem (2.1.3) to have a sense is to assume that it has finite value. This is precisely our definition of admissibility.

Definition 2.1** (Admissibility).**

Let $\mu\in\mathcal{P}(\mathcal{X}),\ \nu\in\mathcal{P}(\mathcal{Y})$ and a l.s.c. cost function $c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}\cup\{\infty\}$. Then, a set of test functions $\left((\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N}\right)\in L^{1}(\mathcal{X},\mu;\mathbb{R})^{N}\times L^{1}(\mathcal{Y},\nu;\mathbb{R})^{N}$ is said to be admissible for $(\mu,\nu,c)$ if

[TABLE]

Thanks to Tchakaloff’s theorem, the admissibility can be checked on finite probability measure as stated in the next Lemma.

Lemma 2.2**.**

Let $\mu\in\mathcal{P}(\mathcal{X})$, $\nu\in\mathcal{P}(\mathcal{Y})$ and $c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}\cup\{+\infty\}$ a l.s.c. function. A set $\left((\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N}\right)\in L^{1}(\mathcal{X},\mu;\mathbb{R})^{N}\times L^{1}(\mathcal{Y},\nu;\mathbb{R})^{N}$ is admissible for $(\mu,\nu,c)$ if, and only if, there exist weights $w_{1},\dots,w_{2N+1}\geq 0$ and points $(x_{1},y_{1}),\dots,(x_{2N+1},y_{2N+1})\in\mathcal{X}\times\mathcal{Y}$ such that

[TABLE]

In particular, if $c$ is finite valued (i.e. $c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}$), any set of test functions $\left((\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N}\right)\in L^{1}(\mathcal{X},\mu;\mathbb{R})^{N}\times L^{1}(\mathcal{Y},\nu;\mathbb{R})^{N}$ is admissible for $(\mu,\nu,c)$ in the sense of Definition 2.1.

Proof.

Let $\Lambda:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}^{2N+1}$ be defined by $\Lambda_{m}(x,y)=\phi_{m}(x)$ and $\Lambda_{m+N}(x,y)=\psi_{m}(y)$ for $m\in\{1,\dots,N\}$, $\Lambda_{2N+1}(x,y)=1$. Let $A=\{(x,y)\in\mathcal{X}\times\mathcal{Y}:c(x,y)=+\infty\}$. Since the set of test function is admissible, there exists $\gamma\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ such that $\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\gamma(x,y)<\infty$. In particular, $\gamma(A)=0$. We can thus apply Proposition 2.1, which gives the implication. The reciprocal result is obvious.

Last, when $c$ is finite valued ($A=\emptyset$), any $\gamma\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ satisfies $\gamma(A)=0$ and the claim follows by using again Proposition 2.1. ∎

3 Existence of discrete minimizers for MCOT problems

3.1 General case

When Definition 2.1 is satisfied, in order to have the existence of a minimizer for the MCOT problem, we make two further assumptions.

• •

We assume that the test function are continuous.

• •

We add to the MCOT problem (2.1.3) a moment inequality constraint.

The additional moment constraint will ensure the tightness of a minimizing sequence satisfying the moment equality and inequality constraints, while the continuity of the test functions will ensure that any limit satisfies the moment constraints. Our main result is stated in Theorem 3.1 thereafter.

Theorem 3.1**.**

Let $\mu\in\mathcal{P}(\mathcal{X})$, $\nu\in\mathcal{P}(\mathcal{Y})$ and $c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}\cup\{+\infty\}$ a l.s.c. function. Let $\Sigma_{\mu}\subset\mathcal{X},\Sigma_{\nu}\subset\mathcal{Y}$ be Borel sets such that $\mu(\Sigma_{\mu})=\nu(\Sigma_{\nu})=1$. Let $N\in\mathbb{N}^{*}$ and let $\left((\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N}\right)\in L^{1}(\mathcal{X},\mu;\mathbb{R})^{N}\times L^{1}(\mathcal{Y},\nu;\mathbb{R})^{N}$ be an admissible set of test functions for $(\mu,\nu,c)$ in the sense of Definition 2.1. We assume that

For all $n\in\{1,\dots,N\}$, the functions $\phi_{n}$ and $\psi_{n}$ are continuous. 2. 2.

There exist $\theta_{\mu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ and $\theta_{\nu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ two non-negative non-decreasing continuous functions such that $\displaystyle\theta_{\mu}(r)\mathop{\longrightarrow}_{r\to+\infty}+\infty$ and $\displaystyle\theta_{\nu}(r)\mathop{\longrightarrow}_{r\to+\infty}+\infty$, and such that there exist $C>0$ and $0<s<1$ such that for all $1\leq n\leq N$, and all $(x,y)\in\mathcal{X}\times\mathcal{Y}$,

[TABLE]

For all $A>0$, let us introduce

[TABLE]

Then, there exists $A_{0}>0$ such that for all $A\geq A_{0}$, $I^{N}_{A}$ is finite and is a minimum. Moreover, for all $A\geq A_{0}$, there exists a minimizer $\pi^{N}_{A}$ for Problem (3.1.2) such that $\pi^{N}_{A}=\sum_{k=1}^{K}p_{k}\delta_{x_{k},y_{k}}$, for some $0<K\leq 2N+2$, with $p_{k}\geq 0$, $x_{k}\in\Sigma_{\mu}$ and $y_{k}\in\Sigma_{\nu}$ for all $1\leq k\leq K$.

Remark 3.1**.**

When $I$ defined by (2.1.1) is finite and

[TABLE]

we have for all $A\geq A^{\prime}_{0}$, $I^{N}_{A}\leq I<\infty$.

Remark 3.2**.**

When the functions $\phi_{m}$ and $\psi_{n}$ are bounded continuous (which holds automatically when $\mathcal{X}$ and $\mathcal{Y}$ are compact), Assumption (3.1.1) is obviously satisfied. Besides, when $\mathcal{X}$ and $\mathcal{Y}$ are compact sets, we can then take $C=\max_{1\leq n\leq N}(\max(\|\phi_{n}\|_{\infty},\|\psi_{n}\|_{\infty}))$ and $\theta_{\mu}=\theta_{\nu}=0$, and therefore we get for all $A>0$, $I^{N}_{A}=I^{N}$ with

[TABLE]

Proof of Theorem 3.1.

Let us introduce the function

[TABLE]

and let us denote by $\Lambda_{i}$ the $i^{th}$ component of $\Lambda$ for all $1\leq i\leq 2N+2$. By assumption there exists $\gamma\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ such that $\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\gamma(x,y)<\infty$. We apply Proposition 2.1 with $N_{0}=2N+2$ and get that there exist $K\in\{1,\dots,2N+2\}$, $x_{1},...,x_{K}\in\mathcal{X}$, $y_{1},...,y_{K}\in\mathcal{Y}$ and weights $w_{1},...,w_{K}\in\mathbb{R}^{*}_{+}$ such that

[TABLE]

Denoting by $\widetilde{\gamma}:=\sum_{k=1}^{K}w_{k}\delta_{x_{k},y_{k}}$, we have that

[TABLE]

We thus get that, for all $A\geq A_{0}:=\int_{\mathcal{X}\times\mathcal{Y}}(\theta_{\mu}(|x|)+\theta_{\nu}(|y|)){\mathrm{d}}\tilde{\gamma}(x,y)$, $I^{N}_{A}$ is finite, since we have $\tilde{\gamma}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$.

Let us now assume that $A\geq A_{0}$ and let us prove that this infimum is a minimum. Let $(\pi_{l})_{l\in\mathbb{N}}$ be a minimizing sequence for the minimization problem (3.1.2). We first prove the tightness of this sequence. For $M_{1},M_{2}>0$, let us introduce the compact sets

[TABLE]

Then, we have

[TABLE]

which implies the tightness of the sequence $(\pi_{l})_{l\in\mathbb{N}}$. We can thus extract a subsequence that weakly converges. For notational simplicity, we still denote $(\pi_{l})_{l\in\mathbb{N}}$ this subsequence, and there exists $\pi_{\infty}\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ such that $\pi_{l}\xrightharpoonup[l\to\infty]{}\pi_{\infty}$.

By Skorokhod’s representation theorem (see e.g. Theorem 4.30 [26]), there exists a probability space $(\Omega,\mathcal{F},\mathbb{P})$ and random variables $(X_{l},Y_{l})_{l\in\mathbb{N}\cup\{\infty\}}$ on this probability space such that $(X_{l},Y_{l})$ is distributed according to $\pi_{l}$ and $(X_{l},Y_{l})\to(X_{\infty},Y_{\infty})$, $\mathbb{P}$-a.s. From $\sup_{l\in\mathbb{N}}\mathbb{E}[\theta_{\mu}(|X_{l}|)+\theta_{\nu}(|Y_{l}|)]\leq A$ and (3.1.1), we deduce that the families $(\phi_{m}(X_{l}))_{l\in\mathbb{N}}$ and $(\psi_{n}(X_{l}))_{l\in\mathbb{N}}$ uniformly integrable. Therefore, we get from the continuity of $\phi_{m}$ and $\psi_{n}$ that $\mathbb{E}[\phi_{m}(X_{\infty})]=\lim_{l\rightarrow\infty}\mathbb{E}[\phi_{m}(X_{l})]=\mu_{m}$ and $\mathbb{E}[\psi_{n}(X_{\infty})]=\lim_{l\rightarrow\infty}\mathbb{E}[\psi_{n}(X_{l})]=\nu_{n}$ (see e.g. Lemma 4.11 [26]). Fatou’s lemma also gives $\mathbb{E}[\theta_{\mu}(|X_{\infty}|)+\theta_{\nu}(|Y_{\infty}|)]\leq A$, which shows that $\pi_{\infty}$ satisfies the constraints of Problem (3.1.2) and thus

[TABLE]

We now show the converse inequality. Since $c$ is l.s.c, $\liminf_{l\to\infty}c(X_{l},Y_{l})\geq c(X_{\infty},Y_{\infty})$, $\mathbb{P}\mathrm{-a.s.}$ and Fatou’s lemma yields that

[TABLE]

As $(\pi_{l})_{l\in\mathbb{N}}$ is a minimizing sequence associated to Problem (3.1.2), it holds that $\mathbb{E}\left[c(X_{l},Y_{l})\right]\xrightarrow[l\to\infty]{}I^{N}_{A}$. Then, by using (3.1.5), one gets

[TABLE]

Thus, $\pi_{\infty}$ is a minimizer of Problem (3.1.2).

Last, we apply Proposition 2.1 to the measure $\pi_{\infty}$ and the application $\Lambda$ defined in (3.1.3) and get the existence of $\pi^{N}_{A}$. ∎

Example 3.1 below shows that the MCOT problem may not be a minimum if we remove the constraint $\int_{\mathcal{X}\times\mathcal{Y}}(\theta_{\mu}(|x|)+\theta_{\nu}(|y|)){\mathrm{d}}\pi(x,y)\leq A$.

Example 3.1**.**

Let

[TABLE]

where for $r\in\mathbb{R}_{+}$, $\varphi(r)=\mathbf{1}_{0\leq r\leq 1}r+\mathbf{1}_{1<r}e^{1-r}$. Let us consider the problem

[TABLE]

The sequence defined for $l\in\mathbb{N}^{*}$ by $\pi_{l}=\left(1-\frac{1}{l}\right)\delta_{(0,0)}+\frac{1}{l}\delta_{(l,l)}$ is a minimizing sequence since $\int_{\mathbb{R}\times\mathbb{R}}x{\mathrm{d}}\pi_{l}(x,y)=\int_{\mathbb{R}\times\mathbb{R}}y{\mathrm{d}}\pi_{l}(x,y)=1$, $c\geq 0$ and

[TABLE]

Hence, $I=0$. Now, since $\varphi(r)>0$ for $r>0$, the only probability measure $\pi\in\mathcal{P}(\mathbb{R}\times\mathbb{R})$ such that $\int c{\mathrm{d}}\pi=0$ is $\delta_{(0,0)}$. Since this probability measure does not satisfy the constraints ($\int_{\mathbb{R}\times\mathbb{R}}x{\mathrm{d}}\delta_{(0,0)}(x,y)=\int_{\mathbb{R}\times\mathbb{R}}y{\mathrm{d}}\delta_{(0,0)}(x,y)=0$), this shows that $I$ is not a minimum.

Let us also note here that the test functions $(\phi_{m})_{1\leq m\leq N}$ and $(\psi_{n})_{1\leq n\leq N}$ are needed to be continuous to guarantee the existence of a minimum in Theorem 3.1 as Example 3.2 shows.

Example 3.2**.**

Let $\mathcal{X}=\mathcal{Y}=[0,1]$, $\,{\mathrm{d}}\nu(x)=\left(\frac{1}{2}\mathbf{1}_{(0,\frac{1}{2})}(x)+\frac{3}{2}\mathbf{1}_{(\frac{1}{2},1)}(x)\right)\,{\mathrm{d}}x$, $\,{\mathrm{d}}\mu(x)=\,{\mathrm{d}}x$ and $c(x,y)=(y-x)^{2}$. Let $N=4$, $\phi_{1}=\mathbf{1}_{[0,\frac{1}{4}]}$, $\phi_{m}=\mathbf{1}_{(\frac{m-1}{4},\frac{m}{4}]}$ for $2\leq m\leq 4$ and $\psi_{m}=\phi_{m}$ for $1\leq m\leq 4$, so that

[TABLE]

For $l\in\mathbb{N}$, $l>4$, let

[TABLE]

For all $l>4$, $\gamma_{l}$ satisfies the constraints of the MCOT problem, and

[TABLE]

Thus, the infimum value of the associated MCOT problem is [math]. Now, let $\pi\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ be such that $\int c{\mathrm{d}}\pi=0$. We have $\pi(\{(x,y)\in\mathcal{X}\times\mathcal{Y}:y=x\})=1$ and thus

[TABLE]

Therefore, we cannot have the left hand side equal to $\mu_{m}$ and the right hand side equal to $\nu_{m}$, which shows that there does not exist any minimizer to the MCOT problem.

3.2 Compactly supported test functions

An alternative statement of Theorem 3.1 that avoids to impose the constraint $\int_{\mathcal{X}\times\mathcal{Y}}(\theta_{\mu}(|x|)+\theta_{\nu}(|y|)){\mathrm{d}}\pi(x,y)\leq A$ can be obtained under stronger assumptions on the test functions and the cost. In all Subsection 3.2, we consider the case

[TABLE]

for some $d_{x},d_{y}\in\mathbb{N}^{*}$, and assume that the cost $c$ is continuous and satisfies:

[TABLE]

This condition is satisfied for example when $d_{x}=d_{y}$ and $c(x,y)=H(|x-y|)$, with $H$ continuous satisfying $H(0)=0$ and $H(r)\underset{r\to+\infty}{\to}+\infty$. We assume also that the test functions $\phi_{m},\psi_{n}$, $1\leq m,n\leq N$ are continuous with compact support, and define their compact support as follows

[TABLE]

Let $M=\max_{x,y\in\mathcal{S}_{\mathcal{X}}\times\mathcal{S}_{\mathcal{Y}}}c(x,y)$ and let us define

[TABLE]

together with

[TABLE]

It can be easily seen that $\tilde{\mathcal{S}}_{\mathcal{X}}$ (resp. $\tilde{\mathcal{S}}_{\mathcal{Y}}$) is a compact set that contains $\mathcal{S}_{\mathcal{X}}$ (resp. ${\mathcal{S}}_{\mathcal{Y}}$), and thus the set $\mathcal{K}$ is compact. Then, from (3.2.2), we take an arbitrary point $(\bar{x},\bar{y})\notin\mathcal{K}$ such that $c(\bar{x},\bar{y})=0$, and we define

[TABLE]

Lemma 3.2**.**

Let $K\in\mathbb{N}^{*}$, and for all $1\leq k\leq K$, $x_{k}\in\mathcal{X}$, $y_{k}\in\mathcal{Y}$, $p_{k}\geq 0$ such that $\sum_{k=1}^{K}p_{k}=1$. If $\gamma=\sum_{k=1}^{K}p_{k}\delta_{x_{k},y_{k}}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ then there exist $K$ points $(\tilde{x}_{k},\tilde{y}_{k})\in\bar{\mathcal{K}}$ for $1\leq k\leq K$ such that the discrete probability measure $\tilde{\gamma}=\sum_{k=1}^{K}p_{k}\delta_{\tilde{x}_{k},\tilde{y}_{k}}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ and

[TABLE]

Proof.

Consider a measure $\gamma=\sum_{k=1}^{K}p_{k}\delta_{x_{k},y_{k}}$ satisfying the assumptions of Lemma 3.2. We construct $\tilde{\gamma}=\sum_{k=1}^{K}p_{k}\delta_{\tilde{x}_{k},\tilde{y}_{k}}$ using the following procedure.

• Case 1:

If $(x_{k},y_{k})\in\mathcal{K}$, then we define $(\tilde{x}_{k},\tilde{y}_{k})=(x_{k},y_{k})$.

• Case 2:

If $x_{k}\notin{\mathcal{S}}_{\mathcal{X}}$ and $y_{k}\notin{\mathcal{S}}_{\mathcal{Y}}$, then we define $(\tilde{x}_{k},\tilde{y}_{k})=(\bar{x},\bar{y})$.

• Case 3:

Let us suppose $x_{k}\in\mathcal{S}_{\mathcal{X}}$ and $y_{k}\notin\tilde{\mathcal{S}}_{\mathcal{Y}}$ (the case $y_{k}\in\mathcal{S}_{\mathcal{Y}}$ and $x_{k}\not\in\tilde{\mathcal{S}}_{\mathcal{X}}$ is treated symmetrically). By definition of $\tilde{\mathcal{S}}_{\mathcal{Y}}$, it holds that

[TABLE]

In particular, we have $c(x_{k},y_{k})>M+1$. Let $y^{*}\in\mathcal{S}_{\mathcal{Y}}$. Then,

[TABLE]

Let $y_{\lambda}:=(1-\lambda)y^{*}+\lambda y_{k}$ for $\lambda\in[0,1]$. As $c$ is continuous, there exists $\lambda^{*}$ such that $c(x_{k},y_{\lambda^{*}})=\frac{2M+1}{2}$. Then, $y_{\lambda^{*}}\notin\mathcal{S}_{\mathcal{Y}}$ because $\frac{2M+1}{2}>M$, and $y_{\lambda^{*}}\in\widetilde{\mathcal{S}}_{\mathcal{Y}}$. Then, we define $(\tilde{x}_{k},\tilde{y}_{k})=(x_{k},y_{\lambda^{*}})$.

This construction preserves the points in the supports $\mathcal{S}_{\mathcal{X}}$ and $\mathcal{S}_{\mathcal{Y}}$, and the points outside the supports are replaced by other points outside the supports. Thus, we have

[TABLE]

and the moment constraints are satisfied by $\tilde{\gamma}$. In addition, it is clear that the cost does not change in Case 1 and is lowered in Cases 2 and 3. ∎

Proposition 3.3**.**

Let us assume that $\mathcal{X}=\mathbb{R}^{d_{x}}$, $\mathcal{Y}=\mathbb{R}^{d_{y}}$ and $c:\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{y}}\rightarrow\mathbb{R}_{+}$ is continuous and satisfies (3.2.1), (3.2.2).

Let us assume that for all $1\leq m,n\leq N$, $\phi_{m}$ and $\psi_{n}$ are compactly supported real-valued continuous functions defined on $\mathbb{R}^{d}$. Then, there exists at least one minimizer to the minimization problem

[TABLE]

Moreover, there exists $K\in\mathbb{N}$ such that $K\leq 2N+2$, and for all $1\leq k\leq K$, $(x_{k},y_{k})\in\bar{\mathcal{K}}$, $p_{k}\geq 0$ such that $\sum_{k=1}^{K}p_{k}=1$ such that $\widetilde{\pi}:=\sum_{k=1}^{K}p_{k}\delta_{x_{k},y_{k}}$ is a minimum.

Proof.

Let us consider a minimizing sequence $(\pi_{l})_{l\in\mathbb{N}}$ for Problem (3.2.6). For all $l\in\mathbb{N}$, we will denote by $\gamma_{l}$ a finite discrete measure which has the same cost and same moments than $\pi_{l}$, with at most $2N+2$ points, which exists thanks to Proposition 2.1, and the fact that the test functions are compactly supported. Then, using Lemma 3.2, for all $l\in\mathbb{N}$, one can define a measure $\tilde{\gamma}_{l}$ which satisfies the moment constraints, has a support contained in the set $\bar{\mathcal{K}}$ defined in (3.2.5), and such that,

[TABLE]

Thus, $(\tilde{\gamma}_{l})_{l\in\mathbb{N}}$ is a minimizing sequence. Besides, $(\tilde{\gamma}_{l})_{l\in\mathbb{N}}$ is tight since the support of $\widetilde{\gamma}_{l}$ is included in the compact set $\bar{\mathcal{K}}$ for all $l\in\mathbb{N}$. Then, following the same lines as in the proof of Theorem 3.1, one can extract a weakly converging subsequence, the cost of the limit of which is equal to $I^{N}$. The fact that there exists a finite discrete measure charging at most $K\leq 2N+2$ points can be deduced from Proposition 2.1, following the same lines as in the proof of Theorem 3.1. ∎

3.3 Multimarginal and martingale OT problem

In this section, two important extensions of the previous problem are introduced, the multimarginal problem and the martingale problem. Alike Problem (3.1.2), several formulations and refinements can be established. We only keep here the more general ones for conciseness.

3.3.1 Multimarginal problem

The propositions introduced until now for two marginal laws can be extended to an arbitrary (finite) number of marginal laws. The proof can be straightforwardly adapted from the one of Theorem 3.1. For all $1\leq i\leq M$, we consider $\mathcal{X}_{i}=\mathbb{R}^{d_{i}}$ with $d_{i}\in\mathbb{N}^{*}$ or more generally a $G_{\delta}$-set $\mathcal{X}_{i}\subset\mathbb{R}^{d_{i}}$. We consider $M$ probability measures $\mu_{1}\in\mathcal{P}(\mathcal{X}_{1}),...,\mu_{M}\in\mathcal{P}(\mathcal{X}_{M})$ and an l.s.c. cost function $c:\mathcal{X}_{1}\times...\times\mathcal{X}_{M}\to\mathbb{R}_{+}\cup\{\infty\}$. We consider the following multimarginal optimal transport problem

[TABLE]

where $\Pi(\mu_{1},...,\mu_{M})=\{\pi\in\mathcal{P}(\mathcal{X}_{1}\times...\times\mathcal{X}_{M})\,\text{s.t.}\forall 1\leq i\leq M,\,\int_{\mathcal{X}_{i}}{\mathrm{d}}\pi={\mathrm{d}}\mu_{i}\}$.

In order to build the moments constrained optimal transport problem, we introduce, for each $i$, $N_{i}\in\mathbb{N}^{*}$ test functions $(\phi^{i}_{n})_{1\leq n\leq N_{i}}\in L^{1}(\mathcal{X}_{i},\mu_{i};\mathbb{R})^{N_{i}}$. We say that this set of test functions is admissible for $(\mu_{1},\dots,\mu_{M},c)$ if there exists $\gamma\in\mathcal{P}(\mathcal{X}_{1}\times...\times\mathcal{X}_{M})$ such that

[TABLE]

and $\int_{\mathcal{X}_{1}\times...\times\mathcal{X}_{M}}c(x_{1},\dots,x_{M}){\mathrm{d}}\gamma(x_{1},\dots,x_{M})<\infty$. We can now state the analogous of Theorem 3.1 for the multimarginal case.

Proposition 3.4**.**

For $i\in\{1,\dots,M\}$, let $\mu_{i}\in\mathcal{P}(\mathcal{X}_{i})$ and $\Sigma_{\mu_{i}}\subset\mathcal{X}_{i}$ a Borel set such that $\mu_{i}(\Sigma_{\mu_{i}})=1$. We assume that $c:\mathcal{X}_{1}\times...\times\mathcal{X}_{M}\to\mathbb{R}^{+}\cup\{\infty\}$ is a l.s.c. cost function, and that the set of test functions $\phi^{i}_{n}\in L^{1}(\mathcal{X}_{i},\mu_{i};\mathbb{R})$ for $i\in\{1,\dots,M\}$ and $n\in\{1,\dots,N_{i}\}$ is admissible for $(\mu_{1},\dots,\mu_{M},c)$. We make the following assumptions.

For all $i$ and $n$, the function $\phi^{i}_{n}$ is continuous. 2. 2.

For all $i$, there exists $\theta_{i}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ a non-decreasing continuous function such that $\theta_{i}(r)\underset{r\to+\infty}{\to}+\infty$ and such that there exist $C>0$ and $0<s<1$ such that for all $1\leq n\leq N_{i}$, we have

[TABLE]

We note $\mathbf{N}=(N_{1},\dots,N_{M})$, $\mathcal{X}=\mathcal{X}_{1}\times\dots\times\mathcal{X}_{M}$ and consider the following problem

[TABLE]

Then, there exists $A_{0}>0$ such that for all $A\geq A_{0}$, $I^{\mathbf{N}}_{A}$ is finite and is a minimum. Moreover, for all $A\geq A_{0}$, there exists a minimizer $\pi^{\mathbf{N}}_{A}$ for the problem (3.1.2) such that $\pi^{\mathbf{N}}_{A}=\sum_{k=1}^{K}p_{k}\delta_{x^{k}_{1},\dots,x^{k}_{M}}$, for some $0<K\leq\sum_{i=1}^{M}N_{i}+2$, with $p_{k}\geq 0$ and $x^{k}_{i}\in\Sigma_{\mu_{i}}$ for all $1\leq i\leq M$ and $1\leq k\leq K$.

An interesting point to remark in Proposition 3.4 is that the number of weighted points of the discrete measure $\pi^{\mathbf{N}}_{A}$ is linear with respect to the number of moment constraints. In particular, if we take the same number of moments $N_{i}=N$ for each marginal, the number of weighted points is equal to $2+MN$ and thus grows linearly with respect to $M$. Since each points has $dM$ coordinates, the dimension of the discrete measure is in $O(M^{2})$. For this reason, the development of algorithms for minimizing $\pi^{\mathbf{N}}_{A}$ by using finite discrete measures may be a way to avoid the curse of dimensionality when $M$ is getting large.

We make here a specific focus on the multimarginal optimal transport problem which arises in quantum chemistry applications [33, 16]. In this particular case, the multi-marginal optimal transport of interest reads as (3.3.1), with $\mathcal{X}_{1}=\cdots\mathcal{X}_{M}=\mathbb{R}^{3}$, $N_{1}=\cdots=N_{M}=N$ for some $N\in\mathbb{N}^{*}$, $\mu_{1}=\cdots=\mu_{M}=\mu$ for some $\mu\in\mathcal{P}(\mathbb{R}^{3})$ and $c$ is given by the Coulomb cost

[TABLE]

The integer $M$ represents here the number of electrons in the system of interest. The inherent symmetries of the system yield interesting consequences on the MCOT problem (3.3.3), which are summarized in the following proposition.

Proposition 3.5**.**

Let $M\in\mathbb{N}^{*}$, $N\in\mathbb{N}^{*}$, $\mu\in\mathcal{P}(\mathcal{X})$ and $\Sigma_{\mu}\subset\mathcal{X}$ a Borel set such that $\mu(\Sigma_{\mu})=1$. We assume that $c:\mathcal{X}^{M}\to\mathbb{R}^{+}\cup\{\infty\}$ is a symmetric l.s.c. cost function. More precisely, we denote by $\mathcal{S}_{M}$ the set of permutations of the set $\{1,\cdots,M\}$ and assume that

[TABLE]

For all $1\leq n\leq N$, let $\phi_{n}\in L^{1}(\mathcal{X},\mu;\mathbb{R})$. We define $\phi_{n}^{i}:=\phi_{n}$ for all $1\leq i\leq M$ and assume the set of test functions $\phi_{n}^{i}$ for $n\in\{1,\dots,N\}$ and $i\in\{1,\cdots,M\}$ is admissible for $(\mu,\dots,\mu,c)$. We make the following assumptions.

For all $n$, the function $\phi_{n}$ is continuous. 2. 2.

There exists $\theta:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ a non-decreasing continuous function such that $\theta(r)\underset{r\to+\infty}{\mathop{\longrightarrow}}+\infty$ and such that there exist $C>0$ and $0<s<1$ such that for all $1\leq n\leq N$, we have

[TABLE]

We consider the following problem

[TABLE]

Then, it holds that

[TABLE]

and there exists $A_{0}>0$ such that for all $A\geq A_{0}$, $I^{N}_{A}$ is finite and is a minimum. Moreover, for all $A\geq A_{0}$, there exists a minimizer $\pi^{N}_{A}$ for the problem (3.3.6) such that $\pi^{N}_{A}=\sum_{k=1}^{K}p_{k}\delta_{x^{k}_{1},\dots,x^{k}_{M}}$, for some $0<K\leq N+2$, with $p_{k}\geq 0$ and $x^{k}_{i}\in\Sigma_{\mu}$ for all $1\leq i\leq M$ and $1\leq k\leq K$. Besides, the symmetric measure

[TABLE]

is a minimizer to (3.3.5) and (3.3.6).

Proof.

It is obvious that the right hand side of (3.3.6) is smaller than the right hand side of (3.3.5). By using the same arguments as in the proof of Theorem 3.1, there exists $A_{0}>0$ such that for all $A\geq A_{0}$ the infimum of (3.3.6) is finite, is a minimum that is attained by some discrete probability measure $\pi^{N}_{A}=\sum_{k=1}^{K}p_{k}\delta_{x^{k}_{1},\dots,x^{k}_{M}}$, for some $0<K\leq N+2$ with $x^{k}_{i}\in\Sigma_{\mu}$ for all $1\leq i\leq M$ and $1\leq k\leq K$. Then, since $c$ is symmetric and the set of constraints is also symmetric, we get that $\pi_{\rm sym,A}^{N}$ also realizes the minimum. Besides, it satisfies $\int_{\mathcal{X}^{M}}\phi_{n}(x_{i}){\mathrm{d}}\pi_{\rm sym,A}^{N}(x_{1},...,x_{M})=\int_{\mathcal{X}}\phi_{n}(x){\mathrm{d}}\mu(x)$ for all $n,i$, which shows that it is also the minimizer of (3.3.5). ∎

Proposition 3.5 is particularly interesting for the design of numerical schemes for the resolution of multimarginal optimal transport problems with Coulomb cost arising in quantum chemistry applications. Indeed, the latter read as (3.3.5) and the number of charged points of the minimizer $\pi^{N}_{A}$ of (3.3.6) only scales at most like $N+2$, and that the dimension of the optimal discrete measure is in $dM(N+2)$. This result states that such formulation of the multimarginal optimal transport problem does not suffer from the curse of dimensionality. Let us mention that this result is close in spirit to the recent work [20], where multimarginal optimal transport problems with Coulomb cost are studied on finite state spaces.

3.3.2 Martingale OT problem

In this paragraph, we assume $\mathcal{X}=\mathcal{Y}=\mathbb{R}^{d}$ with $d\in\mathbb{N}^{*}$, and consider two probability measures $\mu,\nu\in\mathcal{P}(\mathbb{R}^{d})$ such that

[TABLE]

and $\mu$ is lower than $\nu$ for the convex order, i.e.

[TABLE]

for any convex function $\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ non-negative or integrable with respect to $\mu$ and $\nu$. This latter condition is equivalent, by Strassen’s theorem [35], to the existence of a martingale coupling between $\mu$ and $\nu$, i.e.

[TABLE]

The original martingale optimal transport consists then in the minimization problem

[TABLE]

with $c:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}\cup\{\infty\}$ being a l.s.c. cost function. This problem has recently got a great attention in mathematical finance since the work of Beiglböck et al. [6], because it is related to the calculation of model-independent option price bounds on an arbitrage free market.

We consider a set of test functions $(\phi_{m})_{1\leq m\leq N}\in L^{1}({\mathbb{R}^{d}},\mu;\mathbb{R})^{N}$ and $(\psi_{n})_{1\leq n\leq N}\in L^{1}({\mathbb{R}^{d}},\nu;\mathbb{R})^{N}$, and the following problem:

[TABLE]

Suppose for simplicity that there exist some minimizer to this problem $\pi^{*}$. Then, by using Theorem 5.1 in Beiglböck and Nutz [8] that is an extension of Tchakaloff’s theorem to the martingale case, there exists a probability measure $\tilde{\pi}^{*}$ weighting at most $(d+2N+2)^{2}$ points such that $\tilde{\pi}^{*}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$,

[TABLE]

and

[TABLE]

However, the minimization problem for $I^{N}$ still has the martingale constraints. To get a problem that is similar to the MCOT, we then relax in addition the martingale constraint. This constraint is equivalent to have

[TABLE]

for all bounded measurable functions $f:\mathbb{R}^{d}\to\mathbb{R}$, and also for all function $f:\mathbb{R}^{d}\to\mathbb{R}$ such that $\int_{\mathbb{R}^{d}}|xf(x)|{\mathrm{d}}\mu(x)<\infty$. Then, it is natural to consider $N^{\prime}$ test functions $\chi_{l}:\mathbb{R}^{d}\to\mathbb{R}$, $1\leq l\leq N^{\prime}$, such that

[TABLE]

and then to consider the following minimization problem

[TABLE]

We will say that the test functions $(\phi_{m})_{1\leq m\leq N}$, $(\psi_{n})_{1\leq n\leq N}$ and $(\chi_{l})_{1\leq l\leq N^{\prime}}$ are admissible for the martingale problem of $(\mu,\nu,c)$ if $I^{N,N^{\prime}}<\infty$. Similarly to Theorem 3.1, we get the following result.

Proposition 3.6**.**

Let $\mu\in\mathcal{P}(\mathbb{R}^{d})$, $\nu\in\mathcal{P}(\mathbb{R}^{d})$ and $c:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}\cup\{+\infty\}$ a l.s.c. function. Let $\Sigma_{\mu},\Sigma_{\nu}\subset\mathbb{R}^{d}$ be Borel sets such that $\mu(\Sigma_{\mu})=\nu(\Sigma_{\nu})=1$. Let $N\in\mathbb{N}^{*}$ and let $(\phi_{m})_{1\leq m\leq N}\in L^{1}(\mathbb{R}^{d},\mu;\mathbb{R})^{N}$, $(\psi_{n})_{1\leq n\leq N}\in L^{1}(\mathbb{R}^{d},\nu;\mathbb{R})^{N}$ and $(\chi_{l})_{1\leq l\leq N^{\prime}}$ satisfying (3.3.10) be an admissible set of test functions for the martingale problem of $(\mu,\nu,c)$. We make the following assumptions.

For all $n\in\{1,\dots,N\}$, $l\in\{1,\dots,N^{\prime}\}$, the functions $\phi_{n}$, $\psi_{n}$ and $\chi_{l}$ are continuous. 2. 2.

There exist $\theta_{\mu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ and $\theta_{\nu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ two non-negative non-decreasing continuous functions such that $\displaystyle\theta_{\mu}(r)\mathop{\longrightarrow}_{r\to+\infty}+\infty$ and $\displaystyle\theta_{\nu}(r)\mathop{\longrightarrow}_{r\to+\infty}+\infty$, and such that there exist $C>0$ and $0<s<1$ such that for all $1\leq n\leq N$, $1\leq l\leq N^{\prime}$, and all $(x,y)\in{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}$,

[TABLE]

For all $A>0$, let us introduce

[TABLE]

Then, there exists $A_{0}>0$ such that for all $A\geq A_{0}$, $I^{N,N^{\prime}}_{A}$ is finite and is a minimum. Moreover, for all $A\geq A_{0}$, there exists a minimizer $\pi^{N,N^{\prime}}_{A}$ for Problem (3.3.13) such that $\pi^{N,N^{\prime}}_{A}=\sum_{k=1}^{K}p_{k}\delta_{x_{k},y_{k}}$, for some $0<K\leq 2N+N^{\prime}+2$, with $p_{k}\geq 0$, $x_{k}\in\Sigma_{\mu}$ and $y_{k}\in\Sigma_{\nu}$ for all $1\leq k\leq K$.

The proof of Proposition 3.6 follows exactly the same line as the proof of Theorem 3.1, since the relaxation of the martingale moment constraints only brings new moment constraints. Let us stress that the minimizer $\pi^{N,N^{\prime}}_{A}$ does not satisfy in general the martingale constraint. Also, we do not impose in Proposition 3.6 to have (3.3.8), i.e. $\mu$ smaller than $\nu$ for the convex order. In fact, the admissibility condition already ensures that $I^{N,N^{\prime}}<\infty$ and thus, by using Proposition 2.1 that $I^{N,N^{\prime}}_{A}<\infty$ for $A$ large enough. Nonetheless, if we assume in addition that $\mu$ smaller than $\nu$ for the convex order and that $I$, the infimum of Problem 3.3.9, is finite, then we have $I^{N,N^{\prime}}_{A}<\infty$ and $I^{N,N^{\prime}}_{A}\leq I$ for any $A\geq\int_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}(\theta_{\mu}(|x|)+\theta_{\nu}(|y|)){\mathrm{d}}\pi^{1}(x,y)$, where $\pi^{1}\in\Pi(\mu,\nu)$ is such that $\int_{\mathbb{R}^{d}}y{\mathrm{d}}\pi^{1}(x,y)=x$ and $\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}c(x,y){\mathrm{d}}\pi^{1}(x,y)\leq I+1$.

4 Convergence of the MCOT problem towards the OT problem

The aim of this section is to prove that when the number of test functions $N\to+\infty$, the minimizer of the MCOT problem converges towards a minimizer of the OT problem, under appropriate assumptions and up to the extraction of a subsequence.

4.1 Continuous test functions on unbounded domains

Let us consider two sequences of continuous real-valued test functions $(\phi_{m})_{m\in\mathbb{N}^{*}}$ and $(\psi_{n})_{n\in\mathbb{N}^{*}}$ defined on $\mathcal{X}$ (resp. $\mathcal{Y}$) and make the following assumptions.

Let us first assume that there exist continuous non-decreasing non-negative functions $\theta_{\mu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ and $\theta_{\nu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ such that

[TABLE]

and

[TABLE]

In the sequel, we set

[TABLE]

We assume moreover that there exist $(s^{\mu}_{m})_{m\in\mathbb{N}^{*}},(s^{\nu}_{n})_{n\in\mathbb{N}^{*}}\in(0,1)^{\mathbb{N}^{*}}$ and $(C^{\mu}_{m})_{m\in\mathbb{N}^{*}},(C^{\nu}_{n})_{n\in\mathbb{N}^{*}}\in(\mathbb{R}_{+}^{*})^{\mathbb{N}^{*}}$ such that

[TABLE]

Last, we assume that the probability measures $\mu$ and $\nu$ are fully characterized by their moments:

[TABLE]

We consider the optimal cost for the OT problem (2.1.1) that we restate here for convenience

[TABLE]

and for all $N\in\mathbb{N}^{*}$, we define the $N^{th}$ MCOT problem,

[TABLE]

Theorem 4.1**.**

Let $\mu\in\mathcal{P}(\mathcal{X})$ and $\nu\in\mathcal{P}(\mathcal{Y})$ satisfying (4.1.2) for some continuous non-decreasing functions $\theta_{\mu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ and $\theta_{\nu}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ satisfying (4.1.1). Let $c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}\cup\{+\infty\}$ a l.s.c. function. Let $(\phi_{m})_{m\in\mathbb{N}^{*}}\subset L^{1}(\mathcal{X},\mu;\mathbb{R})$ and $(\psi_{n})_{n\in\mathbb{N}^{*}}\subset L^{1}(\mathcal{Y},\nu;\mathbb{R})$ be continuous functions satisfying (4.1.4), (4.1.5), (4.1.6) and (4.1.7). Let us finally assume that $I$, defined by (4.1.8) is finite.

Then, for all $N\in\mathbb{N}^{*}$, there exist at least one minimizer for Problem (4.1.9) and

[TABLE]

Besides, from every sequence $(\pi^{N})_{N\in\mathbb{N}^{*}}$ such that for all $N$, $\pi^{N}\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ is a minimizer for (4.1.9), one can extract a subsequence which converges towards a minimizer $\pi^{\infty}\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ to problem (4.1.8).

Proof.

From Theorem 3.1 and Remark 3.1, We know that there exists at least one minimizer $\pi^{N}\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ to (4.1.9). Since we have

[TABLE]

and (4.1.1), we get that the sequence $(\pi^{N})_{N\in\mathbb{N}^{*}}$ is tight. Thus, up to the extraction of a subsequence, still denoted $(\pi^{N})_{N\in\mathbb{N}^{*}}$ for the sake of simplicity, there exists a measure $\pi^{\infty}\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ such that $\pi^{N}\xrightharpoonup[N\to\infty]{}\pi^{\infty}$. With the same argument as in the proof of Theorem 3.1, we get that for all $m,n\in\mathbb{N}^{*}$,

[TABLE]

Then, Properties (4.1.6) and (4.1.7) give $\pi^{\infty}\in\Pi(\mu,\nu)$. Therefore,

[TABLE]

We now establish that $\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\pi^{\infty}(x,y)\leq I$ which concludes the proof. The Skorokhod representation theorem states that there exist a space $(\Omega,\mathcal{F},\mathbb{P})$ and random variables $(X_{N},Y_{N})_{N\in\mathbb{N}\cup\{\infty\}}$ such that $(X_{N},Y_{N})$ is distributed according to $\pi^{N}$ and $(X_{N},Y_{N})\to(\tilde{X},\tilde{Y})$ $\mathbb{P}$-a.s.

Furthermore, as $c$ is l.s.c,

[TABLE]

and by Fatou’s lemma we get $\liminf_{N\to\infty}\mathbb{E}\left(c({X}_{N},{Y}_{N})\right)\geq\mathbb{E}\left(c(X_{\infty},Y_{\infty})\right)$, i.e.

[TABLE]

Furthermore, note that $(I^{N})_{N\in\mathbb{N}}$ is a non-decreasing sequence and that for all $N\in\mathbb{N}^{*}$, $I^{N}\leq I$. Thus, there exists $I^{\infty}\leq I$ such that $I^{N}\xrightarrow[N\to\infty]{}I^{\infty}$. Recall that

[TABLE]

then (4.1.11) implies,

[TABLE]

Using equations (4.1.10) and (4.1.12), one gets

[TABLE]

Thus, as $\pi^{\infty}\in\Pi(\mu,\nu)$, $\pi^{\infty}\in\operatorname*{arg\,min}_{\pi\in\Pi(\mu,\nu)}\left\{\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\pi(x,y)\right\}$, and $I^{N}\xrightarrow[N\to\infty]{}I$. ∎

4.2 Bounded test functions on compact sets

We now assume that $\mathcal{X}$ and $\mathcal{Y}$ are compact subsets of $\mathbb{R}^{d_{x}}$ and $\mathbb{R}^{d_{y}}$. We state a result analogous to Theorem 4.1 that holds without the additional moment constraint and for possibly discontinuous test functions. We consider two sequences of bounded measurable real-valued test functions $(\phi_{m})_{m\in\mathbb{N}^{*}}\subset L^{\infty}(\mathcal{X})$ and $(\psi_{n})_{n\in\mathbb{N}^{*}}\subset L^{\infty}(\mathcal{Y})$ that satisfy

[TABLE]

and

[TABLE]

It is easy then to see that the properties (4.1.6) and (4.1.7) are satisfied for any $\mu\in\mathcal{P}(\mathcal{X})$ and $\nu\in\mathcal{P}(\mathcal{Y})$. For any $N\geq 1$, we consider the following MCOT problem:

[TABLE]

Proposition 4.2**.**

Let us assume that $\mathcal{X}$ and $\mathcal{Y}$ are compact sets and let $\mu\in\mathcal{P}(\mathcal{X})$ and $\nu\in\mathcal{P}(\mathcal{Y})$. Let $(\phi_{m})_{m\in\mathbb{N}^{*}}\subset L^{\infty}(\mathcal{X})$ and $(\psi_{n})_{n\in\mathbb{N}^{*}}\subset L^{\infty}(\mathcal{Y})$ satisfying (4.2.1) and (4.2.2). Let us assume that $I<+\infty$. Then, it holds that $I^{N}\leq I$ and

[TABLE]

Moreover, from every sequence $(\pi^{N})_{N\in\mathbb{N}}$ such that for all $N\in\mathbb{N}^{*}$, $\pi^{N}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ satisfies

[TABLE]

with $\epsilon_{N}\underset{n\to+\infty}{\longrightarrow}0$, one can extract a subsequence which converges towards a measure $\pi^{\infty}\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ which is a minimizer to Problem (4.1.8).

Remark 4.1**.**

From Proposition 2.1, there exists $0\leq K_{N}\leq 2N+2$, $x_{1},\cdots,x_{K_{N}}\in\mathcal{X}$, $y_{1},\cdots,y_{K_{N}}\in\mathcal{Y}$ and $w_{1},\cdots,w_{K_{N}}\geq 0$ such that $\gamma^{N}:=\sum_{k=1}^{K_{N}}w_{k}\delta_{(x_{k},y_{k})}\in\Pi(\mu,\nu;(\phi_{m})_{1\leq m\leq N},(\psi_{n})_{1\leq n\leq N})$ and

[TABLE]

In other words, any sequence $(\pi^{N})_{N\in\mathbb{N}^{*}}$ satisfying the assumptions of Proposition 4.2 can be chosen as a discrete measure charging at most $2N+2$ points.

Proof of Proposition 4.2.

Since $\mathcal{X}$ and $\mathcal{Y}$ are compact, the sequence $(\pi^{N})$ is tight and we can assume, up to the extraction of a subsequence, that it weakly converges to $\pi^{\infty}$. For $N\in\mathbb{N}^{*}\cup\{\infty\}$, we denote the marginal laws of $\pi^{N}$ respectively by ${\mathrm{d}}\mu^{N}(x):=\int_{\mathcal{Y}}{\mathrm{d}}\pi^{N}(x,y)$ and ${\mathrm{d}}\nu^{N}(y):=\int_{\mathcal{X}}{\mathrm{d}}\pi^{N}(x,y)$. For $f\in C^{0}(\mathcal{X})$, it holds that

[TABLE]

Let $\epsilon>0$. Using the density condition (4.2.1), one can find $M\in\mathbb{N}^{*}$ and $\lambda_{1},...,\lambda_{M}\in\mathbb{R}$ such that $\sup_{x\in\mathcal{X}}\left|f(x)-\sum_{i=1}^{M}\lambda_{i}\phi_{i}(x)\right|\leq\epsilon$. Thus,

[TABLE]

and for $K>M$, $\left|\int_{\mathcal{X}}f{\mathrm{d}}\mu^{K}-\sum_{i=1}^{M}\lambda_{i}\int_{\mathcal{X}}\phi_{i}{\mathrm{d}}\mu^{K}\right|\leq\epsilon$, i.e.

[TABLE]

Then, (4.2.6) and (4.2.7) imply that $\left|\int_{\mathcal{X}}f{\mathrm{d}}\mu^{K}-\int_{\mathcal{X}}f{\mathrm{d}}\mu\right|\leq 2\epsilon$, and taking $K\to\infty$ leads to

[TABLE]

As (4.2.8) holds for any $\epsilon>0$, one gets that for any $f\in C^{0}(\mathcal{X})$,

[TABLE]

which yields that $\mu^{\infty}=\mu$. Similarly, it holds that $\nu^{\infty}=\nu$. Therefore, $\pi^{\infty}\in\Pi(\mu,\nu)$ and

[TABLE]

Now, we use the same arguments as in the proof of Proposition 4.1 to get $\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\pi^{\infty}(x,y)\leq I$, which gives the result. ∎

4.3 Convergence for Martingale Optimal Transport problems

In this subsection, we study the convergence of $I^{N,N^{\prime}}_{A}$ defined by (3.3.11) when the number of test functions for the martingale condition $N^{\prime}\rightarrow+\infty$ towards the following minimization problem:

[TABLE]

This convergence is particularly interesting for the practical application in finance: the marginal laws $\mu,\nu$ are in general not observed and market data only provide some moments. For $d=1$, market data give the prices of European put (or call) options that corresponds to $\phi_{m}(x)=(K_{m}-x)^{+}$ and $\psi_{n}(y)=(K^{\prime}_{m}-y)^{+}$. Let us assume for simplicity zero interest rates. Then, by taking $\theta_{\mu}(|x|)=\theta_{\nu}(|x|)=|x|$, we have from the martingale assumption $\int_{\mathcal{X}\times\mathcal{Y}}(|x|+|y|){\mathrm{d}}\pi(x,y)=2S_{0}$, where $S_{0}>0$ is the current price of the underlying asset. Then, a natural choice would be to take $A_{0}=2S_{0}$. Therefore, the convergence stated in Proposition 4.3 gives a way to approximate option price bounds by taking into account that only some moments are known, while the few existing numerical methods for Martingale Optimal Transport in the literature assume that the marginal laws are known [2, 3, 23].

Proposition 4.3**.**

Let $\mu\in\mathcal{P}(\mathbb{R}^{d})$ lower than $\nu\in\mathcal{P}(\mathbb{R}^{d})$ for the convex order and $c:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}\cup\{+\infty\}$ a l.s.c. function. We assume $|x|\leq\theta_{\mu}(|x|)$, $|y|\leq\theta_{\nu}(|y|)$ and suppose $A_{0}<\infty$ with $A_{0}$ defined by (4.1.3). We assume that the test functions $(\chi_{l},l\in\mathbb{N}^{*})$ are bounded and such that for any function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ continuous with compact support, we have

[TABLE]

Let the assumptions of Proposition 3.6 hold for any $N^{\prime}\geq 1$. Then, we have $I_{A_{0}}^{N,N^{\prime}}\underset{N^{\prime}\to+\infty}{\longrightarrow}I^{N,mg}_{A_{0}}<\infty$.

Proof.

Since $A_{0}<\infty$, any martingale coupling between $\mu$ and $\nu$ satisfies the constraints of $I_{A_{0}}^{N,N^{\prime}}$. By using Tchakaloff’s theorem and the fact that $c$ is finite-valued, we get that $I^{N,N^{\prime}}_{A_{0}}$ is finite for any $N^{\prime}$ and is attained by a measure denoted by $\pi^{N^{\prime}}$ according to Proposition 3.6. Similarly, using Tchakaloff’s theorem for the martingale case, Theorem 5.1 [8], we get that $I^{N,mg}_{A_{0}}<\infty$. Note that from the inclusion of the constraints, we clearly have $I^{N,N^{\prime}_{1}}_{A_{0}}\leq I^{N,N^{\prime}_{2}}_{A_{0}}\leq I^{N,mg}_{A_{0}}$ for $N^{\prime}_{1}\leq N^{\prime}_{2}$. We can then repeat the arguments in the proof of Theorem 4.1 to get that $(\pi^{N^{\prime}})$ is tight and any limit $\pi^{\infty}$ of a weakly converging subsequence satisfies $I^{N,mg}_{A_{0}}=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}c(x,y){\mathrm{d}}\pi^{\infty}(x,y)$.

The only thing to prove is that $\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}(y-x)f(x)d\pi^{\infty}(x,y)=0$ for any function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ continuous with compact support. Let $\epsilon>0$. By assumption, there exists $M\in\mathbb{N}^{*}$ and $\lambda_{1},\dots,\lambda_{M}\in\mathbb{R}$ such that $\sup_{x\in\mathbb{R}^{d}}|f(x)-\sum_{l=1}^{M}\lambda_{l}\chi_{l}(x)|\leq\epsilon$. Therefore, for $N^{\prime}\geq M$, we have

[TABLE]

by using the triangle inequality and the fact that $|x|\leq\theta_{\mu}(|x|)$, $|y|\leq\theta_{\nu}(|y|)$. We conclude then easily letting $N^{\prime}\rightarrow\infty$. ∎

Let us mention that we can obtain using similar arguments that $I_{A_{0}}^{N,mg}$ and $I_{A_{0}}^{N,N^{\prime}}$ converge towards (3.3.9) as $N$ and $N^{\prime}$ go to infinity.

5 Rates of convergence for particular sets of test functions

Throughout this section, we assume that

[TABLE]

and for all $N\in\mathbb{N}^{*}$, we define the intervals

[TABLE]

We investigate in this section the rate of convergence of $I^{N}$ defined by (4.2.3) towards $I$ defined by (4.1.8), when the test functions are piecewise constant (resp. piecewise linear) on $T^{N}_{m}$. We obtain, under suitable assumptions a convergence rate of $O(1/N)$ (resp. $O(1/N^{2})$). This shows, as one may expect, the importance of the choice of test functions to approximate the Optimal Transport problem.

5.1 Piecewise constant test functions on compact sets

In this section, we assume that the cost function $c:[0,1]^{2}\to\mathbb{R}_{+}$ is Lipschitz:

[TABLE]

We define, for $\pi\in\mathcal{P}([0,1]^{2})$, $I(\pi)=\int_{\mathcal{X}\times\mathcal{Y}}c(x,y){\mathrm{d}}\pi(x,y)$ and

[TABLE]

We introduce the piecewise constant test functions

[TABLE]

and consider the MCOT problem:

[TABLE]

Then, Proposition 5.1 establishes the rate of convergence of the sequence $(I^{N})_{N\in\mathbb{N}^{*}}$ to $I$ as $N$ increases. problem $N$ increases.

Proposition 5.1**.**

Let $\mu,\nu\in\mathcal{P}([0,1])$ and $c:[0,1]^{2}\to\mathbb{R}_{+}$ a Lipschitz function with Lipschitz constant $K>0$. Then, for all $N\in\mathbb{N}^{*}$,

[TABLE]

Remark 5.1**.**

Let us note that we are not exactly in the framework of Section 4, since the test functions depends on $N$. However, we have

[TABLE]

and thus Proposition 4.2 gives for any $L\in\mathbb{N}^{*}$,

[TABLE]

Before proving Proposition 5.1, we state a result that bounds the distance between an MCOT optimizer and the minimizer of the OT problem (5.1.2). We define that for $p\geq 1$, the $W_{p}$-Wasserstein distance between $\eta_{1},\eta_{2}\in\mathcal{P}(\mathbb{R}^{d})$ as $W^{p}_{p}(\eta_{1},\eta_{2})=\inf_{\pi\in\Pi(\eta_{1},\eta_{2})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x_{1}-x_{2}\|_{p}^{p}{\mathrm{d}}\pi(x_{1},x_{2})$, i.e. we take the $\|\|_{p}$-norm for $W_{p}$.

Proposition 5.2**.**

Let $p>1$. Let $\mu\in\mathcal{P}([0,1])$. If $\mu^{N}\in\mathcal{P}([0,1])$ is such that $\int_{0}^{1}\phi^{N}_{m}(x){\mathrm{d}}\mu^{N}(x)=\int_{0}^{1}\phi^{N}_{m}(x){\mathrm{d}}\mu(x)$ for all $m\in\{1,\dots,N\}$, then

[TABLE]

*Let us assume besides that the cost function satisfies $c(x,y)=H(y-x)$ with $H:\mathbb{R}\rightarrow\mathbb{R}_{+}$ strictly convex. There exists then a unique minimizer of (5.1.2) that we denote $\pi^{*}$.

Let $\pi^{N}\in\Pi(\mu,\nu;(\phi^{N}_{m})_{1\leq m\leq N},(\phi^{N}_{n})_{1\leq n\leq N})$, $\mu^{N}$ and $\nu^{N}$ the marginal laws of $\pi^{N}$ and assume that*

[TABLE]

Then, we have $W_{p}(\pi^{N},\pi^{*})\leq\frac{2^{1/p}}{N},$ where $W_{p}$ is defined using the $\|\|_{p}$ norm on $\mathbb{R}^{2}$.

Proof.

For $\eta\in\mathcal{P}(\mathbb{R})$, we define $F_{\eta}^{-1}(u)=\inf\{x\in\mathbb{R}:\eta((-\infty,x])\geq u\}$, that coincides with the usual inverse when $F_{\eta}$ is increasing. Let $p>1$. By Theorem 2.9 [31], we have

[TABLE]

If $F_{\mu}\left(\frac{m}{N}\right)=F_{\mu}\left(\frac{m-1}{N}\right)$, we clearly have $\int_{F_{\mu}\left(\frac{m-1}{N}\right)}^{F_{\mu}\left(\frac{m}{N}\right)}|F_{\mu}^{-1}(u)-F_{\mu^{N}}^{-1}(u)|^{p}du=0$. Otherwise, we have $F_{\mu^{N}}\left(\frac{m-1}{N}\right)=F_{\mu}\left(\frac{m-1}{N}\right)<F_{\mu}\left(\frac{m}{N}\right)=F_{\mu^{N}}\left(\frac{m}{N}\right)$, and therefore

[TABLE]

This gives $|F_{\mu}^{-1}(u)-F_{\mu^{N}}^{-1}(u)|\leq 1/N$. Since $F_{\mu}(0)=F_{\mu^{N}}(0)$, we get that $F_{\mu}^{-1}(u)=F_{\mu^{N}}^{-1}(u)=0$ for $u\in(0,F_{\mu}\left(0\right))$. We finally get $W^{p}_{p}(\mu,\mu^{N})\leq N^{-p}$.

Now, let $U\sim\mathcal{U}([0,1])$ be a uniform random variable on $[0,1]$. Still by Theorem 2.9 [31], we have $(F_{\mu}^{-1}(U),F_{\nu}^{-1}(U))\sim\pi^{*}$ and $(F_{\mu^{N}}^{-1}(U),F_{\nu^{N}}^{-1}(U))\sim\pi^{N}$. This gives a coupling between $\pi^{*}$ and $\pi^{N}$, and thus

[TABLE]

∎

In order to prove Proposition 5.1, let us introduce the following auxiliary problem. For all $N\in\mathbb{N}^{*}$, let us define

[TABLE]

and

[TABLE]

Let us introduce the following applications:

[TABLE]

and

[TABLE]

Lemma 5.3**.**

Let $N\in\mathbb{N}^{*}$. It holds that

[TABLE]

Besides, we have

[TABLE]

Proof of Lemma 5.3.

Let $\pi\in\mathcal{P}([0,1]^{2})$. Then, we write

[TABLE]

and get $\left|I(\pi)-J(D(\pi))\right|\leq\frac{K}{2N}$ since $|c(x,y)-c\left(\frac{m-\frac{1}{2}}{N},\frac{n-\frac{1}{2}}{N}\right)|\leq\frac{K}{2N}$ for $(x,y)\in T_{m}^{N}\times T_{n}^{N}$.

Let $N\in\mathbb{N}^{*}$. For all $\bar{\pi}\in\bar{\Pi}(\mu,\nu)$, defining $\pi:=\sum_{m,n=1}^{N}\bar{\pi}_{mn}\delta_{\frac{m-\frac{1}{2}}{N},\frac{n-\frac{1}{2}}{N}}$, one obtains that $\pi\in\mathcal{P}([0,1]^{2})$, $D(\pi)=\bar{\pi}$ and $I(\pi)=J(\bar{\pi})$; this implies that $I^{N}\leq J^{N}$.

Conversely, if $\pi\in\Pi(\mu,\nu;(\phi^{N}_{m})_{1\leq m\leq N},(\phi^{N}_{n})_{1\leq n\leq N})$ is chosen to satisfy $I(\pi)\leq I^{N}+\epsilon$ for some $\epsilon>0$, one gets $J^{N}\leq J(D(\pi))\leq I(\pi)+\frac{K}{2N}=I^{N}+\frac{K}{2N}+\epsilon$. Letting $\epsilon\to 0$ provides the wanted result. ∎

We also need the following auxiliary lemma.

Lemma 5.4**.**

For all $\bar{\pi}\in\bar{\Pi}^{N}(\mu,\nu)$, there exists $\bar{\pi}^{*}\in\Pi(\mu,\nu)$ such that $\bar{\pi}=D(\bar{\pi}^{*})$.

Proof of Lemma 5.4.

Let $\bar{\pi}\in\bar{\Pi}(\mu,\nu)$. We define $\bar{\pi}^{*}$ by

[TABLE]

Since $\sum_{n^{\prime}=1}^{N}\bar{\pi}_{m,n^{\prime}}=\mu(T^{N}_{m})$ and $\sum_{m=1}^{N}\bar{\pi}_{m,n}=\nu(T^{N}_{n})$, we have

[TABLE]

Also, we have $\int_{\mathcal{Y}}{\mathrm{d}}\bar{\pi}^{*}(x,y)={\mathrm{d}}\mu(x)\sum_{m=1}^{N}\mathbf{1}_{T^{N}_{m}}(x)\sum_{n=1}^{N}\frac{\bar{\pi}_{m,n}}{\sum_{n^{\prime}=1}^{N}\bar{\pi}_{m,n^{\prime}}}={\mathrm{d}}\mu(x),$ which gives $\bar{\pi}^{*}\in\Pi(\mu,\nu)$. Last, we have

[TABLE]

which precisely gives $\bar{\pi}=D(\bar{\pi}^{*})$. ∎

We are now in position to give the proof of Proposition 5.1.

Proof of Proposition 5.1.

The inclusion $\Pi(\mu,\nu;(\phi^{N}_{m})_{1\leq m\leq N},(\phi^{N}_{n})_{1\leq n\leq N})\subset\Pi(\mu,\nu)$ gives $I^{N}\leq I$.

Lemma 5.4 implies that for all $\bar{\pi}\in\bar{\Pi}^{N}(\mu,\nu)$, there exists $\bar{\pi}^{*}\in\Pi(\mu,\nu)$ such that $D(\bar{\pi}^{*})=\bar{\pi}$, and we get by Lemma 5.3 $\left|J(\bar{\pi})-I(\bar{\pi}^{*})\right|\leq\frac{K}{2N}$. Let now $\bar{\pi}\in\bar{\Pi}^{N}(\mu,\nu)$ such that $J(\bar{\pi})\leq J^{N}+\epsilon$ for some $\epsilon>0$. Then one gets that $J^{N}+\frac{K}{2N}+\epsilon\geq I\left(\bar{\pi}^{*}\right)\geq I$. Letting $\epsilon$ go to zero yields that

[TABLE]

Furthermore, Lemma 5.3 gives $J^{N}\leq I^{N}+\frac{K}{N}$ and thus $I\leq I^{N}+\frac{K}{N}$. ∎

Remark 5.2**.**

Proposition 5.1 can be easily extended to higher dimensions and in the multi-marginal case. Let us assume that $c:([0,1]^{d})^{M}\to\mathbb{R}_{+}$ be such that

[TABLE]

For $N\in\mathbb{N}^{*}$ and ${\bf m}\in\{1,\dots,N\}^{d}=:\mathcal{E}_{N}$, we consider the test function $\phi^{N}_{\bf m}(x)=\prod_{i=1}^{d}\phi^{N}_{m_{i}}(x_{i})$ for $x\in[0,1]^{d}$. Then, with $I=\inf_{\pi\in\Pi(\mu_{1},\dots,\mu_{M})}\int_{([0,1]^{d})^{M}}c(x_{1},\dots,x_{M}){\mathrm{d}}\pi(x_{1},\dots,x_{M})$ and

[TABLE]

we get similarly (it is straightforward to generalize Proposition 5.2 and we can extend the result of Lemma 5.4 by induction on $M$)

[TABLE]

Since the number of moments (i.e. of test functions) is $MN^{d}$, we see that there is a curse of dimensionality with respect to $d$, not with respect to $M$.

5.2 Piecewise affine test functions in dimension 1 on a compact set

The test functions considered are discontinuous piecewise affine functions, identical on each space. For all $N\in\mathbb{N}^{*}$ and all $1\leq m\leq N$, let us define the following discontinuous piecewise affine functions

[TABLE]

and for all $i=1,2$,

[TABLE]

Lemma 5.5**.**

Let $\mu_{1},\mu_{2}\in\mathcal{P}([0,1])$. Let $N\in\mathbb{N}^{*}$ and let us assume that for all $1\leq m\leq N$ and $i=1,2$,

[TABLE]

Then, denoting by $F_{1}:[0,1]\to[0,1]$ (resp. $F_{2}:[0,1]\to[0,1]$) the cumulative distribution function of $\mu_{1}$ (resp. $\mu_{2}$), one gets that

[TABLE]

and

[TABLE]

Proof.

We have $\phi_{m,1}+\phi_{m,2}=\mathbf{1}_{T^{N}_{m}}$ and thus, for $2\leq m\leq N$, $F_{1}\left(\frac{m}{N}\right)-F_{1}\left(\frac{m-1}{N}\right)=F_{2}\left(\frac{m}{N}\right)-F_{2}\left(\frac{m-1}{N}\right)$. Since $F_{1}(1)=F_{2}(1)=1$, this gives (5.2.2). Now, let $l=1,2$. An integration by parts yields for $1\leq m\leq N$

[TABLE]

Using (5.2.2), this gives (5.2.1). ∎

Let us remark that we may have $F_{1}(0)\not=F_{2}(0)$ under the assumptions of Lemma 5.5, since $\mu_{1}$ and $\mu_{2}$ may charge differently [math].

Let us now explain with a rough calculation why considering these test functions may lead to a convergence rate of $O(1/N^{2})$ when $c$ is $C^{1}$ with a Lipschitz gradient. Let

[TABLE]

We have $I^{N}\leq I$ and, for any $\pi\in\Pi(\mu,\nu;(\phi^{N}_{m,i}),(\phi^{N}_{n,i}))$,

[TABLE]

Thus, we have

[TABLE]

with $\pi^{1}_{mn}=\pi(T_{m}^{N}\times T_{n}^{N})$, $N\pi^{2}_{mn}=\int_{T_{m}^{N}\times T_{n}^{N}}\phi^{N}_{m,1}(x){\mathrm{d}}\pi(x,y)$ and $N\pi^{3}_{mn}=\int_{T_{m}^{N}\times T_{n}^{N}}\phi^{N}_{m,1}(y){\mathrm{d}}\pi(x,y)$. We can thus consider the linear programming problem of minimizing the right-hand-side of (5.2.4) under the constraints $\sum_{n}\pi^{1}_{mn}=\mu^{N}_{m,1}+\mu^{N}_{m,2}$, $\sum_{m}\pi^{1}_{mn}=\nu^{N}_{m,1}+\nu^{N}_{m,2}$, $\sum_{n}\pi^{2}_{mn}=\mu^{N}_{m,1}/N$, $\sum_{m}\pi^{3}_{mn}=\nu^{N}_{m,1}/N$ and $\pi^{i}_{mn}\geq 0$. Suppose for simplicity that we can find an minimum $(\pi^{*,i}_{mn})$ to this discrete problem. If we could find (similarly as Lemma 5.4) $\pi^{*}\in\Pi(\mu,\nu)$ such that $\pi^{*,1}_{mn}=\pi^{*}(T_{m}^{N}\times T_{n}^{N})$, $N\pi^{*,2}_{mn}=\int_{T_{m}^{N}\times T_{n}^{N}}\phi^{N}_{m,1}(x){\mathrm{d}}\pi^{*}(x,y)$ and $N\pi^{*,3}_{mn}=\int_{T_{m}^{N}\times T_{n}^{N}}\phi^{N}_{m,1}(y){\mathrm{d}}\pi^{*}(x,y)$, we would get then

[TABLE]

Unfortunately, such kind of a result is not obvious. Besides, we see from this derivation that the smoothness of the cost function plays an important role.

Let us recall that for $p\geq 1$, the $W_{p}$-Wasserstein distance at the power $p$, $W^{p}_{p}(\mu,\nu)$, corresponds to the cost function $c(x,y)=|x-y|^{p}$. In the following, we prove convergence results with rate $O(1/N^{2})$ for $c(x,y)=|x-y|$ and $c(x,y)=|x-y|^{2}$. In the first case, the cost function is not smooth on the diagonal, and we need to impose an extra condition on $\mu$ and $\nu$ to get this rate. We first state a first result, which is already interesting, but will be not sufficient to prove the desired convergence. Its proof is postponed in Appendix A.

Proposition 5.6**.**

Let $\mu_{1},\mu_{2}\in\mathcal{P}([0,1])$ be two probability measures with cumulative distribution functions $F_{1}$ and $F_{2}$, respectively, such that $F_{1},F_{2}\in C^{2}([0,1])$. Let us assume that for all $1\leq m\leq N$ and $i=1,2$,

[TABLE]

Then,

[TABLE]

In addition, let $m_{1}:=\min_{u\in[0,1]}F_{1}^{\prime}(u)$ and $m_{2}=\min_{u\in[0,1]}F^{\prime}_{2}(u)$ and let us assume that $m_{1}>0$ and $m_{2}>0$. Then, for all $p>1$, it holds that

[TABLE]

Remark 5.3**.**

The result of Proposition 5.6 can be extended through the triangle inequality in order to treat regular measures with different piecewise affine moments. Indeed, for $p\geq 1$:

[TABLE]

thus

[TABLE]

Thus, using Proposition 5.6, one gets that for $\mu$, $\nu$ two measures with cumulative distribution functions $F$ and $G$, respectively, such that $F,G\in C^{2}([0,1])$ and $\tilde{\mu}$, $\tilde{\nu}$ two measures with cumulative distribution functions $\tilde{F}$ and $\tilde{G}$, respectively, such that $\tilde{F},\tilde{G}\in C^{2}([0,1])$; If $\mu$ and $\tilde{\mu}$ (respectively $\nu$ and $\tilde{\nu}$) have the same $2N$ piecewise affine moments of step $1/N$, then

[TABLE]

Besides, if $m_{\mu}=\min_{u\in[0,1]}F^{\prime}(u)$, $m_{\tilde{\mu}}=\min_{u\in[0,1]}\tilde{F}^{\prime}(u)$, $m_{\nu}=\min_{u\in[0,1]}G^{\prime}(u)$ and $m_{\tilde{\nu}}=\min_{u\in[0,1]}\tilde{G}^{\prime}(u)$, are positive, one has for all $p\in\mathbb{N}^{*}$,

[TABLE]

Unfortunately, Proposition 5.6 can not be extended to non-smooth measures, as Example 5.1 below shows. However, the $O(1/N^{2})$ convergence obtained in Remark 5.3 may stay true even for non-smooth measures $\tilde{\mu}$ and $\tilde{\nu}$. This is important in our context to treat the case where $\tilde{\mu}$ and $\tilde{\nu}$ are not smooth since the solution of the MCOT problem may typically be a discrete measure that match respectively the moments of $\mu$ and $\nu$. We tackle this issue for $W_{1}$ and $W_{2}$ in the two following paragraphs.

Example 5.1**.**

In Proposition 5.6, if one of the measures (let us say $\tilde{\mu}$) is not regular enough, then the convergence in $O(1/N^{2})$ may not be true, as shown thereafter.

We consider $\mu\sim\mathcal{U}([0,1])$ and

[TABLE]

Then, for all $1\leq m\leq N$, we have

[TABLE]

and

[TABLE]

However, we have

[TABLE]

5.2.1 Convergence speed for $W_{1}$

Proposition 5.7**.**

Let $\mu,\nu,\tilde{\mu},\tilde{\nu}\in\mathcal{P}([0,1])$. Let us assume that $\mu$ and $\nu$ are absolutely continuous with respect to the Lebesgue measure and let us denote by $\rho_{\mu}$ and $\rho_{\nu}$ their density probability functions. Let us denote by $F_{\mu}$, $F_{\nu}$, $F_{\tilde{\mu}}$ and $F_{\tilde{\nu}}$ the cumulative distribution functions of $\mu$, $\nu$, $\tilde{\mu}$ and $\tilde{\nu}$ respectively. Let $N\in\mathbb{N}^{*}$. Let us assume that

[TABLE]

Let us assume in addition that the function $F_{\mu}-F_{\nu}$ changes sign at most $Q$ times for some $Q\in\mathbb{N}$. More precisely, denoting by $G:=F_{\mu}-F_{\nu}$, we assume that there exist $x_{0}=0<x_{1}<x_{2}<\cdots<x_{Q}<x_{Q+1}=1\in[0,1]$ such that for all $1\leq q\leq Q+1$,

[TABLE]

and for all $1\leq q\leq Q$,

[TABLE]

Let us also assume that $\rho_{\mu}-\rho_{\nu}\in L^{\infty}([0,1],{\mathrm{d}}x;\mathbb{R})$. Then,

[TABLE]

The key thing to notice is that we only assume regularity on the measures $\mu,\nu$, not on $\tilde{\mu},\tilde{\nu}$. The assumption that $F_{\mu}-F_{\nu}$ changes sign at most $Q$ times is related to the fact that $c(x,y)=|x-y|$ is not smooth on the diagonal: an optimal coupling is given by the inverse transform coupling, and $F_{\mu}^{-1}-F_{\nu}^{-1}$ changes sign at most $Q$ times as well. Last, remarkably, we do not need for this result to assume $F_{\mu}(m/N)=F_{\tilde{\mu}}(m/N)$ and $F_{\nu}(m/N)=F_{\tilde{\nu}}(m/N)$. Thus, it is sufficient to work with continuous piecewise affine test functions.

More precisely, for all $N\in\mathbb{N}^{*}$, let us define

[TABLE]

and for all $2\leq m\leq N$,

[TABLE]

We can check by integration by parts that $\int_{[0,1]}\psi^{N}_{1}(x){\mathrm{d}}\mu(x)=N\int_{T^{N}_{1}}F_{\mu}(x){\mathrm{d}}x$ and $\int_{[0,1]}\psi^{N}_{m}(x){\mathrm{d}}\mu(x)=N\int_{T^{N}_{m}}F_{\mu}(x){\mathrm{d}}x-N\int_{T^{N}_{m-1}}F_{\mu}(x){\mathrm{d}}x$ for $2\leq m\leq N$. Therefore, we get

[TABLE]

Last, let us remark that $\psi^{N}_{1}=\phi^{N}_{1,2}$ and $\psi^{N}_{m}=\phi^{N}_{m-1,1}-\phi^{N}_{m,2}$ for $2\leq m\leq N$ so that $\mathrm{Span}\left\{\psi^{N}_{n},\,1\leq n\leq N\right\}\subset\mathrm{Span}\left\{\phi^{N}_{n,1},\phi^{N}_{n,2},\,1\leq n\leq N\right\}$ and

[TABLE]

Corollary 5.8**.**

Let $\mu,\nu\in\mathcal{P}([0,1])$. Let us assume that $\mu$ and $\nu$ are absolutely continuous with respect to the Lebesgue measure and let us denote by $\rho_{\mu}$ and $\rho_{\nu}$ their density probability functions. Let $F_{\mu}$ and $F_{\nu}$ be their cumulative distribution functions. For all $N\in\mathbb{N}^{*}$, let us define

[TABLE]

There exists a minimizer for (5.2.15). Let us assume in addition that the function $F_{\mu}-F_{\nu}$ changes sign at most $Q$ times for some $Q\in\mathbb{N}$ (in the sense of Proposition 5.7) and that $\rho_{\mu}-\rho_{\nu}\in L^{\infty}([0,1],{\mathrm{d}}x;\mathbb{R})$. Then,

[TABLE]

In fact, looking at the proof of Proposition 5.7, it even is sufficient to assume that $\rho_{\mu}-\rho_{\nu}$ is bounded on a neighborhood of the points at which $F_{\mu}-F_{\nu}$ changes sign. For simplicity of statements, we have assumed in Proposition 5.7 and Corollary 5.8 that $\rho_{\mu}-\rho_{\nu}$ is bounded on $[0,1]$.

Proof of Corollary 5.8.

From the inclusion $\Pi(\mu,\nu;(\psi^{N}_{m})_{1\leq m\leq N},(\psi^{N}_{n})_{1\leq n\leq N})\subset\Pi(\tilde{\mu},\tilde{\nu})$, we clearly have $I^{N}\leq W_{1}(\mu,\nu)$. Using Theorem 3.1 together with Remark 3.2, since the functions $\psi_{m}^{N}$ are continuous on $[0,1]$ for all $1\leq m\leq N$, there exists $\pi^{N}\in\Pi(\mu,\nu;(\psi^{N}_{m})_{1\leq m\leq N},(\psi^{N}_{n})_{1\leq n\leq N})$ which is a minimizer to Problem (5.2.15). Let us denote by $\tilde{\mu}$ and $\tilde{\nu}$ the marginal laws of $\pi^{N}$. First, we remark that

[TABLE]

and thus

[TABLE]

Besides, it holds that for all $1\leq m\leq N$,

[TABLE]

and we therefore get (5.2.10) from (5.2.14). We can thus apply Proposition 5.7 and get the desired result. ∎

Proof of Proposition 5.7.

Let $1\leq m\leq N$. If for all $1\leq q\leq Q$, $x_{q}\notin T_{m}^{N}$, then $F_{\mu}-F_{\nu}$ remains non-negative or non-positive on $T_{m}^{N}$. Thus, using (5.2.10), it holds that

[TABLE]

where $\epsilon=1$ if $F_{\mu}-F_{\nu}\geq 0$ on $T_{m}^{N}$ and $\epsilon=-1$ if $F_{\mu}-F_{\nu}\leq 0$ on $T_{m}^{N}$. On the other hand, if there exists $1\leq q\leq Q$, such that $x_{q}\in T_{m}^{N}$, one gets

[TABLE]

since for $x\in T^{N}_{m}$, $F_{\mu}(x)-F_{\nu}(x)=\int_{x_{q}}^{x}\rho_{\mu}-\rho_{\nu}$ and $|x-x_{q}|\leq 1/N$.

Thus, as there are at most $Q$ intervals of that last type, we get

[TABLE]

i.e. $W_{1}(\mu,\nu)\leq W_{1}(\tilde{\mu},\tilde{\nu})+2\|\rho_{\mu}-\rho_{\nu}\|_{\infty}\frac{Q}{N^{2}}$. ∎

5.2.2 Convergence speed for $W_{2}$

Proposition 5.9**.**

Let $\mu,\nu,\tilde{\mu},\tilde{\nu}\in\mathcal{P}([0,1])$. Let us assume that $\mu({\mathrm{d}}x)=\rho_{\mu}(x){\mathrm{d}}x$ and $\nu({\mathrm{d}}x)=\rho_{\nu}(x){\mathrm{d}}x$ with $\rho_{\mu},\rho_{\nu}\in L^{\infty}([0,1],{\mathrm{d}}x;\mathbb{R}_{+})$. Let us denote by $F_{\mu}$, $F_{\nu}$, $F_{\tilde{\mu}}$ and $F_{\tilde{\nu}}$ the cumulative distribution functions of $\mu$, $\nu$, $\tilde{\mu}$ and $\tilde{\nu}$ respectively. Let $N\in\mathbb{N}^{*}$. Let us assume that

[TABLE]

Then,

[TABLE]

This proposition plays the same role as Proposition 5.7 for $W_{1}$. Again, the important point is that no regularity assumption is made on $\tilde{\mu}$ and $\tilde{\nu}$. We note that we no longer have restriction on the number of points where $F_{\mu}-F_{\nu}$ changes sign, which is related to the fact that $c(x,y)=(x-y)^{2}$ is smooth. Contrary to Proposition 5.7, we need here the condition (5.2.17).

Corollary 5.10**.**

Let $\mu,\nu\in\mathcal{P}([0,1])$. Let us assume that $\mu({\mathrm{d}}x)=\rho_{\mu}(x){\mathrm{d}}x$ and $\nu({\mathrm{d}}x)=\rho_{\nu}(x){\mathrm{d}}x$ with $\rho_{\mu},\rho_{\nu}\in L^{\infty}([0,1],{\mathrm{d}}x;\mathbb{R}_{+})$. Let $F_{\mu}$ and $F_{\nu}$ be their cumulative distribution functions. For all $N\in\mathbb{N}^{*}$, let us define

[TABLE]

Then,

[TABLE]

We omit the proof of Corollary 5.10 since it follows the same line as the one of Corollary 5.8. The only difference is that we do not know here if the infimum is a minimum and have to work for an arbitrary $\epsilon>0$ with $\pi\in\Pi(\mu,\nu;(\phi^{N}_{m,l})_{\scriptsize\begin{subarray}{c}1\leq m\leq N\\ 1\leq l\leq 2\end{subarray}},(\phi^{N}_{m,l})_{\scriptsize\begin{subarray}{c}1\leq m\leq N\\ 1\leq l\leq 2\end{subarray}})$ such that $\int_{[0,1]^{2}}(x-y)^{2}{\mathrm{d}}\pi(x,y)\leq I^{N}+\epsilon$. Let us also mention here that we can use Proposition 5.2 to bound the distance between an MCOT minimizer and an OT minimizer since $\phi_{m,1}+\phi_{m,2}=\mathbf{1}_{T^{N}_{m}}$.

Proof of Proposition 5.9.

From Lemma B.3 [25], we have

[TABLE]

The two terms $[F_{\mu}(x)-F_{\nu}(y)]^{+}$ and $[F_{\nu}(x)-F_{\mu}(y)]^{+}$ can be analyzed in the same way by exchanging $\mu$ and $\nu$, and we focus on the first one. Thus, we consider for $k\leq l$ the term $\alpha_{kl}:=\int_{T^{N}_{k}}\int_{T^{N}_{l}}\mathbf{1}_{x<y}[F_{\mu}(x)-F_{\nu}(y)]^{+}{\mathrm{d}}x{\mathrm{d}}y$ and denote $\tilde{\alpha}_{kl}=\int_{T^{N}_{k}}\int_{T^{N}_{l}}\mathbf{1}_{x<y}[F_{\tilde{\mu}}(x)-F_{\tilde{\nu}}(y)]^{+}{\mathrm{d}}x{\mathrm{d}}y$.

$\bullet$ If $F_{\mu}(k/N)\leq F_{\nu}((l-1)/N)$, then from (5.2.17), we have also $F_{\tilde{\mu}}(k/N)\leq F_{\tilde{\nu}}((l-1)/N)$ (note that if $l=1$, $F_{\tilde{\nu}}(0)\geq 0=F_{\nu}(0)$). Thus, $\alpha_{kl}=\tilde{\alpha}_{kl}=0$.

$\bullet$ If $F_{\nu}(l/N)\leq F_{\mu}((k-1)/N)$, then from (5.2.17), we have also $F_{\tilde{\nu}}(l/N)\leq F_{\tilde{\mu}}((k-1)/N)$, and using (5.2.18) we get for $k<l$

[TABLE]

For $k=l$, we have by using (5.2.18) and Lemma A.1 for the inequality

[TABLE]

$\bullet$ We now consider the case where $F_{\mu}(k/N)>F_{\nu}((l-1)/N)$ and $F_{\nu}(l/N)>F_{\mu}((k-1)/N)$. We can thus find $x_{0}\in T^{N}_{k}$ and $y_{0}\in T^{N}_{l}$ such that $F_{\mu}(x_{0})=F_{\nu}(y_{0})$. We then have $\forall x\in T^{N}_{k},y\in T^{N}_{l},|F_{\mu}(x)-F_{\nu}(y)|\leq|F_{\mu}(x)-F_{\mu}(x_{0})|+|F_{\nu}(y_{0})-F_{\nu}(y)|\leq\|\rho_{\mu}\|_{\infty}|x-x_{0}|+\|\rho_{\nu}\|_{\infty}|y-y_{0}|$, and thus using that $\int_{T^{N}_{k}}|x-x_{0}|{\mathrm{d}}x\leq\frac{1}{2N^{2}}$,

[TABLE]

For $1\leq k\leq N$, we note $\{l_{k},l_{k}+1,,\dots,l_{k}+n_{k}-1\}\subset\{1,\dots,N\}$ the set of $l$ such that $F_{\mu}((k-1)/N)<F_{\nu}(l/N)$ and $F_{\mu}(k/N)>F_{\nu}((l-1)/N)$. We necessarily have $l_{k+1}\geq l_{k}+n_{k}-1$ since $F_{\nu}((l_{k}+n_{k}-2)/N)<F_{\mu}(k/N)<F_{\nu}(l_{k+1}/N)$. Therefore, there is at most one element overlap between two consecutive sets, and thus $\sum_{k=1}^{N}n_{k}\leq 2N$.

Combining all cases, and taking into account the contribution of the symmetric term $[F_{\nu}(x)-F_{\mu}(y)]^{+}$ in the integral, we finally get

[TABLE]

which gives (5.2.19) ∎

6 Numerical algorithms to approximate optimal transport problems

This section presents the implementation of two algorithms for the approximation of the Optimal Transport cost. Both algorithms rely on Theorem 3.1, i.e. that the optimum of the MCOT problem is attained by a finite discrete measure $\sum_{k=1}^{2N+2}p_{k}\delta_{(x_{k},y_{k})}$. The two algorithms corresponds to the following choices:

piecewise constant test functions, 2. 2.

(regularized) piecewise linear test functions.

In the first case, the precise positions $(x_{k},y_{k})$ are useless to satisfy the moment constraints: only matters in which cell $(x_{k},y_{k})$ belongs. Thus the optimization problem is essentially discrete on a (large) finite space, for which Metropolis-Hastings algorithms are relevant. In the second case, we implement a penalized gradient algorithm to optimize the positions $(x_{k},y_{k})$ and the weights $p_{k}$.

The goal of these numerical tests is only illustrative to see the potential relevance of this approach. We do not claim that these algorithms are more efficient than other existing methods in the literature, and the improvement of our algorithms is left for future research.

6.1 Metropolis-Hastings algorithm on a finite state space

We expose in the following the principles of the Metropolis-Hastings algorithm used to compute an approximation of the OT cost. For simplicity, we do so in the case of two uni-dimensional marginal laws. However, the algorithm principles can be adapted to solve a Multimarginal MCOT problem with marginal laws defined on spaces of any finite dimension.

6.1.1 Description of the algorithm

For this algorithm, we consider the framework of Subsection 5.1, i.e. $N$ piecewise constant functions $\phi^{N}_{m}=\mathbf{1}_{T^{N}_{m}}$, $1\leq m\leq N$. and the MCOT problem (5.1.3). As mentioned above, if $(x_{k},y_{k})$ belongs to the cell $T^{N}_{i}\times T^{N}_{j}$, its position in this cell is useless regarding the moment constraints. We can therefore assume that the position minimizes the cost in this cell. For $c(x,y)=|y-x|^{2}$, this amounts to take

[TABLE]

We consider then $2N+2$ distinct cells $T^{N}_{i_{k}}\times T^{N}_{j_{k}}$, $k\in\{1,\dots,2N+2\}$. The weights associated to each cell is determined as the solution of the linear optimization of the cost associated under the constraint that the weights satisfy the moments constraints:

[TABLE]

Note that this set of constraints may be void. To start with an initial configuration $(i_{k},j_{k})_{1\leq k\leq 2N+2}$ that allows the existence of weights which satisfy the constraints, we use the inverse transform sampling between the distributions given by $(\mu_{k})_{1\leq k\leq N}$ and $(\nu_{k})_{1\leq k\leq N}$ on $\{1,\dots,N\}$. This gives in fact the optimal solution $(p_{k},(i_{k},j_{k}))_{1\leq k\leq 2N+2}$ for (5.1.3) that satisfy thus in particular the constraints. Since we want here to test the relevance of the Metropolis-Hastings algorithm in this framework, we do not want to start from the optimal solution: thus, we consider a random permutation $\sigma$ on $\{1,\dots,N\}$ and then the inverse transform sampling between the distributions given by $(\mu_{k})_{1\leq k\leq N}$ and $(\nu_{\sigma(k)})_{1\leq k\leq N}$ on $\{1,\dots,N\}$. This gives a configuration that satisfy the constraints and is not a priori optimal.

We now have to specify how the Markov chain defining the algorithm moves from one state $(i_{k},j_{k})_{1\leq k\leq 2N+2}$ to another. We denote by $N(i_{k},j_{k})=\{(i_{k}+1,j_{k}),(i_{k}-1,j_{k}),(i_{k},j_{k}+1),(i_{k},j_{k}-1)\}$ the neighboring cells of $(i_{k},j_{k})$ and

[TABLE]

the neighboring cells that are free, i.e. that are not in the current configuration. We choose randomly and uniformly a cell $l\in\{1,\dots,2N+2\}$. If $FN(i_{l},j_{l})=\emptyset$, we pick randomly another one. This rejection method amounts to choose randomly and uniformly a cell $l$ among those such that $FN(i_{l},j_{l})\not=\emptyset$. Then, we select $(i^{\prime}_{l},j^{\prime}_{l})$ uniformly on $FN(i_{l},j_{l})$ and set $(i^{\prime}_{k},j^{\prime}_{k})=(i_{k},j_{k})$ for $k\not=l$, and we accept the new configuration $(i^{\prime}_{k},j^{\prime}_{k})_{1\leq k\leq 2N+2}$ only if it allows to satisfy the constraints and with an acceptance ratio described in Algorithm 6.1. In practice, we run this Algorithm with $K\geq 2N+2$ cells, in order to increase the chance that the new configuration is compatible with the constraints.

The state space of the Markov Chain describing Algorithm 6.1 is the set of $K$ distinct elements of $\{1,\dots,N\}^{2}$. Note that we can go from any points $(i,j)$ to $(i^{\prime},j^{\prime})$ with at most $2N-2$ moves (a move consists in adding or removing one to one of the coordinate). If we ignore the problem of satisfying the constraints, we can therefore go from a configuration $(i_{k},j_{k})_{1\leq k\leq 2N+2}$ to another one $(i^{\prime}_{k},j^{\prime}_{k})_{1\leq k\leq 2N+2}$ with at most $K(2N-2)$ moves, which let think that the Doeblin condition may be satisfied. This would ensure theoretically the convergence of the algorithm converges towards the infimum

[TABLE]

and that the convergence is exponentially fast (see e.g. Section 2 of [18]).

6.1.2 Numerical examples

We tested the algorithm for the marginal laws with probability density functions

[TABLE]

We consider a number of particles $K=3N+2$ in order at each step to have more freedom among the configurations which satisfy the constraints. We present two minimizations:

• •

$N=20$ and $\beta=0.000075$

• •

$N=60$ and $\beta=0.00002$.

The evolution of the configurations through the iterations are represented for $N=20$ and $N=60$ in Figure 6.1.2 and 6.1.3. The darker the cell, the more weight it has. In green (Figures 2(f) and 3(f)) are represented the optimal configuration for the given number of moment constraints. The convergence of the numerical cost for each minimization is represented in Figure 6.1.1. The pink line represents the cost of the Optimal Transport problem we approximate, the dark blue line the one of the cost of the current configuration and the light blue one the minimum numerical cost encountered during the minimization. The green line is the cost of the optimal configuration for the given number of moment constraints, that we aim to compute.

6.2 Gradient on a penalized functional

6.2.1 Principles

We make use of Theorem 3.1 by searching optima of the MCOT problem with $N$ test functions on each space by looking for an optimal probability measure which is finitely supported on at most $2N+2$ points (note that in the multimarginal case, we can look similarly for measures supported on $DN+2$ points). This algorithm consists in penalizing moments constraints of the MCOT problem for $N$ differentiable test functions on each space ($(\phi_{m})_{1\leq m\leq N}$ and $(\psi_{n})_{1\leq n\leq N}$) and then using a gradient-type algorithm to compute the optimum.

For sake of simplicity, we consider the case of two marginal laws where the cost function $c$ is assumed to be differentiable. Let us write the position of the $2N+2$ particles by $((x_{k},y_{k}))_{1\leq k\leq 2N+2}$ and their weights by $(p_{k})_{1\leq k\leq 2N+2}$. Then, it is natural to consider the minimization of

[TABLE]

for some small parameter $\eta>0$ and under the constraints $p_{k}\geq 0$, $\sum_{k=1}^{2N+2}p_{k}=1$. To avoid the handling of these latter constraints, we prefer to consider weights $p_{k}=\frac{e^{a_{k}}}{\sum_{k=1}^{2N+2}e^{a_{k}}}$ for some $a_{k}\in\mathbb{R}$. Since the positions $(x_{k},y_{k})$ are not assumed to be different from each other, the previous minimization problem is equivalent to minimize

[TABLE]

For a fixed value of $\eta>0$, we use a projected gradient algorithm (see e.g. Algorithm 1.3.16 of [30]), to ensure that $(x_{k},y_{k})\in[0,1]^{2}$ for all $k$, together with a line search method. We implement alternated gradient steps as follows: first, a gradient step is performed on the coefficients $(a_{k})_{1\leq k\leq 2N+2}$ with $(x_{k},y_{k})_{1\leq k\leq 2N+2}$ fixed; second, a gradient step is done on the positions $(x_{k})_{1\leq k\leq 2N+2}$ with the other variables fixed; lastly, a gradient step is done on the positions $(y_{k})_{1\leq k\leq 2N+2}$ with the other variables fixed. This procedure is repeated until the norm of the projected gradient is below some error threshold. The convergence of this algorithm is ensured by Wolfe theorem (see Theorem 1.2.21 of [30]).

The example computations exposed thereafter use regularized continuous piecewise affine test functions. Remark that we do not use discontinuous piecewise affine test functions, for which we have rates of convergence for $W_{1}$ and $W_{2}$. We make this choice because the gradient algorithm that we describe above has better numerical properties for continuously differentiable test functions.

In the MCOT formulation (2.1.3) with $M=N$, minimizers of MCOT problems are the same if we consider test functions $(\bar{\phi}_{m})_{1\leq m\leq N}$ and $(\bar{\psi}_{m})_{1\leq m\leq N}$ such that ${\rm Span}((\bar{\phi}_{m})_{1\leq m\leq N})={\rm Span}(({\phi}_{m})_{1\leq m\leq N})$ and ${\rm Span}((\bar{\psi}_{m})_{1\leq m\leq N})={\rm Span}(({\psi}_{m})_{1\leq m\leq N})$. However, in the penalized version of the problem (6.2.1), the choice of the test functions has a strong impact on the convergence of the gradient algorithms. It appears that considering positive part functions (which are convex functions) greatly improves the efficiency of the procedure with respect to classical hat functions, even if both spans are identical.

Thus, for the numerical examples in 1D, we use the functions for $\epsilon>0$ and for all $N\in\mathbb{N}^{*}$,

[TABLE]

and for all $1\leq m\leq N$,

[TABLE]

which are a regularization of the functions, for all $N\in\mathbb{N}^{*}$, and $1\leq m\leq\mathbb{N}$,

[TABLE]

The vector space spanned by the restriction to $[0,1]$ of these functions is the same as the one spanned by the classical continuous piecewise affine functions (i.e. the functions $\psi^{N}_{m}$ introduced in Section 5.2.1).

For the example in dimension $2$, for $N\in\mathbb{N}^{*}$, we use the following $(N+1)^{2}$ test functions defined as follows: for all $1\leq m,n\leq N$ and $(x,y)\in[0,1]^{2}$,

[TABLE]

where for all $q\in\mathbb{Z}$,

[TABLE]

and

[TABLE]

For $1\leq m,n\leq N$, we set

[TABLE]

Those functions are a regularization of the functions $G^{N}_{m,n}(x,y)=\left(\min\left(x-\frac{m-1}{N},y-\frac{n-1}{N}\right)\right)^{+}$ with $1\leq m,n\leq N$, $G^{N}_{0,0}(x,y)=\left(\min\left(\frac{1}{N}-x,\frac{1}{N}-y\right)\right)^{+}$, $G^{N}_{0,m}(x,y)=\left(\min\left(x-\frac{m-1}{N},\frac{1}{N}-y\right)\right)^{+}$, $G^{N}_{n,0}(x,y)=\left(\min\left(\frac{1}{N}-x,y-\frac{n-1}{N}\right)\right)^{+}$. The vector space spanned by the restriction to $[0,1]^{2}$ of these functions is the same as the one spanned by the classical continuous piecewise affine functions associated to the mesh illustrated in Figure 6.2.1.

6.2.2 1D numerical example

Convergence of the algorithm

We tested the algorithm for the marginal laws with densities

[TABLE]

the quadratic cost function $c(x,y)=|y-x|^{2}$ and a fixed penalization coefficient. The exact optimal transport map between $\mu$ (abscissa) and $\nu$ (ordinate) is represented by the red line on Figures 6.2.3 and 6.2.4. We present two minimizations:

• •

$N=10$ and $1/\eta=100$

• •

$N=40$ and $1/\eta=25$.

Once each minimization process has converged, the cost for $N=10$ is 0.17805 and the one for $N=40$ is 0.17785. The evolution of the configurations through the iterations are represented for $N=10$ and $N=40$ in Figure 6.2.3 and 6.2.4. The darker the particle $(x_{k},y_{k})$, the larger its weight $p_{k}$.

And the convergence of the numerical cost for each one in Figure 6.2.2 the pink line represents the cost of the exact Optimal Transport problem that we approximate.

We note on these examples that the particles $(x_{k},y_{k})$ tend to cluster in some places. This is due to the fact that the cost function is convex and that the test functions are (up to the regularization) locally linear.

6.2.3 2D numerical example

We consider two normal marginal laws in $\mathbb{R}^{2}$: $\mu\sim\mathcal{N}_{2}(\mu,\Sigma_{\mu})$ and $\nu\sim\mathcal{N}_{2}(\nu,\Sigma_{\nu})$, with

[TABLE]

and the quadratic cost function. In this case, it is known that the optimal cost is given by $|m_{\mu}-m_{\nu}|^{2}+\text{Tr}(\Sigma_{\mu}+\Sigma_{\nu}-2(\Sigma_{\mu}^{1/2}\Sigma_{\nu}\Sigma_{\mu}^{1/2})^{1/2})$ and the optimal transport map is given by $x\mapsto m_{\nu}+\Sigma_{\mu}^{-1/2}(\Sigma_{\mu}^{1/2}\Sigma_{\nu}\Sigma_{\mu}^{1/2})^{1/2}\Sigma_{\nu}^{-1/2}$, see e.g. [19]. In Figures 5(a) and 6.2.6, the red density is $\mu$’s one and the blue one $\nu$’s. We consider regularized piecewise linear test functions on $[-4,4]^{2}$ obtained by rescaling the functions (6.2.4), (6.2.5), (6.2.6) and (6.2.7) on $[0,1]^{2}$.

We represent several iterations of the optimization for $N=36$ and $1/\eta=2$ in Figure 6.2.6, where the green arrows represent the transport map computed by the algorithm from $\mu$ (red) to $\nu$ (blue). The greener the arrow, the more weight it has.

We represent the configuration of particles at convergence on Figure 5(a) where each particles consists in a red dot linked to a blue dot. The bigger, the more mass it transports. The green dots represent the location where the red dot would have been transported if the particle was on the transport plan. The convergence of the cost is represented in Figure 5(b) where the pink line represents the cost of the Optimal Transport problem we approximate.

6.2.4 Martingale Optimal Transport numerical example

We tested the algorithm for the marginal laws $\mu$ and $\nu$ being respectively the uniform random variables on $[\frac{1}{4},\frac{3}{4}]$ and $[0,1]$, with the cost $c(x,y)=|y-x|^{3}$. Note that $\int|y-x|^{2}{\mathrm{d}}\pi(x,y)=\int|y|^{2}{\mathrm{d}}\nu(y)-\int|x|^{2}{\mathrm{d}}\mu(x)=1/16$ for any martingale coupling $\pi$. By Jensen’s inequality, we have $\int|y-x|^{3}{\mathrm{d}}\pi(x,y)\geq(1/16)^{3/2}=(1/4)^{3}$ and therefore ${\mathrm{d}}\pi(x,y)={\mathrm{d}}\mu(x)(\frac{1}{2}{\mathrm{d}}\delta_{x+1/4}+\frac{1}{2}{\mathrm{d}}\delta_{x-1/4})$ is an optimal martingale coupling and the equality condition in Jensen’s inequality shows that this is the unique optimal martingale coupling.

The two lines $y=x+1/4$ and $y=x-1/4$ characterizing the optimal martingale coupling are represented by the red lines on Figure 6.2.8. We have made one minimization with $N=10$ and $1/\eta=60$, and $N^{\prime}=10$ continuous piecewise affine moment constraints for the martingale constraint, see Problem (3.3.11). The evolution of the configurations through the iterations are represented in Figure 6.2.8. The darker the particle $(x_{k},y_{k})$, the larger the value of its weight $p_{k}$. The convergence of the numerical cost is illustrated in Figure 6.2.7, where the pink line represents the cost of the exact martingale Optimal Transport problem we approximate.

Appendix A Technical proofs of Section 5

Proof of Proposition 5.6.

Proof.

Let us first prove (5.2.5). Lemma 5.5 implies that

[TABLE]

and

[TABLE]

Then, using a Taylor expansion, as $F_{1},F_{2}\in C^{2}([0,1])$, it holds that for all $1\leq m\leq N$, all $u\in T_{m}^{N}$, and all $l=1,2$,

[TABLE]

Integrating over $T_{m}^{N}$, one gets

[TABLE]

This implies, using (A.0.1) and (A.0.2), that

[TABLE]

Thus, using (A.0.3), for all $l=1,2$ and $u\in T_{m}^{N}$,

[TABLE]

where $|\varphi_{l}(u)|\leq\frac{\|F_{l}^{\prime\prime}\|_{\infty}}{2}\left(u-\frac{m}{N}\right)^{2}$. Then, using (A.0.4), one gets that for all $u\in T_{m}^{N}$,

[TABLE]

Integrating over $T_{m}^{N}$ yields that

[TABLE]

Using the fact that

[TABLE]

we obtain that

[TABLE]

Let us now prove (5.2.6). The main result needed is the expression of the Wasserstein distance in term of the cumulative distribution functions (cdf) and not their inverse (see [25] Lemma B.3), which holds for $p>1$,

[TABLE]

because the reasoning of the beginning of this proof introduced a control on the norm between the cdf of the marginal law and the cdf of a marginal law satisfying the same moments.

Then, one can proceed with the following induction. Suppose that we know for $p\in\mathbb{N}^{*}$ that

[TABLE]

which holds for $p=1$. Then,

[TABLE]

Let us treat the first term of the sum, as the second one can be treated symmetrically. If $F(x)\geq\tilde{F}(x)$, we can define $y_{x}=\tilde{F}^{-1}\left(F(x)\right)$ and because $\tilde{F}:[0,1]\rightarrow[\tilde{F}(0),1]$ is continuous increasing, and we have

[TABLE]

Thus, by using (A.0.3), we get

[TABLE]

[TABLE]

where we used the formula bounding the difference between the cdf (A.0.5).

Therefore, as $m_{\mu}>0$ and $m_{\tilde{\mu}}>0$, and using the symmetry of the the formula (A.0.7), one gets

[TABLE]

Hence, using the induction hypothesis (A.0.8), we obtain that (A.0.8) holds for $p+1$, which gives the claim. ∎

Lemma A.1**.**

Let $\mu\in\mathcal{P}([0,1])$ and $F_{\mu}$ its cumulative distribution function. Let $N\in\mathbb{N}^{*}$. Then, for any $1\leq m\leq N$, we define $x^{N}_{m}\in T^{N}_{m}$ by

[TABLE]

and $\hat{\mu}^{N}=F_{\mu}(0)\delta_{0}+\sum_{m=1}^{N}(F_{\mu}\left(\frac{m}{N}\right)-F_{\mu}\left(\frac{m-1}{N}\right))\delta_{x_{m}}$. Then, we have for all $1\leq m\leq N$,

[TABLE]

Besides, if $\mu({\mathrm{d}}x)=\rho_{\mu}(x){\mathrm{d}}x$ with a $\rho_{\mu}\in L^{\infty}([0,1],{\mathrm{d}}x;\mathbb{R}_{+})$ and $\tilde{\mu}\in\mathcal{P}([0,1])$ is such that $F_{\tilde{\mu}}\left(\frac{m}{N}\right)=F_{\mu}\left(\frac{m}{N}\right)$ and $\int_{T^{N}_{m}}F_{\tilde{\mu}}=\int_{T^{N}_{m}}F_{\mu}$, we have

[TABLE]

Proof.

If $F_{\mu}\left(\frac{m}{N}\right)>F_{\mu}\left(\frac{m-1}{N}\right)$, we have

[TABLE]

since $F_{\mu}$ is non-decreasing and right-continuous. Therefore, there is a unique $x^{N}_{m}\in T^{N}_{m}$ such that

[TABLE]

which is precisely the definition of $x^{N}_{m}$. By construction, we have $F_{\hat{\mu}^{N}}\left(\frac{m}{N}\right)=F_{\mu}\left(\frac{m}{N}\right)$ and the previous equation gives

[TABLE]

when $F_{\mu}\left(\frac{m}{N}\right)>F_{\mu}\left(\frac{m-1}{N}\right)$ (this identity is obvious if $F_{\mu}\left(\frac{m}{N}\right)=F_{\mu}\left(\frac{m-1}{N}\right)$). Last, since for $x\in T^{N}_{m}$, $F_{\mu}\left(\frac{m-1}{N}\right)\leq F_{\mu}(x)\leq F_{\mu}\left(\frac{m}{N}\right)$, we get that $x\mapsto\int_{\frac{m-1}{N}}^{x}(F_{\mu}-F_{\hat{\mu}^{N}})$ is non-decreasing on $[\frac{m-1}{N},x^{N}_{m}]$, non-increasing on $[x^{N}_{m},\frac{m}{N}]$ and vanishes for $x\in\{\frac{m-1}{N},\frac{m}{N}\}$: it is therefore non-negative on $T^{N}_{m}$.

Now, let us assume that $\mu$ has a bounded density probability function $\rho_{\mu}$. We have

[TABLE]

and therefore

[TABLE]

Now, we observe that we either have $\int_{T^{N}_{m}}\left(\int_{\frac{m-1}{N}}^{x}F_{\mu}\right){\mathrm{d}}x\leq\int_{T^{N}_{m}}\left(\int_{\frac{m-1}{N}}^{x}F_{\tilde{\mu}}\right){\mathrm{d}}x$ or $\int_{T^{N}_{m}}\left(\int_{\frac{m-1}{N}}^{x}F_{\mu}\right){\mathrm{d}}x\geq\int_{T^{N}_{m}}\left(\int_{\frac{m-1}{N}}^{x}F_{\tilde{\mu}}\right){\mathrm{d}}x\geq\int_{T^{N}_{m}}\left(\int_{\frac{m-1}{N}}^{x}F_{\hat{\mu}^{N}}\right){\mathrm{d}}x$. In the first case, the claim is obvious. In the second one, we then have

[TABLE]

and we get the result using (A.0.11). ∎

Acknowledgements

The Labex Bézout is acknowledged for funding the PhD thesis of Rafaël Coyaud. Aurélien Alfonsi benefited from the support of the “Chaire Risques Financiers”, Fondation du Risque. We are very grateful to Luca Nenna for stimulating discussions.