# Sample complexity of optimal transport barycenters with discrete support

Léo Portales, Edouard Pauwels, Elsa Cazelles

► To cite this version:

Léo Portales, Edouard Pauwels, Elsa Cazelles. Sample complexity of optimal transport barycenters with discrete support. 2025. hal-05087164v2

HAL Id: hal-05087164

https://hal.science/hal-05087164v2

Preprint submitted on 18 Nov 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sample complexity of optimal transport barycenters with discrete support

Léo Portales

Edouard Pauwels

IRIT, TSE, INP Toulouse

TSE, Université Toulouse Capitole

leo.portales@irit.fr

edouard2.pauwels@ut-capitole.fr

Elsa Cazelles

CNRS, IRIT, Université de Toulouse

elsa.cazelles@irit.fr

October 2025

## Abstract

Computational implementation of optimal transport barycenters for a set of target probability measures requires a form of approximation, a widespread solution being empirical approximation of measures. We provide an $O(\sqrt{N/n})$ statistical generalization bound for the empirical sparse optimal transport barycenters problem, where $N$ is the maximum cardinality of the barycenter (sparse support) and $n$ is the sample size of the target measures empirical approximation. Our analysis includes various optimal transport divergences including Wasserstein, Sinkhorn and Sliced-Wasserstein. We discuss the optimality of our bound as well as the application of our results to specific settings such as K-means, constrained K-means, free and fixed support Wasserstein barycenters.

## 1 Introduction

Optimal transport barycenters provide a geometrically structured notion of mean for a set of probability measures. They are minimizers of functionals of the form

$$ \nu \mapsto \frac{1}{L} \sum_{\ell=1}^{L} \mathbb{D}(\mu^{\ell}, \nu), \qquad (1.1) $$

where $\nu$ is a probability measure, $\mathbb{D}$ is an optimal transport divergence, and $\{\mu^\ell\}_{\ell=1}^L$ is a collection of target probability measures in $\mathbb{R}^d$. This problem was first introduced for $\mathbb{D} = W_2^2$ [1], as a generalization of McCann's interpolation [42]. This was also proposed independently for texture synthesis and mixing [52]. Following fast implementation of $W_2$-barycenters (e.g. [19]), they became a central tool in machine learning and statistics with various applications such as sample blending [57], texture mixing [52], multisource domain adaptation [45] or music genre recognition [44]. Various extensions based on different transport divergences $\mathbb{D}$ in (1.1) exist: entropic [16], Bures [14], sliced [9] and Gromov Wasserstein barycenters [12].

**Sample averaged approximation and sparse barycenters** In many applications, the target probability measures in (1.1) are computationally intractable, and one needs to resort to discrete approximations. Typical schemes include histogram approximation (e.g. images), sample average approximation (SAA) [45, 17] or online stochastic approximation (SA) [23, 22, 6]. We focus on sample average approximation arising in two main use cases: the statistical setting where distributions are only accessible through empirical samples [45] and numerical Monte-Carlo approximation of known densities [51]. In both cases, all measures, including the resulting barycenter, are discrete, possibly supported on a very large number of atoms. Consequently, most optimal transport barycenter algorithms restrict the number of atoms of the barycenter, similarly to optimal quantization of measures [48]. Our analysis includes this support constraint. In other words, the probability measures $\nu$ in (1.1), is subject to the constraint $|\text{Supp}(\nu)| \le N$ for some integer $N \ge 1$, where $\text{Supp}(\nu)$ denotes the support of $\nu$. This problem is known as optimal transport barycenters with sparse support [61, 10].
Many modern machine learning problems leverage this sparse barycenter framework, including applications in LLM [2], flow cytometry [28], graph clustering [25] or cross-machine fault diagnosis [38].

Our main problem of interest is that of statistical guarantees for the optimal transport barycenters sample averaged approximation under support cardinality constraints, with various divergences $D$: (sliced) Wasserstein distances, (debiased) Sinkhorn divergences. The collection of target probability measures $\{\mu^\ell\}_{\ell=1}^L$ in $\mathbb{R}^d$ each relate to $n$ i.i.d random variables $X_1^\ell, \dots, X_n^\ell \stackrel{\text{i.i.d}}{\sim} \mu^\ell$ with empirical measures by $\mu_n^\ell := \frac{1}{n} \sum_{i=1}^n \delta_{X_i^\ell}$, and the empirical sparse barycenters are given by

$$ \underset{(Y, \pi) \in (\mathbb{R}^d)^N \times \Delta_N}{\operatorname{argmin}} \quad \frac{1}{L} \sum_{\ell=1}^{L} \mathbb{D} \left( \mu_n^\ell, \sum_{i=1}^{N} \pi_i \delta_{y_i} \right) \qquad (1.2) $$

where $Y := (y_1, \dots, y_N)$ is a point cloud in $\mathbb{R}^d$ and $\Delta_N$ is the $(N-1)$-probability simplex.

**Statistical guarantees for optimal transport barycenters** For general population measures, statistical rates of the form $O(n^{-1/3d})$ for empirical barycenters were described in [13], suffering from the curse of dimensionality. On the other hand, for discrete target measures, an $O(n^{-1/2})$ rate is known [34]. This is reminiscent of the statistical complexity of optimal transport adapting to the least complex measure [35]. Indeed, empirical optimal transport converges as $O(n^{-1/d})$ [27], while if at least one of the two measures is discrete, the rate is $O(n^{-1/2})$ [20].

**Contribution** We provide uniform rates in expectation for the empirical sparse optimal transport barycenter functional (1.2). They are of the form $O(\sqrt{N/n})$, where the barycenter has $N$ support points, and $n$ is the number of samples per empirical target measure. We consider several divergences bridging the gap between various results in the literature. We deduce sample complexity estimates for the optimal transport barycenters with sparse support. Our proof techniques combine tools from empirical risk minimization theory and estimations on dual variables arising in the dual formulations of transport divergences.

**Notations** Throughout this article we will denote respectively $M_1(\Omega)$ and $M_1^N(\Omega)$, for some $\Omega \subset \mathbb{R}^d$, the set of probability measures on $\Omega$ and the set of probability measures supported on at most $N$ points in $\Omega$. We will denote $\mathbb{B}_R$ the closed ball of $\mathbb{R}^d$ centered in 0 and of radius $R > 0$. Additionally, $\| \cdot \|$ is the Euclidean norm on $\mathbb{R}^d$, $\Delta_N$ the $N-1$ probability simplex and $\mathbb{S}^{d-1}$ the $(d-1)$-sphere in $\mathbb{R}^d$. For a measure $\mu$ and a measurable map $T$ we denote by $T_\mu \mu$ the pushforward (or image measure) of $\mu$ by $T$.

**Organization of the paper** We introduce technical material on optimal transport divergences in Section 2. Computational aspects are discussed in Section 2.3. Our main results are presented in Section 3, with a discussion in Section 4. Dual formulations of the divergences considered, together with technical details and proofs, are postponed to the appendices.

# 2 Introduction to optimal transport and barycenters

## 2.1 Optimal transport divergences

We refer to [49] for a computational point of view, and to [58] for a more theoretical presentation. The optimal transport between two measures $\mu, \nu \in M_1(\mathbb{R}^d)$ is a solution to the optimization problem over measures:

$$ \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} c(x, y) d\gamma(x, y), \qquad (2.1) $$

where $c: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is a lower semi-continuous cost function and $\Pi(\mu, \nu)$ is the set of couplings: probability measures with marginals $\mu$ and $\nu$. In the Euclidean case, if $c: (x, y) \mapsto \|x - y\|^p$ for an integer $p \ge 1$, this is the Wasserstein distance $W_p$ on $\mathbb{R}^d$:

$$ W_p^p(\mu, \nu) = \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x - y\|^p d\gamma(x, y). \qquad (2.2) $$
The semi-discrete setting corresponds to one measure being finitely supported. We focus on this setting which encompasses the sparse barycenter problem. The infinite dimensional problem over measures can be reformulated as a finite-dimensional maximization problem based on optimal transport duality (see Appendix A). Various optimal transport-based divergences were introduced based on variations of the Wasserstein distance.

The Sinkhorn divergences [18, 49] add an entropic regularization term, for $\mu, \nu \in \mathcal{M}_1(\mathbb{R}^d)$:

