Wasserstein metric

Let ${\displaystyle \mu _{1}=\delta _{a_{1}}}$ and ${\displaystyle \mu _{2}=\delta _{a_{2}}}$ be two degenerate distributions (i.e. Dirac delta distributions) located at points ${\displaystyle a_{1}}$ and ${\displaystyle a_{2}}$ in ${\displaystyle \mathbb {R} }$. There is only one possible coupling of these two measures, namely the point mass ${\displaystyle \delta _{(a_{1},a_{2})}}$ located at ${\displaystyle (a_{1},a_{2})\in \mathbb {R} ^{2}}$. Thus, using the usual absolute value function as the distance function on ${\displaystyle \mathbb {R} }$, for any ${\displaystyle p\geq 1}$, the ${\displaystyle p}$-Wasserstein distance between ${\displaystyle \mu _{1}}$ and ${\displaystyle \mu _{2}}$ is

$${\displaystyle W_{p}(\mu _{1},\mu _{2})=|a_{1}-a_{2}|.}$$

By similar reasoning, if ${\displaystyle \mu _{1}=\delta _{a_{1}}}$ and ${\displaystyle \mu _{2}=\delta _{a_{2}}}$ are point masses located at points ${\displaystyle a_{1}}$ and ${\displaystyle a_{2}}$ in ${\displaystyle \mathbb {R} ^{n}}$, and we use the usual Euclidean norm on ${\displaystyle \mathbb {R} ^{n}}$ as the distance function, then

$${\displaystyle W_{p}(\mu _{1},\mu _{2})=\|a_{1}-a_{2}\|_{2}.}$$

If ${\displaystyle P}$ is an empirical measure with samples ${\displaystyle X_{1},\ldots ,X_{n}}$ and ${\displaystyle Q}$ is an empirical measure with samples ${\displaystyle Y_{1},\ldots ,Y_{n}}$, the distance is a simple function of the order statistics:

$${\displaystyle W_{p}(P,Q)=\left({\frac {1}{n}}\sum _{i=1}^{n}\|X_{(i)}-Y_{(i)}\|^{p}\right)^{1/p}.}$$

If ${\displaystyle P}$ and ${\displaystyle Q}$ are empirical distributions, each based on ${\displaystyle n}$ observations, then

$${\displaystyle W_{p}(P,Q)=\inf _{\pi }\left({\frac {1}{n}}\sum _{i=1}^{n}\|X_{i}-Y_{\pi (i)}\|^{p}\right)^{1/p},}$$

where the infimum is over all permutations ${\displaystyle \pi }$ of ${\displaystyle n}$ elements. This is a linear assignment problem, and can be solved by the Hungarian algorithm in cubic time.

Let ${\displaystyle \mu _{1}={\mathcal {N}}(m_{1},C_{1})}$ and ${\displaystyle \mu _{2}={\mathcal {N}}(m_{2},C_{2})}$ be two non-degenerate Gaussian measures (i.e. normal distributions) on ${\displaystyle \mathbb {R} ^{n}}$, with respective expected values ${\displaystyle m_{1}}$ and ${\displaystyle m_{2}\in \mathbb {R} ^{n}}$ and symmetric positive semi-definite covariance matrices ${\displaystyle C_{1}}$ and ${\displaystyle C_{2}\in \mathbb {R} ^{n\times n}}$. Then,[3] with respect to the usual Euclidean norm on ${\displaystyle \mathbb {R} ^{n}}$, the 2-Wasserstein distance between ${\displaystyle \mu _{1}}$ and ${\displaystyle \mu _{2}}$ is

$${\displaystyle W_{2}(\mu _{1},\mu _{2})^{2}=\|m_{1}-m_{2}\|_{2}^{2}+\mathop {\mathrm {trace} } {\bigl (}C_{1}+C_{2}-2{\bigl (}C_{2}^{1/2}C_{1}C_{2}^{1/2}{\bigr )}^{1/2}{\bigr )}.}$$

