Transportation theory (mathematics)

In mathematics and economics, transportation theory or transport theory is a name given to the study of optimal transportation and allocation of resources. The problem was formalized by the French mathematician Gaspard Monge in 1781.^[1]

In the 1920s A. N. Tolstoi was one of the first to study the transportation problem mathematically. In 1930, in the collection Transportation Planning Volume I for the National Commissariat of Transportation of the Soviet Union, he published a paper "Methods of Finding the Minimal Kilometrage in Cargo-transportation in space".^[2][3]

Major advances were made in the field during World War II by the Soviet mathematician and economist Leonid Kantorovich.^[4] Consequently, the problem as it is stated is sometimes known as the Monge–Kantorovich transportation problem.^[5] The linear programming formulation of the transportation problem is also known as the Hitchcock–Koopmans transportation problem.^[6]

Motivation

Mines and factories

Two one-dimensional distributions ${\displaystyle \mu }$ and ${\displaystyle \nu }$, plotted on the ${\displaystyle x}$ and ${\displaystyle y}$. The two distributions can be pictured as two piles of dirt, one before moving, and one after moving. The heatmap in the center is a transport plan and denotes where each atom of dirt would be moved to.

Suppose that we have a collection of ${\displaystyle m}$ mines mining iron ore, and a collection of ${\displaystyle n}$ factories which use the iron ore that the mines produce. Suppose for the sake of argument that these mines and factories form two disjoint subsets ${\displaystyle M}$ and ${\displaystyle F}$ of the Euclidean plane ${\displaystyle \mathbb {R} ^{2}}$. Suppose also that we have a cost function ${\displaystyle c:\mathbb {R} ^{2}\times \mathbb {R} ^{2}\to [0,\infty )}$, so that ${\displaystyle c(x,y)}$ is the cost of transporting one shipment of iron from ${\displaystyle x}$ to ${\displaystyle y}$. For simplicity, we ignore the time taken to do the transporting. We also assume that each mine can supply only one factory (no splitting of shipments) and that each factory requires precisely one shipment to be in operation (factories cannot work at half- or double-capacity). Having made the above assumptions, a transport plan is a bijection ${\displaystyle T:M\to F}$. In other words, each mine ${\displaystyle m\in M}$ supplies precisely one target factory ${\displaystyle T(m)\in F}$ and each factory is supplied by precisely one mine. We wish to find the optimal transport plan, the plan ${\displaystyle T}$ whose total cost

${\displaystyle c(T):=\sum _{m\in M}c(m,T(m))}$

is the least of all possible transport plans from ${\displaystyle M}$ to ${\displaystyle F}$. This motivating special case of the transportation problem is an instance of the assignment problem. More specifically, it is equivalent to finding a minimum weight matching in a bipartite graph.

This can be generalized to the continuous case, where there are infinitely many mines and factories distributed on the real line, or generally in any metric space. This case is usually pictured as "changing the shape of a pile of dirt", and thus called the earth mover's problem.

Moving books: the importance of the cost function

The following simple example illustrates the importance of the cost function in determining the optimal transport plan. Suppose that we have ${\displaystyle n}$ books of equal width on a shelf (the real line), arranged in a single contiguous block. We wish to rearrange them into another contiguous block, but shifted one book-width to the right. Two obvious candidates for the optimal transport plan present themselves:

1. move all ${\displaystyle n}$ books one book-width to the right ("many small moves");
2. move the left-most book ${\displaystyle n}$ book-widths to the right and leave all other books fixed ("one big move").

If the cost function is proportional to Euclidean distance (${\displaystyle c(x,y)=\alpha \|x-y\|}$ for some ${\displaystyle \alpha >0}$) then these two candidates are both optimal. If, on the other hand, we choose the strictly convex cost function proportional to the square of Euclidean distance (${\displaystyle c(x,y)=\alpha \|x-y\|^{2}}$ for some ${\displaystyle \alpha >0}$), then the "many small moves" option becomes the unique minimizer.

Note that the above cost functions consider only the horizontal distance traveled by the books, not the horizontal distance traveled by a device used to pick each book up and move the book into position. If the latter is considered instead, then, of the two transport plans, the second is always optimal for the Euclidean distance, while, provided there are at least 3 books, the first transport plan is optimal for the squared Euclidean distance.

Hitchcock problem

The following transportation problem formulation is credited to F. L. Hitchcock:^[7]