$$W_{\epsilon,p}^{p}(\mu, \nu) = \inf_{\gamma \in \Pi(\mu, \nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x - y\|^p d\gamma(x, y) + \epsilon \text{KL}(\gamma|\mu \otimes \nu) \quad (2.3)$$

where KL($\gamma|\xi$) := $\int_{\mathbb{R}^{2d} \times \mathbb{R}^{2d}} (\log(\frac{d\gamma}{d\xi}(x, y)) - 1)d\gamma(x, y)$ denotes the Kullback-Leibler divergence and $\epsilon > 0$ is a regularization parameter. For discrete measures, the optimization problem in (2.3) can then be efficiently solved using the Sinkhorn algorithm [55, 18]. Similarly as Wasserstein distances, in the semi-discrete setting, there is a finite dimensional dual problem [30] of particular use in our analysis, it is presented in Appendix A.

Considering projections onto lines leads to sliced Wasserstein distances [52, 37]. For any $\theta \in \mathbb{S}^{d-1}$, set $P_\theta$ : $x \in \mathbb{R}^d \mapsto \langle\theta, x\rangle \in \mathbb{R}$, understood as the projection on the line directed by $\theta$. For $\mu, \nu \in \mathcal{M}_1(\mathbb{R}^d)$, one obtains two collections of pushforward univariate measures $(P_{\theta_\sharp\mu})_{\theta\in\mathbb{S}^{d-1}}$ and $(P_{\theta_\sharp\nu})_{\theta\in\mathbb{S}^{d-1}}$. Sliced distances combine divergences between these univariate measures for all directions $\theta \in \mathbb{S}^{d-1}$. We obtain the Sliced and max-Sliced Wasserstein distances respectively defined as

$$\begin{aligned} SW_p^p(\mu, \nu) &= \int_{\mathbb{S}^{d-1}} W_p^p(P_{\theta_\sharp\mu}, P_{\theta_\sharp\nu}) d\sigma(\theta) \\ \max SW_p^p(\mu, \nu) &= \max_{\theta \in \mathbb{S}^{d-1}} W_p^p(P_{\theta_\sharp\mu}, P_{\theta_\sharp\nu}), \end{aligned} \quad (2.4)$$

where $\sigma$ denotes the uniform probability measure over $\mathbb{S}^{d-1}$. In the semi-discrete setting, these also admits finite dimensional dual problems as discussed in Appendix A.

## 2.2 Optimal transport barycenters with sparse support

Given a collection of target probability measures $\{\mu^\ell\}_{\ell=1}^L$, an unconstrained barycenter with respect to divergence $\mathbb{D}$ is a minimizer of (1.1), it is typically absolutely continuous [1, 11] if at least one of the target measures is absolutely continuous. A common strategy is approximation by discrete probability measures, restricting the optimization set in (1.1) to $\mathcal{M}_1^N(\mathbb{R}^d)$. This is referred to as optimal transport barycenter problem with sparse support, with underlying measures of the form $\nu^* := \sum_{i=1}^N \pi_i \delta_{y_i}$ where $Y := (y_1, \dots, y_N) \in (\mathbb{R}^d)^N$ and $\pi \in \Delta_N$. The problem is reformulated as follows:

$$\min_{(Y, \pi) \in \mathcal{A}} \mathrm{F}_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) := \min_{(Y, \pi) \in \mathcal{A}} \frac{1}{L} \sum_{\ell=1}^L \mathrm{D} \left( \mu^\ell, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \quad (2.5)$$

where $\mathcal{A}$ is a closed nonempty subset of $(\mathbb{R}^d)^N \times \Delta_N$ which integrates possibly additional constraints. Formulation (2.5) covers various settings described in the literature [19, 10]:

<table><thead><tr><th>A</th><th>L = 1</th><th>L > 1</th></tr></thead><tbody><tr><td>(R<sup>d</sup>)<sup>N</sup> × Δ<sub>N</sub></td><td>Optimal quantization</td><td>Barycenter with sparse support</td></tr><tr><td>(R<sup>d</sup>)<sup>N</sup> × {π&#773;}</td><td>Constrained quantization</td><td>Free support barycenter</td></tr><tr><td>{Y&#773;} × Δ<sub>N</sub></td><td>—</td><td>Fixed support barycenter</td></tr></tbody></table>

## 2.3 Discretization and computational complexity

Even for continuous target measures, barycenter algorithms operate on discrete data. A common approach is to approximate the target measures by empirical averages, and to compute their barycenter (1.1) with dedicated methods [19, 39, 34]. Most algorithms actually add support constraints as in (2.5), effectively computing a sparse barycenters [61, 10, 43]. For discrete measures supported on $N$, the optimal transport problem has an $O(N^3 \log(N))$ computational cost [49]. The discrete Wasserstein barycenter problem inherits this complexity [4, 10], making it tractable only for measures supported on a small number of atoms [3]. The complexity of Sinkhorn divergences $W_{\epsilon,p}^p$ (2.3) is
$O(N^2/\epsilon)$[24], which allows for a trade-off between optimal transport accuracy and algorithmic efficiency [16]. Sliced-Wasserstein distances are based on univariate closed form quantile expressions [59, Chapter 2.2], making them also computationally relevant [9].

### 3 Sample complexity for sparse optimal transport barycenters

This section presents our main results: statistical generalization bounds for optimal transport barycenters. We draw inspiration from [33, Section 4.4], which treats the K-means functional (see Section 4.4). The following provides uniform generalization bounds for the empirical barycenter cost function, which most barycenter algorithms seek to minimize in practice (see Section 2.3).

**Theorem 3.1 (Generalization error bound).** Let $\mu^1, \dots, \mu^L \in \mathcal{M}_1(\mathbb{B}_R)$, for $R > 0$. For some integer $n$, let $\mu_n^1, \dots, \mu_n^L$ be empirical measures supported over $n$ i.i.d random variables of respective law $\mu^\ell$, for all $\ell \in [1, L]$. Let $\mathbb{D} = W_p^p, W_{\epsilon,p}^p$ (for some $\epsilon > 0$), $SW_p^p$ or $\max-SW_p^p$. Then for any integer $N \ge 1$, and $F_{\mathbb{D}}$ the function defined in (2.5), we have

$$ \mathbb{E} \left[ \sup_{(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N} \left| F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) - F_{\mathbb{D}} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right| \right] \le C_{p,R} \sqrt{\frac{C_{d,N}}{n}}, $$

where the finite constant values $C_{p,R}$ and $C_{d,N}$ depend on the divergence $\mathbb{D}$ and are summarized in Table 1.

<table><caption>Table 1: Explicit constants arising in the generalization error bound for the divergences $\mathbb{D}$.</caption><thead><tr><th>$\mathbb{D}$</th><th>$C_{p,R}$</th><th>$C_{d,N}$</th></tr></thead><tbody><tr><td>$W_p^p$</td><td>$8\sqrt{2} \int_0^{2(2R)^p} \sqrt{\log\left(2 + \frac{32pR^p}{\tau}\right)} d\tau$</td><td>$N(d + 1)$</td></tr><tr><td>$\max-SW_p^p$</td><td>$8\sqrt{2} \int_0^{2(2R)^p} \sqrt{\log\left(2 + \frac{48pR^p}{\tau}\right)} d\tau$</td><td>$N(d + 1) + d$</td></tr><tr><td>$W_{\epsilon,p}^p$</td><td>$8\sqrt{2} \int_0^{(4p+1)(2R)^p} \sqrt{\log\left(2 + \frac{64p(2R)^p}{\tau}\right)} d\tau$</td><td>$N(d + 1)$</td></tr></tbody></table>

**Remark 3.2 (On the hypothesis of Theorem 3.1).** We consider empirical samples of equal size $n$ for each measure, as well as uniform weights $\frac{1}{L}$ for simplicity. Our results hold for different sample sizes, replacing $n$ by the minimal sample size, and any weights $\lambda \in \Delta_L$.

*Sketch of proof.* Complete details are provided in Appendix C. We treat here the case $\mathbb{D} = W_p^p$ for simplicity. Let $(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N$, the goal is to control, for each $\ell = 1, \dots, L$, the variation between empirical and population $W_p^p$ divergences. This is achieved by resorting to the dual formulation of $W_p^p$, in which the dual variables are vectors $w \in \mathbb{R}^N$. Lemma D.2 ensures that there exist functions $f_Y^w$ such that for each $\ell$ and any $\pi \in \Delta_N$

$$ \left| W_p^p \left( \mu_n^\ell, \sum_{i=1}^N \pi_i \delta_{y_i} \right) - W_p^p \left( \mu^\ell, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right| \le \sup_{w \in \mathbb{R}^N} \left| \mathbb{E}_{X \sim \mu_n^\ell} [f_Y^w(X)] - \mathbb{E}_{X \sim \mu^\ell} [f_Y^w(X)] \right|. $$

Proposition B.1 provides a bound on the dual variables, $w \in [-K, K]^N$ for some constant $K$ independent of $Y$ and $\pi$. This allows us to invoke standard results from empirical process theory (e.g., [32, Theorem 3.5.1]). $\square$

**Remark 3.3.** In the specific case where $\mathbb{D}$ is a distance, using the triangular inequality we have $|\mathbb{D}(\mu_n^\ell, \nu) - \mathbb{D}(\mu^\ell, \nu)| \le \mathbb{D}(\mu_n^\ell, \mu^\ell)$. Statistical results for $\mathbb{D}$ [27] could then be used to derive error estimates. However, for $\mathbb{D} = W_1$ for example, the estimate is of order $n^{-1/d}$ [15, Proposition 2.1], compared with $n^{-1/2}$ in Theorem 3.1.

**Corollary 3.4 (Sparse optimal transport barycenters).** Let $\mu^1, \dots, \mu^L \in \mathcal{M}_1(\mathbb{B}_R)$, for $R > 0$. For some integer $n$, let $\mu_n^1, \dots, \mu_n^L$ be empirical measures supported over $n$ i.i.d random variables of respective law $\mu^\ell$, for all $\ell \in [1, L]$.
Let also $N \ge 1$ be some integer, A a closed nonempty subset of $\mathbb{B}_R^N \times \Delta_N$, and let $\mathbb{D} = W_p^p$, $W_{\epsilon,p}^p$ (for some $\epsilon > 0$), $SW_p^p$ or max-$SW_p^p$. We denote by $(Y_n, \pi_n)$ an argmin in A of the empirical barycenter cost $(Y, \pi) \mapsto F_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i})$ where $F_{\mathbb{D}}$ is given in (2.5). Then we have:

$$ \mathbb{E} \left[ F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i^n \delta_{y_i}^n \right) - \min_{(Y, \pi) \in A} F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right] \le 2C_{p,R} \sqrt{\frac{C_{d,N}}{n}}, \quad (3.1) $$

where the constant $C_{p,R}$ and $C_{d,N}$ are summarized in the Table 1.

**Remark 3.5.** For $\mathbb{D} = W_2^2$ and $L = 1$, the estimation error in Corollary 3.4 is upper bounded by $O(N/\sqrt{n})$ [8], independently of the dimension $d$, to be compared to our $O(\sqrt{dN/n})$ results.

*Proof.* This leverages classical statistical learning theory arguments. By continuity of the functionals $(Y, \pi) \mapsto F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i})$ and $(Y, \pi) \mapsto F_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i})$, both admit respective global minimizers in A, $(Y^*, \pi^*)$ and $(Y_n^*, \pi_n^*)$, with associated discrete measures $\nu^* := \sum_{i=1}^N \pi_i \delta_{y_i}$ and $\nu_n^* := \sum_{i=1}^N \pi_i^n \delta_{y_i}^n$. Then we have

$$ \begin{aligned} & F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu_n^*) - F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu^*) \\ & \le F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu_n^*) - F_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \nu_n^*) + F_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \nu^*) - F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu^*) \\ & \le 2 \sup_{(Z, \tau) \in A} \left| F_{\mathbb{D}} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \tau_i \delta_{z_i} \right) - F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \tau_i \delta_{z_i} \right) \right|. \end{aligned} $$

Taking the expectation on both sides and applying Theorem 3.1 allows to conclude. □

Theorem 3.1 and the sample complexity results of Sinkhorn divergences in [29] also yield an upper bound for the debiased Sinkhorn barycenters [31, 36], although with constants which have an exponential dependence in the dimension. The proof of the next corollary is also postponed in Appendix C.

**Corollary 3.6 (Debiased Sinkhorn barycenters).** Let $\mu^1, \dots, \mu^L \in \mathcal{M}_1(\mathbb{B}_R)$, for $R > 0$. For some integer $n$, let $\mu_n^1, \dots, \mu_n^L$ be empirical measures supported over $n$ i.i.d random variables of respective law $\mu^\ell$, for all $\ell \in [1, L]$. Let also $N \ge 1$ be an integer, A a closed nonempty subset of $\mathbb{B}_R^N \times \Delta_N$ and for any $\epsilon > 0$, set $\overline{W}_{\epsilon,p}^p(\mu, \nu) = W_{\epsilon,p}^p(\mu, \nu) - \frac{1}{2}(W_{\epsilon,p}^p(\mu, \mu) + W_{\epsilon,p}^p(\nu, \nu))$. We denote by $(Y_n, \pi_n)$ an argmin in A of the empirical barycenter cost $(Y, \pi) \mapsto F_{\overline{W}_{\epsilon,p}^p}(\mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i})$ where $F_{\overline{W}_{\epsilon,p}^p}$ is given in (2.5). Then we have

$$ \mathbb{E} \left[ F_{\overline{W}_{\epsilon,p}^p} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i^n \delta_{y_i}^n \right) - \min_{(Y, \pi) \in A} F_{\overline{W}_{\epsilon,p}^p} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right] \le \frac{C_{\overline{W}}}{\sqrt{n}} \quad (3.2) $$

with $C_{\overline{W}} = C_{p,R}N(d+1) + 2Ke^{\frac{K}{\epsilon}}(1 + e^{-[d/2]})$ with $K = 2\sqrt{2p}(2R)^{p-1}R + (2R)^p$, $K$ is a constant depending only on $R$ and $p$, and $C_{p,R}$ is the constant related to $W_{\epsilon,p}^p$ in Theorem 3.1.

# 4 Discussion of the results

## 4.1 Sample complexity of optimal transport divergences

Empirical optimal transport suffers from the curse of dimensionality [27], but it is now well established that its statistical convergence depends on the intrinsic complexity of the target measures [35]. If one of them is discrete, the convergence rate is of actually of order $O(1/\sqrt{n})$ [56, 29, 35, 46, 40]. Theorem 3.1 allows to recover convergence rates for the quantity $\mathbb{E}|\mathbb{D}(\mu_n, \nu) - \mathbb{D}(\mu, \nu)|$, for any $\nu \in \mathcal{M}_1(\mathbb{B}_R)$, as well as for the empirical divergences themselves. These do not improve upon the literature.
## 4.2 Barycenters with sparse support

For a sparse barycenter with at most $N$ support points and sample size $n$, Theorem 3.1 provides generalization bounds of order $\sqrt{N/n}$, classical in semi-discrete optimal transport [20, Theorem 2.10]. Interestingly, this rate is uniform in the number of measures $L$, in the regularization parameter $\epsilon$ for the Sinkhorn divergence, without exponential dependency on the ambient dimension $d$.

For $\mathbb{D} = W_2^2$, assuming that each measure $\mu^\ell$ is supported on $n^\ell$ points, and setting $N \ge \sum_{\ell=1}^L n^\ell - L + 1$, problems (1.1) and (1.2) coincide: an unconstrained barycenter turns out to be sparse [5, Theorem 2]. In this setting, [34] obtains the same rate of convergence $\sqrt{N/\min_{\ell}(n^\ell)}$ for the unconstrained barycenter, this analysis being limited to large $N$.

## 4.3 Tightness of the result

In this section, we discuss the optimality of the estimation error bound in Corollary 3.4 with respect to the sample size $n$ and the sparsity level $N$.

### 4.3.1 Lower bound of the estimation error

We aim to establish a lower bound for the estimation error, i.e., the left-hand side term in Corollary 3.4, in a worst-case sense over the choice of target measures. We rely on [7, Theorem 1] which provides minimax lower bounds for the optimal quantization cost ($L = 1$ and $\mathbb{D} = W_2^2$). A direct extension of this result provides a lower bound for the sparse 2-Wasserstein barycenter problem.

**Lemma 4.1.** Let $\mathbb{D} = W_2^2$, the number of measures $L \in \mathbb{N}^*$ and the dimension $d \in \mathbb{N}^*$. Assume $N \ge 3$ and $n \ge 3435N$ and let $\sum_{i=1}^N \pi_i^n \delta_{y_i^n} \in \mathcal{M}_1^N(\mathbb{R}^d)$ be a random measure depending on $nL$ random variables in $\mathbb{B}_R$. Then there exists $\mu^1, \dots, \mu^L \in \mathcal{M}_1(\mathbb{R}^d)$ such that the estimation error of Corollary 3.4 verifies

$$ \mathbb{E} \left[ F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i^n \delta_{y_i^n} \right) - \min_{(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N} F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right] \ge C_d \sqrt{\frac{N^{1-\frac{4}{d}}}{n}} \quad (4.1) $$

where $C_d$ is an explicit constant which only depends on $d$.

The proof of Lemma 4.1 can be found in D.2. For a fixed $N$, Lemma 4.1 implies that the $O(1/\sqrt{n})$ convergence rate for barycenter functional in 3.4 cannot be improved without further assumptions. However the dependency in $N$ differs.

### 4.3.2 Dependency in $N$ in the estimation error

There is an $O(N^{-2/d})$ factor between the lower bound in (4.1) and the upper bound in (3.1). This raises the question of the optimal exponent in $N$. We exhibit a situation where the lower bound is tight.

