# On the potential benefits of entropic regularization for smoothing Wasserstein estimators

Jérémie Bigot, Paul Freulon, Boris P. Hejblum, Arthur Leclaire

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jérémie Bigot, Paul Freulon, Boris P. Hejblum, Arthur Leclaire. On the potential benefits of entropic regularization for smoothing Wasserstein estimators. Electronic Journal of Statistics, 2025, 19 (2), pp.3867-3894. 10.1214/25-EJS2430. hal-03908229

HAL Id: hal-03908229

https://hal.science/hal-03908229v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the potential benefits of entropic regularization for smoothing Wasserstein estimators

Jérémie Bigot¹,², Paul Freulon *¹,²,⁵, Boris P. Hejblum¹,³,⁴, and Arthur Leclaire¹,²,⁶

¹Université de Bordeaux, Bordeaux, 33000, France.

²Institut de Mathématiques de Bordeaux et CNRS (UMR 5251), 33400 Talence, France.

³Bordeaux Population Health Research Center Inserm U1219, Inria SISTM, 33000 Bordeaux, France.

⁴Vaccine Research Institute (VRI), 94010 Créteil, France.

⁵EPFL, Institute of Mathematics, CH-1015 Lausanne, Switzerland.

⁶LTCI, Télécom Paris, IP Paris, 19 place Marguerite Perey, 91120 Palaiseau, France.

## Abstract

This paper focuses on the study of entropic regularization in optimal transport as a smooth- ing method for Wasserstein estimators, through the prism of the classical trade-off between approximation and estimation errors in statistics. Wasserstein estimators are defined as solu- tions of variational problems whose objective function involves the use of an optimal transport cost between probability measures. Such estimators can be regularized by replacing the opti- mal transport cost by its regularized version using an entropy penalty on the transport plan. The use of such a regularization has a potentially significant smoothing effect on the resulting estimators. In this work, we investigate its potential benefits on the approximation and esti- mation properties of regularized Wasserstein estimators. Our main contribution is to discuss how entropic regularization may reach, at a lower computational cost, statistical performances that are comparable to those of un-regularized Wasserstein estimators in statistical learning problems involving distributional data analysis. To this end, we present new theoretical re- sults on the convergence of regularized Wasserstein estimators. We also study their numerical performances using simulated and real data in the supervised learning problem of proportions estimation in mixture models using optimal transport.

## 1 Introduction

Wasserstein estimators are defined as solutions of variational problems whose objective function involves the use of an optimal transport (OT) cost between probability measures. Such estimators typically arise in statistical problems involving the minimization of a Wasserstein distance (or more generally an OT cost) between the empirical measure of the data and a distribution belonging to a parametric model (see [6]), and this class of estimators has found important applications in generative adversarial models for image processing (see e.g. [3]). Wasserstein estimators also represent an important class of inference methods in the field of statistical optimal transport for distributional data analysis where the observations at hand can be modeled as a set of histograms (see e.g. [7, 31, 33] for recent reviews).

Despite the appealing geometric properties of Wasserstein distances for comparing probability distributions, the computational burden required to evaluate an optimal transport cost is an im- portant limitation for its application in data analysis. The seminal paper [15] has opened a breach in the computational complexity of optimal transport by the addition of an entropic regularizing

<footnote>*correspondence: paul.freulon@epfl.ch</footnote>
term in the OT Kantorovich's formulation. In the last years, the benefit of this regularization has been to allow the use of OT based methods in statistics and machine learning with a time complexity that scales quadratically in the number of data using the Sinkhorn algorithm. This represents a significant improvement over the computational cost of un-regularized OT that scales cubically in the number of observations using linear programming. However, regularized OT has been mainly used so far as a fast numerical method to approximate un-regularized OT.

In this paper, we advocate the use of entropic regularization in computational OT as a smooth- ing method for un-regularized Wasserstein estimators. These estimators are obtained by replacing the standard OT cost in a variational problem by its entropy regularized version. The use of such a regularization has a beneficial smoothing effect on the resulting estimators as shown in [8] for the specific problem of computing a smooth Wasserstein barycenter from a set of discrete probability measures. In this paper, we investigate the impact of this smoothing effect of regularized Wasser- stein estimators through the prism of the tradeoff between approximation and estimation errors in statistics which is reminiscent of the classical bias versus variance tradeoff). Our main contribution is to discuss how entropic regularization yields estimators that may reach, at a lower computational cost, statistical performances that are comparable to those of un-regularized Wasserstein estima- tors in statistical learning problems involving distributional data analysis. To this end, we present new theoretical results on the convergence of regularized Wasserstein estimators. We also study their numerical performances using simulated and real data in the supervised learning problem of proportions estimation in mixture models using optimal transport.

1.1 Proportions estimation in mixture models using optimal transport

The motivation of this work comes from the active research field of automated analysis of flow cytometry measurements, see [1]. Flow cytometry is a high-throughput biotechnology used to characterize a large amount of m cells from a biological sample (with m ≥ 10 5 ) that produces a data set X$_{1}$, . . . , X$_{m}$ where each observation X$_{i}$ ∈ R d corresponds to a vector of d biomarkers of each single cell. Automated approaches in flow cytometry aim at clustering the data in order to estimate cellular population proportions in the biological sample. In [19], we have considered that such a data set can be represented with a discrete probability distribution 1 m ∑ m i $_{=1}$δ$_{X}$$_{i}$ with support in R $^{d}$, and we have introduced a new supervised algorithm based on regularized OT to estimate the different cell population proportions from a biological sample characterized with flow cytometry measurements. This approach optimally re-weights class proportions in a mixture model between a source data set (with known segmentation into cell sub-populations) to fit a target data set with unknown segmentation.

Most automated methods in flow cytometry cluster the observations [11]. However, the relevant clinical information is rather the class proportions, i.e. the cell population relative abundance. For instance, when monitoring the immune system of patient, clinicians are more interested by the proportion of CD4+ T-cells than to know whether each cell is a CD4+ T-cell.

To be more precise, let us denote by Y$_{1}$, . . . , Y$_{n}$ , the dataset from the target sample, and by X$_{1}$, . . . , X$_{m}$ the observations from the source biological sample. Thanks to the knowledge of a clustering of the source dataset into K classes C$_{1}$, . . . , C$_{K}$ , the empirical measure ˆ µ = 1 m ∑ m i $_{=1}$δ$_{X}$$_{i}$ can be decomposed as the following mixture of probability measures,

$$
\hat{\mu} = \sum_{k=1}^{K} \frac{m_k}{m} \left( \sum_{i:X_i \in C_k} \frac{1}{m_k} \delta_{X_i} \right) = \sum_{k=1}^{K} \hat{\pi}_k \hat{\mu}_k, \text{ where } \hat{\pi}_k = \frac{m_k}{m}, \quad (1.1)
$$

and each component ˆ µ$_{k}$ = ∑ i : X$_{i}$ ∈ C$_{k}$ 1 m$_{k}$ δ$_{X}$$_{i}$ corresponds to a known sub-population of cells with m$_{k}$ = # C$_{k}$ . Then, the method proposed in [19] aims at modifying the weights ( ˆ π$_{k}$ )$_{1}$$_{≤}$$_{k}$$_{≤}$$_{K}$ in such a way that the re-weighted source measure minimizes a regularized OT cost with respect to the target measure 1 n ∑ n j $_{=1}$δ$_{Y}$$_{j}$ . Then, the resulting weights yield an estimation of the proportions of sub-population of cells in the target sample. However, despite the efficiency of the method for the analysis of flow cytometry data, the work in [19] opens questions on the influence of the regularization, and we set to answer some of them in this work.
Let us now formalize the problem of proportions estimation in mixture models using regularized OT. We denote by $\mu = \sum_{k=1}^{K} \pi_k \mu_k$ a probability measure that can be decomposed as a mixture of K probability measures $\mu_1, \dots, \mu_K$. For $\theta \in \Sigma_K$, where $\Sigma_K = \{(\theta_1, \dots, \theta_K) \in \mathbb{R}_+^K : \sum_{k=1}^{K} \theta_k = 1\}$ is the K-dimensional simplex, we define $\mu_\theta$ as the re-weighted version of $\mu$ that is defined as

$$ \mu_{\theta} = \sum_{k=1}^{K} \theta_k \mu_k. \qquad (1.2) $$

Let $\nu$ be another probability measure. Proportions estimation in mixture models using OT is defined as the problem of finding $\theta^* \in \Sigma_K$ that minimizes an OT cost between $\mu_\theta$ and $\nu$. Denoting by $W_0(\mu, \nu)$ the un-regularized OT cost between $\mu$ and $\nu$ (we shall focus on the squared Wasserstein metric associated to the quadratic cost), the optimal vector of class proportions that we are targeting is:

$$ \theta^* \in \arg \min_{\theta \in \Sigma_K} W_0(\mu_\theta, \nu). $$

In practice, one only has access to independent samples from $\mu$ and $\nu$ denoted by $X_1, \dots, X_m$ (with a known clustering) and $Y_1, \dots, Y_n$ respectively. Therefore, estimators of $\theta^*$ will be obtained from the empirical versions of $\mu_\theta$ and $\nu$ denoted by

$$ \hat{\mu}_{\theta} = \sum_{k=1}^{K} \theta_k \hat{\mu}_k \quad \text{and} \quad \hat{\nu} = \frac{1}{n} \sum_{j=1}^{n} \delta_{Y_j}. $$

The computational cost to numerically evaluate $W_0(\mu, \nu)$ can be prohibitive, which led the author of [19] to consider its regularized version denoted by $W_\lambda(\mu, \nu)$ where $\lambda > 0$ represents the amount of entropic penalty that is put on the transport plan in the primal formulation of OT. Here, this regularized version of the OT cost is computed using the Sinkhorn algorithm, an iterative procedure whose convergence properties are now well understood (see [14] for a recent overview). However, after $\ell$ iterations of the Sinkhorn algorithm, it should be noted that one only has an approximation of the regularized OT cost that we will denote by $W_\lambda^{(\ell)}(\mu, \nu)$. In this work, we focus on the study of the convergence rate of the following estimator towards the optimal vector of class proportions $\theta^*$:

$$ \hat{\theta}_{\lambda}^{(\ell)} = \arg \min_{\theta \in \Sigma_K} W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (1.3) $$

This takes into account both the effect of entropic regularization and the influence of the number of iterations of the Sinkhorn algorithm. Our theoretical results shed some light on how the parameters $\lambda$ and $\ell$ influence the performance of the estimator $\hat{\theta}_\lambda^{(\ell)}$. We demonstrate the practical efficiency of our method and the impact of the regularization parameter $\lambda$ on simulated and real data (flow cytometry measurements).

## 1.2 Related works based on regularized optimal transport

Aside from the computational benefits of entropic regularization mentioned previously, recent developments have studied the statistical properties of a regularized OT cost computed from empirical measures. Indeed, in most cases, $\mu$ and $\nu$ are not available, and one has only access to their empirical versions $\hat{\mu}$ and $\hat{\nu}$ respectively built from $X_1, \dots, X_n$ sampled from $\mu$ and $Y_1, \dots, Y_n$ sampled from $\nu$. In this setting, it is natural to investigate the convergence rate of the plug-in estimator $W_0(\hat{\mu}, \hat{\nu})$ towards $W_0(\mu, \nu)$. This question is addressed in [18] where the authors proved that the resulting estimation error decays to zero at the rate $n^{-2/d}$ when using the quadratic cost in high dimension $d$. Due to its attractive computational efficiency, it is obviously interesting to examine the statistical efficiency of the regularized Wasserstein plug-in estimator naturally defined as $W_\lambda(\hat{\mu}, \hat{\nu})$. This issue as well as the approximation error induced by the regularization parameter is studied in [21]. These questions are thoroughly pursued in [14] as well as the effect of substituting $W_\lambda(\mu, \nu)$ by its debiased counterpart $S_\lambda(\mu, \nu)$. Since the first iteration of our work was posted
on ArXiv, some new results on the rate of convergence of $W_λ(\hat{\mu}, \hat{\nu})$ toward $W_λ(\mu, \nu)$ have been established. For instance, in [37, 25], some structural assumptions on the measures studied enable the authors to derive faster rate of convergence than what was previously known. Other works like [9], and more recently [24], study the limiting distribution of entropic optimal transport estimators. The OT loss functions $W_0, W_λ$ and $S_λ$ also constitute efficient tools for statistical estimation. For instance, a framework of parametric estimation where regularized OT acts as a loss function in learning problems is considered in [4]. Regularized Wasserstein losses are also considered in [23, 35, 28] for the design of generative models in image processing. In a more applied context, the use of regularized OT is investigated in [27, 19] to tackle estimation problems in biostatistics. The influence of the regularization parameter $λ$ for the purpose of computing smooth Wasserstein barycenters is also analyzed in [8, 13].

## 1.3 Organization of the paper

In Section 2 we recall the mathematical aspects of regularized OT needed to derive our results, and we detail the problem of optimal class proportions estimation in mixture models using OT. In Section 3, we introduce the various parametric Wasserstein estimators used to estimate the optimal class proportions. We also give the main results of this paper on a theoretical comparison of the convergence rates of regularized and un-regularized Wasserstein estimators. The influence of the number of iterations of the Sinkhorn algorithm on these convergence rates is also discussed. In Section 4, we propose an alternative bound for the estimation error studied. We also compare our results with related works. Section 5 is focused on numerical experiments that highlight the potential benefits of regularized Wasserstein estimators over un-regularized ones for appropriate choices of the entropic regularization parameter. Section 6 contains a conclusion and some perspectives. In the Appendix A, we detail the main arguments to obtain the convergence rates of regularized and un-regularized Wasserstein estimators.

# 2 Background on optimal transport and the problem of class proportions estimation

In this section, we introduce the notion of entropy regularized OT, and we present some of its mathematical properties needed to derive our results. Then, we describe the main application of this work on class proportions estimation in mixture models using OT. Finally, we discuss some identifiability issues in such models.

## 2.1 The OT problem and its regularized counterpart

**Notations.** In the whole paper, we will work in the space $\mathbb{R}^d$ equipped with the quadratic cost $c(x, y) = \|x - y\|^2$, where $\|x\| = \sqrt{\sum x_i^2}$ is the Euclidean norm. Let $\mathcal{X}$ and $\mathcal{Y}$ be two subsets of $\mathbb{R}^d$ that are assumed to be compact and included in $B(0, R) = \{x \in \mathbb{R}^d : \|x\| \le R\}$ throughout the paper. We denote by $\mathcal{M}_+(\mathcal{X})$ and $\mathcal{M}_+(\mathcal{Y})$ the sets of probability measures on $\mathcal{X}$ and $\mathcal{Y}$ respectively. For $Y_1, \dots, Y_n \sim \nu$, we denote by $\hat{\nu}$ the empirical counterpart of $\nu$ defined by $\hat{\nu} = \frac{1}{n} \sum_{i=1}^n \delta_{Y_i}$. The notation $\lesssim$ means inequality up to a multiplicative universal constant. For $\mu \in \mathcal{M}_+(\mathcal{X})$ and $\nu \in \mathcal{M}_+(\mathcal{Y})$, we let $\Pi(\mu, \nu)$ be the set of probability measures on $\mathcal{X} \times \mathcal{Y}$ with marginals $\mu$ and $\nu$. The problem of entropic optimal transportation between $\mu \in \mathcal{M}_+(\mathcal{X})$ and $\nu \in \mathcal{M}_+(\mathcal{Y})$ is then defined as follows.

**Definition 2.1 (Primal formulation).** For any $(\mu, \nu) \in \mathcal{M}_+(\mathcal{X}) \times \mathcal{M}_+(\mathcal{Y})$, the Kantorovich formulation of the regularized optimal transport between $\mu$ and $\nu$ is the following convex minimization problem

$$W_{\lambda}(\mu, \nu) = \min_{\pi \in \Pi(\mu, \nu)} \int_{\mathcal{X} \times \mathcal{Y}} \|x - y\|^2 d\pi(x, y) + \lambda \mathrm{KL}(\pi|\mu \otimes \nu), \quad (2.1)$$
where $\|x - y\|$ is the Euclidean distance between $x$ and $y$, $\lambda \ge 0$ is a regularization parameter, and KL stands for the Kullback-Leibler divergence, between $\pi$ and a positive measure $\xi$ on $\mathcal{X} \times \mathcal{Y}$, up to the additive term $\int_{\mathcal{X} \times \mathcal{Y}} d\xi(x, y)$, namely

$$KL(\pi|\xi) = \int_{\mathcal{X} \times \mathcal{Y}} \log\left(\frac{d\pi}{d\xi}(x, y)\right) d\pi(x, y),$$

in the case $\pi$ absolutely continuous w.r.t. $\xi$, otherwise $KL(\pi|\xi) = +\infty$. For $\lambda = 0$, the quantity $W_0(\mu, \nu)$ is the standard (un-regularized) OT cost, and for $\lambda > 0$, we refer to $W_\lambda(\mu, \nu)$ as the regularized OT cost between $\mu$ and $\nu$. Note that the continuity of $c$ and the compactness of $\mathcal{X}$ and $\mathcal{Y}$ imply that $W_\lambda(\mu, \nu)$ is finite for any value of $\lambda \ge 0$. Let us now introduce the dual and semi-dual formulations (see e.g. [36, 22]) of the minimization problem (2.1).

**Theorem 2.1 (Dual formulation).** *Strong duality holds for the primal problem (2.1) in the sense that*

$$W_{\lambda}(\mu, \nu) = \sup_{\substack{\varphi \in L^{\infty}(\mathcal{X}), \\ \psi \in L^{\infty}(\mathcal{X})}} \int \varphi(x)d\mu(x) + \int \psi(y)d\nu(y) \\ - \int m_{\lambda}(\varphi(x) + \psi(y) - \|x - y\|^2)d\mu(x)d\nu(y) \qquad (2.2)$$

where $L^\infty(\mathcal{X})$ denotes the space of essentially bounded functions quotiented by a.e. equality, and

$$m_{\lambda}(t) = \begin{cases} +\infty & \text{if } \lambda = 0 \\ \lambda(e^{\frac{t}{\lambda}} - 1) & \text{if } \lambda > 0 \end{cases}$$

A solution $(\varphi, \psi)$ of the dual problem (2.2) is called a pair of Kantorovich potentials. Besides, since $\mathcal{X}, \mathcal{Y}$ are compact and $c$ is continuous, it follows that the dual problem admits a solution $(\varphi, \psi) \in \mathcal{C}_b(\mathcal{X}) \times \mathcal{C}_b(\mathcal{Y})$. Moreover, when $\lambda > 0$, there exists solutions $\varphi, \psi$ to the dual problem (2.2) which are uniquely defined almost everywhere, up to an additive constant. The solutions of this regularized dual problem have the specific structure of $c$-transform functions. For the quadratic cost $c(x, y) = \|x - y\|^2$, the regularized $c$-transform is defined as in [17]: for $\lambda > 0$, we set

$$\forall x \in \mathbb{R}^d, \quad \psi_{\nu}^{c, \lambda}(x) = -\lambda \log \int e^{-\frac{\|x-y\|^2}{\lambda}} \psi(y) d\nu(y), \qquad (2.3)$$

and for $\lambda = 0$, the $c$-transform simply reads

$$\forall x \in \mathbb{R}^d, \quad \psi^c(x) = \min_{y \in \mathcal{Y}} \{\|x - y\|^2 - \psi(y)\}. \qquad (2.4)$$

We also define the analogous operators for the $y$-variable (and for simplicity, we use the same notation for $c$-transform of $x$-functions or $y$-functions). Notice that the operation used in (2.3) can be understood as a smoothed minimum that depends on $\nu$. Therefore, when $\lambda > 0$ we will stick to the notation $\psi_{\nu}^{c, \lambda}$ to keep in mind the possible dependence on $\nu$ of the regularized $c$-transform. Notice also that, even if $\psi_{\nu}^{c, \lambda}$ will be integrated only on $\mathcal{X}$, the formulae allow to extrapolate the $c$-transform on the whole space $\mathbb{R}^d$. In the sequel of this paper we extrapolate the $c$-transform on $B(0, R)$ to manipulate functions defined on a convex subset of $\mathbb{R}^d$ without imposing the convexity of $\mathcal{X}$ and $\mathcal{Y}$. Moreover, considering the $c$-transform amounts to optimizing one of the two potential, thus leading to an optimization problem with respect to one single potential, called the semi-dual problem [22]:

**Corollary 2.1 (Semi-dual formulation).** *The dual problem (2.2) is equivalent to the following semi-dual problem in the sense that*

$$W_{\lambda}(\mu, \nu) = \sup_{\psi \in L^{\infty}(\mathcal{Y})} \int \psi_{\nu}^{c, \lambda}(x) d\mu(x) + \int \psi(y) d\nu(y). \qquad (2.5)$$
A solution $\psi$ of the semi-dual problem is called a Kantorovich potential. In other words, $\psi$ is a Kantorovich potential if and only if $(\psi_{\nu}^{c,\lambda}, \psi)$ is a pair of Kantorovich potentials. By symmetry, we can also formulate a semi-dual problem on the dual variable $\varphi$. For discrete probability distributions, the iterative Sinkhorn algorithm, as defined below, (see e.g. [15]) allows to approximate the regularized OT cost $W_{\lambda}(\mu, \nu)$ as follows.

