Semi-discrete and Wasserstein-1

This chapter develops three computational consequences of duality. Eliminating one potential gives the semi-dual; for discrete measures, this viewpoint leads to auction algorithms, while a continuous source and discrete target lead to Laguerre-cell geometry. The final part specializes duality to $\Wass_1$, where Lipschitz functions and flow fields replace convex potentials. The material connects auction and network-flow methods Bertsekas, 1992Bertsekas & Eckstein, 1988, computational geometry Aurenhammer et al., 1998Mérigot, 2011Mérigot, 2013, and the Kantorovich--Rubinstein and Beckmann formulations Kantorovich & Rubinstein, 1958Beckmann, 1952.

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))

Semi-dual¶

The semi-dual eliminates one potential by an exact $c$-transform. It preserves concavity while removing the explicit pointwise inequality constraint.

General Measure Semi-dual¶

For arbitrary measures, partial maximization converts the constrained two-potential dual into an unconstrained optimization over one function.

Denote the extended full-dual objective by

Thus $\mathcal L_c(\alpha,\beta)=\max_{f,g}\mathcal E_0(f,g)$. For fixed $g$, feasibility is equivalent to $f\leq g^{\bar c}$. Since $\alpha$ is nonnegative, the largest admissible choice $f=g^{\bar c}$ maximizes the objective and gives

Partial maximization preserves concavity. Moreover, $\mathcal E(g+s)=\mathcal E(g)$ because $(g+s)^{\bar c}=g^{\bar c}-s$ and both measures have unit mass. Potentials are therefore defined only up to an additive constant, while the optimization is unconstrained.

Discrete Semi-dual¶

For two discrete measures

with common total mass $M$, use the same notation for vectors:

Eliminating the source vector gives

where

The function $\mathcal E$ is concave, piecewise affine, and invariant under $g\mapsto g+s\mathbf 1$. If ties are resolved by choosing $\sigma_g(i)\in\arg\min_j(\C_{ij}-g_j)$, then a supergradient is

It is therefore the mismatch between the desired target mass and the mass currently attracted by each target coordinate. At ties, splitting source mass among active targets describes the full superdifferential.

Auction Algorithm¶

The auction algorithm is derived from coordinate maximization of the semi-dual of the linear assignment problem. Its practical form uses bidder-specific dual-weight updates that cross the selected row’s next indifference threshold by $\varepsilon$; this controlled relaxation prevents jamming at nonsmooth ties. We follow the account of Mérigot and Thibert Mérigot & Thibert, 2020, itself based on Bertsekas’ auction algorithm and its $\varepsilon$-scaling refinement Bertsekas, 1981Bertsekas & Eckstein, 1988Bertsekas, 1992. In this section, $\C\in\mathbb R^{n\times n}$ and $n\ge2$. A permutation matrix $P$ represents the probability coupling $P/n$.

Coordinate Ascent and Discrete Laguerre Cells¶

Write $g_j$ for the target Kantorovich potential. Specializing the discrete semi-dual to uniform weights gives

Thus $(g^{\bar\C},g)$ is dual feasible. This is the discrete $\bar C$-transform of Remark Remark: Discrete $c$-transform, with the same sign convention as in the general and semi-discrete formulations.

The discrete Laguerre cells are

They can overlap at ties. They are the finite counterparts of the general semi-discrete Laguerre cells in Definition Definition: Laguerre Cells and Power Diagrams: the displayed cell is that definition restricted to row indices. If every row has a unique minimizer, then

An overfull cell therefore calls for a decrease of its dual weight. Proposition Proposition: Dual Certificate for an Assignment shows that a perfect matching in the contact graph $i\in\operatorname{Lag}^{\mathrm D}_j(g)$ certifies optimality of $g$; conversely, assignment duality and complementary slackness produce such a matching at every maximizer.

For $i\in\operatorname{Lag}^{\mathrm D}_j(g)$, define

When the cell is nonempty, the largest maximizing decrement along the negative $j$th coordinate is the largest such bid Mérigot & Thibert, 2020. At a tie, even this largest maximizing decrement can vanish, so naive coordinate ascent may jam before reaching a dual maximizer.

Bids and Relaxed Contacts¶

Auction avoids jamming by moving an unassigned row $\varepsilon$ beyond its next indifference point. Let $j_0$ and $j_1$ be its best and second-best targets:

The winning dual weight is updated by

Unlike exact coordinate maximization, this update uses the selected row’s bid rather than the maximum bid over the whole cell. Afterward, the reduced cost of $j_0$ is exactly $\varepsilon$ above that of the unchanged alternative $j_1$. The row nevertheless takes $j_0$; its former owner, if any, becomes unassigned. Because of this overshoot, a bid need not increase the nonsmooth semi-dual. Coordinate ascent motivates the update, but $\varepsilon$-complementary slackness is the invariant used in the convergence proof.

A bid gives the newly assigned row this property. Decreasing $g_{j_0}$ can only make $j_0$ less attractive to rows assigned elsewhere, and the previous owner of $j_0$ is removed. Every iteration therefore preserves the condition.

Figure Div shows actual iterates on the planar point clouds used in Div. For the current target potential, define the auction reduced costs

Exact zeros are the discrete Laguerre contacts of row $i$, whereas an owned edge only needs to satisfy $r_{ij}(g)\leq\varepsilon$ by $\varepsilon$-complementary slackness.

At an intermediate state, the current ownership matching is $M=\{(i,j):P_{ij}=1\}$; the labels in the figure report its cardinality.

show_book_figure("dual-auction-progression", width=850)

Geometric progression of the unit-mass transportation auction. Thick violet segments form the current partial ownership matching, while thin translucent segments show the $2n$ unmatched edges with smallest auction reduced costs. The labels report matching cardinality and cumulative bid count. With $\varepsilon=0.002$, the final assignment is reached after 505 bids and coincides with the exact squared-distance optimum.

Fixed Tolerance¶

For a prescribed tolerance, relaxed contact gives both an optimality certificate and a finite bound on the number of bids.

$\varepsilon$-Scaling¶

The cold-start estimate above exposes the limitation of a single tolerance. A large $\varepsilon$ is fast but gives only a coarse certificate, whereas the small $\varepsilon$ required for high accuracy can produce a bid count proportional to $R_{\C}/\varepsilon$. Continuation first learns a rough dual-potential landscape and then sharpens it instead of restarting from zero.

Algorithm Algorithm: Auction With $\varepsilon$-Scaling starts at the cost scale $\max\{R_{\C},\eta\}$ and halves the tolerance until it reaches the requested value $\eta$. Each phase rebuilds the ownership matrix but retains the target potential from the preceding phase. The previous complete assignment already certifies approximate contact for this potential, so the next auction refines an existing dual landscape.

A cold-started $\eta$-auction has bound $O(n^2(1+R_{\C}/\eta))$, while scaling uses a logarithmic number of $O(n^3)$ phases. Thus scaling is a high-accuracy guarantee, not an unconditional speedup: its bound improves the cold-start estimate when $R_{\C}/\eta$ is large compared with $n(1+\log_+(R_{\C}/\eta))$. For integer costs, choosing $\eta<1/n$ gives an exact assignment in $O(n^3(1+\log_+(nR_{\C})))$ operations.

Semi-discrete¶

The semi-discrete case is the setting where dual potentials become weights of Laguerre cells. This gives both geometry and algorithms for quantization and density fitting.

Discrete Targets and Laguerre Cells¶

Consider the case where

has distinct atoms and positive weights; zero-weight atoms can be removed. The same construction applies if $\alpha$ is discrete, after exchanging the roles of $\alpha$ and $\beta$. Restricting the minimization in Definition Definition: $c$-Transform to the support of $\beta$, equivalently applying that definition with the discrete target space $\Y=\{y_j\}_{j=1}^m$ and identifying a vector $g\in\RR^m$ with the function $g:\Y\to\RR$ defined by $g(y_j)=g_j$, gives the discrete $\bar c$-transform

This maps a vector $g$ to a continuous function because it is the minimum of finitely many continuous functions. Using this transform when $\beta$ is discrete yields the finite-dimensional semi-dual

The objective is invariant under $g\mapsto g+s\mathbf 1$, so one may impose the gauge $\sum_j g_j=0$.

The geometric object encoded by the dual weights is a weighted nearest-neighbor diagram: each source point is assigned to the target atom that realizes the discrete $\bar c$-transform.

For quadratic costs, varying the dual weights moves the walls between adjacent cells while keeping them parallel. This is the geometric mechanism by which the cell masses are adjusted.

