Transport

1Moving mass at minimum cost

Imagine a pile of earth that must be moved to fill a hole. The question is one of logistics: which shovelful is to be carried where, and at what total effort. If the effort of moving a unit of earth from one place to another is proportional to the distance covered, then there is a definite minimum effort for the job, and a definite assignment of source locations to target locations that achieves it. That is the problem Gaspard Monge posed in 1781 [1], in a memoir written for the French Academy of Sciences. He gave the problem in physical language because he was concerned with the actual movement of soil during fortification work, but he treated it as a mathematical question, the first systematic treatment of what is now called optimal transport.

The modern statement keeps Monge's geometry and adds one generalisation. Let the source pile be described not by a finite list of locations but by a probability measure on a space, written μ; let the destination pile be described by a second probability measure on another space, written ν. The two measures have unit total mass by convention. A transport map is a function T from the source space to the target space that, in a sense made precise below, pushes μ onto ν: the mass μ assigns to a region of the source must equal the mass ν assigns to the image of that region under T. The cost of moving a unit of mass from x to T(x) is given by a function c, often the squared distance for reasons that will become clear. The total cost of the transport is the integral of c with respect to μ. The Monge problem is to find the map T that pushes μ onto ν while making this integral as small as possible.

The difficulty of the Monge problem is not merely numerical. There are pairs of source and target for which no map T pushes μ onto ν at all. The simplest example is a source consisting of a single point and a target consisting of two points: no single-valued function from one point can produce two. More subtly, when the source is a continuous distribution and the target is a discrete one, mass must split, which a map cannot do. The Monge problem is well posed in enough cases to be interesting and ill posed in enough cases to need a generalisation. The next four sections build that generalisation and the algorithms that solve it.

Transport plan

A joint probability measure π on the product of source and target spaces, whose two marginal distributions are μ and ν. The plan specifies, for every pair of locations, how much mass moves from the first to the second. A transport map T induces a plan, namely the one supported on the graph of T. Plans are strictly more general: they allow mass at a single source to split between multiple targets, which is the technical step that lifts Monge's problem out of its existence difficulties.

The cost function in Monge's formulation is part of the data of the problem, not part of its answer, and the choice has a substantive effect on the optimum. The original problem used the Euclidean distance itself, and the optima that result are characterised by the geometric condition that the transport rays do not cross. Squared Euclidean distance, the choice that returns repeatedly in this document, gives optima whose structure is governed by convex potentials rather than crossing geometry. Costs that are not powers of distance arise in matching markets, where the price of pairing two agents need not have any relation to a distance at all; in image processing, where the cost might come from a perceptual metric; and in economics, where it is a utility differential. The mathematics is the same in each case, but the answers, and the algorithms that recover them, depend on the choice.

2The discrete assignment problem

The cleanest concrete instance of optimal transport is the assignment problem: a finite number of source points must be matched, one to one, with the same number of target points, with each pairing incurring a cost and the total cost to be minimised. With five sources and five targets the set of admissible assignments is the set of permutations of five items, which contains one hundred and twenty elements. With ten sources it contains roughly three and a half million. The set is finite for any fixed size, and the optimum is therefore attained, but its enumeration is impractical even at modest sizes, and the central computational question of the next century became how to find the optimum without examining the whole permutation set.

The figure below carries out the small case in two ways. With five source points and five target points it enumerates all hundred and twenty permutations and reports the exact optimal cost. It also provides a local-search step: the Improve button scans pairs of assignments and swaps them when the swap lowers the cost. The local search converges to a configuration in which no single swap reduces the cost, but this two-opt local optimum can fall short of the global one. The contrast between local and global optimality is one of the themes the linear-programming formulation of the next section settles definitively.

For larger sizes brute force is not available. The standard exact algorithm for the assignment problem is the Hungarian method, published by Kuhn in 1955 and improving on Egerváry and Kőnig's combinatorial work from the 1930s; it runs in time cubic in the number of points and solves the problem to optimality. The Hungarian method is one of the early triumphs of combinatorial optimisation, and it remains the standard exact solver for assignments at sizes where the linear program is small enough to fit in memory.