**Definition 2.2.** For $\mu = \sum_{i=1}^{n} \mu_i \delta_{x_i}$ and $\nu = \sum_{j=1}^{m} \nu_j \delta_{y_j}$ two discrete distributions on $\mathbb{R}^d$, the approximation of the regularized OT cost returned by the Sinkhorn approximation after $\ell$ iterations equals

$$W_{\lambda}^{(\ell)}(\mu, \nu) = \sum_{i=1}^{n} \mu_i \varphi_i^{(\ell)} + \sum_{j=1}^{m} \nu_j \psi_j^{(\ell)}. \qquad (2.6)$$

The variables $\varphi^{(\ell)}$ and $\psi^{(\ell)}$ being the dual variables returned after $\ell$ iterations of the Sinkhorn algorithm. Starting from $\psi^0 = 0_m \in \mathbb{R}^m$, the Sinkhorn $\ell^{th}$ iteration is defined by the update of the dual variables:

$$\begin{aligned} \varphi_i^{(\ell)} &= -\lambda \log \left( \sum_{j=1}^{m} \exp \left( -\frac{\|x_i - y_j\|^2 - \psi_j^{(\ell-1)}}{\lambda} \right) \nu_j \right) \\ \psi_j^{(\ell)} &= -\lambda \log \left( \sum_{i=1}^{n} \exp \left( -\frac{\|x_i - y_j\|^2 - \varphi_i^{(\ell)}}{\lambda} \right) \mu_i \right). \end{aligned} \qquad (2.7)$$

**Remark 2.1** (Convergence of Sinkhorn algorithm). Convergence guarantees of Sinkhorn algorithm are established for instance in [14, Prop.2]. That is, when the number of iterations $\ell$ goes to infinity, $W_{\lambda}^{(\ell)}(\mu, \nu)$ converges toward $W_{\lambda}(\mu, \nu)$.

A de-biased version of the regularized OT cost is also applied in [23] and further studied in [17] and [14]. This regularized OT cost is referred to as the Sinkhorn divergence, and defined as follows.

**Definition 2.3** (Sinkhorn Divergence [17]). For $\lambda > 0$, the Sinkhorn divergence between two probability measures $\mu$ and $\nu$ is defined by the formula

$$S_{\lambda}(\mu, \nu) := W_{\lambda}(\mu, \nu) - \frac{1}{2} (W_{\lambda}(\mu, \mu) + W_{\lambda}(\nu, \nu)). \qquad (2.8)$$

In this article, we focus on the classical regularized transport cost $W_{\lambda}$. At the price of stronger assumptions, we also extend our results to the Sinkhorn divergence $S_{\lambda}$ in the Appendix D. In the experimental Section 5, we will study the performance of the estimators based on the Sinkhorn divergence.

## 2.2 An alternative dual problem

We now introduce an alternative dual formulation of regularized OT that is specific to the quadratic cost. This alternative dual problem is restricted to a class of Kantorovich potentials that are concave and Lipschitz functions, which proves useful to derive some of the convergence rates given in Section 3. The relation between those dual problems has already been explicit for un-regularized OT (for example in [14]), and we extend it to the regularized case. Let $\lambda \ge 0$. By expanding the squared Euclidean cost, we have for any $\pi \in \Pi(\mu, \nu)$,

$$\begin{aligned} \int_{\mathcal{X} \times \mathcal{Y}} \|x - y\|^2 d\pi(x, y) + \lambda \text{KL}(\pi|\mu \otimes \nu) &= \int_{\mathcal{X}} \|x\|^2 d\mu(x) + \int_{\mathcal{Y}} \|y\|^2 d\nu(y) \\ &\quad - 2 \int_{\mathcal{X} \times \mathcal{Y}} \langle x, y \rangle d\pi(x, y) + \lambda \text{KL}(\pi|\mu \otimes \nu). \end{aligned} \qquad (2.9)$$

The above decomposition leads us to consider the new regularized transportation problem

$$W_{\lambda}^{s}(\mu, \nu) = \min_{\pi \in \Pi(\mu, \nu)} \int_{\mathcal{X} \times \mathcal{Y}} s(x, y) d\pi(x, y) + \lambda \text{KL}(\pi|\mu \otimes \nu), \qquad (2.10)$$
with $s(x, y) = -2\langle x, y \rangle$. First, we remark that the standard regularized Wasserstein distance $W_\lambda(\mu, \nu)$ and the alternative regularized Wasserstein distance $W_\lambda^s(\mu, \nu)$ are link through the relation

$$W_\lambda(\mu, \nu) = \int_X \|x\|^2 d\mu(x) + \int_Y \|y\|^2 d\nu(y) + W_\lambda^s(\mu, \nu). \quad (2.11)$$

A dual formulation associated to the problem (2.10) is given by the next proposition.

**Proposition 2.1.** The dual problem associated to (2.10) writes as

$$W_\lambda^s(\mu, \nu) = \sup_{\substack{\varphi \in L^\infty(\mathcal{X}) \\ \psi \in L^\infty(\mathcal{Y})}} \int \varphi(x)d\mu(x) + \int \psi(y)d\nu(y) \\ - \int m_\lambda(\varphi(x) + \psi(y) + 2\langle x, y \rangle)d\mu(x)d\nu(y), \quad (2.12)$$

where $m_\lambda$ is defined in Theorem 2.1.

*Proof.* The key argument is to remark that (2.10) is a regularized optimal transportation problem with cost function $s(x, y) = -2\langle x, y \rangle$. Hence, as $\mathcal{X}$ and $\mathcal{Y}$ are assumed to be compact and $s$ continuous, it follows that strong duality holds (see e.g. [22, 5]) in the sense of equation (2.12). $\square$

Fort the cost function $s(x, y) = -2\langle x, y \rangle$, we can also define a $s$-transform and a semi-dual problem as follows. For the cost $s(x, y) = -2\langle x, y \rangle$ and for $\varphi \in L^\infty(\mathcal{X})$ the $s$-transform is defined as

$$\forall y \in \mathbb{R}^d, \quad \varphi_\mu^{s,\lambda}(y) = \begin{cases} -\lambda \log \left( \int \exp \left( \frac{\varphi(x) + 2\langle x, y \rangle}{\lambda} \right) d\mu(x) \right), & \text{for } \lambda > 0, \\ -\max_{x \in \mathcal{X}} (\varphi(x) + 2\langle x, y \rangle) & \text{for } \lambda = 0. \end{cases} \quad (2.13)$$

By the above $s$-transform in the dual problem (2.12) we obtain the following semi-dual formulation

$$\sup_{\varphi \in L^\infty(\mathcal{X})} \int \varphi(x) d\mu(x) + \int \varphi_\mu^{s,\lambda}(y) d\nu(y). \quad (2.14)$$

We conclude this section by studying some properties of this $s$-transform. While already established in [12], we give an elementary proof for completeness.

**Proposition 2.2.** For $\lambda \ge 0$, the $s$-transform $\varphi_\mu^{s,\lambda}$ is concave and R-Lipschitz on $B(0, R)$.

*Proof.* We start with the concavity of $\varphi_\mu^{s,\lambda}$. For $\lambda = 0$, it follows from the fact that a maximum of convex functions is convex. Now, for $\lambda > 0$, $y_1, y_2 \in \mathcal{Y}$ and $t \in (0, 1)$, we have

$$\int_{\mathcal{X}} \exp \left( \frac{\varphi(x) + \langle x, ty_1 + (1-t)y_2 \rangle}{\lambda} \right) d\mu(x) \\ = \int_{\mathcal{X}} \exp \left( t \frac{\varphi(x) + 2\langle x, y_1 \rangle}{\lambda} \right) \exp \left( (1-t) \frac{\varphi(x) + 2\langle x, y_2 \rangle}{\lambda} \right) d\mu(x) \\ \le \left( \int_{\mathcal{X}} \exp \left( \frac{\varphi(x) + 2\langle x, y_1 \rangle}{\lambda} \right) d\mu(x) \right)^t \left( \int_{\mathcal{X}} \exp \left( \frac{\varphi(x) + 2\langle x, y_2 \rangle}{\lambda} \right) d\mu(x) \right)^{1-t},$$

thanks to Hölder inequality with exponents $p = 1/t$ and $p' = 1/(1-t)$. Applying $-\lambda \log$ on both sides gives directly

$$\varphi_{\mu}^{s,\lambda} (ty_1 + (1-t)y_2) \geq t \varphi_{\mu}^{s,\lambda} (y_1) + (1-t) \varphi_{\mu}^{s,\lambda} (y_2). \quad (2.15)$$

Now, we will see as in [17] that the regularized $s$-transform inherits the Lipschitz constant of the cost. For $y_1, y_2 \in \mathcal{Y}$ and $x \in \mathcal{X}$, we have $|\langle x, y_1 - y_2 \rangle| \le R \|y_1 - y_2\|$ thanks to Cauchy-Schwarz inequality, and thus

$$\varphi(x) + 2\langle x, y_1 \rangle \le R \|y_1 - y_2\| + \varphi(x) + 2\langle x, y_2 \rangle.$$
Taking $\lambda \log \int_X \exp(\dot{\lambda})d\mu(x)$ for $\lambda > 0$ (resp. $\max_{x \in X}$ for $\lambda = 0$) on both sides gives

$$ \varphi_{\mu}^{s,\lambda}(y_2) \leq R \|y_1 - y_2\| + \varphi_{\mu}^{s,\lambda}(y_1). $$

By symmetry, we get $|\varphi_{\mu}^{s,\lambda}(y_1) - \varphi_{\mu}^{s,\lambda}(y_2)| \leq R\|y_1 - y_2\|$.

□

## 2.3 Definition of the problem and quantity of interest

Let $\mu = \sum_{k=1}^K \pi_k \mu_k$ be a probability measure that can be decomposed as a mixture of K probability measures $\mu_1, \dots, \mu_K$ in $M_+^1(\mathcal{X})$. For $\theta \in \Sigma_K = \{(\theta_1, \dots, \theta_K) \in \mathbb{R}_+^K : \sum_{k=1}^K \theta_k = 1\}$, the re-weighted version of $\mu$ is defined as

$$ \mu_{\theta} = \sum_{k=1}^{K} \theta_{k} \mu_{k}. \qquad (2.16) $$

Let $\nu$ be another probability measure in $M_+^1(\mathcal{Y})$ referred to as the target distribution. The problem of class proportions estimation consists in estimating an optimal weighting vector

$$ \theta^* \in \arg \min_{\theta \in \Sigma_K} W_0(\mu_{\theta}, \nu) \qquad (2.17) $$

from empirical versions of the $\mu_1, \dots, \mu_K$ and $\nu$. In what follows, we discuss some properties of the optimisation problem (2.17).

First, this minimization problem is motivated by the implicit assumption that representing the target measure $\nu$ as a mixture of K probability measures is relevant. To illustrate this point, we first state a result showing that one can recover the true class proportions in the ideal setting where the target distribution $\nu$ is also a mixture of $\mu_1, \dots, \mu_K$.

**Lemma 2.1.** Suppose that $\mu_{\theta}$ and $\nu$ are mixtures of probability measures with the same components $\mu_1, \dots, \mu_K$ but with different class proportions, respectively denoted by $\theta \in \Sigma_K$ and by $\tau \in \Sigma_K$. If the model $\{\mu_{\theta} = \sum_{k=1}^K \theta_k \mu_k \mid \theta \in \Sigma_K\}$ is identifiable (in the sense that the mapping $\theta \mapsto \mu_{\theta}$ is injective), then the solution of optimization problem (2.17) is unique and one has that $\theta^* = \tau$.

*Proof.* The non-negativity property of $W_0$ ensures that for all $\theta \in \Sigma_K$, $W_0(\mu_{\theta}, \nu) \geq 0$. Next, for $\theta \in \Sigma_K$,

$$ W_0(\mu_{\theta}, \nu) = 0 \Leftrightarrow \sum_{k=1}^{K} \theta_k \mu_k = \sum_{k=1}^{K} \tau_k \mu_k \\ \Leftrightarrow \theta = \tau, $$

where the last equivalence comes from the assumption that the model $\{\mu_{\theta} \mid \theta \in \Sigma_K\}$ is identifiable. From this result, we deduce $\arg \min_{\theta \in \Sigma_K} W_0(\mu_{\theta}, \nu) = \{\tau\}$. □

Notice that the injectivity of $\theta \mapsto \mu_{\theta}$ relates to the affine independence of $\{\mu_1, \dots, \mu_K\}$. It is satisfied for example when the measures $\mu_1, \dots, \mu_K$ have disjoint supports. If all the scenarios are not as friendly as the one considered in Lemma 2.1, in numerous applications (for instance when the data can be clustered into sub-populations), it is relevant to approximate $\nu$ by a mixture model. The next result is about the smoothness of the minimization problem (2.17).

**Lemma 2.2.** Suppose that $\mu_{\theta}$ is defined as in (1.2). Then, the function

$$ F: \left\{ \begin{array}{rcl} \Sigma_K & \rightarrow & \mathbb{R}_+ \\ \theta & \mapsto & W_0(\mu_{\theta}, \nu) \end{array} \right. \quad \text{is continuous on } \Sigma_K. $$
Proof. Let $\theta \in \Sigma_K$ and $(\theta^{(n)})$ a sequence in $\Sigma_K$ that converges to $\theta$. Then, the probability sequence $(\mu_{\theta^{(n)}})$ converges weakly toward $\mu_{\theta}$. Indeed, for any bounded continuous function $\varphi$, one has that $\int \varphi d\mu_{\theta^{(n)}} = \sum_{k=1}^K \theta_k^{(n)} \int \varphi d\mu_k$. As $\theta^{(n)} \to_{n \to \infty} \theta$, it follows that

$$ \sum_{k=1}^{K} \theta_{k}^{(n)} \int \varphi d\mu_{k} \to \sum_{k=1}^{K} \theta_{k} \int \varphi d\mu_{k}. = \int \varphi d\mu_{\theta}. $$

That is, $(\mu_{\theta^{(n)}})$ weakly converges towards $\mu_{\theta}$. As we work under the assumption that $\mathcal{X}$ and $\mathcal{Y}$ are compact, weak convergence is equivalent to Wasserstein convergence [36][Thm. 5.10]. Hence, $W_0(\mu_{\theta^{(n)}}, \nu)$ converges to $W_0(\mu_{\theta}, \nu)$ when $n$ goes to infinity; which shows the continuity of $F$. $\square$

Since the set $\Sigma_K$ is compact, the existence of a minimizer of the optimization problem (2.17) follows from Lemma 2.2. We now give sufficient conditions that ensure the strict convexity of the objective function $\theta \mapsto W_0(\mu_{\theta}, \nu)$.

**Lemma 2.3.** Assume that $\nu$ is absolutely continuous with respect to the Lebesgue measure. Then, if the model $\{\mu_{\theta} | \theta \in \Sigma_K\}$ is identifiable (in the sense that the mapping $\theta \mapsto \mu_{\theta}$ is injective), the function $F: \begin{cases} \Sigma_K \to \mathbb{R}_+ \\ \theta \mapsto W_0(\mu_{\theta}, \nu) \end{cases}$ is strictly convex.

Proof. Thanks to the assumption that $\nu$ is absolutely continuous, Proposition 7.19 in [36] ensures the strict convexity of the functional $\mu \mapsto W_0(\mu, \nu)$. Let $\theta_0, \theta_1 \in \Sigma_K$ with $\theta_0 \neq \theta_1$ and $t \in (0, 1)$. Then, we have that $F(t\theta_0 + (1-t)\theta_1) = W_0(\mu_{t\theta_0+(1-t)\theta_1}, \nu)$, and $\mu_{t\theta_0+(1-t)\theta_1} = t\mu_{\theta_0} + (1-t)\mu_{\theta_1}$.

Since $\theta_0 \neq \theta_1$ and the model $\{\mu_{\theta} | \theta \in \Sigma_K\}$ is supposed to be identifiable, we have that $\mu_{\theta_0} \neq \mu_{\theta_1}$. Therefore, the strict convexity of $\mu \mapsto W_0(\mu, \nu)$ yields

$$ W_0(t\mu_{\theta_0} + (1-t)\mu_{\theta_1}, \nu) < tW_0(\mu_{\theta_0}, \nu) + (1-t)W_0(\mu_{\theta_1}, \nu). $$

which proves the strict convexity of $F: \theta \mapsto W_0(\mu_{\theta}, \nu)$. $\square$

# 3 Convergence rates for the expected excess risk of parametric Wasserstein estimators

In this section, we present the regularized and un-regularized parametric Wasserstein estimators that are considered in this paper, and we compare their convergence rates.

## 3.1 Definition of the estimators

We aim at estimating $\theta^*$ when the distributions $\mu$ and $\nu$ are only observed through samples. Hence, we assume given the following empirical measures (as defined in Section 1.1)

$$ \hat{\mu} = \sum_{k=1}^{K} \hat{\pi}_k \hat{\mu}_k, \text{ where } \hat{\pi}_k = \frac{m_k}{m}, \quad \text{and} \quad \hat{\nu} = \frac{1}{n} \sum_{j=1}^{n} \delta_{Y_j}, $$

where each component $\hat{\mu}_k$ corresponds to a known sub-population of cells of size $m_k$ in the source sample $X_1, \dots, X_m$.

Moreover, we recall that $\hat{\mu}_{\theta} = \sum_{k=1}^{K} \theta_k \hat{\mu}_k$. denotes the empirical version of the re-weighted measure $\mu_{\theta}$.

We can now define the various Wasserstein estimators whose convergence properties are discussed in Section 3.2. Depending on the regularization parameter chosen, and using the empirical measures $\hat{\mu}_{\theta}$ and $\hat{\nu}$, a family of estimators $(\hat{\theta}_{\lambda})_{\lambda \ge 0}$ of the class proportions can be defined as follows:

$$ \text{for } \lambda \ge 0, \quad \hat{\theta}_{\lambda} \in \hat{\Theta}_{\lambda} := \arg \min_{\theta \in \Sigma_K} W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}), \qquad (3.1) $$
When considering entropy regularized OT, that is when $\lambda > 0$, we also propose to study the estimators that are obtained with the Sinkhorn algorithm on the sample distributions after a limited number of $\ell$ iterations, that are

$$ \hat{\theta}_{\lambda}^{(\ell)} \in \hat{\Theta}_{\lambda}^{(\ell)} := \arg \min_{\theta \in \Sigma_K} W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (3.2) $$

As pointed previously in Remark 2.1, the convergence of $W_{\lambda}^{(\ell)}$ towards $W_{\lambda}$ allows to interpret the estimator $\hat{\theta}_{\lambda}$ as a limiting case of $\hat{\theta}_{\lambda}^{(\ell)}$ when $\ell$ goes to infinity. Beside studying the estimators $\hat{\theta}_{\lambda}$ and $\hat{\theta}_{\lambda}^{(\ell)}$; we also extend our result when substituting the regularized transport cost $W_{\lambda}$ in equation (3.1) by the Sinkhorn divergence $S_{\lambda}$ that is defined by formula (2.8). Due to space constraint, we have gathered theoretical results related to the Sinkhorn divergence to the Appendix D. In this paper, to assess the performance of a given estimator $\hat{\theta}$ of $\theta^{*}$ based on $n$ samples, we shall consider the following expected excess risk defined as

$$ r_n(\mu_{\hat{\theta}}, \nu) = \mathbb{E}[W_0(\mu_{\hat{\theta}}, \nu) - W_0(\mu_{\theta^{*}}, \nu)]. \qquad (3.3) $$

**Remark 3.1.** In our context of parametric Wasserstein estimation, we can interpret the excess risk as the representation error of $\nu$ induced by the estimator. Indeed, $\mu_{\theta^{*}}$ defined in equation (2.17) is the best representation of $\nu$ in the model $\{\mu_{\theta} | \theta \in \Sigma_K\}$ w.r.t. the Wasserstein distance. And, $W_0^{1/2}$ being a distance, under the assumption that the function $\theta \mapsto W_0(\mu_{\theta}, \nu)$ is bounded on $\Sigma_K$, we can write

$$ 0 \le W_0(\mu_{\hat{\theta}}, \nu) - W_0(\mu_{\theta^{*}}, \nu) \lesssim W_0^{1/2}(\mu_{\hat{\theta}}, \mu_{\theta^{*}}). $$

This equation shows that the excess risk is closely related to Wasserstein distance between the best representation of $\nu$ in the model that is $\mu_{\theta^{*}}$ and its estimated version $\mu_{\hat{\theta}}$.

**Remark 3.2.** Instead of controlling the excess risk (3.3) we would have preferred to work directly on the weights. That is, upper bounding the quantity $\|\hat{\theta} - \theta^{*}\|$. It would have been possible to derive such a result if the function $\theta \mapsto W_0(\mu_{\theta}, \nu)$ had been strongly convex. However, we can find elementary counter-examples where $\theta \mapsto W_0(\mu_{\theta}, \nu)$ is not strongly convex. Indeed, on the real line, let us consider $\mu_{\theta} = \theta\delta_0 + (1 - \theta)\delta_1$ and $\nu = (\delta_0 + \delta_1)/2$. Then one can show that for every $\theta \in [0, 1]$,

$$ W_0(\mu_{\theta}, \nu) = |1/2 - \theta|, $$

