# Low-rank Optimal Transport:
Approximation, Statistics and Debiasing

Meyer Scetbon
CREST, ENSAE
<EMAIL_ADDRESS>&Marco Cuturi
Apple and CREST, ENSAE
<EMAIL_ADDRESS>
###### Abstract

The matching principles behind optimal transport (OT) play an increasingly
important role in machine learning, a trend which can be observed when OT is
used to disambiguate datasets in applications (e.g. single-cell genomics) or
used to improve more complex methods (e.g. balanced attention in transformers
or self-supervised learning). To scale to more challenging problems, there is
a growing consensus that OT requires solvers that can operate on millions, not
thousands, of points. The low-rank optimal transport (LOT) approach advocated
in Scetbon et al. (2021) holds several promises in that regard, and was shown
to complement more established entropic regularization approaches, being able
to insert itself in more complex pipelines, such as quadratic OT. LOT
restricts the search for low-cost couplings to those that have a low-
nonnegative rank, yielding linear time algorithms in cases of interest.
However, these promises can only be fulfilled if the LOT approach is seen as a
legitimate contender to entropic regularization when compared on properties of
interest, where the scorecard typically includes theoretical properties
(statistical complexity and relation to other methods) or practical aspects
(debiasing, hyperparameter tuning, initialization). We target each of these
areas in this paper in order to cement the impact of low-rank approaches in
computational OT.

##  1 Introduction

Optimal transport (OT) is used across data-science to put in correspondence
different sets of observations. These observations may come directly from
datasets, or, in more advanced applications, depict intermediate layered
representations of data. OT theory provides a single grammar to describe and
solve increasingly complex matching problems (linear, quadratic, regularized,
unbalanced, etc…), making it gain a stake in various areas of science such as
as single-cell biology Schiebinger et al. (2019); Yang et al. (2020); Demetci
et al. (2020), imaging Schmitz et al. (2018); Heitz et al. (2020); Zheng et
al. (2020) or neuroscience Janati et al. (2020); Koundal et al. (2020).

Regularized approaches to OT.  Solving OT problems at scale poses, however,
formidable challenges. The most obvious among them is computational: the
Kantorovich (1942) problem on discrete measures of size $n$ is a linear
program that requires $O(n^{3}\log n)$ operations to be solved. A second and
equally important challenge lies in the estimation of OT in high-dimensional
settings, since it suffers from the curse-of-dimensionality Fournier and
Guillin (2015). The advent of regularized approaches, such as entropic
regularization (Cuturi, 2013), has pushed these boundaries thanks for faster
algorithms (Chizat et al., 2020; Clason et al., 2021) and improved statistical
aspects (Genevay et al., 2018a). Despite these clear strengths, regularized OT
solvers remain, however, costly as they typically scale quadratically in the
number of observations.

Scaling up OT using low-rank couplings.  While it is always intuitively
possible to reduce the size of measures (e.g. using $k$-means) prior to
solving an OT between them, a promising line of work proposes to combine both
(Forrow et al., 2019; Scetbon et al., 2021, 2022). Conceptually, these low-
rank approaches solve simultaneously both an optimal clustering/aggregation
strategy with the computation of an effective transport. This intuition rests
on an explicit factorization of couplings into two sub-couplings. This has
several computational benefits, since its computational cost becomes linear in
$n$ if the ground cost matrix seeded to the OT problem has itself a low-rank.
While these computational improvements, mostly demonstrated empirically, hold
several promises, the theoretical properties of these methods are not yet well
established. This stands in stark contrast to the Sinkhorn approach, which is
comparatively much better understood.

Our Contributions.  The goal of this paper is to advance our knowledge,
understanding and practical ability to leverage low-rank factorizations in OT.
This paper provides five contributions, targeting theoretical and practical
properties of LOT: (i) We derive the rate of convergence of the low-rank OT to
the true OT with respect to the non-nnegative rank parameter. (ii) We make a
first step towards a better understanding of the statistical complexity of LOT
by providing an upper-bound of the statistical error, made when estimating LOT
using the plug-in estimator; that upper-bound has a parametric rate
$\mathcal{O}(\sqrt{1/n})$ that is independent of the dimension. (iii) We
introduce a debiased version of LOT: as the Sinkhorn divergence (Feydy et al.,
2018), we show that debiased LOT is nonnegative, metrizes the weak
convergence, and that it interpolates between the maximum mean discrepancy
(Gretton et al., 2012) and OT. (iv) We exhibit links between the bias induced
by the low-rank factorization and clustering methods. (v) We propose practical
strategies to tune the step-length and the initialization of the algorithm in
(Scetbon et al., 2021).

Notations. We consider $(\mathcal{X},d_{\mathcal{X})}$ and
$(\mathcal{Y},d_{\mathcal{Y}})$ two nonempty compact Polish spaces and we
denote $\mathcal{M}_{1}^{+}(\mathcal{X})$ (resp.
$\mathcal{M}_{1}^{+}(\mathcal{Y})$) the space of positive Radon probability
measures on $\mathcal{X}$ (resp. $\mathcal{Y}$). For all $n\geq 1$, we denote
$\Delta_{n}$ the probability simplex of size $n$ and $\Delta_{n}^{*}$ the
subset of $\Delta_{n}$ of positive histograms. We write
$\mathbf{1}_{n}\triangleq(1,\dots,1)^{T}\in\mathbb{R}^{n}$ and we denote
similarly $\|\cdot\|_{2}$ the Euclidean norm and the Euclidean distance
induced by this norm depending on the context.

##  2 Background on Low-rank Optimal Transport

Let $\mu\in\mathcal{M}_{1}^{+}(\mathcal{X})$,
$\nu\in\mathcal{M}_{1}^{+}(\mathcal{Y})$ and
$c:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}_{+}$ a nonnegative and
continuous function. The Kantorovitch formulation of optimal transport between
$\mu$ and $\nu$ is defined by

$\displaystyle\text{OT}_{c}(\mu,\nu)\triangleq\min_{\pi\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi(x,y)\,,$
(1)

where the feasible set is the set of distributions over the product space
$\mathcal{X}\times\mathcal{Y}$ with marginals $\mu$ and $\nu$:

$\Pi(\mu,\nu)\triangleq\left\\{\pi\in\mathcal{M}_{1}^{+}(\mathcal{X}\times\mathcal{Y})~{}~{}\text{s.t.}~{}~{}P_{1\\#\pi}=\mu\text{,}~{}P_{2\\#\pi}=\nu\right\\}\;,$

with $P_{1\\#\pi}$ (resp. $P_{2\\#\pi}$), the pushforward probability measure
of $\pi$ using the projection maps $P_{1}(x,y)=x$ (resp. $P_{2}(x,y)=y$). When
there exists an optimal coupling solution of (1) supported on a graph of a
function, we call such function a Monge map. In the discrete setting, one can
reformulate the optimal transport problem as a linear program over the space
of nonnegative matrices satisfying the marginal constraints. More precisely,
let $a$ and $b$ be respectively elements of $\Delta_{n}^{*}$ and
$\Delta_{m}^{*}$ and let also $\mathbf{X}\triangleq\\{x_{1},\dots,x_{n}\\}$
and $\mathbf{Y}\triangleq\\{y_{1},\dots,y_{m}\\}$ be respectively two subsets
of $\mathcal{X}$ and $\mathcal{Y}$. By denoting
$\mu_{a,\mathbf{X}}\triangleq\sum_{i=1}^{n}a_{i}\delta_{x_{i}}$ and
$\nu_{b,\mathbf{Y}}\triangleq\sum_{j=1}^{m}b_{j}\delta_{y_{j}}$ the two
discrete distributions associated and writing
$C\triangleq[c(x_{i},y_{j})]_{i,j}$, the discrete optimal transport problem
can be formulated as

$\displaystyle\text{OT}_{c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})=\min_{P\in\Pi_{a,b}}\langle
C,P\rangle\,~{}~{}\text{where}~{}~{}\Pi_{a,b}\triangleq\\{P\in\mathbb{R}_{+}^{n\times
m}\text{ s.t. }P\mathbf{1}_{m}=a,P^{T}\mathbf{1}_{n}=b\\}\,.$ (2)

Scetbon et al. (2021) propose to constrain the discrete optimal transport
problem to couplings that have a low-nonnegative rank:

######  Definition 1.

Given $M\in\mathbb{R}_{+}^{n\times m}$, the nonnegative rank of $M$ is defined
by:
$\textup{rk}_{+}(M)\triangleq\min\\!\left\\{q|M=\sum_{i=1}^{q}R_{i},\forall
i,\textup{rk}(R_{i})=1,R_{i}\geq 0\right\\}.$

Note that for any $M\in\mathbb{R}_{+}^{n\times m}$, we always have that
$\textup{rk}_{+}(M)\leq\min(n,m)$. For $r\geq 1$, we consider the set of
couplings satisfying marginal constaints with nonnegative-rank of at most $r$
as $\Pi_{a,b}(r)\triangleq\\{P\in\Pi_{a,b}\text{,}~{}~{}\textup{rk}_{+}(P)\leq
r\\}$. The discrete Low-rank Optimal Transport (LOT) problem is defined by:

$\displaystyle\text{LOT}_{r,c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})\triangleq\min_{P\in\Pi_{a,b}(r)}\langle
C,P\rangle\,.$ (3)

To solve this problem, Scetbon et al. (2021) show that Problem (3) is
equivalent to

$\displaystyle\min_{(Q,R,g)\in\mathcal{C}_{1}(a,b,r)\cap\mathcal{C}_{2}(r)}\langle
C,Q\operatorname{diag}(1/g)R^{T}\rangle\,,$ (4)

where
$\mathcal{C}_{1}(a,b,r)\triangleq\Big{\\{}(Q,R,g)\in\mathbb{R}_{+}^{n\times
r}\times\mathbb{R}_{+}^{m\times r}\times(\mathbb{R}_{+}^{*})^{r}\text{ s.t.
}Q\mathbf{1}_{r}=a,R\mathbf{1}_{r}=b\Big{\\}}$ and
$\mathcal{C}_{2}(r)\triangleq\Big{\\{}(Q,R,g)\in\mathbb{R}_{+}^{n\times
r}\times\mathbb{R}_{+}^{m\times r}\times\mathbb{R}^{r}_{+}\text{ s.t.
}Q^{T}\mathbf{1}_{n}=R^{T}\mathbf{1}_{m}=g\Big{\\}}\;.$ They propose to solve
it using a mirror descent scheme and prove the non-asymptotic stationary
convergence of their algorithm. While Scetbon et al. (2021) only focus on the
discrete setting, we consider here its extension for arbitrary probability
measures. Following (Forrow et al., 2019), we define the set of rank-$r$
couplings satisfying marginal constraints by:

$\Pi_{r}(\mu,\nu)\triangleq\\{\pi\in\Pi(\mu,\nu):\exists(\mu_{i})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{X})^{r},~{}(\nu_{i})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{Y})^{r},~{}\lambda\in\Delta_{r}^{*}\text{~{}s.t.~{}}\pi=\sum_{i=1}^{r}\lambda_{i}\mu_{i}\otimes\nu_{i}\\}\;.$

This more general definition of LOT between
$\mu\in\mathcal{M}_{1}^{+}(\mathcal{X})$ and
$\nu\in\mathcal{M}_{1}^{+}(\mathcal{Y})$ reads:

$\displaystyle\text{LOT}_{r,c}(\mu,\nu)\triangleq\inf_{\pi\in\Pi_{r}(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi(x,y)\;.$
(5)

Note that this definition of $\text{LOT}_{r,c}$ is consistent as it coincides
with the one defined in (3) on discrete probability measures. Observe also
that $\Pi_{r}(\mu,\nu)$ is compact for the weak topology and therefore the
infimum in (5) is attained. See Appendix A for more details.

##  3 Approximation Error of LOT to original OT as a function of rank

Our goal in this section is to obtain a control of the error induced by the
low-rank constraint when trying to approximate the true OT cost. We provide
first a control of the approximation error in the discrete setting. The proof
is given in Appendix B.1.

######  Proposition 1.

Let $n,m\geq 2$,
$\mathbf{X}\triangleq\\{x_{1},\dots,x_{n}\\}\subset\mathcal{X}$,
$\mathbf{Y}\triangleq\\{y_{1},\dots,y_{m}\\}\subset\mathcal{Y}$ and
$a\in\Delta_{n}^{*}$ and $b\in\Delta_{m}^{*}$. Then for $2\leq
r\leq\min(n,m)$, we have that

$\displaystyle|\mathrm{LOT}_{r,c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})-\mathrm{OT}_{c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})|\leq\|C\|_{\infty}~{}\text{ln}(\min(n,m)/(r-1))$

######  Remark 1.

Note that this result improves the control obtained in (Liu et al., 2021),
where they obtain that
$|\mathrm{LOT}_{r,c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})-\mathrm{OT}_{c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})|\lesssim\|C\|_{\infty}~{}\sqrt{nm}(\min(n,m)-r)$
as we have for any $z,z^{\prime}\geq 1$,
$|\ln(z)-\ln(z^{\prime})|\leq|z-z^{\prime}|$.

It is in fact possible to obtain another control of the approximation error by
partitioning the space where the measures are supported. For that purpose let
us introduce the notion of entropy numbers.

######  Definition 2.

Let $(\mathcal{Z},d)$ a metric space, $\mathcal{W}\subset\mathcal{Z}$ and
$k\geq 1$ an integer. Then by denoting
$B_{\mathcal{Z}}(z,\varepsilon)\triangleq\\{y\in\mathcal{Z}:~{}d(z,y)\leq\varepsilon\\}$,
we define the $k$-th (dyadic) entropy number of $\mathcal{W}$ as

$\displaystyle\mathcal{N}_{k}(\mathcal{W},d)\triangleq\inf\\{\varepsilon~{}\text{s.t.}~{}\exists~{}z_{1},\dots,z_{2^{k}}\in\mathcal{Z}:~{}~{}\mathcal{W}\subset\cup_{i=1}^{2^{k}}B_{\mathcal{Z}}(z_{i},\varepsilon)\\}\,.$

For example, any compact set $\mathcal{W}$ of $\mathbb{R}^{d}$ admits finite
entropy numbers, and by denoting $R\triangleq\sup_{w\in\mathcal{W}}\|w\|_{2}$,
we have $\mathcal{N}_{k}(\mathcal{W},\|\cdot\|_{2})\leq 4R/{2^{k/d}}$. We
obtain next a control of the approximation error of $\text{LOT}_{r,c}$ to the
true OT cost using entropy numbers (see proof in Appendix B.2).

######  Proposition 2.

Let $\mu\in\mathcal{M}_{1}^{+}(\mathcal{X})$,
$\nu\in\mathcal{M}_{1}^{+}(\mathcal{Y})$ and assume that $c$ is $L$-Lipschitz
w.r.t. $x$ and $y$. Then for any $r\geq 1$, we have

$\displaystyle|\mathrm{LOT}_{r,c}(\mu,\nu)-\mathrm{OT}_{c}(\mu,\nu)|\leq
2L\max(\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{X},d_{\mathcal{X}}),\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{Y},d_{\mathcal{Y}}))$

This results in the following bound for the $p$-Wasserstein distance for any
$p\geq 1$ on $\mathbb{R}^{d}$.

######  Corollary 1.

Let $d\geq 1$, $p\geq 1$, $\mathcal{X}$ a compact subspace of $\mathbb{R}^{d}$
and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$. By denoting
$R\triangleq\sup_{x\in\mathcal{X}}\|x\|_{2}$, we obtain that for any $r\geq
1$,

$\displaystyle|\mathrm{LOT}_{r,\|\cdot\|_{2}^{p}}(\mu,\nu)-\mathrm{OT}_{\|\cdot\|_{2}^{p}}(\mu,\nu)|\leq
4dp\frac{(8R^{2})^{p}}{r^{p/2d}}\;.$

As per the Proof of Proposition 2 we can provide a tighter control, assuming a
Monge map exists.

######  Corollary 2.

Under the same assumptions of Proposition 2 and by assuming in addition that
there exists a Monge map solving $\text{OT}_{c}(\mu,\nu)$, we obtain that for
any $r\geq 1$,

$\displaystyle|\mathrm{LOT}_{r,c}(\mu,\nu)-\mathrm{OT}_{c}(\mu,\nu)|\leq
L\mathcal{N}_{\lfloor\log_{2}(r)\rfloor}(\mathcal{Y},d_{\mathcal{Y}})\,.$

When $\mathcal{X}=\mathcal{Y}$ are a subspaces of $\mathbb{R}^{d}$, a
sufficient condition for a Monge map to exists is that either $\mu$ or $\nu$
is absolutely continuous with respect to the Lebesgue measure and that $c$ is
of the form $h(x-y)$ where $h:\mathcal{X}\rightarrow\mathbb{R}_{+}$ is a
strictly convex function (Santambrogio, 2015, Theorem 1.17). Therefore if
$\mu$ is absolutely continuous with respect to the Lebesgue measure, we obtain
for any $r\geq 1$ and $p>1$