Setting $N = n$, $L = 1$ and $\pi = (1/n, \dots, 1/n)$, the upper bound (3.1) of Corollary 3.4 remains constant as $n$ increases. However, we will argue that the error actually decreases like $O(n^{-2/d})$. Indeed, in this setting, the minimizer of $\mu \mapsto W_2^2(\mu_n, \mu)$ over $\mathcal{M}_1^N(\mathbb{R}^d)$, is clearly $\mu_n$ and thus the estimation error writes

$$ \mathbb{E} \left[ W_2^2(\mu, \mu_n) - \min_{\nu \in \mathcal{M}_1^N(\mathbb{B}_R)} W_2^2(\mu, \nu) \right] \le \mathbb{E} [W_2^2(\mu, \mu_n)] = O(n^{-2/d}), $$

see [27]. Hence Corollary 3.4 is not optimal in this particular case. The above argument suggests that the optimal exponent in $N$ is related to estimates of the quantity $W_2(\mu, \mu_n)$, a notoriously difficult problem in statistical optimal transport [15, Section 2]. The curse of dimensionality is in fact hidden in the sparsity level $N$, as remarked in [20, Theorem 2.10].

### 4.3.3 Beyond the $W_2^2$ case

We expect the same conclusions for $\mathbb{D} = SW_p^p$ and max-$SW_p^p$ since the $1/\sqrt{n}$ rate is also known to be optimal for these divergences [60]. The bias of the Sinkhorn divergence $\mathbb{D} = W_{\epsilon,p}^p$ benefits from a more favorable $1/n$ convergence rate [21] suggesting that a sharper rate is attainable for barycenters as well.
## 4.4 K-means and constrained K-means

Choosing $\mathbb{D} = W_2^2$ and $L = 1$, Corollary 3.4 provides error bounds for the celebrated K-means problems. Let $X_1, \dots, X_n$ be $n$ i.i.d samples with distribution $\mu \in \mathcal{M}_1(\mathbb{B}_R)$. It is well known minimizing the K-means functional is actually an optimal quantization problem for the empirical measure $\mu_n$ (see Appendix D.4):

$$ \min_{Y \in (\mathbb{B}_R)^N} \frac{1}{n} \sum_{i=1}^{n} \min_{j=1, \dots, N} \|X_i - y_j\|^2 = \min_{\nu \in \mathcal{M}_1^N(\mathbb{B}_R)} W_2^2(\mu_n, \nu). \quad (4.2) $$

Corollary 3.4 then yields a convergence rate of order $O(\sqrt{N/n})$ for the K-means problem, which yields an improvement by a factor $\sqrt{\log(N)}$ compared to [33, Section 4.4]. The Wasserstein functional formulation is closer to the measure approximation view point. It also allows to consider the constrained K-means problem [47], fixing the mass associated to each centroid with a weight vector $\pi \in \Delta_N$ and minimizing the functional $Y \mapsto W_2^2(\mu_n, \sum_{i=1}^N \pi_i \delta_{y_i})$. For uniform weights $\pi_i = 1/N$ for all $i \in [1, N]$, this is the optimal uniform quantization problem to which Corollary 3.4 apply as well.

# 5 Numerical experiments

We illustrate<fnref n="1" /> Corollary 3.4 for $\mathbb{D} = W_2^2$ and $W_{\epsilon, 2}^2$.

**Dependency in $n$** We consider $L = 3$ target Gaussian distributions in dimension $d = 2$ with random means and covariances. We compute the $W_2^2$ and $W_{\epsilon, 2}^2$ ($\epsilon = 1$) empirical barycenter for $N = 10$, $N = 50$ and for $n$ ranging between 1000 and 10000 using the POT library [26]. Figure 1 represents the estimation error in (3.1) as $n$ increases, matching the rate in Corollary 3.4. The presence of hollows in the curve for $\mathbb{D} = W_{\epsilon, p}^p$ is due to numerical instabilities.

<figure><img src="image_2.png" /><figcaption>Figure 1: The graphs show the evolution of the estimation error (in log-scale) as the sample size $n$ increases. Each experiment is performed 40 times. The mean value of the estimation errors is plotted along with the 95% confidence interval. The theoretical upper bound $n \mapsto 1/\sqrt{n}$ is plotted in red.</figcaption></figure>

**Dependency in $N$** We consider the same setting as above in dimension $d = 6$ for which the theoretical rates predict larger errors as $N$ increases. We fix the sample size $n = 500$ and consider $N$ ranging between 10 and 500. Figure 2 shows the estimation error as a function of $N$ together with the upper bound of Corollary 3.4 and lower bound of Lemma 4.1 for $\mathbb{D} = W_2^2$. This experiment suggests that the dependency in $N$ given by the worst case scenario in Lemma 4.1 for $W_2^2$ matches practical observations on the mild gaussian distributions setting.

<footnote n="1">¹The code to reproduce the experiments is available in Python at https://github.com/LeoPtls.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: The graphs show the estimation error (in log-scale) as the support cardinality $N$ increases, for $n = 500$. Each experiments are performed 10 times. The mean value of the estimation errors is plotted along with the 95% confidence interval. The theoretical bounds (Corollary 3.4 and Lemma 4.1) are plotted in red.</figcaption></figure>

# 6 Conclusion

The expected convergence rate of the empirical sparse optimal transport barycenters is shown to be of order $O(\sqrt{N/n})$ for several types of optimal transport divergences. The rate is min-max optimal for the Wasserstein distance $D = W_2^2$. Investigating optimality for sliced distances $D = SW_p^p$, max-$SW_p^p$ and Sinkhorn divergences $D = W_{\epsilon,p}^p$ is left for future work.

# A Semi-dual formulation of optimal transport

The Kantorovich dual formulation of the Wasserstein distances $W_p$ can be recast as a semi-dual problem using the c-transform. Similarly, the Sinkhorn divergences $W_{\epsilon,p}^p$ admit a semi-dual formulation based on the smoothed c-transform. In the semi-discrete setting (i.e. between a continuous and a discrete measure), these formulations turn an infinite-dimensional program (corresponding to the primal formulations (2.1) and (2.3)) into a finite-dimensional one, which is a key point in our proofs. In this section, we present these semi-dual formulations in the semi-discrete setting, with a particular focus on the case where the discrete measure has zero weights or repeated atoms. Throughout this section, we consider a measure $\mu \in M_1(\mathbb{R}^d)$.

## A.1 Semi-dual formulation outside the generalized diagonal

The generalized diagonal plays a central role throughout our proofs, it denotes the set of point clouds with at least two identical components:

$$D_N = \{Y := (y_1, \dots, y_N) \in (\mathbb{R}^d)^N \mid \exists i \neq j \text{ such that } y_i = y_j\}. \quad (A.1)$$

Let us then define $\nu := \sum_{i=1}^N \pi_i \delta_{y_i}$ a discrete measure supported on $Y := (y_1, \dots, y_N) \in (\mathbb{R}^d)^N \setminus D_N$ and strictly positive weights $\pi \in \text{int}(\Delta_N)$. For a bounded continuous function $\phi \in C_b(\mathbb{R}^d)$, let us also denote $\phi^c \in C_b(\mathbb{R}^d)$ its c-transform, defined in our case for all $x \in \mathbb{R}^d$ as $\phi^c(x) = \min_{i=1,\dots,N}(\|x - y_i\|^p - \phi(y_i))$ [53, Section 1.3]. Then by Kantorovich's duality ([58, Theorem 5.10]) and c-convexity of the optimal potentials we have the following formulation for the p-Wasserstein distance between $\mu$ and $\nu$:

$$W_p^p(\mu, \nu) = \max_{\phi \in \mathbb{R}^N} \int_{\mathbb{R}^d} \min_{i=1,\dots,N} \{\|x - y_i\|^p - \phi_i\} d\mu(x) + \sum_{i=1}^N \pi_i \phi_i. \quad (A.2)$$

In (A.2), with a slight abuse of notation, the dual potential $\phi \in C_b(\mathbb{R}^d)$ is identified with its value on the support of $\nu$, which consists of $N$ distinct points: $\phi_i = \phi(y_i)$, for $i = 1, \dots, N$. For this dual formulation to hold, the measure $\nu$ must be supported outside of the generalized diagonal $D_N$, with strictly positive weights; see the discussion in the next Section A.2.
Entropic optimal transport (2.3) between $\mu$ and $\nu$ also admits both a dual and a semi-dual form (see [30, Proposition 2.1]), which are respectively given by the following optimization problems:

$$ W_{\epsilon,p}^{p}(\mu, \nu) = \max_{u,v \in C(\mathbb{R}^d) \times \mathbb{R}^N} \int_{\mathbb{R}^d} u(x) d\mu(x) + \sum_{i=1}^{N} \pi_i v_i - \epsilon \sum_{i=1}^{N} \pi_i \int_{\mathbb{R}^d} \exp \left( \frac{u(x) + v_i - \|x - y_i\|^p}{\epsilon} \right) d\mu(x) \quad (A.3) $$

and

$$ W_{\epsilon,p}^{p}(\mu, \nu) = \max_{w \in \mathbb{R}^N} -\epsilon \int_{\mathbb{R}^d} \log \left( \sum_{i=1}^{N} \pi_i \exp \left( \frac{w_i - \|x - y_i\|^p}{\epsilon} \right) \right) d\mu(x) + \sum_{i=1}^{N} \pi_i w_i - \epsilon. \quad (A.4) $$

These constructions also necessitate that $Y := (y_1, \dots, y_N) \in (\mathbb{R}^d)^N \setminus D_N$ and $\pi \in \text{int}(\Delta_N)$.

## A.2 The case of support points on the generalized diagonal $D_N$ and/or with null weights

As previously mentioned, it is mandatory that $(Y, \pi) \in (\mathbb{R}^d)^N \setminus D_N \times \text{int}(\Delta_N)$ for a functional of the form $(Y, \pi) \mapsto \mathbb{D}(\mu, \sum_{i=1}^N \pi_i \delta_{y_i})$ to admit the duality formulations described above. Indeed, suppose that $y_i = y_j$, for some $i \neq j$, then the dual formulation of the problem should include the constraint $\phi_i = \phi_j$, as both $\phi_i$ and $\phi_j$ represent the evaluation of a continuous potential in the same argument. In other words, we would have $\phi_i = \phi(y_i) = \phi(y_j) = \phi_j$. While this may seem like a minor technicality that can often be overlooked, it cannot be omitted in our analysis. This is particularly relevant in the case of the sliced Wasserstein distance, where two distinct points in $\mathbb{R}^d$ may project onto the same point on a line. Similarly, it is important to ensure that the weights of the discrete measure remain in the interior of the probability simplex $\text{int}(\Delta_N)$.

This technicality has been handled in [50, Section 5.1] with the following proposition.

**Proposition A.1.** Let $(Y, \pi) \in (\mathbb{R}^d)^N \times \Delta_N$, then there exists $M \in [1, N]$ and $(\tilde{Y}, \tilde{\pi}) \in (\mathbb{R}^d)^M \setminus D_M \times \text{int}(\Delta_M)$ such that for all $\mu \in M_1(\mathbb{R}^d)$ and any $p \ge 1$:

$$ W_p \left( \mu, \sum_{i=1}^{N} \pi_i \delta_{y_i} \right) = W_p \left( \mu, \sum_{i=1}^{M} \tilde{\pi}_i \delta_{\tilde{y}_i} \right) \quad \text{and} \quad W_{\epsilon,p} \left( \mu, \sum_{i=1}^{N} \pi_i \delta_{y_i} \right) = W_{\epsilon,p} \left( \mu, \sum_{i=1}^{M} \tilde{\pi}_i \delta_{\tilde{y}_i} \right). $$

## A.3 Semi-duality formulation for semi-discrete sliced divergences

Starting from the definition of the sliced distances in (2.4), the dual formulation (A.2) can be applied to each one-dimensional distance $W_p(P_{\theta_\sharp\mu}, P_{\theta_\sharp\nu})$, $\theta \in \mathbb{S}^{d-1}$, addressing the ties in the support point of the measures $P_{\theta_\sharp\nu}$. Let $(Y, \pi) \in (\mathbb{R}^d)^N \setminus D_N \times \text{int}(\Delta_N)$, then, for each $\theta \in \mathbb{S}^{d-1}$ there exists an integer $M := M(\theta) \in [1, N]$ such that $\langle y_i - y_j, \theta \rangle \neq 0$ and $\pi_i > 0$ for all $i, j \in [1, M(\theta)]$, $i \neq j$. Up to a reordering of components, we can write $\pi_i$ as $\pi_i^\theta := \pi_i + \pi_j$ whenever $\langle y_i - y_j, \theta \rangle = 0$ for all $i, j \in [1, N]$, $i \neq j$. Following this construction, we obtain a probability measure $\nu^\theta := \sum_{i=1}^{M(\theta)} \pi_i^\theta \delta_{\langle y_i, \theta \rangle}$ with no ties in its support points. Then, the Kantorovich duality formulation (A.2) leads to the following expressions of the sliced distances:

$$ SW_p^p(\mu, \nu) = \int_{\mathbb{S}^{d-1}} \left[ \max_{w \in \mathbb{R}^{M(\theta)}} \int_{\mathbb{R}^d} \min_{i=1, \dots, M(\theta)} \{|\langle x - y_i, \theta \rangle|^p - w_i\} d\mu(x) + \sum_{i=1}^{M(\theta)} \pi_i w_i \right] d\sigma(\theta) \quad (A.5) $$