Figure Div follows this adjustment from unweighted Voronoi cells to a power diagram whose cell masses match the prescribed discrete target weights.

show_book_figure("semidiscrete-laguerre-cells")

Laguerre cells for semi-discrete quadratic transport. The red contours show a continuous source density $\alpha$ given by a three-component Gaussian mixture on the right. The twenty-one colored circular sites are the atoms of the discrete target $\beta$, sampled from a compact cloud on the left; each site color matches its Laguerre cell. Starting from ordinary Voronoi cells, semi-dual weight updates deform the cells so that the $\alpha$-mass captured by each cell approaches the prescribed target mass.

The interactive demo exposes the dual-weight mechanism directly. Increase the number of weight updates to watch cells with too little mass expand and cells with too much mass shrink.

Mass Balance¶

The semi-dual energy can be rewritten as

The first-order optimality condition says that solving the semi-discrete dual amounts to choosing weights $g$ so that

The gradient components sum to zero, consistently with the gauge invariance. Conversely, balanced cells define the piecewise-constant map $T(x)=y_j$ on $\mathcal{L}_j(g)$. Its graph lies in the contact set $g^{\bar c}(x)+g_j=c(x,y_j)$, so continuous complementary slackness proves that both the map and the weights are optimal. For the quadratic cost, uniqueness follows from Brenier’s theorem when $\alpha$ has a density.

The sign of the gradient has a direct geometric interpretation. Increasing $g_j$ lowers the corresponding power distance and expands $\mathcal L_j(g)$; decreasing $g_j$ shrinks it. The dotted outline marks the balanced cell, so semi-dual ascent can be read as a mass-balancing procedure on a power diagram.

Figure Div makes the sign of this gradient geometric.

show_book_figure("semidiscrete-weight-gradient-cells", width=760)

Dual weights control Laguerre cell masses in the semi-discrete quadratic problem. The same blue target sites and red Gaussian source density are used in all panels; only the highlighted violet weight is changed. The dotted violet outline is the balanced cell. If the highlighted cell has too little source mass, then $b_j-\alpha(\mathcal L_j(g))>0$ and the ascent update increases the weight, expanding the cell outward. If it has too much mass, the update decreases the weight, shrinking it inward. At balance, the cell mass matches the prescribed target mass and the first-order update vanishes.

Quadratic power diagrams have polyhedral cells and can be computed efficiently using computational-geometry algorithms Aurenhammer, 1987Aurenhammer et al., 1998Mérigot, 2011. Expanding the cost shows that a cell minimizes $x\mapsto-2\langle x,y_j\rangle+\norm{y_j}^2-g_j$. The lower envelope of these affine functions gives the power diagram, while the lower hull of the lifted sites $(y_j,\norm{y_j}^2-g_j)$ gives its dual regular triangulation. For a planar source, this is a three-dimensional hull and Chan’s output-sensitive algorithm costs $O(m\log Q)$ for $Q$ hull vertices Chan, 1996. A three-dimensional source lifts to four dimensions and is not covered by that particular bound.

Stochastic Optimization¶

The semi-discrete formulation is useful because the objective is an expectation with respect to $\alpha$:

Away from cell boundaries, the stochastic gradient of the integrand is

an unbiased estimator of $\nabla\mathcal{E}(g)$ when cell boundaries have $\alpha$-measure zero. One can therefore maximize the semi-dual without first discretizing $\alpha$: the measure is used as a black box from which independent samples are drawn, a natural setup in high-dimensional statistics and machine learning.

Starting from $g^{(0)}=0$, stochastic gradient ascent draws $x_\ell\sim\alpha$ and performs

The stochastic supergradient has zero coordinate sum and preserves the gauge. For almost-sure stochastic-approximation convergence, one typically imposes

For example, one may use $\tau_\ell=\tau_0(1+\ell/\ell_0)^{-q}$ with $1/2<q\leq1$. The standard finite-horizon rate instead concerns averaged iterates.

This stochastic viewpoint is one of the main algorithmic advantages of the semi-discrete formulation Mérigot, 2011Genevay et al., 2016.

Optimal Quantization¶

Optimal quantization asks for the best discrete approximation of a measure by $m$ codepoints. It is the geometric core of vector quantization, compression and $k$-means clustering.