which is not strongly convex. This result can be established thanks to the formula that links the quantile functions to the optimal transport cost on the real line (see e.g., [36, Prop. 2.17])

**Remark 3.3.** We believe that our results could be extended to other parametric families of the form $\{\mu_{\theta}, \theta \in \Theta\}$, provided compact supports assumptions are satisfied, as well as continuity assumptions on the map $\theta \mapsto \mu_{\theta}$. For instance, further research could aim to extend existing results by [20] on the statistical analysis of un-regularized Wasserstein Generative Adversarial Networks (WGAN) to the case of entropy regularized WGAN considered by [35].

In Section 3.2, we present upper bounds on the above expected risk for the proposed estimators. When the regularization parameter $\lambda$ is involved, we also propose a decreasing choice $(\lambda_n)_{n \ge 0}$ of its value to ensure that the resulting estimator has an expected excess risk that goes to zero when $n \to +\infty$.

## 3.2 Convergence rates for the expected excess risk

This section contains the main results of this paper. We study the rate of convergence of the family of estimators $(\hat{\theta}_{\lambda}^{(\ell)})_{\lambda \ge 0}$ depending on the parameters $\lambda \ge 0$ and $\ell \in \mathbb{N}^{*}$. In the following results, the notation $\lesssim$ means inequality up to a multiplicative universal constant.
As classically done in nonparametric statistics, we decompose the excess risk of an estimator into an estimation error and an approximation error that need to be balanced to derive an optimal choice of the regularization parameter $\lambda \to 0$ as the number of observations tends to infinity. For example, the excess risk of the estimator $\hat{\theta}_{\lambda}$ defined in equation (3.1) is upper bounded as follows:

$$ W_0(\mu_{\hat{\theta}_{\lambda}}, \nu) - W_0(\mu_{\theta^*}, \nu) \leq 2 \sup_{\theta \in \Sigma_K} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \nu)| + 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \nu)|. \quad (3.4) $$

As the introduction of entropic penalty term in the optimal transport problem was motivated by computational improvement [15], it is also useful to take into account the algorithmic error. Therefore, we substitute in equation (3.4) the estimator $\hat{\theta}_{\lambda}$ by its version computed with Sinkhorn algorithm:

$$ \hat{\theta}_{\lambda}^{(\ell)} = \arg \min_{\theta \in \Sigma_K} W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}). $$

In such a case, we provide an upper bound in the next lemma.

**Lemma 3.1.** The excess risk of the estimator $\hat{\theta}_{\lambda}^{(\ell)}$ is upper bounded as follows:

$$ \begin{aligned} 0 \le W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) \le 2 \sup_{\theta \in \Sigma_K} & |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \nu)| \\ & + 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \nu)| \\ & + 2 \sup_{\theta \in \Sigma_K} |W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \nu)|. \end{aligned} \quad (3.5) $$

The computations leading to Lemma 3.1 are gathered in Section A.1 of the Appendix. The main theorem of this article is based on a new bound for the control of the estimation error, given in the following proposition.

**Proposition 3.1.** Set $\lambda \ge 0$, and suppose that the probability measures $\mu_1, \dots, \mu_K$, and $\nu$ have compact supports included in $B(0, R)$. If for all components $\mu_k$ as well as for $\nu$ at least $n$ observations are available, Then the following inequality holds true:

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \nu)| \right] \lesssim K \mathcal{E}(n, d), \quad (3.6) $$

where the upper bound $\mathcal{E}(n, d)$ is defined by

$$ \mathcal{E}(n, d) := \begin{cases} R^2 n^{-1/2} & \text{if } d < 4, \\ R^2 n^{-1/2} \log(n) & \text{if } d = 4, \\ R^2 n^{-2/d} & \text{if } d > 4. \end{cases} \quad (3.7) $$

The proof of Proposition 3.1 is deferred to Section A.2 of the Appendix. This proof relies on [14, Lemma 4] where the maximum of an empirical process is upper bounded. We point out that the upper bound in Proposition 3.1 is independent of the regularization parameter $\lambda$.

**Theorem 3.1.** Set $\lambda > 0$ and suppose that all probability measures $\mu_1, \dots, \mu_K$ and $\nu$ have compact supports. Assume that for all the components $\mu_k$ as well as for $\nu$, at least $n$ observations are available. Then, the expected excess risk of the estimator $\hat{\theta}_{\lambda}^{(\ell)}$ introduced in equation (3.2) is upper bounded as follows:

$$ \mathbb{E} \left[ W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) \right] \lesssim \mathcal{E}(n, d) + \lambda \log \left( \frac{1}{\sqrt{d\lambda}} \right) + \frac{1}{\lambda \ell}, $$
where the quantity $\mathcal{E}(n, d)$ is defined by formula (3.7), and the implicit multiplicative constant depends only on $K$, $d$ and $R$.

The detailed proof of Theorem 3.1 can be found in Section A.3 of the Appendix. We mention that it is based on the upper bound (3.5) where each term of the right-hand side is controlled by the appropriate bound. The expected estimation error is upper bounded thanks to Proposition 3.1. For the remaining terms, we collect results established in the literature. More precisely, we exploit the works of Genevay et al. [21] and of Chizat et al. [14] to control the approximation and the algorithm errors.

**Corollary 3.1.** Suppose that every probability measure has compact support and that $d > 4$. If for all the components $\mu_k$ as well as for $\nu$, at least $n$ observations are available, then the following non-asymptotic rates of convergence hold for the estimator $\hat{\theta}_{\lambda_n}^{(\ell_n)}$ computed with the Sinkhorn algorithm:

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda_n}^{(\ell_n)}, \nu}) - W_0(\mu_{\theta^*}, \nu)] \lesssim n^{-\frac{2}{d}} \log(n), \quad \text{with} \quad \begin{cases} \lambda_n = n^{-2/d}, \\ \ell_n \ge 32R^4 n^{4/d}. \end{cases} $$

*Proof.* Setting the regularization parameter to $\lambda = n^{-2/d}$, and the number of iterations to $\ell = 32R^4 n^{4/d}$ in inequality (3.2) yields the announced rate of convergence. $\square$

**Remark 3.4 (Extension to the Sinkhorn divergence $S_\lambda$).** The estimators analyzed in Theorem 3.1 and Corollary 3.1 are defined as solutions of variational problems based on the regularized transport cost $W_\lambda$. Under ad hoc assumptions, our results can be extended to the Sinkhorn divergence $S_\lambda$ whose definition is reminded in equation (2.8). In such a case, we define the collection of estimators $(\hat{\theta}_\lambda^S)_{\lambda>0}$ as follows:

$$ \hat{\theta}_{\lambda}^{S} \in \arg \min_{\theta \in \Sigma_K} S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (3.8) $$

Provided stronger assumptions are made, the approximation error of the Sinkhorn divergence is smaller than the approximation error $|W_\lambda - W_0|$. However, the estimators $\hat{\theta}_\lambda^S$ and $\hat{\theta}_\lambda$ have estimation errors of the same magnitude $n^{-2/d}$. Therefore, when tuning the parameter $\lambda$ depending on the number of observations and the dimension, we reach the rate $n^{-2/d}$. This same rate, up to logarithm factor, as for the estimator $\hat{\theta}_\lambda$ that we study throughout this article. We can also take into account the algorithm error for the estimator $\hat{\theta}_\lambda^S$ depending on the number of iterations $\ell$. All the results related to the estimator (3.8) introduced in this remark can be found in Section D of the Appendix.

When studying the sample complexity of the regularized OT cost, that is $\mathbb{E}[|W_\lambda(\mu, \nu) - W_\lambda(\hat{\mu}, \hat{\nu})|]$ as done in [21, 30], or when estimating the standard optimal transport cost $W_0(\mu, \nu)$ as in [14], bounds related to the control of the estimation error have been proved. These results give a control of $\mathbb{E}[|W_\lambda(\mu, \nu) - W_\lambda(\hat{\mu}, \hat{\nu})|]$ that is of order $C_\lambda/\sqrt{n}$ with $C_\lambda$ a constant that depends on the regularizing parameter. In the following section we give a similar result adapted to our context of weights estimation, and discuss why we favored the upper bound given in Proposition 3.1.

# 4 Alternative bound on the estimation error, and relation to state of the art

Proposition 3.1 gives a control of the estimation error that is independent of the regularization parameter $\lambda$. We now give a bound much closer to what is known in the literature, where a small regularization parameter severely impacts the rate of convergence.

**Proposition 4.1.** Let $\lambda > 0$, and suppose that all probability measures have compact supports included in $B(0, R)$. If for all components $\mu_k$ as well as for $\nu$ at least $n$ observations are available, the estimation error can be upper bounded as follows:

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_\lambda(\hat{\mu}_\theta, \hat{\nu}) - W_\lambda(\mu_\theta, \nu)| \right] \lesssim \frac{KM_\lambda}{\sqrt{n}}, \qquad (4.1) $$
With $M_λ := M_d \max \left( R^2, \frac{R^{\lfloor d/2 \rfloor + 1}}{\lambda^{\lfloor d/2 \rfloor}} \right)$.

The proof of this last proposition 4.1 can be found in Section B of the Appendix.

This last upper bound (4.1) seems appealing because independent of the dimension $d$ and going to zero at the same rate of $n^{-1/2}$. However, the constant $M_λ$ depends on the dimension $d$ and the regularization parameter $λ$. Thus, when one tries to exploit this bound (4.1), instead of $\mathcal{E}(n, d)$ like it is done in Theorem 3.1, one reaches the following upper bound on the expected risk of the estimator $\hat{\theta}_λ$:

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda_n}}, \nu) - W_0(\mu_{\theta^*}, \nu)] \lesssim \frac{1}{\lambda^{\lfloor d/2 \rfloor} \sqrt{n}} + \lambda \log \left( \frac{1}{\sqrt{d\lambda}} \right). \quad (4.2) $$

In the last inequality we have not taken into account the algorithm error. Balancing the two terms of the right-hand side of (4.2) leads to a regularization parameter $\lambda_n = n^{1/(2\lfloor d/2 \rfloor + 2)}$. Finally, under the assumptions of Theorem 3.1, using the estimation error (4.1) gives a slower rate of convergence than in Corollary 3.1. Indeed the expected excess risk of the estimator $\hat{\theta}_λ$ is upper bounded by

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda_n}}, \nu) - W_0(\mu_{\theta^*}, \nu)] \lesssim n^{-1/(2\lfloor d/2 \rfloor + 2)} \log(n), \quad \text{with} \quad \lambda_n = n^{1/(2\lfloor d/2 \rfloor + 2)}. $$

**Remark 4.1** (Estimation of $W_0(\mu, \nu)$). We can adapt the results established in Corollary 3.1 to estimate the optimal transport cost with a regularized transport cost. Indeed we can use a regularized plug-in estimators $W_λ(\hat{\mu}, \hat{\nu})$. This question is for instanced investigated in [14] where the estimation of $W_0(\mu, \nu)$ is based on the Sinkhorn divergence defined in (2.8). If $n$ samples are available from each measure $\mu$ and $\nu$, thanks to the estimator $S_λ(\hat{\mu}, \hat{\nu})$, they reach the rate of convergence (see [14, Prop. 4])

$$ \mathbb{E}[|S_\lambda(\hat{\mu}, \hat{\nu}) - W_0(\mu, \nu)|] \lesssim n^{-2/(2\lfloor d/2 \rfloor + 4)}, \quad (4.3) $$

for some well chosen regularization parameter $\lambda_n$ that depends on $n$. However, based on the results we established $\hat{\theta}_λ$ we can derive faster rates of convergence for $W_{\lambda_n}(\hat{\mu}, \hat{\nu})$ toward $W_0(\mu, \nu)$.

**Remark 4.2.** When setting $K = 1$ in Proposition 4.1, we recover a result close to what was previously established by other authors. To the best of our knowledge, the first result of this flavor was proven in [21, Thm. 3], but with a factor $e^{R^2/\lambda}$ in the multiplicative constant. Soon after, the result [30, Cor. 3] improved the sample complexity control by removing the exponential factor, and substituting a compact support assumption by a sub-Gaussian one. In Proposition 4.1, exploiting the compact supports assumption enables us to provide a multiplicative constant that scale with $R^{\lfloor d/2 \rfloor + 1} \lambda^{-\lfloor d/2 \rfloor}$. All these controls over the sample complexity dramatically deteriorate in high dimension when $\lambda$ is close to zero. To mitigate this high dimensional phenomenon, the recent article [37, Thm. 2] provides bounds where the impact of the ambient dimension $d$ is substituted by an intrinsic dimensional quantity of the measures compared.

**Proposition 4.2.** Suppose that $\mu$ and $\nu$ have their supports included in $B(0, R)$ and that the dimension $d > 4$. If $n$ samples are available for each probability measure, then the regularized plug-in estimator reaches the rate of convergence

$$ \mathbb{E}[|W_{\lambda_n}(\hat{\mu}, \hat{\nu}) - W_0(\mu, \nu)|] \lesssim n^{-2/d} \log(n) \quad \text{with} \quad \lambda_n = n^{-2/d}. \quad (4.4) $$

*Proof.* We have

$$ \mathbb{E}[|W_\lambda(\hat{\mu}, \hat{\nu}) - W_0(\mu, \nu)|] \le \mathbb{E}[|W_\lambda(\hat{\mu}, \hat{\nu}) - W_\lambda(\mu, \nu)|] + |W_\lambda(\mu, \nu) - W_0(\mu, \nu)|. $$

The first term on the right-hand side is upper bounded by $\mathcal{E}(n, d)$ thanks to Proposition 3.1 applied in the case $K = 1$ and $\mu_1 = \mu$. For the second term, the result established in [21, Thm. 1] gives a control of order $\lambda \log(1/\lambda)$ when $\lambda$ goes to zero. Hence, assuming that $d > 4$, and choosing $\lambda = n^{-2/d}$, we recover the rate of convergence claimed in equation (4.4). $\square$
Hence, in the case $d > 4$, the expected error of the estimator $W_{\lambda_n}(\hat{\mu}, \hat{\nu})$ goes to zero faster than when considering $S_{\lambda_n}(\hat{\mu}, \hat{\nu})$. While establishing inequality (4.4) only requires the measures $\mu$ and $\nu$ to have compact support, the previously known inequality (4.3) requires stronger assumptions on the measures $\mu$ and $\nu$.

**Remark 4.3.** To compare our results with the state-of-the-art, we set $K = 1$ in Proposition 4.2 to recover the sample complexity problem; that is controlling $\mathbb{E}[|W_{\lambda}(\hat{\mu}, \hat{\nu}) - W_{\lambda}(\mu, \nu)|]$. After the first version of our manuscript was posted on ArXiv, the preprint [25] established more general results than what we derive when taking $K = 1$ in Proposition 4.1. Indeed, the result [25, Thm. 3] provides a sample complexity control that holds for more general costs than the squared Euclidean distance we considered in the present article. Moreover, the authors of [25] provide some adaptive results where dimension-dependent bounds are substituted by intrinsic dimensional properties of the measures compared. To stress the difference between our work and [25], we also mention that in this article the authors study the Entropic Gromov-Wasserstein cost. Conversely, our main results, that are Theorem 3.1 and its Corollary 3.1, study an estimation framework where entropic optimal transport costs are exploited as losses function. Our estimation scenario is thus different from [25] where the quantities targeted are entropic optimal transport costs.

**Remark 4.4** (Near minimax-rate for the estimation of $W_0(\mu, \nu)$). It has been shown in [29, Thm. 21] that the minimax rate of convergence of $W_0(\mu, \nu)$ is lower bounded by $(n \log(n))^{-2/d}$ when $n$ observations from each measure are available. Up to a logarithmic factor, as shown in [14, Thm. 2], this rate is reached by the plug-in estimator $W_0(\hat{\mu}, \hat{\nu})$. An application of our work is to show that, up to another logarithmic factor, the regularized plug-in estimator $W_{\lambda_n}(\hat{\mu}, \hat{\nu})$ also reaches this rate of convergence. However, in some cases, the computation of $W_{\lambda_n}(\hat{\mu}, \hat{\nu})$ might be faster than $W_0(\hat{\mu}, \hat{\nu})$.

**Remark 4.5** (Computational cost of $W_{\lambda_n}(\hat{\mu}, \hat{\nu})$). One iteration of Sinkhorn algorithm requires $O(n^2)$ arithmetic operations [14, Page 5]. And we compute an approximation of $W_{\lambda_n}(\hat{\mu}, \hat{\nu})$ with $W_{\lambda_n}^{(\ell_n)}(\hat{\mu}, \hat{\nu})$ where $\ell_n = 32R^4 n^{4/d}$. Hence, the global cost of computing of $W_{\lambda_n}^{(\ell_n)}(\hat{\mu}, \hat{\nu})$ is of $O(n^{2+4/d})$ arithmetic operations. On the other hand computing $W_0(\hat{\mu}, \hat{\nu})$ with a linear programming algorithm requires $O(n^3)$ [32, 15]. Hence, as soon as $d > 4$, approximating $W_0(\mu, \nu)$ based on $n$ samples is faster with Sinkhorn algorithm than with a linear program. From our understanding, the computational advantages of entropic optimal transport are in high dimension. This is due to the fact that in high dimension the estimation error is large enough to allow for a choice of large $\lambda$, and thus a fast convergence of Sinkhorn algorithm.

# 5 Numerical experiments

In this section, using simulated and real data from flow cytometry, we analyze the numerical performances of the estimators introduced in Section 3. These numerical experiments have been designed to demonstrate how the parameters $\lambda$ and $\ell$ impact the performance of regularized Wasserstein estimators. Moreover, these experiments show that an appropriate choice of the parameters allows regularized estimators to reach the performance of estimators based on the standard OT cost $W_0$. For the results reported here, the parameter $\lambda$ ranges in a finite grid $\Lambda \subset \mathbb{R}_+^*$ from 0.01 to 1. Sinkhorn algorithm is either limited to $\ell = 5$ on simulated data, or to $\ell = 10$ on flow cytometry data. To simulate the setting where Sinkhorn in unlimited, our stopping criterion is based the difference between two consecutive outputs. More specifically, it stops if $|W_{\lambda}^{(\ell)} - W_{\lambda}^{(\ell-1)}| < 10^{-9}$.

For the estimators based on the transport cost $W_{\lambda}$ with $\lambda > 0$, we follow the protocol described thereafter. Given a data set $\mathbf{X}$ classified into $K$ classes, and an unclassified data set $\mathbf{Y}$ where we want to estimate the class proportions, we compute the empirical distributions $\hat{\mu}_1, \dots, \hat{\mu}_K$ and $\hat{\nu}$. Then, we compute the estimator $\hat{\theta}_{\lambda}$ of the class proportions by solving the optimization problem

$$ \hat{\theta}_{\lambda} = \arg \min_{\theta \in \Sigma_K} W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (5.1) $$
To solve this problem, we apply a gradient descent algorithm to the function $\theta \mapsto W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})$. In order to move from a constrained problem to an unconstrained one, we re-parameterize the simplex $\Sigma_K$ with a soft-max function $\chi : \mathbb{R}^K \to \Sigma_K$, where the $l$th component of $\chi(z)$ is defined by

$$ \chi(z)_l = \frac{\exp(z_l)}{\sum_{k=1}^{K} \exp(z_k)}. $$

Then, we introduce the linear operator $\Gamma : \Sigma_K \to \Sigma_n$ that maps the weights associated to each component $\hat{\mu}_k$ to the weights associated to each observations. That is

$$ \forall (i, k) \in \{1, \dots, n\} \times \{1, \dots, K\}, \quad \Gamma_{i,k} := \begin{cases} 1/n_k & \text{if } X_i \sim \mu_k, \\ 0 & \text{otherwise.} \end{cases} $$

Thus, our objective function reads $F_{\lambda} = W_{\lambda, \mathbf{X}}(\cdot, \hat{\nu}) \circ \Gamma \circ \chi$, where for $a \in \Sigma_n$, $W_{\lambda, \mathbf{X}}(a, \hat{\nu})$ denotes the transport cost between the measure with weights $a$ and support $\mathbf{X} = (X_1, \dots, X_n)$, and the measure $\hat{\nu}$. From [34, Prop. 9.1], we know that the gradient of $W_{\lambda, \mathbf{X}}(\cdot, \hat{\nu})$ at point $a$ is given by the unique dual potential $\varphi$ associated to the measure $\sum_{i=1}^n a_i \delta_{X_i}$ such that $\sum_{i=1}^n \varphi_i = 0$. From the chain rule of differentiation, we derive that the gradient of the objective function is given by

$$ \nabla_z F_{\lambda}(z) = J_{\chi}(z)^T \Gamma^T \varphi_z, \qquad (5.2) $$

where $J_{\chi}$ is the Jacobian matrix of $\chi$ and $\varphi_z$ is the optimal potential with respect to $W_{\lambda}(\hat{\mu}_{\chi(z)}, \hat{\nu})$. Our approximation of $\varphi_z$ computed is with the Sinkhorn algorithm.

**Algorithm 1:** Approximation of $\hat{\theta}_{\lambda}^{(\ell)} = \arg \min_{\theta \in \Sigma_K} W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})$.