Suppose there are ${\displaystyle m}$ sources ${\displaystyle x_{1},\ldots ,x_{m}}$ for a commodity, with ${\displaystyle a(x_{i})}$ units of supply at ${\displaystyle x_{i}}$ and ${\displaystyle n}$ sinks ${\displaystyle y_{1},\ldots ,y_{n}}$ for the commodity, with the demand ${\displaystyle b(y_{j})}$ at ${\displaystyle y_{j}}$. If ${\displaystyle c(x_{i},\ y_{j})}$ is the unit cost of shipment from ${\displaystyle x_{i}}$ to ${\displaystyle y_{j}}$, find a flow that satisfies demand from supplies and minimizes the flow cost. This challenge in logistics was taken up by D. R. Fulkerson^[8] and in the book Flows in Networks (1962) written with L. R. Ford Jr.^[9]

Tjalling Koopmans is also credited with formulations of transport economics and allocation of resources.

Abstract formulation of the problem

Monge and Kantorovich formulations

The transportation problem as it is stated in modern or more technical literature looks somewhat different because of the development of Riemannian geometry and measure theory. The mines-factories example, simple as it is, is a useful reference point when thinking of the abstract case. In this setting, we allow the possibility that we may not wish to keep all mines and factories open for business, and allow mines to supply more than one factory, and factories to accept iron from more than one mine.

Let ${\displaystyle X}$ and ${\displaystyle Y}$ be two separable metric spaces such that any probability measure on ${\displaystyle X}$ (or ${\displaystyle Y}$) is a Radon measure (i.e. they are Radon spaces). Let ${\displaystyle c:X\times Y\to [0,\infty )}$ be a Borel-measurable function. Given probability measures ${\displaystyle \mu }$ on ${\displaystyle X}$ and ${\displaystyle \nu }$ on ${\displaystyle Y}$, Monge's formulation of the optimal transportation problem is to find a transport map ${\displaystyle T:X\to Y}$ that realizes the infimum

${\displaystyle \inf \left\{\left.\int _{X}c(x,T(x))\,\mathrm {d} \mu (x)\right|T_{*}(\mu )=\nu \right\},}$

where ${\displaystyle T_{*}(\mu )}$ denotes the push forward of ${\displaystyle \mu }$ by ${\displaystyle T}$. A map ${\displaystyle T}$ that attains this infimum (i.e. makes it a minimum instead of an infimum) is called an "optimal transport map".

Monge's formulation of the optimal transportation problem can be ill-posed, because sometimes there is no ${\displaystyle T}$ satisfying ${\displaystyle T_{*}(\mu )=\nu }$: this happens, for example, when ${\displaystyle \mu }$ is a Dirac measure but ${\displaystyle \nu }$ is not.

We can improve on this by adopting Kantorovich's formulation of the optimal transportation problem, which is to find a probability measure ${\displaystyle \gamma }$ on ${\displaystyle X\times Y}$ that attains the infimum

${\displaystyle \inf \left\{\left.\int _{X\times Y}c(x,y)\,\mathrm {d} \gamma (x,y)\right|\gamma \in \Gamma (\mu ,\nu )\right\},}$

where ${\displaystyle \Gamma (\mu ,\nu )}$ denotes the collection of all probability measures on ${\displaystyle X\times Y}$ with marginals ${\displaystyle \mu }$ on ${\displaystyle X}$ and ${\displaystyle \nu }$ on ${\displaystyle Y}$.

Cost duality

Example c-duality transform, where c(x, y) = 2(cos(3x) + 1)|y − x|² + (4 − 2(cos(3x) + 1))|y − x|⁴ and ${\displaystyle \psi (x)=-e^{-x^{2}}}$.

Given a cost function ${\displaystyle c(x,y)}$, it produces a duality transformation ${\displaystyle \psi \mapsto \psi ^{c}}$ defined by$${\displaystyle \psi ^{c}(y):=\inf _{x}(c(x,y)-\psi (x))}$$This generalizes Legendre transformation, which is the case where ${\displaystyle c(x,y)=-xy}$ with a sign flip.

The c-convexification of a curve in the case where ${\displaystyle c(x,y)=|x-y|}$.

${\displaystyle \leq 1}$.