$\displaystyle|\text{LOT}_{r,\|\cdot\|_{2}^{p}}(\mu,\nu)-\text{OT}_{\|\cdot\|_{2}^{p}}(\mu,\nu)|\leq
2dp\frac{(8R^{2})^{p}}{r^{p/d}}\;.$

##  4 Sample Complexity of LOT

We now focus on the statistical performance of the plug-in estimator for LOT.
In the following we assume that $\mathcal{X}=\mathcal{Y}$ for simplicity.
Given $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$, we denote the empirical
measures associated
$\hat{\mu}_{n}\triangleq\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ and
$\hat{\nu}_{n}\triangleq\frac{1}{n}\sum_{i=1}^{n}\delta_{Y_{i}}$, where
$(X_{i},Y_{i})_{i=1}^{n}$ are sampled independently from $\mu\otimes\nu$. We
consider the plug-in estimator defined as
$\text{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})$, and we aim at quantifying the
rate at which it converges towards the true low-rank optimal transport cost
$\text{LOT}_{r,c}(\mu,\nu)$. Before doing so, in the next Proposition we show
that this estimator is consistent on compact spaces. The proof is given in
Appendix B.3.

######  Proposition 3.

Let $r\geq 1$ and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$, then
$\mathrm{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})\xrightarrow[n\rightarrow+\infty]{}\mathrm{LOT}_{r,c}(\mu,\nu)~{}~{}a.s.$

Next we aim at obtaining the convergence rates of our plug-in estimator. In
the following Proposition, we obtain a non-asymptotic upper-bound of the
statistical error. See Appendix B.4 for the proof.

######  Proposition 4.

Let $r\geq 1$ and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$. Then, there
exists a constant $K_{r}$ such that for any $\delta>0$ and $n\geq 1$, we have,
with a probability of at least $1-2\delta$, that

$\displaystyle\mathrm{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})\leq\mathrm{LOT}_{r,c}(\mu,\nu)+11\|c\|_{\infty}\sqrt{\frac{r}{n}}+K_{r}\|c\|_{\infty}\left[\sqrt{\frac{\log(40/\delta)}{n}}+\frac{\sqrt{r}\log(40/\delta)}{n}\right]\,.$

This result is, to the best of our knowledge, the first attempt at providing a
statistical control of low-rank optimal transport. We provide an upper-bound
of the plug-in estimator which converges towards $\text{LOT}_{r,c}$ at a
parametric rate and which is independent of the dimension on general compact
metric spaces. While we fall short of providing a lower bound that could match
that upper bound, and therefore provide a complete statistical complexity
result, we believe this result might provide a first explanation on why, in
practice, $\text{LOT}_{r,c}$ displays better statistical properties than
unregularized OT and its curse of dimensionality (Dudley, 1969). In addition,
that upper bound compares favorably to known results on entropic optimal
transport. The rate of entropy regularized OT does not depend on the ambient
dimension with respect to $n$, but carries an exponential dependence in
dimension with respect to the regularization parameter $\varepsilon$ (Mena and
Niles-Weed, 2019). By contrast, the term associated with the nonnegative rank
$r$ in our bound has no direct dependence on dimension.

Our next aim is to obtain an explicit rate with respect to $r$ and $n$. In
Proposition 4, we cannot control explicitly $K_{r}$ in the general setting.
Indeed, in our proof, we obtain that $K_{r}\triangleq
14/\min_{i}\lambda_{i}^{*}$ where
$(\lambda_{i}^{*})_{i=1}^{r}\in\Delta_{r}^{*}$ are the weights involved in the
decomposition of one optimal solution of the true $\text{LOT}_{r,c}(\mu,\nu)$.
Therefore the control of $K_{r}$ requires additional assumptions on the
optimal solutions of $\text{LOT}_{r,c}(\mu,\nu)$. In the following
Proposition, we obtain an explicit upper-bound of the plug-in estimator with
respect to $r$ and $n$ in the asymptotic regime.

######  Proposition 5.

Let $r\geq 1$, $\delta>0$ and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$.
Then there exists a constant $N_{r,\delta}$ such that if $n\geq N_{r,\delta}$
then with a probability of at least $1-2\delta$, we have

$\displaystyle\mathrm{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})\leq\mathrm{LOT}_{r,c}(\mu,\nu)+11\|c\|_{\infty}\sqrt{\frac{r}{n}}+77\|c\|_{\infty}\sqrt{\frac{\log(40/\delta)}{n}}\;.$

Note that one cannot recover the result obtained in Proposition 5 from the one
obtained in Proposition 4 as we have that $K_{r}\geq
14r\xrightarrow[r\rightarrow+\infty]{}+\infty$. In order to prove the above
result, we use an extension of the McDiarmid’s inequality when differences are
bounded with high probability (Kutin, 2002). See proof in Appendix B.5 for
more details.

##  5 Debiased Formulation of $\mathrm{LOT}$

We introduce here the debiased formulation of $\mathrm{LOT}_{r,c}$ and show
that it is able to distinguish two distributions, metrize the convergence in
law and can be used as a new objective in order to learn distributions. We
focus next on the debiasing terms involving measures with themselves
$\mathrm{LOT}_{r,c}(\mu,\mu)$ in this new divergence, and show that they can
be interpreted as defining a new clustering method generalizing $k$-means for
any geometry.

###  5.1 On the Proprieties of the Debiased Low-rank Optimal Transport

When it comes to learn (or generate) a distribution in ML applications given
samples, it is crucial to consider a divergence that is able to distinguish
between two distributions and metrize the convergence in law. In general,
$\mathrm{LOT}_{r,c}(\mu,\mu)\neq 0$ and the minimum of
$\mathrm{LOT}_{r,c}(\nu,\mu)$ with respect to $\nu$ will not necessarily
recover $\mu$. In order to alleviate this issue we propose a debiased version
of $\mathrm{LOT}_{r,c}$ defined for any
$\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$ as

$\displaystyle\mathrm{DLOT}_{r,c}(\mu,\nu)\triangleq\mathrm{LOT}_{r,c}(\mu,\nu)-\frac{1}{2}[\mathrm{LOT}_{r,c}(\mu,\mu)+\mathrm{LOT}_{r,c}(\nu,\nu)]\;.$

Note that $\mathrm{DLOT}_{r,c}(\nu,\nu)=0$. In the next Proposition, we show
that, as the Sinkhorn divergence (Genevay et al., 2018b; Feydy et al., 2018),
$\mathrm{DLOT}_{r,c}$ interpolates between the Maximum Mean Discrepancy (MMD)
and OT. See proof in Appendix B.6.

######  Proposition 6.

Let $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$. Let us assume that $c$ is
symmetric, then we have

$\mathrm{DLOT}_{1,c}(\mu,\nu)=\frac{1}{2}\int_{\mathcal{X}^{2}}-c(x,y)d[\mu-\nu]\otimes
d[\mu-\nu](x,y)\;.$

If in addition we assume the $c$ is Lipschitz w.r.t to $x$ and $y$, then we
have

$\displaystyle\mathrm{DLOT}_{r,c}(\mu,\nu)\xrightarrow[r\rightarrow+\infty]{}\mathrm{OT}_{c}(\mu,\nu)\;.$

Next, we aim at showing some useful properties of the debiased low-rank OT for
machine learning applications. For that purpose, let us first recall some
definitions.

######  Definition 3.

We say that the cost $c:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}_{+}$
is a semimetric on $\mathcal{X}$ if for all $x,x^{\prime}\in\mathcal{X}$,
$c(x,x^{\prime})=c(x^{\prime},x)$ and $c(x,x^{\prime})=0$ if and only if
$x=x^{\prime}$. In addition we say that $c$ has a negative type if $\forall
n\geq 2$, $x_{1},\dots,x_{n}\in\mathcal{X}$ and
$\alpha_{1},\dots,\alpha_{n}\in\mathbb{R}$ such that
$\sum_{i=1}^{n}\alpha_{i}=0$,
$\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}c(x_{i},x_{j})\leq 0$. We say also that
$c$ has a strong negative type if for all
$\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$,
$\mu\neq\nu\implies\int_{\mathcal{X}^{2}}c(x,y)d[\mu-\nu]\otimes[\mu-\nu]<0$.

Note that if $c$ has a strong negative type, then $c$ has a negative type too.
For example, all Euclidean spaces and even separable Hilbert spaces endowed
with the metric induced by their inner products have strong negative type.
Also, on $\mathbb{R}^{d}$, the squared Euclidean distance has a negative type
(Sejdinovic et al., 2013).

We can now provide stronger geometric guarantees for $\mathrm{DLOT}_{r,c}$. In
the next Proposition, we show that for a large class of cost functions,
$\mathrm{DLOT}_{r,c}$ is nonnegative, able to distinguish two distributions,
and metrizes the convergence in law. The proof is given in Appendix B.8.

######  Proposition 7.

Let $r\geq 1$, and let us assume that $c$ is a semimetric of negative type.
Then for all $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$, we have that

$\mathrm{DLOT}_{r}(\mu,\nu)\geq 0\;.$

In addition, if $c$ has strong negative type then we have also that

$\displaystyle\mathrm{DLOT}_{r,c}(\mu,\nu)=0$
$\displaystyle\iff\mu=\nu~{}~{}\text{and}$
$\displaystyle\mu_{n}\rightarrow\mu$
$\displaystyle\iff\mathrm{DLOT}_{r,c}(\mu_{n},\mu)\rightarrow 0\;.$

where the convergence of the sequence of probability measures considered is
the convergence in law.

Observe that when $c$ has strong negative type,
$\nu\rightarrow\mathrm{DLOT}_{r,c}(\nu,\mu)\geq 0$ and it admits a unique
global minimizer at $\nu=\mu$. Therefore, $\mathrm{DLOT}_{r,c}$ has desirable
properties to be used as a loss. It is also worth noting that, in order to
obtain the metrization of the convergence in law, we show the following
Proposition. See proof in Appendix B.7.

######  Proposition 8.

Let $r\geq 1$ and $(\mu_{n})_{n\geq 0}$ and $(\nu_{n})_{n\geq 0}$ two
sequences of probability measures such that $\mu_{n}\rightarrow\mu$ and
$\nu_{n}\rightarrow\nu$ with respect to the convergence in law. Then we have
that

$\mathrm{LOT}_{r,c}(\mu_{n},\nu_{n})\rightarrow\mathrm{LOT}_{r,c}(\mu,\nu)\;.$

###  5.2 Low-Rank Transport Bias and Clustering

We turn next to the debiasing terms appearing in DLOT and exhibit links
between LOT and clustering methods. Indeed, in the discrete setting, the low-
rank bias of a probability measure $\mu$ defined as
$\text{LOT}_{k,c}(\mu,\mu)$ can be seen as a generalized version of the
$k$-means method for any geometry. In the next Proposition we obtain a new
formulation of $\text{LOT}_{k,c}(\mu,\mu)$ viewed as a general clustering
method on arbitrary metric space. See proof in Appendix B.9.

######  Proposition 9.

Let $n\geq k\geq 1$,
$\mathbf{X}\triangleq\\{x_{1},\dots,x_{n}\\}\subset\mathcal{X}$ and
$a\in\Delta_{n}^{*}$. If $c$ is a semimetric of negative type, then by
denoting $C=(c(x_{i},x_{j}))_{i,j}$, we have that

$\displaystyle\mathrm{LOT}_{k,c}(\mu_{a,\mathbf{X}},\mu_{a,\mathbf{X}})=\min_{Q}\langle
C,Q\text{diag}(1/Q^{T}\mathbf{1}_{n})Q^{T}\rangle~{}~{}\text{s.t.}~{}~{}Q\in\mathbb{R}^{n\times
k}_{+}~{}~{}\text{,}~{}~{}Q\mathbf{1}_{k}=a\;.$ (6)

Let us now explain in more details the link between (6) and $k$-means. When
$\mathcal{X}$ is a subspace of $\mathbb{R}^{d}$, $c$ is the squared Euclidean
distance and $a=\mathbf{1}_{n}$, we recover exactly the $k$-means algorithm.

######  Corollary 3.

Let $n\geq k\geq 1$ and
$\mathbf{X}\triangleq\\{x_{1},\dots,x_{n}\\}\subset\mathbb{R}^{d}$. We have
that

$\displaystyle\mathrm{LOT}_{k,\|\cdot\|_{2}^{2}}(\mu_{\mathbf{1}_{n},\mathbf{X}},\mu_{a,\mathbf{X}})=2\min_{Q,z_{1},\dots,z_{k}}\sum_{i=1}^{n}\sum_{q=1}^{k}Q_{i,q}\|x_{i}-z_{q}\|_{2}^{2}~{}~{}\text{s.t.}~{}~{}Q\in\\{0,1\\}^{n\times
k},~{}~{}Q\mathbf{1}_{k}=\mathbf{1}_{n}\;.$

In the general setting, solving
$\text{LOT}_{k,c}(\mu_{a,\mathbf{X}},\mu_{a,\mathbf{X}})$ for a given geometry
$c$, and a prescribed histrogram $a$ offers a new clustering method where the
assignment of the points to the clusters is determined by the matrix $Q^{*}$
solution of (6).

##  6 Computing LOT: Adaptive Stepsizes and Better Initializations

We target in this section practical issues that arises when using (Scetbon et
al., 2021, Algo.3) to solve (4). Scetbon et al. (2021) propose to apply a
mirror descent scheme with respect to the Kullback-Leibler divergence which
boils down to solve at each iteration $k\geq 0$ the following convex problem
using the Dykstra’s Algorithm (Dykstra, 1983):

$(Q_{k+1},R_{k+1},g_{k+1})\triangleq\\!\\!\operatorname*{argmin}_{\bm{\zeta}\in\mathcal{C}_{1}(a,b,r)\cap\mathcal{C}_{2}(r)}\\!\\!\text{KL}(\bm{\zeta},\bm{\xi}_{k})\,.$
(7)

where $(Q_{0},R_{0},g_{0})\in\mathcal{C}_{1}(a,b,r)\cap\mathcal{C}_{2}(r)$,
$\bm{\xi}_{k}\triangleq(\xi_{k}^{(1)},\xi_{k}^{(2)},\xi_{k}^{(3)})$,
$\xi_{k}^{(1)}\triangleq
Q_{k}\odot\exp(-\gamma_{k}CR_{k}\operatorname{diag}(1/g_{k}))$,
$\xi_{k}^{(2)}\triangleq
R_{k}\odot\exp(-\gamma_{k}C^{T}Q_{k}\operatorname{diag}(1/g_{k}))$,
$\xi_{k}^{(3)}\triangleq g_{k}\odot\exp(\gamma_{k}\omega_{k}/g_{k}^{2})$ with
$[\omega_{k}]_{i}\triangleq[Q_{k}^{T}CR_{k}]_{i,i}$ for all
$i\in\\{1,\dots,r\\}$,
$\mathrm{KL}(\mathbf{w},\mathbf{r})\triangleq\sum_{i}w_{i}\log(w_{i}/r_{i})$
and $(\gamma_{k})_{k\geq 0}$ is a sequence of positive step sizes. In the
general setting, each iteration of their algorithm requires $\mathcal{O}(nmr)$
operations and when the ground cost matrix $C$ admits a low-rank factorization
of the form $C=AB^{T}$ where $A\in\mathbb{R}^{n\times q}$ and
$B\in\mathbb{R}^{m\times q}$ with $q\ll\min(n,m)$, then the total complexity
per iteration becomes linear $\mathcal{O}((n+m)rq)$. Note that for the squared
Euclidean cost on $\mathbb{R}^{d}$, we have that $q=d+2$. In the following we
investigate two practical aspects of the algorithm: the choice of the step
sizes and the initialization.

#### Adaptive choice of $\gamma_{k}$.

Scetbon et al. (2021) show experimentally that the choice of
$(\gamma_{k})_{k\geq 0}$ does not impact the solution obtained upon
convergence, but rather the speed at which it is attained. Indeed the larger
$\gamma_{k}$ is, the faster the algorithm will converge. As a result, their
algorithm simply relies on a fixed $\gamma$ schedule. However, the range of
admissible $\gamma$ depends on the problem considered and it may vary from one
problem to another. Indeed, the algorithm might fail to converge as one needs
to ensure at each iteration $k$ of the mirror descent scheme that the kernels
$\bm{\xi}_{k}$ do not admit 0 entries in order to solve (7) using the
Dykstra’s Algorithm. Such a situation can occur when the terms involved in the
exponentials become too large which may depend on the problem considered.
Therefore, it may be of particular interest for practitioners to have a
generic range of admissible values for $\gamma$ independently of the
considered problem, in order to alleviate parameter tuning issues. We propose
to consider instead an adaptive choice of $(\gamma_{k})_{k\geq 0}$ along
iterations. D’Orazio et al. (2021); Bayandina et al. (2018) have proposed
adaptive mirror descent schemes where, at each iteration, the step-size is
normalized by the squared dual-norm of the gradient. Applying such a strategy
in our case amounts to consider at each iteration