Free Masses and Prescribed Weights¶

The classical problem optimizes both codepoint positions and their probabilities. For a measure $\alpha$, define

This problem is classical in approximation theory and information theory Graf & Luschgy, 2000Lloyd, 1982.

The equal-weight case, $b_j=1/m$, prescribes the weights and is treated at the end of this section.

This deterministic rate mirrors the empirical optimal-transport sample-complexity rate: both are governed by the spacing $m^{-1/d}$ of points in dimension $d$. Quantization is best-case and deterministic, while empirical OT is random, but both display the same curse of dimensionality. Zador’s theorem further identifies the sharp asymptotic constant and limiting codepoint density Graf & Luschgy, 2000.

For fixed codepoints $Y$, the powered cost $b\mapsto\Wass_p^p(\alpha,\sum_jb_j\delta_{y_j})$ is convex. Its $p$-th root $\Wass_p$ need not be convex. The dependence on $Y$ is nonconvex and is generally computationally hard. The rest of this section distinguishes the free-mass Lloyd reduction from the fixed-weight geometry underlying finite-particle $\mathcal W_2$ gradient flows.

Lloyd Algorithm¶

The computational appeal of quantization comes from splitting the nonconvex search over sites into two elementary operations. For fixed sites, the optimal assignment is purely local: each point is sent to one of its nearest sites, and the resulting cells are Voronoi cells. This is the assignment step behind Lloyd’s algorithm and the $k$-means method.

Consequently, the quantization energy can be written in nearest-centroid form:

At a differentiability point of this energy, any local minimizer with nonempty cells satisfies the centroid condition

For the squared Euclidean cost, this becomes

Lloyd’s algorithm, also known as the $k$-means algorithm, iterates this fixed point: assign points to nearest sites, then replace each site by the centroid of its cell Lloyd, 1982. The assignment and centroid steps each minimize the appropriate block, so the objective cannot increase. Nonconvexity means this does not guarantee a global minimizer. For a finite data set with squared Euclidean loss, $k$-means++ gives an expected logarithmic approximation guarantee Arthur & Vassilvitskii, 2007.

Continuous Lloyd Flow¶

There is also an infinitesimal version of Lloyd’s fixed point, but it should first be understood on finite labelled configurations. Assume that $c(x,y)=\norm{x-y}^2$ and that $\alpha$ does not charge Voronoi boundaries. For a configuration $Y$, define, on nonempty cells,

as the cell mass and centroid. Empty cells are singular points of the vector field; one either freezes them, as in the algorithm below, or reseeds them. The relaxed step

is an explicit-Euler step for the cell-mass preconditioned gradient flow of the quantization energy $\mathcal F$,

Indeed, at differentiability points of $\mathcal F$, the envelope theorem gives

so that

Equivalently, this is the gradient flow of $\mathcal F$ for the site metric $g_Y(U,V)=2\sum_j a_j(Y)\langle u_j,v_j\rangle$; it is not the unweighted Euclidean gradient flow unless the masses are absorbed into the time step. Along smooth portions of the flow,

If $\eta_t=\sum_j w_j\delta_{y_j(t)}$ carries fixed positive weights, independent of the Voronoi masses $a_j(Y_t)$, this labelled particle ODE is equivalently a weak continuity equation,

in the sense of the measure evolutions introduced in Chapter Paragraph. The weights $w_j$ in this transport equation are auxiliary weights for the moving labelled particles; they are not the Voronoi masses used to define the quantization energy. If one records instead the free-weight projection

then, formally,

Thus the free-weight quantizer evolves by a balance equation, not by pure transport. This is why the construction is intrinsically finite-dimensional: Voronoi cells, centroids and labels define the velocity, and a canonical extension to arbitrary measures is not obtained by replacing $Y$ with the support of a measure. Indeed, any measure with dense support would have zero support-distance quantization error.

Mean-Field Limit and Ultrafast Diffusion¶

