Statistical Optimal Transport

Optimal transport is rarely evaluated on population measures directly. In machine learning and statistics, the inputs are usually empirical laws, histograms, simulated particles or minibatches, and the central question is therefore no longer only how to compute OT, but how OT behaves as a random estimator. This chapter studies this statistical layer: qualitative consistency of empirical measures, non-asymptotic sample-complexity rates, and asymptotic bias--variance decompositions for exact and regularized transport costs.

This statistical convergence is conceptually different from the algorithmic convergence studied in the previous chapter. There the marginals and the temperature were fixed and one asked how Sinkhorn iterates approach a regularized optimizer. Here the number of samples grows, the empirical measures themselves move, and the regularization parameter may either remain fixed or be sent to zero. The resulting picture explains why exact OT is statistically expensive in high intrinsic dimension, why fixed-temperature Sinkhorn has smoother parametric fluctuations, and why approximating exact OT with entropy always involves a bias--variance tradeoff.

from pathlib import Path
import sys

from IPython.display import Image as DisplayImage
from IPython.display import display

here = Path.cwd()
myst_dir = None
for candidate in [here, here.parent, here / "myst", here.parent / "myst", here.parent.parent / "myst"]:
    if (candidate / "ot4ml_web.py").exists():
        myst_dir = candidate.resolve()
        sys.path.insert(0, str(myst_dir))
        break

if myst_dir is None:
    raise RuntimeError("Could not locate myst/ot4ml_web.py")

repo_root = myst_dir.parent
thumbnails = repo_root / "notebooks-figures" / "thumbnails"

def show_book_figure(name, width=760):
    display(DisplayImage(filename=str(thumbnails / f"{name}.png"), width=width))

Law of Large Numbers and Central Limit Theorem¶

Before discussing sample complexity in Section Sample Complexity, it is useful to separate consistency from rates. If $X_1,\ldots,X_n$ are i.i.d. samples with common law $\alpha$, the associated empirical measure is the random probability measure

The ordinary law of large numbers says that empirical averages converge to expectations. In measure language this means that $\hat\alpha_n$ converges weakly toward $\alpha$, because testing $\hat\alpha_n$ against a bounded continuous function $\varphi$ gives the sample average $n^{-1}\sum_i\varphi(X_i)$. Wasserstein distances strengthen this statement by also recording moment convergence. Thus, if $\alpha$ has a finite $p$-th moment, the empirical law converges to $\alpha$ in $\Wass_p$, almost surely and in $p$-th mean in the sense that $\mathbb E\Wass_p(\hat\alpha_n,\alpha)^p\to0$. This is the qualitative consistency statement behind empirical OT plug-in estimators: empirical transport distances converge to their population counterparts once the sampled laws themselves converge in Wasserstein distance. It says nothing yet about the speed.

Central-Limit Fluctuations¶

The previous proposition is a law-of-large-numbers statement: a random empirical measure converges to the law that generated it. The central limit theorem describes a different, fluctuation-scale limit. As recalled in Remark Remark: Central limit theorem, if $(X_i)_{i\geq1}$ are centered i.i.d. random vectors with identity covariance, the law of $n^{-1/2}\sum_i X_i$ converges weakly toward a Gaussian. Equivalently, if $\alpha$ is the common law of the $X_i$, this law is the rescaled convolution $(D_{1/\sqrt n})_\sharp\alpha^{*n}$. Wasserstein distances make this qualitative convergence quantitative. The next result is a $\Wass_1$ form of the Berry--Esseen theorem: it controls the error uniformly over all 1-Lipschitz test functions. The qualitative Bernoulli example was visualized earlier in Figure Div.

The universal $n^{-1/2}$ bound is sharp over broad classes of input laws, but it need not describe the asymptotic behavior of a fixed law. For symmetric inputs, the smooth skewness correction vanishes. What replaces it depends on whether the law is lattice-valued or has a density.

The absence of a lattice correction reveals the next smooth Edgeworth term. A density automatically satisfies Cramér’s non-lattice condition because its characteristic function vanishes at infinity.