$\displaystyle\gamma_{k}=\frac{\gamma}{\|\left(CR\operatorname{diag}(1/g),C^{T}Q\operatorname{diag}(1/g),-\mathcal{D}(Q^{T}RC)/g^{2}\right)\|_{\infty}^{2}},$
(8)

where the initial $\gamma>0$ is fixed. By doing so, we are able to guarantee a
lower-bound of the exponential terms involved in the expression of the kernels
$\bm{\xi}_{k}$ at each iteration and prevent them from having 0 entries. We
recommend to set such as global $\gamma\in[1,10]$, and observe that this range
works whatever the problem considered.

#### On the choice of the initialization.

As $\text{LOT}_{r,c}$ (4) is a non-convex optimization problem, the question
of choosing an efficient initialization arises in practice. Scetbon et al.
(2021) show experimentally that the convergence of the algorithm does not
depend on the initalization chosen if no stopping criterion is used. Indeed,
their experimental findings support that only well behaved local minimas are
attractive. However, in practice one needs to use a stopping criterion in
order to terminate the algorithm. We do observe in many instances that using
trivial initializers may result in spurious local minima, which trigger the
stopping criterion early on and prevent the algorithm to reach a good
solution. Based on various experimentations, we propose to consider a novel
initialization of the algorithm. Our initialization aims at being close to a
well-behaved local minimum by clustering the input measures. When the measures
are supported on Euclidean space, we propose to find $r$ centroids
$(z_{i})_{i=1}^{r}$ of one of the two input discrete probability measures
using $k$-means and to solve the following convex barycenter problem:

$\displaystyle\min_{Q,R}\langle C_{X,Z},Q\rangle+\langle
C_{Y,Z},R\rangle-\varepsilon H(Q)-\varepsilon
H(R)~{}~{}\text{s.t.}~{}~{}Q\mathbf{1}_{n}=a,~{}~{}R\mathbf{1}_{n}=b,~{}~{}Q^{T}\mathbf{1}_{r}=R^{T}\mathbf{1}_{r}\,,$
(9)

where $C_{X,Z}=(c(x_{i},z_{j}))_{i,j}$, $C_{Y,Z}=(c(y_{i},z_{j}))_{i,j}$, and
$H(P)=-\sum_{i,j}P_{i,j}(\log(P_{i,j}-1)$. In practice we fix
$\varepsilon=1/10$ and we then initialize $\text{LOT}_{r,c}$ using $(Q,R)$
solution of (9) and $g\triangleq Q^{T}\mathbf{1}_{r}(=R^{T}\mathbf{1}_{r})$.
Note that $(Q,R,g)$ is an admissible initialization and finding the centroids
as well as solving (9) requires $\mathcal{O}((n+m)r)$ algebraic operations.
Therefore such initialization does not change the total complexity of the
algorithm.

In the general (non-Euclidean) case, we propose to initialize the algorithm by
applying our generalized $k$-means approach defined in (6) on each input
measure where we fix the common marginal to be $g=\mathbf{1}_{r}/r$. More
precisely, by denoting $C_{X,X}=(c(x_{i},x_{j}))_{i,j}$ and
$C_{Y,Y}=(c(y_{i},y_{j}))_{i,j}$, we initialize the algorithm by solving:

$\displaystyle\begin{aligned} Q&\in\operatorname*{argmin}_{Q}\langle
C_{X,X},Q\text{diag}(1/Q^{T}\mathbf{1}_{n})Q^{T}\rangle~{}~{}\text{s.t.}~{}~{}Q\in\mathbb{R}^{n\times
k}_{+}~{}~{}\text{,}~{}~{}Q\mathbf{1}_{k}=a,~{}~{}Q^{T}\mathbf{1}_{n}=\mathbf{1}_{r}/r\,.\\\
R&\in\operatorname*{argmin}_{R}\langle
C_{Y,Y},R\text{diag}(1/R^{T}\mathbf{1}_{m})R^{T}\rangle~{}~{}\text{s.t.}~{}~{}R\in\mathbb{R}^{m\times
k}_{+}~{}~{}\text{,}~{}~{}R\mathbf{1}_{k}=b,~{}~{}R^{T}\mathbf{1}_{n}=\mathbf{1}_{r}/r\;.\end{aligned}$
(10)

Note that again the $(Q,R,g)$ obtained is an admissible initialization and the
complexity of solving (10) is of the same order as solving (4), thus the total
complexity of the algorithm remains the same.

##  7 Experiments

Figure 1: In this experiment, we consider a mixture of 10 anisotropic
Gaussians supported on $\mathbb{R}^{d}$ and we plot the value of
$\text{DLOT}_{r,c}$ between two independent empirical measures associated to
this mixture when varying the number of samples $n$ and the dimension $d$ for
multiple ranks $r$. The ground cost considered is the squared Euclidean
distance. Note that $\text{LOT}_{r}(\mu,\mu)\neq 0$ and therefore we use
$\text{DLOT}_{r,c}(\mu,\mu)$ instead to evaluate the rates. Each point has
been obtained by repeating 10 times the experiment. We compare the empirical
rates obtained with the theoretical one derived in Proposition 4 for $r=1$. We
observe that our theoretical results match the empirical ones and, as
expected, the rates do not depend on $d$.

In this section, we illustrate experimentally our theoretical findings and
show how our initialization provide practical improvements. For that purpose
we consider 3 synthetic problems and one real world dataset to: (i) provide
illustrations on the statistical rates of $\text{LOT}_{r,c}$, (ii) exhibit the
gradient flow of the debiased formulation $\text{DLOT}_{r,c}$, (iii) use the
clustering method induced by $\text{LOT}_{r,c}$, and (iv) show the effect of
the initialization. All experiments were run on a MacBook Pro 2019 laptop.

Statistical rates. We aim at showing the statistical rates of the plug-in
estimator of $\text{LOT}_{r,c}$. As $\text{LOT}_{r,c}(\mu,\mu)\neq 0$ and as
we do not have access to this value given samples from $\mu$, we consider
instead the debiased version of the low-rank optimal transport,
$\text{DLOT}_{r,c}$. In figure 1, we show that the empiricial rates match the
theoretical bound obtained in Proposition 4. In particular, we show that that
these rates does not depend on the dimension of the ground space. Note also
that we recover our theoretical dependence with respect to the rank $r$: the
higher the rank, the slower the convergence.

Figure 2: We compare the gradient flows $(\mu_{t})_{t\geq 0}$ (in red)
starting from a Gaussian distribution, $\mu_{0}$, to a moon shape distribution
(in blue), $\nu$, in 2D when minimizing either
$L(\mu)\triangleq\text{DLOT}_{r,c}(\mu,\nu)$ or
$L(\mu)\triangleq\text{LOT}_{r,c}(\mu,\nu)$. The ground cost is the squared
Euclidean distance and we fix $r=100$. We consider 1000 samples from each
distribution and and we plot the evolution of the probability measure obtained
along the iterations of a gradient descent scheme. We also display in green
the vector field in the descent direction. We show that the debiased version
allows to recover the target distribution while $\text{LOT}_{r,c}$ is learning
a biased version with a low-rank structure.

Gradient Flows using DLOT. We illustrate here a practical use of DLOT for ML
application. In figure 6, we consider $Y_{1},\dots,Y_{n}$ independent samples
from a moon shape distribution in 2D, and by denoting $\hat{\nu}_{n}$ the
empirical measure associated, we show the iterations obtained by a gradient
descent scheme on the following optimization problem:

$\displaystyle\min_{\mathbf{X}\in\mathbb{R}^{n\times
2}}\text{DLOT}_{r,c}(\mu_{\mathbf{1}_{n}/n,\mathbf{X}},\hat{\nu}_{n})\;.$

We initialize the algorithm using $n=1000$ samples drawn from a Gaussian
distribution. We show that the gradient flow of our debiased version is able
to recover the target distribution. We also compare it with the gradient flow
of the biased version (LOT) and show that it fails to reproduce the target
distribution as it is learning a biased one with a low-rank structure.

Figure 3: In this experiment, we draw 1000 samples from multiple distributions
from the python package scikit-learn (Pedregosa et al., 2011) and we apply the
method proposed in (6) for two different costs: in the top row we consider the
squared Euclidean distance while in the bottom row, we consider the shortest
path distance on the graph associated with the ground cost $c(x,y)=1-k(x,y)$
where $k$ is a Gaussian kernel. In the two first problem (starting from the
left), we fix $r=2$, in the next three problem we fix $r=3$ and in the last
one we fix $r=4$. We observe that the flexibility of our method allows to
recover the clustering for a well chosen ground cost.

Application to Clustering. In this experiment we show some applications of the
clustering method induced by $\text{LOT}_{r,c}$. In figure 3, we consider 6
datasets with different structure and we aim at recovering the clusters using
(6) for some well chosen costs. We compare the clusters obtained when
considering either the squared Euclidean cost (which amounts at applying the
$k$-means) and the shortest-path distance on the data viewed as a graph. We
show that our method is able to recover the clusters on these settings for
well chosen costs and therefore the proposed algorithm in Scetbon et al.
(2021) can be seen as a new alternative in order to clusterize data.

Figure 4: In this experiment, we consider the Newsgroup20 dataset (Pedregosa
et al., 2011) constituted of texts and we embed them into distributions in 50D
using the same pre-processing steps as in (Cuturi et al., 2022). We compare
different initialization when applying the algorithm of (Scetbon et al., 2021)
to compare random texts viewed as distributions for multiple choices of rank
$r$. The ground cost considered in the squared Euclidean distance. We repeat
the experiments 50 times by sampling randomly multiple problems of similar
size ($\simeq 250$ samples). We normalize the cost matrix by its maximum value
in order to have comparable LOT cost. We consider 4 different initialization:
the one using $k$-means algorithm (9), the one using the generalized $k$-means
(10), the rank-2 initialization (Scetbon et al., 2021) and a random
initialization where $Q,R$ and $g$ are drawn from Gaussians. We compare both
the cost value and the criterion value ($\Delta_{k}$) along the iterations of
the MD scheme. Note that the curves obtained do not start at the same point in
time as we start plotting the curves after obtaining the initial point which
in some case requires more algebraic operations (e.g. kmeans methods). First
we observe that whatever the initialization considered, the algorithm
converges toward the same value. In addition, we observe that both $k$-means
and general $k$-means are able to initialize well the algorithm by avoiding
bad local minima at initialization while the two other initialization are
close to spurious local minima at initialization.

Effect of the Initialization. Our goal here is to show the effect of the
initialization. In figure 4, we display the evolution of the cost as well as
the value of the stopping criterion along the iterations of the MD scheme
solving (4) when considering different initialization. The $x$-axis
corresponds to the total number of algebraic operations. This number is
computed at each iteration of the outer loop of the algorithm proposed in
Scetbon et al. (2021) and is obtained by computing the complexity of all the
operations involved in their algorithm to reach it. We consider this notion of
time instead of CPU/GPU time as we do not want to be architecture/machine
dependent. Recall also that the stopping criterion introduced in (Scetbon et
al., 2021) is defined for all $k\geq 1$ by

$\displaystyle\Delta_{k}\triangleq\frac{1}{\gamma_{k}^{2}}(\mathrm{KL}((Q_{k},R_{k},g_{k}),(Q_{k-1},R_{k-1},g_{k-1}))+\mathrm{KL}((Q_{k-1},R_{k-1},g_{k-1}),(Q_{k},R_{k},g_{k}))),$

where $((Q_{k},R_{k},g_{k}))_{k\geq 0}$ is the sequence solution of (7).
First, we show that whatever the initialization chosen, the algorithm manages
to converge to an efficient solution if no stopping criterion is used.
However, the choice of the initialization may impact the termination of the
algorithm as some initialization might be too close to some spurious local
minima. Indeed, the initial points obtained using a “rank 2” or random
initialization can be close to spurious and non-attractive local minima, which
may trigger the stopping criterion too early and prevent the algorithm from
continuing to run in order to converge towards an attractive and well behaved
local minimum. We show also that the initialization we propose in (9) and (10)
are sufficiently far away from bad local minima and allow the algorithm to
converge directly toward the desired solution.

The right figure of Fig.4 shows two main observations: (i) that the initial
point obtained using a “rank 2” or random initialization can be close to
spurious and non-attractive local minima, which may trigger the stopping
criterion too early and prevent the algorithm from continuing to run in order
to converge towards an attractive and well behaved local minimum. (ii) When
initialiazing the algorithm using kmeans methods, we show that our stopping
criterion is a decreasing function of time meaning that the algorithm
converges directly towards the desired solution.

#### Conclusion.

We assembled in this work theoretical and practical arguments to support low-
rank factorizations for OT. We have presented two controls: one concerning the
approximation error to the true optimal transport and another concerning the
statistical rates of the plug-in estimator. The latter is showed to be
independent of the dimension, which is of particular interest when studying OT
in ML settings. We have motivated further the use of LOT as a loss by
introducing its debiased version and showed that it possesses desirable
properties: positivity and metrization of the convergence in law. We have also
presented the links between the bias induced by such regularization and
clustering methods, and studied empirically the effects of hyperparameters
involved in the practical estimation of LOT. The strong theoretical
foundations provided in this paper motivate further studies of the empirical
behaviour of LOT estimator, notably on finding suitable local minima and on
improvements on the convergence of the MD scheme using other adaptive choices
for step sizes.

Acknowledgements. This work was supported by a "Chaire d’excellence de l’IDEX
Paris Saclay". The authors would also like to thank Gabriel Peyré and Jaouad
Mourtada for enlightening conversations on the topics discussed in this work.

## References

  * Bayandina et al. (2018) Anastasia Bayandina, Pavel Dvurechensky, Alexander Gasnikov, Fedor Stonyakin, and Alexander Titov.  Mirror descent and convex optimization problems with non-smooth inequality constraints.  In _Large-scale and distributed optimization_ , pages 181–213. Springer, 2018.
  * Chizat et al. (2020) Lenaic Chizat, Pierre Roussillon, Flavien Léger, François-Xavier Vialard, and Gabriel Peyré.  Faster wasserstein distance estimation with the sinkhorn divergence.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Clason et al. (2021) Christian Clason, Dirk A. Lorenz, Hinrich Mahler, and Benedikt Wirth.  Entropic regularization of continuous optimal transport problems.  _Journal of Mathematical Analysis and Applications_ , 494(1):124432, 2021.  ISSN 0022-247X.  doi: https://doi.org/10.1016/j.jmaa.2020.124432.
  * Cuturi (2013) Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Advances in neural information processing systems_ , pages 2292–2300, 2013.
  * Cuturi et al. (2022) Marco Cuturi, Laetitia Meng-Papaxanthos, Yingtao Tian, Charlotte Bunne, Geoff Davis, and Olivier Teboul.  Optimal transport tools (ott): A jax toolbox for all things wasserstein.  _arXiv preprint arXiv:2201.12324_ , 2022.
  * Demetci et al. (2020) Pinar Demetci, Rebecca Santorella, Björn Sandstede, William Stafford Noble, and Ritambhara Singh.  Gromov-wasserstein optimal transport to align single-cell multi-omics data.  _BioRxiv_ , 2020.
  * D’Orazio et al. (2021) Ryan D’Orazio, Nicolas Loizou, Issam Laradji, and Ioannis Mitliagkas.  Stochastic mirror descent: Convergence analysis and adaptive variants via the mirror stochastic polyak stepsize.  _arXiv preprint arXiv:2110.15412_ , 2021.
  * Dudley (1969) Richard Mansfield Dudley.  The speed of mean glivenko-cantelli convergence.  _The Annals of Mathematical Statistics_ , 40(1):40–50, 1969.
  * Dykstra (1983) Richard L Dykstra.  An algorithm for restricted least squares regression.  _Journal of the American Statistical Association_ , 78(384):837–842, 1983.
  * Feydy et al. (2018) Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-Ichi Amari, Alain Trouvé, and Gabriel Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  _arXiv preprint arXiv:1810.08278_ , 2018.
  * Forrow et al. (2019) Aden Forrow, Jan-Christian Hütter, Mor Nitzan, Philippe Rigollet, Geoffrey Schiebinger, and Jonathan Weed.  Statistical optimal transport via factored couplings.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 2454–2465. PMLR, 2019.
  * Fournier and Guillin (2015) Nicolas Fournier and Arnaud Guillin.  On the rate of convergence in Wasserstein distance of the empirical measure.  _Probability Theory and Related Fields_ , 162(3-4):707–738, 2015.
  * Genevay et al. (2018a) Aude Genevay, Lénaic Chizat, Francis Bach, Marco Cuturi, and Gabriel Peyré.  Sample complexity of sinkhorn divergences.  _arXiv preprint arXiv:1810.02733_ , 2018a.
  * Genevay et al. (2018b) Aude Genevay, Gabriel Peyre, and Marco Cuturi.  Learning generative models with sinkhorn divergences.  In _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_ , Proceedings of Machine Learning Research. PMLR, 09–11 Apr 2018b.
  * Gretton et al. (2012) Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(Mar):723–773, 2012.
  * Heitz et al. (2020) Matthieu Heitz, Nicolas Bonneel, David Coeurjolly, Marco Cuturi, and Gabriel Peyré.  Ground metric learning on graphs.  _Journal of Mathematical Imaging and Vision_ , pages 1–19, 2020.
  * Janati et al. (2020) Hicham Janati, Thomas Bazeille, Bertrand Thirion, Marco Cuturi, and Alexandre Gramfort.  Multi-subject meg/eeg source imaging with sparse multi-task regression.  _NeuroImage_ , page 116847, 2020.
  * Kantorovich (1942) Leonid Kantorovich.  On the transfer of masses (in russian).  _Doklady Akademii Nauk_ , 37(2):227–229, 1942\.
  * Koundal et al. (2020) Sunil Koundal, Rena Elkin, Saad Nadeem, Yuechuan Xue, Stefan Constantinou, Simon Sanggaard, Xiaodan Liu, Brittany Monte, Feng Xu, William Van Nostrand, et al.  Optimal mass transport with lagrangian workflow reveals advective and diffusion driven solute transport in the glymphatic system.  _Scientific reports_ , 10(1):1–18, 2020.
  * Kutin (2002) Samuel Kutin.  Extensions to mcdiarmid’s inequality when differences are bounded with high probability.  _Dept. Comput. Sci., Univ. Chicago, Chicago, IL, USA, Tech. Rep. TR-2002-04_ , 2002.
  * Liu et al. (2021) Weijie Liu, Chao Zhang, Nenggan Zheng, and Hui Qian.  Approximating optimal transport via low-rank and sparse factorization.  _arXiv preprint arXiv:2111.06546_ , 2021.
  * Mena and Niles-Weed (2019) Gonzalo Mena and Jonathan Niles-Weed.  Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Santambrogio (2015) Filippo Santambrogio.  Optimal transport for applied mathematicians.  _Birkäuser, NY_ , 55(58-63):94, 2015.
  * Scetbon et al. (2021) Meyer Scetbon, Marco Cuturi, and Gabriel Peyré.  Low-rank sinkhorn factorization, 2021.
  * Scetbon et al. (2022) Meyer Scetbon, Gabriel Peyré, and Marco Cuturi.  Linear-time gromov wasserstein distances using low rank couplings and costs.  _ICML_ , 2022.
  * Schiebinger et al. (2019) Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, Vidya Subramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin, Peter Berube, et al.  Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming.  _Cell_ , 176(4):928–943, 2019.
  * Schmitz et al. (2018) Morgan A Schmitz, Matthieu Heitz, Nicolas Bonneel, Fred Ngole, David Coeurjolly, Marco Cuturi, Gabriel Peyré, and Jean-Luc Starck.  Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning.  _SIAM Journal on Imaging Sciences_ , 11(1):643–678, 2018.
  * Sejdinovic et al. (2013) Dino Sejdinovic, Bharath Sriperumbudur, Arthur Gretton, and Kenji Fukumizu.  Equivalence of distance-based and rkhs-based statistics in hypothesis testing.  _The Annals of Statistics_ , pages 2263–2291, 2013.
  * Yang et al. (2020) Karren Dai Yang, Karthik Damodaran, Saradha Venkatachalapathy, Ali C Soylemezoglu, GV Shivashankar, and Caroline Uhler.  Predicting cell lineages using autoencoders and optimal transport.  _PLoS computational biology_ , 16(4):e1007828, 2020.
  * Zheng et al. (2020) Xinye Zheng, Jianbo Ye, James Z Wang, and Jia Li.  Scott: Shape-location combined tracking with optimal transport.  _SIAM Journal on Mathematics of Data Science_ , 2(2):284–308, 2020.

## Supplementary materials

##  Appendix A On the Definition of $\text{LOT}_{r,c}$

Let $(\mathcal{X},d_{\mathcal{X})}$ and $(\mathcal{Y},d_{\mathcal{Y}})$ two
nonempty compact Polish spaces, $\mu\in\mathcal{M}_{1}^{+}(\mathcal{X})$,
$\nu\in\mathcal{M}_{1}^{+}(\mathcal{Y})$ two probability measures on these
spaces and $c:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}_{+}$ a
nonnegative and continuous function. We define the generalized low-rank
optimal transport between $\mu$ and $\nu$ as