We say that a function ${\displaystyle \psi }$ is c-convex iff ${\displaystyle \psi =\varphi ^{c}}$ for some ${\displaystyle \varphi }$. Note that because ${\displaystyle \varphi ^{ccc}=\varphi ^{c}}$, we can always assume that ${\displaystyle \varphi }$ is c-convex. The c-convexification of a function ${\displaystyle \psi }$ is ${\displaystyle \psi ^{cc}}$. Equivalently, it is the smallest c-convex function ${\displaystyle \psi '}$ such that ${\displaystyle \psi '\geq \psi }$ pointwise.^{[10]: Prop. 5.8} Like in the case of convex transformation, ${\displaystyle \psi }$ is c-convex iff ${\displaystyle \psi =\psi ^{cc}}$.

If ${\displaystyle \psi =\varphi ^{c}:X\to \mathbb {R} }$ is c-convex, then the set of c-subdifferential of ${\displaystyle \psi }$ at ${\displaystyle x\in X}$ is the set of ${\displaystyle y\in Y}$ such that ${\displaystyle \psi (x)=c(x,y)-\varphi (y)}$. Similarly for ${\displaystyle Y}$.

When ${\displaystyle X=Y}$, the graph ${\displaystyle (y,\psi ^{c}(y))}$ can be constructed as follows: Take the graph of ${\displaystyle \psi }$, and flip it upside down. At each point ${\displaystyle (x,-\psi (x))}$, construct a graph of ${\displaystyle y\mapsto c(x,y)}$ apexed at ${\displaystyle (x,-\psi (x))}$. That is, it is the graph of ${\displaystyle y\mapsto c(x,y)-\psi (x)}$. We obtain a whole set of such graphs. Their lower-edge envelope is the graph of ${\displaystyle \psi ^{c}}$.

In the same image, we can see what it means for a function ${\displaystyle \varphi (y)}$ to be c-convex. It is c-convex iff its entire graph can be "touched" by a "tipped tool" that is moving and shape-shifting. When the tipped tool is at ${\displaystyle x}$, it has a shape of ${\displaystyle y\mapsto c(x,y)}$ and is raised to a height of ${\displaystyle -\psi (x)}$. The graph of the c-convexification ${\displaystyle \varphi ^{cc}(y)}$ is constructed by running the tipped tool so that it is lowered as much as possible, while still touching graph of ${\displaystyle \varphi (y)}$ on the upper side. The lower envelope swept out by the tipped tool is the graph of ${\displaystyle \varphi ^{cc}(y)}$.^{[10]: Fig. 5.2}

For example, if ${\displaystyle X=Y=\mathbb {R} ^{n}}$ is a metric space and ${\displaystyle c(x,y)=\|x-y\|}$, then ${\displaystyle \varphi :X\to \mathbb {R} }$ is c-convex iff it is 1-Lipschitz. This is used in the definition of 1-Wasserstein distance. If ${\displaystyle c(x,y)=\|x-y\|^{2}}$, then ${\displaystyle \varphi }$ is c-convex iff its graph could be touched from above by a tipped tool with the shape of a paraboloid.

Existence and uniqueness

Under fairly permissive assumptions, optimal transport plan exists.

If

• ${\displaystyle (X,\mu _{X}),(Y,\mu _{Y})}$ are Polish probability spaces,
• ${\displaystyle c:X\times Y\to \mathbb {R} \cup \{\infty \}}$ is lower semicontinuous,
• and there exists some upper semicontinuous functions ${\displaystyle a\in L^{1}(\mu _{X}),b\in L^{1}(\mu _{Y})}$ of type ${\displaystyle a:X\to \mathbb {R} \cup \{-\infty \},\;b:Y\to \mathbb {R} \cup \{-\infty \}}$ such that ${\displaystyle c(x,y)\geq a(x)+b(y)}$,

then an optimal transport plan exists. That is, exists ${\displaystyle \gamma ^{*}\in \Gamma (\mu _{X},\mu _{Y})}$ such that it reaches the infimum.^{[10]: Thm. 4.1}

Note that the infimum could be infinite if all transport plans turn out to be infinite. For example, if ${\displaystyle X=Y=\mathbb {R} ,c(x,y)=|x-y|,\;\mu _{X}}$ is the Cauchy distribution, and ${\displaystyle \mu _{Y}=\delta _{0}}$.

If