These two results explain the different behaviors in Figure Div. For the symmetric Bernoulli law, the maximal span is $h=2$, hence $\Wass_1(\alpha_n,\gamma)=1/(2\sqrt n)+o(n^{-1/2})$. The density proposition cannot be used: this law has no density and its characteristic function is $\mathbb E[e^{itX}]=\cos t$, whose modulus returns to one at arbitrarily large frequencies. Thus symmetry removes the smooth skewness term but not the lattice sawtooth. Conversely, $X\sim\operatorname{Unif}[-\sqrt3,\sqrt3]$ has a density, characteristic function $\sin(\sqrt3t)/(\sqrt3t)\to0$, and $\kappa_4=-6/5$. Since

one obtains

Figure Div confronts these equivalents with exact one-dimensional computations. For a centered, unit-variance input law $\alpha$, define

so that $\alpha_1=\alpha$. For Bernoulli input, the distance is evaluated by the exact quantile formula and atom masses are divided by the current lattice spacing in the density display. For continuous-uniform input, the normalized convolution is an affine image of the Irwin--Hall distribution, whose density is, up to rescaling, a cardinal B-spline of degree $n-1$, and the absolute CDF difference is integrated numerically. Neither computation uses Monte Carlo sampling.

show_book_figure("statistical-berry-esseen-w1", width=920)

Sharp lattice and density central-limit asymptotics in $\Wass_1$. The two left panels show $\alpha_1$, $\alpha_2$, and $\alpha_6$ for symmetric Bernoulli and continuous-uniform inputs; the gray curve is the standard Gaussian density. The right panel compares the exact numerical distances (solid) with the sharp equivalents in (17) (dashed). The Bernoulli curve follows its lattice rate $1/(2\sqrt n)$, while the continuous-uniform curve approaches $(1+4e^{-3/2})/(10\sqrt{2\pi}\,n)$.

Sample Complexity¶

This section compares four statistical regimes. Exact OT resolves geometry at all scales and pays rates controlled by the intrinsic dimension. MMD replaces transport by a kernel mean embedding and has a parametric Monte-Carlo rate for bounded kernels. Fixed-temperature Sinkhorn divergences smooth the dual potentials and also recover parametric fluctuations, at the price of regularization bias. Sliced Wasserstein inherits one-dimensional empirical rates through its projected laws.

Unregularized OT¶

The previous section proves qualitative convergence of empirical laws. The sample-complexity question is how fast this convergence happens. The sample complexity of unregularized OT suffers from the curse of dimensionality, but the relevant dimension is geometric. If the distributions are supported on a regular lower-dimensional set, for instance a $d'$-dimensional submanifold of $\RR^d$, the empirical rate is governed by $d'$ rather than by the ambient dimension $d$. In the high-dimensional regime for $\Wass_p$, namely $d'>2p$, this gives the characteristic rate $n^{-1/d'}$. Exact OT is therefore dimension-adaptive: it sees the intrinsic dimension of the data support through its covering numbers Dudley, 1969Weed & Bach, 2019. Related two-sample-testing viewpoints are developed in Ramdas et al., 2017.

For $n\geq2$, $p\geq1$, and $d\geq1$, define the worst-case empirical scale

The following proposition turns this scale into uniform one-sample and two-sample bounds for empirical OT.

The following multiscale estimate makes the mechanism behind the $\Wass_1$ bound explicit.

Figure Div gives a numerical overview of these regimes: exact OT exhibits dimension-dependent empirical fluctuations, whereas MMD and fixed-temperature Sinkhorn divergences lie much closer to the parametric $n^{-1/2}$ scale.

show_book_figure("sinkhorn-bias-variance-tradeoff")

Empirical fluctuations in dimensions three and six. For each sample size $n$, two independent empirical measures are drawn from the same standard Gaussian law. Exact OT follows a slower dimension-dependent scale, while MMD and the fixed-$\epsilon$ Sinkhorn divergence behave closer to the parametric $n^{-1/2}$ guide. This is a statistical illustration, not a solver benchmark.