```c
z ← 1_K
for N ← 1 to N_out do
    /* Sinkhorn Algorithm to compute the dual potentials of W_λ(μ_χ(z), ν) */
    ϕ ← 0_n
    for l ← 1 to ℓ do
        /* One Sinkhorn algorithm iteration */
        ψ ← ϕ^{c,λ}
        ϕ ← ψ^{c,λ}
    end
    /* Approximation of the gradient of z ↦→ W_λ(μ_{θ_z}, ν) where θ_z = χ(z) */
    ω(z) ← (ΓJ_χ(z))^T ϕ
    z ← z − ηω(z)
end
/* Approximation of the estimator θ_λ */
return θ_λ = χ(z)
```

When relying on the other transport costs studied, that are $W_0, S_{\lambda}, W_{\lambda}^{(\ell)}$ or $S_{\lambda}^{(\ell)}$, we apply the same protocol as for $W_{\lambda}$; apart from the gradient formula (5.2). Indeed, denoting by $\mathcal{L}$ an optimal transport cost among $W_0, S_{\lambda}, W_{\lambda}^{(\ell)}$ and $S_{\lambda}^{(\ell)}$, the problem we are trying to solve (after parameterization with the soft-max function $\chi$) is

$$ \min_{z \in \mathbb{R}^K} \mathcal{L}(\hat{\mu}_{\chi(z)}, \hat{\nu}). \qquad (5.3) $$

We can rewrite the objective function $F: z \to F(z) = \mathcal{L}_{\mathbf{X}}(\cdot, \hat{\nu}) \circ \Gamma \circ \chi(z)$. Here, $\mathcal{L}_{\mathbf{X}}(a, \hat{\nu})$ denotes the transport cost criterion between the measure with weights $a \in \Sigma_n$ and support $\mathbf{X} = (X_1, \dots, X_n)$, and the measure $\hat{\nu}$. Then, differentiating this function $F$ gives the gradient

$$ \nabla_z F(z) = J_{\chi}(z)^T \Gamma^T \nabla \mathcal{L}_{\mathbf{X}}(\Gamma \chi(z), \hat{\nu}). $$
Finally, depending on the loss $\mathcal{L}$, we substitute $\nabla \mathcal{L}_{\mathbf{X}}(\Gamma\chi(z), \hat{\nu})$ by its value. For the unregularized case, we have $\nabla W_{0,\mathbf{X}}(\Gamma\chi(z), \hat{\nu}) = \varphi_z$ with $\varphi_z$ the first potential associate to $W_0(\hat{\mu}_{\chi(z)}, \hat{\nu})$. We rely on a linear programming algorithm to approximate this dual vector $\varphi_z$, which in this case is a sub-gradient [34, Prop. 9.1]. For $W_{\lambda}^{(\ell)}$, as in Algorithm 1, we rely on Sinkhorn algorithm to compute $\varphi^{(\ell)}$, the dual potential after $\ell$ iterations. For $S_{\lambda}$, the gradient is given by the formula $\nabla S_{\lambda,\mathbf{X}}(\Gamma\chi(z), \hat{\nu}) = \varphi^{\mu,\nu} - (\varphi^{\mu,\mu} + \psi^{\mu,\mu})/2$ [9, eq. 2.12], where $\varphi^{\mu,\nu}$ is the first potential associated to $W_{\lambda}(\hat{\mu}_{\chi(z)}, \hat{\nu})$, and $\varphi^{\mu,\mu}, \psi^{\mu,\mu}$ are the two potentials associated to $W_{\lambda}(\hat{\mu}_{\chi(z)}, \hat{\mu}_{\chi(z)})$. For $S_{\lambda}^{(\ell)}$, its gradient is given by the same formula as $S_{\lambda}$ while substituting the potentials by their approximations after $\ell$ steps of the Sinkhorn algorithm.

**Remark 5.1.** The algorithm described in the present article is fairly similar to the numerical scheme exploited in [19]. However, in the present work, the approximation of the dual potential required to compute the gradient of $W_{\lambda}$ is based on Sinkhorn algorithm. While in the previous work [19], the authors applied the stochastic optimization algorithm studied in [22, 5]. Relying on Sinkhorn algorithm enables us to incorporate the algorithmic error in our theoretical study.

For each setting, that is choosing a loss among $W_{\lambda}$, $S_{\lambda}$, $W_{\lambda}^{(\ell)}$ or $S_{\lambda}^{(\ell)}$; and setting the parameters $\lambda$ and $\ell$, we sample (or sub-sample when experimenting on real data) $N = 50$ couples of datasets $(\mathbf{X}^{[1]}, \mathbf{Y}^{[1]}), \dots, (\mathbf{X}^{[N]}, \mathbf{Y}^{[N]})$. Then, for each couple $(\mathbf{X}^{[r]}, \mathbf{Y}^{[r]})$, we compute an estimator $\hat{\theta}^{[r]}$ of the class proportions in $\mathbf{Y}^{[r]}$. We thus obtain $N$ realizations $\hat{\theta}^{[1]}, \dots, \hat{\theta}^{[N]}$ of a given estimator of the class proportions. Then, we choose to evaluate performance of the estimator, by computing the quadratic errors $\|\hat{\theta}^{[1]} - \theta^{*}\|^2, \dots, \|\hat{\theta}^{[N]} - \theta^{*}\|^2$. We display these error with box plots as in Figure 2, where circles are the errors $\|\hat{\theta}^{[r]} - \theta^{*}\|^2$ beyond 1.5 times the inter-quartile range. When experimenting on synthetic data, $\theta^{*}$ is known as $\nu = \sum_{k=1}^{K} \theta_{k}^{*} \mu_{k}$ as ensured by Lemma 2.1. In experiments on cytometry data, $\theta^{*}$ is unknown because all probability measures underlying the observations are unknown. In this case, we substitute $\theta^{*}$ by the true proportions in the unclassified data set $\mathbf{Y}$, to which we actually have access.

We also approximate the expected quadratic risk $\mathbb{E}[\|\hat{\theta} - \theta^{*}\|^2]$ by Monte-Carlo repetitions as classically done in statistical experiments:

$$ \mathbb{E}[\|\hat{\theta} - \theta^{*}\|^2] \approx \frac{1}{N} \sum_{r=1}^{N} \|\hat{\theta}^{[r]} - \theta^{*}\|^2. \qquad (5.4) $$

We plot this approximated average error on Figure 3, when considering for instance the losses $W_{\lambda}$ and $S_{\lambda}$. This protocol is repeated for each value of $\lambda$ in the grid $\Lambda$ and each loss function.

**Remark 5.2.** In these numerical experiments, we have chosen to focus on the expected error $\mathbb{E}[\|\hat{\theta} - \theta^{*}\|^2]$ rather than the expected excess risk $r_n(\mu_{\hat{\theta}}, \nu)$ as in flow cytometry the relevant quantity is an accurate estimation of class proportions in the target dataset. Also, notice that the risk $r_n(\mu_{\hat{\theta}}, \nu)$ cannot be computed exactly because it involves the quantity $W_0(\mu_{\hat{\theta}}, \nu)$ for which we have no closed-form formula.

## 5.1 Simulated data

We first simulated two Gaussian mixtures of dimension $d = 6$ with the same $K = 5$ components but with different class proportions. Thus, a source data set $\mathbf{X}$ corresponds to random vectors $X_1, \dots, X_n$ sampled with respect to $\mu$ and a target data set $\mathbf{Y}$ corresponds to random vectors $Y_1, \dots, Y_n$ sampled with respect to the distribution $\nu$, where $\mu$ and $\nu$ are defined below:

$$ \mu = \sum_{k=1}^{5} \pi_k \mathcal{N}(\rho_k, \sigma^2 I_d), \qquad \nu = \sum_{k=1}^{5} \theta_k^* \mathcal{N}(\rho_k, \sigma^2 I_d). \qquad (5.5) $$
Because the vector of proportions $\pi$ and $\theta^*$ are not assumed to be equal, we exploit the known
classes at the source in order to estimate the class proportions $\theta^*$ at the target, based on empirical
versions of $\mu_1, \dots, \mu_K$ and $\nu$.

We have same number of samples $m_k = n$ from each source components $\mu_k$ than samples from
the target distribution $\nu$. This experimentation setting matches the presentation of our theoretical
results given in Section 3.2. To ease the simulation study, we constrain the number of observations
to $m_k = 50$ observations for each class of the source data set. In the target data set, we also
constrain the number of observations per class with $n_1 = 20, n_2 = 5, n_3 = 8, n_4 = 7, n_5 = 10$, so
$n = 50$ in total. We display in Figure 1 two-dimensional projections of one dataset from the source
measure and one dataset from the target measure with their respective clustering. Note that the
clustering of the target dataset is then assumed to be unknown.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** 2D projections of a simulated source data set and a target data set with their clustering.</figcaption></figure>

5.1.1 Unlimited number of Sinkhorn iterations

Through a first series of experiments, we compare the performances of the estimators computed
with the losses $W_0$, $W_\lambda$ and $S_\lambda$. In Figure 2, using a boxplot we display the behavior of the error
$\|\hat{\theta}_\lambda - \theta^*\|^2$ for each value of the regularization parameter $\lambda \in \Lambda$. In Figure 3, we also display
the estimation of $\mathbb{E}[\|\hat{\theta}_\lambda - \theta\|^2]$ using the Monte-Carlo estimator (5.4). For small values of $\lambda$, the
regularized losses $W_\lambda$ and $S_\lambda$ yield competitive estimators compared to the one obtained with $W_0$.
Notice that the regularization parameter advised from Corollary 3.1 is $\lambda_n := n^{-2/d}$. In this first
series of experiments $n = 50$ and $d = 6$, that gives $\lambda_n \approx 0.27$. This parameter $\lambda$ is slightly larger,
than suggested by our empirical results from Figure 2. This gap between theory and practice
might be explained by the fact that we did not take into account the multiplicative constant in
the approximation error. According to [21], this constant is of order $2d\lambda \log(1/\lambda)$. Thus, taking
this constant into account would give a regularization parameter $\tilde{\lambda}_n = (2d)^{-1}\lambda \log(1/\lambda)$, which is
closer to the parameters that perform best in these experiments.

We also point out that the computational complexity of the Sinkhorn algorithm is highly
dependent on the regularization parameter $\lambda$ as discussed in [16], [2]. To illustrate this fact, we
display in Figure 4 the time (in seconds) required to compute $N = 50$ samples of $\hat{\theta}_\lambda$ depending
on the value of $\lambda$. As $\nabla_\theta S_\lambda(\mu_\theta, \nu) = \nabla_\theta W_\lambda(\mu_\theta, \nu) - \frac{1}{2} \nabla_\theta W_\lambda(\mu_\theta, \mu_\theta)$, computing the gradient of
$S_\lambda(\mu_\theta, \nu)$, requires to solve the dual problem associated to $W_\lambda(\mu_\theta, \mu_\theta)$ in addition to the dual
problem associated to $W_\lambda(\mu_\theta, \nu)$. But as noticed in [17], Sinkhorn algorithm converges much
faster for the symmetric term $W_\lambda(a, a)$ than in the general case when computing $W_\lambda(a, b)$. We
have observed in our experiment that the number of iterations before reaching convergence when
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Estimation results on simulated data without limitation on the number of iterations of the Sinkhorn algorithm. We display the error $||\hat{\theta}_{\lambda} - \theta^{*}||^2$ using either the loss $W_{\lambda}$ (left) or $S_{\lambda}$ (right). The black line is the median error of the un-regularized estimator $\hat{\theta}_{0}$ using the standard optimal transport cost $W_{0}$, while the dotted lines are the first and third quartiles of the errors of estimation $||\hat{\theta}_{0} - \theta^{*}||^2$. Circles are the errors $||\hat{\theta}_{\lambda} - \theta^{*}||^2$ beyond 1.5 times the inter-quartile range.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Average error on simulated data of the estimators $\hat{\theta}_{\lambda}$ (orange) and $\hat{\theta}_{\lambda}^{S}$ (blue) as a function of the regularization parameter $\lambda$. Sinkhorn algorithm runs until convergence is reached. That is, when the difference between two iterates $|W^{(\ell)} - W^{(\ell-1)}|$ is less than $10^{-9}$. The black dotted line is the average error of the un-regularized estimator $\hat{\theta}_{0}$.</figcaption></figure>
computing $W_\lambda(a, a)$ does not seem to be a monotonic function with respect to the regularization parameter $\lambda$. This partially accounts for the slightly longer time of computation for $\lambda = 0.02$ in comparison to $\lambda = 0.01$ on the right side of Figure 4, that is when using $S_\lambda$ as loss function.

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Time required to compute $N = 50$ estimators $\hat{\theta}_{\lambda}$ (left) and $\hat{\theta}_{\lambda}^{S}$ (right) when the number of iterations is unlimited.</figcaption></figure>

### 5.1.2 Limited number of Sinkhorn iterations

Figure 2 and Figure 4 presents results questioning the trade-off between the computational cost of regularized OT and the quality of statistical estimation. We have repeated the experiments of Section 5.1.1 by now constraining the number iterations of the Sinkhorn algorithm to be equal to $\ell = 5$ for any value $\lambda$. In other words, we compute the estimators $\hat{\theta}_{\lambda}^{(\ell)}$ and $\hat{\theta}_{\lambda}^{S(\ell)}$ with $\ell = 5$, thus fixing the computational budget. Figure 5 and Figure 6 both present the performances of those estimators: by limiting the number of Sinkhorn iterations, the accuracy of the estimation deteriorates for small values of $\lambda$. This degradation comes from $\ell = 5$ being too small a number of iterations for the Sinkhorn algorithm to converge for small values of $\lambda$. Yet Figure 5 points to some values of $\Lambda$ as offering a nice trade-off between the computational cost of small $\lambda$ and the approximation error induced by larger $\lambda$. For such values, the performances of the regularized estimators $\hat{\theta}_{\lambda}^{(\ell)}$ and $\hat{\theta}_{\lambda}^{S(\ell)}$ are seen to be comparable to those of the un-regularized estimator $\hat{\theta}_0$. However, we must grant a minor divergence between our theoretical findings of Corollary 3.1. Indeed, our theoretical results suggest that $\ell$ should be set of order $\ell_n = n^{4/d}$, which gives $\ell_n \approx 14$ in this context. We suspect two reasons for this gap. Firstly, some constants in the estimation errors $\hat{\theta}$ and $\hat{\theta}^S$ are unknown. In such a case, allowing a larger algorithm error by choosing $\ell$ smaller would not reduce the performance of estimator considered. A second source of error in these experiments is that we are not exactly under the assumptions of Corollary 3.1. Indeed, this corollary requires measures to have compact support. While in our experiments, the probability measures are Gaussian variables, which do not have compact supports.

## 5.2 Flow cytometry data

We now apply our method of class proportions estimation on flow cytometry data. We demonstrate that the regularization parameter $\lambda$ has also a significant impact on the estimation of class proportions on real data. As an illustrative example, we apply our technique to flow cytometry data sets from the T-cell panel of the Human Immunology Project Consortium (HIPC) – publicly available on ImmuneSpace [10]. We arbitrarily chose two data sets that comes from cytometry measurements performed in the “Stanford” laboratory center. One data set, that acts as the source
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Estimation results on simulated data when the number of iterations of the Sinkhorn algorithm is limited to $\ell = 5$. We display the error $\|\hat{\theta}_{\lambda}^{(\ell)} - \theta^*\|^2$ using either the loss $W_{\lambda}^{(\ell)}$ (left) or $S_{\lambda}^{(\ell)}$ (right). The black line is the median error of the un-regularized estimator $\hat{\theta}_0$ using the standard optimal transport cost $W_0$, while the dotted lines are the first and third quartiles of $\|\hat{\theta}_0 - \theta^*\|^2$.</figcaption></figure>

measure, is built from observations measured from a biological sample of a certain patient. Another second data set, acting as the target measure, is built from the observations obtained from a biological sample that comes from another patient. After performing cytometry measurements the observations were manually gated into 10 cell populations: CD4 Effector (CD4 E), CD4 Naive (CD4 N), CD4 Central memory (CD4 CM), CD4 Effector memory (CD4 EM), CD4 Activated (CD4 A), CD8 Effector (CD8 E), CD8 Naive (CD8 N), CD8 Central memory (CD8 CM), CD8 Effector memory (CD8 EM) and CD8 Activated (CD8 A). Hence, for these data sets, a manual clustering is at our disposal to evaluate the performances of our method. In this context $\theta^*$ is defined as the class proportions defined thanks to the manual gating. For each cell, seven biological markers have been measured, and it thus leads to observations $X_i$ and $Y_j$ that belong to $\mathbb{R}^d$ with $d = 7$. A two-dimensional projection of these datasets is displayed in Figure 7 with the resulting manual clustering.

### 5.2.1 Unlimited Sinkhorn iterations

We reproduce the protocol that we have considered in the case of simulated data. To build an empirical distribution of the source distribution when analyzing flow cytometry data, we sub-sample 50 observations from each class of the source data set in order to construct the empirical measures $\hat{\mu}_1, \dots, \hat{\mu}_K$, and to define $\hat{\mu}_\theta = \sum_{k=1}^K \theta_k \hat{\mu}_k$ for $\theta \in \Sigma_K$. Figure 8 shows two sub-samples from the source and target distributions displayed in Figure 7.

We recall that the clustering of the target dataset is not used in the estimation procedure.

The numerical performances of the estimators computed with the three loss functions $W_0$, $W_\lambda$ and $S_\lambda$ are displayed on Figure 9 and Figure 10. In the context of flow cytometry data, the underlying distributions $\mu$ and $\nu$ are obviously unknown, and the quantity $\min W_0(\mu_\theta, \nu)$ is thus not accessible. Therefore, we define the optimal vector $\theta^*$ of class proportions to be the one in the fully observed (not sub-sampled) target dataset that is displayed in Figure 7. Those results on real data are consistent with the results of simulated data. Indeed, one can observe that for small values of $\lambda \in \Lambda$ the accuracy of the estimation obtained with the loss functions $W_\lambda$ and $S_\lambda$ is very similar to the one obtained using $W_0$.
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Average error on simulated data of the estimators $\hat{\theta}_{\lambda}^{(\ell)}$ (orange) and $\hat{\theta}_{\lambda}^{S(\ell)}$ (blue) as a function of the regularization parameter $\lambda$ with a limitation on the number of iterations. For all values of $\lambda$, Sinkhorn algorithm is limited to $\ell = 5$ iterations. The black dotted line is the average error of the un-regularized estimator $\hat{\theta}_0$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Two-dimensional projection of the flow cytometry datasets used in these numerical experiments with a clustering of the cells into 10 sub-populations. Note that the true dimension of the data is $d = 7$, thus limiting the readability of such 2D projections.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8:** Sub-sample of the source and target flow cytometry data sets. In the source sub-sample, $m_k = m = 50$ elements of each class have been sampled. In the target sub-sample, $n = 100$ observations have been randomly chosen.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Results on HIPC data without imposing limitations on the number of Sinkhorn iterations. We display the error $||\hat{\theta}_{\lambda} - \theta^{*}||^2$ using either the loss $W_{\lambda}$ (left) or $S_{\lambda}$ (right). The black line is the median error of the un-regularized estimator $\hat{\theta}_{0}$ using the loss $W_{0}$, while the dotted lines are the first and third quartiles of $||\hat{\theta}_{0} - \theta^{*}||^2$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10:** Average error on HIPC data of the estimators $\hat{\theta}_{\lambda}$ (orange) and $\hat{\theta}_{\lambda}^{S}$ (blue) as a function of the regularization parameter $\lambda$. There is no limitation on the number of iterations, Sinkhorn algorithm runs until convergence is reached. The black dotted line is the average error of the un-regularized estimator $\hat{\theta}_{0}$.</figcaption></figure>

### 5.2.2 Limited Sinkhorn iterations

In order to reduce the computational cost of our estimation method, we limit the number of Sinkhorn iterations to $\ell = 10$. Once again, the results displayed in Figure 11 and Figure 12 show that it is possible to propose a competitive alternative to $W_0$ at a lower computational cost.

<figure><img src="image_13.png" /><figcaption>**Figure 11:** Results on HIPC data when the number of Sinkhorn iterations is limited to $\ell = 10$. We display boxplots of the error $\|\hat{\theta}_{\lambda} - \theta^{*}\|^{2}$ using either the loss $W_{\lambda}$ (left) or $S_{\lambda}$ (right). The black line is the median error of the un-regularized estimator $\hat{\theta}_{0}$ using the loss $W_{0}$, while the dotted lines are the first and third quartiles of $\|\hat{\theta}_{0} - \theta^{*}\|^{2}$.</figcaption></figure>

## 6 Conclusion and discussion

In this work, we have presented a thorough study of Wasserstein estimators based on regularized OT with an emphasis on the influence of the regularization parameter $\lambda$. This study was carried out through the example of a mixture model and weights estimation. We derived upper bounds on the risk of Wasserstein estimators in terms of an estimation error and an approximation error.
<figure><img src="image_14.png" /><figcaption>**Figure 12:** Average error on HIPC data of the estimators $\hat{\theta}_{\lambda}^{(\ell)}$ (orange) and $\hat{\theta}_{\lambda}^{S(\ell)}$ (blue) as a function of the regularization parameter $\lambda$ with a limitation on the number of iterations. For all values of $\lambda$, Sinkhorn algorithm is limited to $\ell = 10$ iterations. The black dotted line is the average error of the un-regularized estimator $\hat{\theta}_0$</figcaption></figure>

