Code
renderShapeDemo(mu1_s, sigma1_s, mu2_s, sigma2_s, lambda_s)Lecture 12: Wasserstein Geometry — Wasserstein Barycenters
Fréchet means in Wasserstein space
1 Learning Goals
By the end of this lecture, learners should be able to:
- Define Wasserstein barycenters as Fréchet means in \((\mathcal{P}_2(\mathcal{X}), W_2)\) and explain how they generalize Euclidean weighted averages.
- Derive the closed-form quantile-averaging formula for 1D Wasserstein barycenters.
- Compute the Wasserstein barycenter of Gaussian measures via the Bures–Wasserstein fixed-point equation.
- State the existence and consistency theorems for sample Wasserstein barycenters.
- Summarize convergence rates in one dimension (\(n^{-1/2}\)) and for Gaussian barycenters.
- Describe the three main computational approaches: linear programming, fixed-point iteration, and entropic regularization.
2 Definition
Wasserstein barycenters are Fréchet means in Wasserstein space. They are the appropriate notion of center when the observations are probability measures and the loss is squared optimal-transport distance.
Let \((\mathcal{X}, d)\) be a complete separable metric space and let \(\mathcal{P}_2(\mathcal{X})\) be equipped with \(W_2\). If \(\Lambda\) is a probability distribution on \(\mathcal{P}_2(\mathcal{X})\) satisfying
\[ \int_{\mathcal{P}_2(\mathcal{X})} W_2^2(\mu, \mu_0)\, d\Lambda(\mu) < \infty \]
for some, hence every, \(\mu_0 \in \mathcal{P}_2(\mathcal{X})\), then a population Wasserstein barycenter of \(\Lambda\) is any minimizer
\[ \bar{\mu} \in \arg\min_{\nu \in \mathcal{P}_2(\mathcal{X})} \Phi(\nu), \qquad \Phi(\nu) = \int_{\mathcal{P}_2(\mathcal{X})} W_2^2(\nu, \mu)\, d\Lambda(\mu). \]
Equivalently, if \(M \sim \Lambda\) is a random probability measure, then \(\bar{\mu}\) minimizes \(\mathbb{E}\, W_2^2(\nu, M)\) over \(\nu \in \mathcal{P}_2(\mathcal{X})\).
In the finite weighted case, where \(\Lambda = \sum_{i=1}^n \lambda_i \delta_{\mu_i}\) with \(\lambda_i > 0\) and \(\sum_i \lambda_i = 1\):
\[ \bar{\mu} \in \arg\min_{\nu \in \mathcal{P}_2(\mathcal{X})} \sum_{i=1}^n \lambda_i W_2^2(\nu, \mu_i). \]
Taking \(\lambda_i=1/n\) gives the Fréchet mean of the measures \(\mu_1, \ldots, \mu_n\).
Note
If \(\mu_i = \delta_{x_i}\) are point masses in Euclidean space, then the barycenter is \(\delta_{\bar{x}}\) with \(\bar{x} = \sum_i \lambda_i x_i\), so the construction extends the ordinary weighted mean.
Barycenter vs. Density Average. A key insight behind Wasserstein barycenters is that they often preserve the shape of the input distributions, whereas a naive pointwise average of densities does not in general. The figure below illustrates this with two 1D Gaussians.
TipWhat to observe
- Barycenter (solid green): Always a single Gaussian — the Wasserstein barycenter preserves the shape of the input family. The mean and standard deviation are the weighted averages of the inputs.
- Density average (dotted purple): When the two Gaussians are well-separated (|μ₁ − μ₂| ≫ σ₁ + σ₂), the pointwise average of the densities becomes bimodal — it has two peaks and does not look like either input distribution.
- Why this matters: The Wasserstein barycenter interpolates between distributions by averaging their quantile functions, which respects the geometry of the space. The density average interpolates pointwise on the x-axis, which ignores the underlying metric structure and can produce distributions that lie far from any individual input in Wasserstein distance.
3 One-Dimensional Formula
When \(\mathcal{X} = \mathbb{R}\), the quantile representation gives an explicit closed-form solution. Let \(Q_\mu = F_\mu^{-1}\) denote the quantile function of \(\mu\). Since the map \(\mu \mapsto Q_\mu\) embeds \(\mathcal{P}_2(\mathbb{R})\) isometrically into the closed convex cone of nondecreasing functions in \(L^2(0, 1)\),
\[ W_2^2(\mu, \nu) = \int_0^1 \{Q_\mu(u) - Q_\nu(u)\}^2\, du. \]
The barycenter is therefore unique and has quantile function
\[ Q_{\bar{\mu}}(u) = \int_{\mathcal{P}_2(\mathbb{R})} Q_\mu(u)\, d\Lambda(\mu), \qquad 0 < u < 1, \]
provided the right side is square-integrable. In the finite case this becomes
\[ Q_{\bar{\mu}}(u) = \sum_{i=1}^n \lambda_i Q_{\mu_i}(u). \]
The pointwise average of nondecreasing quantile functions is again nondecreasing, so it defines a valid probability distribution.
TipWhy 1D is special
This identity is the main reason univariate Wasserstein means are much simpler than their higher-dimensional analogues: the barycenter is simply the pointwise quantile average. This underlies many statistical procedures for distribution-valued data (Petersen and Müller 2016; Panaretos and Zemel 2020).
4 Gaussian Barycenters
If \(\mu_i = N(m_i, \Sigma_i)\) are Gaussian measures on \(\mathbb{R}^D\), then their Wasserstein barycenter is again Gaussian, \(\bar{\mu} = N(\bar{m}, \bar{\Sigma})\), with
\[ \bar{m} = \sum_{i=1}^m \lambda_i m_i, \]
and covariance matrix determined by the Bures–Wasserstein fixed-point equation
\[ \bar{\Sigma} = \sum_{i=1}^m \lambda_i \bigl(\bar{\Sigma}^{1/2} \Sigma_i \bar{\Sigma}^{1/2}\bigr)^{1/2}. \]
This example shows that Wasserstein barycenters average both locations and distributional shapes, not just pointwise density values (Agueh and Carlier 2011; Panaretos and Zemel 2020).
NoteOptimal transport between Gaussians
For nondegenerate Gaussians \(N(m_0, \Sigma_0)\) and \(N(m_1, \Sigma_1)\), the optimal transport map is affine:
\[ T_{0 \to 1}(x) = m_1 + \Sigma_0^{-1/2} \bigl(\Sigma_0^{1/2} \Sigma_1 \Sigma_0^{1/2}\bigr)^{1/2} \Sigma_0^{-1/2}(x - m_0). \]
5 Existence and Uniqueness
Theorem 1 Let \(p \ge 1\) and let \((E, d)\) be a separable locally compact geodesic space. Let \(\Lambda\) be a probability measure on \(\mathcal{W}_p(E)\) such that
\[ \int_{\mathcal{W}_p(E)} W_p^p(\mu, \mu_0)\, d\Lambda(\mu) < \infty \]
for some, hence every, \(\mu_0 \in \mathcal{W}_p(E)\). Then there exists at least one barycenter \(\bar{\mu}_\Lambda \in \mathcal{W}_p(E)\).
Uniqueness is more subtle. Wasserstein barycenters may not be unique in general — \(\mathcal{P}_2(\mathbb{R}^d)\) with \(d \ge 2\) has nonnegative Alexandrov curvature rather than the nonpositive curvature that would force strict convexity of squared distance.
- In one dimension, uniqueness follows from the Hilbert-space quantile embedding.
- In \(\mathbb{R}^d\), a standard sufficient condition is that at least one input measure is absolutely continuous with respect to Lebesgue measure (Agueh and Carlier 2011).
6 Consistency of Sample Barycenters
Let \(M_1, M_2, \ldots\) be i.i.d. random probability measures with law \(\Lambda\), and define the empirical law
\[ \Lambda_n = \frac{1}{n}\sum_{i=1}^n \delta_{M_i}. \]
A sample Wasserstein barycenter is any minimizer
\[ \hat{\mu}_n \in \arg\min_{\nu \in \mathcal{P}_2(\mathcal{X})} \frac{1}{n}\sum_{i=1}^n W_2^2(\nu, M_i). \]
Let \(\mathsf{W}_2\) denote the Wasserstein distance on \(\mathcal{W}_2(E)\) itself (using \(W_2\) as the ground metric). By the strong law of large numbers and Villani’s Theorem 6.9,
\[ \mathsf{W}_2(\Lambda_n, \Lambda) \to 0 \qquad \text{a.s.} \]
Let \(\operatorname{Bar}_2(\Lambda_n)\) denote the set of empirical 2-Wasserstein barycenters, i.e.,
\[ \operatorname{Bar}_2(\Lambda_n) \coloneqq \arg\min_{\nu \in \mathcal{P}_2(\mathcal{X})} \frac{1}{n}\sum_{i=1}^n W_2^2(\nu, M_i). \]
Theorem 2 If the population barycenter \(\bar{\mu}_\Lambda\) is unique and sample barycenters \(\hat{\mu}_n \in \operatorname{Bar}_2(\Lambda_n)\), then
\[ W_2(\hat{\mu}_n, \bar{\mu}_\Lambda) \to 0 \qquad \text{a.s.} \]
In \(\mathbb{R}^d\), the uniqueness condition holds if \(\Lambda\) assigns positive probability to absolutely continuous ground-space distributions (Le Gouic and Loubes 2017, Proposition 6).
7 Convergence Rates
7.1 One Dimension
In \(\mathcal{W}_2(\mathbb{R})\), the quantile representation gives
\[ Q_{\hat{\mu}_n}(u) = \frac{1}{n}\sum_{i=1}^n Q_{M_i}(u), \qquad Q_{\bar{\mu}}(u) = \mathbb{E}\{Q_M(u)\}. \]
Hence
\[ \mathbb{E}\, W_2^2(\hat{\mu}_n, \bar{\mu}) = \frac{1}{n}\int_0^1 \operatorname{Var}\{Q_M(u)\}\, du, \]
so the squared error is \(O(n^{-1})\) and the distance is \(O(n^{-1/2})\). This parametric rate reflects the Hilbert-space structure of the quantile embedding — the barycenter is simply a sample mean in \(L^2(0,1)\).
7.2 Higher Dimensions
\(\alpha\)-Strong Convexity and \(\beta\)-Smoothness. The convergence rate of sample Wasserstein barycenters depends on the geometric regularity of the transport maps pushing \(\bar{\mu}\) to each input measure. The key structural conditions are strong convexity and smoothness of the potentials whose gradients realize these maps.
Definition 1 Let \(\alpha > 0\). A differentiable function \(f: \mathbb{R}^d \to \mathbb{R}\) is \(\alpha\)-strongly convex if for all \(x, y \in \mathbb{R}^d\),
\[ f(y) \ge f(x) + \langle\nabla f(x),\, y - x\rangle + \frac{\alpha}{2}\|y - x\|^2. \]
Equivalently, \(\nabla^2 f(x) \succeq \alpha I_d\) (in the sense of Loewner order) wherever the Hessian exists, meaning all eigenvalues of the Hessian are bounded below by \(\alpha\).
Strong convexity quantitatively strengthens ordinary convexity (\(f(y) \ge f(x) + \langle\nabla f(x), y - x\rangle\)). The quadratic penalty \(\frac{\alpha}{2}\|y-x\|^2\) ensures the function curves at least as sharply as \(\frac{\alpha}{2}\|x\|^2\). Geometrically, the gradient map \(\nabla f\) is strictly expanding: \(\langle\nabla f(y) - \nabla f(x), y - x\rangle \ge \alpha\|y - x\|^2\), which guarantees injectivity of \(\nabla f\) and a well-behaved inverse.
Definition 2 Let \(\beta > 0\). A differentiable function \(f: \mathbb{R}^d \to \mathbb{R}\) is \(\beta\)-smooth if its gradient is \(\beta\)-Lipschitz continuous:
\[ \|\nabla f(x) - \nabla f(y)\| \le \beta\,\|x - y\|, \qquad \forall x, y \in \mathbb{R}^d. \]
Equivalently, \(\nabla^2 f(x) \preceq \beta I_d\) wherever the Hessian exists, and
\[ f(y) \le f(x) + \langle\nabla f(x),\, y - x\rangle + \frac{\beta}{2}\|y - x\|^2. \]
Whereas strong convexity provides a lower quadratic bound, smoothness provides an upper quadratic bound. Together, the two conditions sandwich \(f\) between quadratics with curvatures \(\alpha\) and \(\beta\). The ratio \(\beta/\alpha \ge 1\) is the condition number of \(f\).
NoteConnection to optimal transport
When each \(\mu \in \operatorname{supp}(\Lambda)\) is the pushforward of \(\bar{\mu}\) by \(T_\mu = \nabla\phi_\mu\) — the gradient of a convex potential — the regularity of \(\phi_\mu\) controls how much the geometry of \(\mathcal{P}_2(\mathbb{R}^d)\) near \(\bar{\mu}\) resembles a Hilbert space. Brenier’s theorem (Panaretos and Zemel 2020, sec. 2.3) guarantees that optimal transport maps between absolutely continuous measures are exactly of this gradient-of-convex-potential form.
TipTwo notions of strong convexity
The \(\alpha\)-strong convexity defined here is the classical Euclidean notion — a lower quadratic bound on the function via its gradient. This is distinct from the geodesic \(\lambda\)-strong convexity introduced in Lecture 2, which applies to functions on arbitrary geodesic metric spaces:
\[ f(\gamma(t)) \le (1-t)f(\gamma(0)) + t f(\gamma(1)) - \lambda\,t(1-t)\,d^2(\gamma(0), \gamma(1)), \]
where \(\gamma: [0,1] \to \mathcal{M}\) is a geodesic. The two notions coincide when \(\mathcal{M} = \mathbb{R}^d\) with the Euclidean metric and \(f\) is differentiable — in that case the geodesic inequality reduces to the gradient inequality with \(\lambda = \alpha\).
In \(d \ge 2\), the Wasserstein space lacks the flat Hilbert geometry of the 1D case, and convergence rates depend on the curvature of the transport maps from the barycenter. Le Gouic et al. (2023) showed that if each \(\mu \in \operatorname{supp}(\Lambda)\) is the pushforward of \(\bar{\mu}\) by the gradient of an \(\alpha\)-strongly convex and \(\beta\)-smooth potential \(\phi_\mu\), i.e.
\[ \mu = (\nabla\phi_\mu)_{\#}\bar{\mu}, \]
then the sample barycenter attains the parametric rate, with a constant governed by the gap \(\beta - \alpha\).
Theorem 3 If \(\beta - \alpha < 1\), then
\[ \mathbb{E}\, W_2^2(\hat{\mu}_n, \bar{\mu}) \le \frac{4\sigma^2}{(1 - \beta + \alpha)^2\, n}, \]
where \(\sigma^2 = \int_{\mathcal{P}_2(\mathbb{R}^d)} W_2^2(\mu, \bar{\mu})\, d\Lambda(\mu)\) is the population variance in Wasserstein space.
7.3 Gaussian Case
For Gaussian measures, the transport maps are affine and the regularity parameters are determined by the eigenvalue spread of the covariance matrices. If all input covariances have eigenvalues in \([\kappa_0, \kappa_1]\) and we set \(\kappa = \kappa_1 / \kappa_0 \ge 1\), then:
\[ \mathbb{E}\, W_2^2(\hat{\mu}_n, \bar{\mu}) \le \frac{4\sigma^2}{(1 - \kappa + \kappa^{-1})^2\, n}, \]
provided \(\kappa - \kappa^{-1} < 1\). This is exactly the higher-dimensional bound with the transport-map regularity gap \(\beta - \alpha\) replaced by the eigenvalue condition-number gap \(\kappa - \kappa^{-1}\). The condition \(\kappa - \kappa^{-1} < 1\) (equivalently \(\kappa < \frac{1+\sqrt{5}}{2} \approx 1.618\)) requires the covariance matrices to be moderately well-conditioned — extreme anisotropy in the input measures can degrade the convergence rate.
8 Computation
Computing Wasserstein barycenters is substantially harder than computing Euclidean means. Three main paradigms exist, each with different trade-offs between accuracy, speed, and scalability. For a comprehensive treatment, see Peyré and Cuturi (2019), Chapters 4 and 9.
8.1 Linear Programming on Discrete Supports
For discrete inputs, the barycenter problem is a linear program (LP) once the barycenter support locations are fixed. Suppose each \(\mu_i\) is supported on \(\{x_{i1},\ldots,x_{i n_i}\}\) with probability vector \(a_i\), and restrict the barycenter to \(\{z_1,\ldots,z_L\}\) with unknown weight vector \(b\). Let \(C_{i\ell k} = \|z_\ell - x_{ik}\|^2\) be the transport cost. The fixed-support barycenter LP is:
\[ \begin{aligned} \min_{b,\;\pi_1,\ldots,\pi_m} \quad & \sum_{i=1}^m \lambda_i \sum_{\ell=1}^{L} \sum_{k=1}^{n_i} C_{i\ell k}\,\pi_{i\ell k} \\ \text{subject to} \quad & \pi_i \mathbf{1}_{n_i} = b, \qquad \pi_i^\mathsf{T} \mathbf{1}_L = a_i, \qquad \pi_i \ge 0, \qquad i=1,\ldots,m, \end{aligned} \]
where \(b \ge 0\), \(\sum_\ell b_\ell = 1\), and \(\pi_i\) is the coupling between the barycenter and \(\mu_i\). The constraints enforce that each \(\pi_i\) has the barycenter weights \(b\) as its first marginal and the input weights \(a_i\) as its second marginal — all input measures couple to the same \(b\), which is the discrete barycenter.
Complexity. The formulation has \(L+\sum_i L n_i\) variables before eliminating \(b\), and \(\sum_i(L+n_i)\) displayed marginal equalities (with linear dependencies). Thus even writing the dense coupling variables and costs requires \(O(\sum_i L n_i)\) storage. General-purpose interior-point methods have polynomial worst-case complexity (Nesterov and Nemirovskii 1994), while transportation and network-flow structure can be exploited by specialized solvers (Ahuja and Orlin 1992). The actual running time depends strongly on sparsity and the solver, so there is no universal cubic running-time formula in \(mL\bar n\). Exact LP is consequently most useful for modest supports or as a reference solution.
Choice of the barycenter support. The LP formulation requires fixing \(\{z_\ell\}\) in advance. A grid over the convex hull or the union of the input supports gives a restricted, generally approximate barycenter. For quadratic cost, an exact discrete barycenter can be sought on the much larger candidate set of weighted averages \(\sum_i\lambda_i x_{i k_i}\), but that set can contain as many as \(\prod_i n_i\) points (Peyré and Cuturi 2019, sec. 9.2). Alternatively, one can alternate between optimizing the weights/couplings and moving a prescribed number of support locations; this is a free-support, nonconvex optimization rather than a single LP. A coarse or otherwise misspecified candidate support introduces discretization error even when the restricted LP is solved exactly.
8.2 Fixed-Point Iteration
When measures are absolutely continuous (or approximated by large discrete samples), the fixed-point iteration of Álvarez-Esteban et al. (2016) avoids discretizing the barycenter support in advance by iteratively updating the barycenter via averaged transport maps.
NoteAlgorithm: Fixed-Point Barycenter Iteration
Input: Measures \(\mu_1,\ldots,\mu_m\), weights \(\lambda_i\), initial guess \(\nu^{(0)}\).
Repeat for \(t = 0, 1, 2, \ldots\):
Compute transport maps. For each \(i\), compute the optimal transport map \(T_i^{(t)}\) from the current iterate \(\nu^{(t)}\) to \(\mu_i\). In \(\mathbb{R}^d\), if \(\nu^{(t)}\) is absolutely continuous, Brenier’s theorem guarantees \(T_i^{(t)} = \nabla\phi_i^{(t)}\) for a convex potential \(\phi_i^{(t)}\).
Form the barycentric map. Define the weighted average \[ T^{(t)}(x) = \sum_{i=1}^m \lambda_i\, T_i^{(t)}(x). \]
Push forward. Update the iterate: \[ \nu^{(t+1)} \leftarrow T^{(t)}_{\#}\nu^{(t)}. \]
Until \(\|T^{(t)}(x) - x\|\) is sufficiently small (in \(L^2(\nu^{(t)})\)).
Why it works. At the true barycenter \(\bar{\mu}\), the first-order optimality condition is
\[ \sum_{i=1}^m \lambda_i\,T_i(x) = x \qquad \bar{\mu}\text{-a.e.}, \]
where \(T_i\) is the optimal map from \(\bar{\mu}\) to \(\mu_i\) under the absolute-continuity hypotheses above (Agueh and Carlier 2011, Proposition 3.8 and Remark 3.9). Hence \(\bar{\mu}\) is a fixed point of the update \(\nu \mapsto (\sum \lambda_i T_i^\nu)_\#\nu\). The algorithm is a generalized Procrustes procedure: at each step it computes the best way to align the current guess with each input, averages the alignments, and uses the result as the new guess.
Convergence. The rigorous result is more conditional than the fixed-point intuition suggests. Álvarez-Esteban et al. (2016), Theorem 3.6, assume that all targets are absolutely continuous, at least one has a bounded density, and the initial measure is absolutely continuous. They prove tightness of the iterates and show that every subsequential limit is a fixed point of the update. Convergence in \(W_2\) to the barycenter follows if the update operator has a unique fixed point. A fixed point need not be a barycenter in complete generality, so compact support alone does not give the convergence claim. In the Gaussian and, more generally, location-scatter settings treated in that paper, the iteration does converge to the barycenter.
Each exact step requires solving \(m\) optimal-transport problems, which can be done in parallel.
8.3 Entropic Regularization and Sinkhorn Scaling
Entropic regularization (Cuturi and Doucet 2014) makes the discrete coupling subproblems strictly convex and gives them a matrix-scaling structure, enabling fast parallel implementations.
Regularized fixed-support problem. To match the iterative Bregman-projection algorithm below, define the discrete entropic transport cost using entropy relative to a fixed counting reference:
\[ \operatorname{OT}_\varepsilon(b,a_i) = \min_{\pi_i\mathbf 1=b,\;\pi_i^\mathsf T\mathbf 1=a_i} \left\{ \langle C_i,\pi_i\rangle +\varepsilon\sum_{\ell,k}\pi_{i\ell k} \bigl(\log\pi_{i\ell k}-1\bigr) \right\}, \qquad \min_{b\in\Delta_L}\sum_{i=1}^m\lambda_i \operatorname{OT}_\varepsilon(b,a_i). \]
Changing the reference measure changes terms that depend on the unknown \(b\) and therefore changes the regularized barycenter. In particular, the frequently used penalty \(\mathrm{KL}(\pi_i\|b\otimes a_i)\) is not interchangeable with the entropy above when \(b\) is being optimized.
The optimal coupling has the scaling form
\[ \pi_{i\ell k} = u_{i\ell}\, K_{i\ell k}\, v_{ik}, \]
where \(K_{i\ell k} = \exp(-\|z_\ell - x_{ik}\|^2 / \varepsilon)\) is the Gibbs kernel and \((u_i, v_i)\) are positive scaling vectors. This factorized form is the key computational advantage.
NoteAlgorithm: Sinkhorn Barycenter (Iterative Bregman Projections)
Input: Discrete measures \(\mu_i\) with supports \(\{x_{ik}\}\) and weights \(a_i\), target support \(\{z_\ell\}\), regularization \(\varepsilon > 0\), weights \(\lambda_i\).
Initialise: \(u_i\leftarrow\mathbf 1_L\) and \(v_i\leftarrow\mathbf 1_{n_i}\) for all \(i\).
Repeat:
Project onto the known input marginals. For each \(i\), \[ v_i \leftarrow \frac{a_i}{K_i^\mathsf{T}u_i}. \]
Project onto a common barycenter marginal. Compute the weighted geometric mean of the current first marginals, \[ b \leftarrow \prod_{i=1}^m \bigl[u_i\odot(K_i v_i)\bigr]^{\lambda_i}, \qquad u_i \leftarrow \frac{b}{K_i v_i}, \] where products, powers, and divisions are componentwise.
Until the marginal residuals (or the change in \(b\)) are sufficiently small (Benamou et al. 2015, sec. 3.2).
This is an instance of iterative Bregman projections: Step 1 projects all couplings onto their prescribed input marginals, and Step 2 projects them onto the constraint that their first marginals agree. Alternating KL projections converge to the regularized solution when the Gibbs kernels are positive (Benamou et al. 2015). A dense projection sweep costs \(O(\sum_i L n_i)\) arithmetic operations and has the same order of storage if all kernels are materialized. The number of sweeps is problem- and tolerance-dependent and typically increases as \(\varepsilon\) decreases; a general linear rate does not follow from the alternating-projection result.
Regularization bias. For a fixed finite problem, regularized minimizers approach unregularized minimizers as \(\varepsilon\downarrow0\) (with the usual qualification that the unregularized minimizer may not be unique). Smaller \(\varepsilon\) reduces this bias but makes the kernels more ill-conditioned and matrix scaling slower. There is no universal \(W_2(\bar\mu_\varepsilon,\bar\mu)=O(\varepsilon^{1/2})\) barycenter bound, nor a universal \(O(1/\varepsilon)\) iteration count, under only compact support. With the entropy convention above, large \(\varepsilon\) favors diffuse couplings and a high-entropy barycenter; on a fixed finite support the barycenter tends toward uniform weights as the entropy term dominates.
Practical considerations:
Log-domain stabilisation. For small \(\varepsilon\), the kernel entries \(K_{i\ell k} = \exp(-\|z_\ell - x_{ik}\|^2/\varepsilon)\) underflow to zero in floating point. The standard fix is to run Sinkhorn in log-space using the
logsumexpoperation (Schmitzer 2019; Peyré and Cuturi 2019, sec. 4.4).Debiasing. Entropic OT has a nonzero self-cost and commonly produces overly diffuse or blurred barycenters. The Sinkhorn divergence corrects both arguments symmetrically: \[ S_\varepsilon(\mu,\nu) = \operatorname{OT}_\varepsilon(\mu,\nu) - \frac12\operatorname{OT}_\varepsilon(\mu,\mu) - \frac12\operatorname{OT}_\varepsilon(\nu,\nu). \] A debiased Sinkhorn barycenter minimizes \(\sum_i\lambda_iS_\varepsilon(\nu,\mu_i)\). Its important extra term is \(-\frac12\operatorname{OT}_\varepsilon(\nu,\nu)\), which depends on the candidate barycenter. Subtracting only the input self-costs would add a constant and could not alter the minimizer. Debiased barycenters require a modified scaling algorithm (Janati et al. 2020).
Scalability and GPU. The Sinkhorn algorithm uses matrix-vector products and elementwise operations, which parallelize well on GPUs. Dense kernels still require \(O(\sum_iLn_i)\) memory and work per sweep; very large supports require additional structure such as convolutional kernels on grids, low-rank approximations, sparsity, or lazy kernel evaluation.
Free-support barycenters. The Sinkhorn algorithm as stated fixes the barycenter support \(\{z_\ell\}\). When the support should also be learned, one can alternate regularized coupling solves with updates of the support locations (Cuturi and Doucet 2014, sec. 4). The joint problem is nonconvex, so this procedure generally guarantees only a stationary/local solution.
8.4 Summary of Methods
| Method | Accuracy | Speed | Support size | Key reference |
|---|---|---|---|---|
| Fixed-support LP | Exact for the chosen support (to solver tolerance) | General LP solve; coupling storage \(O(\sum_iLn_i)\) | Limited by LP variables and sparsity | Peyré and Cuturi (2019), Ch. 9 |
| Fixed-point iteration | Exact only under its continuous-map assumptions; discrete projection is approximate | Requires \(m\) OT solves per outer iteration | No fixed spatial grid; atom count chosen by the user | Álvarez-Esteban et al. (2016) |
| Entropic Sinkhorn | Solves a regularized fixed-support problem; approaches exact OT as \(\varepsilon\downarrow0\) | Dense sweep \(O(\sum_iLn_i)\); highly parallel | Limited by kernel storage unless structure is exploited | Cuturi and Doucet (2014); Benamou et al. (2015) |
9 Interactive Exploration: 1D Barycenters via Quantile Averaging
Code
renderBarycenterDemo(bres, bary_controls_view, lambda1)
TipTry these experiments
- Equal weights (λ₁ = 0.5): The barycenter sits halfway between the two distributions in Wasserstein space.
- λ₁ = 0 or 1: The barycenter collapses to one of the input distributions.
- Different distribution types: Mix a normal with a uniform — the barycenter quantile averages them pointwise, preserving the nondecreasing property.
- Wasserstein vs. linear averaging: The Wasserstein barycenter of densities is NOT the pointwise average of densities — it’s the average of quantile functions. Try comparing to what a density average would look like.
10 Key Takeaways
- Wasserstein barycenters are Fréchet means in \((\mathcal{P}_2(\mathcal{X}), W_2)\) — minimizers of expected squared \(W_2\) distance.
- In 1D, the barycenter is the pointwise quantile average: \(Q_{\bar{\mu}}(u) = \sum \lambda_i Q_{\mu_i}(u)\).
- Gaussian barycenters have closed-form means and a fixed-point equation for covariances (Bures–Wasserstein).
- Existence holds under moment and compactness conditions; uniqueness requires absolute continuity or 1D structure.
- Sample barycenters are strongly consistent, with \(n^{-1/2}\) rate in 1D and structured rates in higher dimensions.
- Computation ranges from exact LP (small discrete), fixed-point iteration (smooth), to entropic Sinkhorn (large-scale regularized).
11 Exercises
Barycenter of point masses. Show that the Wasserstein barycenter of \(\delta_{x_1}, \ldots, \delta_{x_m}\) with weights \(\lambda_i\) is \(\delta_{\bar{x}}\) where \(\bar{x} = \sum \lambda_i x_i\) is the Euclidean mean. 📝 Show Solution
1D quantile formula. Derive the 1D barycenter formula \(Q_{\bar{\mu}}(u) = \sum \lambda_i Q_{\mu_i}(u)\) from the isometry \(\mu \mapsto Q_\mu\) into \(L^2(0,1)\). 📝 Show Solution
Gaussian fixed-point. For two univariate Gaussians \(N(m_1, \sigma_1^2)\) and \(N(m_2, \sigma_2^2)\), solve the Bures–Wasserstein equation to find the barycenter variance. 📝 Show Solution
Entropic bias-variance. Explain how the entropic regularization parameter \(\varepsilon\) trades bias for variance in barycenter computation. What happens as \(\varepsilon \to 0\) and \(\varepsilon \to \infty\)? 📝 Show Solution
12 Further Reading
- Le Gouic and Loubes (2017) — Existence and consistency theory for Wasserstein barycenters.
- Agueh and Carlier (2011) — Gaussian barycenters and the Bures–Wasserstein metric.
- Le Gouic et al. (2023) — Sharp convergence rates under convexity/smoothness of transport maps.
- Cuturi and Doucet (2014); Benamou et al. (2015) — Entropic regularization and Sinkhorn scaling.
- Álvarez-Esteban et al. (2016) — Fixed-point computation of Wasserstein barycenters.
- Panaretos and Zemel (2020) — Statistical aspects of optimal transport.
13 Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
References
Agueh, M., and G. Carlier. 2011. “Barycenters in the Wasserstein Space.” SIAM Journal on Mathematical Analysis 43 (2): 904–24.
Ahuja, Ravindra K., and James B. Orlin. 1992. “The Scaling Network Simplex Algorithm.” Operations Research 40 (1-supplement-1): S5–13. https://doi.org/10.1287/opre.40.1.S5.
Álvarez-Esteban, Pedro C., Eustasio del Barrio, Juan A. Cuesta-Albertos, and Carlos Matrán. 2016. “A Fixed-Point Approach to Barycenters in Wasserstein Space.” Journal of Mathematical Analysis and Applications 441 (2): 744–62. https://doi.org/10.1016/j.jmaa.2016.04.045.
Benamou, Jean-David, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré. 2015. “Iterative Bregman Projections for Regularized Transportation Problems.” SIAM Journal on Scientific Computing 37 (2): A1111–38. https://doi.org/10.1137/141000439.
Cuturi, Marco, and Arnaud Doucet. 2014. “Fast Computation of Wasserstein Barycenters.” Proceedings of the 31st International Conference on Machine Learning, Proceedings of machine learning research, vol. 32: 685–93.
Janati, Hicham, Marco Cuturi, and Alexandre Gramfort. 2020. “Debiased Sinkhorn Barycenters.” Proceedings of the 37th International Conference on Machine Learning, Proceedings of machine learning research, vol. 119: 4692–701. https://proceedings.mlr.press/v119/janati20a.html.
Le Gouic, Thierry, and Jean-Michel Loubes. 2017. “Existence and Consistency of Wasserstein Barycenters.” Probability Theory and Related Fields 168 (3–4): 901–17. https://doi.org/10.1007/s00440-016-0727-z.
Le Gouic, Thierry, Quentin Paris, Philippe Rigollet, and Alexander J. Stromme. 2023. “Fast Convergence of Empirical Barycenters in Alexandrov Spaces and the Wasserstein Space.” Journal of the European Mathematical Society 25 (6): 2229–50. https://doi.org/10.4171/JEMS/1234.
Nesterov, Yurii, and Arkadii Nemirovskii. 1994. Interior-Point Polynomial Algorithms in Convex Programming. Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9781611970791.
Panaretos, Victor M., and Yoav Zemel. 2020. An Invitation to Statistics in Wasserstein Space. SpringerBriefs in Probability and Mathematical Statistics. Springer. https://doi.org/10.1007/978-3-030-38438-8.
Petersen, A., and H.-G. Müller. 2016. “Functional Data Analysis for Density Functions by Transformation to a Hilbert Space.” The Annals of Statistics 44 (1): 183–218.
Peyré, Gabriel, and Marco Cuturi. 2019. “Computational Optimal Transport.” Foundations and Trends in Machine Learning 11 (5–6): 355–607. https://doi.org/10.1561/2200000073.
Schmitzer, Bernhard. 2019. “Stabilized Sparse Scaling Algorithms for Entropy Regularized Transport Problems.” SIAM Journal on Scientific Computing 41 (3): A1443–81. https://doi.org/10.1137/16M1106018.