Lower Bounds and Minimax Optimality¶

The high-dimensional rate is not merely a defect of the empirical plug-in estimator. A hypercube of localized density perturbations gives the following distribution-estimation obstruction Singh & Póczos, 2018Weed & Bach, 2019Chewi et al., 2025.

Leveraging Smoothness¶

For a Sobolev-smooth density, convolving the empirical measure at bandwidth $h$ creates a bias of order $h^{s+1}$ but reduces the empirical fluctuation to order $n^{-1/2}h^{1-d/2}$. Balancing these terms selects $h_n\asymp n^{-1/(d+2s)}$, which improves on the empirical $n^{-1/d}$ rate when $d>2$ Niles-Weed & Berthet, 2022Divol, 2022.

The restriction $s\leq1$ permits a nonnegative symmetric kernel and hence a genuine probability estimator. Exploiting higher smoothness requires higher-order, generally signed kernels followed by a positivity correction. Dimensions one and two have separate parametric or logarithmic regimes Divol, 2022.

This statistical acceleration has a numerical price. For a Gaussian kernel, $\tilde\alpha_n$ is a mixture of $n$ Gaussians, and OT between two such mixtures has no closed form in general. A direct particle implementation draws $M$ independent vectors $\xi_{i,j}\sim\mathcal N(0,\mathrm{Id}_d)$ around each observation and defines

Here addition is understood modulo $\mathbb Z^d$. Conditionally on the observations, convexity of $W_1$ under mixing, together with translation, scaling and the 1-Lipschitz projection onto the torus, gives

The last estimate is the moment form of the empirical OT rate Fournier & Guillin, 2015. At the optimal bandwidth $h_n\asymp n^{-1/(d+2s)}$, keeping this extra error below the statistical rate is guaranteed by $M\gtrsim n^{sd/(d+2s)}$. This componentwise strategy therefore uses $N_{\rm comp}=nM\gtrsim n^{1+sd/(d+2s)}$ particles per measure.

Approximating every Gaussian separately is conservative. One can instead sample or quantize the whole smoothed mixture using $N$ points. In dimension $d>2$, the corresponding $W_1$ error is of order $N^{-1/d}$, up to dimension-dependent constants Graf & Luschgy, 2000. Thus target accuracy $\delta$ requires $N\gtrsim\delta^{-d}$, whereas the componentwise construction uses $N_{\rm comp}\gtrsim\delta^{-d-2s/(s+1)}$. The global strategy is cheaper, but its particle count remains exponential in $d$ at fixed accuracy. Forming and storing a dense $N\times N$ transport kernel costs $O(N^2)$ work and memory, and each direct Sinkhorn scaling costs another $O(N^2)$ operations. Fast kernel summation or low-rank approximations may reduce this matrix cost, but without additional structure they do not remove the underlying $N^{-1/d}$ discretization barrier.

MMD¶

MMD contains no transport optimization: its square is a combination of expectations of kernel evaluations, estimated empirically by sums of kernel values. Ordinary Monte Carlo averaging therefore gives a dimension-free parametric rate for bounded kernels. This sharply contrasts with exact empirical OT, whose value is defined through an optimization problem and whose worst-case $W_1$ plug-in bias in dimension $d>2$ decays only as $n^{-1/d}$, corresponding to a sample requirement of order $\delta^{-d}$ at accuracy $\delta$. The constants for MMD still depend on the kernel and its bandwidth, and changing the kernel changes the discrepancy being estimated.

Entropic OT¶

Entropic regularization keeps the transport geometry while smoothing the dual potentials. The sharp fixed-temperature statement below has polynomial, rather than exponential, dependence on the inverse temperature Genevay et al., 2019Mena & Niles-Weed, 2019.

