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, 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, which is the topic of Sample Complexity.

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_0$, define

so that $\alpha_1=\alpha_0$. 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 and Vallender’s CDF formula 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_0$, $\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)$.

Empirical-Process Fluctuations¶

The fluctuation object used by statistical OT is the empirical law itself, not only a normalized sum in Euclidean space. There is no canonical Gaussian law on an unrestricted infinite-dimensional space of signed measures, so the precise finite-dimensional statement is made through test functions. Uniform convergence over an infinite class requires a Donsker condition Vaart & Wellner, 1996.

This linear CLT is the input to the delta-method arguments in Bias and Variance of OT. Exact OT may have only directional derivatives, whereas fixed-temperature entropic OT is smooth in the relative interiors of finite probability simplices.

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.

Figure Div gives a numerical overview: 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¶

Smooth densities permit a better estimator: first remove fine empirical oscillations with a wavelet or kernel smoother, then compute OT between the smoothed laws Niles-Weed & Berthet, 2022.

This gain is statistical rather than automatically computational: a grid at bandwidth $h$ has $O(h^{-d})$ degrees of freedom, so direct smoothing remains expensive in high dimension. Sum-of-squares relaxations offer theoretical dimension-free exponents but are still costly in practice Vacher et al., 2021.

MMD¶

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.

Sample Complexity of Estimating OT Maps¶

The preceding estimates concern transport values. For map estimation, one solves the empirical dual problem 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 rates require stronger stability assumptions on Brenier potentials. For unregularized map estimation, Hütter--Rigollet, Deb--Ghosal--Sen, and Manole--Balakrishnan--Niles-Weed--Wasserman obtain, under smoothness and curvature assumptions, squared $L^2(\alpha)$ risk of order

up to logarithms, and this rate is minimax optimal Hütter & Rigollet, 2021Deb et al., 2021Manole et al., 2024. For the Sinkhorn-computable entropic barycentric estimator, Pooladian--Niles-Weed prove

with $\bar s=s\wedge3$ Pooladian & Niles-Weed, 2021. The first term is the statistical fluctuation of the entropic map, while the last two terms are regularization-bias terms. Thus fixed $\epsilon$ gives a parametric empirical fluctuation for the regularized map, whereas $\epsilon\simeq n^{-1/(d+\bar s+1)}$ gives a convergent estimator of the Brenier map. The unregularized case can achieve sharper minimax rates when combined with enough smooth extension machinery; the regularized estimator is instead directly Sinkhorn-computable.

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.

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. Statistically, it inherits one-dimensional empirical rates rather than ambient-dimensional matching rates Nadjahi et al., 2019Nadjahi et al., 2020Manole et al., 2022.

Bias and Variance of OT¶

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. Rather than only estimating the size of the error, we seek the first expansion of the plug-in estimator. Given empirical laws $\hat\alpha_n$ and $\hat\beta_m$, and a scalar OT value $\mathcal V$, this expansion separates the leading statistical bias

and the centered fluctuation

for values such as $\mathcal V(\alpha,\beta)=\Wass_p^p(\alpha,\beta)$, $\mathcal V(\alpha,\beta)=\MK_c(\alpha,\beta)$, or their entropic analogues. This goes beyond Proposition Proposition: Finite-Dimensional Empirical-Process CLT: the empirical fluctuation is now pushed through a nonlinear, often nonsmooth, transport value. The resulting law is governed by local differentiability and by the geometry of the optimal dual face.

Literature Map¶

The shape of the limit depends sharply on the analytic regularity of the transport value. For exact OT on a finite or countable space, the map from empirical weights to optimal cost is convex and piecewise affine. Sommerfeld and Munk Sommerfeld & Munk, 2018, and the countable-space extensions of Tameling--Sommerfeld--Munk Tameling et al., 2019, show that the natural limit is a directional delta-method limit: often a support function of a Gaussian process over an optimal dual face, rather than a Gaussian random variable. For Euclidean costs, del Barrio--Loubes and collaborators obtain central-limit theorems for empirical transportation costs under regularity and uniqueness assumptions on the Kantorovich potentials Barrio & Loubes, 2019Barrio et al., 2021. Entropic regularization makes the problem smoother. For fixed $\epsilon>0$, the value is differentiable with respect to the marginals, and the limiting variance is read directly from the entropic dual potentials. This is the point of view developed by Bigot--Cazelles--Papadakis, Klatt--Tameling--Munk, Hundrieser--Klatt--Munk, and Mena--Niles-Weed Bigot et al., 2019Klatt et al., 2020Hundrieser et al., 2021Mena & Niles-Weed, 2019. More recent results clarify weak limits beyond smooth costs and second-order null limits for Sinkhorn divergences González-Sanz & Hundrieser, 2023Goldfeld et al., 2024. These asymptotic statements do not replace the dimension-dependent bounds of Sample Complexity; they explain what the leading random object is, what constants appear, and where the bias comes from.