The assignment problem has a structural property that the next section turns into the basis of the general theory: although the objective is linear in the assignment, the feasible set is a finite, highly non-convex collection of permutation matrices. Replacing this combinatorial set with its convex hull, which is the set of doubly stochastic matrices by a theorem of Birkhoff, converts the problem into a linear program whose vertices are exactly the permutations. Optimisation over a polytope, with a linear objective, attains its extremum at a vertex; so the LP solution recovers an assignment, and the general theory recovers Monge's problem in a form for which an optimum is guaranteed to exist.

3A history of optimal transport: a literature review

Optimal transport is a subject whose modern shape was assembled over two centuries. Each major step removed a definite obstruction to the previous one, and the account that follows is deliberately causal: in each case the question is which earlier failure the next idea was built to answer.

Monge's 1781 memoir [1] stated the problem in geometric language and gave a partial geometric characterisation of the optimum for the cost given by Euclidean distance. The treatment was incomplete in two ways that mattered for the next century and a half. The existence of an optimal map was not guaranteed under Monge's formulation, because, as the first section noted, no map at all may push a given source onto a given target. And the characterisation, in terms of straight lines that do not cross, was qualitative rather than constructive. Methods for actually finding the optimum, beyond a handful of small examples, did not exist.

Leonid Kantorovich, writing in 1942 in the wartime Soviet Union [2], transformed the problem. He replaced Monge's deterministic map with a transport plan, a joint probability measure on the product of the source and target spaces, whose two marginals are the given μ and ν. A plan can split mass at a single source between multiple targets, which a map cannot, and the set of plans for given marginals is always non-empty. Kantorovich observed that the total cost is linear in the plan and the constraints on the plan are linear in its values, so the problem becomes a linear program. The existence of an optimum becomes a standard fact in linear programming, and the characterisation of the optimum becomes a question for the methods of linear duality. Kantorovich was led to this formulation through his earlier work on economic planning, and the dual potentials of the next section coincide with the prices in his economic interpretation. He later shared the 1975 Nobel Memorial Prize in Economic Sciences for this line of work.

The computational treatment of the linear program followed quickly from Dantzig's invention of the simplex method [3], and a specialised variant for the structure of the transportation problem, the network simplex, gave the first practical algorithm for large instances. By the mid-1950s, optimal transport was a workable problem of operations research. The mathematical question of when the Kantorovich optimum coincides with a Monge map, however, remained open in essential generality.

The decisive step on that question was taken by Yann Brenier in 1991 [4]. For the squared Euclidean distance cost and a source measure that is absolutely continuous with respect to Lebesgue measure, Brenier proved that the optimal plan is supported on the graph of a unique map, and that this map is the gradient of a convex function. This theorem connected optimal transport to convex analysis and to the Monge–Ampère equation, and it answered the existence question Monge had not been able to answer: under mild conditions the Monge problem does have a solution, and that solution is encoded in the same convex potential that the Kantorovich dual produces. Earlier independent work by Knott and Smith had given a similar result in different language; Brenier's framing made the connection to convex analysis explicit, and the resulting theory of monotone maps quickly became a standard tool in analysis.

The third major reframing came from Robert McCann in 1997 [5], who studied not the optimal map at a single moment but the interpolation between two distributions. He defined a geodesic in the space of probability measures whose endpoints are μ and ν, and whose intermediate values are the displacement interpolants, obtained by moving each mass fraction along its assigned line at constant speed. The space of probability measures, equipped with the Wasserstein-2 distance, was revealed by these geodesics to be a metric space in which differential geometry could be done. Felix Otto's 2001 paper [6] made that geometry explicit: he described a formal Riemannian metric on the space of measures whose geodesics are precisely McCann's, and whose gradient flows are partial differential equations of the kind that arise in physics.

Jordan, Kinderlehrer and Otto had already opened the connection between optimal transport and partial differential equations in 1998 [7] by showing that the Fokker–Planck equation, which describes the time evolution of a probability density under drift and diffusion, is the gradient flow of the free energy in the Wasserstein metric. Their JKO scheme made this constructive: it solves a sequence of finite-time minimisation problems on the metric space of measures, and in the limit of small time steps recovers the partial differential equation. The construction was introduced because no constructive method existed for gradient flows on the space of measures; once the metric was available, the scheme followed. Benamou and Brenier [8] developed a complementary fluid-dynamic interpretation in 2000, characterising the Wasserstein-2 distance as the minimum kinetic energy of any time evolution carrying μ to ν, and turning the static linear program into a partial-differential-equation-constrained problem amenable to numerical schemes.