where ${\displaystyle C^{1/2}}$ denotes the principal square root of ${\displaystyle C}$. Note that the second term (involving the trace) is precisely the (unnormalised) Bures metric between ${\displaystyle C_{1}}$ and ${\displaystyle C_{2}}$. This result generalises the earlier example of the Wasserstein distance between two point masses (at least in the case ${\displaystyle p=2}$), since a point mass can be regarded as a normal distribution with covariance matrix equal to zero, in which case the trace term disappears and only the term involving the Euclidean distance between the means remains.

Let ${\displaystyle \mu _{1},\mu _{2}\in P_{p}(\mathbb {R} )}$ be probability measures on ${\displaystyle \mathbb {R} }$, and denote their cumulative distribution functions by ${\displaystyle F_{1}(x)}$ and ${\displaystyle F_{2}(x)}$. Then the transport problem has an analytic solution: Optimal transport preserves the order of probability mass elements, so the mass at quantile ${\displaystyle q}$ of ${\displaystyle \mu _{1}}$ moves to quantile ${\displaystyle q}$ of ${\displaystyle \mu _{2}}$. Thus, the ${\displaystyle p}$-Wasserstein distance between ${\displaystyle \mu _{1}}$ and ${\displaystyle \mu _{2}}$ is

$${\displaystyle W_{p}(\mu _{1},\mu _{2})=\left(\int _{0}^{1}\left|F_{1}^{-1}(q)-F_{2}^{-1}(q)\right|^{p}\,\mathrm {d} q\right)^{1/p},}$$

where ${\displaystyle F_{1}^{-1}}$ and ${\displaystyle F_{2}^{-1}}$ are the quantile functions (inverse CDFs). In the case of ${\displaystyle p=1}$, a change of variables leads to the formula

$${\displaystyle W_{1}(\mu _{1},\mu _{2})=\int _{\mathbb {R} }\left|F_{1}(x)-F_{2}(x)\right|\,\mathrm {d} x.}$$

It can be shown that W_p satisfies all the axioms of a metric on P_p(M). Furthermore, convergence with respect to W_p is equivalent to the usual weak convergence of measures plus convergence of the first pth moments.[10]

The following dual representation of W₁ is a special case of the duality theorem of Kantorovich and Rubinstein (1958): when μ and ν have bounded support,

$${\displaystyle W_{1}(\mu ,\nu )=\sup \left\{\left.\int _{M}f(x)\,\mathrm {d} (\mu -\nu )(x)\right|{\text{continuous }}f:M\to \mathbb {R} ,\operatorname {Lip} (f)\leq 1\right\},}$$

where Lip(f) denotes the minimal Lipschitz constant for f. This form shows that W₁ is an integral probability metric.

Compare this with the definition of the Radon metric:

$${\displaystyle \rho (\mu ,\nu ):=\sup \left\{\left.\int _{M}f(x)\,\mathrm {d} (\mu -\nu )(x)\right|{\text{continuous }}f:M\to [-1,1]\right\}.}$$

If the metric d of the metric space (M,d) is bounded by some constant C, then

$${\displaystyle 2W_{1}(\mu ,\nu )\leq C\rho (\mu ,\nu ),}$$

and so convergence in the Radon metric (identical to total variation convergence when M is a Polish space) implies convergence in the Wasserstein metric, but not vice versa.

The following is an intuitive proof which skips over technical points. A fully rigorous proof is found in.[11]

Discrete case: When ${\displaystyle M}$ is discrete, solving for the 1-Wasserstein distance is a problem in linear programming:

$${\displaystyle {\begin{cases}\min _{\gamma }\sum _{x,y}c(x,y)\gamma (x,y)\\\sum _{y}\gamma (x,y)=\mu (x)\\\sum _{x}\gamma (x,y)=\nu (y)\\\gamma \geq 0\end{cases}}}$$

where ${\displaystyle c:M\times M\to [0,\infty )}$ is a general "cost function".

By carefully writing the above equations as matrix equations, we obtain its dual problem:[12]

$${\displaystyle {\begin{cases}\max _{f,g}\sum _{x}\mu (x)f(x)+\sum _{y}\nu (y)g(y)\\f(x)+g(y)\leq c(x,y)\end{cases}}}$$

and by the duality theorem of linear programming, since the primal problem is feasible and bounded, so is the dual problem, and the minimum in the first problem equals the maximum in the second problem. That is, the problem pair exhibits strong duality.

