One-dimensional optimal transport is often introduced as a sorting problem: line up the source and target quantiles, and match equal ranks. This gives an exact formula for Wasserstein distance. But sometimes a cumulative distribution function is easy to evaluate while its inverse is inconvenient. Gaussian mixtures are one example; data distributed across many machines are another.
There is an exact alternative for the quadratic distance which uses only CDF values. The formula has a short geometric proof, works with atoms, and naturally suggests integration and sampling algorithms. The identity is recorded by Bobkov and Ledoux, who credit Eustasio del Barrio [1]. Its recent use for data-parallel sliced-Wasserstein estimation is due to Vauthier, Mérigot, and Korba [2]. The formula also appears as Proposition 2.45 in my book [3].
From ranks to cumulative masses
For a probability law \(\alpha\) on \(\mathbb R\), write
\[ F_\alpha(x)=\alpha((-\infty,x]), \qquad Q_\alpha(r)=\inf\{x:F_\alpha(x)\ge r\}. \]The square of the quadratic Wasserstein distance is the least mean squared displacement over all couplings. If both laws have finite second moment, monotone matching gives
\[ W_2^2(\alpha,\beta) =\int_0^1|Q_\alpha(r)-Q_\beta(r)|^2\,dr. \]The aim is not to change the optimal coupling. It is to integrate the cost of that coupling in a different order. Define the nonnegative function, for \(x<y\),
\[ H(x,y) =(F_\alpha(x)-F_\beta(y))_+ +(F_\beta(x)-F_\alpha(y))_+, \]where \(a_+=\max(a,0)\). The first term measures how many quantile levels have their source position to the left of \(x\) and target position to the right of \(y\). The second records motion in the opposite direction.
Theorem 1 (the CDF-only formula). For \(\alpha,\beta\in\mathcal P_2(\mathbb R)\),
\[ W_2^2(\alpha,\beta) =2\int_{\mathbb R}\int_x^\infty H(x,y)\,dy\,dx. \]If both laws are supported in \([a,b]\), the integration region reduces to the triangle \(a<x<y<b\). No continuity or density assumption is needed.
Proof: integrate the triangle swept out by each matched pair
For \(u<v\), the region \(u\le x<y<v\) is a triangle of area \((v-u)^2/2\). Thus, for arbitrary \(u,v\),
\[ |u-v|^2 =2\iint_{x\lt y} \left[ \mathbf1_{\{u\le x,\ y\lt v\}} +\mathbf1_{\{v\le x,\ y\lt u\}} \right]dx\,dy. \]Substitute \(u=Q_\alpha(r)\), \(v=Q_\beta(r)\), and integrate over \(r\in(0,1)\). The generalized-inverse identities give, up to irrelevant endpoints,
\[ \begin{aligned} Q_\alpha(r)\le x&\iff r\le F_\alpha(x),\\ Q_\beta(r)>y&\iff r>F_\beta(y). \end{aligned} \]The set of \(r\)'s satisfying both has length \((F_\alpha(x)-F_\beta(y))_+\). The reverse direction gives the other positive part. All integrands are nonnegative, so Tonelli's theorem permits the change in integration order, even before finiteness is established. The second moments make the resulting Wasserstein value finite. \(\square\)
The triangular region matters: simply replacing it by the full square would change the answer. So does the use of two different thresholds \(x,y\). The expression \(\int|F_\alpha(x)-F_\beta(x)|^2\,dx\) is a Cramér-type discrepancy, not \(W_2^2\). In contrast, \(W_1=\int|F_\alpha-F_\beta|\) does use a single threshold.
An exactly solvable example
Take the uniform law \(\beta\) on \([0,1]\), and let \(\alpha\) have density
\[ \rho_\alpha(x)=1+\lambda\sin(2\pi x), \qquad |\lambda|<1. \]Its CDF is explicit:
\[ F_\alpha(x)=x+\frac{\lambda}{2\pi}(1-\cos(2\pi x)), \qquad 0\le x\le1. \]Although its inverse has no comparably simple expression, the exact transport cost is
\[ W_2^2(\alpha,\beta)=\frac{3\lambda^2}{8\pi^2}. \]Proof: for a uniform target, the source CDF is the transport map
The density is strictly positive, so \(F_\alpha\) is increasing and sends \(\alpha\) to the uniform law. It is therefore the optimal map. Consequently,
\[ W_2^2(\alpha,\beta) =\int_0^1|F_\alpha(x)-x|^2 (1+\lambda\sin(2\pi x))\,dx. \]The integral of \((1-\cos(2\pi x))^2\) is \(3/2\). Its product with \(\sin(2\pi x)\) integrates to zero, either by symmetry under \(x\mapsto1-x\) or by direct differentiation. This gives the claimed value. \(\square\)
This provides an independent reference for the numerical CDF integral. With \(\lambda=0.6\), the value is \(0.135/\pi^2\). Refining the integration grid converges to this value, although individual errors need not decrease monotonically with every choice of grid size. The identity itself has no numerical error.
Sampling the triangle
For compactly supported laws, the formula also gives a particularly simple unbiased Monte Carlo estimator. Draw two independent uniforms on \([a,b]\), sort them into \(X<Y\), and evaluate \(H(X,Y)\). The sorted pair is uniform on the triangle, whose area is \((b-a)^2/2\).
Proposition 1 (an unbiased CDF estimator). For independent triangle samples \((X_\ell,Y_\ell)\), set
\[ \widehat D_L =\frac{(b-a)^2}{L}\sum_{\ell=1}^L H(X_\ell,Y_\ell). \]Then
\[ \mathbb E\widehat D_L=W_2^2(\alpha,\beta), \qquad \operatorname{Var}(\widehat D_L) \le\frac{(b-a)^4}{4L}. \]Proof: use the area factor and a bounded random variable
Theorem 1 and the triangle area give \(W_2^2=(b-a)^2\mathbb EH(X,Y)\). For \(x<y\), the two positive parts cannot both be positive. Otherwise monotonicity of the CDFs would imply
\[ F_\alpha(x)>F_\beta(y)\ge F_\beta(x) >F_\alpha(y)\ge F_\alpha(x), \]a contradiction. Hence \(0\le H\le1\). A random variable in \([0,1]\) has variance at most \(1/4\), and averaging \(L\) independent copies divides the variance by \(L\). \(\square\)
This unbiasedness assumes exact CDF evaluations. If the CDFs are estimated from data, the estimator is conditionally unbiased for the distance between those empirical laws, not automatically for the distance between the population laws: the positive part is nonlinear.
Errors in evaluating genuine CDFs can nevertheless be controlled. If two approximating laws on \([a,b]\) satisfy \(\|\widehat F_\alpha-F_\alpha\|_\infty\le\delta_\alpha\) and \(\|\widehat F_\beta-F_\beta\|_\infty\le\delta_\beta\), then
\[ \left|W_2^2(\widehat\alpha,\widehat\beta) -W_2^2(\alpha,\beta)\right| \le (b-a)^2(\delta_\alpha+\delta_\beta). \]Indeed, the two positive parts cannot be simultaneously positive, so \(H\) is the maximum of zero and the two CDF differences. Changing each difference by at most \(\delta_\alpha+\delta_\beta\) changes their maximum by at most that amount. Integrate over a triangle of area \((b-a)^2/2\). This bound concerns CDF approximation; quadrature or Monte Carlo error must still be added separately.
When avoiding inversion helps
A Gaussian-mixture CDF is a weighted sum of Gaussian CDFs. It can be evaluated directly, whereas evaluating its inverse requires root finding. For unbounded laws, restricting the formula to a finite triangle has a precise interpretation: it computes the distance between the laws clipped to that interval.
Controlled truncation: what a finite integration interval actually computes
For \(I=[a,b]\), let \(P_I(x)=\max(a,\min(b,x))\), \(\alpha_I=(P_I)_\#\alpha\), and \(\beta_I=(P_I)_\#\beta\). Their CDFs agree with the original ones in the interior of \(I\). Therefore
\[ D_I:=2\int_a^b\int_x^b H(x,y)\,dy\,dx =W_2^2(\alpha_I,\beta_I). \]Define the root-mean-square clipping displacement by
\[ \delta_I(\alpha) =\left(\int\operatorname{dist}(x,I)^2\,d\alpha(x)\right)^{1/2}. \]Clipping both coordinates of a coupling cannot increase its displacement, because \(P_I\) is 1-Lipschitz. Hence \(W_2(\alpha_I,\beta_I)\le W_2(\alpha,\beta)\). The triangle inequality and the explicit couplings \(x\mapsto P_I(x)\) give the reverse estimate with a tail error:
\[ 0\le W_2(\alpha,\beta)-\sqrt{D_I} \le\delta_I(\alpha)+\delta_I(\beta). \]Thus finite-window integration underestimates the true distance by a controlled amount. The relevant tail quantity is a squared displacement, not just the probability of being outside the interval.
For distributed empirical data, each worker can contribute its local count below a threshold. Adding these weighted counts gives the global CDF. Crucially, one must aggregate the CDFs before taking the positive parts; averaging local Wasserstein distances solves a different problem.
The same idea applies after projecting high-dimensional measures onto a line. For a direction \(\theta\), compute the CDFs of \(\theta^\top X\) and \(\theta^\top Y\), apply the formula, and average the resulting squared distances over uniformly sampled directions. This estimates the square of the sliced Wasserstein distance; take a square root to obtain the distance itself. It is a different metric from the full high-dimensional Wasserstein distance. The recent paper [2] develops estimators tailored to this setting.
There is a computational tradeoff: a two-dimensional quadrature is not automatically cheaper than one-dimensional quantile integration or sorting a sample. The benefit is access to a different computational primitive—CDF evaluation—which can be particularly convenient for structured distributions and distributed data. For the broader algorithmic picture, see [4].
References
- Sergey Bobkov and Michel Ledoux. One-Dimensional Empirical Measures, Order Statistics, and Kantorovich Transport Distances. Memoirs of the American Mathematical Society 261(1259), 2019. Author-hosted PDF.
- Christophe Vauthier, Quentin Mérigot, and Anna Korba. Highly Data Parallelizable Estimation of the Sliced-Wasserstein Distance Using Cumulative Distribution Functions. arXiv:2606.30310, 2026.
- Gabriel Peyré. Optimal Transport for Machine Learners, 2025–2026. Proposition 2.45, p. 31.
- Gabriel Peyré and Marco Cuturi. Computational Optimal Transport. Foundations and Trends in Machine Learning 11(5–6), 355–607, 2019.