The most recent wave of work has been computational and was driven by machine learning. Solving the exact linear program for transport between two empirical distributions of size ten thousand is already beyond the simplex method, and a million-point source against a million-point target was unreachable. Sinkhorn and Knopp had shown, long before any of this was a computational issue, that a positive matrix can be made doubly stochastic by alternately scaling its rows and columns [9]. Marco Cuturi observed in 2013 [10] that the entropically regularised optimal transport problem, in which a small multiple of the entropy of the plan is subtracted from the objective, has a unique optimum that is exactly a doubly scaled version of the Gibbs kernel exp(-C/ε), and that the Sinkhorn iteration converges to it. This made approximate optimal transport parallelisable and differentiable, two properties indispensable for modern machine learning, and within a few years the Sinkhorn method had become the standard way to compute large-scale optimal transport. Villani's two textbooks [11] consolidated the analytic theory; the computational synthesis is set out by Peyré and Cuturi [12].

Two threads run through this history. The first is a movement from deterministic maps toward probabilistic plans and from static assignments toward dynamic flows of measures: optimal transport began as a question about moving a pile of earth and became a question about the geometry of probability measures and the equations governing their evolution. The second is a movement from exact, expensive solvers toward approximate, scalable ones: the progression from the simplex through interior-point methods to entropic regularisation and Sinkhorn iteration was driven by the rising size of the problems that practitioners wanted to solve.

4The Kantorovich relaxation

To turn Monge's problem into one for which an optimum always exists, Kantorovich enlarged the set of admissible objects. Where Monge asked for a deterministic map sending each source point to a unique target point, Kantorovich admitted a transport plan: a joint probability measure on source and target whose two marginals are the given source and target distributions. The plan can be visualised as a mass distribution on the product space. The amount of mass located at a particular source–target pair is the amount of source mass at the first coordinate that is sent to the second.

Three properties of this relaxation matter. First, the feasible set is always non-empty: the product measure μ ⊗ ν is a plan with the right marginals, although in general not the optimal one. Second, the objective, which is the integral of the cost against the plan, is linear in the plan, and the marginal constraints are linear in its values. The Kantorovich problem is therefore a linear program in infinite or finite dimension according to whether μ and ν are continuous or discrete. Third, the relaxation is exact in the following precise sense: when the optimal Monge map exists, the Kantorovich optimum is supported on its graph, so passing to plans does not change the optimal cost. The relaxation is genuinely a relaxation, not a different problem.

The figure below illustrates a case in which the Monge problem has no solution but the Kantorovich problem does. Three source points with unequal weights must be transported to four target points, again with unequal weights. No one-to-one map can connect three locations to four. The optimal Kantorovich plan resolves this by splitting mass: a single source contributes to several targets, fractionally. The line widths in the figure encode the amount of mass moved along each pair, and the configuration of the plan changes continuously as the source point is moved.

The relaxation reveals a structural fact: although the feasible set is the set of all plans with the given marginals, the optimum is attained at an extreme point of that set, by general results in linear programming. For uniform marginals on n source and n target points the extreme points are precisely the permutation matrices, which is the Birkhoff–von Neumann theorem; the LP optimum is therefore an assignment. For non-uniform marginals the extreme points are the sparse plans with at most n + m − 1 non-zero entries. The Kantorovich relaxation is convex, but its optima are as sparse as the constraints allow.

The Birkhoff polytope, named for the 1946 theorem of Garrett Birkhoff that the set of n by n doubly stochastic matrices is the convex hull of the n by n permutation matrices, is the geometric object underlying this fact in the uniform-marginal case. A doubly stochastic matrix is a non-negative matrix whose rows and columns each sum to one, which is precisely the discrete Kantorovich constraint when the source and target are both uniform on n points. The Birkhoff polytope therefore has n! vertices, exactly the assignments, and the Kantorovich linear program is the minimisation of a linear functional over this polytope. The fundamental theorem of linear programming says that the optimum is attained at a vertex, recovering the discrete assignment from the convex relaxation. The same picture holds, with the obvious changes, for non-uniform marginals, where the polytope is more general and its vertices are the sparse plans.