For fixed $\epsilon>0$, the empirical fluctuation is therefore parametric, while the prefactor deteriorates only polynomially as $\epsilon\downarrow0$. Compact support is a special case: if both supports lie in a common ball $B(z,R)$, translating them by $-z$ leaves the quadratic transport problem unchanged, and the subgaussian assumption holds with $\sigma=R/\sqrt{2d\log 2}=O(R)$. Genevay--Chizat--Bach--Cuturi--Peyré established the compactly supported result, and Mena--Niles-Weed extended it to the subgaussian setting above Genevay et al., 2019Mena & Niles-Weed, 2019.

The interactive demo below is only a scaling guide: change the dimension to see the exact-OT exponent flatten, and change $\epsilon$ to move the Sinkhorn bias floor against its parametric fluctuation term.

Sample Complexity of Estimating OT Maps¶

The preceding estimates concern transport values. For map estimation, one solves the empirical discrete dual problem (2) or its entropic counterpart (47) for potentials $(\mathbf f_{n,m}^\epsilon,\mathbf g_{n,m}^\epsilon)\in\RR^n\times\RR^m$, with $\epsilon=0$ denoting the unregularized Kantorovich problem, and then uses the soft $c$-transform as an out-of-sample extrapolator.

For the quadratic cost $c_2(x,y)=\norm{x-y}^2/2$, let $\P_{n,m}^\epsilon=(\P_{ij}^\epsilon)$ be the empirical entropic coupling and let $\mathbf g_{n,m}^\epsilon$ be the target-side dual potential. Define

Equivalently, subtracting the quadratic marginal terms turns $c_2$ into the bilinear cost $c_{\rm ip}(x,y)=-\langle x,y\rangle$. With the sign convention used here, the associated soft convex potential is

Differentiating this log-sum-exp gives the explicit soft barycentric form

At a sampled point $X_i$, Sinkhorn normalization gives

which is the barycentric projection of the entropic plan. At $\epsilon=0$, the hard $c$-transform is generally nondifferentiable: away from Laguerre boundaries its gradient selects one target site, whereas the barycentric projection of an optimal unregularized plan may average several target sites when a row splits. The consistency statement is therefore formulated directly for barycentric projections. It uses uniform empirical weights; deterministic nonuniform weights require the corresponding entropy bound.

Quantitative map convergence requires stronger regularity than consistency. The following result gives a direct guarantee for the Sinkhorn-based estimator defined above Pooladian & Niles-Weed, 2021.

This estimator is directly computable from Sinkhorn potentials. Under stronger modeling and smoothing assumptions, unregularized estimators can attain sharper minimax rates Hütter & Rigollet, 2021Deb et al., 2021Manole et al., 2024.

Sliced Wasserstein¶

For $P_\theta(x)=\langle\theta,x\rangle$ and normalized surface measure $\sigma$ on $\mathbb S^{d-1}$, define

The construction is studied in Sliced Wasserstein Distances. For $p=1$, it inherits a particularly clean one-dimensional empirical-CDF rate rather than an ambient-dimensional matching rate Nadjahi et al., 2019Nadjahi et al., 2020Manole et al., 2022.

Bias and Variance of OT¶

Bias, Variance and Approximation Errors¶

The previous section answered a coarse but essential question: how many samples are needed before an empirical OT quantity is accurate, typically up to universal constants. We now ask for a finer statistical description of the KL-normalized entropic OT value $\MK_c^\epsilon$ defined in (20), with the convention $\MK_c^0\eqdef\MK_c$ for unregularized OT. Fixing $c$, $\alpha$ and $\beta$, define the statistical bias

and the centered fluctuation

If the target is the exact cost but the statistic is entropically regularized, these quantities isolate the three sources of error:

For the debiased Sinkhorn divergence, define $\bar B_{n,m}^\epsilon$ and $\bar Z_{n,m}^\epsilon$ by replacing $\MK_c^\epsilon$ with $\bar\MK_c^\epsilon$ above; the same decomposition then holds with bars. It is usually read after choosing a temperature $\epsilon=\epsilon_n$. At fixed $\epsilon>0$ and on finite supports, the statistical bias is typically $O(n^{-1}+m^{-1})$, while the centered fluctuation is $O_{\mathbb P}(n^{-1/2}+m^{-1/2})$. The regularization bias disappears only when $\epsilon\to0$, whereas the fluctuation constants can deteriorate in this limit. This is the bias--variance tradeoff illustrated in Div.