$\displaystyle\text{LOT}_{r,c}(\mu,\nu)\triangleq\inf_{\pi\in\Pi_{r}(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi(x,y)\;.$

where

$\Pi_{r}(\mu,\nu)\triangleq\\{\pi\in\Pi(\mu,\nu):\exists(\mu_{i})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{X})^{r},~{}(\nu_{i})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{Y})^{r},~{}\lambda\in\Delta_{r}^{*}\text{~{}s.t.~{}}\pi=\sum_{i=1}^{r}\lambda_{i}\mu_{i}\otimes\nu_{i}\\}\;.$

As $\mathcal{X}$ and $\mathcal{Y}$ are compact, $\Pi_{r}(\mu,\nu)$ is tight,
then Prokhorov’s theorem applies and the closure of $\Pi_{r}(\mu,\nu)$ is
sequentially compact. Let us now show that $\Pi_{r}(\mu,\nu)$ is closed.
Indeed, Let $(\pi_{n})_{n\geq 0}$ a sequence of $\Pi_{r}(\mu,\nu)$ converging
towards $\pi_{*}$. Then by definition there exists for all $k\in[|1,r|]$,
$(\mu_{n}^{(k)})_{n\geq 0}$, $(\nu_{n}^{(k)})_{n\geq 0}$ and
$(\lambda_{n}^{(k)})_{n\geq 0}$ such that for all $n\geq 0$

$\displaystyle\pi_{n}=\sum_{i=1}^{r}\lambda_{n}^{(k)}\mu_{n}^{(k)}\otimes\nu_{n}^{(k)}\;.$

However, $(\mu_{n}^{(k)})_{n\geq 0}$ and $(\nu_{n}^{(k)})_{n\geq 0}$ are also
tight, and Prokhorov’s theorem applies, therefore we can extract a common
subsequence such that for all $k$,

$\displaystyle\mu_{n}^{(k)}\rightarrow\mu_{*}^{(k)}~{}~{}\text{and
}\nu_{n}^{(k)}\rightarrow\nu_{*}^{(k)}$

In addition as $(\lambda_{n})_{n\geq 0}$ live in the simplex $\Delta_{r}$, we
can also extract a sub-sequence, such that
$\lambda_{n}\rightarrow\lambda_{*}\in\Delta_{r}$. Finally by unicity of the
limit we obtain that

$\displaystyle\pi_{*}=\sum_{k=1}^{r}\lambda_{*}^{(k)}\mu_{*}^{(k)}\otimes\nu_{*}^{(k)}\;.$

Finally, by denoting $I\triangleq\\{k:\lambda_{*}^{(k)}>0\\}$, and by
considering $i^{*}\in I$, we obtain that

$\displaystyle\pi_{*}=\sum_{i\in
I\backslash\\{{i^{*}}\\}}^{r}\lambda_{*}^{(i)}\mu_{*}^{(i)}\otimes\nu_{*}^{(i)}+\sum_{j=1}^{r-|I|+1}\frac{\lambda_{*}^{(i^{*})}}{r-|I|+1}\mu_{*}^{(i^{*})}\otimes\nu_{*}^{(i^{*})}\;.$

from which follows that $\pi_{*}\in\Pi_{r}(\mu,\nu)$.

##  Appendix B Proofs

###  B.1 Proof of Proposition 1

######  Proposition.

Let $n,m\geq 2$,
$\mathbf{X}\triangleq\\{x_{1},\dots,x_{n}\\}\subset\mathcal{X}$,
$\mathbf{Y}\triangleq\\{y_{1},\dots,y_{m}\\}\subset\mathcal{Y}$ and
$a\in\Delta_{n}^{*}$ and $b\in\Delta_{m}^{*}$. Then for $2\leq
r\leq\min(n,m)$, we have that

$\displaystyle|\mathrm{LOT}_{r,c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})-\mathrm{OT}_{c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})|\leq\|C\|_{\infty}~{}\text{ln}(\min(n,m)/(r-1))$

###### Proof.

Let $P\in\operatorname*{argmin}_{P\in\Pi_{a,b}}\langle C,P\rangle$. As $P$ is
a nonnegative matrix, its nonnegative rank cannot exceed $\min(n,m)$. Assume
for simplicity, that $n=m$, then there exists $(R_{i})_{i=1}^{n}$ nonnegative
matrices of rank 1 such that

$\displaystyle P=\sum_{i=1}^{n}R_{i}\;.$

As for all $i\in[|1,n|]$, $R_{i}$ is a rank 1 matrix, there exist
$\tilde{q}_{i},\tilde{r}_{i}\in\mathbb{R}_{+}^{n}$ such that
$R_{i}=\tilde{q}_{i}\tilde{r}_{i}^{T}$. Then by denoting
$q_{i}=\tilde{q}_{i}/|\tilde{q}_{i}|$, $r_{i}=\tilde{r}_{i}/|\tilde{r}_{i}|$
and $\lambda_{i}=|\tilde{q}_{i}||\tilde{r}_{i}|$ where for any
$h\in\mathbb{R}^{n}$ $|h|\triangleq\sum_{i=1}^{n}h_{i}$, we obtain that

$P=\sum_{i=1}^{n}\lambda_{i}q_{i}r_{i}^{T}\;.$

Without loss of generality, we can consider the case where
$\lambda_{1}\geq\dots\geq\lambda_{n}$. Let us now denote
$\mathbf{\lambda}:=(\lambda_{1},\dots,\lambda_{n})$, and by using the fact the
$P$ is a coupling we obtain that $\mathbf{\lambda}\in\Delta_{n}$. Also, by
definition of $\lambda$, we have that for all $k\in[|1,n|]$, $\lambda_{k}\leq
1/k$. Let us now define

$\displaystyle\tilde{P}\triangleq\sum_{i=1}^{r-1}\lambda_{i}q_{i}r_{i}^{T}+\left(\sum_{i=r}^{n}\lambda_{i}\right)\alpha_{r}\beta_{r}^{T}$

where

$\displaystyle\alpha_{r}$
$\displaystyle\triangleq\frac{\sum_{i=r}^{n}\lambda_{i}q_{i}}{\sum_{i=r}^{n}\lambda_{i}}$
$\displaystyle\beta_{r}$
$\displaystyle\triangleq\frac{\sum_{i=r}^{n}\lambda_{i}r_{i}}{\sum_{i=r}^{n}\lambda_{i}}$

Remark that $\tilde{P}\in\Pi_{a,b}(r)$, therefore we obtain that

$\displaystyle|\text{LOT}_{r,c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})-\text{OT}_{c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})|$
$\displaystyle=\text{LOT}_{r,c}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})-\text{OT}_{x}(\mu_{a,\mathbf{X}},\nu_{b,\mathbf{Y}})$
$\displaystyle\leq\langle C,\tilde{P}\rangle-\langle C,P\rangle$
$\displaystyle\leq\langle
C,\left(\sum_{i=r}^{n}\lambda_{i}\right)\alpha_{r}\beta_{r}^{T}\rangle-\langle
C,\sum_{i=r}^{n}\lambda_{i}q_{i}r_{i}^{T}\rangle$ $\displaystyle\leq\langle
C,\left(\sum_{i=r}^{n}\lambda_{i}\right)\alpha_{r}\beta_{r}^{T}\rangle$
$\displaystyle\leq\|C\|_{\infty}\sum_{i=r}^{n}\lambda_{i}\leq\|C\|_{\infty}\sum_{i=r}^{n}\frac{1}{i}\leq\|C\|_{\infty}~{}\text{ln}(n/(r-1))$

∎

###  B.2 Proof of Proposition 2

######  Proposition 10.

Let $\mu\in\mathcal{M}_{1}^{+}(\mathcal{X})$,
$\nu\in\mathcal{M}_{1}^{+}(\mathcal{Y})$ and let us assume that $c$ is
$L$-Lipschitz w.r.t. $x$ and $y$ . Then for any $r\geq 1$, we have

$\displaystyle|\text{LOT}_{r,c}(\mu,\nu)-\text{OT}_{c}(\mu,\nu)|\leq
2L\max(\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{X},d_{\mathcal{X}}),\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{Y},d_{\mathcal{Y}}))$

###### Proof.

As $\mathcal{X}$ and $\mathcal{Y}$ are compact,
$\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{X},d),\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{Y},d)<+\infty$
and then by denoting
$\varepsilon_{\mathcal{X}}\triangleq\mathcal{N}_{\lfloor\log_{2}(\lfloor\sqrt{r}\rfloor)\rfloor}(\mathcal{X},d_{\mathcal{X}})$,
there exists $x_{1},\dots,x_{\lfloor\sqrt{r}\rfloor}\in\mathcal{X}$, such that
$\mathcal{X}\subset\bigcup_{i=1}^{r}\mathcal{B}_{\mathcal{X}}(x_{i},\varepsilon)$
from which we can extract a partition
$(S_{i,\mathcal{X}})_{i=1}^{\lfloor\sqrt{r}\rfloor}$ of $\mathcal{X}$ such
that for all $i\in[|1,\lfloor\sqrt{r}\rfloor|]$, and $x,y\in
S_{i,\mathcal{X}}$, $d_{\mathcal{X}}(x,y)\leq\varepsilon_{\mathcal{X}}$.
Similarly we can build a partition
$(S_{i,\mathcal{Y}})_{i=1}^{\lfloor\sqrt{r}\rfloor}$ of $\mathcal{Y}$. Let us
now define for all $k\in[|1,\lfloor\sqrt{r}\rfloor|]$,

$\displaystyle\mu_{k}\triangleq\frac{\mu|_{S_{k,\mathcal{X}}}}{\mu(S_{k,\mathcal{X}})}~{}~{}\text{and}~{}~{}\nu_{k}\triangleq\frac{\nu|_{S_{k,\mathcal{Y}}}}{\nu(S_{k,\mathcal{Y}})}$

with the convention that $\frac{0}{0}=0$, we can define

$\displaystyle\pi_{r}\triangleq\sum_{i,j=1}^{\lfloor\sqrt{r}\rfloor}\pi^{*}(S_{i,\mathcal{X}}\times
S_{j,\mathcal{Y}})\nu_{j}\otimes\mu_{i}\;.$

First remarks that $\pi_{r}\in\Pi_{r}(\mu,\nu)$. Indeed we have for any
measurable set $B$

$\displaystyle\pi_{r}(\mathcal{X}\times B)$
$\displaystyle=\sum_{j=1}^{\lfloor\sqrt{r}\rfloor^{2}}\nu_{j}(B)\sum_{i=1}^{r}\pi^{*}(S_{i,\mathcal{X}}\times
S_{j,\mathcal{Y}})$
$\displaystyle=\sum_{j=1}^{\lfloor\sqrt{r}\rfloor}\nu_{j}(B)\nu(S_{j,\mathcal{Y}})$
$\displaystyle=\sum_{j=1}^{\lfloor\sqrt{r}\rfloor}\nu|_{S_{j,\mathcal{X}}}(B)$
$\displaystyle=\nu(B)\;,$

similarly $\pi_{r}(A\times\mathcal{Y})=\mu(A)$ and we have that
$\lfloor\sqrt{r}\rfloor^{2}\leq r$. Therefore we obtain that

$\displaystyle|\text{LOT}_{r,c}(\mu,\nu)-\text{OT}_{c}(\mu,\nu)|$
$\displaystyle=\text{LOT}_{r,c}(\mu,\nu)-\text{OT}_{c}(\mu,\nu)$
$\displaystyle\leq\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi_{r}(x,y)-\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi^{*}(x,y)$
$\displaystyle\leq\sum_{i,j=1}^{\lfloor\sqrt{r}\rfloor}\int_{S_{i,\mathcal{X}}\times
S_{j,\mathcal{Y}}}c(x,y)d[\pi_{r}(x,y)-\pi^{*}(x,y)]$
$\displaystyle\leq\sum_{i,j=1}^{\lfloor\sqrt{r}\rfloor}\pi^{*}(S_{i,\mathcal{X}}\times
S_{j,\mathcal{Y}})$ $\displaystyle\times[\sup_{(x,y)\in
S_{i,\mathcal{X}}\times S_{j,\mathcal{Y}}}c(x,y)-\inf_{(x,y)\in
S_{i,\mathcal{X}}\times S_{j,\mathcal{Y}}}c(x,y)]$ $\displaystyle\leq
L[\varepsilon_{\mathcal{X}}+\varepsilon_{\mathcal{Y}}]$