5The cost matrix, the coupling and the dual problem

A finite Kantorovich problem is determined by three objects: the source weights, the target weights, and the cost matrix C, whose entry C(i, j) gives the cost of moving a unit of mass from the i-th source to the j-th target. The choice of cost is the modeller's, and almost every property of the problem depends on it. Squared Euclidean distance is the canonical choice for two reasons. It is smooth and convex, which makes the analysis tractable, and it makes the Wasserstein-2 distance defined in the eighth section a true metric on the space of probability measures with finite second moment. Other costs are useful in other contexts: the Euclidean-distance cost gives the Wasserstein-1 distance and a dual theory involving Lipschitz functions; costs that are not powers of distance arise in matching markets and in image processing.

The figure above shows the cost matrix and the optimal coupling matrix side by side as heat grids, computed for three source and three target points whose positions can be moved. The two grids encode different objects, but they are tightly linked. The coupling is concentrated on the cells where the cost is small, subject to the marginal constraints. In the limit of vanishing regularisation, the coupling is supported on a set of at most five cells out of nine; with the small regularisation the figure uses, the support is slightly diffused but the dominant pattern is the same.

Underlying this concentration is a duality theorem due to Kantorovich. Each linear program has a dual, and for the transportation problem the dual is a maximisation over pairs of functions f on the source and g on the target. The objective is the sum of the integrals of f against μ and g against ν, and the constraint is the inequality f(x) + g(y) ≤ c(x, y) for every source point x and every target point y. Strong duality holds, so the primal minimum equals the dual maximum, and the optimal coupling is supported on the set where the dual constraint is tight: f(x) + g(y) equals exactly c(x, y). This is the complementary slackness of the linear program restated in the language of transport.

The dual potentials f and g carry an economic interpretation that Kantorovich had in mind from the start. Think of f(x) as a price a producer charges for releasing a unit of mass at the source point x, and g(y) as a price a consumer pays for receiving a unit at the target point y. The constraint says that no one would ever buy at x and sell at y if doing so were cheaper than transporting between them. The market clears when the constraint is tight on the locations actually traded between, and the resulting prices are the dual potentials. For the special case c(x, y) = |x − y|, the dual can be written in terms of a single Lipschitz function f with g = −f, giving the Kantorovich–Rubinstein duality used widely in probability and statistics, and now central to certain generative modelling architectures.

A geometric condition characterises the support of an optimal plan in any dimension and goes by the name of c-cyclic monotonicity. For the squared Euclidean cost in two dimensions it reduces to the following: if (x₁, y₁) and (x₂, y₂) are two pairs in the support of the optimal plan, then the sum |x₁ − y₁|² + |x₂ − y₂|² is at most |x₁ − y₂|² + |x₂ − y₁|², which is to say that swapping the two assignments cannot lower the cost. This is the continuous analogue of the local-optimality condition the two-opt search of the first figure tested, and it is the multidimensional generalisation of the one-dimensional monotonicity that the seventh section will use to derive the closed-form solution. Rüschendorf showed in the 1980s that c-cyclic monotonicity, in fact, characterises optimality globally as well as locally for a wide class of costs, and that result is one of the structural facts behind Brenier's later proof that the optimal map for quadratic cost is the gradient of a convex function.

6Sinkhorn and entropic regularisation

The linear program of the Kantorovich problem has nm variables in the discrete case of n source and m target points, and on contemporary machine-learning data those values can each be in the millions. The standard simplex and interior-point solvers scale poorly in that regime, both in time and in their inability to use parallel hardware. The remedy that has come to dominate modern optimal-transport computation is entropic regularisation, introduced to this problem in machine learning by Cuturi [10].