We assessed the influence of the chosen OT-based loss (among $W_{\lambda}$, $S_{\lambda}$ and $W_0$) on the decay of the estimation and approximation terms. We have also proposed an optimal decay of the regularization parameter $\lambda = \lambda_n$ based on these upper bounds to achieve decreasing rate of $n^{-2/d}$ for the expected excess risk. Secondly, motivated by the sensitive question of the computational cost of regularized OT, we have studied the algorithmic error induced by limiting the number of iterations in the Sinkhorn algorithm. This study resulted in a principled strategy to set the number of Sinkhorn iterations $\ell = \ell_n$ in order to maintain the algorithm error below the statistical error. We have also demonstrated with numerical experiments that an appropriate choice of $\lambda$ and a limited number of Sinkhorn iterations $\ell$ allow to equal the performances of the un-regularized estimator at a reduced computational cost.

We now present a few perspectives for future research. For an estimator $\hat{\theta}_n$ of $\theta^*$, we have derived a control on the excess risk, that is $W_0(\mu_{\hat{\theta}_n}, \nu) - W_0(\mu_{\theta^*}, \nu)$. However, a direct control of the weights estimator, i.e. of the quantity $\|\hat{\theta}_n - \theta^*\|$, would be even more valuable. For instance, a control of $\|\hat{\theta}_n - \theta^*\|$ may allow to develop statistical tests on the estimator $\hat{\theta}_n$. An other possible extension of this work is suggested by our numerical experiments. Figure 2 and Figure 5 indicate that limiting the number of iterations for Sinkhorn algorithm could improve statistical performance. These better results with limited iterations are not accounted by the present work. Hence, further investigation on this observation is an other direction for research.

## References

[1] N. Aghaeepour, G. Finak, H. Hoos, T. Mosmann, R. Brinkman, R. Gottardo, R. Scheuermann, F. Consortium, D. Consortium, et al. Critical assessment of automated flow cytometry data analysis techniques. *Nature methods*, 10(3):228, 2013.

[2] J. Altschuler, J. Niles-Weed, and P. Rigollet. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration. *Advances in neural information processing systems*, 30, 2017.

[3] M. Arjovsky, S. Chintala, and L. Bottou. Wasserstein generative adversarial networks. In *International Conference on Machine Learning*, pages 214–223, 2017.
[4] M. Ballu, Q. Berthet, and F. Bach. Stochastic optimization for regularized wasserstein estimators. In *International Conference on Machine Learning*, pages 602–612. PMLR, 2020.

[5] B. Bercu and J. Bigot. Asymptotic distribution and convergence rates of stochastic algorithms for entropic optimal transportation between probability measures. *The Annals of Statistics*, 49(2):968–987, 2021.

[6] E. Bernton, P. E. Jacob, M. Gerber, and C. Robert. On parameter estimation with the Wasserstein distance. *Information and Inference: A Journal of the IMA*, 8(4):657–676, 2019.

[7] J. Bigot. Statistical data analysis in the Wasserstein space. *ESAIM: ProcS*, 68:1–19, 2020.

[8] J. Bigot, E. Cazelles, and N. Papadakis. Data-driven regularization of Wasserstein barycenters with an application to multivariate density registration. *Information and Inference: A Journal of the IMA*, 8, 04 2018.

[9] J. Bigot, E. Cazelles, and N. Papadakis. Central limit theorems for entropy-regularized optimal transport on finite spaces and statistical applications. *Electronic Journal of Statistics*, 13(2):5120–5150, 2019.

[10] V. Brusic, R. Gottardo, S. Kleinstein, M. Davis, D. Hafler, H. Quill, A. Palucka, G. Poland, B. Pulendran, E. Reinherz, et al. Computational resources for high-dimensional immune analysis from the human immunology project consortium. *Nature biotechnology*, 32(2):146, 2014.

[11] M. Cheung, J. J. Campbell, L. Whitby, R. J. Thomas, J. Braybrook, and J. Petzing. Current trends in flow cytometry automated data analysis software. *Cytometry Part A*, 99(10):1007–1021, 2021.

[12] S. Chewi and A.-A. Pooladian. An entropic generalization of Caffarelli's contraction theorem via covariance inequalities. *arXiv preprint arXiv:2203.04954*, 2022.

[13] L. Chizat. Doubly Regularized Entropic Wasserstein Barycenters, 2023.

[14] L. Chizat, P. Roussillon, F. Léger, F. Vialard, and G. Peyré. Faster Wasserstein Distance Estimation with the Sinkhorn Divergence. In *Proc. NeurIPS'20*, 2020.

[15] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In *Advances in neural information processing systems*, pages 2292–2300, 2013.

[16] P. Dvurechensky, A. Gasnikov, and A. Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent is better than by sinkhorn algorithm. In *International conference on machine learning*, pages 1367–1376. PMLR, 2018.

[17] J. Feydy, T. Séjourné, F. Vialard, S. Amari, A. Trouve, and G. Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In *The 22nd International Conference on Artificial Intelligence and Statistics*, pages 2681–2690, 2019.

[18] N. Fournier and A. Guillin. On the rate of convergence in Wasserstein distance of the empirical measure. *Probability Theory and Related Fields*, 162(3):707–738, 2015.

[19] P. Freulon, J. Bigot, and B. P. Hejblum. Cytopt: Optimal transport with domain adaptation for interpreting flow cytometry data. *The Annals of Applied Statistics*, 17(2):1086–1104, 2023.

[20] G. Biau and M. Sangnier and U. Tanielian. Some Theoretical Insights into Wasserstein GANs. *Journal of Machine Learning Research*, 22(119):1–45, 2021.

[21] A. Genevay, L. Chizat, F. Bach, M. Cuturi, and G. Peyré. Sample complexity of sinkhorn divergences. In *The 22nd International Conference on Artificial Intelligence and Statistics*, pages 1574–1583. PMLR, 2019.
[22] A. Genevay, M. Cuturi, G. Peyré, and F. Bach. Stochastic optimization for large-scale optimal transport. In *Advances in neural information processing systems*, pages 3440–3448, 2016.

[23] A. Genevay, G. Peyré, and M. Cuturi. Learning generative models with sinkhorn divergences. In A. Storkey and F. Perez-Cruz, editors, *Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics*, volume 84 of *Proceedings of Machine Learning Research*, pages 1608–1617. PMLR, 2018.

[24] A. González-Sanz and S. Hundrieser. Weak limits for empirical entropic optimal transport: Beyond smooth costs. *arXiv preprint arXiv:2305.09745*, 2023.

[25] M. Groppe and S. Hundrieser. Lower complexity adaptation for empirical entropic optimal transport. *arXiv preprint arXiv:2306.13580*, 2023.

[26] M. Hardy. Combinatorics of Partial Derivatives. *The Electronic Journal of Combinatorics*, pages R1–R1, 2006.

[27] G. Huizing, G. Peyré, and L. Cantini. Optimal transport improves cell-cell similarity inference in single-cell omics data. *bioRxiv*, 2021.

[28] D. Liu, M. T. Vu, S. Chatterjee, and L. K. Rasmussen. Entropy-regularized optimal transport generative models. In *ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 3532–3536. IEEE, 2019.

[29] T. Manole and J. Niles-Weed. Sharp convergence rates for empirical optimal transport with smooth costs. *arXiv preprint arXiv:2106.13181*, 2021.

[30] G. Mena and J. Niles-Weed. Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem. *Advances in Neural Information Processing Systems*, 32, 2019.

[31] V. M. Panaretos and Y. Zemel. Statistical Aspects of Wasserstein Distances. *Annual Reviews of Statistics and its Applications*, 6:405–431, 2018.

[32] O. Pele and M. Werman. Fast and robust earth mover’s distances. In *2009 IEEE 12th international conference on computer vision*, pages 460–467. IEEE, 2009.

[33] A. Petersen, C. Zhang, and P. Kokoszka. Modeling Probability Density Functions as Data Objects. *Econometrics and Statistics*, 21(C):159–178, 2022.

[34] G. Peyré and M. Cuturi. Computational optimal transport. *Foundations and Trends® in Machine Learning*, 11(5-6):355–607, 2019.

[35] M. Sanjabi, J. Ba, M. Razaviyayn, and J. D. Lee. On the convergence and robustness of training gans with regularized optimal transport. In *Advances in Neural Information Processing Systems*, pages 7091–7101, 2018.

[36] F. Santambrogio. Optimal transport for applied mathematicians. *Birkhäuser, NY*, 55(58-63):94, 2015.

[37] A. J. Stromme. Minimum intrinsic dimension scaling for entropic optimal transport. In *International Conference on Soft Methods in Probability and Statistics*, pages 491–499. Springer, 2024.

[38] A. van der Vaart and J. Wellner. *Weak Convergence and Empirical Processes. With Applications to Statistics*. New York: Springer, 1996.

[39] R. van Handel. *Probability in High Dimension*. Princeton University, 2016.
# A Proofs of the main results

The goal of this section is to derive the rate of convergence of regularized Wasserstein estimators. That is, when considering $W_\lambda$ as a loss function with $\lambda \ge 0$. To be more specific, we investigate in this section the estimators

$$ \hat{\theta}_{\lambda} := \arg \min_{\theta \in \Sigma_K} W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) \quad \text{for } \lambda \ge 0, \qquad (A.1) $$

or

$$ \hat{\theta}_{\lambda}^{(\ell)} \in \arg \min_{\theta \in \Sigma_K} W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}) \quad \text{for } \lambda > 0, \qquad (A.2) $$

when taking into account the algorithm error.

## A.1 Decomposition of the excess risk

We first detail how the excess risk of $\hat{\theta}_{\lambda}^{(\ell)}$ defined in equation (A.2) can be upper bounded by the sum of three terms. They represent a tradeoff between an approximation error, estimation error and an algorithm error.

**Lemma A.1.** Set $\lambda \ge 0$. the excess risk of the estimator $\hat{\theta}_{\lambda}^{(\ell)}$ defined by (A.2) is bounded as follows:

$$ \begin{aligned} 0 \le W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) \le & 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - W_\lambda(\mu_{\theta}, \nu)| \\ & + 2 \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_{\theta}, \nu) - W_\lambda(\hat{\mu}_{\theta}, \hat{\nu})| \\ & + 2 \sup_{\theta \in \Sigma_K} |W_\lambda(\hat{\mu}_{\theta}, \hat{\nu}) - W_\lambda^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|. \end{aligned} \qquad (A.3) $$

*Proof.* We begin with the decomposition

$$ \begin{aligned} W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) = & W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_\lambda(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) + W_\lambda(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_\lambda(\hat{\mu}_{\theta^*}, \hat{\nu}) \\ & + W_\lambda(\hat{\mu}_{\theta^*}, \hat{\nu}) - W_\lambda^{(\ell)}(\hat{\mu}_{\theta^*}, \hat{\nu}) + W_\lambda^{(\ell)}(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_\lambda^{(\ell)}(\mu_{\theta^*}, \nu) \\ & + W_\lambda^{(\ell)}(\mu_{\theta^*}, \hat{\nu}) - W_\lambda(\mu_{\theta^*}, \hat{\nu}) + W_\lambda(\mu_{\theta^*}, \hat{\nu}) - W_\lambda(\mu_{\theta^*}, \nu) \\ & + W_\lambda(\mu_{\theta^*}, \nu) - W_0(\mu_{\theta^*}, \nu). \end{aligned} \qquad (A.4) $$

Let us focus on the right hand side of this last equation (A.4). The first and the last differences are controlled by the approximation error $\sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \nu) - W_0(\mu_\theta, \nu)|$. The second and sixth differences can be upper bounded by the estimation error $\sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \nu) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})|$. The third and fifth differences are upper bounded by the algorithm error $\sup_{\theta \in \Sigma_K} |W_\lambda^{(\ell)}(\hat{\mu}_\theta, \hat{\nu}) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})|$.

It only remains to control $W_\lambda^{(\ell)}(\hat{\mu}_{\hat{\theta}_{\lambda}^{(\ell)}}, \hat{\nu}) - W_\lambda^{(\ell)}(\mu_{\theta^*}, \hat{\nu})$. However, $\hat{\theta}_{\lambda}^{(\ell)}$ minimizes the function $\theta \mapsto W_\lambda^{(\ell)}(\hat{\mu}_{\theta^*}, \hat{\nu})$. Hence $W_\lambda^{(\ell)}(\hat{\mu}_{\hat{\theta}_{\lambda}^{(\ell)}}, \hat{\nu}) - W_\lambda^{(\ell)}(\mu_{\theta^*}, \hat{\nu}) \le 0$.

Going back to equation (A.4) and substituting every difference of the right hand side by its appropriate bound we derive

$$ \begin{aligned} W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) \le & 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - W_\lambda(\mu_{\theta}, \nu)| \\ & + 2 \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_{\theta}, \nu) - W_\lambda(\hat{\mu}_{\theta}, \hat{\nu})| \\ & + 2 \sup_{\theta \in \Sigma_K} |W_\lambda(\hat{\mu}_{\theta}, \hat{\nu}) - W_\lambda^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|, \end{aligned} $$

which is the result claimed in Lemma A.1.

□
## A.2 Control of the estimation error

To control the estimation error, we split it into two terms:

$$ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \nu) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})| \le \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \hat{\nu}) - W_\lambda(\mu_\theta, \nu)| \\ + \sup_{\theta \in \Sigma_K} |W_\lambda(\hat{\mu}_\theta, \hat{\nu}) - W_\lambda(\mu_\theta, \hat{\nu})|. \quad (A.5) $$

Hence, controlling the expected estimation error boils down to controlling the (closely related) quantities

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_\lambda(\mu_\theta, \nu) - W_\lambda(\hat{\mu}_\theta, \hat{\nu}) \right| \right], \quad \text{and} \quad \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_\lambda(\hat{\mu}_\theta, \hat{\nu}) - W_\lambda(\mu_\theta, \hat{\nu}) \right| \right]. $$

The upper bounds we will derive are based on the supremum of an empirical process that has been previously studied by Chizat et al. in [14]. More precisely, we rely on the following Lemma.

**Lemma A.2.** [14, Lemma 4 and proof of Theorem 2] Assume that all probability measures have compact supports included in $B(0, R)$, and that $n$ samples for $\nu$ are available. Then,

$$ \mathbb{E} \left[ \sup_{\varphi \in \mathcal{F}_R} \left| \int \varphi d(\nu - \hat{\nu}) \right| \right] \lesssim \begin{cases} R^2 n^{-1/2} & \text{if } d < 4, \\ R^2 n^{-1/2} \log(n) & \text{if } d = 4, \\ R^2 n^{-2/d} & \text{if } d > 4, \end{cases} \quad (A.6) $$

where $\lesssim$ hides a constant that depends only on $d$, and $\mathcal{F}_R$ denotes the class of concave and R-Lipschitz functions on $B(0, R)$. In the same paper, the authors established that

$$ \mathbb{E} \left[ \left| \int_{\mathcal{Y}} \|y\|^2 d(\nu - \hat{\nu})(y) \right| \right] \le 4R^2 n^{-1/2}. \quad (A.7) $$

As we make a repeating use of the upper bound in equation (A.6) we denote it by $\mathcal{E}(d, n)$. From now on

$$ \mathcal{E}(d, n) := \begin{cases} R^2 n^{-1/2} & \text{if } d < 4, \\ R^2 n^{-1/2} \log(n) & \text{if } d = 4, \\ R^2 n^{-2/d} & \text{if } d > 4. \end{cases} \quad (A.8) $$

The next proposition gives an upper bound of the estimation error that is independent of the regularization parameter $\lambda$.

**Proposition A.1.** Let $\lambda \ge 0$. Suppose that every probability measure considered has compact support included in $B(0, R)$.

(i) If $n$ samples from $\nu$ are available, then it holds that

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \nu) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})| \right] \lesssim \mathcal{E}(n, d). \quad (A.9) $$

(ii) If for each distribution $\mu_k$, $m_k$ samples are available, then

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_\lambda(\mu_\theta, \hat{\nu}) - W_\lambda(\hat{\mu}_\theta, \hat{\nu}) \right| \right] \lesssim K \mathcal{E}(\underline{m}, d), \quad (A.10) $$

where $\underline{m} = \min(m_1, \dots, m_K)$.
Proof. The key point is to exploit the alternative dual formulation of regularized OT that has been introduced in Section 2.2 of the article. Using relation (2.11) of the article, we remark that for any $\theta \in \Sigma_K$,

$$ 
\begin{aligned}
W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu}) &= \int_{\mathcal{Y}} \|y\|^2 d\nu(y) - \int_{\mathcal{Y}} \|y\|^2 d\hat{\nu}(y) + W_{\lambda}^s(\mu_{\theta}, \nu) - W_{\lambda}^s(\mu_{\theta}, \hat{\nu}) \\
&= \int_{\mathcal{Y}} \|y\|^2 d(\nu - \hat{\nu})(y) + W_{\lambda}^s(\mu_{\theta}, \nu) - W_{\lambda}^s(\mu_{\theta}, \hat{\nu}). \quad (A.11)
\end{aligned}
$$

Now, let us denote by $\varphi$ and $\hat{\varphi}$ two optimal dual potentials respectively associated to $W_{\lambda}^s(\mu_{\theta}, \nu)$ and $W_{\lambda}^s(\mu_{\theta}, \hat{\nu})$ when exploiting the semi-dual formulation 2.14 of the article. We can thus write

$$ 
\begin{aligned}
W_{\lambda}^{s}(\mu_{\theta}, \nu) - W_{\lambda}^{s}(\mu_{\theta}, \hat{\nu}) &= \int \varphi(x) d\mu_{\theta}(x) + \int \varphi^{s}(y) d\nu(y) \\
&\quad - \left( \int \hat{\varphi}(x) d\mu_{\theta}(x) + \int \hat{\varphi}^{s}(y) d\hat{\nu}(y) \right) \\
&= \int \varphi^{s}(y) d\nu(y) - \int \varphi^{s}(y) d\hat{\nu}(y) \\
&\quad + \int \varphi(x) d\mu_{\theta}(x) + \int \varphi^{s}(y) d\hat{\nu}(y) \\
&\quad - \left( \int \hat{\varphi}(x) d\mu_{\theta}(x) + \int \hat{\varphi}^{s}(y) d\hat{\nu}(y) \right) \\
&\le \int \varphi^{s}(y) d(\nu - \hat{\nu})(y),
\end{aligned}
$$

where the last inequality derives from the optimality of $\hat{\varphi}$ for the semi-dual formulation of $W_{\lambda}^s(\mu_{\theta}, \hat{\nu})$. A similar reasoning yields

$$ W_{\lambda}^{s}(\mu_{\theta}, \hat{\nu}) - W_{\lambda}^{s}(\mu_{\theta}, \nu) \le \int \hat{\varphi}^{s}(y) d(\hat{\nu} - \nu)(y). $$

As $\varphi^s$ and $\hat{\varphi}^s$ are both s-transform, Proposition 2.2 in the article ensures that both $\varphi^s$ and $\hat{\varphi}$ are concave and R-Lipschitz on $B(0, R)$. We deduce the upper bound

$$ |W_{\lambda}^{s}(\mu_{\theta}, \nu) - W_{\lambda}^{s}(\mu_{\theta}, \hat{\nu})| \le \sup_{\varphi \in \mathcal{F}_R} \left| \int_{\mathcal{Y}} \varphi d(\nu - \hat{\nu}) \right|, \quad (A.12) $$

where $\mathcal{F}_R$ denotes the class of concave and R-Lipschitz functions on $B(0, R)$. Taking the expectation of inequality (A.12), point (i) of Proposition A.1 follows from Lemma A.2.

Point (ii) of Proposition A.1 can be obtained with a similar reasoning. That is, we exploit the relation between $W_{\lambda}$ and $W_{\lambda}^s$. Moreover, Proposition 3.1 ensures that the optimal potentials associated to $W_{\lambda}^s$ can be chosen R-Lipschitz and concave. Performing the same computations as in point (i), and decomposing $\mu_{\theta}$ into a convex combination of $\mu_1, \dots, \mu_K$, we derive

$$ 
\begin{aligned}
|W_{\lambda}(\mu_{\theta}, \hat{\nu}) - W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| &\le \sup_{\varphi \in \mathcal{F}_R} \left| \int_{\mathcal{X}} \varphi d(\mu_{\theta} - \hat{\mu}_{\theta}) \right| + \left| \int_{\mathcal{X}} \|x\|^2 d(\mu_{\theta} - \hat{\mu}_{\theta}) \right| \\
&\le \sum_{k=1}^{K} \theta_k \left( \sup_{\varphi \in \mathcal{F}_R} \left| \int_{\mathcal{X}} \varphi d(\mu_k - \hat{\mu}_k) \right| + \left| \int_{\mathcal{X}} \|x\|^2 d(\mu_k - \hat{\mu}_k) \right| \right) \\
&\le \sum_{k=1}^{K} \left( \sup_{\varphi \in \mathcal{F}_R} \left| \int_{\mathcal{X}} \varphi d(\mu_k - \hat{\mu}_k) \right| + \left| \int_{\mathcal{X}} \|x\|^2 d(\mu_k - \hat{\mu}_k) \right| \right),
\end{aligned}
$$
using the crude upper bound $\theta_k \le 1$. This last upper bound being independent of $\theta$, we can write