There is nevertheless a precise high-resolution continuum theory when the number $m$ of codepoints tends to infinity. If $\alpha=\rho\,\d x$, the limiting Eulerian variable is the density $\sigma$ of sites, meaning heuristically that $m\sigma(x)\,\d x$ codepoints lie in $\d x$. Thus the limit is $m\to+\infty$, not a limit in the exponent of the PDE. In one dimension, Caglioti, Golse and Iacobelli embed the ordered particle configuration in $L^2(0,1)$ and prove quantitative convergence of the discrete gradient flow toward a limiting flow Caglioti et al., 2015. A perturbative two-dimensional analysis around the optimal hexagonal lattice is developed in Caglioti et al., 2018. For $p$-quantization in dimension $d$, set $r=p/d$, so $r=2/d$ for the quadratic cost used in this section, and write $\mathcal F_p(Y)=\int\min_j\norm{x-y_j}^p\,\d\alpha(x)$. For well-prepared configurations whose empirical site distributions converge to $\sigma\,\d x$, the rescaled energy is described, up to a universal cell-shape constant, by

Formally,

Since the first variation is $-r\rho\sigma^{-r-1}$, the formal $\mathcal W_2$-gradient flow is the weighted ultrafast diffusion equation

with periodic or no-flux boundary conditions. Iacobelli studies the associated one-dimensional very-fast-diffusion equation and its convergence to equilibrium Iacobelli, 2019; Iacobelli, Patacchini and Santambrogio then use the JKO scheme and Wasserstein-gradient-flow tools to prove well-posedness, regularity estimates and convergence for a multidimensional weighted version Iacobelli et al., 2019. When $\rho>0$, set $\omega=\rho^{1/(r+1)}$ and $u=\sigma/\omega$. The same equation becomes

which makes the negative exponent, hence the ultrafast-diffusion character, explicit. Its stationary site density is proportional to $\rho^{d/(d+p)}$.

Figure Div shows the relaxed Lloyd flow in a two-dimensional toy problem.

show_book_figure("semidiscrete-lloyd-flow-mixtures")

Relaxed Lloyd flow from a source Gaussian-mixture initialization toward a different target Gaussian-mixture density. The blue contours and shading show the target density $\alpha$, while the colored disks are the moving codepoints initialized from the source mixture. The faint curves trace the labelled sites under the explicit-Euler Lloyd ODE. The right panel displays the relative quantization energy, illustrating the monotone decay of the objective along the relaxed iterations.

Figure Div follows the associated Voronoi cells and generators through Lloyd iterations, making the decrease of the quantization energy visible.

show_book_figure("semidiscrete-lloyd-quantization")

Lloyd quantization for the same continuous density and twenty-one initial sites as the Laguerre-cell figure. The red contours show the density $\alpha$, while the colored disks are the current codepoints and have the same colors as their Voronoi cells. The iterations move the initially left-located sites toward the high-density region and reshape the cells according to centroidal Voronoi geometry.

The interactive demo separates the nonconvex geometry from the fixed-point update: increase the iteration counter and watch sites migrate toward the density before settling into a local centroidal configuration.

Quantization with Fixed Equal Weights¶

The free-mass formulation above optimizes the positions and the weights of the atoms. A different problem is obtained by prescribing the weights. In the equal-weight case, set

and minimize only over the positions $Y=(y_j)_j$. Assume that the sites are distinct, that $\alpha$ has a density, and that cell boundaries have zero $\alpha$-mass. Let $C_j(Y)$ be the Laguerre cell transported to $y_j$, so that $\alpha(C_j(Y))=1/m$, and define its centroid

At differentiability points, the envelope theorem gives

Locally, while labels remain optimally matched, the $\mathcal W_2$ metric on equal-weight empirical measures induces the particle metric $g_Y(U,V)=m^{-1}\sum_j\langle u_j,v_j\rangle$. Hence the associated $\mathcal W_2$ gradient flow is the coupled system

Equivalently, $\nu_{Y_t}$ satisfies a continuity equation with velocity $v_t(y_j(t))=\bar x_j(Y_t)-y_j(t)$. This is the so-called finite-particle $\mathcal W_2$ gradient-flow viewpoint developed more systematically in Chapter Paragraph; the fixed-weight cells are Laguerre cells rather than the free-mass Voronoi cells used by Lloyd’s method.

Equal-Weight Quantization on the Line¶

The following classical scalar quantization result gives the precise form of the inverse-CDF rule for equal-weight quadratic quantization Graf & Luschgy, 2000. The atoms are not exactly the midpoint quantiles in general; they are the averages of the quantile function over equal mass bins. Midpoint inverse-CDF samples are nevertheless asymptotically equivalent and are often the most convenient rule in numerical examples.