from which the result follows. ∎

######  Corollary.

Under the same assumptions of Proposition 2 and by assuming in addition that
there exists a Monge map solving $\text{OT}_{c}(\mu,\nu)$, we obtain that for
any $r\geq 1$,

$\displaystyle|\mathrm{LOT}_{r,c}(\mu,\nu)-\mathrm{OT}_{c}(\mu,\nu)|\leq
L\mathcal{N}_{\lfloor\log_{2}(r)\rfloor}(\mathcal{Y},d_{\mathcal{Y}})$

###### Proof.

Let us denote $T$ a Monge map solution of $\text{OT}_{c}(\mu,\nu)$ and as in
the proof above, let us consider a partition of
$(S_{i,\mathcal{Y}})_{i=1}^{r}$ of $\mathcal{Y}$ such that for all
$i\in[|1,r|]$, and $x,y\in S_{i,\mathcal{Y}}$,
$d_{\mathcal{Y}}(x,y)\leq\varepsilon_{\mathcal{Y}}$ with
$\varepsilon_{\mathcal{Y}}\triangleq\mathcal{N}_{\lfloor\log_{2}(r)\rfloor}(\mathcal{Y},d_{\mathcal{Y}})$.
Let us now define for all $k\in[|1,\lfloor\sqrt{r}\rfloor|]$,

$\displaystyle\mu_{k}\triangleq\frac{\mu|_{T^{-1}(S_{k,\mathcal{Y}})}}{\mu(T^{-1}(S_{k,\mathcal{Y}}))}~{}~{}\text{and}~{}~{}\nu_{k}\triangleq\frac{\nu|_{S_{k,\mathcal{Y}}}}{\nu(S_{k,\mathcal{Y}})}$

with the convention that $\frac{0}{0}=0$, we can define

$\displaystyle\pi_{r}\triangleq\sum_{k=1}^{r}\pi^{*}(T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}})\nu_{k}\otimes\mu_{k}\;.$

Again we have that $\pi_{r}\in\Pi_{r}(\mu,\nu)$, and we obtain that

$\displaystyle|\text{LOT}_{r,c}(\mu,\nu)-\text{OT}_{c}(\mu,\nu)|$
$\displaystyle=\text{LOT}_{r,c}(\mu,\nu)-\text{OT}_{c}(\mu,\nu)$
$\displaystyle\leq\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi_{r}(x,y)-\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)d\pi^{*}(x,y)$
$\displaystyle\leq\sum_{k=1}^{r}\pi^{*}(T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}})\int_{T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}}}c(x,y)d\mu_{k}(y)\otimes\nu_{k}(y)$
$\displaystyle-\sum_{k=1}^{r}\int_{T^{-1}(S_{k,\mathcal{Y}})}c(x,T(x))d\mu(x)$
$\displaystyle\leq\sum_{k=1}^{r}\pi^{*}(T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}})\int_{T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}}}c(x,y)d\mu_{k}(y)\otimes\nu_{k}(y)$
$\displaystyle-\sum_{k=1}^{r}\pi^{*}(T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}})\int_{T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}}}c(x,T(x))d\mu_{k}(x)\otimes\nu_{k}(y)$
$\displaystyle\leq\sum_{k=1}^{r}\pi^{*}(T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}})\int_{T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}}}[c(x,y)-c(x,T(x))]d\mu_{k}(x)\otimes\nu_{k}(x)$
$\displaystyle\leq L\varepsilon_{\mathcal{Y}}$

from which the result follows. Note that to obtain the above inequalities, we
use the fact that $\pi^{*}$ is supported on the graph of $T$, and therefore we
have have for all $k\in[|1,r|]$,

$\displaystyle\pi^{*}(T^{-1}(S_{k,\mathcal{Y}})\times
S_{k,\mathcal{Y}})=\mu(T^{-1}(S_{k,\mathcal{Y}}))=\nu(S_{k,\mathcal{Y}}).$

∎

###  B.3 Proof of Proposition 3

######  Proposition.

Let $r\geq 1$ and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$, then
$\mathrm{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})\xrightarrow[n\rightarrow+\infty]{}\mathrm{LOT}_{r,c}(\mu,\nu)~{}~{}a.s.$

###### Proof.

Let $\pi^{*}$ solution of $\text{LOT}_{r,c}(\mu,\nu)$. Then there exists
$\lambda^{*}\in\Delta_{r}^{*}$,
$(\mu_{i}^{*})_{i=1}^{r},(\nu_{i}^{*})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{X})^{r}$
such that

$\pi^{*}=\sum_{i=1}^{r}\lambda_{i}^{*}\mu_{i}^{*}\otimes\nu_{i}^{*}.$

Note that by definition, we have that

$\displaystyle\mu=\sum_{i=1}^{r}\lambda_{i}^{*}\mu_{i}^{*}\text{~{}~{}and~{}~{}}\nu=\sum_{i=1}^{r}\lambda_{i}^{*}\nu_{i}^{*}\;.$

Let us now define $\pi_{\mu}$ and $\pi_{\mu}$ both elements of
$\mathcal{M}_{1}^{+}(\mathcal{X}\times[|1,r|])$ as follows:

$\displaystyle\pi_{\mu}(A\times\\{k\\})\triangleq\lambda_{k}\mu_{k}(A)\text{~{}~{}and~{}~{}}\pi_{\nu}(A\times\\{k\\})\triangleq\lambda_{k}\nu_{k}(A)~{}~{}\text{for
any measurable set $A$ and $k\in[|1,r|]$}\;.$

Observe that the right marginals of $\pi_{\mu}$ and $\pi_{\nu}$ is the same
and we will denote it $\rho$. We can now define for all $x,y\in\mathcal{X}$
the family of kernels
$(k_{\mu}(\cdot,x))_{x\in\mathcal{X}}\in\mathcal{M}_{1}^{+}([|1,r|])^{\mathcal{X}}$
and
$(k_{\nu}(\cdot,y))_{y\in\mathcal{X}}\in\mathcal{M}_{1}^{+}([|1,r|])^{\mathcal{X}}$
corresponding to the disintegration with respect to the projection of
respectively $\mu$ and $\nu$. Let us now consider $n$ independent samples
$(Z^{\mu}_{i})_{i=1}^{n}$ and $(Z^{\nu}_{i})_{i=1}^{n}$ such that for all
$i\in[|1,n|]$, $Z_{i}^{\mu}\sim k_{\mu}(\cdot,X_{i})$ and $Z_{i}^{\nu}\sim
k_{\nu}(\cdot,Y_{i})$ and let us define for all $k\in[|1,r|]$

$\displaystyle\tilde{\mu}_{k}\triangleq\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{Z^{\mu}_{i}=k}\delta_{X_{i}}~{}~{}\text{and}~{}~{}\tilde{\nu}_{k}\triangleq\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{Z^{\nu}_{i}=k}\delta_{Y_{i}}\;.$

Let us now define

$\displaystyle\tilde{\pi}$
$\displaystyle\triangleq\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\tilde{\mu}_{k}\otimes\tilde{\nu}_{k}$
$\displaystyle+\frac{1}{1-\sum_{k=1}^{r-1}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}\left[\hat{\mu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\tilde{\mu}_{k}\right]\otimes\left[\hat{\nu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\tilde{\nu}_{k}\right]$

with the convention that $\frac{0}{0}=0$. Now it is easy to check that
$\tilde{\pi}\in\Pi_{r}(\hat{\mu},\hat{\nu})$, indeed we have that

$\displaystyle\tilde{\pi}(A\times\mathcal{X})$
$\displaystyle=\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\tilde{\mu}_{k}(A)$
$\displaystyle+\frac{1}{1-\sum_{k=1}^{r-1}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}\left[\hat{\mu}(A)-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\tilde{\mu}_{k}(A)\right]\left[1-\sum_{k=1}^{r-1}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\right]$
$\displaystyle=\hat{\mu}(A)$

in addition by construction we have that

$\left|\hat{\mu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\tilde{\mu}_{k}\right|=\left|\hat{\nu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\tilde{\nu}_{k}\right|=1-\sum_{k=1}^{r-1}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)$

and both
$\hat{\mu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\tilde{\mu}_{k}$
and
$\hat{\nu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\tilde{\nu}_{k}$
are positive measures. Therefore we obtain that

$\displaystyle\text{LOT}_{r,c}(\hat{\mu},\hat{\nu})\leq\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\pi}(x,y)$

Now we aim at showing at
$\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\pi}(x,y)\rightarrow\text{LOT}_{r,c}(\mu,\nu)~{}~{}a.s.$.
Indeed first observe that from the law of large numbers we have that for all
$k\in[|1,r|]$, $|\tilde{\mu}_{k}|\rightarrow\lambda_{k}^{*}$ and similarly
$|\tilde{\nu}_{k}|\rightarrow\lambda_{k}^{*}$. In addition, for all $k,q$ we
have that almost surely, $\tilde{\mu}_{k}\otimes\tilde{\nu}_{q}$ converges
weakly towards $\lambda_{k}^{*}\lambda_{q}^{*}\mu_{k}\otimes\nu_{q}$. Indeed
one can consider the following algebra
$\mathcal{F}\triangleq\left\\{(x,y)\in\mathcal{X}^{2}\rightarrow
f(x)g(y)~{}~{}f,g\in\mathcal{C}(\mathcal{X})\right\\}$, and then by Stone-
Weierstrass, one obtains by density the desired result. Now remark that

$\displaystyle\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\pi}(x,y)$
$\displaystyle=\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\mu}_{k}\otimes\tilde{\nu}_{k}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\int_{\mathcal{Z}^{2}}c(x,y)d\tilde{\mu}_{r}\otimes\tilde{\nu}_{r}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\right)\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\mu}_{r}\otimes\tilde{\nu}_{k}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\mu}_{k}\otimes\tilde{\nu}_{r}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\sum_{k,q=1}^{r-1}\int_{\mathcal{X}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\left(1-\frac{\min(|\tilde{\mu}_{q}|,|\tilde{\nu}_{q}|)}{|\tilde{\nu}_{q}|}\right)c(x,y)d\tilde{\mu}_{k}(x)d\tilde{\nu}_{q}(y)$

from which follows directly that
$\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\pi}(x,y)\rightarrow\text{LOT}_{r,c}(\mu,\nu)~{}~{}a.s.$
Let us now denote for all $n\geq 1$, $\pi_{n}$ a solution of
$\text{LOT}_{r,c}(\hat{\mu},\hat{\nu})$. Let $\omega\in\Omega$ an element of
the probability space where live the random variables $(X_{i})_{i\geq 0}$ and
$(Y_{i})_{i\geq 0}$ such that
$\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\pi}^{(\omega)}(x,y)\rightarrow\text{LOT}_{r,c}(\mu,\nu)$.
As $\mathcal{X}$ is compact Thanks to Prokhorov’s Theorem, we can extract a
sequence such that $(\pi_{n}^{(\omega)})_{n\geq 0}$ converge weakly towards
$\pi^{(\omega)}\in\Pi_{r}(\mu,\nu)$. In addition we have that for all $n\geq
1$

$\displaystyle\int_{\mathcal{X}^{2}}c(x,y)d\pi_{n}^{(\omega)}(x,y)\leq\int_{\mathcal{X}^{2}}c(x,y)d\tilde{\pi}^{(\omega)}(x,y)$

And by considering the limit we obtain that

$\displaystyle\int c(x,y)d\pi^{(\omega)}(x,y)\leq\text{LOT}_{r,c}(\mu,\nu)$

However $\pi^{(\omega)}\in\Pi_{r}(\mu,\nu)$ and by optimality we obtain that

$\displaystyle\int c(x,y)d\pi^{(\omega)}(x,y)=\text{LOT}_{r,c}(\mu,\nu)$

This holds for an arbitrary subsequence of $(\pi_{n}^{(\omega)})_{n\geq 0}$,
from which follows that $\int
c(x,y)d\pi_{n}^{(\omega)}(x,y)\rightarrow\text{LOT}_{r,c}(\mu,\nu)$. Finally
this holds almost surely and the result follows. ∎

###  B.4 Proof of Proposition 4

######  Proposition.

Let $r\geq 1$ and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$. Then, there
exists a constant $K_{r}$ such that for any $\delta>0$ and $n\geq 1$, we have,
with a probability of at least $1-2\delta$, that

$\displaystyle\mathrm{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})-\mathrm{LOT}_{r,c}(\mu,\nu)\leq
11\|c\|_{\infty}\sqrt{\frac{r}{n}}+K_{r}\|c\|_{\infty}\left[\sqrt{\frac{\log(40/\delta)}{n}}+\frac{\sqrt{r}\log(40/\delta)}{n}\right]$

###### Proof.

We reintroduce the same notation as in the proof of Proposition 3. Let
$\pi^{*}$ solution of $\text{LOT}_{r,c}(\mu,\nu)$. Then there exists
$\lambda^{*}\in\Delta_{r}^{*}$,
$(\mu_{i}^{*})_{i=1}^{r},(\nu_{i}^{*})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{Z})^{r}$
such that

$\pi^{*}=\sum_{i=1}^{r}\lambda_{i}^{*}\mu_{i}^{*}\otimes\nu_{i}^{*}.$

As before let us also consider $\pi_{\mu}$ and $\pi_{\mu}$ defined as
$\pi_{\mu}(A\times\\{k\\})\triangleq\lambda_{k}\mu_{k}(A)\text{~{}~{}and~{}~{}}\pi_{\nu}(A\times\\{k\\})\triangleq\lambda_{k}\nu_{k}(A)~{}~{}\text{for
any measurable set $A$ and $k\in[|1,r|]$}$ and denote $\rho$ their common
right marginal. We also consider $n$ independent samples
$(Z^{\mu}_{i})_{i=1}^{n}$ and $(Z^{\nu}_{i})_{i=1}^{n}$ such that for all
$i\in[|1,n|]$, $Z_{i}^{\mu}\sim k_{\mu}(\cdot,X_{i})$ and $Z_{i}^{\nu}\sim
k_{\nu}(\cdot,Y_{i})$ and we denote for all $k\in[|1,r|]$

$\displaystyle\tilde{\mu}_{k}\triangleq\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{Z^{\mu}_{i}=k}\delta_{X_{i}}~{}~{}\text{and}~{}~{}\tilde{\nu}_{k}\triangleq\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{Z^{\nu}_{i}=k}\delta_{Y_{i}}$

Let us now define

$\hat{\pi}\triangleq\sum_{i=1}^{r}\frac{1}{\lambda_{k}^{*}}\tilde{\mu}_{k}\otimes\tilde{\nu}_{k}\;.$

Our goal is to control the following quantity:

$\displaystyle\left|\text{LOT}_{r,c}(\mu,\nu)-\int_{\mathcal{Z}^{2}}c(x,y)d\hat{\pi}(x,y)\right|,$

First observe that

$\displaystyle\mathbb{E}\left[\int_{\mathcal{Z}^{2}}c(x,y)d\hat{\pi}(x,y)\right]$
$\displaystyle=\sum_{i=1}^{r}\frac{1}{\lambda_{k}^{*}}\mathbb{E}\left[\int_{\mathcal{Z}^{2}}c(x,y)d\tilde{\mu}_{k}(x)d\tilde{\nu}_{k}(y)\right]$
$\displaystyle=\sum_{i=1}^{r}\frac{1}{\lambda_{k}^{*}n^{2}}\times\sum_{i,j}\mathbb{E}\left[c(X_{i},Y_{j})\mathbf{1}_{Z^{\mu}_{i}=k}\mathbf{1}_{Z^{\nu}_{j}=k}\right]$

Moreover, we have that

$\displaystyle\mathbb{E}\left[c(X_{i},Y_{j})\mathbf{1}_{Z^{\mu}_{i}=k}\mathbf{1}_{Z^{\nu}_{j}=k}\right]$
$\displaystyle=\int_{(\mathcal{Z}\times[|1,r|])^{2}}c(x,y)\mathbf{1}_{z=k}\mathbf{1}_{z^{\prime}=k}d\pi_{\mu}(x,z)d\pi_{\nu}(y,z^{\prime})$
$\displaystyle=\int_{(\mathcal{Z}\times[|1,r|])^{2}}c(x,y)\mathbf{1}_{z=k}\mathbf{1}_{z^{\prime}=k}d\mu_{z}(x)d\nu_{z^{\prime}}(y)d\rho(z)d\rho(z^{\prime})$
$\displaystyle=\lambda_{k}^{2}\int_{\mathcal{Z}^{2}}c(x,y)d\mu_{k}(x)d\nu_{k}(y)$

from which follows that