Bias Versus Centered Fluctuation¶

The decomposition

is elementary, but in OT it is genuinely diagnostic because the two terms may live on different scales. For the self-distance $\mathcal V(\hat\alpha_n,\alpha)=\Wass_p^p(\hat\alpha_n,\alpha)$, the population value is zero and the deterministic bias $\mathbb E\Wass_p^p(\hat\alpha_n,\alpha)$ is the leading term. On a regular $d'$-dimensional support, the high-dimensional matching scale is $n^{-p/d'}$ when $d'>2p$, with a critical logarithmic correction when $d'=2p$. Below this threshold the rate is controlled by empirical-process fluctuations and can be faster; in one dimension, smooth positive densities give the familiar quantile-process asymptotics for $\Wass_p^p$, of order $n^{-p/2}$ for $p>1$, and in particular $n^{-1}$ for $\Wass_2^2$. Thus $p/d'$ is the high-dimensional matching exponent, not a universal rate formula Dudley, 1969Fournier & Guillin, 2015Weed & Bach, 2019. By contrast, when $\mathcal V(\alpha,\beta)$ is nonzero and the OT value is differentiable at $(\alpha,\beta)$, the centered fluctuation is often governed by an ordinary $n^{-1/2}$ central limit theorem. A ``slow OT statistic’’ can therefore have three distinct causes: a large empirical bias, a nonsmooth directional limit, or a regularization bias $\MK_c^\epsilon(\alpha,\beta)-\MK_c(\alpha,\beta)$ that still has to be removed.

The same argument gives the independent two-sample limit. If $N/(N+M)\to\lambda\in(0,1)$, then, with $r_{N,M}=\sqrt{NM/(N+M)}$,

where $G_a$ and $G_b$ are independent multinomial Gaussian limits. When both dual potentials are unique, the first-order bias is negligible and the asymptotic variance of the unscaled estimator is

If the optimal face contains several dual potentials, the limit is the supremum of a Gaussian process over that face and need not itself be Gaussian. The leading bias is then the expectation of this supremum divided by the square-root sample size. This nonsmoothness is the basic reason why inference for exact OT can remain delicate even on a finite space.

For the entropy-only convention $\MKD_\C^\epsilon$ used in Entropic Regularization for Discrete Measures,

The marginal entropy derivatives described after Proposition: First variations of entropic OT must therefore be added to the potentials. For the debiased Sinkhorn divergence $\bar\MK_c^\epsilon$, the derivative is instead the difference between cross and self entropic potentials. At the null $\alpha=\beta$, this first derivative can vanish for the debiased statistic, and second-order limits then become relevant Goldfeld et al., 2024. This is a useful warning: the asymptotic law is determined not only by a rate, but also by the local geometry of the functional.

Three Error Terms When Entropy Estimates Exact OT¶

If the target is $\MK_c(\alpha,\beta)$ but the statistic uses $\MK_c^\epsilon$, three errors should be kept separate:

Although algebraic, this identity is important statistically. It is usually read after choosing $\epsilon=\epsilon_n$. For fixed $\epsilon$, the statistical bias is typically $O(n^{-1}+m^{-1})$ in finite dimension, and the centered fluctuation is $O_{\mathbb P}(n^{-1/2}+m^{-1/2})$; boundary events where an empirical finite histogram has a zero entry have exponentially small probability when the population weights are positive. The regularization bias disappears only when $\epsilon\to0$, whereas the variance constants usually deteriorate as $\epsilon\downarrow0$. This is the asymptotic form of the bias--variance tradeoff illustrated in Div.

What Changes in Continuous Spaces¶

The finite-dimensional formulas above are not merely toy models; they are the cleanest template for the general mechanism. Whenever an OT value is Hadamard differentiable and admits a unique sufficiently regular dual potential $f^\star$, the empirical-process CLT and the delta method give a first-order Gaussian limit with one-sample variance $\operatorname{Var}_\alpha(f^\star(X))$. When differentiability fails, as can happen for exact OT because the dual optimizer is not unique or because one studies the degenerate self-distance $\MK_c(\alpha,\alpha)=0$, the limit may be non-Gaussian or may occur at a slower non-parametric scale. Entropic regularization smooths the dual potentials and typically restores a conventional first-order CLT for fixed $\epsilon$, even under weak assumptions on the cost. If the target is unregularized OT, however, this statistical expansion must still be balanced against the bias $\MK_c^\epsilon(\alpha,\beta)-\MK_c(\alpha,\beta)$. Thus Sample Complexity gives dimension-dependent magnitudes, while the present section identifies the local asymptotic constants and covariance formulas needed for inference.

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. For data points $(x_i)_{i=1}^n$, the dense Gram matrix $K=(k(x_i,x_j))_{ij}$ is replaced by $\widetilde K=\Phi\Phi^\top$, where $\Phi\in\RR^{n\times rp}$ is the sampled feature matrix with entries $\Phi_{i,(\ell,q)}=r^{-1/2}\phi_q(x_i,z_\ell)$, $1\leq \ell\leq r$, $1\leq q\leq p$. The number $r$ of sampled latent variables is thus a second sample size, distinct from the number of data points. 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.