• ${\displaystyle (X,\mu _{X}),(Y,\mu _{Y})}$ are Polish probability spaces,
• ${\displaystyle c:X\times Y\to \mathbb {R} }$ is lower semicontinuous,
• there exists some upper semicontinuous functions ${\displaystyle a\in L^{1}(\mu _{X}),b\in L^{1}(\mu _{Y})}$ of type ${\displaystyle a:X\to \mathbb {R} ,\;b:Y\to \mathbb {R} }$ such that ${\displaystyle c(x,y)\geq a(x)+b(y)}$,
• there exists a finite-cost transport plan,
• and for any c-convex function ${\displaystyle \psi :X\to \mathbb {R} \cup \{\infty \}}$, for ${\displaystyle \mu _{X}}$-almost all ${\displaystyle x\in X}$, ${\displaystyle \psi }$ has a unique c-subdifferential at ${\displaystyle x}$

then an optimal transport map exists.^{[10]: Thm. 5.30}

A restriction of an optimal transport plan is still optimal. That is, suppose ${\displaystyle \gamma \in \Gamma (\mu ,\nu )}$ is optimal, and ${\displaystyle 0<\gamma '<\gamma }$, and define the normalized transport plan ${\displaystyle {\bar {\gamma }}':=\gamma '/\gamma '(X\times Y)}$, then ${\displaystyle {\bar {\gamma }}'}$ is an optimal transport plan between its own marginals.^{[10]: Thm. 4.6} If ${\displaystyle {\bar {\gamma }}'}$ isn't optimal, then there exists an improvement of it, which then translates back to an improvement of the original ${\displaystyle \gamma }$.

Kantorovich duality

The Kantorovich duality states that:^{[10]: Thm. 5.10}

If ${\displaystyle (X,\mu _{X}),(Y,\mu _{Y})}$ are Polish probability spaces, ${\displaystyle c:X\times Y\to \mathbb {R} \cup \{\infty \}}$ is lower semicontinuous, and there exists some upper semicontinuous functions ${\displaystyle a\in L^{1}(\mu _{X}),b\in L^{1}(\mu _{Y})}$ of type ${\displaystyle a:X\to \mathbb {R} ,\;b:Y\to \mathbb {R} }$ such that ${\displaystyle c(x,y)\geq a(x)+b(y)}$, then$${\displaystyle \inf _{\gamma \in \Gamma (\mu ,\nu )}\left(\int _{X\times Y}c(x,y)\,\mathrm {d} \gamma (x,y)\right)=\sup _{\varphi {\text{ is }}c{\text{-convex}}}\left(\int _{X}\varphi (x)\,\mathrm {d} \mu (x)+\int _{Y}\varphi ^{c}(y)\,\mathrm {d} \nu (y)\right)}$$If furthermore, ${\displaystyle c}$ only takes real values, there exists a transport plan with finite cost, and there exists some functions ${\displaystyle a'\in L^{1}(\mu _{X}),b'\in L^{1}(\mu _{Y})}$ such that ${\displaystyle c(x,y)\leq a'(x)+b'(y)}$, then$${\displaystyle \min _{\gamma \in \Gamma (\mu ,\nu )}\left(\int _{X\times Y}c(x,y)\,\mathrm {d} \gamma (x,y)\right)=\max _{\varphi {\text{ is }}c{\text{-convex}}}\left(\int _{X}\varphi (x)\,\mathrm {d} \mu (x)+\int _{Y}\varphi ^{c}(y)\,\mathrm {d} \nu (y)\right)}$$

Consider the second case, where we can actually arrive at an exactly optimal plan, instead of merely getting closer and closer. In this case, an optimal transport plan ${\displaystyle \gamma \in \Gamma (\mu ,\nu )}$, constrains the form of an optimal pricing pair ${\displaystyle (\varphi ,\psi )}$, and vice versa.

Given such an optimal pricing pair ${\displaystyle (\varphi ,\psi )}$,^{[10]: Remark 5.13}

• given an arbitrary transport plan ${\displaystyle \gamma \in \Gamma (\mu ,\nu )}$, if all ${\displaystyle (x,y)\in \operatorname {supp} (\gamma )}$ satisfies the exact equality ${\displaystyle c(x,y)=\varphi (x)+\psi (y)}$, then ${\displaystyle \gamma }$ is an optimal plan;
• given an optimal transport plan ${\displaystyle \gamma \in \Gamma (\mu ,\nu )}$, any ${\displaystyle (x,y)\in \operatorname {supp} (\gamma )}$ must satisfy the exact equality ${\displaystyle c(x,y)=\varphi (x)+\psi (y)}$.

More succinctly, a transport plan is optimal iff it is supported on the set of c-subdifferential pairs of ${\displaystyle (\varphi ,\psi )}$.

Stability