$$ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \hat{\nu}) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})| \le \sum_{k=1}^K \left( \sup_{\varphi \in \mathcal{F}_R} \left| \int_X \varphi d(\mu_k - \hat{\mu}_k) \right| + \left| \int_X \|x\|^2 d(\mu_k - \hat{\mu}_k) \right| \right). $$

Taking the expectation of last inequality yields

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \hat{\nu}) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})| \right] \le \sum_{k=1}^K \mathbb{E} \left[ \left( \sup_{\varphi \in \mathcal{F}_R} \left| \int_X \varphi d(\mu_k - \hat{\mu}_k) \right| + \left| \int_X \|x\|^2 d(\mu_k - \hat{\mu}_k) \right| \right) \right]. $$

Next, applying Lemma A.2 to the probability distribution $\mu_k$, we obtain

$$ \mathbb{E} \left[ \sup_{\varphi \in \mathcal{F}_R} \left| \int_X \varphi d(\mu_k - \hat{\mu}_k) \right| \right] + \mathbb{E} \left[ \left| \int_X \|x\|^2 d(\mu_k - \hat{\mu}_k) \right| \right] \lesssim \mathcal{E}(d, m_k), $$

where $\mathcal{E}(d, m_k)$ is defined in equation (A.8). It follows that

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \hat{\nu}) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})| \right] \lesssim \sum_{k=1}^K \mathcal{E}(d, m_k). $$

Finally, introducing the notation $\underline{m} = \min(m_1, \dots, m_K)$, we derive

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \hat{\nu}) - W_\lambda(\hat{\mu}_\theta, \hat{\nu})| \right] \lesssim K \mathcal{E}(d, \underline{m}), $$

which gives the last inequality of Proposition A.1. $\square$

We now gather the pieces to prove Proposition 3.1. Taking the expectation of equation (A.2) and using Proposition A.1 gives

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_\lambda(\mu_\theta, \nu) - W_\lambda(\mu_\theta, \hat{\nu}) \right| \right] \lesssim \mathcal{E}(n, d) + K \mathcal{E}(\underline{m}, d). $$

Under the assumption of Proposition 3.1 from the main article that we have access to $n$ samples from each probability distribution $\mu_1, \dots, \mu_K$, as well as for $\nu$; last inequality reads

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_\lambda(\mu_\theta, \nu) - W_\lambda(\mu_\theta, \hat{\nu}) \right| \right] \lesssim K \mathcal{E}(n, d), $$

which is the result announced in Proposition 3.1 of the main article.

## A.3 Collecting existing results to prove the main Theorem

**Approximation error** Thanks to [21, Theorem 1] adapted to the squared Euclidean cost $c(x, y) = \|x - y\|^2$ (which is $R$-Lipschitz on $B(0, R)$ w.r.t. both its variables), we can control the impact of entropic regularization on the approximation of the value of the un-regularized OT cost.

**Proposition A.2.** [21, Theorem 1] Assume that $\mathcal{X}, \mathcal{Y}$ are compact subsets of $B(0, R)$. Then, it holds that

$$ 0 \le W_\lambda(\mu, \nu) - W_0(\mu, \nu) \le 2d\lambda \log \left( \frac{8 \exp(2) R^2}{\sqrt{d\lambda}} \right), \qquad (A.13) $$

and consequently

$$ \sup_{\theta \in \Sigma_K} |W_0(\mu_\theta, \nu) - W_\lambda(\mu_\theta, \nu)| \le B(\lambda) \quad \text{where} \quad B(\lambda) = 2d\lambda \log \left( \frac{8 \exp(2) R^2}{\sqrt{d\lambda}} \right) \qquad (A.14) $$

Notice that $B(\lambda)$ goes to zero when $\lambda \to 0$ at the speed

$$ B(\lambda) \sim_{\lambda \to 0} 2d\lambda \log (1/\lambda). $$
**Algorithm error** For $\mu = \sum_{i=1}^{n} a_i \delta_{x_i}$, and $\nu = \sum_{j=1}^{m} b_j \delta_{y_j}$ two discrete distributions, we denote by

$$W_{\lambda}^{(\ell)}(\mu, \nu) = \sum_{i=1}^{n} a_i \varphi_i^{(\ell)} + \sum_{j=1}^{m} b_j \psi_j^{(\ell)}, \qquad (A.15)$$

the approximation of the regularized OT cost $W_{\lambda}(\mu, \nu)$ that is returned by the Sinkhorn approximation after $\ell$ iterations. The variables $\varphi^{(\ell)}$ and $\psi^{(\ell)}$ denote the dual variables after $\ell$ iterations of the Sinkhorn algorithm. We thus consider the estimator used in our numerical experiments that is defined as

$$\hat{\theta}_{\lambda}^{(\ell)} = \arg \min_{\theta \in \Sigma_{\kappa}} W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (A.16)$$

The computational complexity of Sinkhorn algorithm has been studied in [14] and we remind the error after $\ell$ iterations of the Sinkhorn algorithm with respect to the regularized OT cost.

**Proposition A.3.** [14, Proposition 2]. Assume that $\lambda > 0$. For $\mu = \sum_{i=1}^{n} a_i \delta_{x_i}$ and $\nu = \sum_{j=1}^{m} b_j \delta_{y_j}$ two discrete distributions and a ground cost set to $c(x, y) = \|x - y\|^2$ on $\mathbb{R}^d$. The approximation of the regularized OT cost after $\ell$ iterations of the Sinkhorn algorithm satisfies:

$$|W_{\lambda}^{(\ell)}(\mu, \nu) - W_{\lambda}(\mu, \nu)| \le \frac{\|c\|_{\infty}^{2}}{\lambda \ell} \qquad (A.17)$$

where $\|c\|_{\infty} = \max_{(i,j)} \|x_i - y_j\|^2$.

**Remark A.1.** If the discrete distributions $\mu$ and $\nu$ both have their supports subsets of $B(0, R)$, it implies that $\max_{(i,j)} \|x_i - y_j\|^2 \le 4R^2$. And then, the quantity $\|c\|_{\infty}^2$ can be upper bounded by $\|c\|_{\infty}^2 \le 16R^4$.

**Proof of the main result** We now conclude the proof of our main result; that is Theorem 3.1 of the main paper. From Lemma A.1 we have the inequality

$$W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) \le 2 \sup_{\theta \in \Sigma_{\kappa}} |W_0(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \nu)| \\ + 2 \sup_{\theta \in \Sigma_{\kappa}} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \\ + 2 \sup_{\theta \in \Sigma_{\kappa}} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|.$$

Taking the expectation on both sides of this last inequality yields

$$\mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu)] \le 2 \sup_{\theta \in \Sigma_{\kappa}} |W_0(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \nu)| \\ + 2 \mathbb{E} \left[ \sup_{\theta \in \Sigma_{\kappa}} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right] \\ + 2 \sup_{\theta \in \Sigma_{\kappa}} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|.$$

Exploiting Proposition A.2, we have that

$$\sup_{\theta \in \Sigma_{\kappa}} |W_0(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \nu)| \le 2d\lambda \log \left( \frac{8 \exp(2) R^2}{\sqrt{d\lambda}} \right).$$

Under the assumption that $n$ samples are available from each probability measure, Proposition A.1 ensures the control

$$\mathbb{E} \left[ \sup_{\theta \in \Sigma_{\kappa}} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right] \lesssim K \mathcal{E}(n, d).$$
Finally, the algorithm error is upper bounded thanks to Proposition A.3 as follows:

$$ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})| \le \frac{16R^4}{\lambda\ell}. \quad (A.18) $$

Gathering the pieces together, we derive

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda}^{(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu)] \lesssim \lambda \log \left(\frac{1}{\lambda}\right) + \mathcal{E}(n, d) + \frac{1}{\lambda\ell}, $$

as claimed in Theorem 3.1 of the main paper.

## B Alternative bound for the estimation error of $W_{\lambda}$

This section aims at controlling the estimation error

$$ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| \quad (B.1) $$

with an upper bound of the magnitude $M_{\lambda}n^{-1/2}$, where $M_{\lambda}$ is a constant that depends on $\lambda$. The arguments that we use are very much inspired by the works [21, 14]. We will also see how the constants $M_{\lambda}$ involved in this upper bound depend on the regularizing parameter $\lambda$ with a power that depends on the dimension $d$. We now introduce the space of functions that we exploit in our analysis. For a bounded subset $\mathcal{Z}$ of $\mathbb{R}^d$, we shall denote by $\mathcal{C}^{\mathcal{K}}(\mathcal{Z})$ the set of $\mathcal{C}^{\mathcal{K}}$ functions on $\mathcal{Z}$ equipped with the norm $\|f\|_{\mathcal{K}} = \max_{|\kappa| \le \mathcal{K}} \|\partial^{\kappa} f\|_{\infty}$ and

$$ \mathcal{C}_{M}^{\mathcal{K}}(\mathcal{Z}) = \{f \in \mathcal{C}^{\mathcal{K}}(\mathcal{Z}) \mid \|f\|_{\mathcal{K}} \le M\}, \quad (B.2) $$

where the notation $|\kappa| \le \mathcal{K}$ denotes any multi-index $\kappa$ of differentiation of length $|\kappa|$ at most $\mathcal{K}$.

**Proposition B.1.** Suppose that every probability measure considered has compact support included in $B(0, R)$. Then, we can bound the estimation error (B.1) as follows:

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu}) \right| \right] \lesssim \frac{M_{\lambda}}{\sqrt{n}}, \quad (B.3) $$

with

$$ M_{\lambda} = M_d \max \left( R^2, \frac{R^{[d/2]+1}}{\lambda^{[d/2]}} \right), \quad (B.4) $$

and $M_d$ is a constant that depends only on $d$.

This proof is built upon the following lemmas. A first step is to study the regularity of the optimal potentials of the dual formulation (2.2) of the paper. To this end, we adapt the analysis of [21] to the setting where the measure $\mu$ belongs to the parametric model $\{\mu_{\theta} \mid \theta \in \Sigma_K\}$. This result is adapted from [21] and relies on the fact that an optimal potential can be chosen as the c-transform of $\varphi \in L^{\infty}(\mathcal{X})$. The choice of the squared Euclidean cost allows for a very clear description of the regularity of an optimal potential. We recall below the definition of the c-transform of $\varphi \in L^{\infty}(\mathcal{X})$ that we denote by $\varphi_{\mu}^{c, \lambda}$ as $\lambda > 0$. With the aim of manipulating functions defined on a convex and bounded subset of $\mathbb{R}^d$, the c-transforms are defined on $B(0, R)$. Hence, even if integrated only against $\mu$ or $\nu$ that have support $\mathcal{X}$ and $\mathcal{Y}$ respectively, the c-transform are defined on $B(0, R)$. The expression of $\varphi_{\mu}^{c, \lambda} \in L^{\infty}(B(0, R))$ is given by

$$ \forall y \in B(0, R), \quad \varphi_{\mu}^{c, \lambda}(y) = -\lambda \log \int e^{-\frac{\|x-y\|^2}{\lambda}} d\mu(x). \quad (B.5) $$
**Lemma B.1.** Suppose that every probability measure has compact support included in $B(0, R)$. Then, for all $\theta \in \Sigma_K$, there exists a couple of dual potentials $(\varphi, \psi)$ with respect to $W_\lambda(\mu_\theta, \nu)$, that satisfies $\psi(0) = 0$, and $\psi = \varphi_{\mu_\theta}^{c,\lambda}$. Moreover, $\psi$ belongs to $\mathcal{C}^\infty(B(0, R))$, and for each $\mathcal{H} > 0$, there exists a constant $M_{\mathcal{H}} > 0$ that depends only on $\mathcal{H}$ such that

$$ \|\psi\|_{\mathcal{H}} \le M_{\mathcal{H}} \max \left( R^2, \frac{R^{\mathcal{H}}}{\lambda^{\mathcal{H}-1}} \right). \qquad (B.6) $$

The sup norm is taken over $B(0, R)$.

This lemma will be proved in Section C. Two observations on this Lemma B.1 will reveal useful in the sequel.

**Remark B.1.** We stress that $\psi = \varphi_{\mu_\theta}^{c,\lambda}$ is a regularized c-transform with respect to $\mu_\theta$ (with $\theta \in \Sigma_K$), and that the constant $M_{\mathcal{H}} > 0$ that appears in Lemma B.1 does not depend on $\theta$.

**Remark B.2.** If we denote by $\varphi$ the dual potential of Lemma B.1 such that $\psi = \varphi_{\mu_\theta}^{c,\lambda}$ where $\psi$ meets the requirements of the Lemma B.1, this variable $\varphi$ can be chosen as an optimal potential for the semi-dual formulation given in equation (2.5) of the paper.

In the sequel of this section, we make a repeating use of the constant of Lemma B.1. Thus, we introduce the notation $M_{\lambda, \mathcal{H}}$ to refer to the upper bound of equation (B.6) which is defined as

$$ M_{\lambda, \mathcal{H}} := M_{\mathcal{H}} \max \left( R^2, \frac{R^{\mathcal{H}}}{\lambda^{\mathcal{H}-1}} \right), \qquad (B.7) $$

and

$$ M_{\mathcal{H}} > 0 \text{ is a constant that depends only on } \mathcal{H}. \qquad (B.8) $$

The next proposition links the estimation error $\sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \nu) - W_\lambda(\mu, \hat{\nu})|$ to the regularity of the c-transforms established in Lemma B.1.

**Proposition B.2.** Suppose that all probability measures have compact supports included in $B(0, R)$. Then, for $\mathcal{H} > 0$, we have the following upper bound

$$ \sup_{\theta \in \Sigma_K} |W_\lambda(\mu_\theta, \nu) - W_\lambda(\mu_\theta, \hat{\nu})| \le \sup_{\psi \in \mathcal{C}_{M_{\lambda, \mathcal{H}}}^{\mathcal{H}}(B(0, R))} \int \psi d(\nu - \hat{\nu}), \qquad (B.9) $$

with $M_{\lambda, \mathcal{H}} > 0$ a constant defined in equation (B.7).

**Proof.** Let $\theta \in \Sigma_K$, $\mathcal{H} > 0$ and introduce $\varphi, \psi$ (resp. $\hat{\varphi}, \hat{\psi}$) two optimal potentials for the dual formulation of $W_\lambda(\mu_\theta, \nu)$ (resp. $W_\lambda(\mu_\theta, \hat{\nu})$) chosen as in Lemma B.1. In particular $\psi, \hat{\psi} \in \mathcal{C}_{M_{\lambda, \mathcal{H}}}^{\mathcal{H}}(B(0, R))$, and the potentials $\varphi$ and $\hat{\varphi}$ are respectively optimal potentials for the semi-dual formulation of $W_\lambda(\mu_\theta, \nu)$ and $W_\lambda(\mu_\theta, \hat{\nu})$ as precised in Remark B.2. We can thus write

$$ \begin{aligned} W_\lambda(\mu_\theta, \nu) - W_\lambda(\mu_\theta, \hat{\nu}) &= \int \varphi d\mu_\theta + \int \psi d\nu - \int \hat{\varphi} d\mu_\theta - \int \hat{\psi} d\hat{\nu} \\ &= \int \psi d\nu - \int \psi d\hat{\nu} \\ &\quad + \underbrace{\left( \int \varphi d\mu_\theta + \int \psi d\hat{\nu} - \int \hat{\varphi} d\mu_\theta - \int \hat{\psi} d\hat{\nu} \right)}_{\le 0}. \end{aligned} $$

By optimality of $\hat{\varphi}$ for the semi dual formulation of $W_\lambda(\mu_\theta, \hat{\nu})$, the last term in the above parenthesis is non-positive. Using a symmetric optimality argument for $W_\lambda(\mu_\theta, \nu)$ and its optimal potential $\varphi$ for the semi-dual formulation, we get

$$ \int \hat{\psi} d(\nu - \hat{\nu}) \le W_\lambda(\mu_\theta, \nu) - W_\lambda(\mu_\theta, \hat{\nu}) \le \int \psi d(\nu - \hat{\nu}). \qquad (B.10) $$
As $\psi$ and $\hat{\psi}$ belong to $\mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))$, we can write

$$ 
\begin{aligned}
|W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| &\le \sup_{\psi \in \mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))} \left|\int \psi d(\nu - \hat{\nu})\right| \\
&= \sup_{\psi \in \mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))} \int \psi d(\nu - \hat{\nu}).
\end{aligned}
\quad (B.11) $$

The set $\mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))$ being independent of $\theta$, we get

$$ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| \le \sup_{\psi \in \mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))} \int \psi d(\nu - \hat{\nu}). \quad (B.12) $$

□

Therefore, the search for a control over $\sup_{\theta \in \Sigma_K} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu, \hat{\nu})|$ leads us to the study of the following empirical process

$$ \sup_{\psi \in \mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))} Y_{\psi} \quad \text{with} \quad Y_{\psi} = \int \psi d(\hat{\nu} - \nu) = \frac{1}{n} \sum_{i=1}^{n} \psi(Y_i) - \int \psi d\nu. \quad (B.13) $$

In order to bound the empirical process (B.13), we will need several ingredients. First, we show that this empirical process has a sub-Gaussian behavior (see the definition in [39]).

**Lemma B.2.** Under the assumption that all probability measures have compact supports included in $B(0, R)$, if $Y = \{Y_1, \dots, Y_n\}$ where $Y_1, \dots, Y_n$ are independent random samples from $\nu$, the empirical process $(Y_{\psi})_{\psi \in \mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))}$ defined in equation (B.13) has zero mean and is subgaussian w.r.t. $2n^{-\frac{1}{2}} \| \cdot \|_{\infty}$. In other terms, for all $\varphi, \psi \in \mathcal{C}_{M_{\lambda, \mathcal{X}}}(B(0, R))$ we have

$$ \forall s \in \mathbb{R}, \quad \mathbb{E}[e^{s(Y_{\varphi} - Y_{\psi})}] \le e^{\frac{2s^2}{n} \|\varphi - \psi\|_{\infty}^2} = e^{\frac{s^2}{2} (2n^{-\frac{1}{2}} \|\varphi - \psi\|_{\infty})^2}. \quad (B.14) $$

This lemma in an application of Azuma-Hoeffding inequality [39][Corollary 3.9]. We can now use Dudley's entropy integral inequality, which we recall below.

**Theorem B.1.** [39][Corollary 5.25] Let $(Y_{\varphi})_{\varphi \in \Phi}$ be a zero mean stochastic process which is sub-Gaussian with respect to the distance induced by a norm $\| \cdot \|$ on the indexing set $\Phi$. Then

$$ \mathbb{E} \left[ \sup_{\varphi \in \Phi} Y_{\varphi} \right] \le 12 \int_{0}^{\infty} \sqrt{\log(N(\varepsilon, \Phi, \| \cdot \|))} d\varepsilon, $$

where $N(\varepsilon, \Phi, \| \cdot \|)$ is the covering number of $\Phi$ by balls of radius $\varepsilon$ with respect to the norm $\| \cdot \|$.

A classical bound on the covering number for smooth functions (see e.g. [38, Theorem 2.7.1]) will prove highly valuable.

**Theorem B.2.** If $\mathcal{Z}$ is a bounded convex subset of $\mathbb{R}^d$ with nonempty interior, then there exists a constant $L(\mathcal{K}, d)$ such that

$$ \forall \varepsilon > 0, \quad \log N(\varepsilon, \mathcal{C}_1^{\mathcal{K}}(\mathcal{Z}), \| \cdot \|_{\infty}) \le L(\mathcal{K}, d) |\mathcal{Z} + B(0, 1)| \frac{1}{\varepsilon^{d/\mathcal{K}}}. \quad (B.15) $$

where $N(\varepsilon, \mathcal{C}_1^{\mathcal{K}}(\mathcal{Z}), \| \cdot \|_{\infty})$ denotes the covering number of $\mathcal{C}_1^{\mathcal{K}}(\mathcal{Z})$ (by balls of radius $\varepsilon$) with respect to the $\ell_{\infty}$ norm, and where $|\mathcal{Z} + B(0, 1)|$ is the Lebesgue measure of $\mathcal{Z} + B(0, 1)$

We now have all the ingredients to bound the expectation of the empirical process (B.13).
**Proposition B.3.** Suppose that all probability measures have compact supports included in $B(0, R)$. Then, we have the following upper bound