$\displaystyle\mathbb{E}\left[\int_{\mathcal{Z}^{2}}c(x,y)d\hat{\pi}(x,y)\right]$
$\displaystyle=\sum_{i=1}^{r}\lambda_{k}^{*}\int_{\mathcal{Z}^{2}}c(x,y)d\mu_{k}(x)d\nu_{k}(y)=\text{LOT}_{r,c}(\mu,\nu)$

Now let us define for all
$(x_{i},z_{i})_{i=1}^{n},(y_{i},z_{i}^{\prime})\in(\mathcal{Z}\times[|1,r|])^{n}$,

$\displaystyle
g((x_{1},z_{1}),\dots,(x_{n},z_{n}),(y_{1},z_{1}^{\prime}),\dots,(y_{n},z_{n}^{\prime}))\triangleq\sum_{q=1}^{r}\frac{1}{\lambda_{q}^{*}n^{2}}\sum_{i,j}c(x_{i},y_{j})\mathbf{1}_{z_{i}=q}\mathbf{1}_{z_{j}^{\prime}=q}\;,$

since $\mathcal{Z}$ is compact and $c$ is continuous, we have that

$\displaystyle|g(\dots,(x_{k},z_{k}),\dots)-g(\dots,(\tilde{x}_{k},\tilde{z}_{k}),\dots)|$
$\displaystyle=\left|\sum_{q=1}^{r}\frac{1}{\lambda_{q}^{*}n^{2}}\sum_{j}[c(x_{k},y_{j})\mathbf{1}_{z_{k}=q}-c(\tilde{x}_{k},y_{j})\mathbf{1}_{\tilde{z}_{k}=q}]\mathbf{1}_{z_{j}^{\prime}=q}\right|$
$\displaystyle=\left|\frac{1}{\lambda_{z_{k}}^{*}n^{2}}\sum_{j=1}^{n}c(x_{k},y_{j})\mathbf{1}_{z_{j}^{\prime}=z_{k}}-\frac{1}{\lambda_{\tilde{z}_{k}}^{*}n^{2}}\sum_{j=1}^{n}c(\tilde{x}_{k},y_{j})\mathbf{1}_{z_{j}^{\prime}=\tilde{z}_{k}}\right|$
$\displaystyle\leq\frac{\|c\|_{\infty}}{n^{2}}\left[\frac{\sum_{j=1}^{n}\mathbf{1}_{z_{j}^{\prime}=z_{k}}}{\lambda_{z_{k}}^{*}}+\frac{\sum_{j=1}^{n}\mathbf{1}_{z_{j}^{\prime}=\tilde{z}_{k}}}{\lambda_{\tilde{z}_{k}}^{*}}\right]$
$\displaystyle\leq\frac{2\|c\|_{\infty}}{\min\limits_{1\leq q\leq
r}\lambda_{q}^{*}}\frac{1}{n}$

Then by applying the McDiarmid’s inequality we obtain that for $\delta>0$,
with a probability at least of $1-\delta$, we have

$\displaystyle\left|\text{LOT}_{r,c}(\mu,\nu)-\int_{\mathcal{Z}^{2}}c(x,y)d\hat{\pi}(x,y)\right|\leq\frac{2\|c\|_{\infty}}{\min\limits_{1\leq
q\leq r}\lambda_{q}^{*}}\sqrt{\frac{\log(2/\delta)}{n}}$

Now we aim at building a coupling $\tilde{\pi}\in\Pi_{r}(\hat{\mu},\hat{\nu})$
from $\hat{\pi}$. Let us consider the same as the one introduce in the proof
of Proposition B.3, that is

$\displaystyle\tilde{\pi}$
$\displaystyle\triangleq\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\tilde{\mu}_{k}\otimes\tilde{\nu}_{k}$
$\displaystyle+\frac{1}{1-\sum_{k=1}^{r-1}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}\left[\hat{\mu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\tilde{\mu}_{k}\right]\otimes\left[\hat{\nu}-\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\tilde{\nu}_{k}\right]$

with the convention that $\frac{0}{0}=0$. Let us now expand the above
expression, and by denoting
$\tilde{\lambda}_{r}=1-\sum_{k=1}^{r-1}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)$
we obtain that

$\displaystyle\tilde{\pi}$
$\displaystyle=\sum_{k=1}^{r-1}\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\tilde{\mu}_{k}\otimes\tilde{\nu}_{k}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\tilde{\mu}_{r}\otimes\tilde{\nu}_{r}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\tilde{\mu}_{r}\otimes\left[\sum_{k=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\right)\tilde{\nu}_{k}\right]$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\left[\sum_{k=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\tilde{\mu}_{k}\right]\otimes\tilde{\nu}_{r}$
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\left[\sum_{k=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\tilde{\mu}_{k}\right]\otimes\left[\sum_{k=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\right)\tilde{\nu}_{k}\right]$

Now we aim at controlling the following quantity
$\left|\int_{\mathcal{Z}^{2}}c(x,y)d\hat{\pi}(x,y)-\int_{\mathcal{Z}^{2}}c(x,y)d\tilde{\pi}(x,y)\right|$
and we observe that

$\displaystyle\int_{\mathcal{Z}^{2}}c(x,y)d[\hat{\pi}(x,y)-\tilde{\pi}(x,y)]=\sum_{k=1}^{r-1}\int_{\mathcal{Z}^{2}}c(x,y)\left[\frac{1}{\lambda_{k}^{*}}-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\right]d\tilde{\mu}_{k}(x)\tilde{\nu}_{k}(y)$
(11)
$\displaystyle+\int_{\mathcal{Z}^{2}}c(x,y)\left[\frac{1}{\lambda_{r}^{*}}-\frac{1}{\tilde{\lambda}_{r}}\right]d\tilde{\mu}_{r}(x)\tilde{\nu}_{r}(y)$
(12)
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\right)c(x,y)d\tilde{\mu}_{r}(x)d\tilde{\nu}_{k}(y)$
(13)
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)c(x,y)d\tilde{\mu}_{k}(x)d\tilde{\nu}_{r}(y)$
(14)
$\displaystyle+\frac{1}{\tilde{\lambda}_{r}}\sum_{k,q=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\left(1-\frac{\min(|\tilde{\mu}_{q}|,|\tilde{\nu}_{q}|)}{|\tilde{\nu}_{q}|}\right)c(x,y)d\tilde{\mu}_{k}(x)d\tilde{\nu}_{q}(y)$
(15)

Let us now control each term of the RHS of the above equality. Let us first
consider the term in Eq. 11, remark that we have

$\displaystyle\left|\int_{\mathcal{Z}^{2}}c(x,y)\left[\frac{1}{\lambda_{k}^{*}}-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\right]d\tilde{\mu}_{k}(x)\tilde{\nu}_{k}(y)\right|$
$\displaystyle\leq\left|\left[\frac{1}{\lambda_{k}^{*}}-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\right]\right|\|c\|_{\infty}|\tilde{\mu}_{k}||\tilde{\nu}_{k}|$
$\displaystyle\leq\left|\left[\frac{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}{\lambda_{k}^{*}}-\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\right]\right|\|c\|_{\infty}$
$\displaystyle\leq\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\left|\frac{\max(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{\lambda_{k}^{*}}-1\right|\|c\|_{\infty}$
$\displaystyle\leq\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{\lambda_{k}^{*}}|\max(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)-\lambda_{k}^{*}|\|c\|_{\infty}$
$\displaystyle\leq\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{\lambda_{k}^{*}}\max(\|\tilde{\mathbf{\lambda}}_{\mu}-\mathbf{\lambda}^{*}\|_{\infty},\|\tilde{\mathbf{\lambda}}_{\nu}-\mathbf{\lambda}^{*}\|_{\infty})\|c\|_{\infty}$
$\displaystyle\leq\|c\|_{\infty}\max\left(\Big{\|}\frac{\tilde{\lambda}_{\mu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty},\Big{\|}\frac{\tilde{\lambda}_{\nu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty}\right)\max(\|\tilde{\mathbf{\lambda}}_{\mu}-\mathbf{\lambda}^{*}\|_{\infty},\|\tilde{\mathbf{\lambda}}_{\nu}-\mathbf{\lambda}^{*}\|_{\infty})$

where we have denoted
$\tilde{\mathbf{\lambda}}_{\mu}\triangleq(|\tilde{\mu}_{k}|)_{k=1}^{r}$ and
$\tilde{\mathbf{\lambda}}_{\nu}\triangleq(|\tilde{\nu}_{k}|)_{k=1}^{r}$. Now
observe that

$\displaystyle\mathbb{P}\left(\max(\|\tilde{\mathbf{\lambda}}_{\mu}-\mathbf{\lambda}^{*}\|_{\infty},\|\tilde{\mathbf{\lambda}}_{\nu}-\mathbf{\lambda}^{*}\|_{\infty})\geq
t\right)$ $\displaystyle\leq
2\mathbb{P}\left(\|\tilde{\mathbf{\lambda}}_{\mu}-\mathbf{\lambda}^{*}\|_{\infty}\geq
t\right)$
$\displaystyle\leq\mathbb{P}\left(d_{K}(\mathbf{\lambda}^{*},\tilde{\mathbf{\lambda}}_{\mu})\geq\frac{t}{2}\right)$
$\displaystyle\leq 4\exp(-nt^{2}/2)$

where $d_{K}$ is the Kolmogorov distance. In addition we have

$\displaystyle\max\left(\Big{\|}\frac{\tilde{\lambda}_{\mu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty},\Big{\|}\frac{\tilde{\lambda}_{\nu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty}\right)\leq
1+\frac{1}{\min\limits_{1\leq i\leq
r}\lambda_{i}^{*}}\max\left(\|\tilde{\lambda}_{\mu}-\mathbf{\lambda^{*}}\|_{\infty},\|\tilde{\lambda}_{\nu}-\mathbf{\lambda^{*}}\|_{\infty}\right)$

Combining the two above controls, we obtain that for all $\delta>0$, with a
probability of at least $1-\delta$,

$\displaystyle\left|\int_{\mathcal{Z}^{2}}c(x,y)\left[\frac{1}{\lambda_{k}^{*}}-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}||\tilde{\nu}_{k}|}\right]d\tilde{\mu}_{k}(x)\tilde{\nu}_{k}(y)\right|\leq\|c\|_{\infty}\sqrt{{\frac{2\ln{8/\delta}}{n}}}+\frac{\|c\|_{\infty}}{n}\frac{2\ln{8/\delta}}{\min\limits_{1\leq
i\leq r}\lambda_{i}^{*}}$

Let us now consider the term in Eq. 12, we have that

$\displaystyle\left|\int_{\mathcal{Z}^{2}}c(x,y)\left[\frac{1}{\lambda_{r}^{*}}-\frac{1}{\tilde{\lambda}_{r}}\right]d\tilde{\mu}_{r}(x)\tilde{\nu}_{r}(y)\right|$
$\displaystyle\leq\frac{|\tilde{\mu}_{r}||\tilde{\nu}_{r}|}{\lambda_{r}^{*}\tilde{\lambda}_{r}}\left|1-\sum_{i=1}^{r}\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)-\lambda_{r}\right|\|c\|_{\infty}$
$\displaystyle\leq\max\left(\Big{\|}\frac{\tilde{\lambda}_{\mu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty},\Big{\|}\frac{\tilde{\lambda}_{\nu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty}\right)\sum_{k=1}^{r-1}\left|\lambda_{k}^{*}-\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\right|\|c\|_{\infty}$
$\displaystyle\leq\max\left(\Big{\|}\frac{\tilde{\lambda}_{\mu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty},\Big{\|}\frac{\tilde{\lambda}_{\nu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty}\right)\|c\|_{\infty}(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1}+\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\nu}\|_{1})$
$\displaystyle\leq
2\|c\|_{\infty}\max\left(\Big{\|}\frac{\tilde{\lambda}_{\mu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty},\Big{\|}\frac{\tilde{\lambda}_{\nu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty}\right)\max(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1},\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\nu}\|_{1})$

However we have that

$\displaystyle\mathbb{P}\left(\max(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1},\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\nu}\|_{1})\geq
t\right)\leq
2\mathbb{P}\left(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1}\geq
t\right)$

In addition we have that
$\mathbb{E}(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1})\leq\sqrt{\frac{r}{n}}$
and by applying the McDiarmid’s Inequality, we obtain that for all $\delta>0$,
with a probability of $1-\delta$

$\displaystyle\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1}\leq\sqrt{\frac{r}{n}}+\sqrt{\frac{2\ln(2/\delta)}{n}}$

Therefore we obtain that with a probability of at least $1-\delta$,

$\displaystyle\left|\int_{\mathcal{Z}^{2}}c(x,y)\left[\frac{1}{\lambda_{r}^{*}}-\frac{1}{\tilde{\lambda}_{r}}\right]d\tilde{\mu}_{r}(x)\tilde{\nu}_{r}(y)\right|\leq
2\|c\|_{\infty}\left[\sqrt{\frac{r}{n}}+\sqrt{\frac{2\ln(8/\delta)}{n}}+\frac{2\ln(8/\delta)+\sqrt{2r\ln(8/\delta)}}{n\times\min\limits_{1\leq
i\leq r}\lambda_{i}^{*}}\right]$

For the term in Eq. 13 and 14, we obtain that

$\displaystyle\left|\frac{1}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\right)c(x,y)d\tilde{\mu}_{r}(x)d\tilde{\nu}_{k}(y)\right|$
$\displaystyle\leq\frac{|\tilde{\mu}_{r}|}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\left(|\tilde{\nu}_{k}|-\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\right)\|c\|_{\infty}$
$\displaystyle\leq\frac{|\tilde{\mu}_{r}|}{\tilde{\lambda}_{r}}[\tilde{\lambda}_{r}-|\tilde{\nu}_{r}|]\|c\|_{\infty}$
$\displaystyle\leq[|\tilde{\lambda}_{r}-\lambda_{r}^{*}|+|\lambda_{r}^{*}-\tilde{\nu}_{r}|]\|c\|_{\infty}$
$\displaystyle\leq
3\|c\|_{\infty}\max(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1},\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\nu}\|_{1})$

Therefore we obtain that with a probability of at least $1-\delta$,

$\displaystyle\left|\frac{1}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\nu}_{k}|}\right)c(x,y)d\tilde{\mu}_{r}(x)d\tilde{\nu}_{k}(y)\right|\leq
3\|c\|_{\infty}\left[\sqrt{\frac{r}{n}}+\sqrt{\frac{2\ln(2/\delta)}{n}}\right]$

Finally the last term in Eq. 15 can be controlled as the following:

$\displaystyle\left|\frac{1}{\tilde{\lambda}_{r}}\sum_{k,q=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\left(1-\frac{\min(|\tilde{\mu}_{q}|,|\tilde{\nu}_{q}|)}{|\tilde{\nu}_{q}|}\right)c(x,y)d\tilde{\mu}_{k}(x)d\tilde{\nu}_{q}(y)\right|$
$\displaystyle\leq\frac{\|c\|_{\infty}}{\tilde{\lambda}_{r}}\sum_{k,q=1}^{r-1}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\left(1-\frac{\min(|\tilde{\mu}_{q}|,|\tilde{\nu}_{q}|)}{|\tilde{\nu}_{q}|}\right)|\tilde{\mu}_{k}||\tilde{\nu}_{q}|$
$\displaystyle\leq\frac{\|c\|_{\infty}}{\tilde{\lambda}_{r}}\sum_{k=1}^{r-1}\left(|\tilde{\mu}_{k}|-\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\right)\sum_{k=1}^{r-1}\left(|\tilde{\nu}_{k}|-\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)\right)$
$\displaystyle\leq
3\|c\|_{\infty}\max(\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\mu}\|_{1},\|\mathbf{\lambda}^{*}-\tilde{\mathbf{\lambda}}_{\nu}\|_{1})$

and we obtain that with a probability of at least $1-\delta$,

$\displaystyle\left|\frac{1}{\tilde{\lambda}_{r}}\sum_{k,q=1}^{r-1}\int_{\mathcal{Z}^{2}}\left(1-\frac{\min(|\tilde{\mu}_{k}|,|\tilde{\nu}_{k}|)}{|\tilde{\mu}_{k}|}\right)\left(1-\frac{\min(|\tilde{\mu}_{q}|,|\tilde{\nu}_{q}|)}{|\tilde{\nu}_{q}|}\right)c(x,y)d\tilde{\mu}_{k}(x)d\tilde{\nu}_{q}(y)\right|$
$\displaystyle\leq
3\|c\|_{\infty}\left[\sqrt{\frac{r}{n}}+\sqrt{\frac{2\ln(2/\delta)}{n}}\right]$

Then by applying a union bound we obtain that with a probability of at least
$1-\delta$

$\displaystyle\left|\int_{\mathcal{Z}^{2}}c(x,y)d[\hat{\pi}(x,y)-\tilde{\pi}(x,y)]\right|\leq\|c\|_{\infty}\left[11\sqrt{\frac{r}{n}}+12\sqrt{{\frac{2\ln{40/\delta}}{n}}}+\frac{6\ln(40/\delta)+2\sqrt{2r\ln(40/\delta)}}{n\times\min\limits_{1\leq
i\leq r}\lambda_{i}^{*}}\right]$