The regularised problem subtracts from the Kantorovich objective a small multiple of the entropy of the plan, written H(π) = −Σ π log π. The regularised problem is strictly convex, has a unique optimum, and its solution takes a closed form. The optimal regularised plan factors as π(i, j) = u(i) K(i, j) v(j), where K(i, j) = exp(−C(i, j) / ε) is the Gibbs kernel and u and v are positive vectors fixed by the marginal constraints. The two constraints are linear in u and v separately, and the standard way to enforce them is to alternate: rescale the rows of the current π to match the source marginal, then rescale the columns to match the target marginal, and repeat. This is the Sinkhorn iteration, named for Richard Sinkhorn, who showed in 1967 with Knopp [9] that the procedure converges from any positive starting point.

The figure above runs the iteration step by step. The Step button performs one row scaling and one column scaling; the Run button continues until the marginals match to numerical precision; the Epsilon slider sets the regularisation. Two qualitative facts can be read off the figure. The matrix evolves from the initial Gibbs kernel, with mass spread roughly evenly over all source–target pairs, toward a concentrated plan whose entries are largest where the cost is smallest. The smaller the regularisation, the sharper the final plan and the closer it is to the unregularised Kantorovich optimum; but smaller regularisation also slows the convergence, because the kernel becomes more singular and the alternating scaling has more work to do.

Two practical remarks complete the picture. First, the regularised plan is not exactly the unregularised optimum; the regularisation biases the cost slightly upward and biases the entropy of the plan slightly upward as well. For machine-learning applications this bias is usually acceptable, and where it is not, the Sinkhorn divergence, a debiased combination of three Sinkhorn distances, removes it. Second, at very small regularisation the kernel exp(−C/ε) underflows to zero in double precision, and the multiplicative iteration becomes unstable. Production implementations carry the iteration in the log domain, working with log u and log v and the log-sum-exp trick, and the resulting scheme is numerically robust at regularisations many orders of magnitude smaller than the naive implementation can support.

A second refinement removes the bias the regularisation introduces. The regularised optimum is not the unregularised optimum; it pays a small premium for spreading mass more diffusely than the unregularised plan would. The Sinkhorn divergence, defined by Genevay and colleagues as the regularised Sinkhorn distance minus half of the Sinkhorn distance from each measure to itself, vanishes when the two measures agree and recovers the unregularised optimal transport in the limit of vanishing regularisation. Used together with the log-domain Sinkhorn iteration, the Sinkhorn divergence is the form in which entropic optimal transport now most commonly appears in production code, and it is the divergence used to fit generative models and to perform domain adaptation when the empirical samples number in the hundreds of thousands.

7The one-dimensional case and the monotone map

Optimal transport in one dimension is exceptional. For any cost of the form c(x, y) = h(y − x) with h strictly convex, the optimal Monge map exists, is unique, and has a closed-form expression: it is the unique non-decreasing map sending the source to the target. The monotonicity is what makes the problem trivial in one dimension and delicate in higher dimensions, where monotonicity in one direction says nothing about behaviour transverse to that direction.

The non-decreasing map between two one-dimensional measures has a familiar name. If F_μ is the cumulative distribution function of μ and F_ν is the cumulative distribution function of ν, the optimal Monge map T satisfies T(x) = F_ν⁻¹(F_μ(x)), the composition of the source CDF with the target quantile function. To send mass at the source point x to a target point, look up its source quantile, then send it to the matching target quantile. This is the quantile transport map, and it is the unique non-decreasing function pushing μ onto ν.

The figure above shows two Gaussian densities on a number line and draws the quantile map as a family of arrows from points of the source to their images in the target. The probe arrow, drawn in full, can be moved along the source axis, and its target is computed by inverting the target's cumulative distribution numerically. The figure also displays the closed-form Wasserstein-2 distance between the two Gaussians, which for one-dimensional Gaussians equals the square root of (μ₁ − μ₂)² + (σ₁ − σ₂)². The closed form is a rare luxury and a useful diagnostic: a numerical Sinkhorn or LP solver applied to one-dimensional empirical data should reproduce it.

The one-dimensional case is more than a special example. Two facts extend its reach. First, optimal transport on the real line, for the quadratic cost, has the same algebra as optimal transport between the marginals of higher-dimensional measures along a one-dimensional projection; that observation gives the sliced Wasserstein distance, which averages the one-dimensional Wasserstein distances over random projection directions. Second, the monotone-rearrangement argument that establishes the one-dimensional optimal map generalises, in higher dimensions, to the c-monotonicity that characterises optimal plans, the multivariable analogue of monotonicity that Rüschendorf used to develop the theory in the eighties.