and

$$ \max - SW_p^p(\mu, \nu) = \max_{\theta \in \mathbb{S}^{d-1}} \max_{w \in \mathbb{R}^{M(\theta)}} \int_{\mathbb{R}^d} \min_{i=1, \dots, M(\theta)} \{|\langle x - y_i, \theta \rangle|^p - w_i\} d\mu(x) + \sum_{i=1}^{M(\theta)} \pi_i w_i, \quad (A.6) $$

where $\sigma$ denotes the uniform measure over the $d$-sphere $\mathbb{S}^{d-1}$. Note that actually, one can replace $M(\theta)$ by $N$ in the maximization over the variable $w$ in (A.5), without taking into account any dependency on $\theta$. Indeed, let $M(\theta) := \text{card}(\{i \in [1, N] : \langle y_i - y_j, \theta \rangle \neq 0, \forall j \neq i\})$. We thus have $M(\theta) < N$ if and only if there exist indices
i ≠ j such that ⟨yᵢ - yⱼ, θ⟩ = 0, meaning that θ lies in a union of sets of the form {θ ∈ ℝ<sup>d-1</sup> : ⟨yᵢ - yⱼ, θ⟩ = 0}, for i, j ∈ [1, N] and i ≠ j. Since Y ∈ (ℝ<sup>d</sup>)<sup>N</sup> \ D<sub>N</sub> this is a finite union of hyperplanes intersected with the sphere, which is σ-negligible. It follows that M(θ) = N for almost every θ ∈ ℝ<sup>d-1</sup>.

B Boundedness of the dual variables

In this section, we consider a measure μ ∈ M₁(ℝ<sub>R</sub>), R > 0 and a discrete measure ν ∈ M<sub>1</sub><sup>N</sup>(ℝ<sub>R</sub>). We prove that the dual variables associated to the optimal transport problem D(μ, ν) for the Wasserstein distances D = W<sub>p</sub><sup>p</sup> and the Sinkhorn divergences D = W<sub>ε,p</sub><sup>p</sup> are bounded.

**Proposition B.1.** Let μ ∈ M₁(ℝ<sub>R</sub>) for some R > 0 and D = W<sub>p</sub><sup>p</sup>, or W<sub>ε,p</sub><sup>p</sup>. Then for any (Y, π) ∈ (ℝ<sub>R</sub><sup>N</sup> \ D<sub>N</sub>) × int(Δ<sub>N</sub>), we have

$$
W_p^p \left( \mu, \sum_{i=1}^{N} \pi_i \delta_{y_i} \right) = \max_{\substack{w \in \mathbb{R}^N \\ |w_i| \le K_1}} \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \{ \|x - y_i\|^p - w_i \} d\mu(x) + \sum_{i=1}^{N} \pi_i w_i,
$$

and

$$
W_{\epsilon,p}^{p} \left( \mu, \sum_{i=1}^{N} \pi_i \delta_{y_i} \right) = \max_{\substack{w \in \mathbb{R}^N \\ |w_i| \le K_2}} -\epsilon \int_{\mathbb{R}^d} \log \left( \sum_{i=1}^{N} \pi_i \exp \left[ \frac{w_i - \|x - y_i\|^p}{\epsilon} \right] \right) d\mu(x) + \sum_{i=1}^{N} \pi_i w_i - \epsilon
$$

where $K_1 = 2pR^p$ and $K_2 = 4p(2R)^p$.

**Remark B.2.** In all generality, the dual variables $w \in \mathbb{R}^N$ are not bounded since adding a constant $c \in \mathbb{R}$ to $w$ yields another valid solution. Therefore, to fix a dual solution, we set one coordinate of $w$ to 0 and under this constraint the dual variable $w$ remains bounded uniformly in $Y$ and $\pi$.

*Proof of Proposition B.1.* In the following, we prove both cases independently.

**Case** $\mathbb{D} = W_p^p$. Let $(Y, \pi) \in (\mathbb{R}_R^N \setminus D_N) \times \text{int}(\Delta_N)$, then by Kantorovich duality (A.2) between $\mu$ and the discrete measure $\sum_{i=1}^N \pi_i \delta_{y_i}$, we have

$$
W_p^p \left( \mu, \sum_{i=1}^{N} \pi_i \delta_{y_i} \right) = \max_{w \in \mathbb{R}^N} \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \{ \|x - y_i\|^p - w_i \} d\mu(x) + \sum_{i=1}^{N} \pi_i w_i. \quad (B.1)
$$

As discussed in Remark (B.2), we may select a choice of dual variables $w \in \mathbb{R}^N$ optimal in (B.1) and verifying $w_N = 0$. From [53, Proposition 5.11] and [53, Remark 1.13] that follows, we may also select a c-convave function $w^c$ which verifies for all $i \in [1, N]: w_i = (w^c(y_i))^c$. We then have

$$
\begin{align*}
|w_i - w_N| &= |(w^c(y_i))^c - (w^c(y_N))^c| \\
&= \left| \min_{x \in \text{Supp}(\mu)} \{\|x - y_i\|^p - w^c(x)\} - \min_{x \in \text{Supp}(\mu)} \{\|x - y_N\|^p - w^c(x)\} \right| \\
&\le \max_{x \in \mathbb{B}_R} \|x - y_i\|^p - \|x - y_N\|^p \tag{B.2}
\end{align*}
$$

$$
\le pR^{p-1} \|z_i - y_i\| \le 2pR^p, \qquad (B.3)
$$

where inequalities (B.2) and (B.3) come from Lemmas D.2 and D.1. Using that $w_N = 0$, we get the result: the dual variables for the *p*-Wasserstein distances are bounded, and $|w_i| \le 2pR^p$.

**Case** $\mathbb{D} = W_{\epsilon,p}^p$. Let $(Y, \pi) \in (\mathbb{R}_R^N \setminus D_N) \times \text{int}(\Delta_N)$. Let also $u \in C_b(\mathbb{R}_R)$ be an optimal solution in the dual (A.3) formulation of the entropic *p*-Wasserstein distance between $\mu$ and $\sum_{i=1}^N \pi_i \delta_{y_i}$. By [29, Proposition 1], $u$ is $\kappa$-Lipshitz, where $\kappa$ is the Lipschitz constant of $(a, b) \in \mathbb{B}_R^2 \mapsto \|a - b\|^p$. Then since $\text{Supp}(\mu) \subset \mathbb{B}_R$, we have by Lemma D.1 that $\kappa$ is upper bounded by $p\sqrt{2}(2R)^{p-1}$. Let $w$ be a solution of (A.4) for which we can suppose
$w_N = 0$ (Remark (B.2)). As in [30], one can show with a first order condition on (A.3) that for all $y$ a coordinate in $Y: w(y) = -\epsilon \log \left( \int_{\mathbb{R}^d} \exp \left( \frac{u(x)-\|x-y\|^p}{\epsilon} \right) d\mu(x) \right)$. We thus have for all $i \neq N$:

$$
\begin{aligned}
|w(y_i) - w(y_N)| &= \left| \epsilon \log \left( \frac{\int_{\mathbb{R}^d} \exp \left( \frac{u(x)-\|x-y_i\|^p}{\epsilon} \right) d\mu(x)}{\int_{\mathbb{R}^d} \exp \left( \frac{u(x)-\|x-y_N\|^p}{\epsilon} \right) d\mu(x)} \right) \right| \\
&\le \left| \epsilon \log \left( \frac{\int_{\mathbb{R}^d} \exp \left( \frac{u(\bar{x}_i)-\|\bar{x}_i-y_i\|^p}{\epsilon} \right) d\mu(x)}{\int_{\mathbb{R}^d} \exp \left( \frac{u(\bar{x}_N)-\|\bar{x}_N-y_N\|^p}{\epsilon} \right) d\mu(x)} \right) \right| \\
&= |u(\bar{x}_i) - \|\bar{x}_i - y_i\|^p - (u(\underline{x}_N) - \|\underline{x}_N - y_N\|^p)| \\
&\le |u(\bar{x}_i) - u(\underline{x}_N)| + \|\bar{x}_i - y_i\|^p - \|\underline{x}_N - y_N\|^p \\
&\le p\sqrt{2}(2R)^{p-1}\|\bar{x}_i - \underline{x}_N\| + \|\bar{x}_i - y_i\|^p + \|\underline{x}_N - y_N\|^p \quad (B.4) \\
&\le p\sqrt{2}(2R)^p + 2(2R)^p \le 4p(2R)^p,
\end{aligned}
$$

where $\bar{x}_i$ and $\underline{x}_N$ denote respectively a maximizer of $x \mapsto \exp\left(\frac{u(x)-\|x-y_i\|^p}{\epsilon}\right)$ and a minimizer of $x \mapsto \exp\left(\frac{u(x)-\|x-y_N\|^p}{\epsilon}\right)$ over $\mathbb{B}_R$. Using that $w(y_N) = 0$ then yields the result. $\square$

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

# C Proofs of the main results

In this section, we prove the sample complexity of optimal transport barycenters with sparse support in Theorem 3.1. For $\ell \in [1, L]$, we denote $\mathbb{P}_n^\ell$ the empirical process $f \mapsto \mathbb{P}_n^\ell f := \frac{1}{n} \sum_{i=1}^n f(X_i^\ell)$ for $n$ i.i.d random variables $X_1^\ell, \dots, X_n^\ell$ of law $\mu^\ell$. For any subset $\mathcal{F}$ of a normed space $(\Theta, \|\cdot\|)$ and $\tau > 0$, we denote $\mathcal{N}(\tau, \mathcal{F}, \|\cdot\|)$ the $\tau$-covering number of $\mathcal{F}$ with respect to $\|\cdot\|$, which is the minimum number of balls of radius $\tau$, with respect to the norm $\|\cdot\|$, needed to cover $\mathcal{F}$. We refer the reader to [54, Definition 2.2] for more details on the covering number. Finally, we also denote $\|\cdot\|_{\mathbb{L}^2(\mu)}$ the $\mathbb{L}^2$ norm with respect to the measure $\mu$ and $\bar{B}_d(0, R)$ the $d$-dimensional closed ball of radius $R$. The proofs for Wasserstein (Section C.1), Sliced Wasserstein (Section C.2) and Sinkhorn (Section C.3) barycenters all follow the same approach : we first derive the covering number for a class of functions associated to the dual formulation of the divergences $\mathcal{D}$, and then apply tools from empirical process theory to obtain a generalization error bound for the barycenters.

## C.1 Wasserstein barycenter

**Lemma C.1.** Let

$$ \mathcal{F}_p = \left\{ \left\{ f_Y^w : x \in \mathbb{B}_R \mapsto \min_{i=1,\dots,N} \|x - y_i\|^p - w_i \right\}, Y \in \mathbb{B}_R^N, w \in \bar{B}_1(0, 2pR^p)^N \right\}, \quad (C.1) $$

then for all $\tau > 0$ and $\mu \in \mathcal{M}_1(\mathbb{B}_R)$, $\mathcal{N}(\tau, \mathcal{F}_p, \|\cdot\|_{\mathbb{L}^2(\mu)}) \le \left(1 + \frac{4pR^p}{\tau}\right)^{dN} \left(1 + \frac{8pR^p}{\tau}\right)^N$. Also, the subset $\mathcal{F}_{p,\mathcal{Q}}$ of functions of $\mathcal{F}_p$ with rationals parameters $(Y, w)$, is dense in $\mathcal{F}_p$.

**Proof.** Let $(Y, w)$, $(Z, v) \in \mathbb{B}_R^N \times [-2pR^p, 2pR^p]^N$. Then for all $x \in \mathbb{B}_R$

$$
\begin{aligned}
|f_Y^w(x) - f_{Z,v}(x)| &= \left| \min_{i=1,\dots,N} \|x - y_i\|^p - w_i - \min_{i=1,\dots,N} \|x - z_i\|^p - v_i \right| \\
&\le \max_{i=1,\dots,N} \|x - y_i\|^p - w_i - (\|x - z_i\|^p - v_i) \quad (C.2)
\end{aligned}
$$

$$ \le p R^{p-1} \max_{i=1,\dots,N} \|y_i - z_i\| + \max_{i=1,\dots,N} |w_i - v_i|, \quad (C.3) $$
where inequalities (C.2) and (C.3) are due to Lemmas D.2 and D.1. Let us denote $T = pR^{p-1}$, and let $\tau > 0$. We then consider $z^1, \dots, z^{n(\tau)}$ a $\frac{\tau}{2T}$-covering of $\mathbb{B}_R$ and $v^1, \dots, v^{m(\tau)}$ a $\frac{\tau}{2}$-covering of $\mathbb{B}_1(0, 2pR^p)$. For any vector of indices $i = (i_1, \dots, i_N) \in [1, n(\tau)]^N$ and $k = (k_1, \dots, k_N) \in [1, m(\tau)]^N$, we define for all $x \in \mathbb{B}_R$ the function $f_{i,k}(x) = \min_{j=1,\dots,N} \|x - z_{ij}^\ell\|^p - v^{k_j}$. Let $f := f_Y^w \in \mathcal{F}_p$, then by (C.3) there exists $i, k$ such that