The optimal transportation is stable in the following sense:^{[10]: Thm. 5.20}

Assume that ${\displaystyle (X,\mu ),(Y,\nu )}$ are Polish probability spaces, ${\displaystyle c:X\times Y\to \mathbb {R} }$ is continuous, and ${\displaystyle \inf c}$ is finite. Given a sequence of continuous functions ${\displaystyle c:X\times Y\to \mathbb {R} }$ converging uniformly to ${\displaystyle c}$ over ${\displaystyle X\times Y}$, a sequence ${\displaystyle \mu _{k}\to \mu }$ weakly, a sequence ${\displaystyle \nu _{k}\to \nu }$ weakly, and a sequence of optimal transport plans ${\displaystyle \gamma _{k}\in \Gamma (\mu _{k},\nu _{k})}$. If the transport costs ${\displaystyle \int c_{k}d\pi _{k}}$ satisfy ${\displaystyle \int c_{k}d\pi _{k}<+\infty ,\;\forall k}$ and ${\displaystyle \liminf _{k}\int c_{k}d\pi _{k}<+\infty }$, then ${\displaystyle \gamma _{k}}$ converges weakly to some ${\displaystyle \gamma }$, and ${\displaystyle \gamma }$ is an optimal transport plan from ${\displaystyle \mu }$ to ${\displaystyle \nu }$.

Similarly, the optimal transport map is also stable.^{[10]: Cor. 5.23}

Assume that ${\displaystyle (X,\mu ),(Y,\nu )}$ are Polish probability spaces, ${\displaystyle X}$ is locally compact, ${\displaystyle c:X\times Y\to \mathbb {R} }$ is lower semicontinuous, and ${\displaystyle \inf c}$ is finite. Given a sequence of lower semicontinuous functions ${\displaystyle c:X\times Y\to \mathbb {R} }$ converging uniformly to ${\displaystyle c}$ over ${\displaystyle X\times Y}$, a sequence ${\displaystyle \nu _{k}\to \nu }$ weakly,

Economic interpretation

The optimal transport problem has an economic interpretation.^[11] Cédric Villani recounts the following interpretation from Luis Caffarelli:^[12]

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.

In the interpretation, the duality transformation transforms a loading cost function ${\displaystyle \varphi (x)}$ into the optimal (for the shipper) unloading cost function ${\displaystyle \psi (y)=\varphi ^{c}(y)}$. If the unloading cost function ${\displaystyle \psi (y)}$ were any higher at any point, then there would be some route ${\displaystyle x\to y}$ on which ${\displaystyle c(x,y)<\varphi (x)+\psi (y)}$, meaning that there is some route on which you would rather ship yourself. But if the unloading cost function were any lower at any point, then the shipper could have earned more money by raising the price there. Therefore, the shipper should always choose ${\displaystyle \psi =\varphi ^{c}}$. The same argument applied again then states that the shipper should always choose ${\displaystyle \varphi =\psi ^{c}}$, and therefore we obtain the lower bound half of the duality formula:$${\displaystyle \inf _{\gamma \in \Gamma (\mu ,\nu )}\left(\int _{X\times Y}c(x,y)\,\mathrm {d} \gamma (x,y)\right)\geq \sup _{\varphi {\text{ is }}c{\text{-convex}}}\left(\int _{X}\varphi (x)\,\mathrm {d} \mu (x)+\int _{Y}\varphi ^{c}(y)\,\mathrm {d} \nu (y)\right),}$$The Kantorovich duality states that it is in fact an equality, i.e. the shipper can make you pay as much as you would pay yourself, though the shipper might never exactly reach the bound (thus the use of infimum and supremum, instead of minimum and maximum).

Assume that the shipper in fact must pay the same cost function and us, and can exactly reach the maximum revenue using ${\displaystyle (\varphi ,\psi )}$ as their pricing chart. Then the shipper must use an optimal plan, at which point the shipper just breaks even with no profit. Conversely, any shipping plan that allows the shipper to exactly break even must be optimal.

Solution of the problem

Optimal transportation on the real line

Optimal transportation matrix

Continuous optimal transport

For ${\displaystyle 1\leq p<\infty }$, let ${\displaystyle {\mathcal {P}}_{p}(\mathbb {R} )}$ denote the collection of probability measures on ${\displaystyle \mathbb {R} }$ that have finite ${\displaystyle p}$-th moment. Let ${\displaystyle \mu ,\nu \in {\mathcal {P}}_{p}(\mathbb {R} )}$ and let ${\displaystyle c(x,y)=h(x-y)}$, where ${\displaystyle h:\mathbb {R} \to [0,\infty )}$ is a convex function.