Thus the common deterministic rule $m^{-1}\sum_i\delta_{Q((i-1/2)/m)}$ should be read as a midpoint approximation of the optimal bin-average formula. Orthogonal projection onto constants shows that it has the same leading squared error under the same smoothness assumptions; for the uniform law on $[0,1]$, both rules coincide and give the regular grid $y_i=(i-1/2)/m$. Random sampling has a different asymptotic regime.

Combining these propositions gives a sharp contrast between optimal placement and random placement on the line. Deterministic equal-weight quantization has squared error of order $m^{-2}$, hence $\mathcal W_2$ error of order $m^{-1}$, while i.i.d. empirical sampling has expected squared error of order $m^{-1}$, hence root-mean-square $\mathcal W_2$ error of order $m^{-1/2}$. This is consistent with broader empirical OT sample-complexity theory Dereich et al., 2013Fournier & Guillin, 2015Weed & Bach, 2019.

Figure Div illustrates both parts of this comparison: optimal atoms are uniform in quantile coordinates, and their error decays one power of $m$ faster than the root-mean-square empirical error.

show_book_figure("semidiscrete-quantile-quantization-rates")

One-dimensional equal-weight quantization in quantile coordinates. Left: for a smooth positive density on $[0,1]$, the colored atoms are bin averages of the inverse CDF over equal quantile intervals, while the gray atoms show one i.i.d. empirical draw with the same number of particles. Right: expected squared $\mathcal W_2$ errors. The deterministic bin averages and midpoint quantiles follow the $m^{-2}$ squared-error law, whereas i.i.d. empirical measures follow the slower $m^{-1}$ expected squared-error law.

Wasserstein-1 Norm¶

The $\Wass_1$ distance has an especially transparent dual: the admissible potentials are exactly 1-Lipschitz test functions. This makes $\Wass_1$ the meeting point between transport, PDE formulations and weak norms on signed measures.

c-Transform for Wasserstein-1¶

Assume that $d$ is a distance on $\X=\Y$ and take the ground cost $c(x,y)=d(x,y)$.

By the preceding proposition, a closed dual pair has the form $(f,-f)$ with $\Lip(f)\le1$. The Kantorovich dual therefore becomes the Kantorovich--Rubinstein formula

This expression depends only on the signed measure $\xi=\alpha-\beta$. On compact $\X$, the same supremum defines the Kantorovich--Rubinstein norm on finite signed measures with zero mass Kantorovich & Rubinstein, 1958. Homogeneity and the triangle inequality are immediate, while definiteness follows because Lipschitz functions separate finite Radon measures. On a noncompact pointed space, one uses normalized Lipschitz functions and measures with finite first moment.

For a discrete signed measure $\alpha-\beta=\sum_k r_k\delta_{z_k}$ with $\sum_k r_k=0$,

This finite-dimensional linear program can be solved by generic interior-point or first-order methods. If $N$ support points are involved, however, it still contains $O(N^2)$ Lipschitz constraints, mirroring the $O(nm)$ coupling variables of the original discrete Kantorovich LP; the dual formulation alone does not remove the all-pairs structure. The gain comes on structured metric spaces where the distance is generated locally: it is then enough to impose Lipschitz inequalities on neighboring pairs, because summing along paths recovers the constraints between arbitrary points. The one-dimensional ordered case is the first example; the graph-geodesic case is described later in Proposition Proposition: $\Wass_1$ and Beckmann Flow on a Graph.

When $d(x,y)=|x-y|$ on $\RR$, ordering the support points $z_1\le z_2\le\cdots$ reduces the constraints to neighboring pairs:

In one dimension this is equivalent to the cumulative formula given in the one-dimensional transport section.

Wasserstein-1 on Euclidean Spaces¶

In Euclidean space, the Lipschitz constraint has a local differential form and its dual variable is a flux. Let $\alpha,\beta$ have finite first moments and set $\xi=\alpha-\beta$. Rademacher’s theorem gives

The flux need not have a Lebesgue density: Dirac-to-Dirac transport already produces a measure concentrated on a segment. Let $m\in\mathcal M(\RR^d;\RR^d)$ be a vector-valued Radon measure, let $|m|$ denote its total variation, and define $\langle\operatorname{div}m,\varphi\rangle=-\int\langle\nabla\varphi,\d m\rangle$.