The sliced Wasserstein distance deserves a closer look, because it is the simplest way to scale optimal transport to large or high-dimensional data. Given two high-dimensional empirical measures, one chooses a random unit direction, projects each measure onto that direction to obtain two one-dimensional empirical measures, computes their exact one-dimensional Wasserstein distance by sorting and pairing the projected samples, and averages the result over many random directions. The complexity per direction is dominated by a sort, which is linear-logarithmic, and the procedure parallelises trivially across directions. The sliced distance is a true metric on the space of measures, and it agrees with the one-dimensional Wasserstein when the data is one-dimensional. It does not, in general, equal the full Wasserstein distance in higher dimensions, but it captures enough of the same geometry to serve as a practical surrogate in applications such as generative modelling and shape comparison, and it is now a standard part of the computational toolkit.

8The Wasserstein distance

The minimum cost of the Kantorovich problem with c(x, y) = |x − y|^p is a function of the source and target measures alone; raising it to the power 1/p yields the Wasserstein-p distance, written W_p(μ, ν). For p = 1 and p = 2 these are the two cases of greatest practical use. The Wasserstein-2 distance is the one associated with squared cost and with the geometry of probability measures discussed in the next section.

The Wasserstein-p distance is a genuine metric on the space of probability measures with finite p-th moment. It is non-negative; it equals zero only when μ equals ν; it is symmetric; and it satisfies the triangle inequality. These four properties separate the Wasserstein distance from the more familiar Kullback–Leibler divergence, which is non-symmetric, fails the triangle inequality, and assigns infinity to any pair of measures with disjoint support. That last failure of the Kullback–Leibler divergence is the one that most often pushes practitioners toward Wasserstein. Two distributions slightly displaced from each other, but with the same shape, are near in Wasserstein distance and at infinite KL distance, and the Wasserstein answer is the one that matches geometric intuition.

The figure above shows two clouds of points in the plane and the Wasserstein-2 distance between them, computed by Sinkhorn applied to the empirical measures. Moving the source cloud changes the distance continuously, and the optimal coupling is drawn faintly as line segments whose widths encode the moved mass. Two structural properties of the Wasserstein distance are visible. Translating the source cloud rigidly changes W₂ by an amount equal to the translation length, because the optimal coupling is also translated rigidly. Widening the source cloud, holding its centre fixed, changes W₂ by an amount that depends on the change in the standard deviation, and for isotropic Gaussian clouds the change is given by the closed form from the previous section.

The Wasserstein distance also has a dynamic interpretation that rounds out its character. By a result of Benamou and Brenier [8], W₂(μ, ν) equals the square root of the minimum kinetic energy of any time-dependent flow of mass that carries μ to ν over a unit time interval. The optimal flow is the one whose velocity field is the gradient of a convex function, and it traces out the geodesic of the Wasserstein space discussed next. Optimal transport is therefore not only a static assignment but a least-energy flow, and the same distance describes both views.

The Benamou–Brenier dynamic reformulation recasts the static linear program as a partial-differential-equation-constrained variational problem. The unknowns are a curve of densities ρ(t, x) joining μ at time zero to ν at time one, together with a velocity field v(t, x), related by the continuity equation that expresses conservation of mass. The objective is the time integral of the kinetic energy carried by the mass, and the constraint is the conservation law alone; the marginals at the two endpoints are the boundary conditions. The minimum of this dynamic problem is exactly the squared Wasserstein-2 distance, and the minimisers are the geodesics in the Wasserstein space. The reformulation is appealing because it turns a problem on the polytope of plans into one on the much smaller space of densities and velocity fields, and it is the starting point for several numerical methods that scale better than the static linear program in high dimensions.

9Gradient flows in Wasserstein space

The Wasserstein-2 distance turns the space of probability measures with finite second moment into a metric space, and McCann's displacement interpolation [5] equips it with geodesics. With a distance and geodesics the space becomes formally Riemannian, and the formal Riemannian structure worked out by Otto [6] provides the machinery of differential geometry on measures. A functional assigning a real number to each measure has, on this Riemannian space, a gradient, and a measure can be evolved in the direction opposite to that gradient. The resulting curve of measures is the gradient flow.