1. If ${\displaystyle \mu }$ has no atom, i.e., if the cumulative distribution function ${\displaystyle F_{\mu }:\mathbb {R} \to [0,1]}$ of ${\displaystyle \mu }$ is a continuous function, then ${\displaystyle F_{\nu }^{-1}\circ F_{\mu }:\mathbb {R} \to \mathbb {R} }$ is an optimal transport map. It is the unique optimal transport map if ${\displaystyle h}$ is strictly convex.
2. We have

${\displaystyle \min _{\gamma \in \Gamma (\mu ,\nu )}\int _{\mathbb {R} ^{2}}c(x,y)\,\mathrm {d} \gamma (x,y)=\int _{0}^{1}c\left(F_{\mu }^{-1}(s),F_{\nu }^{-1}(s)\right)\,\mathrm {d} s.}$

The proof of this solution appears in Rachev & Rüschendorf (1998).^[13]

Discrete version and linear programming formulation

In the case where the margins ${\displaystyle \mu }$ and ${\displaystyle \nu }$ are discrete, let ${\displaystyle \mu _{x}}$ and ${\displaystyle \nu _{y}}$ be the probability masses respectively assigned to ${\displaystyle x\in \mathbf {X} }$ and ${\displaystyle y\in \mathbf {Y} }$, and let ${\displaystyle \gamma _{xy}}$ be the probability of an ${\displaystyle xy}$ assignment. The objective function in the primal Kantorovich problem is then

${\displaystyle \sum _{x\in \mathbf {X} ,y\in \mathbf {Y} }\gamma _{xy}c_{xy}}$

and the constraint ${\displaystyle \gamma \in \Gamma (\mu ,\nu )}$ expresses as

${\displaystyle \sum _{y\in \mathbf {Y} }\gamma _{xy}=\mu _{x},\forall x\in \mathbf {X} }$

and

${\displaystyle \sum _{x\in \mathbf {X} }\gamma _{xy}=\nu _{y},\forall y\in \mathbf {Y} .}$

In order to input this in a linear programming problem, we need to vectorize the matrix ${\displaystyle \gamma _{xy}}$ by either stacking its columns or its rows, we call ${\displaystyle \operatorname {vec} }$ this operation. In the column-major order, the constraints above rewrite as

${\displaystyle \left(1_{1\times |\mathbf {Y} |}\otimes I_{|\mathbf {X} |}\right)\operatorname {vec} (\gamma )=\mu }$ and ${\displaystyle \left(I_{|\mathbf {Y} \|}\otimes 1_{1\times |\mathbf {X} |}\right)\operatorname {vec} (\gamma )=\nu }$

where ${\displaystyle \otimes }$ is the Kronecker product, ${\displaystyle 1_{n\times m}}$ is a matrix of size ${\displaystyle n\times m}$ with all entries of ones, and ${\displaystyle I_{n}}$ is the identity matrix of size ${\displaystyle n}$. As a result, setting ${\displaystyle z=\operatorname {vec} (\gamma )}$, the linear programming formulation of the problem is

${\displaystyle {\begin{aligned}&{\text{Minimize }}&&\operatorname {vec} (c)^{\top }z\\[4pt]&{\text{subject to:}}&&z\geq 0,\\[4pt]&&&{\begin{pmatrix}1_{1\times |\mathbf {Y} |}\otimes I_{|\mathbf {X} |}\\I_{|\mathbf {Y} |}\otimes 1_{1\times |\mathbf {X} |}\end{pmatrix}}z={\binom {\mu }{\nu }}\end{aligned}}}$

which can be readily inputted in a large-scale linear programming solver (see chapter 3.4 of Galichon (2016)^[11]).

Semi-discrete case

In the semi-discrete case, ${\displaystyle X=Y=\mathbb {R} ^{d}}$ and ${\displaystyle \mu }$ is a continuous distribution over ${\displaystyle \mathbb {R} ^{d}}$, while ${\displaystyle \nu =\sum _{j=1}^{J}\nu _{j}\delta _{y_{i}}}$ is a discrete distribution which assigns probability mass ${\displaystyle \nu _{j}}$ to site ${\displaystyle y_{j}\in \mathbb {R} ^{d}}$. In this case, we can see^[14] that the primal and dual Kantorovich problems respectively boil down to:

${\displaystyle \inf \left\{\int _{X}\sum _{j=1}^{J}c(x,y_{j})\,d\gamma _{j}(x),\gamma \in \Gamma (\mu ,\nu )\right\}}$

for the primal, where ${\displaystyle \gamma \in \Gamma (\mu ,\nu )}$ means that ${\displaystyle \int _{X}d\gamma _{j}(x)=\nu _{j}}$ and ${\displaystyle \sum _{j}d\gamma _{j}(x)=d\mu (x)}$, and:

${\displaystyle \sup \left\{\int _{X}\varphi (x)d\mu (x)+\sum _{j=1}^{J}\psi _{j}\nu _{j}:\psi _{j}+\varphi (x)\leq c(x,y_{j})\right\}}$

for the dual, which can be rewritten as:

${\displaystyle \sup _{\psi \in \mathbb {R} ^{J}}\left\{\int _{X}\inf _{j}\left\{c(x,y_{j})-\psi _{j}\right\}d\mu (x)+\sum _{j=1}^{J}\psi _{j}\nu _{j}\right\}}$

which is a finite-dimensional convex optimization problem that can be solved by standard techniques, such as gradient descent.

In the case when ${\displaystyle c(x,y)=|x-y|^{2}/2}$, one can show that the set of ${\displaystyle x\in \mathbf {X} }$ assigned to a particular site ${\displaystyle j}$ is a convex polyhedron. The resulting configuration is called a power diagram.^[15]

Quadratic normal case

Assume the particular case ${\displaystyle \mu ={\mathcal {N}}(0,\Sigma _{X})}$, ${\displaystyle \nu ={\mathcal {N}}(0,\Sigma _{Y})}$, and ${\displaystyle c(x,y)=|y-Ax|^{2}/2}$ where ${\displaystyle A}$ is invertible. One then has

${\displaystyle \varphi (x)=-x^{\top }\Sigma _{X}^{-1/2}\left(\Sigma _{X}^{1/2}A^{\top }\Sigma _{Y}A\Sigma _{X}^{1/2}\right)^{1/2}\Sigma _{X}^{-1/2}x/2}$

${\displaystyle \psi (y)=-y^{\top }A\Sigma _{X}^{1/2}\left(\Sigma _{X}^{1/2}A^{\top }\Sigma _{Y}A\Sigma _{X}^{1/2}\right)^{-1/2}\Sigma _{X}^{1/2}Ay/2}$

${\displaystyle T(x)=(A^{\top })^{-1}\Sigma _{X}^{-1/2}\left(\Sigma _{X}^{1/2}A^{\top }\Sigma _{Y}A\Sigma _{X}^{1/2}\right)^{1/2}\Sigma _{X}^{-1/2}x}$

The proof of this solution appears in Galichon (2016).^[11]

Separable Hilbert spaces

Let ${\displaystyle X}$ be a separable Hilbert space. Let ${\displaystyle {\mathcal {P}}_{p}(X)}$ denote the collection of probability measures on ${\displaystyle X}$ that have finite ${\displaystyle p}$-th moment; let ${\displaystyle {\mathcal {P}}_{p}^{r}(X)}$ denote those elements ${\displaystyle \mu \in {\mathcal {P}}_{p}(X)}$ that are Gaussian regular: if ${\displaystyle g}$ is any strictly positive Gaussian measure on ${\displaystyle X}$ and ${\displaystyle g(N)=0}$, then ${\displaystyle \mu (N)=0}$ also.

Let ${\displaystyle \mu \in {\mathcal {P}}_{p}^{r}(X)}$, ${\displaystyle \nu \in {\mathcal {P}}_{p}(X)}$, ${\displaystyle c(x,y)=|x-y|^{p}/p}$ for ${\displaystyle p\in (1,\infty ),p^{-1}+q^{-1}=1}$. Then the Kantorovich problem has a unique solution ${\displaystyle \kappa }$, and this solution is induced by an optimal transport map: i.e., there exists a Borel map ${\displaystyle r\in L^{p}(X,\mu ;X)}$ such that

${\displaystyle \kappa =(\mathrm {id} _{X}\times r)_{*}(\mu )\in \Gamma (\mu ,\nu ).}$

Moreover, if ${\displaystyle \nu }$ has bounded support, then

${\displaystyle r(x)=x-|\nabla \varphi (x)|^{q-2}\,\nabla \varphi (x)}$

