In 1781, Gaspard Monge, geometer of the French crown, asked how to move a pile of earth into an excavation of the same volume at the least total cost of transport. Every grain must travel, every grain pays some cost for its journey, and the question is which grain goes where. It is an exercise in logistics, and it turned out to be interesting geometrically: Monge extracted some lovely geometry about routes that never cross and then the problem sat, essentially open, for a hundred and sixty years.
It was worth the wait. The mathematics that eventually grew out of the earth-moving question now prices goods in economics, defines curvature on spaces that are too rough for ordinary calculus, trains generative models, reconstructs the histories of embryonic cells and of the early universe, and collects prizes: a Nobel for Kantorovich in 1975 for the idea that cracked it, Fields medals for Villani in 2010 and Figalli in 2018 for what it became. My aim in this essay is to build up the subject from nothing but an undergraduate understanding of math, and to show you where each of the payoffs comes from.
The subjects of this essay lie on the square below. Parts I and II build a distance between probability distributions out of Monge's transport cost, and find a linear program that has associated with it a shadow economy of prices. From that corner the square offers two moves. Add time, and the distance becomes a principle of least action for moving crowds (Part III). Add heat, and the brittle optimal plan melts into a blur that a sixty-year-old scaling algorithm computes at an impressive speed (Part IV). Make both moves at once and we land on a question Schrödinger asked in 1931, about the most probable history of a cloud of particles caught misbehaving, and the other three corners assemble themselves into its answer (Part V). You can click underlined terms for definitions and derivations, and interact with the figures to get intuition for the machinery that this essay introduces.
Let's start where Monge did, with a cost. There is stuff distributed as \(\mu\) (the pile), a desired distribution \(\nu\) (the hole), and a price \(c(x,y)\) for hauling one unit of stuff from \(x\) to \(y\); throughout this essay the price will be the squared distance \(c=\|x-y\|^2\), a choice that looks arbitrary now and become better motivated by Part III. Monge asked for a map: an assignment \(T\) sending each point of the pile to a point of the hole, moving \(\mu\) onto \(\nu\), at least total cost
$$\min_{T:\ T_\#\mu=\nu}\ \int c\big(x,\,T(x)\big)\,\mathrm d\mu(x),$$
where says that the map really does turn the pile into the hole. The cost is a sum over grains, which is as tame as objectives get. All of the trouble lives in the constraint: the set of maps pushing \(\mu\) onto \(\nu\) is a horrible, non-convex thing, small perturbations of a valid \(T\) are invalid, and there is no calculus to do. This is why the problem stalled for so long, and it stalled in good company; the Academy offered a prize for it, and the young Kantorovich did not know of Monge at all when he broke it from the other side.
\(T_\#\mu = \nu\) means \(\nu(A) = \mu(T^{-1}(A))\) for every set \(A\): the amount of mass \(T\) delivers into \(A\) is the amount \(\nu\) demands there. For densities in \(\mathbb R^d\) and nice enough \(T\) this is the change-of-variables formula \(\rho_\mu(x) = \rho_\nu(T(x))\,|\!\det DT(x)|\), whose nonlinearity in \(T\) is the technical face of everything I just complained about.
To see the structure of the question with nothing in the way, let's make everything finite for a moment: \(n\) bakeries, \(n\) cafés, one truckload each, and an \(n\times n\) matrix of haulage prices. A map is now a permutation, the problem is the assignment problem, and you can play with it below. There are \(n!\) permutations, so brute force stops being computable at around \(n=15\), and it is worth knowing that this particular wall was demolished: , solvable in polynomial time, and its study is entangled with the birth of linear programming itself.
Leonid Kantorovich published the relaxation you are about to meet in 1942, in a three-page note, having already invented linear programming in 1939 to schedule Soviet plywood production. George Dantzig built the simplex method in 1947 without knowledge of it; Koopmans connected the two worlds; Kantorovich and Koopmans shared the 1975 Nobel prize in economics. Transport is the continuous linear program, the one the founders kept in their heads as the meaning of the subject.
For the finite case, specialized algorithms beat generic simplex methods: the Hungarian method of Kuhn (1955, assembled from lemmas of Kőnig and Egerváry, hence the name) solves the \(n\times n\) assignment problem in \(O(n^3)\), and transport with unequal weights is a minimum-cost flow problem, one of the workhorse models of operations research. If you have ever seen matching markets analyzed in economics, kidney exchanges, or ad auctions, you have seen this LP wearing work clothes.
The figure below has a toggle between two costs, and they disagree in a way worth understanding before you start dragging points around. One is the plain distance \(c(x,y)=\|x-y\|\), the un-squared one, what a truck's odometer would charge. The other is the squared distance \(c(x,y)=\|x-y\|^2\), which punishes a long haul quadratically and is the cost we'll commit to in this essay. Under one of them the optimal routes provably never cross, under the other they sometimes do, and .
Take two shipments, \(x_1\to y_1\) and \(x_2\to y_2\), and ask whether swapping their destinations would be cheaper. Whatever the cost, the plan can only be optimal if the answer is no, for every pair. Run that test for each of our two costs.
Plain distance. Suppose the two segments cross, at a point \(z\). Then each route's length splits at \(z\), and we can regroup the four pieces:
$$\|x_1-y_1\|+\|x_2-y_2\| \;=\; \underbrace{\big(\|x_1-z\|+\|z-y_2\|\big)}_{\ge\ \|x_1-y_2\|} \;+\; \underbrace{\big(\|x_2-z\|+\|z-y_1\|\big)}_{\ge\ \|x_2-y_1\|},$$
each bracket bounded below by the triangle inequality. So the swapped pairing costs no more, and strictly less unless all four points sit on one line. An optimal matching for the distance cost therefore has no crossings at all, in any dimension. This is the geometry Monge found, and his cost was distance.
Squared distance. Now the same comparison is pure algebra. Expand all four squares and watch everything except the cross terms cancel:
$$\big[\|x_1-y_1\|^2+\|x_2-y_2\|^2\big] - \big[\|x_1-y_2\|^2+\|x_2-y_1\|^2\big] \;=\; -\,2\,(x_1-x_2)\cdot(y_1-y_2).$$
So the swap fails to help precisely when \((x_1-x_2)\cdot(y_1-y_2)\ge 0\): supply points that are far apart in some direction should send their mass far apart in that same direction. That is a condition on an inner product, which is a weaker demand than a condition on the picture. On a line the two demands coincide, because there the inner product is a product of two numbers and "same sign" means "same order" means "no crossing", which is where Part I's quantile formula will come from. In the plane they part company: two segments can cross while their difference vectors still point the same way. About two thirds of the random six-point instances the figure generates do exactly that, so a nonzero crossing counter under squared cost is the theory working, not the solver failing.
The inner-product condition is the one that survives and generalizes. Part II extends it from pairs to cycles of any length (cyclical monotonicity), and Brenier's theorem then identifies the maps satisfying it: gradients of convex functions, the multidimensional meaning of "nondecreasing".
Drag any point; the optimal assignment is recomputed from scratch by pricing all 720 permutations, and the counter at the bottom reports how many of the chosen routes cross. Under plain distance it stays at zero forever. Under squared distance, which is what the toggle starts on, drag for a while and you will catch crossings: circled where they happen, optimal anyway, and explained in the box above.
Where Monge's formulation fails: One bakery, two cafés, each wanting half the bread: no map exists, because a map must send the bakery somewhere, entire. Mass sometimes has to split, and maps cannot split. Kantorovich's move, one of the great acts of lateral thinking in applied mathematics, is to stop asking where each grain goes and ask instead how much goes from here to there. The unknown becomes a plan (formally a coupling between measures): a joint distribution \(\pi(x,y)\) whose first marginal is the supply and whose second is the demand,
$$\min_{\pi\,\in\,\Pi(\mu,\nu)}\ \int c(x,y)\,\mathrm d\pi(x,y), \qquad \Pi(\mu,\nu)=\Big\{\pi\ :\ \textstyle\int_y \pi(\cdot,y)=\mu,\ \int_x \pi(x,\cdot)=\nu\Big\}.$$
Look at what has happened to the geometry of the problem. The objective is linear in \(\pi\), the constraints are linear equalities, the feasible set is convex and never empty (the independent coupling \(\mu\otimes\nu\) always qualifies), and minima of linear functions over convex sets are found at extreme points. Our problem has become a linear program. And the relaxation costs nothing: , so allowing mass to split never helps unless splitting was forced, as at the lone bakery. You give up the map formulation and the map comes back as a theorem.
With \(n\) unit supplies and \(n\) unit demands, plans are doubly stochastic matrices (nonnegative, all row and column sums one). Birkhoff's theorem: the extreme points of this polytope are exactly the permutation matrices, every doubly stochastic matrix being a convex average of them. A linear objective attains its minimum at an extreme point, hence at a permutation. The general principle, that linear programs are honest about combinatorial problems whenever the constraint polytope has integral vertices, grew into the field of polyhedral combinatorics; transport is its founding example.
With the plan formulation in hand, let's define the central object of the essay. For \(p\ge 1\), the Wasserstein distance of order \(p\) between two distributions with finite \(p\)-th moments is
$$W_p(\mu,\nu) \;=\; \Big(\min_{\pi\in\Pi(\mu,\nu)} \int \|x-y\|^p \,\mathrm d\pi\Big)^{1/p},$$
the cheapest transport bill, read as a distance. It , and it is a different animal from the pointwise comparisons of densities you may know. Kullback–Leibler and its relatives compare \(\mu\) and \(\nu\) vertically, ratio of heights at each point, and become infinite once the supports stop overlapping. \(W_p\) compares the two densitities horizontally: two narrow spikes a millimeter apart are close in \(W_p\) because moving one onto the other is cheap, while their KL divergence is infinite. A distance that incorporates the geometry of the underlying space, that metrizes weak convergence, and that is finite between anything and anything is a distance you can do analysis with, and this property is behind most of its adoptions: it is why \(W_1\) was taken for training GANs whose generator starts on a low-dimensional manifold, and why comparing point clouds, images, and shapes is often done with a Wasserstein distance.
Symmetry: transpose the plan. Identity: \(W_p(\mu,\mu)=0\) via the diagonal coupling, and positivity because a zero-cost plan is supported on \(x=y\). The triangle inequality is the one with content: given good plans \(\pi_{12}\) and \(\pi_{23}\), the gluing lemma builds a three-variable coupling with those two as marginals, and Minkowski's inequality on \(\|x-z\|\le\|x-y\|+\|y-z\|\) finishes it. Gluing is intuitively obvious (chain the shipments through the middle distribution) and slightly fussy to prove; e.g. Villani's book (linked at the end) does it in brief.
One case can be solved by hand, and it anchors everything: dimension one. On the line, with \(c=(x-y)^2\), the best plan never lets two routes cross, because . A plan with no crossings preserves order, so it must ship the leftmost fraction of supply to the leftmost fraction of demand, quantile to quantile: with \(F\) and \(G\) the cumulative distribution functions,
Suppose a plan ships \(x_1\to y_2\) and \(x_2\to y_1\) with \(x_1<x_2\), \(y_1<y_2\) (a crossing). Swap the destinations. The change in cost is
$$\big[(x_1-y_1)^2+(x_2-y_2)^2\big]-\big[(x_1-y_2)^2+(x_2-y_1)^2\big] \;=\; -2\,(x_2-x_1)(y_2-y_1)\;<\;0,$$
which is the one-dimensional case of the algebra we did for the figure above, now with the extra information that \(x_2-x_1\) and \(y_2-y_1\) are both positive numbers, so the whole expression is negative. Every crossing is therefore a saving left on the table, and the optimal plan is monotone.
This is where the line is special. Above, in the plane, monotone was weaker than non-crossing; here they are the same condition, because "the difference vectors point the same way" is nothing but "the order is preserved". And a monotone plan is unique once the marginals are fixed, since there is only one order-preserving way to pour \(\mu\) into \(\nu\): the leftmost \(q\) fraction of the pile has to become the leftmost \(q\) fraction of the hole, for every \(q\). That is the quantile coupling, mass at quantile level \(q\) of \(\mu\) going to quantile level \(q\) of \(\nu\), and integrating the squared displacement \(|F^{-1}(q)-G^{-1}(q)|^2\) over all levels \(q\in[0,1]\) gives the formula below.
Keep the swap argument in view. Run it over cycles of every length instead of pairs and it becomes cyclical monotonicity, the subject of Part II; the one-dimensional case is the whole theory in miniature.
$$W_2^2(\mu,\nu) \;=\; \int_0^1 \big|F^{-1}(q) - G^{-1}(q)\big|^2 \,\mathrm dq,$$
which is an explicit formula with no optimization left to be done. Sorting solves transport on the line, which is why the 1D figures in this essay can show you exact optima and not approximations. In higher dimension there is no sorting and no formula, and the honest routes to the optimum are the LP (this part), the dual (next part), or the two deformations of the problem that the rest of the essay is about.
Every linear program has a dual, and the dual of transport is a small economics lesson. Suppose you own the bakeries and the cafés and are dreading the haulage. A shipping company makes an offer: it will buy bread at bakery \(x\) for a price \(\varphi(x)\), handle everything in between, and sell bread at café \(y\) for a price \(\psi(y)\). You will accept only if no route is a ripoff, that is, only if
$$\varphi(x) + \psi(y) \;\le\; c(x,y) \qquad \text{for all pairs } (x,y),$$
since otherwise it'd be cheaper for you to transport bread between that pair yourself at cost \(c\). The company, subject to keeping you on board, maximizes its revenue \(\int \varphi\,\mathrm d\mu + \int \psi\,\mathrm d\nu\). Revenue can never exceed the true haulage cost (integrate the constraint against any feasible plan). The theorem, Kantorovich duality, is that a clever enough company loses nothing:
$$\min_{\pi\in\Pi(\mu,\nu)} \int c\,\mathrm d\pi \;=\; \max_{\varphi(x)+\psi(y)\le c(x,y)} \;\int \varphi\,\mathrm d\mu + \int \psi\,\mathrm d\nu,$$
haulage and pricing are the same number, computed by . The prices \(\varphi,\psi\) are the Lagrange multipliers of the marginal constraints, the shadow value of one extra unit of supply here or demand there, and economists use exactly this machinery on matching markets, from labor to marriage models.
In the finite case this is textbook LP duality: dualize the two families of marginal constraints, and the multipliers are \(\varphi_i\) and \(\psi_j\). In the continuum one route is a minimax swap: write the constrained minimum as \(\min_{\pi\ge0}\max_{\varphi,\psi}\big[\langle c,\pi\rangle + \langle\varphi,\mu-\pi\mathbf 1\rangle + \langle\psi,\nu-\pi^{\mathsf T}\mathbf 1\rangle\big]\) and exchange min and max; making the exchange rigorous is a compactness argument (or Fenchel–Rockafellar duality), and holds in vast generality.
Given \(\psi\), the best feasible \(\varphi\) is forced: \(\varphi(x) = \min_y\, [\,c(x,y)-\psi(y)\,]\), called the c-transform \(\psi^{c}\). Optimal price pairs are fixed points of transforming back and forth, exactly as convex functions are fixed points of the double Legendre transform. File the shape of this update away: in Part IV an algorithm will do nothing except apply a softened c-transform, alternately, until it converges.
The dual is where structure enters, through the complementary slackness of the LP: at the optimum, mass flows only where the price is tight, that is, the optimal plan lives on the set where \(\varphi(x)+\psi(y)=c(x,y)\). On every used route the company breaks even against the haulage; the slack routes carry nothing. From tightness, one differentiation away, comes a geometric fingerprint of optimality: , meaning that rerouting any finite cycle of shipments can only raise the total cost. The figure lets you feel this with your hands.
A set \(S\) of routes is c-cyclically monotone if for every finite cycle \((x_1,y_1),\dots,(x_k,y_k)\in S\),
$$\sum_{i=1}^k c(x_i,y_i) \;\le\; \sum_{i=1}^k c(x_i,y_{i+1}),\qquad y_{k+1}:=y_1,$$
every cyclic reshuffle of destinations is a loss. Optimal plans are supported on such sets (if a profitable cycle existed you could reroute a little mass around it and beat the optimum), and, for reasonable costs, the converse holds too: cyclical monotonicity plus the marginals pins down optimality. The quantifier over all cycle lengths matters. Pairwise swaps alone can be locally unbeatable while a three-cycle still saves money, which you can occasionally verify below by getting stuck.
Click one supply point, then another: their destinations swap. Starting from the optimum, every swap you can make raises the bill, which is cyclical monotonicity for cycles of length two. Scramble and try to walk back down by swaps alone; if you get stuck above the optimum, you have met a longer cycle, and the definition's quantifier over all cycle lengths has earned its keep.
Now let's see why we're using a quadratic cost. Take \(c(x,y)=\tfrac12\|x-y\|^2\); the one half halves every bill and moves no minimizer, it is there to keep the algebra clean. Write out the price constraint we just imposed, expanding the square:
$$\varphi(x)+\psi(y) \;\le\; \tfrac12\|x-y\|^2 \;=\; \tfrac12\|x\|^2 \;-\; x\cdot y \;+\; \tfrac12\|y\|^2 .$$
Both quadratic pieces on the right depend on one variable each, so each can be moved to the left and absorbed into the price living on that side. Define
$$u(x) \;=\; \tfrac12\|x\|^2-\varphi(x), \qquad w(y) \;=\; \tfrac12\|y\|^2-\psi(y),$$
which is a change of unknowns and nothing more: from \(u\) and \(w\) we can read the prices back off, and vice versa. In the new variables the constraint becomes
$$u(x)+w(y) \;\ge\; x\cdot y \quad \text{for all } x,y, \qquad\text{that is,}\qquad u(x) \;\ge\; x\cdot y - w(y) \quad \text{for every } y.$$
Read the second form one \(x\) at a time. Fix a destination \(y\). The function \(x \mapsto x\cdot y - w(y)\) is affine in \(x\): a straight line in one dimension, a hyperplane in general, whose slope is the vector \(y\) and whose height is set by the number \(w(y)\). So the constraint says: the graph of \(u\) lies above each member of a whole family of hyperplanes, one for every possible destination.
Which member is binding, and how high does \(u\) sit? Recall what the company wants. Its revenue grows with \(\varphi\), and \(\varphi = \tfrac12\|x\|^2 - u\), so making the company richer means pushing \(u\) down, as far as the constraint allows. It allows exactly down to the highest hyperplane in the family. At the optimum, therefore,
$$u(x) \;=\; \sup_y\ \big[\,x\cdot y \;-\; w(y)\,\big] \;=\; w^{*}(x),$$
the of \(w\). Two things follow immediately, and they are the whole theorem.
For a function \(f\) on \(\mathbb R^d\), its Legendre transform is
$$f^{*}(y) \;=\; \sup_x\ \big[\,x\cdot y - f(x)\,\big].$$
Here is what the supremum is measuring. Fix a slope \(y\) and slide the linear function \(x\mapsto x\cdot y\) downward until it stops touching the graph of \(f\) from below. How far you had to slide is \(f^{*}(y)\), and the line you end up with is \(x\mapsto x\cdot y - f^{*}(y)\), the highest affine minorant of \(f\) with that slope. So \(f^{*}\) is a catalogue: for every slope, the intercept of the best line with that slope. Convex functions are precisely the ones this catalogue can reconstruct, which is the content of the involution \(f^{**}=f\) (for convex, lower semicontinuous \(f\)).
Two properties are the ones we use above. First, \(f^{*}\) is a supremum of functions that are affine in \(y\), and a supremum of affine functions is convex no matter what \(f\) was, which is why the \(u\) built in the main text is convex without anyone having assumed it. Second, for smooth strictly convex \(f\) the gradients invert each other, \(\nabla f^{*} = (\nabla f)^{-1}\), and the supremum defining \(f^{*}(y)\) is attained at the \(x\) where \(\nabla f(x)=y\). That attainment is exactly the tightness condition the main text turns into the transport map.
Examples worth carrying: \(f(x)=\tfrac12\|x\|^2\) has \(f^{*}=f\), the self-dual case; \(f(x)=\tfrac12 x^{\mathsf T}\!Ax\) has \(f^{*}(y)=\tfrac12 y^{\mathsf T}\!A^{-1}y\); \(f(x)=e^{x}\) has \(f^{*}(y)=y\log y - y\), which is where entropy comes from and why the transform keeps showing up in Part IV.
You have probably met this transform already, wearing other hats. It is what takes a Lagrangian to a Hamiltonian in mechanics, trading velocity for momentum, and what takes energy to free energy in thermodynamics, trading entropy for temperature. Each time the move is the same: describe a convex function not by its values but by its slopes.
And it is the quadratic-cost special case of Part II's c-transform \(\psi^{c}(x)=\min_y[c(x,y)-\psi(y)]\). Replace the family of hyperplanes \(x\cdot y - \text{const}\) by the family of cost functions \(x \mapsto c(x,y) - \text{const}\), and everything above goes through with \"convex\" replaced by \"c-concave\". The quadratic cost is the case where the cost-shaped family is the hyperplane family, which is why it, and only it, hands us convexity.
First, \(u\) is a supremum of affine functions, hence convex, and nobody assumed convexity anywhere. It fell out of the price constraint.
Second, remember complementary slackness from a moment ago: mass moves only along tight routes. Tightness at a used route \((x,y)\) means \(u(x)+w(y)=x\cdot y\), which says that of all the hyperplanes in the family, the one indexed by this particular \(y\) is the one that actually touches the graph of \(u\) at this particular \(x\). A hyperplane touching a convex function from below at a point is called a supporting hyperplane, and its slope is called a subgradient there. So the destination \(y\) of the mass sitting at \(x\) is a subgradient of \(u\) at \(x\).
That last step is where the splitting problem dies. A convex function can have several subgradients at a point, but only at its kinks, and a convex function on \(\mathbb R^d\) has kinks only on a set of measure zero. So if \(\mu\) has a density, then at \(\mu\)-almost every \(x\) the function \(u\) is differentiable, the subgradient is one single vector, and the mass at \(x\) has exactly one destination:
$$y \;=\; \nabla u(x).$$
Kantorovich let mass split; the prices have just told us it does not, and where it goes. This is (1991), the structural heart of the quadratic theory: the optimal plan is a map, and the map is the gradient of a convex function,
$$T \;=\; \nabla u, \qquad u \text{ convex}.$$
Statement: if \(\mu\) has a density (no mass concentrated on small sets, so no bakery is ever forced to split), the quadratic problem has a unique optimal plan, it is induced by a map \(T\) which is the gradient of a convex potential, and conversely any gradient of a convex function pushing \(\mu\) to \(\nu\) is optimal. Uniqueness and structure in one stroke, from bookkeeping about prices.
The main text derived the forward direction. The converse is worth a sentence because it is the one you use in practice, to certify that a map you have built is the right one. Gradients of convex functions satisfy \((\nabla u(x_1)-\nabla u(x_2))\cdot(x_1-x_2)\ge0\), and more generally the cyclic version of that inequality; by the criterion of the previous box, a plan supported on such a set cannot be improved by any reshuffle, so it is optimal. Monotone map, optimal plan, and convex potential are three descriptions of one object.
Where do the hypotheses bite? Drop the density assumption on \(\mu\) and the argument breaks exactly where you would expect: an atom of mass sitting on a kink of \(u\) has several subgradients and does split, which is Kantorovich's relaxation earning its keep. Nothing is required of \(\nu\), which may be discrete; that asymmetric case has its own name, semi-discrete transport, and its own algorithms.
If we insert \(T=\nabla u\) into the change-of-variables formula, the mass balance becomes a PDE for the potential, the Monge–Ampère equation: \(\det D^2u(x) = \rho_\mu(x)\,/\,\rho_\nu(\nabla u(x))\), a fully nonlinear second-order equation, degenerate where densities vanish. Its regularity theory (when is the optimal map smooth? continuous? where does it tear?) is extensive: Caffarelli's work in the 1990s opened it, the map can be discontinuous outright when the target's support is non-convex, and sharpened partial regularity in this family of problems is part of what Figalli's 2018 Fields medal cites.
In one dimension, gradients of convex functions are the nondecreasing functions, and Brenier's theorem collapses to Part I's monotone rearrangement.
Pause on how strange this is. We fed in a shipping problem and a bill of prices; out came the statement that optimal rearrangements of mass are exactly the gradients of convex potentials, a class of maps with no rotation in them, no shear, no swirl, the multidimensional meaning of monotone. Convexity was nowhere in the assumptions. It condensed out of complementary slackness, the way crystal structure condenses out of a cooling liquid, and most of modern transport theory is built on this one identification. The 2D figures that follow can only be exact because Gaussians hand us their Brenier maps .
Between \(\mathcal N(m_0,\Sigma_0)\) and \(\mathcal N(m_1,\Sigma_1)\) the Brenier map is affine, \(T(x) = m_1 + A\,(x-m_0)\) with
$$A \;=\; \Sigma_0^{-1/2}\big(\Sigma_0^{1/2}\Sigma_1\Sigma_0^{1/2}\big)^{1/2}\Sigma_0^{-1/2},$$
symmetric positive definite (hence a gradient of the convex quadratic \(\tfrac12 x^{\mathsf T}\!Ax + \dots\)), and the distance comes out as \(W_2^2 = \|m_0-m_1\|^2 + \operatorname{tr}\!\big(\Sigma_0+\Sigma_1 - 2(\Sigma_0^{1/2}\Sigma_1\Sigma_0^{1/2})^{1/2}\big)\), the Bures metric of quantum information under another name. Figure 3 runs exactly this formula.
Before we move on, one more thing duality provides us. For \(p=1\) the two prices collapse into one: , which is the form in which transport entered deep learning.
For \(c=\|x-y\|\), feasibility forces \(\psi=-\varphi\) with \(\varphi\) 1-Lipschitz, and duality reads \(W_1(\mu,\nu) = \sup_{\mathrm{Lip}(\varphi)\le1} \int\varphi\,\mathrm d\mu - \int\varphi\,\mathrm d\nu\). The Wasserstein GAN (Arjovsky, Chintala, Bottou 2017) trains a neural network to play \(\varphi\), the critic, against a generator minimizing the estimated \(W_1\); the horizontal character of the distance is the point, since a generator supported on a thin manifold gets useful gradients from \(W_1\) where likelihood-style objectives saturate. The Lipschitz constraint is enforced only loosely in practice, and an honest reading is that WGANs are transport-flavored rather than transport-exact; the flavor was enough to matter.
So far the transport plan is a ledger: it says what leaves \(x\) for \(y\) and says nothing about the trip. When I took a class on optimal transport, I used to joke with friends that transport plans are "optimal teleport", not optimal transport. The square in the introduction has a downward arrow labelled "add time", and this part is that arrow: we move from the top-left corner, where the answer is a ledger, to the bottom-left corner, where the answer is a motion. Fix a unit of time and ask for the whole history: a curve of densities \(\rho_t\) with \(\rho_0=\mu\), \(\rho_1=\nu\).
There is a lazy way to interpolate that you should look at once and then discard forever: the blend \((1-t)\,\mu + t\,\nu\), which fades the old density out and the new one in. Nothing travels; mass at the source dwindles while mass at the target swells, as if teleported. The transport answer instead carries each grain along a road. Take the Brenier map \(T\) and walk every grain a fraction \(t\) of its journey:
$$\rho_t \;=\; \big((1-t)\,\mathrm{id} + t\,T\big)_{\#}\,\mu,$$
McCann's displacement interpolation (1997). The figure puts the two side by side, and the difference is not subtle: one is a crossfade, the other is traffic. is a story of its own, and it is the door through which transport entered pure geometry.
Some functionals of a density, the entropy \(\int\rho\log\rho\) foremost, fail to be convex along blends yet are convex along displacement interpolations. McCann introduced the interpolation to prove uniqueness results with exactly this trick, and the idea scaled astonishingly: on a curved space, how convex the entropy is along these paths measures Ricci curvature. Lott, Sturm, and Villani turned that observation into a definition of curvature bounds needing no derivatives at all, valid on metric measure spaces where classical Ricci curvature is meaningless. That line of work is a large part of Villani's 2010 Fields citation, and it started with watching entropy along moving mass.
Two Gaussians, exact formulas on both panels. Drag either blob by its center, in either panel. On the left, mass dies here and is born there, passing through nowhere; scrub slowly and watch the corridor between them stay empty. On the right every particle walks a straight line at constant speed, and the traced dots are ten of them doing it.
Displacement interpolation looks like mechanics, and the resemblance can be stated as a theorem. Benamou and Brenier (2000) proved that the static distance \(W_2\) can be equivalently formulated as the solution of a least-action problem over all conceivable traffic: among every density curve \(\rho_t\) and velocity field \(v_t\) that move \(\mu\) to \(\nu\) while conserving mass,
$$W_2^2(\mu,\nu) \;=\; \min_{(\rho,v)}\ \int_0^1\!\!\int \rho_t(x)\,\|v_t(x)\|^2\,\mathrm dx\,\mathrm dt \qquad \text{s.t.}\quad \partial_t\rho_t + \nabla\!\cdot\!(\rho_t v_t)=0,\quad \rho_0=\mu,\ \rho_1=\nu,$$
the total kinetic energy of the crowd, minimized subject to the . The is two inequalities, and the optimizer is exactly the displacement interpolation: the squared distance of Part I is the least kinetic energy needed to physically execute the move in unit time. This is why the exponent \(p=2\) was the right choice all along; squared distance is what kinetic energy charges for straight-line travel.
\(\partial_t\rho + \nabla\cdot(\rho v)=0\) is local conservation of mass. Take any small region: the amount of stuff inside it can change only by stuff crossing the boundary, and the flux carrying it across is \(\rho v\), density times velocity. Summing the crossings over the boundary and shrinking the region is what the divergence does, and the equation says the accumulation inside is exactly minus that. Nothing is created, nothing is destroyed, and this is the same equation that keeps the books for charge in electromagnetism and for probability in fluid form.
Its role in the formula above is to be the promise that \(v\) really carries \(\rho\), rather than being some unrelated arrow field we drew on top of a movie of densities. One technical remark worth having: the equation is linear in the pair \((\rho,\ m)\) where \(m=\rho v\) is the momentum field, and in those variables the kinetic energy \(\int\!\!\int \|m\|^2/\rho\) is jointly convex. So the whole least-action problem is a convex problem in disguise, which is what let Benamou and Brenier solve it numerically at all.
We want to show that two numbers are equal: the least kinetic energy \(\mathcal K\) over all traffic obeying the constraints, and the static \(W_2^2\) from Part I. As usual with a claimed equality between two optimization problems, we show each is at most the other.
Why the objective is called kinetic energy. A single particle of mass \(m\) and velocity \(v\) carries \(\tfrac12 m\|v\|^2\). Here the "mass" sitting at \(x\) at time \(t\) is \(\rho_t(x)\,\mathrm dx\) and its velocity is \(v_t(x)\), so the crowd's instantaneous energy is \(\tfrac12\int\rho_t\|v_t\|^2\,\mathrm dx\) and the action is its time integral. The factor \(\tfrac12\) is dropped by convention, which is why the formula lands on \(W_2^2\) rather than \(\tfrac12W_2^2\).
Step one: \(\mathcal K \le W_2^2\). Exhibit one legal traffic pattern with that energy. Take the Brenier map \(T\) from Part II and let every grain drive straight from \(x\) to \(T(x)\) at constant speed, arriving at time one. Its position at time \(t\) is \(X_t(x)=(1-t)x+t\,T(x)\) and its velocity is \(T(x)-x\), the same at all times. The density of the crowd at time \(t\) is then \(\rho_t=(X_t)_\#\mu\), which is McCann's displacement interpolation, and the pair \((\rho_t,v_t)\) built this way satisfies the continuity equation for free, since it was made by physically moving mass around. Its energy is
$$\int_0^1\!\!\int \rho_t\|v_t\|^2\,\mathrm dx\,\mathrm dt \;=\; \int_0^1\!\!\int \|T(x)-x\|^2\,\mathrm d\mu(x)\,\mathrm dt \;=\; \int \|T(x)-x\|^2\,\mathrm d\mu(x) \;=\; W_2^2(\mu,\nu).$$
The middle equality is the change of variables \(y=X_t(x)\): integrating a quantity against the crowd at time \(t\) is the same as integrating it against the original crowd, since the mass at \(X_t(x)\) came from \(x\). The time integral then does nothing, because the speed of each grain does not depend on \(t\). So legal traffic with energy \(W_2^2\) exists, and the minimum can only be lower.
Step two: \(\mathcal K \ge W_2^2\). Now take any legal traffic \((\rho_t,v_t)\) and show it pays at least \(W_2^2\). Follow the velocity field: let \(X_t(x)\) solve \(\dot X_t = v_t(X_t)\) with \(X_0(x)=x\). The continuity equation is precisely the statement that this flow carries the crowd, \(\rho_t = (X_t)_\#\mu\), so \(X_1\) pushes \(\mu\) to \(\nu\). Then the pair \(\pi=(X_0,X_1)_\#\mu\), read as "the grain that started at \(x\) ended at \(X_1(x)\)", is a feasible transport plan, and by definition of the minimum in Part I,
$$W_2^2(\mu,\nu) \;\le\; \int \|X_1(x)-x\|^2\,\mathrm d\mu(x).$$
It remains to bound that displacement by the energy. Each grain's net displacement is the integral of its velocity, \(X_1(x)-x=\int_0^1 v_t(X_t(x))\,\mathrm dt\), and squared norms do not commute with averaging: by Jensen's inequality (equivalently Cauchy–Schwarz on the time interval),
$$\Big\|\int_0^1 v_t(X_t)\,\mathrm dt\Big\|^2 \;\le\; \int_0^1 \big\|v_t(X_t)\big\|^2\,\mathrm dt .$$
Integrate over the crowd and change variables back, exactly as in step one, to recognize the right-hand side as the kinetic energy of the traffic we started with. Chaining the three displays gives \(W_2^2\le \mathcal K\), which is what we wanted.
What equality costs. Jensen is an equality only when the thing being averaged is constant, so a traffic pattern achieving the minimum must move every grain at a velocity that never changes in time: straight lines, constant speed. That pins the optimizer down to the displacement interpolation and nothing else, which is the sentence the main text was waiting for. In words, the crowd that spends the least energy is the crowd in which nobody accelerates and nobody swerves.
The paper also provided an algorithm: solved as a convex problem in \((\rho,\rho v)\) by an augmented Lagrangian scheme, it was the standard way to compute large transport problems for a decade, and it earned the authors a central place in OT.
The minimizer has more structure than feasibility requires. At the optimum the velocity is the gradient of a potential, \(v_t=\nabla\phi_t\), and the latter solves the \(\partial_t\phi + \tfrac12\|\nabla\phi\|^2 = 0\); the particles of the optimal crowd feel no force and never collide. Read as a control problem, which is what it is: steer a crowd from \(\mu\) so that it lands exactly on \(\nu\) at time one, paying \(\int\!\!\int\rho\|v\|^2\) for effort, with no randomness anywhere. \(\phi\) is the value function of that problem, its characteristics are the straight roads of Figure 3, and this reading of \(W_2\) as the cost of deterministic optimal control is the lower left corner of the square, which we'll diffuse in Part V.
Fix the traffic density \(\rho_t\) and minimize kinetic energy over velocity fields carrying it: any component of \(v\) orthogonal (in \(L^2(\rho)\)) to gradients stirs mass without changing \(\rho\)'s evolution, pure waste, so optimal velocities are gradients \(v=\nabla\phi\). Inserting \(v=\nabla\phi\) and taking the variation of \(\int\!\!\int\rho\,\|\nabla\phi\|^2\) subject to continuity, the multiplier computation yields \(\partial_t\phi+\tfrac12\|\nabla\phi\|^2=0\) along the optimum (Burgers' equation for \(\nabla\phi\)): each particle keeps its initial velocity, i.e. free streaming, a pressureless gas. The pair (continuity for \(\rho\), Hamilton–Jacobi for \(\phi\)) is the geodesic equation of an infinite-dimensional geometry, which is the next paragraph.
Benamou–Brenier licenses a change of self-image for the whole subject. If \(W_2\) is a least action of the kinetic type, then \(W_2\) is a Riemannian distance: the space of probability densities becomes a manifold whose tangent vectors at \(\rho\) are the velocity fields \(v=\nabla\phi\), with kinetic energy \(\int\rho\|v\|^2\) as the metric, displacement interpolations as its geodesics. Otto worked out this calculus around 2000, and its famous dividend is thermodynamic: . The heat equation, of all things, is the steepest descent of the Boltzmann entropy with respect to the transport metric, and dissipative PDE after dissipative PDE turned out to be a gradient flow of a free energy in the same sense. That single sentence is a bridge between this essay and statistical mechanics, and I will cross it once more, from the other side, at the end of Part V.
Jordan, Kinderlehrer, and Otto (1998) discretized time and showed that the scheme
$$\rho_{k+1} \;=\; \arg\min_{\rho}\ \Big[\ \int \rho\log\rho \;+\; \frac{1}{2\tau}\,W_2^2(\rho,\rho_k)\ \Big]$$
converges, as \(\tau\to0\), to the heat equation \(\partial_t\rho=\Delta\rho\): each step trades entropy gained against transport spent, and the continuum limit of that bargain is diffusion. With an energy \(\int V\rho\) added to the entropy the same scheme produces the Fokker–Planck equation, so the dissipative dynamics of statistical mechanics are gradient flows in the Benamou–Brenier geometry. Otto's formal Riemannian calculus turns this into a computational device (geodesics, Hessians, curvature of the density manifold), and it is the reason "Wasserstein gradient flow" is now a standard phrase in PDE, sampling theory, and the analysis of learning dynamics.
Now back to the static case; which has three practical issues. The LP is expensive: algorithms that are cubic-ish in time stall around a few thousand points and are therefore hopeless for millions of pixels or cells. The optimal plan is brittle: supported on a lower-dimensional set (a permutation, a map), it jumps around discontinuously when the data wiggles even slightly, and estimating so singular an object from samples is statistically cursed. And the LP value is a piecewise-linear function of its inputs, so gradients (i.e. learning) give only very localized information. The repair for all three of these problems is the same and it is thermodynamically flavored: stop demanding the cheapest plan and ask instead for the cheapest warm plan,
$$\pi_\varepsilon \;=\; \arg\min_{\pi\in\Pi(\mu,\nu)}\ \int c\,\mathrm d\pi \;+\; \varepsilon\,\mathrm{KL}\big(\pi\,\|\,\mu\otimes\nu\big),$$
haulage plus \(\varepsilon\) times the against the independent coupling \(\mu\otimes\nu\), the plan that ships from everywhere to everywhere in proportion to supply and demand and pays no attention to distance. Relative entropy measures how far a plan has moved away from that reference, so the second term is a fee for being opinionated, and it fattens whatever the first term wants to sharpen. Turning \(\varepsilon\) down recovers the linear program; turning it up drives the plan back to \(\mu\otimes\nu\).
Calling \(\varepsilon\) a temperature is not decoration. The objective is an energy \(\langle c,\pi\rangle\) plus temperature times negative entropy, which is a free energy, this time over shipping ledgers rather than over microstates, and free energies have Boltzmann minimizers. Setting the derivative to zero, which the box below does line by line, gives the
First the object. For two distributions \(\pi\) and \(R\) over the same set of routes,
$$\mathrm{KL}(\pi\|R) \;=\; \sum_{i,j} \pi_{ij}\,\log\frac{\pi_{ij}}{R_{ij}} \;\ge\; 0,$$
with equality only when \(\pi=R\). The usual reading is in bits: it is the excess description length you pay if you compress data drawn from \(\pi\) using a code designed for \(R\). For our purposes the shape matters more than the interpretation. It is finite only when \(\pi\) puts mass nowhere \(R\) does not, it blows up logarithmically as \(\pi\) concentrates, and, since \(t\mapsto t\log t\) is strictly convex, \(\mathrm{KL}\) is strictly convex in \(\pi\). That last property is why the warm problem has exactly one solution while the linear program can have a flat face full of them.
Now the computation. We minimize
$$\sum_{i,j} c_{ij}\pi_{ij} \;+\; \varepsilon \sum_{i,j}\pi_{ij}\log\frac{\pi_{ij}}{\mu_i\nu_j}$$
over nonnegative \(\pi\) with row sums \(\mu_i\) and column sums \(\nu_j\). Attach a multiplier \(\varphi_i\) to each row constraint and \(\psi_j\) to each column constraint, exactly as in Part II, and write the Lagrangian
$$L(\pi,\varphi,\psi) \;=\; \sum_{i,j} c_{ij}\pi_{ij} + \varepsilon \sum_{i,j}\pi_{ij}\log\frac{\pi_{ij}}{\mu_i\nu_j} \;-\; \sum_i \varphi_i\Big(\sum_j \pi_{ij}-\mu_i\Big) \;-\; \sum_j \psi_j\Big(\sum_i \pi_{ij}-\nu_j\Big).$$
Differentiate with respect to one entry \(\pi_{ij}\). The first term contributes \(c_{ij}\). The entropy term contributes \(\varepsilon\big(\log\frac{\pi_{ij}}{\mu_i\nu_j}+1\big)\), the \(+1\) coming from differentiating \(t\log t\). The two constraint terms contribute \(-\varphi_i\) and \(-\psi_j\). Setting the sum to zero:
$$c_{ij} + \varepsilon\log\frac{\pi_{ij}}{\mu_i\nu_j} + \varepsilon - \varphi_i - \psi_j \;=\; 0 .$$
Solve for \(\pi_{ij}\) by isolating the logarithm and exponentiating:
$$\pi_{ij} \;=\; \mu_i\,\nu_j\ \exp\!\Big(\frac{\varphi_i+\psi_j-c_{ij}}{\varepsilon}-1\Big) \;=\; \mu_i\,\nu_j\ e^{\varphi_i/\varepsilon}\,e^{-c_{ij}/\varepsilon}\,e^{\psi_j/\varepsilon},$$
where the stray \(e^{-1}\) has been absorbed by shifting \(\varphi\) by the constant \(\varepsilon\), which no constraint notices. This is the displayed Gibbs form: a fixed kernel \(e^{-c_{ij}/\varepsilon}\) built from the cost alone, multiplied by one positive number per row and one per column.
Two observations before you leave the box. The exponential is automatically positive, so the constraint \(\pi\ge0\) never has to be enforced; the entropy term has quietly removed the inequality constraints that make the linear program combinatorial. And the pattern here, "linear cost plus \(\varepsilon\) times relative entropy to a reference, minimized under linear constraints, yields an exponential tilt of the reference", is the same computation that produces the Boltzmann distribution in statistical mechanics, with energy playing the role of \(c\) and the multiplier on the normalization playing the role of the free energy. Here the physical system is a shipping ledger; in Part V it will turn out to be an actual gas.
Gibbs form of the solution: with \(K = e^{-c/\varepsilon}\) the kernel of the cost,
$$\pi_\varepsilon \;=\; \operatorname{diag}(a)\ K\ \operatorname{diag}(b) \qquad \text{for some positive vectors } a, b.$$
Sit with what this says about degrees of freedom. An \(n\times n\) plan had \(n^2\) unknowns; the warm plan is a fixed known matrix \(K\) corrected by \(2n\) numbers, one price per bakery and one per café (the multipliers \(a=e^{\varphi/\varepsilon}\), \(b=e^{\psi/\varepsilon}\) are Part II's dual prices, exponentiated). And the two constraints suggest their own algorithm, so obvious it has been discovered at least three times: fix \(b\), then matching the row sums to \(\mu\) determines \(a\) in closed form; fix \(a\), match the column sums, get \(b\); repeat. The two updates are one line each,
$$a_i \;\leftarrow\; \frac{\mu_i}{(Kb)_i}, \qquad\qquad b_j \;\leftarrow\; \frac{\nu_j}{(K^{\mathsf T}a)_j},$$
and that is the entire algorithm. It is , and each half-step turns out to be Part II's c-transform with the hard minimum softened into a log-sum-exp at temperature \(\varepsilon\): the duality of Part II, running as a while-loop. Watch it work:
The updates. We know from the previous box that the answer looks like \(\pi_{ij}=a_i K_{ij} b_j\) with \(K_{ij}=\mu_i\nu_j e^{-c_{ij}/\varepsilon}\) folded together, and we need the row sums to be \(\mu\) and the column sums \(\nu\). Impose the row constraint on that form:
$$\sum_j a_i K_{ij} b_j \;=\; \mu_i \quad\Longleftrightarrow\quad a_i \,(Kb)_i \;=\; \mu_i \quad\Longleftrightarrow\quad a_i \;=\; \frac{\mu_i}{(Kb)_i}.$$
Since \(a_i\) appears in row \(i\) and nowhere else, the row constraints do not interact, and one division fixes all of them at once. The same argument on the columns gives \(b_j=\nu_j/(K^{\mathsf T}a)_j\). Each update makes its own family of constraints exact and, in general, breaks the other family; alternate, and both errors shrink. At the fixed point both hold and you are looking at the unique warm plan.
Why this is a softened c-transform. Undo the exponentials by writing \(a_i=e^{\varphi_i/\varepsilon}\) and \(b_j=e^{\psi_j/\varepsilon}\), so the unknowns are the prices themselves. The row update becomes
$$\varphi_i \;\leftarrow\; -\,\varepsilon\,\log \sum_j \nu_j\, e^{(\psi_j-c_{ij})/\varepsilon}.$$
Look at what happens as \(\varepsilon\to0\). Let \(m_i=\min_j\,[c_{ij}-\psi_j]\), factor \(e^{-m_i/\varepsilon}\) out of the sum, and the remaining terms are \(e^{-(\text{something}\ \ge\ 0)/\varepsilon}\), which die unless they belong to the minimizing \(j\). What survives is \(\varphi_i\to \min_j[c_{ij}-\psi_j]\), which is exactly the c-transform of Part II. So Sinkhorn is the c-transform iteration with \(\min\) replaced by a soft minimum at temperature \(\varepsilon\), and \(\varepsilon\) controls how soft. Implementations run it in this logarithmic form, because \(e^{-c/\varepsilon}\) underflows to zero long before \(\varepsilon\) is as small as one would like.
Why it converges. Each half-step is an exact maximization of the entropic dual over one block of variables with the other held fixed, and that dual is smooth and concave, so the loop is textbook block coordinate ascent and cannot get stuck. The sharp statement is stronger and older: the iteration is a contraction in Hilbert's projective metric, so the error decays geometrically, with a rate that degrades like \(e^{-\mathrm{osc}(c)/\varepsilon}\) as the temperature falls (Franklin–Lorenz 1989). That is the theory behind the iteration counter in the figure climbing as you cool it: sharpness is bought with sweeps, and the exchange rate is exponential.
Three births. Statisticians invented the loop in 1940 as iterative proportional fitting (Deming–Stephan), to force census tables to agree with known margins. Sinkhorn proved the matrix-scaling convergence theorems in the 1960s. And Fortet ran essentially the same iteration in 1940 to prove existence for the Schrödinger system you will meet in Part V, decades before anyone knew these were the same fixed point. The modern era starts with Cuturi (2013), who noticed the updates are two matrix-vector products, embarrassingly parallel, GPU-native, and differentiable end to end, whereupon transport went from a specialist's tool to a layer you can drop into a network.
The square is the plan \(\pi_\varepsilon(x,y)\), computed live by Sinkhorn; the dashed curve is the exact \(\varepsilon=0\) optimum from Part I's quantile formula. Cool the temperature and watch the blur condense onto the curve while the iteration counter, the price of sharpness, climbs. Heat it and the plan forgets the cost and drifts toward the independent product. Then draw your own marginals and watch the row and column sums being alternately snapped into place.
What did the blur buy? Speed, first: the figure's plan on a \(96\times96\) grid resolves in milliseconds, and the same loop runs on millions of points on a GPU, which no exact LP solver approaches. Uniqueness and stability, second: the warm objective is strictly convex, the plan depends smoothly on the data, and gradients of the entropic cost with respect to anything (the marginals, the point positions, parameters upstream) exist and backpropagate. Third, and less advertised, : the exact distance estimated from \(n\) samples converges at the cursed rate \(n^{-1/d}\), while entropic quantities achieve \(n^{-1/2}\) with the dimension demoted into the constant. The blur is regularization in the statistical sense too, and in high dimension it is arguably the only version of the distance whose measurement from data means anything.
Empirical measures of \(n\) samples in \(\mathbb R^d\) satisfy \(\mathbb E\,W_2(\hat\mu_n,\mu) \asymp n^{-1/d}\) for \(d>4\) (Dudley's argument: a sample cloud cannot fill space faster than that), so in dimension 50 the plain distance between empirical measures reflects sample noise long before it reflects \(\mu\) and \(\nu\). For fixed \(\varepsilon>0\) the entropic cost concentrates at the parametric rate \(n^{-1/2}\) with constants growing as \(\varepsilon\) shrinks (Genevay, Chizat, Bach, Cuturi, Peyré 2019), interpolating the two regimes. One wart: the entropic cost of \(\mu\) against itself is not zero; the repair \(S_\varepsilon(\mu,\nu) = \mathrm{OT}_\varepsilon(\mu,\nu) - \tfrac12\mathrm{OT}_\varepsilon(\mu,\mu) - \tfrac12\mathrm{OT}_\varepsilon(\nu,\nu)\), the Sinkhorn divergence, is nonnegative, zero only at equality, and metrizes weak convergence, a fully respectable discrepancy that happens to be computable.
The adoption list is long and I will keep only the vivid entries. Single-cell biology: cells are destroyed by measurement, so a developing tissue is seen as unpaired snapshots at successive times, and entropic transport between consecutive snapshots reconstructs which ancestors likely became which descendants; the Waddington-OT pipeline (Schiebinger et al. 2019) did this for tens of thousands of reprogramming cells and is by now a standard trick in the field. Graphics moves color palettes between photographs along transport plans. Economists fit matching markets with entropic transport because the blur models unobserved heterogeneity, the Gibbs form is their logit choice model. And cosmologists, in the other temporal direction, reconstructed the early universe's density field by solving one gigantic assignment problem between today's galaxy positions and a uniform primordial soup (Frisch et al. 2002), Monge's earth at the largest available scale.
There is one more thing the blur bought, and it is the strangest. Tilt your head and the warm plan stops looking like a numerical compromise: a Gibbs distribution over routes, prices playing potentials, a temperature dial, an algorithm that equilibrates it. Part V shows this is no analogy. The entropically blurred plan is exactly the endpoint statistics of a cloud of physical Brownian particles, and the question it answers was asked by Schrödinger before any of the computational motivations existed.
In 1931, between inventing wave mechanics and acquiring his cat, Schrödinger published a thought experiment about hot gas. Release a huge number of independent Brownian particles, observe the cloud's density to be \(\mu\) at time zero, look away, and look back at time one to find the density is \(\nu\), where \(\nu\) is not what diffusion predicts from \(\mu\). You have witnessed something exponentially unlikely, a conspiracy of the noise. His question: given that the anomaly happened, what is the most probable history of the cloud in between? What did the gas most likely do?
The modern formulation, which makes exact, is a projection in the geometry of information: among all stochastic processes \(P\) whose time-zero law is \(\mu\) and time-one law is \(\nu\), find the one closest in relative entropy to the physics that was supposed to happen,
$$\min_{P:\ P_0=\mu,\ P_1=\nu}\ \mathrm{KL}\big(P \,\|\, W^{\varepsilon}\big),$$
where \(W^\varepsilon\) is the law of Brownian motion with diffusivity \(\varepsilon\) (started from \(\mu\)). This is the fourth corner of the square: the unknown is an entire evolving cloud, as in Part III, and the objective is an entropy relative to a noisy reference, as in Part IV. And it dismantles beautifully, in two moves.
Sanov's theorem: the empirical distribution of \(N\) independent copies of a process with law \(R\) satisfies a large deviation principle, \(\Pr[\hat P_N \approx P] \sim e^{-N\,\mathrm{KL}(P\|R)}\). Conditioning on a rare event (the observed marginals) therefore concentrates, as \(N\to\infty\), on the conditional mode: the \(P\) satisfying the constraints with least \(\mathrm{KL}(P\|R)\). "Most probable history of the cloud" is the informal reading of that concentration, and Föllmer made the connection to Schrödinger's problem precise in the 1980s. Schrödinger, writing before large deviation theory existed, got the right functional by a combinatorial Boltzmann-counting argument, which is very much the same proof in period costume.
First move: a path is an endpoint pair plus a story in between, so split the entropy . The in-between term is minimized to zero by letting particles travel exactly as the reference physics travels when its endpoints are known, which for Brownian motion means Brownian bridges, pinned random walks with variance \(\varepsilon\,t(1-t)\). All actual freedom collapses onto the endpoints: only the coupling \(\pi(x,y)\), who starts where and ends where, remains to be chosen, and its objective is
A path is two pieces of information: where it starts and ends, and what it does in between. Both our laws can be written that way. Write \(P = \pi \otimes P^{x,y}\), meaning "draw endpoints from \(\pi\), then draw a path conditioned on those endpoints from \(P^{x,y}\)", and likewise \(W^\varepsilon = R_{01}\otimes W^{x,y}\), where \(R_{01}\) is the endpoint law of the reference Brownian motion and \(W^{x,y}\) is a Brownian motion conditioned to start at \(x\) and finish at \(y\), which is called a Brownian bridge. The chain rule for relative entropy splits the objective along that seam:
$$\mathrm{KL}(P\|W^\varepsilon) \;=\; \underbrace{\mathrm{KL}(\pi\|R_{01})}_{\text{endpoints}} \;+\; \underbrace{\int \mathrm{KL}\big(P^{x,y}\,\|\,W^{x,y}\big)\,\mathrm d\pi(x,y)}_{\text{the story in between}}.$$
The second term is an average of relative entropies, so it is nonnegative, and it is zero exactly when \(P^{x,y}=W^{x,y}\) for (almost) every pair. Nothing in the constraints touches it: we only constrained the marginals at times \(0\) and \(1\), which are properties of \(\pi\) alone. So the optimal cloud travels by ordinary Brownian bridges, and the only real decision left is the endpoint coupling \(\pi\).
Now identify the problem that \(\pi\) solves. The reference endpoint law is \(R_{01}(x,y)=\mu(x)\,k_\varepsilon(x,y)\) with the heat kernel \(k_\varepsilon(x,y)=(2\pi\varepsilon)^{-d/2}e^{-\|x-y\|^2/2\varepsilon}\): start from \(\mu\), diffuse for one unit of time. Insert it and split the logarithm by multiplying and dividing by \(\nu(y)\):
$$\mathrm{KL}(\pi\|R_{01}) \;=\; \int \log\frac{\pi(x,y)}{\mu(x)\nu(y)}\,\mathrm d\pi \;+\; \int \log\frac{\nu(y)}{k_\varepsilon(x,y)}\,\mathrm d\pi .$$
The first integral is \(\mathrm{KL}(\pi\|\mu\otimes\nu)\), the entropic penalty of Part IV. In the second, \(\int\log\nu(y)\,\mathrm d\pi = \int \log\nu\,\mathrm d\nu\) depends only on the fixed marginal, and \(-\log k_\varepsilon(x,y)=\frac{\|x-y\|^2}{2\varepsilon}+\frac d2\log(2\pi\varepsilon)\) contributes the squared cost plus another constant. Collecting,
$$\mathrm{KL}(\pi\|R_{01}) \;=\; \frac{1}{2\varepsilon}\int\|x-y\|^2\,\mathrm d\pi \;+\; \mathrm{KL}(\pi\|\mu\otimes\nu) \;+\; \text{const},$$
and multiplying by \(\varepsilon\), which changes no minimizer, gives \(\tfrac12\int\|x-y\|^2\,\mathrm d\pi + \varepsilon\,\mathrm{KL}(\pi\|\mu\otimes\nu)\): the entropic transport problem of Part IV, coefficients and all. The regularization strength was never a knob we chose; it is the diffusivity of the gas, and the Gibbs kernel \(e^{-c/\varepsilon}\) was the heat kernel all along.
$$\min_{\pi\in\Pi(\mu,\nu)}\ \int \tfrac12\|x-y\|^2\,\mathrm d\pi \;+\; \varepsilon\,\mathrm{KL}\big(\pi\,\|\,\mu\otimes\nu\big)\quad(+\ \mathrm{const}),$$
Part IV verbatim, with the temperature revealed as the diffusivity of the gas and the Gibbs kernel revealed as the heat kernel. Nothing about the warm plan was a numerical metaphor: entropic transport is the endpoint bookkeeping of Schrödinger's cloud, Sinkhorn's fixed-point system was written down by Schrödinger himself (existence proved by Fortet in 1940, by running the iteration), and the algorithm of Figure 4, applied to the heat-kernel cost, computes most-likely histories of a physical gas. Second move: sample endpoints from the plan, hang Brownian bridges between them. That full object, the Schrödinger bridge, is what Figure 5 shows.
Time runs left to right between the two drawn densities; each blue thread is one particle of the most probable cloud, endpoints sampled from the Sinkhorn plan, interior filled in by a Brownian bridge. The grey lines are the zero-temperature answer, Part III's displacement interpolation. Cool ε and the threads lay themselves down onto the geodesic; heat it and the plan loosens while the bridges thrash.
The figure makes a claim in its corner that deserves its theorem. As \(\varepsilon\to0\) the Schrödinger problem : the fourth corner of the square sits directly above the third and settles onto it as the heat is withdrawn. The precise statement identifies the finite-\(\varepsilon\) problem as least-action too, for a corrected action: the bridge minimizes kinetic energy plus \(\tfrac{\varepsilon^2}{8}\) times the Fisher information of the evolving density. The extra term charges for sharpness, a penalty on \(\int\rho\,\|\nabla\log\rho\|^2\), so the hot crowd keeps itself smooth in a way the cold crowd never bothers to; you can see the reluctance to focus in the figure at high \(\varepsilon\).
Write the bridge's density flow \(\rho_t\) and its current velocity \(v_t\) (the gradient part of the probability flux). A computation going back to Yasue and made systematic by Léonard and by Gentil, Léonard, and Ripani turns \(\varepsilon\,\mathrm{KL}(P\|W^\varepsilon)\), up to boundary terms depending only on \(\mu,\nu\), into
$$\int_0^1\!\!\int \Big(\ \tfrac12\|v_t\|^2 \;+\; \tfrac{\varepsilon^2}{8}\,\big\|\nabla\log\rho_t\big\|^2\ \Big)\rho_t\,\mathrm dx\,\mathrm dt,$$
minimized over solutions of the continuity equation with pinned ends: Benamou–Brenier plus a Fisher-information regularizer. As \(\varepsilon\to0\) the problems \(\Gamma\)-converge to the kinetic action, optimal values converging to \(\tfrac12 W_2^2\) and the bridges to the displacement interpolation (Mikami 2004, Léonard 2012). So "entropic interpolation" deserves the name: it is displacement interpolation for a world that refuses to be perfectly sharp.
What does the optimal cloud do, particle by particle? The solution has a closed structural form: the bridge is a diffusion with the reference noise and a modified drift,
$$\mathrm dX_t \;=\; \varepsilon\,\nabla\log\phi_t(X_t)\,\mathrm dt \;+\; \sqrt{\varepsilon}\;\mathrm dB_t,$$
where , and the density factorizes as \(\rho_t = \phi_t\,\hat\phi_t\), one factor propagating from each endpoint constraint. The drift is worth reading aloud: the gradient of the log of a probability, the score, aimed at futures compatible with the destination. A particle of Schrödinger's cloud drifts along the score of where it must end up, and this mechanism, steering a diffusion by adding a score to its drift, is precisely how modern generative models sample: score-based diffusion models learn \(\nabla\log\) of noised data and run it as a drift, and their guided and bridge-matching variants (Diffusion Schrödinger Bridge and descendants) are numerical Schrödinger problems between the data distribution and something easy, trained with neural surrogates in place of Fortet's iterates. The 1931 gas has a very direct line to the image generators of the 2020s.
Ansatz \(\rho_t=\phi_t\hat\phi_t\) with \(\partial_t\phi = -\tfrac{\varepsilon}{2}\Delta\phi\) (backward) and \(\partial_t\hat\phi = \tfrac{\varepsilon}{2}\Delta\hat\phi\) (forward): the product then solves the Fokker–Planck equation of the displayed SDE, as one checks by differentiating. The boundary conditions \(\phi_0\hat\phi_0=\mu\), \(\phi_1\hat\phi_1=\nu\) constitute the Schrödinger system; its unknowns are two positive functions, exactly the \(a\) and \(b\) of Part IV's Gibbs form promoted to function space, and Fortet's 1940 existence proof iterates them alternately, which is Sinkhorn before Sinkhorn. Probabilistically, the drift correction \(\varepsilon\nabla\log\phi\) is Doob's h-transform: conditioning a Markov process on future information tilts its drift by the score of the conditional probability of that future. Schrödinger noted the eerie resemblance of the \((\phi,\hat\phi)\) pair to quantum mechanics' \((\psi,\bar\psi)\); that resemblance, in imaginary time, is a rabbit hole I am closing the lid on here.
And thermodynamics collects its dividend in both directions. Part III showed diffusion running downhill in the transport geometry; here is the converse trade, transport priced in dissipation: . The Wasserstein distance is the price sheet of finite-time thermodynamics: rushing a transformation costs dissipation in proportion to squared transport per unit time, a refinement of Landauer's bound with the geometry of this essay sitting inside a physical inequality. The square, which entered as a table of contents, exits as physics on every edge: its right edge heats a ledger into a gas, its bottom edge turns a distance into least action, and its corners are joined by limits that commute.
For overdamped Langevin dynamics steered by a time-dependent potential, the total entropy production of a protocol carrying \(\rho_0=\mu\) to \(\rho_\tau=\nu\) equals \(\frac{1}{T}\!\int_0^\tau\!\!\int \rho\,\|v\|^2\), the kinetic action of the realized density flow (in units with unit mobility, temperature \(T\)). Minimizing over protocols is then exactly Benamou–Brenier on a stretched clock, giving \(\Delta S_{\mathrm{tot}}^{\min} = W_2^2(\mu,\nu)\,/\,(T\tau)\) (Aurell, Mejía-Monasterio, Muratore-Ginanneschi 2011). Consequences include finite-time corrections to Landauer erasure and "thermodynamic speed limits": you cannot transform a state quickly and reversibly, and the irreversibility floor is a Wasserstein length. Optimal transport turns out to be the geometry in which nonequilibrium thermodynamics keeps its accounts.
Kantorovich priced the routes, Brenier bent the prices into a gradient, Benamou and Brenier melted the ledger into motion, and Schrödinger had asked, a decade before any of them wrote, what the ledger could mean in a world that jitters. The subject reads as if it were planned backwards from his question. Moving mass optimally and conditioning noise on its destination are the same calculation, carried out at different temperatures.
Three pieces of fine print, in increasing order of importance.
Existence is not regularity. Kantorovich plans always exist; Brenier maps need the source to have a density, and even then the map can be discontinuous, tearing mass apart wherever the target's support has dents. The regularity theory of the Monge–Ampère equation is a subject of its own, and any application that silently assumes a smooth transport map is writing checks that Caffarelli's theorems only sometimes cash.
The curse is real. In high dimension, exact Wasserstein distances between empirical measures mostly measure sample noise (\(n^{-1/d}\) does not forgive), and the honest quantities are the entropic ones, at a temperature you must choose and report. The \(\varepsilon\to0\) limit exists on paper; on data, chasing it re-buys the curse and the iteration count of Figure 4 at once. Practitioners anneal \(\varepsilon\), work in log-space for stability, and debias with Sinkhorn divergences, and all three habits concede the same point, that the blur is doing statistical work and cannot be waved away.
The cost is a modeling decision. Everything downstream of Part II leaned on the quadratic cost: Brenier, the Riemannian picture, the heat-kernel identity of Part V. Change the cost and the structure changes species: \(c=\|x-y\|\) forgets direction along rays, concave costs reward consolidating shipments, general costs on manifolds carry their own geometry. The square is the quadratic square; other costs have other essays.
Doors I pointed at without opening, each one room deep: Wasserstein barycenters, averaging several distributions at once (Agueh–Carlier), which is how one builds atlases of shapes and templates of anatomies; unbalanced transport, where mass may be created and destroyed at a price, indispensable for growing tissues and imperfect measurements; Gromov–Wasserstein, transport between spaces with no common coordinates, matching graphs to graphs and proteins to proteins by their internal geometry alone; semi-discrete transport, densities against point masses, whose optimal plans are Laguerre cell diagrams and whose solvers power the cosmological reconstruction mentioned in Part IV; the multi-marginal and mean-field versions, where many populations negotiate at once; and the ongoing project of putting neural surrogates inside every corner of the square, from input-convex networks playing Brenier potentials to the bridge-matching samplers of Part V.
If this essay did its job, the square is now furniture in your head: a distance, its prices, its motion, its heat. For the longer road, the sources below are the ones I actually recommend.