$$|f(x) - f_{i,k}(x)| \le pR^{p-1} \max_{j=1,\dots,N} \|y_j - z_{ij}^\ell\| + \max_{j=1,\dots,N} |w_j - v^{k_j}| \le T \frac{\tau}{2T} + \frac{\tau}{2} = \tau,$$

so that $\|f - f_{i,k}\|_{L^2(\mu)} \le \tau$. Since there exists $n(\tau)^N \times m(\tau)^N$ possibilities for such functions $f_{i,k}$, it follows from the covering construction that:

$$\begin{aligned} \mathcal{N}(\tau, \mathcal{F}_p, \|\cdot\|_{L^2(\mu)}) &\le \mathcal{N}\left(\frac{\tau}{2T}, \overline{\mathbf{B}_d(0, R)}, \|\cdot\|\right)^N \mathcal{N}\left(\frac{\tau}{2}, \overline{\mathbf{B}_1(0, 2pR^p)}, |\cdot|\right)^N \\ &\le \left(1 + \frac{4RT}{\tau}\right)^{dN} \left(1 + \frac{8pR^p}{\tau}\right)^N = \left(1 + \frac{4pR^p}{\tau}\right)^{dN} \left(1 + \frac{8pR^p}{\tau}\right)^N \end{aligned}$$

where the second inequality is due to [41, Lemma 4.14].

Finally, $\mathcal{F}_{p,\mathbb{Q}}$ is dense in $\mathcal{F}_p$ by a continuity argument and inequality (C.3). $\square$

*Proof of Theorem 3.1 for $\mathbb{D} = W_p^p$. Let $(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N$ and let $M \in [1, N]$ and $(\tilde{Y}, \tilde{\pi}) \in (\mathbb{B}_R^M \setminus D_M) \times \text{int}(\Delta_M)$ be as in Proposition A.1. For all $\ell \in [1, L]$, let $X_1^\ell, \dots, X_n^\ell$ be $n$ i.i.d random variables of law $\mu^\ell$. Using Proposition B.1 and denoting $K := 2pR^p$ we therefore have*

$$\begin{aligned} \text{Gen}_n(Y, \pi) &:= \left| \mathbb{F}_{\mathbb{D}} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) - \mathbb{F}_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right| \\ &\le \frac{1}{L} \sum_{\ell=1}^L \max_{|\tilde{w}_i| \le K} \left| \frac{1}{n} \sum_{k=1}^n \min_{i=1, \dots, M} \{ \|X_k^\ell - \tilde{y}_i\|^p - \tilde{w}_i \} - \int_{\mathbb{R}^d} \min_{i=1, \dots, M} \{ \|x - \tilde{y}_i\|^p - \tilde{w}_i \} d\mu^\ell(x) \right| \\ &= \frac{1}{L} \sum_{\ell=1}^L \max_{\substack{w_i = w_j \\ w_i \le K \text{ if } y_i = y_j}} \left| \frac{1}{n} \sum_{k=1}^n \min_{i=1, \dots, N} \{ \|X_k^\ell - y_i\|^p - w_i \} - \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \{ \|x - y_i\|^p - w_i \} d\mu^\ell(x) \right|. \end{aligned}$$

where the inequality is due to Lemma D.2. The last equality is due to the fact that for all $Y \in \mathbb{B}_R^N$ and $x \in \mathbb{R}^d$, we have $\min_{k=1, \dots, M} \|x - y_k\|^p - w_k = \min_{k=1, \dots, N} \|x - y_k\|^p - w_k$ if $w_i = w_j$ whenever $y_i = y_j$. Notice that the dependency in the vector of weights $\pi$ has disappeared. For any $(x, Y, w) \in \mathbb{B}_R \times \mathbb{B}_R^N \times [-K, K]^N$, let $f_Y^w(x) = \min_{i=1, \dots, N} \|x - y_i\|^p - w_i$. We then have

$$\begin{aligned} \sup_{(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N} \text{Gen}_n(Y, \pi) &\le \frac{1}{L} \sum_{\ell=1}^L \sup_{(Y, w) \in \mathbb{B}_R^N \times [-K, K]^N} \left| \mathbb{P}_n^\ell (f_Y^w) - \mathbb{E}_{X \sim \mu^\ell} [f_Y^w(X)] \right| \\ &= \frac{1}{L} \sum_{\ell=1}^L \sup_{f \in \mathcal{F}_p} \left| \mathbb{P}_n^\ell (f) - \mathbb{E}_{X \sim \mu^\ell} [f(X)] \right| \end{aligned}$$

where $\mathcal{F}_p$ is defined in (C.1). Notice further that for any $\ell \in [1, L]$,

$$\begin{aligned} \sup_{f \in \mathcal{F}_p} \mathbb{P}_n^\ell f^2 &= \sup_{(Y,w) \in \mathbb{B}_R^N \times [-K,K]^N} \frac{1}{n} \sum_{k=1}^n \left( \min_{i=1, \dots, N} \|X_k^\ell - y_i\|^p - w_i \right)^2 \\ &\le \frac{1}{n} \sum_{k=1}^n \left( \sup_{x,y \in \mathbb{B}_R} \|x - y\|^p + \sup_{w \in [-K,K]} |w| \right)^2 \\ &\le ((2R)^p + 2pR^p)^2 \le (2(2R)^p)^2, \end{aligned}$$
which gives $\sqrt{\sup_{f \in \mathcal{F}} \mathbb{P}_n f^2} \le 2(2R)^p$. By Lemma C.1, $\mathcal{F}_p$ is dense in $\mathcal{F}_{p, \mathbb{Q}}$, and we can then apply [32, Theorem 3.5.1] on empirical processes and covering number. Therefore for all $\ell$:

$$ \mathbb{E} \left[ \sqrt{n} \sup_{f \in \mathcal{F}_p} |\mathbb{P}_n^\ell(f) - \mathbb{E}_{X \sim \mu^\ell}[f(X)]| \right] \le 8\sqrt{2}\mathbb{E} \left[ \int_0^{2(2R)^p} \sqrt{\log(2\mathcal{N}(\tau/2, \mathcal{F}_p, \|\cdot\|_{\mathbb{L}^2(\mu_n^\ell)})} d\tau \right]. $$

From Lemma C.1 applied to the measure $\mu_n^\ell$, denoting $\tilde{K} = 2(2R)^p$, we get

$$ \begin{aligned} \int_0^{\tilde{K}} \sqrt{\log(2\mathcal{N}(\tau/2, \mathcal{F}_p, \|\cdot\|_{\mathbb{L}^2(\mu_n^\ell)}))} d\tau &\le \int_0^{\tilde{K}} \sqrt{\log \left( 2 \left[ 1 + \frac{8pR^p}{\tau} \right]^{dN} \left[ 1 + \frac{16pR^p}{\tau} \right]^N \right)} d\tau \\ &\le \sqrt{N(d+1)} \int_0^{\tilde{K}} \sqrt{\log \left( 2 + \frac{32pR^p}{\tau} \right)} d\tau. \end{aligned} $$

Finally, as this last term does not depend on $\ell$, we get

$$ \mathbb{E} \left[ \frac{1}{L} \sum_{\ell=1}^{L} \sup_{f \in \mathcal{F}_p} |\mathbb{P}_n^\ell(f) - \mathbb{E}_{X \sim \mu^\ell}[f(X)]| \right] \le C_{p,R} \sqrt{\frac{N(d+1)}{n}}, $$

$$ \text{where } C_{p,R} = 8\sqrt{2} \int_0^{2(2R)^p} \sqrt{\log\left(2 + \frac{32pR^p}{\tau}\right)} d\tau. \quad \square $$

## C.2 Sliced Wasserstein barycenter

**Lemma C.2.** Let

$$ \mathcal{F}_p^s = \left\{ \left\{ f_{\theta,Y,w} : x \mapsto \min_{i=1,\dots,N} \left| \langle x - y_i, \theta \rangle^p - w_i \right| \right\} , Y \in \mathbb{B}_R^N, w \in \overline{B_1(0, 2pR^p)}^N, \theta \in \mathbb{S}^{d-1} \right\}, $$

then for all $\tau > 0$, $\mu \in \mathcal{M}_1(\mathbb{B}_R)$, $\mathcal{N}(\tau, \mathcal{F}_p^s, \|\cdot\|_{\mathbb{L}^2(\mu)}) \le \left(1 + \frac{12pR^p}{\tau}\right)^d \left(1 + \frac{6pR^p}{\tau}\right)^{dN} \left(\frac{12pR^p}{\tau}\right)^N$. Also, the subset $\mathcal{F}_{p, \mathbb{Q}}^s$ of functions of $\mathcal{F}_p^s$ with rationals parameters $(Y, w, \theta)$, is dense in $\mathcal{F}_p^s$.

**Proof.** Let $K = 2pR^p$, $Y, Z \in \mathbb{B}_R^N$, $w, v \in [-K, K]^N$ and $\theta, \varphi \in \mathbb{S}^{d-1}$. Then by Lemma (D.1), we have for $x \in \mathbb{B}_R$

$$ \begin{aligned} & | |\langle x - y_i, \theta \rangle^p - |\langle x - z_i, \varphi \rangle^p| \le pR^{p-1} |\langle x - y_i, \theta \rangle - \langle x - z_i, \varphi \rangle| \\ & = pR^{p-1} |\langle x, \theta - \varphi \rangle + \langle z_i, \varphi - \theta \rangle + \langle z_i, \theta \rangle - \langle y_i, \theta \rangle| \\ & \le pR^{p-1} (2R \|\theta - \varphi\| + \|z_i - y_i\|). \end{aligned} $$

Therefore, thanks to Lemma (D.2), we get

$$ \begin{aligned} |f_{\theta,Y,w}(x) - f_{\varphi,Z,v}(x)| &\le \max_{i=1,\dots,N} | |\langle x - y_i, \theta \rangle^p - |\langle x - z_i, \varphi \rangle^p| + |w_i - v_i| \\ &\le \max_{i=1,\dots,N} pR^{p-1} (2R \|\theta - \varphi\| + \|z_i - y_i\|) + |w_i - v_i| \\ &\le 2pR^p \|\theta - \varphi\| + pR^{p-1} \max_{i=1,\dots,N} \|z_i - y_i\| + \max_{i=1,\dots,N} |w_i - v_i|. \end{aligned} \quad (C.4) $$

Then, using the same argument as in the proof of Lemma C.1, we can upper bound the covering number $\mathcal{N}(\tau, \mathcal{F}_p^s, \|\cdot\|_{\mathbb{L}^2(\mu)})$ by the following product of covering numbers:

$$ \begin{aligned} \mathcal{N}(\tau, \mathcal{F}_p^s, \|\cdot\|_{\mathbb{L}^2(\mu)}) &\le \mathcal{N} \left[ \frac{\tau}{6pR^p}, \mathbb{S}^{d-1}, \|\cdot\| \right] \mathcal{N} \left[ \frac{\tau}{3pR^{p-1}}, \mathbb{B}_R, \|\cdot\| \right]^N \mathcal{N} \left[ \frac{\tau}{3}, \overline{B_1(0, 2pR^p)}, |\cdot| \right]^N \\ &\le \left(1 + \frac{12pR^p}{\tau}\right)^d \left(1 + \frac{6pR^p}{\tau}\right)^{dN} \left(\frac{12pR^p}{\tau}\right)^N, \end{aligned} $$

which concludes the proof of the upper bound of the covering number. Finally, $\mathcal{F}_{p, \mathbb{Q}}^s$ is dense in $\mathcal{F}_p^s$ due to inequality (C.4).
Proof of Theorem 3.1 for $\mathbb{D} = SW_p^p$ and $\max-SW_p^p$. For brevity, we only prove the result for $\mathbb{D} = SW_p^p$ as the max-$SW_p^p$ case is similar. Let $M \in [1, N]$ and $(\tilde{Y}, \tilde{\pi}) \in (\mathbb{B}_R^M \setminus D_M) \times \text{int}(\Delta_N)$ be as in Proposition A.1. For all $\ell \in [1, L]$, let $X_1^\ell, \dots, X_n^\ell$ be i.i.d random variables of law $\mu^\ell$. Following the construction in Section A.3, we define the function $M: \theta \mapsto \text{Card}(\{i \in [1, M], \langle y_i - y_j, \theta \rangle \neq 0, \forall j \neq i\})$. By Proposition B.1 and for $K := 2pR^p$, we get

$$
\begin{align*}
\text{Gen}_n^s(Y, \pi) &:= \left| \mathbb{F}_{\mathbb{D}} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \tilde{\pi}_i \delta_{\tilde{y}_i} \right) - \mathbb{F}_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \tilde{\pi}_i \delta_{\tilde{y}_i} \right) \right| \\
&\le \frac{1}{L} \sum_{\ell=1}^L \max_{\theta \in \mathbb{S}^{d-1}} \max_{w \in \mathbb{R}^{M(\theta)} \atop |w_i| \le K} \left| \frac{1}{n} \sum_{k=1}^n \min_{i=1, \dots, M(\theta)} \left\{ \left| \langle X_k^\ell - y_i, \theta \rangle \right|^p - w_i \right\} \right| \\
&\quad - \int_{\mathbb{R}^d} \min_{i=1, \dots, M(\theta)} \left| \langle x - y_i, \theta \rangle|^p - w_i \right| d\mu^\ell(x) \\
&= \frac{1}{L} \sum_{\ell=1}^L \max_{\theta \in \mathbb{S}^{d-1}} \max_{w \in \mathbb{R}^N \atop |w_i| \le K} \left| \frac{1}{n} \sum_{k=1}^n \min_{i=1, \dots, N} \left\{ \left| \langle X_k^\ell - y_i, \theta \rangle \right|^p - w_i \right\} \right| \\
&\quad - \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \left| \langle x - y_i, \theta \rangle|^p - w_i \right| d\mu^\ell(x),
\end{align*}
$$