$$ \mathbb{E} \left[ \sup_{\psi \in \mathcal{C}_{M_{\lambda}}^{d'}(B(0,R))} \int \psi d(\hat{\nu} - \nu) \right] \lesssim \frac{M_{\lambda}}{\sqrt{n}} \quad \text{with} \quad M_{\lambda} = M_d \max \left( R^2, \frac{R^{\lfloor d/2 \rfloor + 1}}{\lambda^{\lfloor d/2 \rfloor}} \right), \quad (B.16) $$

and $M_d > 0$ is a constant that depends only on $d$.

*Proof.* Set $\mathcal{H} > 0$. We denote the empirical process under study by

$$ (Y_{\psi})_{\psi \in \mathcal{C}_{M_{\lambda}, \mathcal{H}}^{\mathcal{H}}(B(0,R))} \quad \text{with} \quad Y_{\psi} = \int \psi d(\hat{\nu} - \nu) = \frac{1}{n} \sum_{i=1}^{n} \psi(Y_i) - \int \psi d\nu, $$

and $M_{\lambda, \mathcal{H}} > 0$ the constant defined in equation (B.7). Dudley's inequality with the entropy integral (see Theorem B) gives

$$ \mathbb{E} \left[ \sup_{\psi \in \mathcal{C}_{M_{\lambda}, \mathcal{H}}^{\mathcal{H}}(B(0,R))} Y_{\psi} \right] \le 12 \int_{0}^{\infty} \sqrt{\log N(\varepsilon, \mathcal{C}_{M_{\lambda}, \mathcal{H}}^{\mathcal{H}}(B(0,R)), 2n^{-\frac{1}{2}} \| \cdot \|_{\infty})} d\varepsilon \quad (B.17) $$

$$ \le 12 \int_{0}^{\infty} \sqrt{\log N \left( \frac{1}{2} \sqrt{n} M_{\lambda, \mathcal{H}}^{-1} \varepsilon, \mathcal{C}_{1}^{\mathcal{H}}(B(0, R)), \| \cdot \|_{\infty} \right)} d\varepsilon \quad (B.18) $$

$$ \le \frac{24 M_{\lambda, \mathcal{H}}}{\sqrt{n}} \int_{0}^{\infty} \sqrt{\log N(\varepsilon, \mathcal{C}_{1}^{\mathcal{H}}(B(0, R)), \| \cdot \|_{\infty})} d\varepsilon \quad (B.19) $$

$$ \lesssim \frac{M_{\lambda, \mathcal{H}}}{\sqrt{n}} \int_{0}^{1/2} \varepsilon^{-\frac{d}{2\mathcal{H}}} d\varepsilon. \quad (B.20) $$

This integral is finite as soon as $\mathcal{H} > d/2$. As $M_{\lambda, \mathcal{H}} = M_{\mathcal{H}} \max(R^2, \frac{R^{\mathcal{H}}}{\lambda^{\mathcal{H}-1}})$ and $\lambda$ will be chosen little in the sequel, we set $\mathcal{H} = \lfloor d/2 \rfloor + 1$ in order to have the quantity $M_{\lambda, \mathcal{H}}$ as small as possible. From now on, we denote by $d' := \lfloor d/2 \rfloor + 1$, and with this choice of $\mathcal{H} = d'$, we can substitute the constant $M_{\lambda, \mathcal{H}}$ of Lemma B.1 with a new constant $M_{\lambda} := M_{\lambda, d'}$ that reads

$$ M_{\lambda} = M_{d} \max \left( R^{2}, \frac{R^{\lfloor d/2 \rfloor + 1}}{\lambda^{\lfloor d/2 \rfloor}} \right), $$

where $M_d$ is a constant that depends only on $d$. Finally, we have the following bound for the empirical process under study

$$ \mathbb{E} \left[ \sup_{\psi \in \mathcal{C}_{M_{\lambda}}^{d'}(B(0,R))} Y_{\psi} \right] \lesssim \frac{M_{\lambda}}{\sqrt{n}}. \quad (B.21) $$

□

**Proof of Proposition B.1** Gathering the results established since Lemma B.1, we are in a favorable position to prove Proposition B.1.

*Proof.* Indeed, by combining inequality (B.16) from Proposition B.3 with upper bound (B.9), we derive

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu}) \right| \right] \lesssim \frac{M_{\lambda}}{\sqrt{n}}. \quad (B.22) $$

And the above inequality is the result claimed in Proposition B.1.

□
The second empirical process in (A.2) can be upper bounded in a similar manner, as shown in the next proposition.

**Proposition B.4.** Suppose that all probability measures have compact support included in $B(0, R)$. Make the additional assumption that for all $k \in \{1, \dots, K\}$, $m_k$ samples from $\mu_k$ are available, and denote by $\underline{m} = \min(m_1, \dots, m_K)$. Then, the following inequality holds

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| \right] \lesssim \frac{KM_{\lambda}}{\sqrt{\underline{m}}}, \qquad (B.23) $$

where $M_{\lambda}$ is defined in equation (B.4).

*Proof.* We begin by setting $\theta \in \Sigma_K$. Let us denote by $\varphi$ (resp. $\hat{\varphi}$) an optimal potential chosen as in Lemma B.1 when considering the semi dual formulation of the regularized optimal transport problem between $\mu_{\theta}$ and $\hat{\nu}$ (resp. $\hat{\mu}_{\theta}$ and $\hat{\nu}$). Thus,

$$ \begin{aligned} W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \hat{\nu}) &= \int \hat{\varphi}_{\hat{\nu}}^{c, \lambda} d\hat{\mu}_{\theta} + \int \hat{\varphi} d\hat{\nu} - \left( \int \varphi_{\hat{\nu}}^{c, \lambda} d\mu_{\theta} + \int \varphi d\hat{\nu} \right) \\ &= \int \hat{\varphi}_{\hat{\nu}}^{c, \lambda} d\hat{\mu}_{\theta} - \int \hat{\varphi}_{\hat{\nu}}^{c, \lambda} d\mu_{\theta} \\ &\quad + \underbrace{\left( \int \hat{\varphi}_{\hat{\nu}}^{c, \lambda} d\mu_{\theta} + \int \hat{\varphi} d\hat{\nu} - \left( \int \varphi_{\hat{\nu}}^{c, \lambda} d\mu_{\theta} + \int \varphi d\hat{\nu} \right) \right)}_{\le 0} \end{aligned} $$

The optimality of the variable $\varphi$ with respect to the measures $\mu_{\theta}$ and $\hat{\nu}$ ensures the last term of the previous equation to be non positive. Hence,

$$ W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \hat{\nu}) \le \int \hat{\varphi}_{\hat{\nu}}^{c, \lambda} d(\hat{\mu}_{\theta} - \mu_{\theta}). $$

With a slight modification of the last argument we get

$$ W_{\lambda}(\mu_{\theta}, \hat{\nu}) - W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) \le \int \varphi_{\hat{\nu}}^{c, \lambda} d(\mu_{\theta} - \hat{\mu}_{\theta}). $$

Combining these last inequalities, we have

$$ \begin{aligned} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| &\le \sup_{\psi \in L^{\infty}(Y)} \left| \sum_{k=1}^{K} \theta_k \int \psi_{\hat{\nu}}^{c, \lambda} d(\mu_k - \hat{\mu}_k) \right| \\ &\le \sum_{k=1}^{K} \theta_k \sup_{\psi \in L^{\infty}(Y)} \left| \int \psi_{\hat{\nu}}^{c, \lambda} d(\mu_k - \hat{\mu}_k) \right| \\ &\le \sum_{k=1}^{K} \sup_{\psi \in L^{\infty}(Y)} \left| \int \psi_{\hat{\nu}}^{c, \lambda} d(\mu_k - \hat{\mu}_k) \right|. \end{aligned} $$

This last upper bound being independent of $\theta$, we derive

$$ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| \le \sum_{k=1}^{K} \sup_{\psi \in L^{\infty}(Y)} \left| \int \psi_{\hat{\nu}}^{c, \lambda} d(\mu_k - \hat{\mu}_k) \right|. \qquad (B.24) $$

Next, the application of Lemma B.1 when computing a c-transform w.r.t. $\hat{\nu}$ gives that $\psi_{\hat{\nu}}^{c, \lambda} \in \mathcal{C}_{M_{\lambda}}^{d'}(B(0, R))$ with $d'$ and $M_{\lambda}$ both defined in equation (B.4). Hence, for $k \in \{1, \dots, K\}$, using the
same ingredients as in Proposition B.2, we reach the study of an empirical process indexed by the class of functions $\mathcal{C}_{M_\lambda}^{d'}(B(0, R))$. Finally, the straight application of Proposition B.3 yields

$$ \mathbb{E} \left[ \sup_{f \in \mathcal{C}_{M_\lambda}^{d'}((B(0,R)))} \left| \int fd(\mu_k - \hat{\mu}_k) \right| \right] \lesssim \frac{M_\lambda}{\sqrt{m_k}}, $$

where $m_k$ is the number of observations sampled from distribution $\mu_k$. Then, exploiting inequality (B.24) we derive

$$ \begin{aligned} \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} \left| W_\lambda(\hat{\mu}_\theta, \hat{\nu}) - W_\lambda(\mu_\theta, \hat{\nu}) \right| \right] & \lesssim \sum_{k=1}^K \frac{M_\lambda}{\sqrt{m_k}} \\ & \le \frac{K M_\lambda}{\sqrt{\underline{m}}}, \end{aligned} $$

where $\underline{m} = \min\{m_k : 1 \le k \le K\}$. $\square$

We now gather the results from the previous section to obtain an upper bound on the expected excess risk of our regularized Wasserstein estimators.

**Proposition B.5.** Suppose that all probability measures have compact supports included in $B(0, R)$. If, for all $k \in \{1, \dots, K\}$, $m_k$ samples from $\mu_k$ are available and $n$ samples from $\nu$ are available, denoting $\underline{m} = \min(m_1, \dots, m_K)$, we have

$$ \mathbb{E} \left[ W_0(\mu_{\hat{\theta}_\lambda}, \nu) - W_0(\mu_{\theta^*}, \nu) \right] \lesssim \underbrace{\frac{2KM_\lambda}{\sqrt{\underline{m}}} + \frac{2M_\lambda}{\sqrt{n}}}_{\text{Estimation error}} + \underbrace{4d\lambda \log \left( \frac{8 \exp(2)R^2}{\sqrt{d\lambda}} \right)}_{\text{Approximation error}}, \quad (B.25) $$

where $M_\lambda = M_d \max \left( R^2, \frac{R^{\lfloor d/2 \rfloor + 1}}{\lambda^{\lfloor d/2 \rfloor}} \right)$, and $M_d$ is a constant that depends only on $d$.

*Proof.* Gathering the results on the approximation error for the regularized OT cost in Proposition A.2 and the upper bounds from Proposition B.1 and Proposition B.4 on the empirical processes defined in (A.2), we obtain the convergence rate claimed in (B.25). $\square$

From the upper bound on the expected excess risk of $\hat{\theta}_\lambda$ established in the previous Proposition B.5, we can propose a regularization choice in order to balance the estimation error and the approximation error. This regularization choice and the corresponding rate of convergence are given in the next Corollary.

**Corollary B.1.** Suppose that all probability measures have compact supports included in $B(0, R)$. Make the additional assumption that for all the distributions $\mu_k$ and for $\nu$, at least $n$ samples are available. Then, choosing $\lambda_n = n^{-1/(2\lfloor d/2 \rfloor + 2)}$ we get

$$ \mathbb{E} \left[ W_0(\mu_{\hat{\theta}_\lambda}, \nu) - W_0(\mu_{\theta^*}, \nu) \right] \lesssim n^{-\frac{1}{2\lfloor d/2 \rfloor + 2}} \log(n), \quad (B.26) $$

where $\lesssim$ hides a constant that depends on $R$ and $d$.

*Proof.* In order to drive the approximation term towards zero, the regularization parameter will converge towards 0. And in this case, $\lambda \log \left( \frac{8 \exp(2)R^2}{\sqrt{d\lambda}} \right) \sim_{\lambda \to 0} \lambda \log(\lambda^{-1})$. Next as we have assumed that all the distributions have $n$ samples, the estimation term equals $\frac{M_\lambda}{\sqrt{n}} = \frac{M_d R^{\lfloor d/2 \rfloor + 1}}{\lambda^{\lfloor d/2 \rfloor} \sqrt{n}}$. To balance these two terms we set $\lambda_n = n^{-\frac{1}{2\lfloor d/2 \rfloor + 2}}$, and with this choice of regularization parameter, there exists a constant $M_{R,d} > 0$ such that for $n$ sufficiently large,

$$ \mathbb{E} \left[ W_0(\mu_{\hat{\theta}_\lambda}, \nu) - W_0(\mu_{\theta^*}, \nu) \right] \le M_{R,d} n^{\frac{-1}{2\lfloor d/2 \rfloor + 2}} \log(n). \quad (B.27) $$

$\square$
C Proof of Lemma B.1

In this section we give a precise bound on the derivatives of a c-transform. Some results of the same flavor had already been established, for instance in [21, Lemma 1, Lemma 2]. The specificity of our result is to exploit the particular cost function $c(x, y) = \|x - y\|^2$ to give a precise description of the bound and to ensure that it is independent of the parameter $\theta$.

We precise the notations previously introduced in equation (B.2). For a multi-index $\kappa = (\kappa_1, \dots, \kappa_d) \in \mathbb{N}^d$, we denote by $|\kappa| = \sum_{i=1}^d \kappa_i$ and $D^\kappa$ the differential operator defined as follows

$$D^{\kappa} = \frac{\partial^{|\kappa|}}{\partial x_1^{\kappa_1} \dots \partial x_d^{\kappa_d}}. \qquad (C.1)$$

**Lemma C.1.** Set $\lambda > 0$. We assume that $\mu$ and $\nu$ have compact supports included in $B(0, R)$. Denoting by $\psi := \varphi_{\mu}^{c,\lambda}$ the c-transform of a given function $\varphi \in L^{\infty}(\mathcal{Y})$, for every multi-index $\kappa \in (\kappa_1, \dots, \kappa_d)$, there exists a constant $C_K$ that depends only on $\mathcal{K}$ such that

$$\sup_{y \in B(0, R)} |D^{\kappa} \varphi(y)| \le C_{\mathcal{K}} \frac{R^{\mathcal{K}}}{\lambda^{\mathcal{K}-1}}.$$

Notice that for every $\theta \in \Sigma_K$, the probability measure $\mu_\theta = \sum_{k=1}^K \theta_k \mu_k$ has compact support included $B(0, R)$. Therefore, the upper bound established in the previous Lemma C.1 holds for every measure $\mu_\theta$ with $\theta \in \Sigma_K$.

**Proof.** We proceed by induction on $\mathcal{K} \ge 1$ with inductive hypothesis $H(\mathcal{K})$: there exists a constant $C_K$ such that for every multi-index $\kappa$ with $|\kappa| = \mathcal{K} \sup_{y \in B(0, R)} |D^\kappa \varphi(y)| \le C_K \frac{R^{\mathcal{K}}}{\lambda^{\mathcal{K}-1}}$.

**Base case.** As $\psi = \varphi_{\mu}^{c,\lambda}$,

$$\forall y \in B(0, R), \quad \exp\left(-\frac{\psi(y)}{\lambda}\right) = \int_{\mathcal{X}} \exp\left(\frac{\varphi(x) - c(x, y)}{\lambda}\right) d\mu(x). \qquad (C.2)$$

As $\mu$ has compact support we can differentiate with respect to some $y_j$. From this we deduce that $\psi$ is differentiable with respect to $y_j$ and that for $y \in B(0, R)$,

$$-\frac{1}{\lambda} \frac{\partial \psi}{\partial y_j}(y) \exp\left(-\frac{\psi(y)}{\lambda}\right) = -\frac{1}{\lambda} \int_{\mathcal{X}} \frac{\partial c}{\partial y_j}(x, y) \exp\left(\frac{\varphi(x) - c(x, y)}{\lambda}\right) d\mu(x).$$

Taking the absolute value on both sides of last equality gives

$$\left| \frac{\partial \psi}{\partial y_j}(y) \right| \exp \left( -\frac{\psi(y)}{\lambda} \right) \le \int_{\mathcal{X}} \left| \frac{\partial c}{\partial y_j}(x, y) \right| \exp \left( \frac{\varphi(x) - c(x, y)}{\lambda} \right) d\mu(x).$$

The cost $c$ being the squared euclidean distance, we have $\left| \frac{\partial c}{\partial y_j}(x, y) \right| \le 4R$. As $\psi$ is the c-transform of $\varphi$, $\int_{\mathcal{X}} \exp\left(\frac{\varphi(x)+\psi(y)-c(x,y)}{\lambda}\right) d\mu(x) = 1$. From this we deduce

$$\left| \frac{\partial \psi}{\partial y_j}(y) \right| \le 4R,$$

which concludes the base case.

**Induction case.** Set $\mathcal{K} \ge 2$, and suppose that $H(1), \dots, H(\mathcal{K} - 1)$ hold true. Set $\kappa \in \mathbb{N}^d$ a multi-index with $|\kappa| = \mathcal{K}$. The computation of $D^\kappa$ relies on Faà di Bruno's formula [26] and the correspondence

$$\kappa = (\kappa_1, \dots, \kappa_d) \longleftrightarrow \underbrace{\{1, \dots, 1\}}_{\kappa_1} \underbrace{\{2, \dots, 2\}}_{\kappa_2} \dots \underbrace{\{d, \dots, d\}}_{\kappa_d}. \qquad (C.3)$$
Throughout this induction step, we will always implicitly refer to this one to one relation (C.3). For instance, when mentioning a partition of $κ$, even if not explicitly stated, we refer to the right-hand side of relation (C.3). And when referring to the differential operator $D^B$ with $B$ a subset of $\{1, \dots, 1, 2, \dots, 2, \dots, d, \dots, d\}$ we identify $B$ with the corresponding multi index on the left-hand side of relation (C.3).

Applying Faà di Bruno's formula on the right side of equation (C.2), as well as the differentiation theorem under the integral, gives

$$ 
\begin{aligned}
& D^{\kappa} \int_{\mathcal{X}} \exp \left( \frac{\varphi(x) - c(x, y)}{\lambda} \right) d\mu(x) && (C.4) \\
& = \int_{\mathcal{X}} \left( \sum_{\pi \in \mathcal{P}(\kappa)} \left( \frac{-1}{\lambda} \right)^{|\pi|} \exp \left( \frac{\varphi(x) - c(x, y)}{\lambda} \right) \prod_{B \in \pi} D^B c(x, y) \right) d\mu(x),
\end{aligned}
$$

where $\mathcal{P}(\kappa)$ denotes the collection of partitions of the right hand side of relation (C.3), and $|\pi|$ the number of sets that compose $\pi$. Applying the same formula on the left hand side of equation (C.2) yields

$$ 
\begin{aligned}
D^{\kappa} \exp \left(-\frac{\psi(y)}{\lambda}\right) &= \sum_{\pi \in \mathcal{P}(\kappa)} \left(-\frac{1}{\lambda}\right)^{|\pi|} \exp \left(-\frac{\psi(y)}{\lambda}\right) \prod_{B \in \pi} D^B \psi(y) && (C.5) \\
&= \left(-\frac{1}{\lambda} D^{\kappa} \psi(y) + \sum_{\pi \in \tilde{\mathcal{P}}(\kappa)} \left(-\frac{1}{\lambda}\right)^{|\pi|} \prod_{B \in \pi} D^B \psi(y)\right) \exp \left(-\frac{\psi(y)}{\lambda}\right),
\end{aligned}
$$

with $\tilde{\mathcal{P}}(\kappa)$ the collection of partitions of $κ$, when relying on correspondence (C.3), without the partition composed of the full set. Reminding that both quantities (C.5) and (C.4) are equals enables us to derive

$$ 
\begin{aligned}
D^{\kappa} \psi(y) &= \sum_{\pi \in \tilde{\mathcal{P}}(\kappa)} \left(-\frac{1}{\lambda}\right)^{|\pi|-1} \prod_{B \in \pi} D^B \psi(y) \\
&\quad + \int_{\mathcal{X}} \left( \sum_{\pi \in \mathcal{P}(\kappa)} \left(-\frac{1}{\lambda}\right)^{|\pi|-1} \prod_{B \in \pi} D^B c(x, y) \right) \gamma(x, y) d\mu(x),
\end{aligned}
$$

where $\gamma(x, y) = \exp\left(\frac{\varphi(x)+\psi(y)-c(x,y)}{\lambda}\right)$. Then, taking the absolute value on both sides of last equality gives

$$ 
\begin{aligned}
|D^{\kappa} \psi(y)| &\le \sum_{\pi \in \tilde{\mathcal{P}}(\kappa)} \lambda^{1-|\pi|} \prod_{B \in \pi} |D^B \psi(y)| && (C.6) \\
&\quad + \int_{\mathcal{X}} \left( \sum_{\pi \in \mathcal{P}(\kappa)} \lambda^{1-|\pi|} \prod_{B \in \pi} |D^B c(x, y)| \right) \gamma(x, y) d\mu(x).
\end{aligned}
$$

We begin by controlling the first sum of the right-hand side of last inequality. As $\pi$ is not the partition built from one block, for every block $B \in \pi$, $|B| \le |\kappa| - 1 = \mathcal{K} - 1$. We can thus apply the inductive hypothesis to every factor of $\prod_{B \in \pi} |D^B \psi(y)|$ to derive

$$ \lambda^{1-|\pi|} \prod_{B \in \pi} |D^B \psi(y)| \le \lambda^{1-|\pi|} \prod_{B \in \pi} C_{|B|} \frac{R^{|B|}}{|\lambda^{|B|-1}|}. $$
As $\pi$ is a partition of $\kappa$ composed of $|\pi|$ blocks, we have

$$ \prod_{B \in \pi} \frac{R^{|B|}}{\lambda^{|B|-1}} = \frac{R^{|\kappa|}}{\lambda^{|\kappa|-|\pi|}}. $$

From these previous computations we deduce that up to a multiplicative constant that depends only on $|\kappa| = \mathcal{K}$ the following inequality holds true

$$ \sum_{\pi \in \mathcal{P}(\kappa)} \lambda^{1-|\pi|} \prod_{B \in \pi} |D^B \psi(y)| \lesssim \lambda^{1-|\pi|} \frac{R^{|\kappa|}}{\lambda^{|\kappa|-|\pi|}} = \frac{R^{\mathcal{K}}}{\lambda^{\mathcal{K}-1}}. $$

We now control, the integral on the right-hand side of inequality (C.6). For $\pi$ a partition of $\kappa$ and $B$ a block of $\pi$, as $c$ is the squared euclidean distance, for all $x, y \in B(0, R)$, $|D^B(x, y)| \le 4R$. From this we deduce

$$ \sum_{\pi \in \mathcal{P}(\kappa)} \lambda^{1-|\pi|} \prod_{B \in \pi} |D^B c(x, y)| \le \sum_{\pi \in \mathcal{P}(\kappa)} \lambda^{1-|\pi|} (4R)^{|\pi|}. $$

The biggest term in the sum over the partitions of $\kappa$ is when considering the partition composed of singletons. In this case, $|\pi| = |\kappa| = \mathcal{K}$. Hence, up to a multiplicative constant that depends only on $\mathcal{K}$,

$$ \int_{\mathcal{X}} \left( \sum_{\pi \in \mathcal{P}(\kappa)} \lambda^{1-|\pi|} \prod_{B \in \pi} |D^B c(x, y)| \right) \gamma(x, y) d\mu(x) \lesssim \frac{R^{\mathcal{K}}}{\lambda^{\mathcal{K}-1}} \int_{\mathcal{X}} \gamma(x, y) d\mu(x) = \frac{R^{\mathcal{K}}}{\lambda^{\mathcal{K}-1}}. $$

To derive the last equality, remind that $\gamma(x, y) = \exp\left(\frac{\varphi(x)+\psi(y)-c(x,y)}{\lambda}\right)$; and as a consequence of equality (C.2), $\int_{\mathcal{X}} \gamma(x, y)d\mu(x) = 1$. We now have upper bounded all the terms of the left hand side of inequality (C.6). We thus derive that for all $y \in B(0, R)$,

$$ |D^{\kappa} \psi(y)| \lesssim \frac{R^{\mathcal{K}}}{\lambda^{\mathcal{K}-1}}. $$

As $\kappa$ is an arbitrary multi index with weights $|\kappa| = \mathcal{K}$, it shows $H(\mathcal{K})$. $\square$

## Conclusion of the proof of Lemma B.1

*Proof.* As $\psi = \varphi_{\mu}^{c, \lambda}$, The application of Proposition 12 from [17] ensures that a c-transform inherits the Lipschitz constant of the cost function. The cost function being $c(x, y) = \|x - y\|^2$ with $x, y \in B(0, R)$, we have that $\psi$ is $4R$-Lipschitz. And as $\psi(0) = \rho(0) = 0$, we can write,

$$ \forall y \in B(0, R), \quad \|\psi(y)\| \le 4R \|y\|. $$

Hence $\|\psi\|_{\infty} \le 4R^2$. Finally, Lemma C.1 gives a control of the derivatives of $\psi$. Indeed, for every $\kappa \in \mathbb{N}^d$ such that $|\kappa| = \mathcal{K}$,

$$ \forall y \in B(0, R), \quad |D^{\kappa} \psi(y)| \lesssim \frac{R^{\mathcal{K}}}{\lambda^{1-\mathcal{K}}}. $$

From this we deduce

$$ \|\psi\|_{\mathcal{H}} \lesssim \max \left( R^2, \frac{R^{\mathcal{H}}}{\lambda^{\mathcal{H}-1}} \right), $$

as claimed in Lemma B.1. $\square$
# D Extension of the results to the Sinkhorn divergence $S_\lambda$

We now study the collection of estimators $(\hat{\theta}_\lambda^{S(\ell)})_{\lambda>0}$ defined by

$$ \hat{\theta}_{\lambda}^{S(\ell)} := \arg \min_{\theta \in \Sigma_K} S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (D.1) $$

When not taking into account the computational cost, we use the notation $\hat{\theta}_\lambda^S$ for the estimator defined by

$$ \hat{\theta}_{\lambda}^{S} := \arg \min_{\theta \in \Sigma_K} S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}). \qquad (D.2) $$

In the next Section D.1, we introduce two additional assumptions in order to exploit the approximation result established in [14] between the Sinkhorn divergence $S_\lambda(\mu, \nu)$ and the Wasserstein distance $W_0(\mu, \nu)$. Indeed, Theorem 1 in [14], that we remind in Theorem D.1 of this paper, allows to bound $|S_\lambda(\mu, \nu) - W_0(\mu, \nu)|$ with a constant that depends on $\lambda$, the standard Fisher information $I(\mu)$, $I(\nu)$ of $\mu, \nu$, and $I(\mu, \nu)$ that is the Fisher information of the Wasserstein geodesic between $\mu$ and $\nu$ defined as in [14].

## D.1 Fisher information and approximation error of the Sinkhorn divergence $S_\lambda(\mu, \nu)$

We discuss conditions that ensure a control of the approximation error between the Sinkhorn divergence $S_\lambda(\mu, \nu)$ and $W_0(\mu, \nu)$.

**Theorem D.1.** [14, Theorem 1] Suppose that $\mu$ and $\nu$ have bounded densities and supports. Then, it holds that

$$ |S_{\lambda}(\mu, \nu) - W_{0}(\mu, \nu)| \leq \frac{\lambda^2}{4} \max\{I(\mu, \nu), (I(\mu) + I(\nu))/2\}, \qquad (D.3) $$

where $I(\mu)$ refers to the standard Fisher information of $\mu$, and $I(\mu, \nu)$ is the Fisher information of the Wasserstein geodesic between $\mu$ and $\nu$ as defined in [14].

First, we introduce sufficient conditions to ensure that the Fisher information of $\mu_\theta$ can be upper bounded without dependence on $\theta$.

**Assumption A.1.** The probability distributions $\mu_1, \dots, \mu_K$ have finite Fisher information with respective densities $f_1, \dots, f_K$ w.r.t. the Lebesgue measure, and all the components $\mu_1, \dots, \mu_K$ have disjoint supports $\mathcal{X}_1, \dots, \mathcal{X}_K$.

**Proposition D.1.** Suppose that assumption A.1 holds. Then, one has that

$$ \forall \theta \in \Sigma_K, I(\mu_\theta) \leq \max_{k \in \{1, \dots, K\}} I(\mu_k). \qquad (D.4) $$

*Proof.* For simplicity, we consider the case $d=1$. Let $\theta \in \Sigma_K$. Then, using the assumption that the components $\mu_k$ have disjoint supports, we decompose the Fisher information of $\mu_\theta$ as follows

$$ \begin{aligned} I(\mu_\theta) &= \int_{\mathcal{X}} \left( \frac{f'_\theta(x)}{f_\theta(x)} \right)^2 f_\theta(x) dx = \sum_{k=1}^{K} \theta_k \int_{\mathcal{X}_k} \left( \frac{f'_\theta(x)}{f_\theta(x)} \right)^2 f_k(x) dx \\ &= \sum_{k=1}^{K} \theta_k \int_{\mathcal{X}_k} \left( \frac{\theta_k f'_k(x)}{\theta_k f_k(x)} \right)^2 f_k(x) dx = \sum_{k=1}^{K} \theta_k I(\mu_k) \leq \max_{k \in \{1, \dots, K\}} I(\mu_k), \end{aligned} $$

which proves Inequality (D.4) for $d=1$. The case $d>1$ can be treated analogously. $\square$

Next, in order to bound the Fisher information of the Wasserstein geodesic between $\mu_\theta$ and $\nu$ with a constant independent of $\theta$, we adapt the assumptions of Proposition 1 from [14] to our needs.
**Assumption A.2.** The probability distribution $\nu$ is absolutely continuous with respect to the Lebesgue measure. Moreover, there exist two constants $m > 0$ and $L > 0$ such that for all $\theta \in \Sigma_K$ the Brenier potential $\varphi_\theta$ between $\mu_\theta$ and $\nu$ has a L-Lipschitz continuous Hessian satisfying $mId \le \nabla^2 \varphi_\theta$.

**Proposition D.2.** Suppose that Assumptions A.1 and A.2 hold. Then, we have the following inequality

$$ \forall \theta \in \Sigma_K, I(\mu_\theta, \nu) \le \frac{2}{m} \left( \max_{k \in \{1, \dots, K\}} I(\mu_k) + \frac{L^2}{3m^2} \right). \qquad (D.5) $$

*Proof.* A straight application of Proposition 1 from [14] gives that

$$ \forall \theta \in \Sigma_K, I(\mu_\theta, \nu) \le \frac{2}{m} \left( I(\mu_\theta) + \frac{L^2}{3m^2} \right). $$

Then, using $I(\mu_\theta) \le \max_{k \in \{1, \dots, K\}} I(\mu_k)$ from Proposition D.1 yields inequality (D.5). $\square$

As a consequence of Theorem D.1, we have the following result.

**Corollary D.1.** Suppose that Assumptions A.1 and A.2 hold. Then, we have that there exists a constant $M_I > 0$ such that

$$ \forall \theta \in \Sigma_K, |S_{\lambda}(\mu_{\theta}, \nu) - W_0(\mu_{\theta}, \nu)| \le M_I \lambda^2, \qquad (D.6) $$

where

$$ M_I = \frac{1}{4} \max \left\{ \frac{2}{m} \left( \max_{k \in \{1, \dots, K\}} I(\mu_k) + \frac{L^2}{3m^2} \right), \frac{\max_{k \in \{1, \dots, K\}} I(\mu_k) + I(\nu)}{2} \right\}. \qquad (D.7) $$

*Proof.* The combination of the upper bounds on $I(\mu_\theta)$ and $I(\mu_\theta, \nu)$ established in Proposition D.1 and in Proposition D.2 respectively, as well as the upper bound of Theorem D.1 yields (D.6). $\square$

## D.2 Estimation error of the Sinkhorn divergence $S_\lambda$ and rate of convergence

**Lemma D.1.** For $\lambda > 0$, the excess risk of the estimator $\hat{\theta}_\lambda^{S(\ell)}$ is upper bounded as follows:

$$ \begin{aligned} 0 \le W_0(\mu_{\hat{\theta}_\lambda^{S(\ell)}}, \nu) - W_0(\mu_{\theta^*}, \nu) \le & 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_\theta, \nu) - S_\lambda(\mu_\theta, \nu)| \quad (D.8) \\ & + 2 \sup_{\theta \in \Sigma_K} |S_\lambda(\mu_\theta, \nu) - S_\lambda(\hat{\mu}_\theta, \hat{\nu})| \\ & + 2 \sup_{\theta \in \Sigma_K} |S_\lambda(\hat{\mu}_\theta, \hat{\nu}) - S_\lambda^{(\ell)}(\hat{\mu}_\theta, \hat{\nu})|. \end{aligned} $$

Mutatis mutandis, this is exactly the same proof as when considering the regularized transport cost $W_\lambda$.

*Proof.* We begin with the decomposition

$$ \begin{aligned} W_0(\mu_{\hat{\theta}_\lambda^{S(\ell)}, \nu}) - W_0(\mu_{\theta^*}, \nu) = & W_0(\mu_{\hat{\theta}_\lambda^{S(\ell)}, \nu}) - S_\lambda(\mu_{\hat{\theta}_\lambda^{S(\ell)}, \nu}) + S_\lambda(\mu_{\hat{\theta}_\lambda^{S(\ell)}, \nu}) - S_\lambda(\hat{\mu}_{\hat{\theta}_\lambda^{S(\ell)}, \hat{\nu}}) \\ & + S_\lambda(\hat{\mu}_{\hat{\theta}_\lambda^{S(\ell)}, \hat{\nu}}) - S_\lambda^{(\ell)}(\hat{\mu}_{\hat{\theta}_\lambda^{S(\ell)}, \hat{\nu}}) + S_\lambda^{(\ell)}(\hat{\mu}_{\hat{\theta}_\lambda^{S(\ell)}, \hat{\nu}}) - S_\lambda^{(\ell)}(\hat{\mu}_{\theta^*, \hat{\nu}}) \\ & + S_\lambda^{(\ell)}(\hat{\mu}_{\theta^*, \hat{\nu}}) - S_\lambda(\hat{\mu}_{\theta^*, \hat{\nu}}) + S_\lambda(\hat{\mu}_{\theta^*, \hat{\nu}}) - S_\lambda(\mu_{\theta^*}, \nu) \\ & + S_\lambda(\mu_{\theta^*}, \nu) - W_0(\mu_{\theta^*}, \nu). \end{aligned} \qquad (D.9) $$

Let us focus on the right hand side of this last equation (D.9). The first and the last differences are controlled by the approximation error $\sup_{\theta \in \Sigma_K} |S_\lambda(\mu_\theta, \nu) - W_0(\mu_\theta, \nu)|$. The second and sixth differences can be upper by the estimation error $\sup_{\theta \in \Sigma_K} |S_\lambda(\mu_\theta, \nu) - S_\lambda(\hat{\mu}_\theta, \hat{\nu})|$. The third and
fifth differences are upper bounded by the algorithm error $\sup_{\theta \in \Sigma_K} |S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})|$.