The remainder of this section analyzes the latter two terms. The superscript records their dependence on the regularization strength. The empirical fluctuation is pushed through the nonlinear transport value $\MK_c^\epsilon$, which can be nonsmooth at $\epsilon=0$. Its limiting law is governed by local differentiability and by the geometry of the optimal dual potentials.

For exact OT on finite spaces, empirical central-limit theorems follow from a directional delta method Sommerfeld & Munk, 2018Tameling et al., 2019; fixed positive entropic regularization makes the transport value smooth on the interiors of the marginal simplices Bigot et al., 2019Klatt et al., 2020Hundrieser et al., 2021Mena & Niles-Weed, 2019. The following statement isolates the Gaussian regime common to both settings: at zero temperature, uniqueness of the dual potentials turns the directional derivative into an ordinary linear derivative.

The uniqueness assumption at $\epsilon=0$ is substantive. Without it, the directional delta method gives the generally non-Gaussian limit

where $\mathcal D_c(\alpha,\beta)$ is the set of optimal dual pairs modulo the additive gauge. Its expectation need not vanish, so the exact-OT bias can then have order $n^{-1/2}$ rather than being negligible at the central-limit scale Sommerfeld & Munk, 2018Tameling et al., 2019.

Fixed Support Versus a Continuum¶

Proposition: Finite-Space Bias and CLT for OT keeps the support sizes $N$ and $M$ fixed while $n$ grows. More explicitly, its limiting variable is

The variance does not necessarily diverge with $N,M$: for a uniformly bounded cost, normalized dual potentials have bounded oscillation, hence $\mathsf{v}_\epsilon$ remains bounded. What fails to be uniform is the differentiability argument. If $a_{\min}=\min_i a_i$ and $b_{\min}=\min_j b_j$, then

For nearly uniform weights, making this probability vanish already requires both $n\gg N\log N$ and $n\gg M\log M$. Moreover, exact-dual uniqueness margins can close as the supports become dense, while entropic Hessian bounds can deteriorate with the smallest weights and as $\epsilon\downarrow0$. Thus the constant hidden in the bias estimate may grow even when the first-order variance stays bounded.

There is therefore no distribution-free passage from this multinomial CLT to arbitrary continuous laws. Continuous CLTs nevertheless hold under hypotheses ensuring stable unique potentials, with limiting variance $\operatorname{Var}_\alpha(f_\epsilon^\star)+ \operatorname{Var}_\beta(g_\epsilon^\star)$ in the independent two-sample setting Barrio & Loubes, 2019Barrio et al., 2021Mena & Niles-Weed, 2019González-Sanz & Hundrieser, 2023. The first-order description ceases to be informative in degenerate regimes. For instance, when $\epsilon=0$, $c=\dist^p$, and $\alpha=\beta$, the first derivative vanishes; the two-sample empirical cost $\MK_{\dist^p}(\hat\alpha_n,\hat\beta_n)$ is instead governed by the dimension-dependent matching behavior studied in Sample Complexity. This is where the curse of dimensionality re-enters. The first derivative of the debiased Sinkhorn divergence also vanishes under the null, so its nontrivial limit is of second order Goldfeld et al., 2024.

Sketching Sinkhorn in Linear Time¶

The statistical results above explain why entropic OT is attractive at fixed temperature, but they also expose a computational tension. Sampling the marginals replaces the input laws by empirical measures and reduces the problem to finitely many support points; subsampling or coarsening these supports is governed by the sample-complexity estimates of this chapter. Even after the supports are fixed, however, a Sinkhorn sweep still costs $O(nm)$, because it repeatedly applies the Gibbs kernel to vectors. Kernel sketching attacks this second bottleneck: it keeps the sampled marginals fixed and samples, or otherwise compresses, the feature representation of the kernel used by Sinkhorn. With a rank $R$ factorization, $N_{\rm it}$ scaling sweeps cost $O(N_{\rm it}(n+m)R)$, which is linear in the support sizes only when both the required rank and the iteration count are controlled. The rank generally depends on dimension, temperature, and target accuracy.