where the first inequality uses Lemma D.2 for the maximum with respect to $w$, and then upper bound the integral with respect to $\theta$ by its maximum value. The last equality is due to the fact that for all $\theta \in \mathbb{S}^{d-1}$, $Y \in (\mathbb{R}^d)^N$ and $x \in \mathbb{R}^d$, we have that $\min_{k=1,...,M(\theta)} |\langle x - y_k, \theta \rangle|^p - w_k = \min_{k=1,...,N} |\langle x - y_k, \theta \rangle|^p - w_k$ whenever $w_i = w_j$ if $\langle y_i - y_j, \theta \rangle = 0$.

For all $(x, \theta, Y, w) \in \mathbb{B}_R \times \mathbb{S}^{d-1} \times \mathbb{B}_R^N \times [-K, K]^N$, let $f_{\theta,Y,w}(x) = \min_{i=1,...,N} |\langle x - y_i, \theta \rangle|^p - w_i$. We then have

$$
\begin{align*}
\sup_{(Y,\pi) \in \mathbb{B}_R^N \times \Delta_N} \text{Gen}_n^s(Y, \pi) &\le \frac{1}{L} \sum_{\ell=1}^L \sup_{(\theta,Y,w) \in \mathbb{S}^{d-1} \times \mathbb{B}_R^N \times [-K,K]^N} |\mathbb{P}_n^\ell(f_{\theta,Y,w}) - \mathbb{E}_{X \sim \mu^\ell}[f_{\theta,Y,w}(X)]| \\
&= \frac{1}{L} \sum_{\ell=1}^L \sup_{f \in \mathcal{F}_p^s} |\mathbb{P}_n^\ell(f) - \mathbb{E}_{X \sim \mu^\ell}[f(X)]|,
\end{align*}
$$

where $\mathcal{F}_p^s$ is defined in Lemma C.2. Moreover, we have for any $\ell \in [1, L]$,

$$
\begin{align*}
\sup_{f \in \mathcal{F}_p^s} \mathbb{P}_n^\ell f^2 &= \sup_{(\theta, Y, w) \in \mathbb{S}^{d-1} \times \mathbb{B}_R^N \times [-K, K]^N} \frac{1}{n} \sum_{k=1}^n \left( \min_{i=1, \dots, N} |\langle X_k^\ell - y_i, \theta \rangle|^p - w_i \right)^2 \\
&\le \frac{1}{n} \sum_{k=1}^n \left( \sup_{x, y \in \mathbb{B}_R, \theta \in \mathbb{S}^{d-1}} \|x - y\|^p \|\theta\|^p + \sup_{w \in [-K, K]} w \right)^2 \\
&\le ((2R)^p + 2pR^p)^2 \le (2(2R)^p)^2.
\end{align*}
$$

With the same argument as in the proof of case $\mathbb{D} = W_p^p$, this quantity can be upper bounded in expectation using [32, Theorem 3.5.1]:

$$
\mathbb{E} \left[ \sqrt{n} \sup_{f \in \mathcal{F}_p^s} | \mathbb{P}_n^\ell(f) - \mathbb{E}_{X \sim \mu^\ell}[f(X)] | \right] \le 8\sqrt{2} \mathbb{E} \left[ \int_0^{2(2R)^p} \sqrt{\log (2\mathcal{N}(\tau/2, \mathcal{F}_p^s, \| \cdot \|_{\mathbb{L}^2(\mu_n^\ell)}) )} d\tau \right].
$$
From Lemma C.2 we thus have

$$
\begin{align*}
& \int_0^{2(2R)^p} \sqrt{\log \left( 2\mathcal{N}(\tau/2, \mathcal{F}_p^s, \|\cdot\|_{\mathbb{L}^2(\mu_n^\ell)}) \right)} d\tau \\
& \le \int_0^{2(2R)^p} \sqrt{\log \left( 2 \left( 1 + \frac{24pR^p}{\tau} \right)^d \left( 1 + \frac{12pR^p}{\tau} \right)^{dN} \left( \frac{24pR^p}{\tau} \right)^N \right)} d\tau \\
& \le \sqrt{N(d+1)+d} \int_0^{2(2R)^p} \sqrt{\log \left( 2 + \frac{48pR^p}{\tau} \right)} d\tau.
\end{align*}
$$

Finally, as this last term does not depend on $\ell$, we get

$$
\mathbb{E} \left[ \frac{1}{L} \sum_{\ell=1}^{L} \sup_{f \in \mathcal{F}_p^s} \left| \mathbb{P}_n^\ell(f) - \mathbb{E}_{X \sim \mu^\ell} [f(X)] \right| \right] \le C_{p,R} \sqrt{\frac{N(d+1)+d}{n}},
$$

$$
\text{where } C_{p,R} = \int_{0}^{2(2R)^p} \sqrt{\log \left( 2 + \frac{48pR^p}{\tau} \right)} d\tau.
$$

□

C.3 Entropy regularized barycenter

**Lemma C.3.** Let

$$
\mathcal{F}_p^\epsilon = \left\{ \left\{ g_{Y,\pi,w}^\epsilon : x \in \mathbb{B}_R \mapsto -\epsilon \log \left( \sum_{i=1}^N \pi_i \exp \left( \frac{w_i - \|x - y_i\|^p}{\epsilon} \right) \right) \right\} \middle| \begin{array}{l} Y \in \mathbb{B}_R^N, \pi \in \Delta_N, w \in \bar{B}_1(0, 4p(2R)^p)^N, \end{array} \right\}
$$

then for all $\tau > 0$ and $\mu \in \mathcal{M}_1(\mathbb{B}_R)$, we have $\mathcal{N}(\tau, \mathcal{F}_p^\epsilon, \|\cdot\|_{\mathbb{L}^2(\mu)}) \le (1 + \frac{4pR^p}{\tau})^{dN} (\frac{16p(2R)^p}{\tau})^N$. Also, the subset $\mathcal{F}_{p,\mathbb{Q}}^\epsilon$ of functions of $\mathcal{F}_p^\epsilon$ with rational parameters $(Y, \pi, w)$, is dense in $\mathcal{F}_p^\epsilon$.

*Proof.* Let $(Y, \pi, w)$, $(Z, \gamma, v) \in \mathbb{B}_R^N \times \Delta_N \times [-4p(2R)^p, 4p(2R)^p]^N$. Then for all $x \in \mathbb{B}_R$, we have

$$
\begin{align*}
|g_{Y,\pi,w}^{\epsilon}(x) - g_{Z,\gamma,v}^{\epsilon}(x)| &= \epsilon \left| \log \left( \frac{\sum_{i=1}^{N} \pi_i \exp\left(\frac{w_i - \|x-y_i\|^p}{\epsilon}\right)}{\sum_{i=1}^{N} \gamma_i \exp\left(\frac{v_i - \|x-z_i\|^p}{\epsilon}\right)} \right) \right| \\
&\le \epsilon \left| \log \left( \frac{\max_{i=1,\dots,N} \exp\left(\frac{w_i - \|x-y_i\|^p}{\epsilon}\right) \sum_{i=1}^{N} \pi_i}{\min_{i=1,\dots,N} \exp\left(\frac{v_i - \|x-z_i\|^p}{\epsilon}\right) \sum_{i=1}^{N} \gamma_i} \right) \right| \\
&= \left| \max_{i=1,\dots,N} w_i - \|x-y_i\|^p - \min_{i=1,\dots,N} v_i - \|x-z_i\|^p \right|.
\end{align*}
$$

Denoting $\bar{i}$ (resp. $i$) is an index for which the maximum (resp. minimum) of $i \mapsto w_i - \|x - y_i\|^p$ (resp. $i \mapsto v_i - \|x - z_i\|^p$) is attained, we get by Lemma D.1,

$$
|g_{Y,\pi,w}(x) - g_{Z,\gamma,v}(x)| \le |w_i - v_i| + pR^{p-1} \|y_i - z_i\|. \quad (C.5)
$$

Notice that the right-hand side of the above inequality depends on neither $\pi$ nor $\epsilon$. Let $T = pR^{p-1}$, with the same proof argument as of case $\mathbb{D} = W_p^p$ and by boundedness of the entropic dual variables (Proposition B.1) we thus get

$$
\begin{align*}
\mathcal{N}(\tau, \mathcal{F}_p^\epsilon, \|\cdot\|_{\mathbb{L}^2(\mu)}) &\le \mathcal{N}\left(\frac{\tau}{2T}, \overline{\mathbb{B}}_d(0, R), \|\cdot\|\right)^N \mathcal{N}\left(\frac{\tau}{2}, \overline{\mathbb{B}}_1(0, 4p(2R)^p), |\cdot|\right)^N \\
&\le \left(1 + \frac{4pR^p}{\tau}\right)^{dN} \left(\frac{16p(2R)^p}{\tau}\right)^N
\end{align*}
$$
which concludes the proof of the upper bound of the covering number. Finally, $\mathcal{F}_{p,\mathbb{Q}}^{\epsilon}$ is dense in $\mathcal{F}_{p}^{\epsilon}$ due to inequality (C.5). □

**Proof of Theorem 3.1 for $\mathbb{D} = W_{\epsilon,p}^p$.** We fix $\epsilon > 0$. Let $(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N$ and let $M \in [1, N]$ and $(\tilde{Y}, \tilde{\pi}) \in (\mathbb{B}_R^M \setminus D_M) \times \text{int}(\Delta_M)$ be as in Proposition A.1. For any triplet $(Y, \pi, w) \in \mathbb{B}_R^N \times \Delta_N \times \mathbb{R}^N$ and $x \in \mathbb{B}_R$, we let $g_{Y,\pi,w}^{\epsilon}(x) := -\epsilon \log \left(\sum_{j=1}^N \pi_j \exp \left(\frac{w_j - \|x-y_j\|^p}{\epsilon}\right)\right)$. For all $\ell \in [1, L]$, let $X_1^{\ell}, \dots, X_n^{\ell}$ be i.i.d random variables of law $\mu^{\ell}$. Using Proposition B.1 and the same argument as in the case $\mathbb{D} = W_p^p$, we get for $K := 4p(2R)^p$:

$$ \sup_{(Y,\pi) \in \mathbb{B}_R^N \times \Delta_N} \left| \mathbb{F}_{\mathbb{D}} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) - \mathbb{F}_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right| \\ \le \frac{1}{L} \sum_{\ell=1}^L \sup_{(Y,\pi,w) \in \mathbb{B}_R^N \times \Delta_N \times [-K,K]^N} \left| \mathbb{P}_n^{\ell} g_{Y,\pi,w}^{\epsilon} - \mathbb{E}_{X \sim \mu^{\ell}} [g_{Y,\pi,w}^{\epsilon}(X)] \right| \\ = \frac{1}{L} \sum_{\ell=1}^L \sup_{g \in \mathcal{F}_p^{\epsilon}} \left| \mathbb{P}_n^{\ell}(g) - \mathbb{E}_{X \sim \mu^{\ell}} [g(X)] \right| $$

where $\mathcal{F}_p^{\epsilon}$ is defined in Lemma C.3. Moreover, for all $w \in \overline{\mathbb{B}}_1(0, 4p(2R)^p)^N$, we have

$$ \begin{aligned} -\frac{4p(2R)^p + (2R)^p}{\epsilon} &\le \frac{w_j - \|x - y_j\|^p}{\epsilon} \le \frac{4p(2R)^p}{\epsilon} \\ \iff -\frac{(4p+1)(2R)^p}{\epsilon} &\le \log \left( \sum_{j=1}^N \pi_j \exp \left( \frac{w_j - \|x - y_j\|^p}{\epsilon} \right) \right) \le \frac{4p(2R)^p}{\epsilon} \\ \iff -4p(2R)^p &\le g_{Y,\pi,w}^{\epsilon}(x) \le (4p+1)(2R)^p. \end{aligned} \quad (C.6) $$

Then, by Lemma C.3, the bound obtained in (C.6) and, once again, with the same argument as in the case $\mathbb{D} = W_p^p$, we obtain by [32, Theorem 3.5.1] the following bound