Now observe that

$\displaystyle\text{LOT}_{r,c}(\hat{\mu},\hat{\nu})-\text{LOT}_{r,c}(\mu,\nu)$
$\displaystyle\leq\int_{\mathcal{Z}^{2}}c(x,y)d\tilde{\pi}(x,y)-\int_{\mathcal{Z}^{2}}c(x,y)d\pi^{*}(x,y)$
$\displaystyle\leq\int_{\mathcal{Z}^{2}}c(x,y)d[\tilde{\pi}-\hat{\pi}](x,y)+\int_{\mathcal{Z}^{2}}c(x,y)d[\hat{\pi}-\pi^{*}](x,y)$

and by combining the two control we obtain that with a probability of at least
$1-2\delta$,

$\displaystyle\text{LOT}_{r,c}(\hat{\mu},\hat{\nu})-\text{LOT}_{r,c}(\mu,\nu)$
$\displaystyle\leq\|c\|_{\infty}\left[11\sqrt{\frac{r}{n}}+12\sqrt{{\frac{2\ln{40/\delta}}{n}}}+\frac{1}{\alpha}\left(2\sqrt{\frac{\log(2/\delta)}{n}}+\frac{6\ln(40/\delta)+2\sqrt{2r\ln(40/\delta)}}{n}\right)\right]$
$\displaystyle\leq
11\|c\|_{\infty}\sqrt{\frac{r}{n}}+\frac{14\|c\|_{\infty}}{\alpha}\sqrt{\frac{\log(40/\delta)}{n}}+\frac{2\|c\|_{\infty}\max(6,\sqrt{2r})\log(40/\delta)}{n\alpha}$

where $\alpha\triangleq\min\limits_{1\leq i\leq r}\lambda_{i}^{*}$ and the
result follows. ∎

###  B.5 Proof Proposition 5

######  Proposition.

Let $r\geq 1$, $\delta>0$ and $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$.
Then there exists a constant $N_{r,\delta}$ such that if $n\geq N_{r,\delta}$
then with a probability of at least $1-2\delta$, we have

$\displaystyle\mathrm{LOT}_{r,c}(\hat{\mu}_{n},\hat{\nu}_{n})-\mathrm{LOT}_{r,c}(\mu,\nu)\leq
11\|c\|_{\infty}\sqrt{\frac{r}{n}}+77\|c\|_{\infty}\sqrt{\frac{\log(40/\delta)}{n}}\;.$

###### Proof.

We consider the same notations as in the proof of Proposition 4. In particular
let us define for all
$(x_{i},z_{i})_{i=1}^{n},(y_{i},z_{i}^{\prime})\in(\mathcal{Z}\times[|1,r|])^{n}$,

$\displaystyle
g((x_{1},z_{1}),\dots,(x_{n},z_{n}),(y_{1},z_{1}^{\prime}),\dots,(y_{n},z_{n}^{\prime}))\triangleq\sum_{q=1}^{r}\frac{1}{\lambda_{q}^{*}n^{2}}\sum_{i,j}c(x_{i},y_{j})\mathbf{1}_{z_{i}=q}\mathbf{1}_{z_{j}^{\prime}=q}\;,$

Recall that we have

$\displaystyle|g(\dots,(x_{k},z_{k}),\dots)-g(\dots,(\tilde{x}_{k},\tilde{z}_{k}),\dots)|$
$\displaystyle\leq\frac{\|c\|_{\infty}}{n^{2}}\left[\frac{\sum_{j=1}^{n}\mathbf{1}_{z_{j}^{\prime}=z_{k}}}{\lambda_{z_{k}}^{*}}+\frac{\sum_{j=1}^{n}\mathbf{1}_{z_{j}^{\prime}=\tilde{z}_{k}}}{\lambda_{\tilde{z}_{k}}^{*}}\right]$
$\displaystyle\leq\frac{2\|c\|_{\infty}}{n}\max\left(\Big{\|}\frac{\tilde{\lambda}_{\mu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty},\Big{\|}\frac{\tilde{\lambda}_{\nu}}{\mathbf{\lambda^{*}}}\Big{\|}_{\infty}\right)$
$\displaystyle\leq\frac{2\|c\|_{\infty}}{n}+\frac{2\|c\|_{\infty}}{n\times\min\limits_{1\leq
i\leq
r}\lambda_{i}^{*}}\max\left(\|\tilde{\lambda}_{\mu}-\mathbf{\lambda^{*}}\|_{\infty},\|\tilde{\lambda}_{\nu}-\mathbf{\lambda^{*}}\|_{\infty}\right)$

In fact if we have a control in probability of the bounded difference we can
use an extension of the McDiarmid’s Inequality. For that purpose let us first
introduce the following definition.

######  Definition 4.

Let $(X_{i})_{i=1}^{m}$, $m$ independent random variables and $g$ a measurable
function. We say that g is weakly difference-bounded with respect to
$(X_{i})_{i=1}^{m}$ by $(b,\beta,\delta)$ if

$\displaystyle\mathbb{P}\left(|g(X_{1},\dots,X_{m})-g(X_{1}^{\prime},\dots,X_{m}^{\prime})|\leq\beta\right)\geq
1-\delta$

with $X_{i}^{\prime}=X_{i}$ except for one coordinate $k$ where
$X_{k}^{\prime}$ is an independent copy of $X_{k}$. Furthermore for any
$(x_{i})_{i=1}^{m}$ and $(x_{i}^{\prime})_{i=1}^{m}$ where for all coordinate
except on $x_{j}=x_{j}^{\prime}$

$\displaystyle|g(x_{1},\dots,x_{m})-g(x_{1}^{\prime},\dots,x_{m}^{\prime})|\leq
b\;.$

Let us now introduce an extension of McDiarmid’s Inequality [Kutin, 2002].

######  Theorem 1.

Let $(X_{i})_{i=1}^{m}$, $m$ independent random variables and $g$ a measurable
function which is weakly difference-bounded with respect to
$(X_{i})_{i=1}^{m}$ by $(b,\beta/m,\exp(-Km))$, then if $0<\tau\leq
T(b,\beta,K)$ and $m\geq M(b,\beta,K,\tau)$, then

$\displaystyle\mathbb{P}(|g(X_{1},\dots,X_{m})-\mathbb{E}(g(X_{1},\dots,X_{m}))|\geq\tau)\leq
4\exp\left(\frac{-\tau^{2}m}{8\beta^{2}}\right)$

where

$\displaystyle T(b,\beta,K)$
$\displaystyle\triangleq\min\left(\frac{14c}{2},4\beta\sqrt{K},\frac{\beta^{2}K}{b}\right)$
$\displaystyle M(b,\beta,K,\tau)$
$\displaystyle\triangleq\max\left(\frac{b}{\beta},\beta\sqrt{40},3\left(\frac{24}{K}+3\right)\log\left(\frac{24}{K}+3\right),\frac{1}{\tau}\right)$

Given the above Theroem we can obtain an asymptotic control of the deviation
of $g$ from its mean. Let $\delta^{\prime}>0$ and let us denote

$\displaystyle m$ $\displaystyle\triangleq 2n$ $\displaystyle b$
$\displaystyle\triangleq\frac{2\|c\|_{\infty}}{n\times\min\limits_{1\leq i\leq
r}\lambda_{i}^{*}}$ $\displaystyle K$
$\displaystyle\triangleq\frac{\log(1/\delta^{\prime})}{2n}$
$\displaystyle\beta$ $\displaystyle\triangleq
4\|c\|_{\infty}\left[1+\frac{1}{\min\limits_{1\leq i\leq
r}\lambda_{i}^{*}}\sqrt{\frac{2\log(4/\delta^{\prime})}{n}}\right]$

Observe now that with a probability of at least $1-\exp(-Km)$

$\displaystyle|g(\dots,(x_{k},z_{k}),\dots)-g(\dots,(\tilde{x}_{k},\tilde{z}_{k}),\dots)|\leq\frac{2\|c\|_{\infty}}{n}\left[1+\frac{1}{\min\limits_{1\leq
i\leq r}\lambda_{i}^{*}}\sqrt{\frac{2\log(4/\delta^{\prime})}{n}}\right]$

Let us now fix $\delta>0$ and let us choose $\delta^{\prime}$ such that
$\delta^{\prime}\triangleq 4/n$ and
$\tau\triangleq\beta\sqrt{\frac{4\log(4/\delta)}{n}}$, then we obtain that for
$n$ sufficiently large (such that $n\geq M(b,\beta,K,\tau)/2$ and $\tau\leq
T(b,\beta,K)$), we have that with a probability of at least $1-\delta$

$\displaystyle\left|\text{LOT}_{r,c}(\mu,\nu)-\int_{\mathcal{Z}^{2}}c(x,y)d\hat{\pi}(x,y)\right|$
$\displaystyle\leq 4\|c\|_{\infty}\left[1+\frac{1}{\min\limits_{1\leq i\leq
r}\lambda_{i}^{*}}\sqrt{\frac{2\log(n)}{n}}\right]\sqrt{\frac{4\log(4/\delta)}{n}}$
$\displaystyle\leq
4\|c\|_{\infty}\sqrt{\frac{4\log(4/\delta)}{n}}+\frac{16\sqrt{5}\|c\|_{\infty}\sqrt{\log(n)\log(4/\delta)}}{n\times\min\limits_{1\leq
i\leq r}\lambda_{i}^{*}}$

Recall also from the proof of Proposition 4, that we have with a probability
of at least $1-\delta$

$\displaystyle\left|\int_{\mathcal{Z}^{2}}c(x,y)d[\hat{\pi}(x,y)-\tilde{\pi}(x,y)]\right|\leq\|c\|_{\infty}\left[11\sqrt{\frac{r}{n}}+12\sqrt{{\frac{2\ln{40/\delta}}{n}}}+\frac{6\ln(40/\delta)+2\sqrt{2r\ln(40/\delta)}}{n\times\min\limits_{1\leq
i\leq r}\lambda_{i}^{*}}\right]$

Finally by imposing in addition that

$\displaystyle\sqrt{\frac{n}{\log(n)}}\geq\frac{1}{\min\limits_{1\leq i\leq
r}\lambda_{i}^{*}}\text{~{}~{},~{}~{}}\sqrt{n}\geq\frac{\sqrt{\log(40/\delta)}}{\min\limits_{1\leq
i\leq
r}\lambda_{i}^{*}}\text{~{}~{}and~{}~{}}\sqrt{n}\geq\frac{\sqrt{r}}{\min\limits_{1\leq
i\leq r}}$

we obtain that for $n$ is large enough (such that (such that $n\geq
M(b,\beta,K,\tau)/2$ and $\tau\leq T(b,\beta,K)$) and satysfing the above
inequalities, we have with a probability of at least $1-2\delta$ that

$\displaystyle\text{LOT}_{r,c}(\hat{\mu},\hat{\nu})-\text{LOT}_{r,c}(\mu,\nu)$
$\displaystyle\leq
11\|c\|_{\infty}\sqrt{\frac{r}{n}}+77\|c\|_{\infty}\sqrt{\frac{\log(40/\delta)}{n}}$

∎

###  B.6 Proof Proposition 6

######  Proposition.

Let $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$. Let us assume that $c$ is
symmetric, then we have

$\mathrm{DLOT}_{1,c}(\mu,\nu)=\frac{1}{2}\int_{\mathcal{X}^{2}}-c(x,y)d[\mu-\nu]\otimes
d[\mu-\nu](x,y)\;.$

If in addition we assume the $c$ is Lipschitz w.r.t to $x$ and $y$, then we
have

$\displaystyle\mathrm{DLOT}_{r,c}(\mu,\nu)\xrightarrow[r\rightarrow+\infty]{}\mathrm{OT}_{c}(\mu,\nu)\;.$

###### Proof.

When $r=1$, it is clear that for any
$\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$,
$\Pi_{r}(\mu,\nu)=\\{\mu\otimes\nu\\}$ and thanks to the symmetry of $c$, we
have directly that

$\displaystyle\text{DLOT}_{1,c}(\mu,\nu)=\frac{1}{2}\int_{\mathcal{X}^{2}}-c(x,y)d[\mu-\nu]\otimes
d[\mu-\nu](x,y)=\frac{1}{2}\text{MMD}_{-c}(\mu,\nu)\;.$

The limit is a direct consequence of Proposition 2. ∎

###  B.7 Proof of Proposition 8

######  Proposition.

Let $r\geq 1$ and $(\mu_{n})_{n\geq 0}$ and $(\nu_{n})_{n\geq 0}$ two
sequences of probability measures such that $\mu_{n}\rightarrow\mu$ and
$\nu_{n}\rightarrow\nu$ with respect to the convergence in law. Then we have
that

$\text{LOT}_{r,c}(\mu_{n},\nu_{n})\rightarrow\text{LOT}_{r,c}(\mu,\nu)\;.$

###### Proof.

Let us denote $\pi$ an optimal solution of $\text{LOT}_{r,c}(\mu,\nu)$ and let
us denote $(\mu^{(i)})_{i=1}^{r}$, $(\nu^{(i)})_{i=1}^{r}$ and
$(\lambda^{(i)})_{i=1}^{r}$ the decomposition associated. In the following
Lemma, we aim at building specific decompositions of the sequences
$(\mu_{n})_{n\geq 0}$ and $(\nu_{n})_{n\geq 0}$.

######  Lemma 1.

Let $r\geq 1$, $\mu\in\mathcal{M}_{1}^{+}(\mathcal{X})$ and
$(\mu^{(i)})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{X})$ and
$(\lambda^{(i)})_{i=1}^{r}\in\Delta_{r}^{*}$ such that
$\mu=\sum_{i=1}^{r}\lambda_{i}\mu^{(i)}$. Then for any sequence of probability
measures $(\mu_{n})_{\geq 0}$ such that $\mu_{n}\rightarrow\mu$, there exist
for all $i\in[|1,r|]$ a sequence of nonnegative measures
$(\mu_{n}^{(i)})_{n\geq 0}$ such that

$\displaystyle\mu_{n}^{(i)}\rightarrow\lambda_{i}\mu^{(i)}~{}~{}\text{for all
$i\in[|1,r|]$ and}$
$\displaystyle\sum_{i=1}^{r}\mu_{n}^{(i)}=\mu_{n}~{}~{}\text{for all $n\geq
0$}$

###### Proof.

For $r=1$ the result is clear. Let us now show the result for $r=2$. Let us
denote $(\tilde{\mu}^{(1)}_{n})$ a sequence converging weakly towards
$\lambda_{1}\mu^{(1)}$. Then by denoting
$\mu^{(1)}_{n}\triangleq\mu_{n}-(\mu_{n}-\tilde{\mu}^{(1)}_{n})_{+}$ where
$(\cdot)_{+}$ correspond to the non-negative part of the measure, we have that

$\displaystyle\mu^{(1)}_{n}$ $\displaystyle\geq
0,~{}~{}\mu^{(1)}_{n}\rightarrow\lambda_{1}\mu^{(1)},$
$\displaystyle\mu^{(2)}_{n}$ $\displaystyle\triangleq\mu_{n}-\mu^{(1)}_{n}\geq
0,~{}~{}\mu^{(2)}_{n}\rightarrow\lambda_{2}\mu^{(2)}~{}~{}\text{and}$
$\displaystyle\mu_{n}$
$\displaystyle=\mu^{(1)}_{n}+\mu^{(2)}_{n}~{}~{}\text{for all}~{}~{}n\geq 0$

which is the result. Let $r\geq 2$ and let us assume that the result holds for
all $1\leq k\leq r$. Let us now consider a decomposition of $\mu$ such that
$\mu=\sum_{i=1}^{r+1}\lambda_{i}\mu^{(i)}$. By denoting
$\tilde{\mu}^{(1)}\triangleq\frac{\sum_{i=1}^{r}\lambda_{i}\mu^{(i)}}{\sum_{i=1}^{r}\lambda_{i}}$,
we obtain that

$\displaystyle\mu=\left(\sum_{i=1}^{r}\lambda_{i}\right)\tilde{\mu}^{(1)}+\lambda_{r+1}\mu^{(r+1)}\;.$

Then by recursion we have that there exists sequences of nonnegative measures
$(\tilde{\mu}_{n}^{(1)})$ and $(\mu_{n}^{(r+1)})$ such that

$\displaystyle\tilde{\mu}_{n}^{(1)}\rightarrow\left(\sum_{i=1}^{r}\lambda_{i}\right)\tilde{\mu}^{(1)},~{}~{}\mu_{n}^{(r+1)}\rightarrow\lambda_{r+1}\mu^{(r+1)}~{}~{}\text{and}~{}~{}\mu_{n}=\tilde{\mu}_{n}^{(1)}+\mu_{n}^{(r+1)}~{}~{}\text{for
all}~{}~{}n\geq 0$

Now observe that
$\frac{\tilde{\mu}_{n}^{(1)}}{|\tilde{\mu}_{n}^{(1)}|}\rightarrow\tilde{\mu}^{(1)}=\sum_{i=1}^{r}\frac{\lambda_{i}}{\sum_{i=1}^{r}\lambda_{i}}\mu^{(i)}$.
Therefore applying the recursion on this problem allows us to obtain a
decomposition of $\tilde{\mu}_{n}^{(1)}$ of the form

$\displaystyle\frac{\tilde{\mu}_{n}^{(1)}}{|\tilde{\mu}_{n}^{(1)}|}$
$\displaystyle=\sum_{i=1}^{r}\mu_{n}^{(i)}~{}~{}\text{where}$
$\displaystyle\mu_{n}^{(i)}$ $\displaystyle\geq
0~{}~{}\text{and}~{}~{}\mu_{n}^{(i)}\rightarrow\frac{\lambda_{i}}{\sum_{i=1}^{r}\lambda_{i}}\mu^{(i)}\;.$

Therefore we obtain that

$\displaystyle\mu_{n}$
$\displaystyle=\sum_{i=1}^{r}|\tilde{\mu}_{n}^{(1)}|\mu_{n}^{(i)}+\mu_{n}^{(r+1)}~{}~{}\text{where}~{}~{}$
$\displaystyle\mu_{n}^{(i)}$ $\displaystyle\geq
0,~{}~{}|\tilde{\mu}_{n}^{(1)}|\mu_{n}^{(i)}\rightarrow\lambda_{i}\mu^{(i)}~{}~{}\text{for
all}~{}~{}i\in[|1,r|]~{}~{}\text{and}$ $\displaystyle\mu_{n}^{(r+1)}$
$\displaystyle\geq 0,~{}~{}\mu_{n}^{(r+1)}\rightarrow\lambda_{r+1}\mu^{(r+1)}$

from which follows the result. ∎

Let us now consider such decompositions of $(\mu_{n})_{n\geq 0}$ and
$(\nu_{n})_{n\geq 0}$ such that each factor converges toward the target
decomposition of $\mu$. Now let us build the following coupling:

$\displaystyle\tilde{\pi}_{n}$
$\displaystyle\triangleq\sum_{k=1}^{r-1}\frac{\min(|\mu_{n}^{(k)}|,|\nu_{n}^{(k)}|)}{|\mu_{n}^{(k)}||\nu_{n}^{(k)}|}\mu_{n}^{(k)}\otimes\mu_{n}^{(k)}$
$\displaystyle+\frac{1}{1-\sum_{k=1}^{r-1}\min(|\mu_{n}^{(k)}|,|\nu_{n}^{(k}|)}\left[|\mu_{n}|-\sum_{k=1}^{r-1}\frac{\min(|\mu_{n}^{(k)}|,|\nu_{n}^{(k)}|)}{|\mu_{n}^{(k)}|}\mu_{n}^{(k)}\right]\otimes\left[\nu_{n}-\sum_{k=1}^{r-1}\frac{\min(|\mu_{n}^{(k)}|,|\nu_{n}^{(k)}|)}{|\nu_{n}^{(k)}|}\nu_{n}^{(k)}\right]$

with the convention that $\frac{0}{0}=0$. Now it is easy to check that
$\tilde{\pi}_{n}\in\Pi_{r}(\mu_{n},\nu_{n})$, and we have that

$\displaystyle\text{LOT}_{r,c}(\mu_{n},\nu_{n})\leq\int_{\mathcal{X}^{2}}d(x,y)d\tilde{\pi}_{n}(x,y)\rightarrow\text{LOT}_{r,c}(\mu,\nu)$

and by Prokhorov’s theorem and the optimality of the limit of
$(\tilde{\pi}_{n})_{n\geq 0}$ (up to an extraction) we obtain that
$\text{LOT}_{r,c}(\mu_{n},\nu_{n})\rightarrow\text{LOT}_{r,c}(\mu,\nu)$. ∎

###  B.8 Proof Proposition 7

######  Proposition.

Let $r\geq 1$, and let us assume that $c$ is a semimetric of negative type.
Then for all $\mu,\nu\in\mathcal{M}_{1}^{+}(\mathcal{X})$, we have that

$\text{DLOT}_{r}(\mu,\nu)\geq 0\;.$

In addition, if $c$ has strong negative type then we have also that

$\displaystyle\mathrm{DLOT}_{r,c}(\mu,\nu)=0$
$\displaystyle\iff\mu=\nu~{}~{}\text{and}$
$\displaystyle\mu_{n}\rightarrow\mu$
$\displaystyle\iff\mathrm{DLOT}_{r,c}(\mu_{n},\mu)\rightarrow 0\;.$

where the convergence of the sequence of probability measures considered is
the convergence in law.

###### Proof.

Let $\pi^{*}$ solution of $\text{LOT}_{r,c}(\mu,\nu)$. Then there exists
$\lambda^{*}\in\Delta_{r}^{*}$,
$(\mu_{i}^{*})_{i=1}^{r},(\nu_{i}^{*})_{i=1}^{r}\in\mathcal{M}_{1}^{+}(\mathcal{X})^{r}$
such that

$\pi^{*}=\sum_{i=1}^{r}\lambda_{i}^{*}\mu_{i}^{*}\otimes\nu_{i}^{*}.$

Note that by definition, we have that

$\displaystyle\mu=\sum_{i=1}^{r}\lambda_{i}^{*}\mu_{i}^{*}\text{~{}~{}and~{}~{}}\nu=\sum_{i=1}^{r}\lambda_{i}^{*}\nu_{i}^{*},$

By definition we have also that

$\displaystyle\text{LOT}_{r,c}(\mu,\mu)\leq\sum_{k=1}^{r}\lambda_{k}^{*}\int_{\mathcal{X}^{2}}c(x,y)d\mu_{k}^{*}\otimes\mu_{k}^{*}$

similarly for $\text{LOT}_{r,c}(\nu,\nu)$ we have

$\displaystyle\text{LOT}_{r,c}(\nu,\nu)\leq\sum_{k=1}^{r}\lambda_{k}^{*}\int_{\mathcal{X}^{2}}c(x,y)d\nu_{k}^{*}\otimes\nu_{k}^{*}$

Therefore we have

$\displaystyle\text{DLOT}_{r,c}(\mu,\nu)$
$\displaystyle\geq\sum_{k=1}^{r}\lambda_{k}^{*}\left(\int_{\mathcal{X}^{2}}c(x,y)d\mu_{k}^{*}\otimes\nu_{k}^{*}-\frac{1}{2}\left[\int_{\mathcal{X}^{2}}c(x,y)d\mu_{k}^{*}\otimes\mu_{k}^{*}+\int_{\mathcal{X}^{2}}c(x,y)d\nu_{k}^{*}\otimes\nu_{k}^{*}\right]\right)$
$\displaystyle\geq\sum_{k=1}^{r}\lambda_{k}^{*}\int_{\mathcal{X}^{2}}-c(x,y)d[\mu_{k}^{*}-\nu_{k}^{*}]\otimes[\mu_{k}^{*}-\nu_{k}^{*}]$
$\displaystyle\geq\sum_{k=1}^{r}\frac{\lambda_{k}^{*}}{2}D_{c}(\mu_{k}^{*},\nu_{k}^{*})$

where for any any probability measures $\alpha,\beta$ on $\mathcal{X}$ we
define

$\displaystyle D_{c}(\alpha,\beta)\triangleq
2\int_{\mathcal{X}^{2}}c(x,y)d\alpha\otimes\beta-\int_{\mathcal{X}^{2}}c(x,y)d\alpha\otimes\alpha-\int_{\mathcal{X}^{2}}c(x,y)d\beta\otimes\beta$

However, as $c$ is assumed to have a negative type, we have that

$\displaystyle D_{c}(\mu_{k}^{*},\nu_{k}^{*})\geq 0~{}~{}\forall k$

In addition if we assume that $c$ has a strong negative type, then we obtain
directly that

$\displaystyle\text{DLOT}_{r,c}(\mu,\nu)=0\implies\mu_{k}^{*}=\nu_{k}^{*}~{}~{}\forall
k\;.$

Let us now show that $\text{DLOT}_{r,c}$ metrize the convergence in law. The
direct implication is a direct consequence of the Proposition 8. Conversely,
if $\text{DLOT}_{r,c}(\mu_{n},\mu)\rightarrow 0$, then by compacity of
$\mathcal{X}$ and thanks to the Prokhorov’s theorem we can extract a
subsequence of $\mu_{n}\rightarrow\mu^{*}$, and thanks to Proposition 8, we
also obtain that
$\text{DLOT}_{r,c}(\mu_{n},\mu)\rightarrow\text{DLOT}_{r,c}(\mu^{*},\mu)$.
Finally we deduce that $\text{DLOT}_{r,c}(\mu^{*},\mu)=0$ and $\mu^{*}=\mu$.

∎

###  B.9 Proof Proposition 9

######  Proposition.

Let $n\geq k\geq 1$,
$\mathbf{X}\triangleq\\{x_{1},\dots,x_{n}\\}\subset\mathcal{X}$ and
$a\in\Delta_{n}^{*}$. If $c$ is a semimetric of negative type, then by
denoting $C=(c(x_{i},x_{j}))_{i,j}$, we have that

$\displaystyle\mathrm{LOT}_{k,c}(\mu_{a,\mathbf{X}},\mu_{a,\mathbf{X}})=\min_{Q}\langle
C,Q\text{diag}(1/Q^{T}\mathbf{1}_{n})Q^{T}\rangle~{}~{}\text{s.t.}~{}~{}Q\in\mathbb{R}^{n\times
k}_{+}~{}~{}\text{,}~{}~{}Q\mathbf{1}_{k}=a\;.$ (16)

###### Proof.

First remarks that one can reformulate the $\text{LOT}_{k,c}$ problem as

$\displaystyle\text{LOT}_{k,c}(\mu,\mu)\triangleq\min_{g\in\Delta_{k}^{*}}\min_{(\mathbf{x},\mathbf{y})\in
K_{a,g}^{2}}\sum_{i=1}^{k}\frac{\mathbf{x}_{i}^{T}C\mathbf{y}_{i}}{g_{i}}$

where

$\displaystyle K_{a,g}$
$\displaystyle\triangleq\\{\mathbf{x}\in\mathbb{R}^{nk}~{}\text{s.t.}~{}A\mathbf{x}=[a,g]^{T},~{}\mathbf{x}\geq
0\\}$ $\displaystyle A$
$\displaystyle\triangleq\begin{pmatrix}\begin{matrix}\mathbf{1}_{n}^{T}\otimes\mathbb{I}_{k}\\\
\mathbb{I}_{n}^{T}\otimes\mathbf{1}_{k}\end{matrix}\end{pmatrix}~{}~{}\text{and}$
$\displaystyle\mathbf{x}_{i}$ $\displaystyle\triangleq[x_{(i-1)\times
n+1},\dots,x_{i\times n}]^{T},~{}~{}\mathbf{y}_{i}\triangleq[y_{(i-1)\times
n+1},\dots,y_{i\times n}]^{T}~{}~{}\text{for all}~{}~{}i\in[|1,k|]$

Indeed the above optimization problem is just a reformulation of
$\text{LOT}_{k,c}(\mu,\mu)$ where we have vectorized the couplings in a
column-wise order. Let us now show the following lemma from which the result
will follow.

######  Lemma 2.

Under the same assumption of Proposition 9 we have that for all
$g\in\Delta_{k}^{*}$

$\displaystyle\min_{(\mathbf{x},\mathbf{y})\in
K_{a,g}^{2}}\sum_{i=1}^{k}\frac{\mathbf{x}_{i}^{T}C\mathbf{y}_{i}}{g_{i}}=\min_{\mathbf{x}\in
K_{a,g}}\sum_{i=1}^{k}\frac{\mathbf{x}_{i}^{T}C\mathbf{x}_{i}}{g_{i}}$

###### Proof.

Let $(\mathbf{x}^{*},\mathbf{y}^{*})$ solution of the LHS optimization
problem. Then we have that

$\displaystyle\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{x}^{*}_{i}}{g_{i}}$
$\displaystyle\geq\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}$
$\displaystyle\sum_{i=1}^{k}\frac{(\mathbf{y}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}$
$\displaystyle\geq\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}$