for ${\displaystyle \mu }$-almost all ${\displaystyle x\in X}$ for some locally Lipschitz, ${\displaystyle c}$-concave and maximal Kantorovich potential ${\displaystyle \varphi }$. (Here ${\displaystyle \nabla \varphi }$ denotes the Gateaux derivative of ${\displaystyle \varphi }$.)

By minimizing flows

A gradient descent formulation for the solution of the Monge–Kantorovich problem was given by Sigurd Angenent, Steven Haker, and Allen Tannenbaum.^[16]

Entropic regularization

Consider a variant of the discrete problem above, where we have added an entropic regularization term to the objective function of the primal problem

${\displaystyle {\begin{aligned}&{\text{Minimize }}\sum _{x\in \mathbf {X} ,y\in \mathbf {Y} }\gamma _{xy}c_{xy}+\varepsilon \gamma _{xy}\ln \gamma _{xy}\\[4pt]&{\text{subject to: }}\\[4pt]&\gamma \geq 0\\[4pt]&\sum _{y\in \mathbf {Y} }\gamma _{xy}=\mu _{x},\forall x\in \mathbf {X} \\[4pt]&\sum _{x\in \mathbf {X} }\gamma _{xy}=\nu _{y},\forall y\in \mathbf {Y} \end{aligned}}}$

One can show that the dual regularized problem is

${\displaystyle \max _{\varphi ,\psi }\sum _{x\in \mathbf {X} }\varphi _{x}\mu _{x}+\sum _{y\in \mathbf {Y} }\psi _{y}v_{y}-\varepsilon \sum _{x\in \mathbf {X} ,y\in \mathbf {Y} }\exp \left({\frac {\varphi _{x}+\psi _{y}-c_{xy}}{\varepsilon }}\right)}$

where, compared with the unregularized version, the "hard" constraint in the former dual (${\displaystyle \varphi _{x}+\psi _{y}-c_{xy}\geq 0}$) has been replaced by a "soft" penalization of that constraint (the sum of the ${\displaystyle \varepsilon \exp \left((\varphi _{x}+\psi _{y}-c_{xy})/\varepsilon \right)}$ terms). The optimality conditions in the dual problem can be expressed as

Eq. 5.1: ${\displaystyle \mu _{x}=\sum _{y\in \mathbf {Y} }\exp \left({\frac {\varphi _{x}+\psi _{y}-c_{xy}}{\varepsilon }}\right)~\forall x\in \mathbf {X} }$

Eq. 5.2: ${\displaystyle \nu _{y}=\sum _{x\in \mathbf {X} }\exp \left({\frac {\varphi _{x}+\psi _{y}-c_{xy}}{\varepsilon }}\right)~\forall y\in \mathbf {Y} }$

Denoting ${\displaystyle A}$ as the ${\displaystyle |\mathbf {X} |\times |\mathbf {Y} |}$ matrix of term ${\displaystyle A_{xy}=\exp \left(-c_{xy}/\varepsilon \right)}$, solving the dual is therefore equivalent to looking for two diagonal positive matrices ${\displaystyle D_{1}}$ and ${\displaystyle D_{2}}$ of respective sizes ${\displaystyle |\mathbf {X} |}$ and ${\displaystyle |\mathbf {Y} |}$, such that ${\displaystyle D_{1}AD_{2}1_{|\mathbf {Y} |}=\mu }$ and ${\displaystyle (D_{1}AD_{2})^{\top }1_{|\mathbf {X} |}=\nu }$. The existence of such matrices generalizes Sinkhorn's theorem and the matrices can be computed using the Sinkhorn–Knopp algorithm,^[17] which simply consists of iteratively looking for ${\displaystyle \varphi _{x}}$ to solve Equation 5.1, and ${\displaystyle \psi _{y}}$ to solve Equation 5.2. Sinkhorn–Knopp's algorithm is therefore a coordinate descent algorithm on the dual regularized problem.

Applications

The Monge–Kantorovich optimal transport has found applications in wide range in different fields. Among them are:

• Image registration and warping^[18]
• Reflector design^[19]
• Retrieving information from shadowgraphy and proton radiography^[20]
• Seismic tomography and reflection seismology^[21]
• The broad class of economic modelling that involves gross substitutes property (among others, models of matching and discrete choice).

See also

• Wasserstein metric
• Transport function
• Hungarian algorithm
• Transportation planning
• Earth mover's distance
• Monge–Ampère equation