$$ \begin{aligned} \mathbb{E} \sqrt{n} \sup_{g \in \mathcal{F}_p^{\epsilon}} |\mathbb{P}_n^{\ell}(g) - \mathbb{E}_{X \sim \mu^{\ell}} [g(X)]| &\le 8\sqrt{2} \mathbb{E} \left[ \int_0^{(4p+1)(2R)^p} \sqrt{\log \left( 2N(\tau/2, \mathcal{F}_p^{\epsilon}, \| \cdot \|_{\mathbb{L}^2(\mu_n^{\ell})}) \right)} d\tau \right] \\ &\le 8\sqrt{2} \int_0^{(4p+1)(2R)^p} \sqrt{\log \left( 2 \left( 1 + \frac{8pR^p}{\tau} \right)^{dN} \left( \frac{32p(2R)^p}{\tau} \right)^N \right)} d\tau \\ &\le 8\sqrt{2} \sqrt{N(d+1)} \int_0^{(4p+1)(2R)^p} \sqrt{\log \left( 2 + \frac{64p(2R)^p}{\tau} \right)} d\tau \end{aligned} $$

so that finally

$$ \mathbb{E} \left[ \frac{1}{L} \sum_{\ell=1}^{L} \sup_{g \in \mathcal{F}_p^{\epsilon}} |\mathbb{P}_n^{\ell}(g) - \mathbb{E}_{X \sim \mu^{\ell}} [g(X)]| \right] \le C_{p,R} \sqrt{\frac{N(d+1)}{n}}, $$

$$ \text{where } C_{p,R} = 8\sqrt{2} \int_{0}^{(4p+1)(2R)^p} \sqrt{\log \left( 2 + \frac{64p(2R)^p}{\tau} \right)} d\tau. \quad \square $$

In the following, we prove the result for debiased Sinkhorn barycenters.

**Proof of Corollary 3.6.** Let $(Y^*, \pi^*)$ and $(Y^n, \pi^n)$ be respective minimizers of $(Y, \pi) \mapsto \mathbb{F}_{\mathbb{D}}(\mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i})$ and $(Y, \pi) \mapsto \mathbb{F}_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{y_i})$, which exist by continuity on the compact set A. We denote $\nu^* =$
$\sum_{i=1}^{N} \pi_i^* \delta_{y_i^*}$ and $\nu^n = \sum_{i=1}^{N} \pi_i^n \delta_{y_i^n}$. Then

$$
\begin{align*}
& F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu^n) - F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu^*) \\
& \leq F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu^n) - F_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \nu^n) + F_{\mathbb{D}}(\mu_n^1, \dots, \mu_n^L, \nu^*) - F_{\mathbb{D}}(\mu^1, \dots, \mu^L, \nu^*) \\
& \leq 2 \sup_{(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N} \left| F_{\mathbb{D}} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{Y_i} \right) - F_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{Y_i} \right) \right| \\
& \leq \sup_{(Y, \pi) \in \mathbb{B}_R \times \Delta_N} 2 \left| F_{W_{\epsilon,p}^p} \left( \mu_n^1, \dots, \mu_n^L, \sum_{i=1}^N \pi_i \delta_{Y_i} \right) - F_{W_{\epsilon,p}^p} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{Y_i} \right) \right| \\
& \qquad + \frac{1}{L} \sum_{\ell=1}^L \left| W_{\epsilon,p}^p (\mu_n^\ell, \mu_n^\ell) - W_{\epsilon,p}^p (\mu^\ell, \mu^\ell) \right|.
\end{align*}
$$

In the right-hand side of the inequality, we bound the first term by Theorem 3.1 while the second term can be bounded thanks to the proof of [29, Theorem 3]. In fact, their proof can be adapted almost directly to show that $E|W_{\epsilon,p}^p(\mu_n, \mu_n) - W_{\epsilon,p}^p(\mu, \mu)| = O\left(\frac{\epsilon^{\frac{K}{\epsilon}}}{\sqrt{n}}\left(1 + \frac{1}{\epsilon^{[d/2]}}\right)\right)$. This result is given for arbitrary measures $\alpha, \beta \in \mathcal{M}_1(\mathbb{R}^d)$, their empirical counterparts $\hat{\alpha}_n, \hat{\beta}_n$ and the quantity $W_{\epsilon,p}^p(\hat{\alpha}_n, \hat{\beta}_n)$, but at no point in the proof the independence of the samples upon which $\hat{\alpha}_n$ and $\hat{\beta}_n$ are constructed is mandatory. As a result, one can obtain the same bound with $\beta = \alpha = \mu$. $\square$

# D Additional Lemmas

This section includes elementary, yet useful Lemmas for the proof of our results. Additionally, for completeness, we prove a well known fact about optimal quantization that is used in several instances of our article.

## D.1 Technical lemmas

**Lemma D.1.** For $R > 0$ and $p \ge 1$, we have

(i) the function $f: x \in \mathbb{B}_R \mapsto \|x\|^p$ is $pR^{p-1}$ Lipschitz,

(ii) the function $c: (x, y) \in \mathbb{B}_R \times \mathbb{B}_R \mapsto \|x - y\|^p$ is $\sqrt{2p(2R)^{p-1}}$ Lipschitz.

*Proof.* (i) The case $p = 1$ is just the result of triangular inequality. Second, by mean value theorem and the chain rule, one can show that for $p > 1$,

$$ \nabla f(x) = p \|x\|^{p-2} x $$

whose norm is upper bounded by $pR^{p-1}$.

(ii) Let $x, y \in \mathbb{B}_R$. We first consider the case $p = 1$:

$$
\begin{align*}
\left| \|x - y\| - \|x' - y'\| \right|^2 &= \left| \|x - y\| - \|x - y'\| - (\|x' - y'\| - \|x - y'\|) \right|^2 \\
&\leq \left| \|x - y\| - \|x - y'\| \right|^2 + \left| \|x - y'\| - \|x' - y'\| \right|^2 \\
&\leq \|x - x'\|^2 + \|y - y'\|^2
\end{align*}
$$

by triangular inequality.

Fore the case $p > 1$, we differentiate the functional $c: (x, y) \mapsto \|x - y\|^p$, that is

$$
\begin{align*}
\partial_x c(x, y) &= \frac{p}{2} (\|x - y\|^2)^{\frac{p}{2}-1} \partial_x \|x - y\|^2 \\
&= p \|x - y\|^{p-2} (x - y).
\end{align*}
$$

It results that

$$ \|\partial_x c(x, y)\| = p \|x - y\|^{p-1} $$
and therefore

$$ \begin{aligned} \|\nabla c(x, y)\|^2 &= \|\partial_x c(x, y)\|^2 + \|\partial_y c(x, y)\|^2 \\ &= 2p^2\|x - y\|^{2p-2}. \end{aligned} $$

We then conclude by mean value theorem and the fact that $x, y, \in \mathbb{B}_R$.

□

**Lemma D.2.** For any functions $f, g: \mathbb{R}^d \to \mathbb{R}$ and subset $U$ of $\mathbb{R}^d$

$$ \left| \sup_{x \in U} f(x) - \sup_{x \in U} g(x) \right| \le \sup_{x \in U} |f(x) - g(x)| $$

and

$$ \left| \inf_{x \in U} f(x) - \inf_{x \in U} g(x) \right| \le \sup_{x \in U} |f(x) - g(x)| $$

*Proof.* Let $x \in U$, then

$$ f(x) - \sup_{t \in U} g(t) \le f(x) - g(x) \le \sup_{x \in U} |f(x) - g(x)|. $$

This being true for all $x \in U$ means that

$$ \sup_{x \in U} f(x) - \sup_{t \in U} g(x) \le \sup_{x \in U} |f(x) - g(x)|. $$

By inverting the role of $f$ and $g$ we can then show the first result. As for the second result for all $\bar{x} \in U$ we have

$$ \begin{aligned} \sup_{x \in U} |f(x) - g(x)| &\ge f(\bar{x}) - g(\bar{x}) \\ &\ge \inf_{x \in U} f(x) - g(\bar{x}) \end{aligned} $$

This being true for all $\bar{x} \in U$ means that

$$ \sup_{x \in U} |f(x) - g(x)| \ge \inf_{x \in U} f(x) - \inf_{x \in U} g(x). $$

By inverting the role of $f$ and $g$ we can show the second result.

□

## D.2 Lower bound Lemma

In this subsection we prove that the worst-case lower bound given in [7] also holds for Wasserstein sparse barycenters.

**Lemma D.3.** Let $\mathbb{D} = W_2^2$, the number of measures $L \in \mathbb{N}^*$ and the dimension $d \in \mathbb{N}^*$. Assume $N \ge 3$ and $n \ge 3435N$ and let $\sum_{i=1}^N \pi_i^n \delta_{y_i^n} \in \mathcal{M}_1^N(\mathbb{R}^d)$ be a random measure depending on $nL$ random variables in $\mathbb{B}_R$. Then there exists $\mu^1, \dots, \mu^L \in \mathcal{M}_1(\mathbb{R}^d)$ such that the estimation error of Corollary 3.4 verifies

$$ \mathbb{E} \left[ \mathrm{F}_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i^n \delta_{y_i^n} \right) - \min_{(Y, \pi) \in \mathbb{B}_R^N \times \Delta_N} \mathrm{F}_{\mathbb{D}} \left( \mu^1, \dots, \mu^L, \sum_{i=1}^N \pi_i \delta_{y_i} \right) \right] \ge C_d \sqrt{\frac{N^{1-\frac{4}{d}}}{n}} \quad (\text{D.1}) $$

where $C_d$ is an explicit constant which only depends on $d$.

*Proof.* Let us first consider the case $L = 1$. By [7, Theorem 1], there exists a measure $\mu \in \mathcal{M}_1(\mathbb{R}^d)$ and an optimal quantizer $(Y^*, \pi^*) \in \underset{(Y, \pi) \in (\mathbb{R}^d)^N \times \Delta_N}{\operatorname{argmin}} W_2^2(\mu, \sum_{i=1}^N \pi_i \delta_{y_i})$ such that

$$ \mathbb{E} \left[ \min_{\pi \in \Delta_N} W_2^2 \left( \mu, \sum_{i=1}^N \pi_i \delta_{y_i^n} \right) - W_2^2 \left( \mu, \sum_{i=1}^N \pi_i^* \delta_{y_i^*} \right) \right] \gtrsim \sqrt{\frac{N^{1-\frac{4}{d}}}{n}}. $$
Since almost surely

$$ \min_{\pi \in \Delta_N} W_2^2 \left( \mu, \sum_{i=1}^N \pi_i \delta_{y_i^n} \right) \le W_2^2 \left( \mu, \sum_{i=1}^N \pi_i^n \delta_{y_i^n} \right), $$

we deduce that

$$ \mathbb{E} \left[ W_2^2 \left( \mu, \sum_{i=1}^N \pi_i^n \delta_{y_i^n} \right) - W_2^2 \left( \mu, \sum_{i=1}^N \pi_i^* \delta_{y_i^*} \right) \right] \gtrsim \sqrt{\frac{N^{1-\frac{4}{d}}}{n}}. $$

The result for $L > 1$ in (D.1) then follows by choosing $\mu^1 = \cdots = \mu^L = \mu$. $\square$

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

## D.3 Optimal quantization

For completeness, we recall and prove this well known fact about optimal quantization.

**Lemma D.4.** Let $\mu$ be an arbitrary measure and $Y \in (\mathbb{R}^d)^N \setminus D_N$, then for all $p \ge 1$

$$ \min_{\pi \in \Delta_N} W_p^p \left( \mu, \sum_{i=1}^N \pi_i \delta_{y_i} \right) = \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \|x - y_i\|^p d\mu(x). $$

*Proof.* Consider the functional

$$ (w, \pi) \in \mathbb{R}^N \times \Delta_N \mapsto \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \{\|x - y_i\|^p - w_i\} d\mu(x) + \sum_{i=1}^N \pi_i w_i. $$

It is concave in $w$ and convex in $\pi$, and we can therefore swap the min and the max in the following problem

$$ \min_{\pi \in \Delta_N} \max_{w \in \mathbb{R}^N} \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \{\|x - y_i\|^p - w_i\} d\mu(x) + \sum_{i=1}^M \pi_i w_i. \quad (D.2) $$

Notice that a first order condition with respect to $\pi$ implies $w^* = 0_{\mathbb{R}^N}$ at optimality. Let $\pi^* \in \Delta_N$ be a minimizer of (D.2). In order to use the duality expression (A.2), we consider $\tilde{\pi} \in \text{int}(\Delta_M)$ the vector of weights $\pi^*$ whose null components have been removed. Following this construction and with a slight abuse of notation, we reorder the point cloud $Y := (y_1, \dots, y_N)$ so that the components $y_i$ corresponding to null weights have indexes $i \in [M + 1, N]$.

$$ \begin{align*} \min_{\pi \in \Delta_N} W_p^p \left( \mu, \sum_{i=1}^N \pi_i \delta_{y_i} \right) &= W_p^p \left( \mu, \sum_{i=1}^M \tilde{\pi}_i \delta_{y_i} \right) \\ &= \min_{\pi \in \Delta_M} \max_{w \in \mathbb{R}^M} \int_{\mathbb{R}^d} \min_{i=1, \dots, M} \{\|x - y_i\|^p - w_i\} d\mu(x) + \sum_{i=1}^M \pi_i w_i \\ &= \max_{w \in \mathbb{R}^M} \min_{\pi \in \Delta_M} \int_{\mathbb{R}^d} \min_{i=1, \dots, M} \{\|x - y_i\|^p - w_i\} d\mu(x) + \sum_{i=1}^M \pi_i w_i \\ &= \int_{\mathbb{R}^d} \min_{i=1, \dots, M} \|x - y_i\|^p d\mu(x). \end{align*} $$