Therefore we obtain that

$\displaystyle 0$
$\displaystyle\leq\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{x}^{*}_{i}}{g_{i}}-\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}=\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C(\mathbf{x}^{*}_{i}-\mathbf{y}^{*}_{i})}{g_{i}}$
$\displaystyle 0$
$\displaystyle\leq\sum_{i=1}^{k}\frac{(\mathbf{y}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}-\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}=\sum_{i=1}^{k}\frac{(\mathbf{y}^{*}_{i}-\mathbf{x}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}$

Then by symmetry of $C$, we obtain by adding the two terms that

$\displaystyle\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i}-\mathbf{y}^{*}_{i})^{T}C(\mathbf{x}^{*}_{i}-\mathbf{y}^{*}_{i})}{g_{i}}\geq
0$

However, thanks to the linear constraints, we have that for all $i\in[|1,k|]$,

$\sum_{q=0}^{n-1}x^{*}_{(i-1)\times n+1+q}=\sum_{q=0}^{n-1}y^{*}_{(i-1)\times
n+1+q}=g_{i}$

Therefore $(\mathbf{x}_{i}^{*}-\mathbf{y}_{i}^{*})^{T}\mathbf{1}_{n}=0$ and
thanks to the negativity of the cost function $c$ we obtain that

$\displaystyle(\mathbf{x}_{i}^{*}-\mathbf{y}_{i}^{*})^{T}C(\mathbf{x}_{i}^{*}-\mathbf{y}_{i}^{*})\leq
0$

Therefore we have that

$\displaystyle(\mathbf{x}_{i}-\mathbf{y}_{i})^{T}C(\mathbf{x}_{i}-\mathbf{y}_{i})=0$

from which follows that

$\displaystyle\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{x}^{*}_{i}}{g_{i}}=\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}=\sum_{i=1}^{k}\frac{(\mathbf{y}^{*}_{i})^{T}C\mathbf{y}^{*}_{i}}{g_{i}}$

and the result follows. ∎

As the above result holds for any $g\in\Delta_{k}^{*}$, we obtain that

$\displaystyle\text{LOT}_{k,c}(\mu,\mu)=\min_{g\in\Delta_{k}^{*}}\min_{\mathbf{x}\in
K_{a,g}}\sum_{i=1}^{k}\frac{(\mathbf{x}^{*}_{i})^{T}C\mathbf{x}^{*}_{i}}{g_{i}}$

Then by formulating back this problem in term of matrices, we obtain that

$\displaystyle\text{LOT}_{k,c}(\mu,\mu)=\min_{g\in\Delta_{k}^{*}}\min_{Q\in\Pi_{a,g}}\langle
C,Q\text{diag}(1/g)Q^{T}\rangle$

from which the result follows. ∎

##  Appendix C Additional Experiments

###  C.1 Comparison of the $\gamma$ schedules

Figure 5: In this experiment, we compare two strategies for the choice of the
step-size in the MD scheme proposed by Scetbon et al. [2021] on two different
problems. More precisely, we compare the constant $\gamma$ schedule with the
proposed adaptive one and compare them when the distributions are sampled from
either uniform distributions (_left_) or mixtures of anisotropic Gaussians
(_right_). We show that the range of admissible $\gamma$ when considering a
constant schedule varies from one problem to another. Indeed, in the right
plot, we observe that the algorithm converges only when $\gamma\leq 1$, while
in the left plot, the algorithm manages to converge for $\gamma\leq 100$. We
also observe that our adaptive strategy allows to have a consistent choice of
admissible values for $\gamma$ whatever the problem considered. It is worth
noticing that whatever the $\gamma$ chosen, the algorithm converges towards
the same value, however the larger $\gamma$ is chosen in its admissible range,
the faster the algorithm converges.

###  C.2 Gradient Flows between two Moons

Figure 6: We compare the gradient flows $(\mu_{t})_{t\geq 0}$ (in red)
starting from a moon shape distribution, $\mu_{0}$, to another moon shape
distribution (in blue), $\nu$, in 2D when minimizing either
$L(\mu)\triangleq\text{DLOT}_{r,c}(\mu,\nu)$ or
$L(\mu)\triangleq\text{LOT}_{r,c}(\mu,\nu)$. The ground cost is the squared
Euclidean distance and we fix $r=100$. We consider 1000 samples from each
distribution and and we plot the evolution of the probability measure obtained
along the iterations of a gradient descent scheme. We also display in green
the vector field in the descent direction. We show that the debiased version
allows to recover the target distribution while $\text{LOT}_{r,c}$ is learning
a biased version with a low-rank structure.