For the general case, the dual problem is found by converting sums to integrals:

$${\displaystyle {\begin{cases}\sup _{f,g}\mathbb {E} _{x\sim \mu }[f(x)]+\mathbb {E} _{y\sim \nu }[g(y)]\\f(x)+g(y)\leq c(x,y)\end{cases}}}$$

and the strong duality still holds. This is the Kantorovich duality theorem. Cédric Villani recounts the following interpretation from Luis Caffarelli:[13]

Suppose you want to ship some coal from mines, distributed as ${\displaystyle \mu }$, to factories, distributed as ${\displaystyle \nu }$. The cost function of transport is ${\displaystyle c}$. Now a shipper comes and offers to do the transport for you. You would pay him ${\displaystyle f(x)}$ per coal for loading the coal at ${\displaystyle x}$, and pay him ${\displaystyle g(y)}$ per coal for unloading the coal at ${\displaystyle y}$.

For you to accept the deal, the price schedule must satisfy ${\displaystyle f(x)+g(y)\leq c(x,y)}$. The Kantorovich duality states that the shipper can make a price schedule that makes you pay almost as much as you would ship yourself.

This result can be pressed further to yield:

Theorem (Kantorovich-Rubenstein duality) — When the probability space ${\displaystyle \Omega }$ is a metric space, then for any fixed ${\displaystyle K>0}$,

$${\displaystyle W_{1}(\mu ,\nu )={\frac {1}{K}}\sup _{\|f\|_{L}\leq K}\mathbb {E} _{x\sim \mu }[f(x)]-\mathbb {E} _{y\sim \nu }[f(y)]}$$

where ${\displaystyle \|\cdot \|_{L}}$ is the Lipschitz norm.

Proof

It suffices to prove the case of ${\displaystyle K=1}$. Start with

$${\displaystyle W_{1}(\mu ,\nu )=\sup _{f(x)+g(y)\leq d(x,y)}\mathbb {E} _{x\sim \mu }[f(x)]+\mathbb {E} _{y\sim \nu }[g(y)].}$$

Then, for any choice of ${\displaystyle g}$, one can push the term higher by setting ${\displaystyle f(x)=\inf _{y}d(x,y)-g(y)}$, making it an infimal convolution of ${\displaystyle -g}$ with a cone. This implies ${\displaystyle f(x)-f(y)\leq d(x,y)}$ for any ${\displaystyle x,y}$, that is, ${\displaystyle \|f\|_{L}\leq 1}$.

Thus,

$${\displaystyle {\begin{aligned}W_{1}(\mu ,\nu )&=\sup _{g}\sup _{f(x)+g(y)\leq d(x,y)}\mathbb {E} _{x\sim \mu }[f(x)]+\mathbb {E} _{y\sim \nu }[g(y)]\\&=\sup _{g}\sup _{\|f\|_{L}\leq 1,f(x)+g(y)\leq d(x,y)}\mathbb {E} _{x\sim \mu }[f(x)]+\mathbb {E} _{y\sim \nu }[g(y)]\\&=\sup _{\|f\|_{L}\leq 1}\sup _{g,f(x)+g(y)\leq d(x,y)}\mathbb {E} _{x\sim \mu }[f(x)]+\mathbb {E} _{y\sim \nu }[g(y)].\end{aligned}}}$$

Next, for any choice of ${\displaystyle \|f\|_{L}\leq 1}$, ${\displaystyle g}$ can be optimized by setting ${\displaystyle g(y)=\inf _{x}d(x,y)-f(x)}$. Since ${\displaystyle \|f\|_{L}\leq 1}$, this implies ${\displaystyle g(y)=-f(y)}$.

Infimal convolution of a cone with a curve. Note how the lower envelope has slope ${\displaystyle \leq 1}$, and how the lower envelope is equal to the curve on the parts where the curve itself has slope ${\displaystyle \leq 1}$.

The two infimal convolution steps are visually clear when the probability space is ${\displaystyle \mathbb {R} }$.

For notational convenience, let ${\displaystyle \square }$ denote the infimal convolution operation.