Application to Sinkhorn Kernels¶

Let

Assume here that the weights are positive and that $K_{ij}>0$. The Sinkhorn scaling equations are

where divisions are componentwise. The factorization used by Sinkhorn is rectangular, because the source points $x_i$ and target points $y_j$ need not coincide. Starting from the same sketch $k(x,y)\simeq r^{-1}\sum_{\ell=1}^r \dotp{\phi(x,z_\ell)}{\phi(y,z_\ell)}_{\RR^p}$, set $R=rp$ and flatten the pair $(\ell,q)$ into a single column index. Define

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

When the two supports are the same, one has $\Phi_X=\Phi_Y=\Phi$ and this reduces to the symmetric Gram approximation $\widetilde K=\Phi\Phi^\top$. In the rectangular case, if $K$ is replaced by the rank-$R$ factorization $\widetilde K=\Phi_X\Phi_Y^\top$, 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.

The next proposition is not a substitute for the sharp low-rank Sinkhorn theory. It is an elementary concentration lemma whose role is to isolate the two facts that matter for Sinkhorn: positivity of the approximate kernel and accuracy after taking logarithms. The full algorithmic statements rely on more refined Nyström or Gaussian-kernel approximations together with stability of Sinkhorn scaling Altschuler et al., 2018Altschuler et al., 2019.

This proposition pinpoints the obstruction. When $\epsilon$ is small or the data diameter is large, $\kappa=\min K_{ij}$ can be exponentially small, so an ordinary signed sketch may need a very large $r$ just to keep all entries positive. What Sinkhorn really needs is relative control of the kernel, equivalently additive control of the sketched costs. Positive features do not solve the whole approximation problem, but they remove the most basic failure mode, negative entries in $\widetilde K$. The log-normal positive features below are unbounded, so this Hoeffding proposition does not apply to them without truncation or a different concentration argument.

Available worst-case bounds retain explicit dimension dependence. For the quadratic cost on points contained in a radius-$R$ ball of $\RR^d$, the Nyström analysis of Altschuler, Bach, Rudi and Niles-Weed Altschuler et al., 2019 studies $K_{ij}=\exp(-\eta\norm{x_i-x_j}^2)$ and bounds the effective rank by a quantity of the form

up to universal constants and harmless logarithms. Thus the resulting Sinkhorn computation is near-linear in $n$ once this rank is fixed, but the required rank is polynomial in the inverse kernel bandwidth, with exponent $d$ and only logarithmic dependence on $n$ in this bound. When the same machinery is used to approximate unregularized quadratic OT to additive accuracy $\tau$, one takes $\eta$ of order $\tau^{-1}\log n$, so the complexity contains a factor $\tau^{-O(d)}$, often summarized informally as an $\epsilon^{-d}$-type dependence on the target precision. The same work also explains why this can be much better on data with lower intrinsic dimension.

Positive Features and Total Positivity¶

A positive feature sketch asks for more than positive semidefiniteness. In finite dimension, let

Matrices in $\mathrm{DNN}_n$ are doubly nonnegative, while matrices in $\mathrm{CP}_n$ are completely positive. One always has $\mathrm{CP}_n\subset\mathrm{DNN}_n$, with equality for $n\leq4$ and strict inclusion for $n\geq5$ Berman & Shaked-Monderer, 2003Burer et al., 2009. Checking membership in the completely positive cone is computationally hard: strong membership is NP-hard, and weak membership for the completely positive cone and its copositive dual is NP-hard Dickinson & Gijben, 2014.

For kernels, it is useful to distinguish the two corresponding notions. A kernel $k$ is doubly positive if every finite Gram matrix belongs to $\mathrm{DNN}_n$, i.e. if $k$ is PSD and pointwise nonnegative. It is totally positive in the positive-feature sense if every finite Gram matrix belongs to $\mathrm{CP}_n$. This terminology should not be confused with Karlin’s classical total positivity, which concerns nonnegativity of all minors of an ordered kernel Karlin, 1968; it is a complete-positive-type condition, studied for instance in DeCorte et al., 2022.

The positive-feature construction of Scetbon and Cuturi Scetbon & Cuturi, 2020 makes the reason for this definition direct. If

then every finite Gram matrix is a conic combination of rank-one matrices $vv^\top$ with $v\geq0$, and hence is completely positive. Conversely, on a fixed finite set, $K\in\mathrm{CP}_n$ is exactly the existence of a nonnegative factorization $K=BB^\top$. Thus total positivity is the relevant finite-sample kernel condition behind positive random features. For rectangular source--target matrices, one applies the same statement to the Gram matrix on the union of the sampled source and target points and then extracts its $n\times m$ cross block.