Suppose the input of a learning problem is a whole distribution: a population, a point cloud, or a collection of tokens of arbitrary size. We would like to approximate a continuous functional \(F(\mu)\) of that distribution by something finite. There is a remarkably direct recipe: draw \(n\) points from \(\mu\), evaluate \(F\) on their empirical distribution, and average the answer.
This recipe produces a polynomial of the measure. As \(n\) grows, it approximates every continuous functional, uniformly over all distributions on a compact space. Sampling error becomes approximation error. Bernstein’s probabilistic proof of the Weierstrass theorem [1] already contains the idea; its extension to probability measures belongs to the Bernstein–Schnabl construction [2].
The companion post Polynomials of probability measures explains the Stone–Weierstrass argument and its links to neural networks. Here the emphasis is constructive: which polynomial should we choose, why does it converge uniformly, and what determines its rate?
A polynomial whose variable is a distribution
Let \((K,d)\) be a nonempty compact metric space. Write \(\mathcal P(K)\) for its Borel probability measures. For \(1\le p<\infty\), the Wasserstein distance is
\[ W_p(\mu,\nu)^p=\inf_{\pi\in\Pi(\mu,\nu)} \int d(x,y)^p\,d\pi(x,y), \]where \(\Pi(\mu,\nu)\) is the set of couplings with the prescribed marginals. Compactness makes every finite \(W_p\) induce the weak topology, and makes \(\mathcal P(K)\) compact.
A continuous interaction kernel \(\varphi:K^m\to\mathbb R\) gives a polynomial of order at most \(m\):
\[ P_\varphi(\mu)=\int_{K^m}\varphi(x_1,\ldots,x_m)\, d\mu(x_1)\cdots d\mu(x_m). \]The order counts particles, not powers of their coordinates. For instance, the mean is order one and a pair interaction is order two. Along a mixture \(\mu_t=(1-t)\mu+t\nu\), expanding the \(m\) factors shows that \(P_\varphi(\mu_t)\) is an ordinary polynomial in \(t\), of degree at most \(m\). This is a statement about mixtures, not about Wasserstein geodesics.
Products stay in the class: multiply the two kernels on separate groups of variables. Sums do too: add unused variables to a lower-order kernel, then add the kernels. Constants belong to it, and order-one observables \(\int\phi\,d\mu\) separate measures. These functions form an algebra, so Stone–Weierstrass proves density, as in Cardaliaguet’s mean field game notes [4]. The following operator gives a particular approximating sequence.
Sample, evaluate, average
For independent \(X_1,\ldots,X_n\sim\mu\), define
\[ \widehat\mu_n=\frac1n\sum_{i=1}^n\delta_{X_i}, \qquad B_nF(\mu)=\mathbb E\,F(\widehat\mu_n). \]This is a polynomial of order at most \(n\), with the symmetric kernel
\[ \varphi_n(x_1,\ldots,x_n) =F\!\left(\frac1n\sum_{i=1}^n\delta_{x_i}\right). \]The kernel is continuous: pairing corresponding points gives
\[ W_p\!\left(\frac1n\sum_i\delta_{x_i}, \frac1n\sum_i\delta_{y_i}\right)^p \le\frac1n\sum_i d(x_i,y_i)^p. \]Moreover, \(B_n\) is linear in \(F\), preserves constants and positivity, and satisfies \(\|B_nF\|_\infty\le\|F\|_\infty\). Averaging cannot amplify a uniform error.
It also reproduces every linear observable exactly: if \(F(\mu)=\int\phi\,d\mu\), then \(B_nF(\mu)=\int\phi\,d\mu\), because the sample average is unbiased. Approximation enters when \(F\) depends nonlinearly on the distribution.
Theorem 1 (uniform constructive approximation). If \(K\) is compact and \(F:\mathcal P(K)\to\mathbb R\) is continuous, then
\[ \sup_{\mu\in\mathcal P(K)}|B_nF(\mu)-F(\mu)|\longrightarrow0. \]The supremum is the essential point: one particle count works for every input law, including discrete, singular, and continuous laws.
Proof: reduce the whole space to a finite histogram
Fix \(\delta>0\). Choose a finite \(\delta\)-net with \(M\) points, and a measurable nearest-point map \(q\). Set \(\mu^\delta=q_\#\mu\) and \(\widehat\mu_n^\delta=q_\#\widehat\mu_n\), using the same samples. Each original law is within \(\delta\) of its quantized version in \(W_p\).
Let \(w_j\) and \(\widehat w_j\) be the masses of the true and empirical histograms. Their expected mismatch satisfies
\[ \mathbb E\sum_{j=1}^M|\widehat w_j-w_j| \le\sum_{j=1}^M\sqrt{\frac{w_j(1-w_j)}n} \le\sqrt{\frac Mn}. \]Keep the common mass in place and transport the unmatched mass, which is half this \(\ell^1\) difference, over distance at most \(D=\operatorname{diam}K\). Minkowski’s inequality and the triangle inequality give
\[ \sup_\mu\bigl(\mathbb E W_p(\widehat\mu_n,\mu)^p\bigr)^{1/p} \le 2\delta+D\left(\frac12\sqrt{\frac Mn}\right)^{1/p}. \]First let \(n\to\infty\), then \(\delta\to0\). Thus the expected \(p\)-th transport error tends to zero uniformly in \(\mu\).
Finally \(F\) is bounded and uniformly continuous on the compact space \(\mathcal P(K)\). If \(\omega_F(\varepsilon)\) is its modulus of continuity, splitting according to whether \(W_p(\widehat\mu_n,\mu)\le\varepsilon\) gives
\[ |B_nF(\mu)-F(\mu)| \le\omega_F(\varepsilon) +\frac{2\|F\|_\infty}{\varepsilon^p} \mathbb E W_p(\widehat\mu_n,\mu)^p. \]Let \(n\to\infty\), then \(\varepsilon\to0\). \(\square\)
How dimension enters the rate
Suppose now \(K\subset\mathbb R^d\) is compact and \(F\) is Hölder in \(W_p\): for some \(L\ge0\) and \(0<\eta\le1\),
\[ |F(\mu)-F(\nu)|\le L\,W_p(\mu,\nu)^\eta. \]Theorem 2 (a uniform sampling rate). There is a constant \(C\), depending only on a bounding cube for \(K\), \(d,p,\eta\), such that
\[ \|B_nF-F\|_\infty\le CL\,r_n^\eta, \qquad r_n= \begin{cases} n^{-1/(2p)},&d<2p,\\ n^{-1/(2p)}(\log(1+n))^{1/p},&d=2p,\\ n^{-1/d},&d>2p. \end{cases} \]These are general upper bounds over all laws on \(K\), of the type developed by Fournier and Guillin [5]. A particular law or functional can converge faster. The critical logarithm is part of this bound, not a claim of optimality for every setting.
Proof: Hölder continuity followed by a multiscale histogram estimate
Jensen’s inequality immediately reduces the problem to empirical transport:
\[ |B_nF(\mu)-F(\mu)| \le L\,\mathbb E W_p(\widehat\mu_n,\mu)^\eta \le L\bigl(\mathbb E W_p(\widehat\mu_n,\mu)^p\bigr)^{\eta/p}. \]After rescaling, work in \([0,1]^d\), partitioned at level \(j\) into \(2^{dj}\) dyadic cubes \(\mathcal D_j\). Boundary points are assigned consistently to one cube.
Lemma 1 (transport by scales). For probability laws \(\mu,\nu\) on the cube,
\[ W_p(\mu,\nu)^p \le C_{d,p}\sum_{j\ge1}2^{-pj} \sum_{Q\in\mathcal D_j}|\mu(Q)-\nu(Q)|. \]Proof of Lemma 1: balance mass from the leaves to the root
First replace each measure by its masses at the centers of level-\(J\) cubes. Within each parent cube, match as much of the surplus and deficit of its children as possible, then pass the remaining imbalance to the parent. Continue to the root, whose total imbalance is zero. This constructs a coupling.
The amount matched between children at level \(j\) is at most \(\frac12\sum_{Q\in\mathcal D_j}|\mu(Q)-\nu(Q)|\). Every such pair lies in one parent cube, so its distance is at most \(\sqrt d\,2^{-(j-1)}\). Charge its \(p\)-th power to this level and sum. The bound is independent of \(J\). The center approximations converge to the original laws in \(W_p\), so letting \(J\to\infty\) proves the claim. \(\square\)
For the empirical histogram, the variance calculation from the previous proof gives
\[ \mathbb E\sum_{Q\in\mathcal D_j} |\widehat\mu_n(Q)-\mu(Q)| \le \min\!\left\{2,\sqrt{\frac{2^{dj}}n}\right\}. \]Thus
\[ \mathbb E W_p(\widehat\mu_n,\mu)^p \le C_{d,p}\sum_{j\ge1}2^{-pj} \min\{1,2^{dj/2}n^{-1/2}\}. \]Split at \(J=\lfloor(\log_2 n)/d\rfloor\). The coarse-scale sum is \(n^{-1/2}\sum_{j\le J}2^{j(d/2-p)}\); the fine-scale tail is at most \(C2^{-pJ}\). The coarse sum is bounded, grows like \(J\), or is dominated by its last term according as \(d<2p\), \(d=2p\), or \(d>2p\). These are exactly the three displayed rates for \(\mathbb E W_p^p\), followed by the exponent \(\eta/p\). \(\square\)
The threshold \(d=2p\) compares two effects: the \(n^{-1/2}\) fluctuations of the mass in a coarse cell, and the spatial resolution needed to locate many particles in high dimension. Finer results can exploit the intrinsic geometry of the law rather than just its ambient dimension [6].
The low-dimensional exponent can already be attained on a two-point space. Put \(K=\{0,1\}\), \(\mu_*=(\delta_0+\delta_1)/2\), and \(F(\mu)=W_p(\mu,\mu_*)^\eta\). This is an \(\eta\)-Hölder functional with constant one. At \(\mu_*\), its Bernstein error is of order \(n^{-\eta/(2p)}\), exactly the first rate in Theorem 2.
Why two atoms can attain this rate
Write \(\mu_a=(1-a)\delta_0+a\delta_1\). Unmatched mass travels distance one, so \(W_p(\mu_a,\mu_b)^p=|a-b|\). If \(K_n\) has binomial law with parameters \(n,1/2\), then
\[ B_nF(\mu_*)-F(\mu_*) =\mathbb E\left|\frac{K_n}{n}-\frac12\right|^{\eta/p}. \]The central limit theorem gives \((K_n-n/2)/\sqrt n\Rightarrow Z/2\), where \(Z\) is standard Gaussian. These variables have second moment \(1/4\), so their powers of order \(\eta/p\le1\) are uniformly integrable. Consequently,
\[ n^{\eta/(2p)}B_nF(\mu_*) \longrightarrow\mathbb E|Z/2|^{\eta/p}>0. \]The obstruction is random mass imbalance, not spatial dimension. This concerns the Bernstein sampling operator; it is not a lower bound for every possible polynomial approximation.
The original Bernstein formula reappears
Take \(K=\{0,1\}\), \(\mu_a=(1-a)\delta_0+a\delta_1\), and write \(f(a)=F(\mu_a)\). The number of ones in the sample has binomial law, so
\[ B_nF(\mu_a)= \sum_{k=0}^n f(k/n)\binom nk a^k(1-a)^{n-k}. \]This is precisely the classical Bernstein polynomial. It averages nearby values of \(f\), with a window whose width is \(\sqrt{a(1-a)/n}\).
Three examples make the mechanism visible. For \(f(a)=a\), the approximation is exact. For \(f(a)=a^2\),
\[ B_nf(a)=a^2+\frac{a(1-a)}n. \]More generally, if \(F(\mu)=(\int x\,d\mu)^2\) on a bounded real interval, its bias is exactly \(\operatorname{Var}(\mu)/n\). After subtracting the square of the mean, independence leaves only the diagonal variance terms.
For the corner \(f(a)=|a-\tfrac12|\), the approximation rounds the tip. At \(a=\tfrac12\), for even \(n=2m\), a binomial-sum identity gives
\[ B_{2m}f(1/2) =\frac12\,\frac{\binom{2m}{m}}{4^m} \sim \frac1{\sqrt{2\pi n}}. \]Thus the \(n^{-1/2}\) fluctuation scale is genuinely visible at a nonsmooth observable, even on a two-point space. It is the error of this particular operator; it does not rule out better polynomial approximations chosen by another method.
Numerical experimentation
A polynomial recipe is not a complexity bound
The construction is explicit, but the kernel \(\varphi_n\) lives on \(K^n\), a space of dimension \(dn\) in Euclidean coordinates. Representing a completely general kernel can be expensive. The rate in \(n\) is therefore a rate in interaction order, not automatically a rate in the number of network parameters or operations.
There are two separate approximations if one evaluates \(B_nF\) by averaging \(R\) independent empirical clouds. The bias is \(B_nF-F\), studied above. The additional Monte Carlo fluctuation has root-mean-square size at most \(\|F\|_\infty/\sqrt R\). Increasing \(R\) reduces the second error, but does not remove the first.
Useful architectures exploit additional structure: low-order interactions, a small family of features, or compositional rules acting on the law. This connects with Stochastic Deep Networks [7] and the pooling constructions in the companion post. The Bernstein operator gives a clean baseline against which those more economical representations can be understood.
References
- Sergei N. Bernstein. Démonstration du théorème de Weierstrass fondée sur le calcul des probabilités. Communications of the Kharkov Mathematical Society, second series, 13(1), 1–2, 1912.
- Roman Schnabl. Zur Approximation durch Bernsteinpolynome auf gewissen Räumen von Wahrscheinlichkeitsmaßen. Mathematische Annalen 180, 326–330, 1969.
- Gabriel Peyré. Optimal Transport for Machine Learners, 2025–2026. §4.4, Proposition 4.36, and Chapter 10.
- Pierre Cardaliaguet. Notes on Mean Field Games, from P.-L. Lions’s lectures at the Collège de France, 2012. §5.3, Proposition 5.9.
- Nicolas Fournier and Arnaud Guillin. On the rate of convergence in Wasserstein distance of the empirical measure. Probability Theory and Related Fields 162, 707–738, 2015.
- Jonathan Weed and Francis Bach. Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance. Bernoulli 25(4A), 2620–2648, 2019.
- Gwendoline de Bie, Gabriel Peyré, and Marco Cuturi. Stochastic Deep Networks. ICML, PMLR 97, 1556–1565, 2019.