For the first step, where we used ${\displaystyle f=cone\square (-g)}$, plot out the curve of ${\displaystyle -g}$, then at each point, draw a cone of slope 1, and take the lower envelope of the cones as ${\displaystyle f}$, as shown in the diagram, then ${\displaystyle f}$ cannot increase with slope larger than 1. Thus all its secants have slope ${\displaystyle {\Big |}{\frac {f(x)-f(y)}{x-y}}{\Big |}\leq 1}$.

For the second step, picture the infimal convolution ${\displaystyle cone\square (-f)}$, then if all secants of ${\displaystyle f}$ have slope at most 1, then the lower envelope of ${\displaystyle cone\square (-f)}$ are just the cone-apices themselves, thus ${\displaystyle cone\square (-f)=-f}$.

1D Example. When both ${\displaystyle \mu ,\nu }$ are distributions on ${\displaystyle \mathbb {R} }$, then integration by parts give

$${\displaystyle \mathbb {E} _{x\sim \mu }[f(x)]-\mathbb {E} _{y\sim \nu }[f(y)]=\int f'(x)(F_{\nu }(x)-F_{\mu }(x))dx,}$$

thus

$${\displaystyle f(x)=K\cdot \mathop {sign} (F_{\nu }(x)-F_{\mu }(x)).}$$

Benamou & Brenier found a dual representation of ${\displaystyle W_{2}}$ by fluid mechanics, which allows efficient solution by convex optimization.[14][15]

Given two probability distributions on ${\displaystyle \mathbb {R} ^{n}}$ with density ${\displaystyle p,q}$, then

$${\displaystyle W_{2}(p,q)=\min _{v}\int _{0}^{T}\int _{\mathbb {R} ^{n}}\|v(x,t)\|^{2}\rho (x,t)dxdt}$$

where ${\displaystyle v(x,t)}$

is a velocity field, and ${\displaystyle \rho (x,t)}$ is the fluid density field, such that

$${\displaystyle {\begin{aligned}&{\frac {\partial \rho (\cdot ,t)}{\partial t}}=-\operatorname {div} \left(\rho (\cdot ,t)v\right)\\&\rho (x,0)=p\\&\rho (x,T)=q.\end{aligned}}}$$

That is, the mass should be conserved, and the velocity field should transport the probability distribution ${\displaystyle p}$ to ${\displaystyle q}$ during the time interval ${\displaystyle [0,T]}$.

Under suitable assumptions, the Wasserstein distance ${\displaystyle W_{2}}$ of order two is Lipschitz equivalent to a negative-order homogeneous Sobolev norm.[16] More precisely, if we take ${\displaystyle M}$ to be a connected Riemannian manifold equipped with a positive measure ${\displaystyle \pi }$, then we may define for ${\displaystyle f\colon M\to \mathbb {R} }$ the seminorm

$${\displaystyle \|f\|_{{\dot {H}}^{1}(\pi )}^{2}=\int _{M}|\nabla f(x)|^{2}\,\pi (\mathrm {d} x)}$$

and for a signed measure ${\displaystyle \mu }$ on ${\displaystyle M}$ the dual norm

$${\displaystyle \|\mu \|_{{\dot {H}}^{-1}(\pi )}=\sup {\bigg \{}|\langle f,\mu \rangle |\,{\bigg |}\,\|f\|_{{\dot {H}}^{1}(\pi )}\leq 1{\bigg \}}.}$$

Then any two probability measures ${\displaystyle \mu }$ and ${\displaystyle \nu }$ on ${\displaystyle M}$ satisfy the upper bound

$${\displaystyle W_{2}(\mu ,\nu )\leq 2\|\mu -\nu \|_{{\dot {H}}^{-1}(\mu )}.}$$

In the other direction, if ${\displaystyle \mu }$ and ${\displaystyle \nu }$ each have densities with respect to the standard volume measure on ${\displaystyle M}$ that are both bounded above by some ${\displaystyle 0<C<\infty }$, and ${\displaystyle M}$ has non-negative Ricci curvature, then

$${\displaystyle \|\mu -\nu \|_{{\dot {H}}^{-1}(\pi )}\leq {\sqrt {C}}W_{2}(\mu ,\nu ).}$$

For any p ≥ 1, the metric space (P_p(M), W_p) is separable, and is complete if (M, d) is separable and complete.[17]