The most consequential gradient flow is the one Jordan, Kinderlehrer and Otto identified [7]. They showed that the Fokker–Planck equation, ∂_tρ = ∇·(ρ∇V) + Δρ, which describes a density evolving under an external potential V together with diffusion, is the gradient flow of the free energy F(ρ) = ∫V ρ + ∫ρ log ρ under the Wasserstein-2 metric. The drift term is the gradient of the first integral, the diffusion term is the gradient of the entropy, and the partial differential equation is the assembly of the two pieces. The same identification turns the heat equation into the gradient flow of entropy alone, the porous-medium equation into the gradient flow of a generalised entropy, and other equations of mathematical physics into gradient flows of correspondingly natural functionals.

The figure above runs a discretised gradient flow toward a target cloud. At each step the algorithm computes the optimal coupling between the current source cloud and the target, then moves each source particle a small fraction of the way toward the barycentre of its target assignments. This is the discrete-time analogue of a velocity field given by the gradient of the Wasserstein potential, and over enough steps the source cloud matches the target. The figure honours prefers-reduced-motion by jumping to the final configuration; otherwise it animates the evolution.

The reason this matters beyond aesthetic appeal is that gradient flows give a constructive route to the solutions of the partial differential equations they correspond to. The JKO scheme, implementing the gradient flow as a sequence of minimisations on the space of measures, converges in the small-time-step limit to the underlying equation. For non-linear and degenerate diffusions where classical methods are awkward, the JKO scheme provides both a proof of existence and a numerical method. The optimal-transport point of view, in this sense, has reorganised a substantial part of the analysis of partial differential equations of evolution.

10Conclusion

The narrative arc of optimal transport is the conversion of a physical question of moving earth into a structure on the space of probability measures. The question began with Monge as a problem of maps and partial differential equations; it became, with Kantorovich, a linear program with a dual structure that admits an economic interpretation; it returned, with Brenier, to the language of maps and convex potentials in a form that finally resolved Monge's existence question; it left, with McCann and Otto, the static realm altogether for the geometry of measures and the evolution of densities; and it returned, with Sinkhorn iteration and its later refinements, to the computational world where million-point empirical measures can be transported into one another in seconds on ordinary hardware.

A practical postscript on tuning is worth adding here, because the theory of gradient flows and the practice of computational optimal transport meet on the question of how to take steps. The JKO scheme and its discrete analogues require a time step τ; too large a τ turns the discrete trajectory into a sequence of jumps that bear little relation to the continuous flow, while too small a τ wastes computation by taking infinitesimal steps. The standard practical recipe is an adaptive step size, in which τ is increased when the energy decreases as predicted and decreased when the iterates begin to oscillate. Behind this recipe is the abstract observation that the Wasserstein gradient flow is a particular case of a metric gradient flow on a Polish space, for which a well-developed convergence theory due to Ambrosio, Gigli and Savaré is available.

The applications follow from this breadth. In economics, the dual potentials are equilibrium prices in matching markets, and optimal transport gives the existence and characterisation of those prices. In probability, the Wasserstein distance is the natural metric on empirical distributions, and the laws of large numbers and central limit theorems for the Wasserstein distance are an active area of study. In partial differential equations, the JKO scheme provides a variational characterisation of a large class of evolution equations. In machine learning, the Wasserstein distance and its regularised approximations have given rise to generative models, to domain-adaptation algorithms, and to a wide family of methods for comparing histograms and point clouds. In each of these settings the same five objects appear: a source measure, a target measure, a cost function, a plan or map, and a distance. The story this monograph has told is how each of those objects came to be defined.

The reader who has followed the development is positioned to read further into any of the application areas, because each rests on the same five-piece structure. The methods of Villani [11] develop the analytic side at length; the computational synthesis of Peyré and Cuturi [12] is the reference for the algorithms. The seven figures in this document are not a substitute for either of those treatments, but they make concrete, in the case of small two-dimensional examples, what those treatments handle in generality. The geometry on display in the figures is, in every case, the same geometry that governs the high-dimensional problems.