We can then consider two cases. First, suppose that $\mu$ does not give mass to the Voronoï set $\{x \in \mathbb{R}^d \mid \|x - y_i\|^p \le \|x - y_j\|^p, j \in [1, N]\}$ for all $i \in [M + 1, N]$; this may be the case for example if $\mu$ is discrete. Then we have

$$ \int_{\mathbb{R}^d} \min_{i=1, \dots, M} \|x - y_i\|^p d\mu(x) = \int_{\mathbb{R}^d} \min_{i=1, \dots, N} \|x - y_i\|^p d\mu(x). $$
Second, suppose that there exists $i \in [M + 1, N]$ such that $\mu$ gives mass to a set $\{x \in \mathbb{R}^d | \|x - y_i\|^p \le \|x - y_j\|^p, j \in [1, N]\}$. By contradiction, suppose that $M < N$. Because $Y \notin D_N$ we then have

$$ 
\begin{aligned}
\int_{\mathbb{R}^d} \min_{i=1,\dots,M} \|x - y_i\|^p d\mu(x) &> \int_{\mathbb{R}^d} \min_{i=1,\dots,N} \|x - y_i\|^p d\mu(x) \\
&= \max_{w \in \mathbb{R}^N} \min_{\pi \in \Delta_N} \int_{\mathbb{R}^d} \min_{i=1,\dots,N} \{\|x - y_i\|^p - w_i\} d\mu(x) + \sum_{i=1}^N \pi_i w_i \\
&= \min_{\pi \in \Delta_N} W_p^p \left( \mu, \sum_{i=1}^N \pi_i \delta_{y_i} \right)
\end{aligned}
$$

which contradicts the optimality of the minimum in $\pi$. We thus get the result. $\square$

## References

[1] Martial Agueh and Guillaume Carlier. Barycenters in the Wasserstein space. SIAM Journal on Mathematical Analysis, 43(2):904–924, 2011.

[2] Mengting Ai, Tianxin Wei, Yifan Chen, Zhichen Zeng, Ritchie Zhao, Girish Varatkar, Bita Darvish Rouhani, Xianfeng Tang, Hanghang Tong, and Jingrui He. Resmoe: Space-efficient compression of mixture of experts llms via residual restoration. arXiv preprint arXiv:2503.06881, 2025.

[3] Jason M Altschuler and Enric Boix-Adsera. Wasserstein barycenters can be computed in polynomial time in fixed dimension. Journal of Machine Learning Research, 22(44):1–19, 2021.

[4] Jason M Altschuler and Enric Boix-Adsera. Wasserstein barycenters are np-hard to compute. SIAM Journal on Mathematics of Data Science, 4(1):179–203, 2022.

[5] Ethan Anderes, Steffen Borgwardt, and Jacob Miller. Discrete Wasserstein barycenters: Optimal transport for discrete data. Mathematical Methods of Operations Research, 84:389–409, 2016.

[6] Julio Backhoff, Joaquin Fontbona, Gonzalo Rios, and Felipe Tobar. Stochastic gradient descent for barycenters in Wasserstein space. Journal of Applied Probability, pages 1–29, 2024.

[7] Peter L Bartlett, Tamás Linder, and Gábor Lugosi. The minimax distortion redundancy in empirical quantizer design. IEEE Transactions on Information theory, 44(5):1802–1813, 2002.

[8] Gérard Biau, Luc Devroye, and Gábor Lugosi. On the performance of clustering in hilbert spaces. IEEE Transactions on Information Theory, 54(2):781–790, 2008.

[9] Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and Radon Wasserstein barycenters of measures. Journal of Mathematical Imaging and Vision, 51:22–45, 2015.

[10] Steffen Borgwardt and Stephan Patterson. On the computational complexity of finding a sparse Wasserstein barycenter. Journal of Combinatorial Optimization, 41(3):736–761, 2021.

[11] Camilla Brizzi, Gero Friesecke, and Tobias Ried. p-Wasserstein barycenters. Nonlinear Analysis, 251:113687, 2025.

[12] Luc Brogat-Motte, Rémi Flamary, Céline Brouard, Juho Rousu, and Florence d'Alché Buc. Learning to predict graphs with fused Gromov-Wasserstein barycenters. In International Conference on Machine Learning, pages 2321–2335. PMLR, 2022.

[13] Guillaume Carlier, Alex Delalande, and Quentin Merigot. Quantitative stability of barycenters in the Wasserstein space. Probability Theory and Related Fields, 188(3):1257–1286, 2024.

[14] Sinho Chewi, Tyler Maunu, Philippe Rigollet, and Austin J Stromme. Gradient descent algorithms for Bures-Wasserstein barycenters. In Conference on Learning Theory, pages 1276–1304. PMLR, 2020.
[15] Sinho Chewi, Jonathan Niles-Weed, and Philippe Rigollet. Statistical optimal transport. arXiv preprint arXiv:2407.18163, 3, 2024.

[16] Lénaïc Chizat. Doubly regularized entropic Wasserstein barycenters. arXiv preprint arXiv:2303.11844, 2023.

[17] Sebastian Claici, Edward Chien, and Justin Solomon. Stochastic Wasserstein barycenters. In International Conference on Machine Learning, pages 999–1008. PMLR, 2018.

[18] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.

[19] Marco Cuturi and Arnaud Doucet. Fast computation of Wasserstein barycenters. In International conference on machine learning, pages 685–693. PMLR, 2014.

[20] Eustasio Del Barrio, Alberto González Sanz, and Jean-Michel Loubes. Central limit theorems for semi-discrete Wasserstein distances. Bernoulli, 30(1):554–580, 2024.

[21] Eustasio del Barrio, Alberto González Sanz, Jean-Michel Loubes, and Jonathan Niles-Weed. An improved central limit theorem and fast convergence rates for entropic transportation costs. SIAM Journal on Mathematics of Data Science, 5(3):639–669, 2023.

[22] Darina Dvinskikh. Stochastic approximation versus sample average approximation for population Wasserstein barycenters. arXiv preprint arXiv:2001.07697, 2020.

[23] Pavel Dvurechenskii, Darina Dvinskikh, Alexander Gasnikov, Cesar Uribe, and Angelia Nedich. Decentralize and randomize: Faster algorithm for Wasserstein barycenters. Advances in Neural Information Processing Systems, 31, 2018.

[24] Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's algorithm. In International conference on machine learning, pages 1367–1376. PMLR, 2018.

[25] Lasalle Etienne, Vaudaine Rémi, Vayer Titouan, Borgnat Pierre, Gonçalves Paulo, Gribonval Rémi, and Karsai Márton. Pasco (parallel structured coarsening): an overlay to speed up graph clustering algorithms. Machine Learning, 114(10):212, 2025.

[26] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer. Pot: Python optimal transport. Journal of Machine Learning Research, 22(78):1–8, 2021.

[27] Nicolas Fournier and Arnaud Guillin. On the rate of convergence in Wasserstein distance of the empirical measure. Probability theory and related fields, 162(3):707–738, 2015.

[28] Erell Gachon, Jérémie Bigot, Elsa Cazelles, Aguirre Mimoun, and Jean-Philippe Vial. Low dimensional representation of multi-patient flow cytometry datasets using optimal transport for minimal residual disease detection in leukemia, 07 2024.

[29] Aude Genevay, Lénaïc Chizat, Francis Bach, Marco Cuturi, and Gabriel Peyré. Sample complexity of Sinkhorn divergences. In Kamalika Chaudhuri and Masashi Sugiyama, editors, Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pages 1574–1583. PMLR, 2019.

[30] Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach. Stochastic optimization for large-scale optimal transport. Advances in neural information processing systems, 29, 2016.

[31] Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with Sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pages 1608–1617. PMLR, 2018.
[32] Evarist Giné and Richard Nickl. Mathematical foundations of infinite-dimensional statistical models. Cambridge university press, 2021.

[33] László Györfi. Principles of nonparametric learning. Springer, 2002.

[34] Florian Heinemann, Axel Munk, and Yoav Zemel. Randomized Wasserstein barycenter computation: Resampling with statistical guarantees. SIAM Journal on Mathematics of Data Science, 4(1):229–259, 2022.

[35] Shayan Hundrieser, Thomas Staudt, and Axel Munk. Empirical optimal transport between different measures adapts to lower complexity. In Annales de l'Institut Henri Poincare (B) Probabilites et statistiques, volume 60, pages 824–846. Institut Henri Poincaré, 2024.

[36] Hicham Janati, Marco Cuturi, and Alexandre Gramfort. Debiased Sinkhorn barycenters. In International Conference on Machine Learning, pages 4692–4701. PMLR, 2020.

[37] S. Kolouri, K. Nadjahi, U. Simsekli, R. Badeau, and G. Rohde. Generalized Sliced Wasserstein distances. Advances in neural information processing systems, 32, 2019.

[38] Chun-Yao Lee and Guang-Lin Zhuo. Barycenter-based optimal transport neighborhood learning network for multi-source cross-machine bearing fault diagnosis. IEEE Transactions on Instrumentation and Measurement, 2025.

[39] Lingxiao Li, Aude Genevay, Mikhail Yurochkin, and Justin M Solomon. Continuous regularized Wasserstein barycenters. Advances in Neural Information Processing Systems, 33:17755–17765, 2020.

[40] Tianyi Lin, Zeyu Zheng, Elynn Chen, Marco Cuturi, and Michael I Jordan. On projection robust optimal transport: Sample complexity and model misspecification. In International Conference on Artificial Intelligence and Statistics, pages 262–270. PMLR, 2021.

[41] Pascal Massart. Concentration inequalities and model selection: Ecole d'Eté de Probabilités de Saint-Flour XXXIII-2003. Springer, 2007.

[42] Robert J McCann. A convexity principle for interacting gases. Advances in mathematics, 128(1):153–179, 1997.

[43] Daniel Mimouni, Welington de Oliveira, and Gregorio M Sempere. On the computation of constrained wasserstein barycenters. hal-05216990, 2025.

[44] Eduardo F Montesuma and Fred-Maurice Ngolè Mboula. Wasserstein barycenter transport for acoustic adaptation. In ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3405–3409. IEEE, 2021.

[45] Eduardo Fernandes Montesuma and Fred Maurice Ngole Mboula. Wasserstein barycenter for multi-source domain adaptation. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 16785–16793, 2021.

[46] Kimia Nadjahi, Alain Durmus, Lénaïc Chizat, Soheil Kolouri, Shahin Shahrampour, and Umut Simsekli. Statistical and topological properties of Sliced probability divergences. Advances in Neural Information Processing Systems, 33:20802–20812, 2020.

[47] Michael Kwok-Po Ng. A note on constrained k-means algorithms. Pattern Recognition, 33(3):515–519, 2000.

[48] Gilles Pagès, Huyên Pham, and Jacques Printems. Optimal quantization methods and applications to numerical problems in finance. Handbook of computational and numerical methods in finance, pages 253–297, 2004.

[49] Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.

[50] Léo Portales, Elsa Cazelles, and Edouard Pauwels. On the sequential convergence of Lloyd's algorithms. arXiv preprint arXiv:2405.20744, 2024.
[51] Giovanni Puccetti, Ludger Rüschendorf, and Steven Vanduffel. On the computation of Wasserstein barycenters. *Journal of Multivariate Analysis*, 176:104581, 2020.

[52] Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot. Wasserstein barycenter and its application to texture mixing. In *International conference on scale space and variational methods in computer vision*, pages 435–446. Springer, 2011.

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

[54] Bodhisattva Sen. A gentle introduction to empirical process theory and applications. Lecture Notes, Columbia University, 11:28–29, 2018.

[55] Richard Sinkhorn and Paul Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics, 21(2):343–348, 1967.

[56] Max Sommerfeld and Axel Munk. Inference for empirical Wasserstein distances on finite spaces. Journal of the Royal Statistical Society Series B: Statistical Methodology, 80(1):219–238, 2018.

[57] Sanvesh Srivastava, Volkan Cevher, Quoc Dinh, and David Dunson. Wasp: Scalable Bayes via barycenters of subset posteriors. In Artificial intelligence and statistics, pages 912–920. PMLR, 2015.

[58] Cédric Villani. Optimal Transport: Old and New, volume 338 of Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 2008.

[59] Cédric Villani. Topics in optimal transportation, volume 58. American Mathematical Soc., 2021.

[60] Xianliang Xu and Zhongyi Huang. Central limit theorem for the sliced 1-Wasserstein distance and the max-sliced 1-Wasserstein distance. arXiv preprint arXiv:2205.14624, 2022.

[61] Qingyuan Yang and Hu Ding. Approximate algorithms for $k$-sparse Wasserstein barycenter with outliers. arXiv preprint arXiv:2404.13401, 2024.