It only remains to control $S_{\lambda}^{(\ell)}(\hat{\mu}_{\hat{\theta}^{S(\ell)}}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta^{*}}, \hat{\nu})$. However, $\hat{\theta}^{S(\ell)}$ minimizes the function $\theta \mapsto S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})$. Hence $S_{\lambda}^{(\ell)}(\hat{\mu}_{\hat{\theta}^{S(\ell)}}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta^{*}}, \hat{\nu}) \le 0$.

Going back to equation (D.9) and substituting every difference of the right hand side by its appropriate bound we derive

$$W_0(\mu_{\hat{\theta}^{S(\ell)}}, \nu) - W_0(\mu_{\theta^{*}}, \nu) \le 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \nu)| \\ + 2 \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \\ + 2 \sup_{\theta \in \Sigma_K} |S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|,$$

which is the result claimed in Lemma D.1.

□

**Proposition D.3.** Set $\lambda > 0$ and suppose that every probability measures have compact supports included in $B(0, R)$. Then, denoting by $\mathcal{E}$ the quantity introduced in equation (A.8), the following inequalities hold true.

(i) If n sample are available from $\nu$, then

$$\mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \hat{\nu}) \right] \lesssim \mathcal{E}(n, d).$$

(ii) If for every probability measure $\mu_k$, $m_k$ observations are available, then

$$\mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \hat{\nu}) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) \right] \lesssim K \mathcal{E}(\underline{m}, d),$$

where $\underline{m} = \min(m_1, \dots, m_K)$.

*Proof.* For concision, we only proof point (i). Set $\theta \in \Sigma_K$, and using the definition of $S_{\lambda}$ we write

$$\begin{aligned} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \hat{\nu})| &= |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu}) + \frac{1}{2}(W_{\lambda}(\hat{\nu}, \hat{\nu}) - W_{\lambda}(\nu, \nu))| \\ &\le |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| \\ &\quad + \frac{1}{2}(|W_{\lambda}(\hat{\nu}, \hat{\nu}) - W_{\lambda}(\nu, \hat{\nu})| + |W_{\lambda}(\nu, \hat{\nu}) - W_{\lambda}(\nu, \nu)|). \end{aligned}$$

From last inequality we deduce

$$\mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \hat{\nu})| \right] \le \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |W_{\lambda}(\mu_{\theta}, \nu) - W_{\lambda}(\mu_{\theta}, \hat{\nu})| \right] \quad (D.10) \\ + \frac{1}{2} \left( \mathbb{E} [ |W_{\lambda}(\hat{\nu}, \hat{\nu}) - W_{\lambda}(\nu, \hat{\nu})| ] + \mathbb{E} [ |W_{\lambda}(\nu, \hat{\nu}) - W_{\lambda}(\nu, \nu)| ] \right).$$

The first term is upper bounded, up to a multiplicative constant, by $\mathcal{E}(n, d)$ thanks to point (i) of Proposition A.1. To control the second term, we apply point (ii) of Proposition A.1 in the case $K = 1$ and $\mu_1 = \nu$. Hence the second term is also upper bounded by $\mathcal{E}(n, d)$. The last term is upper bounded by the same quantity with a similar reasoning. From this we deduce

$$\mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \hat{\nu})| \right] \lesssim \mathcal{E}(n, d),$$

as announced in Proposition D.3.

□
Under ad hoc assumptions, we can now control the expected excess risk of the estimator $\hat{\theta}_{\lambda}^{S(\ell)}$ defined by problem (D.1).

**Theorem D.2.** Set $\lambda > 0$. Suppose that all probability measures, that are $\mu_1, \dots, \mu_K$ and $\nu$, have compact supports; and that Assumptions A.1 and A.2 hold true. Assume that for all component $\mu_k$, as well as for $\nu$, at least $n$ observations are available. Then, the expected excess risk of the estimator $\hat{\theta}_{\lambda}^{S(\ell)}$ introduced in equation (D.1) is upper bounded as follows:

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda}^{S(\ell)}, \nu}) - W_0(\mu_{\theta^*}, \nu)] \lesssim \mathcal{E}(n, d) + \lambda^2 + \frac{1}{\lambda \ell}, $$

where the quantity $\mathcal{E}(n, d)$ is defined by formula (A.8).

Thanks to the upper bound established in Theorem D.2, we can propose a choice of the parameters $\lambda$ and $\ell$. For concision, we only address the case $d > 4$.

**Corollary D.2.** Suppose that all probability measures have compact supports included in $B(0, R)$ and that Assumptions A.1 and A.2 hold true. If at least $n$ samples are available for each probability measure, setting $\lambda_n = n^{-1/d}$ and $\ell_n = 64R^4 n^{3/d}$, the estimator $\hat{\theta}_{\lambda_n}^{S(\ell_n)}$ defined in equation D.1 admits the following rate of convergence:

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda}^{S(\ell)}, \nu}) - W_0(\mu_{\theta^*}, \nu)] \lesssim n^{-2/d}. $$

*Proof of Theorem D.2.* From Lemma D.1, we have the inequality

$$ W_0(\mu_{\hat{\theta}_{\lambda}^{S(\ell)}, \nu}) - W_0(\mu_{\theta^*}, \nu) \le 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \nu)| \\ + 2 \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \\ + 2 \sup_{\theta \in \Sigma_K} |S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|. $$

Taking the expectation on both sides of this last inequality yields

$$ \mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda}^{S(\ell)}, \nu}) - W_0(\mu_{\theta^*}, \nu)] \le 2 \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \nu)| \\ + 2 \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right] \\ + 2 \sup_{\theta \in \Sigma_K} |S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})|. $$

Thanks to Assumptions A.1 and A.2, we can exploit corollary D.1 to derive

$$ \sup_{\theta \in \Sigma_K} |W_0(\mu_{\theta}, \nu) - S_{\lambda}(\mu_{\theta}, \nu)| \le M_I \lambda^2. $$

Then, the triangular inequity yields

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right] \le \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right] \\ + \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right]. $$

Under the assumption that $n$ samples are available from each probability measure, Proposition D.3 ensures the control

$$ \mathbb{E} \left[ \sup_{\theta \in \Sigma_K} |S_{\lambda}(\mu_{\theta}, \nu) - S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu})| \right] \lesssim K \mathcal{E}(n, d). $$
Finally, to control the algorithm error, we write

$$
\begin{align*}
|S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})| &= |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - \frac{1}{2}(W_{\lambda}(\hat{\mu}_{\theta}, \hat{\mu}_{\theta}) + W_{\lambda}(\hat{\nu}, \hat{\nu}))| \\
&\quad - W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu}) + \frac{1}{2}(W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\mu}_{\theta}) + W_{\lambda}^{(\ell)}(\hat{\nu}, \hat{\nu}))| \\
&\le |W_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})| \\
&\quad + \frac{1}{2}|W_{\lambda}(\hat{\mu}_{\theta}, \hat{\mu}_{\theta}) - W_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\mu}_{\theta})| + \frac{1}{2}|W_{\lambda}(\hat{\nu}, \hat{\nu}) - W_{\lambda}^{(\ell)}(\hat{\nu}, \hat{\nu})|.
\end{align*}
$$

As we work with the squared euclidean distance as a ground cost function, from Proposition A.3 we deduce that each term in the last inequality is upper bounded by $16R^4(\lambda\ell)^{-1}$. Hence, the algorithm error is controlled as follows:

$$
\sup_{\theta \in \Sigma_K} |S_{\lambda}(\hat{\mu}_{\theta}, \hat{\nu}) - S_{\lambda}^{(\ell)}(\hat{\mu}_{\theta}, \hat{\nu})| \le \frac{32R^4}{\lambda\ell}.
$$

Gathering the pieces together, we derive

$$
\mathbb{E} [W_0(\mu_{\hat{\theta}_{\lambda}^{S(\ell)}}, \nu) - W_0(\mu_{\theta^{*}}, \nu)] \lesssim \lambda^2 + \mathcal{E}(n, d) + \frac{1}{\lambda \ell},
$$

as claimed in Theorem D.2.

<figure><img src="image_15.png" /></figure>