PSD Kernels and Ordinary Feature Sketches¶

The RKHS/MMD section, Dual RKHS Norms and Maximum Mean Discrepancies, already used positive semidefinite kernels to define Hilbertian discrepancies between measures. Recall that a symmetric kernel $k:\Xx\times\Xx\to\RR$ is positive semidefinite if, for every finite family $(x_i)_{i=1}^n$, the Gram matrix $(k(x_i,x_j))_{ij}$ is positive semidefinite. Many scalable kernel methods start from an integral feature representation

where $\xi$ is a latent distribution. Drawing $z_1,\ldots,z_r\sim\xi$ gives

This is the classical random-feature principle. Given two point families $(x_i)_{i=1}^n$ and $(y_j)_{j=1}^m$, define the rectangular cross-kernel matrix $K\in\RR^{n\times m}$ by $K_{ij}=k(x_i,y_j)$. Set $R=rp$, flatten $(\ell,q)$ into a single feature index, and define

Here $\Phi_X\in\RR^{n\times R}$, $\Phi_Y\in\RR^{m\times R}$, and

Although $K$ is generally rectangular and hence is not itself a Gram matrix, it is a cross-block of the Gram matrix evaluated on the union of the two point families. When the families coincide, $\Phi_X=\Phi_Y=\Phi$, and the construction reduces to the symmetric approximation $\widetilde K=\Phi\Phi^\top$. The number $r$ of sampled latent variables is thus a second sample size, distinct from $n$ and $m$. At this stage no positivity of $\phi$ is assumed; signed or oscillatory features are exactly what make Fourier sketches useful.

For translation-invariant kernels on $\RR^d$, Bochner’s theorem gives the standard Fourier sketch. If $k(x,y)=\kappa(x-y)$ and $\kappa(t)=\int e^{i\dotp{\omega}{t}}\,d\Lambda(\omega)$, with $\Lambda$ a positive spectral measure, one may use real features $\phi(x,\omega)=(\cos\dotp{\omega}{x},\sin\dotp{\omega}{x})$. For the Gaussian kernel $k_\sigma(x,y)=\exp(-\norm{x-y}^2/(2\sigma^2))$, this corresponds to $\omega\sim\mathcal N(0,\sigma^{-2}\Id)$. This random Fourier feature construction was introduced by Rahimi and Recht to accelerate large-scale kernel machines Rahimi & Recht, 2007; its statistical role for kernel ridge regression and related supervised methods is analyzed, for instance, in Rudi & Rosasco, 2017Avron et al., 2017.

The following concentration bound assesses the entrywise quality of the rectangular sketch (69). This is the relevant control when the kernel will subsequently be used inside logarithms or componentwise divisions.

Application to Sinkhorn Kernels¶

Let

Specialize the preceding construction to the Gibbs kernel $k_\epsilon(x,y)=e^{-c(x,y)/\epsilon}$, assumed to be PSD and to admit the displayed feature representation. Thus the rectangular matrix defined above is precisely the Sinkhorn kernel, $K_{ij}=k_\epsilon(x_i,y_j)$. Assume that the weights and the entries of $K$ are positive. Recall the Sinkhorn scaling equations (8):

where divisions are componentwise. Replacing $K$ by its rank-$R$ sketch (69), the two matrix-vector products are evaluated as

so each scaling sweep costs $O((n+m)R)$ operations and stores only $O((n+m)R)$ numbers. The difficulty is that Sinkhorn is not a generic kernel method: it divides by $\widetilde K v$ and, in dual variables, applies logarithms. Thus one needs $\widetilde K_{ij}>0$, preferably with a relative error on the entries, because $-\epsilon\log K_{ij}$ is the effective cost. This is a more stringent requirement than the usual spectral or Frobenius approximation of a PSD Gram matrix.