Jean-David Benamou
email: jean-david.benamou@inria.fr

Optimal transportation, modelling and numerical simulation

Keywords: ulation, volume=30, DOI=10.1017/S0962492921000040, journal=Acta Numerica, publisher=Cambridge University Press, author=Benamou, Jean-

We present an overview of the basic theory, modern optimal transportation extensions and recent algorithmic advances. Selected modelling and numerical applications illustrate the impact of optimal transportation in numerical analysis.

applications can be found in the monographs by [START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Villani | Optimal Transport: Old and New[END_REF] and -in a more probabilistic setting - [START_REF] Rachev | Mass Transportation Problems: Applications, Probability and its Applications[END_REF]. They were mostly theory-oriented.

Offering a Riemannian metric on the space of measures defined on general manifolds, optimal transportation has spread far in mathematics and scientific computation. The construction of Wasserstein gradient flows, for instance, introduced by [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], led to a new understanding of nonlinear parabolic equations and is an active research field. Developed on the infinite-dimensional space of measures but also including discrete, possibly finite 'empirical' measures, optimal transportation theory offers an elegant bypass of the discretization process, a familiar but sometime burdensome routine for numerical analysts working on continuous models.

Use of the optimal transportation framework in the natural sciences appeared during the same period. The fundamental paper of [START_REF] Brenier | Polar factorization monotone rearrangement of vector-valued functions[END_REF] was in fact preceded and inspired by an earlier paper, [START_REF] Brenier | The least action principle and the related concept of generalized flows for incompressible perfect fluids[END_REF], constructing weak solutions to Euler geodesics variational interpretation [START_REF] Arnold | Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF]. Other applications began to appear in fields as diverse as medical imaging [START_REF] Haker | Optimal mass transport for registration and warping[END_REF], meteorology [START_REF] Cullen | An extended Lagrangian theory of semigeostrophic frontogenesis[END_REF], astrophysics [START_REF] Frisch | A reconstruction of the initial conditions of the universe by optimal mass transportation[END_REF] and optics [START_REF] Glimm | Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem[END_REF]Oliker 2003, Wang 2004).

Similar developments occurred in parallel in economics, for which optimal transportation tools have been instrumental, particularly in the study of the famous 'principal-agent' problem; see [START_REF] Carlier | A general existence result for the principal-agent problem with adverse selection[END_REF], [START_REF] Ekeland | Notes on optimal transportation[END_REF], [START_REF] Figalli | When is multidimensional screening a convex program[END_REF] and [START_REF] Rochet | sweeping, and multidimensional screening[END_REF]. The role of optimal transportation in economics is presented in depth in [START_REF] Galichon | Optimal Transport Methods in Economics[END_REF]; it will not be covered in this paper. Other important variants of optimal transportation will be omitted, for example transport occurring on discrete Markov chains (Maas 2011) and branched transport and optimal transportation networks; see [START_REF] Pegon | Transport branché et structures fractales[END_REF] and [START_REF] Bernot | Optimal Transportation Networks: Models and Theory[END_REF] for a review of these topics. It continues to pop up in fields as diverse as number theory [START_REF] Steinerberger | On a Kantorovich-Rubinstein inequality[END_REF], quantum mechanics [START_REF] Golse | Quantum and semiquantum pseudometrics and applications[END_REF], relativity [START_REF] Cavalletti | Optimal transport in Lorentzian synthetic spaces, synthetic timelike Ricci curvature lower bounds and applications[END_REF], genomics [START_REF] Lavenant | Unconditional convergence for discretizations of dynamical optimal transport[END_REF] . . . The generic optimal transportation problem is a nonlinear optimization problem. Kantorovich relaxation (which earned Kantorovich the Nobel Prize for Economics) as a linear program was a giant leap forward for the theory but at the heavy computational price of replacing the space of transport maps with the space of couplings, hence squaring the number of unknowns. Combined with the cubic complexity of linear programming solvers, it was never used as such in the abovementioned applications.

The tighter convex relaxation proposed by [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] offered the first convergent and tractable algorithm. It is difficult to assess its computational complexity precisely because it adds a time dimension, but one can safely estimate that it is somewhere between quadratic and cubic. The proposed dynamic vision of optimal transportation has also been inspirational, and has strongly influenced subsequent research in the field.

Significant progress was achieved in the early 2010s. Based on the convexity of the 'semi-dual' formulation of optimal transportation, a damped Newton method was proposed to solve the Monge-Ampère [START_REF] Benamou | Numerical solution of the optimal transportation problem using the Monge-Ampère equation[END_REF] and semi-discrete [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF]) settings, respectively. This approach resulted in the first linear complexity solvers. The 'damping' ensures that the optimization is kept to the domain of strict convexity of the optimal transportation cost. So far this approach has essentially been limited to the most commonly used Euclidean quadratic cost (1.1), but a general analysis of the method is available in [START_REF] Kitagawa | Convergence of a Newton algorithm for semi-discrete optimal transport[END_REF].

Another strategy, based on a convexification of Kantorovich relaxation, was initiated by [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] and in the economics context by [START_REF] Salanié | Cupid's invisible hand: Social surplus and identification in matching models[END_REF] independently. It has been popular ever since for several different reasons. It does not solve the optimal transportation problem but yields its 'Schrödinger entropic approximation', a connection well known in probability and statistical physics [START_REF] Léonard | A survey of the Schrödinger problem and some of its connections with optimal transport[END_REF]. It is simple to implement and applies to virtually any optimal transportation generalization. There is a trade-off between the computational cost and the 'entropic bias', or the precision error, with classical optimal transportation. It is very competitive for estimating the transportation cost but not the transport plan itself. Finally it is easily parallelizable and, while naively quadratic, linear-cost off-the-shelf software can be found in [START_REF] Feydy | Geometric loss functions between sampled measures, images and volumes[END_REF], for instance.

These numerical breakthroughs triggered an acceleration of the use of optimal transportation for numerical modelling and solutions in an increasing number of domains, including the recent and productive machine learning topic. Following all these developments is far beyond my mental abilities and energy. Here is a first disclaimer: this review paper will necessarily be incomplete, and biased towards my personal scientific environment. I apologize for the important contributions I have missed.

The first lecture notes oriented towards PDEs I came across -those of [START_REF] Evans | Partial differential equations and Monge-Kantorovich mass transfer[END_REF] -are still available online and useful! As optimal transportation is now taught at graduate level across the globe, a wealth of such lecture notes is available on the internet. Several books are also out: the book by [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians: Calculus of Variations[END_REF] is a modern self-contained review of the theory with a wide range of applications. As suggested by the title, Optimal Transport for Applied Mathematicians, this is definitely the reference book for applied mathematicians working in this field. In § 8.4.4, for instance, it provides an easy introduction to semi-geostrophic equations [START_REF] Cullen | A Mathematical Theory of Large-Scale Atmosphere/Ocean Flow[END_REF] and early universe reconstruction [START_REF] Frisch | A reconstruction of the initial conditions of the universe by optimal mass transportation[END_REF], two important applications of optimal transportation. The paper by [START_REF] Peyré | Computational optimal transport[END_REF] is a precious resource with precise descriptions of the algorithms and their analysis, in particular entropic regularization. It is nicely complemented by the book chapter by [START_REF] Mérigot | Optimal transport: Discretization and algorithms[END_REF], which presents the semi-discrete approach in more detail along with its link to entropic regularization. Finally, [START_REF] Carlier | Classical and Modern Optimization[END_REF] has provided a didactic presentation of the mathematics and the algorithmics of convex duality and convex optimization, which are ubiquitous in optimal transportation. As we will try to emphasize, convex Fenchel-Rockafellar duality is a cornerstone of optimal transportation. A good knowledge of these techniques, and the associated methods of solution, is compulsory for numerical analysts interested in optimal transportation as a research topic.

This paper has been designed and organized to present, in a unified and compact framework, the basic theory, the motivations and also the insight behind recent optimal transportation developments, made possible by the algorithmic advances mentioned above. It synthesizes a large corpus; each section could be expanded into a longer technical survey paper. Except for the last section, the paper should be browsed or read linearly, at least the first time. Sections present -in a non-technical way -elements of the theory and algorithms together with flagship applications, and often use earlier expositions. The presentation will remain formal with some elements of proof. The chain of logic and the derivation of formulae will often be sketchy, and may be considered as exercises for the reader. I try to indicate precisely references with detailed proofs and rigorous mathematical formulations. In an attempt to keep the presentation self-contained, I have also gathered some required basic tools and results in footnotes. A list of notations is available after this introduction.

By definition, this paper should be useless for optimal transportation experts. It may, however, offer a broader vision of the topic and unseen connections between the many concepts used. For beginners, it cannot replace careful study of the monographs cited above but may instead be used as a study guide; each section roughly corresponds to one or two lectures.

A second disclaimer: my knowledge and understanding is of course the result of many years, even decades, of listening to lectures, discussing with colleagues and reading the literature. The scientific fruits will appear in some form or another in this paper.

Many thanks to Irène Wadspurger, Francis Collino, Guillaume Carlier and Gabriel Peyré for their pre-reading and correction of this manuscript.

I am deeply grateful to the Editorial Board of Acta Numerica for their invitation. Many thanks to Glennis Starling for copy-editing the document and improving my English.

This foreword would not be complete without acknowledging the fantastic group of colleagues and students at MOKAPLAN and our close collaborators. The list is too long to be reproduced here, but if they read these lines they will know they are part of the family. I thank them warmly for their friendship, patience, contributions, advice and vision.

Outline

Section 1 contains in a condensed form the classical presentation of the original Monge optimal transportation, the Kantorovich linear program and basic linear programming and linear assignment methods. This can be found in the general references above. Notations are introduced, and in particular the important idea that we can use the same set of notations for measures with continuous densities and discrete measures.

See Santambrogio (2015, §1), Carlier (2021, §6) and Peyré and Cuturi (2019, §2) for comprehensive presentations.

Section 2 introduces the Kantorovich dual problem and its role in the study of the well-posedness of optimal transportation. The key notion is the 𝑐-transform optimality of the Kantorovich potentials. This can also be used to tighten the linear Kantorovich problem into a concave 'semi-dual' problem. It leads to two modern numerical approaches: the Monge-Ampère and semi-discrete approaches. We present an application to the reflector cost.

As for Section 1, see Santambrogio (2015, §1), Carlier (2021, §6) and Peyré and Cuturi (2019, §2), and also Mérigot and Thibert (2020, §2).

Section 3 gives the basics of the distance on the space of probability measures. It also introduces entropic optimal transportation and the Sinkhorn algorithm and discusses the smoothing bias on the distance introduced by the approximation. This bias may be useful when computing interpolations or barycentres of discrete measures for the optimal transportation metric. These notions naturally lead to the introduction of time and dynamic optimal transportation in the next section.

Rigorous proofs (there are different techniques) are detailed in Santambrogio (2015, §5). See Peyré and Cuturi (2019, §4) for entropic optimal transportation. Section 4 brings the notion of time or dynamics to optimal transportation. This led to the computational fluid dynamics formulation (also known as 'Benamou-Brenier'). This approach seems artificial when considering the classical source-totarget optimal transportation problem but was beneficial from the computational point of view. It also helps us to understand the richer models in Sections 5 and 6.1, and the Riemannian-like distance over the set of probability densities used, in particular, in the Wasserstein gradient flow theory (Section 7.2).

Classical first-order optimization 'proximal splitting' methods are applicable. See Carlier (2021, §7) and the review by [START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF].

Section 5 covers the optimal transportation treatment of Euler geodesics, its connection with dynamic optimal transportation and the use of optimal transportation solvers in this context.

I recommend [START_REF] Daneri | Variational Models for the Incompressible Euler Equations[END_REF] as a complete and precise review. See also Brenier (2020, §2).

Section 6 makes the connection between entropic optimal transportation and the Schrödinger problem in statistical physics. The entropic bias can also be interpreted as adding a diffusion to the transport model. The application of the Sinkhorn algorithm to variational mean field games is a good illustration. It also changes the nature of the transport from deterministic to stochastic and hyperbolic to parabolic. Quite remarkably, using the volatility as a control allows us to extend most of the classical optimal transportation results. A good example is martingale optimal transportation used in finance. [START_REF] Léonard | A survey of the Schrödinger problem and some of its connections with optimal transport[END_REF] is THE reference. It is a probabilistic presentation but the (long) introduction is accessible to all.

Section 7 introduces arguably the most active area of applications. We move one level up and use the optimal transportation distance as a metric (or 'loss' in the machine learning terminology) for variational problems set on the space of probability and non-negative Radon measures. It was a crucial motivation in the development of 'unbalanced' optimal transportation distances and the introduction of a simple but very efficient variation of entropic optimal transportation called Sinkhorn divergence. This is illustrated via gradient flow problems.

See Santambrogio (2015, §8) for Wasserstein gradient flows and Peyré and Cuturi (2019, §7-8) for losses and applications. The reference monograph on the theory of gradient flows is the book by [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF].

Section 8 covers a few subjects and directions, which I have been interested in and are open or partly open.

Introduction

What is optimal transportation?

A general presentation on optimal transport traditionally starts with a reference to Gaspard Monge's problem 'des déblais et des remblais' (1781): minimize the amount of work to shovel some mass from a source to a target 'configuration' (we will define this last term more precisely). The energy is defined as the total distance the shovel, which contains exactly one unit of mass, has been carried over. Gaspard Monge was a military engineer and his motivation concerned fortifications: shovelling a ditch and using the soil to elevate a protective bank. In real life, the actual digging/building process is affected by accessibility/structural stability constraints: shovel transfers are done in sequence, and some sequences are infeasible (imagine you want to start placing the dirt on top of the earth bank with no foundations). One way to keep the Monge fortification analogy in line with modern mathematical formulations of optimal transportation is to use a 'magic' shovel. The magic shovel would be able to dig and carry all units of mass (one shovel) instantly, according to a predefined mapping, from the source to the target configuration.

We therefore look for a 'transport' map 𝑇, in the mathematical sense, between two source/target sets, denoted 𝑋 0 and 𝑋 1 , i.e. 𝑇 : 𝑋 0 → 𝑋 1 . The simplest situation, considered in most of this review, corresponds to 𝑋 0 = 𝑋 1 = 𝑋, a compact metric space. The source and target 'configurations' are real-valued non-negative mass measures on 𝑋 0 and 𝑋 1 ; we will first assume they are non-negative Radon measures, denoted 𝜇 0 and 𝜇 1 . In mathematical notation, 𝜇 0 ∈ M + (𝑋 0 ) and 𝜇 1 ∈ M + (𝑋 1 ). Because the shovel always carries the same amount of mass, 1 Continuous linear forms over the set of continuous functions on the compact space 𝑋.

the transport map is constrained to 'conserve the mass'. This translates into the following mathematical measure-preserving constraint.

For any 𝜇 1 -measurable set 𝐴 1 ⊂ 𝑋 1 , 𝑇 -1 (𝐴 1 ) is a 𝜇 0 -measurable subset of 𝑋 0 and 𝜇 0 (𝑇 -1 (𝐴 1 )) = 𝜇 1 (𝐴 1 ) (in the modern literature this is generally denoted 𝑇#𝜇 0 = 𝜇 1 : we say 𝑇 'pushes forward' 𝜇 0 to 𝜇 1 ). As a direct consequence, the total mass also needs to be preserved: 𝜇 0 (𝑋 0 ) = 𝜇 1 (𝑋 1 ). Setting the total mass to 1, we may restrict the space of 𝜇 0 and 𝜇 1 to probability measures. We will use the notation 𝜇 0 ∈ P(𝑋 0 ) and 𝜇 1 ∈ P(𝑋 1 ). In fact, 𝜇 0,1 will often be discrete empirical measures (sums of Dirac masses) but might also be absolutely continuous measures with respect to Lebesgue measure. In the continuous case and by abuse of notation the probability densities will still be denoted 𝜇 0,1 .

The Monge shovel 'carries a unit of mass' instead of volume because the relevant measure of work is mass times distance. In the simplest Monge configuration, densities are constant and can be scaled to 1, i.e. mass and volume are equivalent. In this case the mass conservation then implies that 𝑋 0 and 𝑋 1 have the same cardinality (a finite number of elements) or volume (infinite number of elements) for the set of admissible maps 𝑇 to be non-empty. If we want to tackle non-constant densities, the shovel has to pile several shovel/mass units on the same target location in order to reach the prescribed mass density and dig several units of mass from the same place, then send them to different locations. The second requirement rules out the map representation as each element in 𝑋 0 may only possess one image 𝑇(𝑥 0 ) ∈ 𝑋 1 .

Kantorovich relaxation

The mathematical framework to deal with this limitation was invented by Leonid Kantorovich in 1954. Instead of a map he proposed looking for couplings. These represent the amount of mass transported between all pairs (𝑥 0 , 𝑥 1 ) ∈ 𝑋 0 ×𝑋 1 . These couplings are called Kantorovich transport plans. They are probability measures, here denoted 𝜋, on the product space 𝑋 0 × 𝑋 1 : 𝜋 ∈ P(𝑋 0 × 𝑋 1 ). For all measurable subsets (𝐴 0 , 𝐴 1 ) ⊂ 𝑋 0 × 𝑋 1 , 𝜋(𝐴 0 , 𝐴 1 ) is the amount of mass transported from 𝐴 0 to 𝐴 1 . The mass conservation is now automatically enforced: you receive exactly what you send. The original formulation of the problem requires the transported mass to reach the exact target configuration. Therefore you cannot send more or less than what is available, and likewise for the reception. This translates into 2 It can also be expressed as

∫ 𝑋 0 ( 𝑓 • 𝑇) 𝜇 0 d𝑥 0 = ∫ 𝑋 1
𝑓 𝜇 1 d𝑥 1 for any continuous test function 𝑓 ∈ C(𝑋). If 𝑇 is smooth and one-to-one, the change of variable 𝑥 1 = 𝑇(𝑥 0 ) gives

∫ 𝑋 0 ( 𝑓 • 𝑇) d𝜇 0 d𝑥 0 = ∫ 𝑋 0 ( 𝑓 • 𝑇) det(𝐷𝑇) (𝜇 1 • 𝑇) d𝑥 0 ,
a variational formulation for the Jacobian equation det(𝐷𝑇) (𝜇 1 • 𝑇) = 𝜇 0 . 3 The (•) 0,1 will systematically denote the ordered pair ((•) 0 , (•) 1 ). 4 That is, mass per volume or how much mass is located at each element of the sets. marginal constraints on the plans: 1 𝑋 1 , 𝜋 𝑋 1 = 𝜇 0 and 1 𝑋 0 , 𝜋 𝑋 0 = 𝜇 1 (see the footnote for the definition of the brackets). This is often written using the canonical projection 𝑃 𝑋 0 : (𝑥 0 , 𝑥 1 ) → 𝑥 0 (not to be mistaken for P(𝑋), the probability space over 𝑋) and the push-forward notation 𝑃 𝑋 0 #𝜋 = 𝜇 0 and 𝑃 𝑋 1 #𝜋 = 𝜇 1 . The last ingredient we need to describe our first formal optimal transportation problem is the displacement cost 𝑐 : (𝑥 0 , 𝑥 1 ) ∈ 𝑋 × 𝑋 → 𝑐(𝑥 0 , 𝑥 1 ), i.e. the cost of moving the shovel from 𝑥 0 to 𝑥 1 . The mathematical and modelling properties of 𝑐 are of course of tantamount importance and we will return later to this point. In line with Monge modelling, this is often linked to the metric associated with 𝑋. The standard case is the square of the Euclidean distance and we will mostly restrict to this case in our presentation:

𝑋 0,1 = compact subsets of R 𝑑 , 𝑐(𝑥 0 , 𝑥 1 ) = 1 2 𝑥 1 -𝑥 0 2 . (1.1)
Kantorovich's problem will be the first formal optimal transportation problem discussed in this paper. Optimizers will always carry a (•) * throughout the paper, that is,

𝜋 * ∈ arg inf 𝜋 ∈Π(𝜇 0 , 𝜇 1 ) 𝑐, 𝜋 𝑋 0 ×𝑋 1 , (1.2)
where Π(𝜇 0 , 𝜇 1 ) is the following set of admissible transport plans:

Π(𝜇 0 , 𝜇 1 ) := {𝜋 ∈ P(𝑋 0 × 𝑋 1 ) : 𝑃 𝑋 1 #𝜋 = 𝜇 0 and 𝑃 𝑋 0 #𝜋 = 𝜇 1 }. (1.3)
We stress that (1.2)-(1.3) is a linear optimization problem with linear constraints in 𝜋. Minimizers exist but may not be unique simply on the condition that 𝑐 ∈ C(𝑋 0 × 𝑋 1 ) and Π(𝜇 0 , 𝜇 1 ) is compact for weak-★ topology of measures.

Let us give a trivial example. Choose 𝜇 0 = 𝛿 𝑎 0 and 𝜇 1 = 𝛿 𝑎 1 , two Dirac masses at 𝑎 0 and 𝑎 1 . We can set the unit of mass to 1 and Monge transport is done in just one shovel displacement between 𝑎 0 and 𝑎 1 . The transport cost is the displacement cost, 𝑐(𝑎 0 , 𝑎 1 ). Generalizing to 𝜇 0 = (1 -𝛼 0 ) 𝛿 𝑎 0 + 𝛼 0 𝛿 𝑏 0 (𝛼 ∈]0, 1[) and the same 𝜇 1 , the transport map is again trivial. There will be (1 -𝛼 0 ) shovel displacement between 𝑎 0 and 𝑎 1 and 𝛼 0 between 𝑏 0 and 𝑎 1 . The cost (1 -𝛼 0 ) 𝑐(𝑎 0 , 𝑎 1 ) + 𝛼 0 𝑐(𝑏 0 , 𝑎 1 ) is bilinear in the weight 𝛼 0 and the distance between points. Let us complicate this slightly to

𝜇 1 = (1 -𝛼 1 ) 𝛿 𝑎 1 + 𝛼 1 𝛿 𝑏 1 . If 𝛼 0 ∉ {𝛼 1 , 1 -𝛼 1 },
there is no measure-preserving Monge map. We can, however, find Kantorovich transport plans. They are 2 × 2 matrices assigning mass between couples in (𝑎 0 , 𝑏 0 ) × (𝑎 1 , 𝑏 1 ). The Kantorovich formulation is more general than the Monge formulation and subsumes Monge. See also Figure 1.1. 5 The notation 𝑓 , 𝛼 Ω will stand for the duality product ∫ Ω 𝑓 d𝛼 between continuous functions 𝑓 ∈ C(Ω) and probability measures 𝛼 ∈ P(Ω); the Ω subscript will often be omitted if we just consider a unit constant function. 6 A sequence of measures 𝜇 𝑛 converges to 𝜇 for the weak-★ topology if lim 𝑛→+∞ 𝑓 , 𝜇 𝑛 𝑋 = 𝑓 , 𝜇 𝑋 for all 𝑓 ∈ C(𝑋).

mass

𝑥 1 𝑥 0 (a) 𝛼 0 = 1 4 , 𝛼 1 = 1 4 mass 𝑥 1 𝑥 0 (b) 𝛼 0 = 1 4 , 𝛼 1 = 3 4 Figure 1.1. Blue, 𝜇 0 = 𝛼 0 𝛿 𝑎 0 + (1 -𝛼 0 ) 𝛿 𝑏 0 ; red, 𝜇 1 = 𝛼 1 𝛿 𝑎 1 + (1 -𝛼 1 ) 𝛿 𝑏 1 ; black, 𝜋 * .
The optimal plan cannot be reduced to a map when 𝛼 0 ≠ (𝛼 1 , 1 -𝛼 1 ).

Linear programming

The first and most obvious numerical method to attack the Kantorovich problem numerically is the linear programming approach. It assumes that the data are given in discrete form, for 𝑘 = 0, 1, 𝜇 𝑘 = 𝑁 𝑘 𝑖=1 𝑤 𝑘,𝑖 𝛿 𝑥 𝑘,𝑖 . The measures 𝜇 𝑘 may be empirical measures set on 𝑁 𝑘 points in 𝑋: 𝑋 𝑘 = {𝑥 𝑘,𝑖 } 𝑖=1, 𝑁 𝑘 or a discretization of probability densities on this discrete set representing a Cartesian grid, for instance. The weights {𝑤 𝑘,𝑖 } 𝑖=1, 𝑁 𝑘 are strictly positive and sum to 1.

Plugging the discrete data and support into (1.2)-(1.3), brackets become sums and we obtain the finite-dimensional linear program

𝜋 * ∈ arg inf 𝜋 ∈Π(𝜇 0 , 𝜇 1 ) ∑ (𝑖, 𝑗)∈𝑁 0 ×𝑁 1 𝑐(𝑥 0,𝑖 , 𝑥 1, 𝑗 ) 𝜋 𝑖, 𝑗 .
(1.4)

Note that Π(𝜇 0 , 𝜇 1 ) can be rewritten as a set of 𝑁 0 × 𝑁 1 matrices with constraints on sums of rows and

Π(𝜇 0 , 𝜇 1 ) := (𝜋 𝑖, 𝑗 ) ∈ M 𝑁 0 ×𝑁 1 (R + ) : ∑ 𝑗 ∈ [1, 𝑁 1 ]
𝜋 𝑖, 𝑗 = 𝑤 0,𝑖 and

∑ 𝑖 ∈ [1, 𝑁 0 ]
𝜋 𝑖, 𝑗 = 𝑤 1, 𝑗 .

(1.5)

Problem (1.4-1.5) is a discrete linear programming problem in 𝑁 0 ×𝑁 1 unknowns and 𝑁 0 + 𝑁 1 constraints. Numerical resolution with linear solvers which have cubic complexity in the number of unknowns is therefore out of reach for reasonable discretizations (typically 𝑁 𝑘 < 100 or even less depending on machine memory). As discussed in the example at the end of the previous subsection, Monge solutions are not guaranteed to exist: they will appear when 𝑁 0 = 𝑁 1 = 𝑁 and a permutation matrix 𝜋 𝑖, 𝑗=𝜎(𝑖) (𝜎 is a permutation of 1, 𝑁 ) is feasible. There is only one non-zero element per line and column and the permutation can be interpreted as a map. This implicitly constrains the weights to be able to match exactly {𝑤 0,𝑖 = 𝑤 1, 𝜎(𝑖) } 𝑖=1, 𝑁 .

When the weights are all equal (to 1/𝑁), the Birkhoff-von Neumann theorem (see Carlier 2021, §6.4, or Peyré andCuturi 2019, §2.3) states that the permutation matrices are the extreme points of Π(𝜇 0 , 𝜇 1 ). Problem (1.4)-(1.5) is then known as the linear assignment problem.

Assuming all weights {𝑤 𝑘,𝑖 } 𝑖=1, 𝑁 are integer multiples of the elemental mass 𝐸 := min 𝑘 {1/𝑁 𝑘 }, it is always possible to rewrite (1.4) as a linear assignment problem by piling Dirac masses at the same point, 𝑤 0,𝑖 𝛿 𝑥 0,𝑖 = 𝐸 𝑤 0,𝑖 /𝐸 𝑙=1 𝛿 𝑥 0,𝑖,𝑙 (and performing the same decomposition on 𝜇 1 ). The new Dirac masses are at the same location but counted 𝑤 0,𝑖 times. The weights are units but the size of the problem and the cost matrix may be larger. Optimal permutations 𝜎 * then characterize the support of 𝜋 * , supp 𝜋 * = {(𝑥 0,𝑖 , 𝑥 1, 𝜎 * (𝑖) )} 𝑖=1, 𝑁 , where 𝑁 ≥ 𝑁 is the new bigger number of unit Dirac masses.

The optimality over permutations translates into an important property called cyclical monotonicity (Carlier 2021, §6.4): for every subset {𝑥 0,𝑖 , 𝑥 1,𝑖 } 𝑖=1,𝐼 ∈ supp 𝜋 * , we have

∑ 𝑖 ∈ [1,𝐼 ] 𝑐(𝑥 0,𝑖 , 𝑥 1,𝑖 ) ≤ ∑ 𝑖 ∈ [1,𝐼 ] 𝑐(𝑥 0,𝑖+1 , 𝑥 1,𝑖 ), (1.6)
where by convention 𝑥 0,𝐼 +1 = 𝑥 0,1 . This property is in fact also a sufficient condition for the optimality.

If we apply (1.6) to a simple pair (𝑖, 𝑖 ) and use the quadratic displacement cost (1.1), we find (after simple computations) (𝑥 0,𝑖 -𝑥 0,𝑖 ) • (𝑥 1, 𝜎 * (𝑖) -𝑥 1, 𝜎 * (𝑖 ) ) ≥ 0.

(1.7) This is a simple monotonicity condition on the transport with the important consequence that the affine trajectories of the planned optimal transport paths 𝑥 0,𝑖 ↦ → 𝑥 1, 𝜎 * (𝑖) and 𝑥 0,𝑖 ↦ → 𝑥 1, 𝜎 * (𝑖 ) cannot intersect. It allows us to reconcile Monge's ordinary (not magic) shovel with the mathematical formulation. The optimality ensures that the map corresponds to a structurally feasible sequence of shovel transport: there will be no obstruction.

Kantorovich duality and the reflector problem

Kantorovich duality and semi-duality

Kantorovich duality is the application of linear programming duality theory to (1.2). For simplicity we assume 𝑋 0 and 𝑋 1 to be compact; for a general formulation, see Santambrogio (2015, Theorem 1.39). The dual Kantorovich problem is

(𝑢 * 0 , 𝑢 * 1 ) ∈ arg sup (𝑢 0 ,𝑢 1 )∈𝐶 𝐷 𝑢 0 , 𝜇 0 𝑋 0 + 𝑢 1 , 𝜇 1 𝑋 1 , (2.1)
where 𝑢 0 ⊕ 𝑢 1 (𝑥 0 , 𝑥 1 ) = 𝑢 0 (𝑥 0 ) + 𝑢 1 (𝑥 1 ) and

𝐶 𝐷 := {(𝑢 0 , 𝑢 1 ) ∈ C(𝑋 0 ) × C(𝑋 1 ) : 𝑢 0 ⊕ 𝑢 1 ≤ 𝑐}. (2.2)
The derivation of (2.1) is classical and can be found for instance in Santambrogio (2015, §1.2) or Carlier (2021, §6.4.1). The dual variables (𝑢 0 , 𝑢 1 ) are the Lagrange multipliers of the marginal constraints in (1.3). There is a nice economics interpretation in terms of optimal prices: 𝜇 0 and 𝜇 1 now represent two different types of goods, their mass and storage location. For some reason it is decided that the goods must be swapped, either via ground transport or by selling locally produced goods instead of moving the already available but misplaced ones. Before deciding on the method, the government needs to evaluate the feasibility of producing the goods locally. This depends on the maximal revenue a (single) producer may expect by selling the goods at price 𝑢 0 and 𝑢 1 respectively; this is given by (2.1). Of course selling will only happen if (2.2) is satisfied, that is, the sale price of a unit of mass on each side 𝑢 0 ⊕ 𝑢 1 is less than the fixed price 𝑐 of swapping a unit of mass through transport. Assuming maximizers (𝑢 * 0 , 𝑢 * 1 ) exist, the maximal revenue is the transport cost,

𝑐, 𝜋 * 𝑋 0 ×𝑋 1 = 𝑢 * 0 , 𝜇 0 𝑋 0 + 𝑢 * 1 , 𝜇 1 𝑋 1 , (2.3)
and the optimal prices saturate the constraint when 𝜋 * > 0:

𝑢 * 0 (𝑥 0 ) + 𝑢 * 1 (𝑥 1 ) = 𝑐(𝑥 0 , 𝑥 1 ) for all (𝑥 0 , 𝑥 1 ) ∈ supp 𝜋 * , 𝑢 * 0 (𝑥 0 ) + 𝑢 * 1 (𝑥 1 ) < 𝑐(𝑥 0 , 𝑥 1 ) else.

(2.4)

In the literature (𝑢 * 0 , 𝑢 * 1 ) are known as the Kantorovich potentials. They are always defined up to a constant as (𝑢 * 0 + 𝐶, 𝑢 * 1 -𝐶) is also a solution for any 𝐶 ∈ R. The optimality condition (2.4) imposes a strong structure on the maximizers, in that they are 𝑐-transforms of one another:

𝑢 * 0 = (𝑢 * 1 ) 𝑐 and 𝑢 * 1 = (𝑢 * 0 ) 𝑐 . (2.5)
As (2.5) is a necessary condition for optimality, the semi-dual formulation idea is to leverage 𝑢 1 = (𝑢 0 ) 𝑐 to optimize only on 𝑢 0 . This is a key part of the theory 8 ⊕ is the direct sum. 9 A quick summary: 𝑓 𝑐 , the 𝑐-transform of a function 𝑓 : 𝑋 0 → R ∪ +∞ is defined as 𝑓 𝑐 (𝑥 1 ) := inf 𝑥 0 ∈𝑋 0 𝑐(𝑥 0 , 𝑥 1 ) -𝑓 (𝑥 0 ) (note that the definition depends on the domain of 𝑓 ). The domain of 𝑓 𝑐 is 𝜕 𝑐 𝑓 (𝑋 0 ) = ∪ 𝑥 0 ∈𝑋 0 {𝜕 𝑐 𝑓 (𝑥 0 )}, where 𝜕 𝑐 𝑓 is the 𝑐-subdifferential of 𝑓 , i.e. the multivalued map 𝜕 𝑐 𝑓 (𝑥 0 ) = {𝑥 1 ∈ 𝑋 : 𝑐(𝑥 0 , 𝑥 1 ) -𝑓 (𝑥 0 ) ≤ 𝑐(𝑥 0 , 𝑥 1 ) -𝑓 (𝑥 0 ) for all 𝑥 0 ∈ 𝑋 0 }, that is, where the 𝑐-transform is defined; this set may be empty. A function 𝑔 is 𝑐-concave if there exists 𝑓 such that 𝑔 = 𝑓 𝑐 and then the domain is given by Dom(𝑔 𝑐 ) = 𝜕 𝑐 𝑓 𝑐𝑐 (𝑋 0 ). If 𝑓 is also 𝑐-concave then 𝑓 𝑐𝑐 = 𝑓 .

as the set of potentials {(𝑢 0 , 𝑢 𝑐 0 )} is compact while 𝐶 𝐷 is not (see Santambrogio 2015, §1.6). The other condition 𝑢 0 = (𝑢 1 ) 𝑐 requires 𝑢 0 to be 𝑐-concave, hence 𝑢 0 = 𝑢 𝑐𝑐 0 = ((𝑢 0 ) 𝑐 ) 𝑐 = (𝑢 1 ) 𝑐 . Finally the domain of (𝑢 0 ) 𝑐 must match 𝑋 1 , the support of 𝜇 1 , where the mass is supported.

The semi-dual optimal transportation problem reads sup

𝑢 0 ∈𝐶 SD SD(𝑢 0 ) := 𝑢 0 , 𝜇 0 𝑋 0 + 𝑢 𝑐 0 , 𝜇 1 𝑋 1 , (2.6)
where

𝐶 SD := {𝑢 0 is 𝑐-concave and Dom(𝑢 𝑐 0 ) = 𝑋 1 }.
(2.7)

We will comment on the set of constraints in the simpler quadratic cost formalism in Section 2.3. Expanding the 𝑐-transform in (2.6), we get

SD(𝑢 0 ) = 𝑢 0 , 𝜇 0 𝑋 0 + 𝑥 1 ↦ → inf 𝑥 0 ∈𝑋 0 𝑐(𝑥 0 , 𝑥 1 ) + 𝑢 0 (𝑥 0 ) , 𝜇 1 𝑋 1 . (2.8)
The concavity of the cost function 𝑢 0 → SD(𝑢 0 ) follows directly from the concavity of 𝑢 0 (•) ↦ → 𝑢 𝑐 0 (•) := inf 𝑥 0 {𝑐(𝑥 0 , •) + 𝑢 0 (𝑥 0 )} and the linearity of (2.6). We move from a linear to a concave program at the price of constraint (2.7). The computation of the critical point of 𝑢 0 ↦ → SD(𝑢 0 ) is the basis of the Monge-Ampère and semi-discrete formulations detailed below.

Existence of Monge map uniqueness of Kantorovich plans

The optimality condition (2.4) is also central to establishing the existence and uniqueness of Monge maps. This is one of the fundamental contributions of [START_REF] Brenier | Polar factorization monotone rearrangement of vector-valued functions[END_REF] in the Euclidean quadratic cost case (1.1) and it was generalized by [START_REF] Gangbo | The geometry of optimal transportation[END_REF]. The idea is straightforward (once you have it!) and relies on the twist condition on the cost: for all 𝑥 0 ∈ 𝑋 0 , the map 𝑥 1 ∈ 𝑋 1 ↦ → 𝐷 𝑥 0 𝑐(𝑥 0 , 𝑥 1 ) ∈ T 𝑥 0 𝑋 0 (2.9) is one-to-one. Using (2.9), (2.4) can be solved as an implicit equation in 𝑥 1 ; the solution for all 𝑥 0 ∈ 𝑋 0 defines the unique optimal Monge map

𝑇 : 𝑥 0 → 𝑥 * 1 := {{𝑥 1 ↦ → 𝐷 𝑥 0 𝑐(𝑥 0 , 𝑥 1 )}} -1 (𝐷 𝑥 0 𝑢 * 0 (𝑥 0 )).
(2.10)

Note that the gradient gets rid of the additive constant underdetermination. Precise theorems can be found in the literature (e.g. Theorem 2.9 of (Dafni, McCann and Stancu 2013, §6)). We give a few important consequences of this result. First, condition (2.9) implicitly assumes that 𝑐 is continuous and at least differentiable with respect to 𝑥 ∈ 𝑋. Defining this differential for all 𝑥 0 ∈ 𝑋 0 rules out discrete measures: 𝜇 0 has to be absolutely continuous with respect to 1 𝑋 , the Lebesgue measure on 𝑋 (remember the Dirac counterexample in Section 1.2). The twist condition (2.9) is obviously satisfied by the Euclidean quadratic cost (1.1), but also by the more complex reflector cost presented at the end of this section; it only needs to hold between 𝑋 0 and 𝑋 1 .

The second important consequence of (2.4) is that the unique transport map is sufficient to characterize the unique Kantorovich solution. For all Borel sets 𝐴 0,1 ⊂ 𝑋 0,1 and because 𝜋 * , while satisfying (1.3), vanishes outside the graph {𝑥 0 , 𝑥 1 = 𝑇(𝑥 0 )}, 𝜋 * (𝐴 0 , 𝐴 1 ) = 𝜋 * (𝐴 0 ∩ 𝑇 -1 (𝐴 1 ), 𝑋 1 ) = 𝜇 0 (𝐴 0 ∩ 𝑇 -1 (𝐴 1 )).

(2.11)

In compact form, (2.11) is usually written 𝜋 * = (Id 𝑋 0 , 𝑇)#𝜇 0 , where Id 𝑋 0 is the identity map on 𝑋 0 and (Id 𝑋 0 , 𝑇) := {𝑥 0 ↦ → (𝑥 0 , 𝑇(𝑥 0 ))} is a map from 𝑋 0 to 𝑋 0 × 𝑋 1 . See Figure 2.1 for an example.

Monge-Ampère solutions

For the quadratic cost (1.1), 𝑐-concavity can be expressed, after a transformation of the potentials, using the classical Legendre-Fenchel transform (denoted (•) ★ and not to be mistaken for the optimizer notation (•) * ):

𝑥 0 ↦ → 1 2 𝑥 0 -. 2 -𝜙(𝑥 0 ) ★ (𝑥 1 ) := sup 𝑥 0 ∈𝑋 𝑥 1 • 𝑥 1 - 1 2 𝑥 0 2 -𝜙(𝑥 0 ) = 1 2 𝑥 1 2 -𝜙 𝑐 (𝑥 1 ).
(2.12) 13 The Legendre-Fenchel transform 𝑓 ★ of 𝑓 for a function 𝑓 : 𝑋 0 → R ∪ +∞ is defined by 𝑓 ★ (𝑥 1 ) := sup 𝑥 0 ∈𝑋 0 (𝑥 0 • 𝑥 1 ) -𝑓 (𝑥 0 ) (note that the definition depends on the domain of 𝑓 ); 𝑓 ★ is always convex. If 𝑓 is convex , ( 𝑓 ★ ) ★ = 𝑓 , 𝑓 and 𝑓 ★ are a.e. differentiable and 𝐷 𝑓

• 𝐷 𝑓 ★ = 𝐷 𝑓 ★ • 𝐷 𝑓 = Id.
The functions (2.13) correspond to the seminal setting of [START_REF] Brenier | Polar factorization monotone rearrangement of vector-valued functions[END_REF]. To optimize (2.6) over the {𝜙 0,1 } instead of the {𝑢 0,1 }, change the maximization into a minimization and add the constant second moments of 𝜇 0 and 𝜇 1 (omitted below): (2.14) The connection between convexity and the Legendre-Fenchel transform are well known (see footnote 13). In this setting (2.7) becomes

𝜙 𝑖 = 𝑥 𝑖 ↦ → 1 2 𝑥 𝑖 2 -𝑢 𝑖 (𝑥) , 𝑖 = 0, 1,
𝜙 * 0 := arg inf 𝜙 0 ∈𝐶 MA SD(𝜙 0 ), SD(𝜙 0 ) := 𝜙 0 , 𝜇 0 𝑋 0 + 𝜙 ★ 0 , 𝜇 1 𝑋 1 .
𝐶 MA = {𝜙 0 convex and 𝜕𝜙 0 (𝑋 0 ) = 𝑋 1 }, (2.15)
where

𝜕𝜙 0 (𝑥 0 ) := {𝑥 1 , 𝑥 0 • 𝑥 1 -𝜙 0 (𝑥 0 ) ≥ 𝑥 0 • 𝑥 1 -𝜙 0 (𝑥 0 ) for all 𝑥 0 ∈ 𝑋 0 } (2.16)
is the usual subgradient of 𝜙 0 at 𝑥 0 (and 𝜕𝜙 0 (𝑋 0 ) = {𝜕𝜙 0 (𝑥 0 )} 𝑥 0 ∈𝑋 0 ). As 𝜙 * 0 is convex, subgradients are gradients almost everywhere. This defines the singlevalued Monge map (2.10)

𝑥 0 ↦ → 𝑥 1 = 𝑇(𝑥 0 ) = 𝐷𝜙 * 0 (𝑥 0 ) = Id 𝑋 0 -𝐷𝑢 * 0 (𝑥 0 ).
(2.17)

Likewise (𝜙 * 0 ) ★ is convex and its domain is constrained to be 𝑋 1 . For almost every 𝑥 1 ∈ 𝑋 1 there is a unique 𝑥 0 = 𝐷(𝜙 * 0 ) ★ (𝑥 1 ) achieving the maximum in the Legendre-Fenchel transform:

(𝜙 * 0 ) ★ (𝑥 1 ) = 𝑥 1 • 𝐷(𝜙 * 0 ) ★ (𝑥 1 ) -𝜙 * 0 • 𝐷(𝜙 * 0 ) ★ (𝑥 1
+ 𝛿𝜙 0 ) -SD(𝜙 * 0 ) = 𝛿𝜙 0 , 𝜇 0 𝑋 0 -𝛿𝜙 0 • 𝐷(𝜙 * 0 ) ★ , 𝜇 1 𝑋 1 + 𝑂( 𝛿𝜙 0 2 ). (2.19)
Because of the change of variables (2.13), the maximization is now a minimization and SD is convex instead of concave. We conclude that the optimal potential satisfies for all 𝛿𝜙 0 the equation

𝛿𝜙 0 , 𝜇 0 𝑋 0 -𝛿𝜙 0 • 𝐷(𝜙 * 0 ) ★ , 𝜇 1 𝑋 1 = 0.
(2.20) 14 Use inf{•} =sup -{•}. 15 We give a simplified formulation. Let Φ : (𝑥 0 , 𝑥 1 ) ∈ 𝑋 0 × 𝑋 1 → R be continuous and convex in 𝑥 1 for all 𝑥 0 and 𝑋 0 be compact. Set Ψ(𝑥 1 ) = sup 𝑥 0 Φ(𝑥 0 , 𝑥 1 ). Then, if {𝑥 * 0 := arg sup 𝑥 0 Φ(𝑥 0 , 𝑥 1 )} is a singleton, Ψ is differentiable at 𝑥 1 and

𝐷 𝑥 1 Ψ(𝑥 1 ) = 𝜕 𝑥 1 Φ(𝑥 * 0 , 𝑥 1 )
. The assumptions on Φ may be considerably relaxed: see Carlier (2021, §5.3).

J.-D. B

The change of variable 𝑥 1 = 𝐷𝜙 * 0 (𝑥 0 ) gives (again formally)

𝛿𝜙 0 , 𝜇 0 𝑋 0 -𝛿𝜙 0 , (𝜇 1 • 𝐷𝜙 * 0 ) det(𝐷 2 𝜙 * 0 ) 𝑋 0 = 0. (2.21)
Therefore the optimal transportation potential, together with the constraint (2.15), may be understood in a weak sense as a solution of the (so-called) second boundary value problem for the Monge-Ampère equation:

det(𝐷 2 𝜙 0 ) = 𝜇 0 𝜇 1 • 𝐷𝜙 0 on 𝑋 0 . (2.22)
The link with the Monge-Ampère operator has triggered an important theoretical activity on the regularity theory of such equations; see [START_REF] Caffarelli | The regularity of mappings with a convex potential[END_REF] and [START_REF] Figalli | The Monge-Ampère Equation and its Applications[END_REF]. A short summary is as follows. Assuming 𝑋 1 is convex (this is imposed by 𝜕𝜙 0 (𝑋 0 ) = 𝑋 1 ) and 𝜇 0,1 bounded above and below away from 0, we gain two orders of regularity in Hölder spaces for 𝜙 0 with respect to the data 𝜇 0,1 , as is classically expected of second-order elliptic equations. The bounds on the densities force the right-hand side of (2.22) to be strictly positive and 𝜙 0 to be strictly convex. This is necessary for the variational formulation (2.21) to allow 𝜙 * 0 + 𝛿𝜙 0 to remain in 𝐶 for all 'small' 𝛿𝜙 0 (s).

The PDE (2.22)-(2.15) may be discretized and the resulting nonlinear system solved with a damped Newton's method yielding, experimentally at least, a linear cost (see [START_REF] Benamou | Numerical solution of the optimal transportation problem using the Monge-Ampère equation[END_REF][START_REF] Benamou | Minimal convex extensions and finite difference discretisation of the quadratic Monge-Kantorovich problem[END_REF], and the references therein ). The discretization is delicate: it has to preserve a notion of strict convexity to keep a well-defined Hessian for the discretization of the convex 𝜙 0 ↦ → SD(𝜙 0 ), i.e. the Jacobian of (2.22). The Newton's method damping ensures the right-hand side of (2.22) remains strictly positive; more on this can be found in Benamou, Collino and Mirebeau (2016d).

Finally let us mention that the Monge-Ampère formulation has been generalized to other costs satisfying (2.9); see Section 2.5 below for an example.

Semi-discrete solutions

When 𝜇 1 ∈ L 1 (𝑋 1 ) is simply integrable and 𝜇 0 = 𝑛∈ [1, 𝑁 ] 𝛼 𝑛 𝛿 𝑥 0,𝑛 is a finite sum of weighted Dirac masses, (2.6) becomes an unconstrained finite-dimensional concave optimization problem called semi-discrete optimal transportation (not to be confused with semi-dual).

We discuss the two components in (2.6) separately. First we have

𝑢 0 , 𝜇 0 𝑋 0 = ∑ 𝑛∈ [1, 𝑁 ]
𝛼 𝑛 𝑢 0 (𝑥 0,𝑛 ), so we only need to set the 𝑁 prices {𝑢 0,𝑛 } := {𝑢 0 (𝑥 0,𝑛 )} where the mass is located. Now, by construction and for all 𝑥 1 ∈ 𝑋 1 ,

𝑢 𝑐 0 (𝑥 1 ) = inf 𝑛∈ [1, 𝑁 ] {𝑐(𝑥 0,𝑛 , 𝑥 1 ) -𝑢 0,𝑛 } (2.23)
(remember that 𝑋 0 = supp(𝜇 0 ) = {𝑥 0,𝑛 }). We could plug this expression directly into (2.6), but it is helpful to introduce the sets of Laguerre cells. They depend on the vector 𝑢 0 = {𝑢 0,𝑛 } ∈ R 𝑁 , for all 𝑚 ∈ [1, 𝑁]:

Lag 𝑚 (𝑢 0 ) := {𝑥 1 ∈ 𝑋 1 : 𝑢 0,𝑚 -𝑐(𝑥 0,𝑚 , 𝑥 1 ) ≤ 𝑢 0,𝑛 -𝑐(𝑥 0,𝑛 , 𝑥 1 ) for all 𝑛 ∈ [1, 𝑁]}.
(2.24)

In our economics illustration, Lag 𝑚 is the subset of 𝑋 1 where the price is set to 𝑢 𝑐 0 (𝑥 1 ) = 𝑐(𝑥 0,𝑚 , 𝑥 1 ) -𝑢 0,𝑚 (the infimum is reached for 𝑚 in (2.23)). We now plug this into (2.6) by using the decomposition

𝑋 1 = ∪ 𝑛∈ [1, 𝑁 ] Lag 𝑛 (𝑢 0 ): 𝑢 0 , 𝜇 0 𝑋 0 + 𝑢 𝑐 0 , 𝜇 1 𝑋 1 = ∑ 𝑛∈ [1, 𝑁 ] 𝛼 𝑛 𝑢 0,𝑛 + 𝑢 𝑐 0 , 𝜇 1 ∪ 𝑛∈ [1, 𝑁 ] Lag 𝑛 = ∑ 𝑛∈ [1, 𝑁 ] 𝛼 𝑛 𝑢 0,𝑛 + ∑ 𝑛∈ [1, 𝑁 ] 𝑐(𝑥 0,𝑛 , •) -𝑢 0 𝑛 , 𝜇 1 Lag 𝑛 = ∑ 𝑛∈ [1, 𝑁 ] (𝛼 𝑛 -𝜇 1 (Lag 𝑛 )) 𝑢 0,𝑛 + ∑ 𝑛∈ [1, 𝑁 ] 𝑐(𝑥 0,𝑛 , •), 𝜇 1 Lag 𝑛 . (2.25)
The last term is constant and can be removed from the maximization. The second term expresses the mass balance of what is swapped between Lag 𝑛 and 𝑥 1,𝑛 and the actual mass to be transported 𝛼 𝑛 . Optimality of this concave program is therefore obtained by solving the set of 𝑁 nonlinear equations

𝛼 𝑛 = 𝜇 1 (Lag 𝑛 ), 𝑛 ∈ [1, 𝑁].
(2.26)

The constraints (2.7) are automatically taken into account by this formulation. First, the 𝑋 1 support constraint is implicitly enforced by the computation of the Laguerre cells (2.24). The 𝑐-concavity condition is more subtle and can be transformed into a local constraint. If

Lag 𝑛 (𝑢 0 ) ≠ ∅, 𝑛 ∈ [1, 𝑁],
(2.27)

then 𝑢 𝑐𝑐 0 (𝑥 0,𝑛 ) = 𝑢 0,𝑛 for all 𝑛 ∈ [1, 𝑁] and we can use the function 𝑢 𝑐𝑐 0 , which is by construction 𝑐-concave.

The formulation, analysis and solution method is explained in great detail in [START_REF] Mérigot | Optimal transport: Discretization and algorithms[END_REF]. A more computationally oriented survey is that of [START_REF] Lévy | Notions of optimal transport theory and how to implement them on a computer[END_REF]. In the case of the Euclidean quadratic cost (1.1), Laguerre cells can be computed in 𝑂(𝑁 log(𝑁)) operations. The integration of the measure 𝜇 0 on each cell is the only numerical approximation in the method. Finally, [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF] has shown that (2.26) was amenable to a damped Newton's method forcing the volume condition (2.27) and yielding a linear cost optimal transportation solver for the Euclidean quadratic cost.

The far field reflector cost

Analysts and differential geometers have gone a long way to extend optimal transportation to Riemannian manifolds with a significant impact on pure mathematics. The place to delve further is [START_REF] Villani | Optimal Transport: Old and New[END_REF]. It is also important for numerical analysts. We illustrate this generalization with an optimal transportation application set on the (𝑑 -1)-dimensional unit sphere S 𝑑-1 . Here 𝑋 0 ⊂ (S 𝑑-1 ) + and 𝑋 1 ⊂ (S 𝑑-1 ) -are connected domains in the northern and southern hemispheres respectively; 𝜇 0 ∈ 𝑋 0 represents a given illuminance from a point source (the R 𝑑 origin); rays carry 𝜇 0 (𝑥 0 ) light intensity in the direction 𝑥 0 . The other prescribed marginal is the illumination: 𝜇 1 (𝑥 1 ) is the light intensity in the direction 𝑥 1 . The reflector 'freeforming' problem is to find a surface in R 𝑑 such that the specular law of reflection maps 𝜇 0 to 𝜇 1 .

This problem has an elegant optimal transportation formulation introduced independently by [START_REF] Glimm | Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem[END_REF] and [START_REF] Wang | On the design of a reflector antenna II[END_REF], as follows.

Let us apply the Monge theory to the displacement cost:

𝑐(𝑥 0 , 𝑥 1 ) = -log(1 -𝑥 0 • 𝑥 1 ), (𝑥 0 , 𝑥 1 ) ∈ (𝑋 0 × 𝑋 1 ).
(2.28)

The TWIST condition (2.9) is satisfied; note that the restriction to 𝑋 0,1 is essential to preserve the regularity of the cost. A transport map exists, given by (2.10).

Taking the exponential of the equation/inequation over the support of the optimal transportation plan in (2.4), we get

e -𝑢 * 1 (𝑇 (𝑥 0 )) 1 -𝑥 0 • 𝑇(𝑥 0 ) = e 𝑢 * 0 (𝑥 0 ) ≤ e -𝑢 * 1 (𝑥 1 )) 1 -𝑥 0 • 𝑥 1 for all 𝑥 1 ∈ 𝑋 1 .
(2.29)

We now notice that

𝑥 0 ∈ 𝑋 0 ⊂ (S 𝑑-1 ) + → Par 𝑥 1 (𝑥 0 ) := 𝑥 0 e -𝑢 * 1 (𝑥 1 )) 1 -𝑥 0 • 𝑥 1 (2.30)
is a family of parabolic reflectors in R 𝑑 with axis 𝑥 1 ∈ 𝑋 1 , focal at the origin and focal length 1 2 e -𝑢 * 1 (𝑥 1 ) . A parabolic mirror reflects all incoming light rays in its axis direction. Equation (2.29) therefore shows that a ray shot in direction 𝑥 0 will first touch the parabola Par 𝑇 (𝑥 0 ) after travelling the distance

e 𝑢 * 0 (𝑥 0 ) = inf 𝑥 1 ∈𝑋 1 e -𝑢 * 1 (𝑥 1 )) 1 -𝑥 0 • 𝑥 1 (2.31)
and then be reflected in the 𝑥 * 1 = 𝑇(𝑥 0 ) direction. The transport map enforces the law of specular reflection of the surface Refl = {𝑥 0 e 𝑢 * 0 (𝑥 0 ) , 𝑥 0 ∈ 𝑋 0 }.

By construction 𝜇 1 = 𝑇#𝜇 0 , so the illuminance and illumination constraints are satisfied. The simplest example corresponds to a single Dirac target 𝜇 1 = 𝛿 𝑎 1 , the solution being the parabolic mirror Par 𝑎 1 : all the light from 𝜇 0 is reflected in the 𝑎 1 direction. The case of two Dirac masses

𝜇 1 = 𝛼 1 𝛿 𝑎 1 + (1 -𝛼 1 )
𝛿 𝑏 1 is sufficient to get the general idea. According to (2.31) rays will hit the closest parabola, and the mirror is the inf-envelope of Par 𝑎 1 and Par 𝑏 1 . The focal lengths, depending on the 'prices' (𝑢 1 (𝑎 1 ), 𝑢 1 (𝑏 1 )), determine their unique intersection point. It partitions the support of 𝜇 0 into the rays reflecting in the 𝑎 1 and 𝑏 1 direction respectively. It is adjusted so that the energy carried splits according to the Dirac weights (𝛼 1 , 1-𝛼 1 ). This is easily generalized to 𝑀 Dirac masses and, as the reflection depends on the local tangent to the reflector, to densities. The corresponding 𝑐-Monge-Ampère equation has been solved numerically using a B-spline collocation method and a multi-scale approach in [START_REF] Brix | Solving the Monge-Ampère equations for the inverse reflector problem[END_REF]. The method is tested by providing an image on a projection plane, computing the corresponding reflector, resimulating via ray-tracing, and finally comparing the two pictures: see Figure 2.2.

Semi-discrete optimal transportation has been applied to the reflector problem in de Castro et al. (2016). Equation (2.31) is simply the 𝑢 1 𝑐-transform. When 𝜇 1 is discrete we have a finite collection of parabolae, each sending rays to one fixed direction 𝑥 1,𝑛 . The Laguerre cells contain the ray directions of the source touching the corresponding parabola. Adjusting the focal length e -𝑢 1 (𝑥 1 ) amounts to moving the parabola closer to or further away from the origin such that the amount of light energy carried by the rays in the Laguerre cell is exactly the weight of the corresponding target direction; see Figure 2.3.

Use of the Sinkhorn algorithm and the entropic regularization method has been investigated in [START_REF] Benamou | An entropic optimal transport numerical approach to the reflector problem[END_REF].

Wasserstein distance and entropic interpolation

Displacement interpolation and W 2 metric

In this section we will restrict ourselves to the classical Euclidean quadratic cost (1.1) and show that the optimal value function, or Monge optimal transport 'work'

W 2 (𝜇 0 , 𝜇 1 ) := √ 𝑐, 𝜋 * 𝑋 0 ×𝑋 1 = √ ∫ 𝑇 * (𝑥 0 ) -𝑥 0 2 𝑑𝜇 0 (𝑥 0 ), (3.1)
defines a distance, in the mathematical sense on P(𝑋), known as the Wasserstein-2 distance. Identity (nothing moves) and symmetry (transport back) are immediate.

The triangle inequality is a consequence of the same property for 𝑐, the displacement cost. Let 𝑇 𝜇 0 →𝜌 and 𝑇 𝜌→𝜇 1 , respectively, be the optimal Monge maps from 𝜇 0 to 𝜌 and 𝜌 to 𝜇 1 , and assume they are all absolutely continuous densities of probabilities.

Then 𝑇 𝜌→𝜇 1 • 𝑇 𝜇 0 →𝜌 is an admissible transport from 𝜇 0 to 𝜇 1 . We have

W 2 (𝜇 0 , 𝜇 1 ) = √ ∫ 𝑇 𝜌→𝜇 1 (𝑇 𝜇 0 →𝜌 (𝑥 0 )) -𝑥 0 2 𝑑𝜇 0 (𝑥 0 ) ≤ √ ∫ 𝑇 𝜌→𝜇 1 (𝑇 𝜇 0 →𝜌 (𝑥 0 )) -𝑇 𝜇 0 →𝜌 (𝑥 0 ) 2 𝑑𝜇 0 (𝑥 0 ) + √ ∫ 𝑇 𝜇 0 →𝜌 (𝑥 0 ) -𝑥 0 2 𝑑𝜇 0 (𝑥 0 ).
The last term is W 2 (𝜇 0 , 𝜌) and the change of variable 𝑥 1 = 𝑇 𝜇 0 →𝜌 (𝑥 0 ) in the first one gives W 2 (𝜌, 𝜇 1 ) (see Santambrogio 2015, §5.1).

We are now going to make a detour through displacement interpolation and an incursion into multi-marginal optimal transportation.

We assume that 𝑋 0 and 𝑋 1 are subsets of a larger 𝑋. The interpolation in P(𝑋) with respect to the transport cost is

𝜌 𝜏 := arg inf 𝜌∈P(𝑋 ) {(1 -𝜏) W 2 2 (𝜇 0 , 𝜌) + 𝜏 W 2 2 (𝜌, 𝜇 1 )}, 𝜏 ∈ (0, 1). (3.2)
The mapping 𝜌 ↦ → W 2 2 (𝜇 0 , 𝜌) is continuous, differentiable and convex on P(𝑋). These properties are further discussed in Section 7.2. The interpolation 𝜌 𝜏 is well-defined. It is the distribution of mass closest to 𝜇 0 and 𝜇 1 in the sense of the weighted average transport cost. We have two nested optimization problems in (3.2), that is, inf

𝜌∈P(𝑋 ) (1 -𝜏) inf 𝜋 0→𝜏 ∈Π(𝜇 0 ,𝜌) 𝑐, 𝜋 0→𝜏 𝑋 0 ×𝑋 + 𝜏 inf 𝜋 𝜏→1 ∈Π(𝜌, 𝜇 1 )
𝑐, 𝜋 𝜏→1 𝑋 ×𝑋 1 .

(3.

3)

The plans 𝜋 0→𝜏 and 𝜋 𝜏→1 specify the quantity of mass flowing from 𝑥 0 ∈ 𝑋 0 to some 𝑥 𝜏 ∈ 𝑋 and then from 𝑥 𝜏 ∈ 𝑋 to 𝑥 1 ∈ 𝑋 1 . We can rewrite this problem in the Kantorovich spirit. Instead of looking for coupling, let us search for 'triplings': the quantity of mass, denoted 𝜋 3 , flowing across triplets (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) ↦ → 𝜋 3 (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) in the bigger space P(𝑋 0 × 𝑋 × 𝑋 1 ):

inf

𝜌∈P(𝑋 ) inf 𝜋 3 ∈Π 3 (𝜇 0 ,𝜌, 𝜇 1 ) 𝑐 3 , 𝜋 3 𝑋 0 ×𝑋 ×𝑋 1 , (3.4)
, 23

where

Π 3 (𝜇 0 , 𝜌, 𝜇 1 ) (3.5) := {𝜋 3 ∈ P(𝑋 0 × 𝑋 × 𝑋 1 ) : 𝑃 𝑋 0 #𝜋 3 = 𝜇 0 , 𝑃 𝑋 1 #𝜋 3 = 𝜇 1 and 𝑃 𝑋 #𝜋 3 = 𝜌}
is the set of 3-marginal plans and 𝑐 3 is given by (3.6) below. For all couples of 2-plans 𝜋 0→𝜏 ∈ Π(𝜇 0 , 𝜌) and 𝜋 𝜏→1 ∈ Π(𝜌, 𝜇 1 ) there exists a 3-plan 𝜋 3 ∈ Π 3 (𝜇 0 , 𝜌, 𝜇 1 ) such that 𝜋 0→𝜏 = 𝑃 𝑋 1 #𝜋 3 and 𝜋 𝜏→1 = 𝑃 𝑋 0 #𝜋 3 (see Lemma 5.5 in Santambrogio 2015) andconversely. Problems (3.3) and(3.4) are equivalent. The intermediate marginal 𝜌 support can be determined with the help of the displacement cost:

𝑐 3 (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) = (1 -𝜏) 1 2 𝑥 𝜏 -𝑥 0 2 + 𝜏 1 2 𝑥 1 -𝑥 𝜏 2 . (3.6)
If 𝜋 3 (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) > 0 and because 𝑐 3 is a strictly convex function of 𝑥 𝜏 , mass travelling from 𝑥 0 to 𝑥 1 necessarily goes through its minimum 𝑥 𝜏 = 𝑥 0 + 𝜏(𝑥 1 -𝑥 0 ). Straight trajectories are optimal. The 3-marginal cost 𝑐 3 therefore simplifies to the 2-marginal cost 𝑐 (1.1) for the mass carried by minimizers 𝜋 * 3 of (3.5)-(3.4). The integration

𝑃 𝑋 0 ×𝑋 1 #𝜋 * 3 = 1, 𝜋 3 𝑋 coincides with the Kantorovich 2- plan 𝜋 * 0→1 solution of (1.2). The 3-plan 𝜋 * 3 is fully determined by 𝜋 * 3 (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) = 𝜋 * 0→1 (𝑥 0 , 𝑥 1 ) if 𝑥 𝜏 = 𝑥 0 + 𝜏(𝑥 1 -𝑥 0 )
and 0 else; remember that optimal transport paths do not cross. The outer minimization in (3.4) can be eliminated and the problem relaxed to the simpler inf

𝜋 3 ∈Π 2 (𝜇 0 , 𝜇 1 ) 𝑐 3 , 𝜋 3 𝑋 0 ×𝑋 ×𝑋 1 , (3.7)
where

Π 2 (𝜇 0 , 𝜇 1 ) := {𝜋 3 ∈ P(𝑋 0 × 𝑋 × 𝑋 1 ) : 𝑃 𝑋 0 #𝜋 3 = 𝜇 0 , 𝑃 𝑋 1 #𝜋 3 = 𝜇 1 }. (3.8)
We are back to the 2-marginal constraints, but optimize on the 3-marginal plans. The solutions of (3.3) are recovered via marginal integrations: Monge level, McCann (1997) has shown that, given the map 𝑥 0 ↦ → 𝑇 0→1 (𝑥 0 ) between 𝜇 0 and 𝜇 1 , the map

𝜌 𝜏 = 𝑃 𝑋 #𝜋 * 3 , 𝜋 * 0→𝜏 = 𝑃 𝑋 0 ×𝑋 #𝜋 * 3 and 𝜋 * 𝜏→1 = 𝑃 𝑋 ×𝑋 1 #𝜋 * 3 . At the
𝑥 0 ↦ → 𝑇 0→𝜏 (𝑥 0 ) := 𝑥 0 + 𝜏 (𝑇 0→1 (𝑥 0 ) -𝑥 0 ) ( = 𝑥 𝜏 above) (3.9)
is the Monge map between 𝜇 0 and 𝜌 𝜏 , the solution of the interpolation problem (3.2). He called this process displacement interpolation. It can be used to generate a smooth dynamic 𝜏 ∈ (0, 1) → 𝜌 𝜏 sequence of probability measures morphing 𝜇 0 to 𝜇 1 continuously. This process will, in particular, displace the support of 𝜌 𝜏 smoothly and 'horizontally'. Assuming 𝜇 0 and 𝜇 1 are integrable, it is interesting to compare with the L 2 interpolation of functions, that is,

𝜌 𝜏 := arg inf 𝜌 {(1 -𝜏) 𝜇 0 -𝜌 2 L 2 + 𝜏 𝜌 -𝜇 1 2 L 2 },
which yields a pointwise interpolated sequence 𝜌 𝜏 = (1 -𝜏) 𝜇 0 + 𝜏 𝜇 1 . The mass will only move 'vertically'. Figure 3.1 illustrates these fundamentally different behaviours.

An important remark is that the curve 𝜏 ∈ (0, 1) ↦ → 𝜌 𝜏 is a geodesic curve for the W 2 distance. For any two times (𝜏 0 , 𝜏 1 ) ∈ [0, 1] 2 , 𝜏 ∈]𝜏 0 , 𝜏 1 [→ 𝜌 𝜏 is the displacement interpolation between 𝜌 𝜏 0 and 𝜌 𝜏 1 . The piece of the map (3.9) 𝑥 𝜏 0 → 𝑥 𝜏 1 is also the Monge map between 𝜌 𝜏 0 and 𝜌 𝜏 1 .

Displacement interpolation is harder to compute than L 2 interpolation, but its ability to model smooth displacements has caught the attention of many researchers in signal processing. All is far from perfect, though. Signals are not necessarily probability distributions. The measure-preserving/mass-conservation property also rules out any 'rigidity' or 'geometric' constraint on supp (𝜌 𝜏 ), but it may be a no-go for many morphing applications. We discuss these restrictions further in Sections 4 and 7.

Finally, let us mention attempts at formulating higher-order optimal transportation interpolations, e.g. Benamou, Gallouët and Vialard (2019c) and [START_REF] Chen | Measure-valued spline curves: An optimal transport viewpoint[END_REF].

Sinkhorn algorithm

A method based on an 'entropic regularization' [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]) was introduced to speed up optimal transportation computations at the price of an approximation.

The method comes with a parallelization-friendly computational method usually known as the Sinkhorn algorithm. It can be traced back to much older works. We will present a statistical physics interpretation in Section 6 but it is also linked to exponential barrier functions in discrete optimization [START_REF] Cominetti | Asymptotic analysis of the exponential penalty trajectory in linear programming[END_REF].

A quick and somewhat dirty derivation of the algorithm (the literature usually starts with the primal problem) starts from the dual problem (2.1) in its equivalent formulation sup

{𝑢 0 ,𝑢 1 } 𝑢 0 , 𝜇 0 𝑋 0 + 𝑢 1 , 𝜇 1 𝑋 1 + 𝜄 {𝑢 0 ⊕𝑢 1 -𝑐 ≤0} , 𝜇 0 ⊗ 𝜇 1 𝑋 0 ×𝑋 1 , (3.10)
where 𝜄 is the characteristic function. The integration of the constraint against the product measure 𝜇 0 ⊗ 𝜇 1 will be discussed for the primal regularized problem in Section 3.3 below. Entropic regularization replaces the hard constraint (2.2) with a smooth differentiable barrier function depending on a small parameter 𝜖:

OT 𝜖 (𝜇 0 , 𝜇 1 ) := sup {𝑢 𝜖 ,0 ,𝑢 𝜖 ,1 } 𝐽 𝜖 (𝑢 𝜖 ,0 , 𝑢 𝜖 ,1 ), (3.11) 𝐽 𝜖 (𝑢 𝜖 ,0 , 𝑢 𝜖 ,1 ) := 𝑢 𝜖 ,0 , 𝜇 0 𝑋 0 + 𝑢 𝜖 ,1 , 𝜇 1 𝑋 1 -𝜖 e (𝑢 𝜖 ,0 ⊕𝑢 𝜖 ,1 -𝑐)/𝜖 , 𝜇 0 ⊗ 𝜇 1 𝑋 0 ×𝑋 1 .
The penalization is strictly concave. As in the semi-dual approach, this is a large gain as this formulation is amenable to the solver arsenal of smooth convex/concave optimization. The Sinkhorn algorithm is the simplest iterative (in 𝑘) coordinate ascent method (we drop the . 𝜖 notation):

𝑢 𝑘 0 = arg sup 𝑢 0 𝑢 0 , 𝜇 0 𝑋 0 -𝜖 e (𝑢 0 ⊕𝑢 𝑘-1 1 -𝑐)/𝜖 , 𝜇 0 ⊗ 𝜇 1 𝑋 0 ×𝑋 1 , 𝑢 𝑘 1 = arg sup 𝑢 1 𝑢 1 , 𝜇 1 𝑋 1 -𝜖 e (𝑢 𝑘 0 ⊕𝑢 1 -𝑐)/𝜖 , 𝜇 0 ⊗ 𝜇 1 𝑋 0 ×𝑋 1 .
(3.12)

The maximization can be expressed in closed form as

𝑢 𝑘 0 = LSE 𝜖 𝜇 1 ,𝑋 1 (𝑢 𝑘-1 1 ), 𝑢 𝑘 1 = LSE 𝜖 𝜇 0 ,𝑋 0 (𝑢 𝑘 0 ), (3.13)
where we introduced log/sum/exp (LSE), a classical and optimized function in modern scientific software (e.g. PyTorch):

LSE 𝜖 𝜌,𝑋 (𝑢) := {𝑥 1 ↦ → -𝜖 log( e (𝑢-𝑐)/𝜖 , 𝜌 𝑋 )}, (3.14)
which sends, alternately, a function or vector defined on 𝑋 0 to one defined on 𝑋 1 , and reciprocally (remember 𝑐 is defined on 𝑋 0 × 𝑋 1 ). At convergence, (3.13) is an 𝜖-smooth version of (2.5) up to a constant shift of the potentials. The formulation (3.12) holds for both continuous and discrete measures. Using the discretization of Section 1.3, one iteration of (3.13) is quadratic in 𝑁. A nice feature of the algorithm is its formal independence on 𝑐, unlike the Monge-Ampère and semidiscrete approaches. Details and references can be found in the survey by [START_REF] Vialard | An elementary introduction to entropic regularization and proximal methods for numerical optimal transport[END_REF] and the recent papers by Di [START_REF] Marino | An optimal transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm[END_REF] and Peyré and Cuturi (2019, §4.2). It gives a linear convergence rate in

𝑂(1 -e 𝐿 𝑐 (Diam 𝜇 0 )/𝜖 ) 𝑂(1 -e 𝐿 𝑐 (Diam 𝜇 1 )/𝜖 ),
where 𝐿 𝑐 is the Lipschitz constant of 𝑐. For the Euclidean quadratic cost (1.1) this is the friendly constant 1, but, as we can see, the dependence on 𝜖 is not so nice.

We finish this section with several important remarks for the numerical analysts interested in using (3.13).

The kernel appearing in the iterations scales like e 𝐿 𝑐 (Diam 𝜇 0,1 )/𝜖 . In finite precision this leads to computer underflow and overflow when 𝜖 is too small or 𝐿 𝑐 Diam 𝜇 0,1 , the scale of the ground displacement, is too large (or both). Several stabilization techniques are available: see [START_REF] Schmitzer | Stabilized sparse scaling algorithms for entropy regularized transport problems[END_REF].

The operation cost of Sinkhorn iterations can be drastically reduced for the Euclidean quadratic cost (1.1). Rewrite e (𝑢-𝑐)/𝜖 , 𝜌 𝑋 in the LSE function as

{𝑦 ↦ → e 𝑢(𝑥)/𝜖 𝛾 𝜖 (𝑦 -𝑥)), 𝜌 𝑋 },
where

𝛾 𝜖 (𝑧) = 1 (2𝜋𝜖) 𝑑/2 e -𝑧 2 /(2𝜖 ) (3.15)
is a Gaussian with variance 𝜖. The normalization constant is not important at this stage (and we will assume in this section it is absorbed by the potentials (𝑢 0 , 𝑢 1 )) but will be useful in Section 6.1. If 𝑋 0 = 𝑋 1 are 𝑑-dimensional regular Cartesian grids,

𝑧 2 = 𝑧 1 2 + 𝑧 1 2 + • • • + 𝑧 𝑑 2
and the number of 𝑧 𝑖 points is 𝑁 1/𝑑 . The exponential matrix tensorizes and the operation cost of one iteration reduces to 𝑂(𝑁 1+1/𝑑 ).

When the problem admits a Monge solution, it is possible to take advantage of the sparsity of the entropic plan to obtain an 𝑂(𝑁 log(𝑁)) algorithm via a strategy that is multiscale in 𝜖; see [START_REF] Schmitzer | Stabilized sparse scaling algorithms for entropy regularized transport problems[END_REF].

Much more is to be found in Peyré and Cuturi (2019, §4).

Entropic bias

The OT 𝜖 solution is, however, biased by the penalization term. The entropic bias is better understood for the Fenchel-Rockafellar primal corresponding to (3.11):

OT 𝜖 (𝜇 0 , 𝜇 1 ) := inf

𝜋 𝜖 ∈Π(𝜇 0 , 𝜇 1 ) 𝑐, 𝜋 𝜖 𝑋 0 ×𝑋 1 + 𝜖 KL(𝜋 𝜖 |𝜇 0 ⊗ 𝜇 1 ), (3.16)
where KL(𝜈|𝜈 0 ) := log(𝜈/𝜈 0 )-1, 𝜈 + 1, 𝜈 0 is the Kullback-Leibler divergence.

The compact notations may be a bit confusing: 𝜈/𝜈 0 implies that 𝜇 is absolutely continuous with respect to the measure 𝜈 0 (continuous or discrete), otherwise KL takes, by convention, an infinite value. So supp 𝜈 ⊂ supp 𝜈 0 and 𝜈/𝜈 0 is a density. The penalization 𝜈 ↦ → KL(𝜈|𝜈 0 ) is strictly convex, positive and vanishes for 𝜈 = 𝜈 0 . The penalization 𝜖 KL(𝜋 𝜖 |𝜇 0 ⊗ 𝜇 1 ) biases the entropic plan 𝜋 * 𝜖 towards the measure 𝜇 0 ⊗ 𝜇 1 . Figure 3.2 illustrates this effect. The choice of this particular reference measure is guided by the following consideration: it is our best direct guess for 𝜋 * 𝜖 in Π(𝜇 0 , 𝜇 1 ) in terms of both support and value. We therefore expect to reduce the bias with the classical optimal transportation solution introduced by the regularization. Note that it is possible to choose a reference measure such that Π(𝜇 0 , 𝜇 1 ) will be empty of plans with the corresponding support. The problem then becomes infeasible.

A few lines of computations rewrite (3.16) as inf

𝜋 𝜖 ∈M(𝑋 0 ×𝑋 1 ) 𝜄 {𝑃 𝑋 0 # 𝜋 𝜖 =𝜇 0 } + 𝜄 {𝑃 𝑋 1 # 𝜋 𝜖 =𝜇 1 } + 𝜖 KL(𝜋 𝜖 |𝜋 0 𝜖 ), (3.17)
where M is the set of Radon measure and (𝑥 0 , 𝑥 1 ) ↦ → 𝜋 0 𝜖 (𝑥 0 , 𝑥 1 ) = 𝛾 𝜖 (𝑥 1 -𝑥 0 )𝜇 0 (𝑥 0 )𝜇 1 (𝑥 1 ).

(3.18)

It is important to note that this formulation is unconstrained. The (strict) positivity of 𝜋 𝜖 is enforced by the KL penalization and its mass sums to 1 because of the marginal constraints. Optimizers are necessarily in P(𝑋 0 × 𝑋 1 ). The equivalence between (3.10) and (3.17) is a direct application of Fenchel-Rockafellar convex duality; see Theorem 6.3 in [START_REF] Carlier | Classical and Modern Optimization[END_REF]. Specifically, let (E, E ) and (F, F ) be a pair of dual normed spaces. Let 𝐴 : E → F be a continuous linear operator and 𝐴 its adjoint. Let 𝐹, 𝐺 be proper, convex and lower semicontinuous functions defined on E and F respectively. If there exists 𝜎 ∈ Dom(𝐹) such that 𝐺 is continuous at 𝐴𝜎, then sup

𝜎 ∈𝐸 -𝐹(-𝜎) -𝐺(𝐴𝜎) = inf 𝑞 ∈𝐹 𝐹 ★ (𝐴 𝑞) + 𝐺 ★ (𝑞) (3.19)
and the infimum is attained. Moreover, if there exists a maximizer 𝜎 * ∈ 𝐸, then there exists 𝑞 * ∈ 𝐹 such that We now apply (3.19) with

𝐴 : 𝜎 := (𝑢 0 , 𝑢 1 ) ∈ C(𝑋 0 ) × C(𝑋 1 ) ↦ → 𝑢 0 + 𝑢 1 ∈ C(𝑋 0 × 𝑋 1 ), 𝐴 : 𝑞 := 𝜋 𝜖 ∈ M(𝑋 0 × 𝑋 1 ) ↦ → (𝑃 𝑋 0 #𝜋 𝜖 , 𝑃 𝑋 1 #𝜋 𝜖 ) ∈ M(𝑋 0 ) × M(𝑋 1 ), 𝐹 : (𝑢 0 , 𝑢 1 ) ∈ C(𝑋 0 ) × C(𝑋 1 ) ↦ → 𝑢 0 , 𝜇 0 𝑋 0 + 𝑢 1 , 𝜇 1 𝑋 1 , 𝐹 ★ : (𝜌 0 , 𝜌 1 ) ∈ M(𝑋 0 ) × M(𝑋 1 ) ↦ → 𝜄 {𝜌 0 =𝜇 0 } + 𝜄 {𝜌 1 =𝜇 1 } , 𝐺 : 𝑢 ∈ C(𝑋 0 × 𝑋 1 ) ↦ → {𝑞 → KL(𝑞|𝜋 0 𝜖 )} ★ (𝑢) = -𝜖 e (𝑢-𝑐)/𝜖 , 𝜇 0 ⊗ 𝜇 1 𝑋 0 ×𝑋 1 , 𝐺 ★ : 𝑞 ∈ M(𝑋 0 × 𝑋 1 ) ↦ → 𝐺 ★ (𝑞) = 𝜖 KL(𝑞|𝜋 0 𝜖 ). (3.21)
Remember that 𝐹 = 𝐹 ★★ and 𝐺 = 𝐺 ★★ are convex though not necessarily continuous. In particular, 𝐹 ★ above takes infinite values. Using this particular choice the primal-dual optimality conditions (3.20) take the form

𝜋 * 𝜖 = e (𝑢 * 𝜖 ,0 ⊕𝑢 * 𝜖 ,1 -𝑐)/𝜖 𝜇 0 ⊗ 𝜇 1 and (𝑃 𝑋 0 #𝜋 𝜖 , 𝑃 𝑋 1 #𝜋 𝜖 ) = (𝜇 0 , 𝜇 1 ). (3.22)
Under the Fenchel-Rockafellar formalism, Sinkhorn iterations are the coordinate ascents

𝑢 𝑘 0 = arg sup 𝑢 0 -𝐹(-(𝑢 0 , 𝑢 𝑘-1 1 )) -𝐺(𝐴 (𝑢 0 , 𝑢 𝑘-1 1 ), 𝑢 𝑘 1 = arg sup 𝑢 1 -𝐹(-(𝑢 𝑘 0 , 𝑢 1 )) -𝐺(𝐴 (𝑢 𝑘 0 , 𝑢 1 )).
(3.23)

Details and generalizations of this formalism are given in Chizat, Peyré, Schmitzer and Vialard (2018a).

Let us restrict again to the Euclidean quadratic cost (1.1) and the discrete case (Section 1.3) to discuss (3.22). The entropic optimal transportation matrix is

𝜋 * 𝜖 ,(𝑖, 𝑗) = (𝑤 0,𝑖 e 𝑢 * 0,𝑖 /𝜖 ) 𝛾 𝜖 (𝑥 1, 𝑗 -𝑥 0,𝑖 ) (𝑤 1,𝑖 e 𝑢 * 1,𝑖 /𝜖 ), (3.24)
where {𝑢 0,𝑖 } := {𝑢 𝜖 ,0 (𝑥 0,𝑖 )} and {𝑢 1, 𝑗 } := {𝑢 𝜖 ,1 (𝑥 1, 𝑗 )} are converged solutions of the (3.12) Sinkhorn iterations. The Gaussian matrix Γ 𝜖 := {𝛾 𝜖 (𝑥 1, 𝑗 -𝑥 0,𝑖 )} (see (3.15)) is strictly positive and decreases exponentially away from the diagonal (if the marginal discretizations are identical). On each side there is a diagonal matrix scaling Γ 𝜖 such that 𝜋 * 𝜖 ∈ Π(𝜇 0 , 𝜇 1 ). Denoting the discretization scale by ℎ, the effective bandwidth of Γ 𝜖 in finite precision admits a maximum admissible displacement. It will decrease as ℎ 2 /𝜖 increases. It may again cause the set of 𝜋 * 𝜖 ∈ Π(𝜇 0 , 𝜇 1 ) in the form (3.24) to be empty. Decreasing the entropic bias therefore involves ℎ and 𝜖; Berman (2020) gives quantitative results. Set ℎ = 1/𝑁 1/𝑑 = 𝜖 and assume the ground space discretization of the smooth densities 𝜇 0,1 is on the scale of 𝜖. After running (3.13) for 𝑘 𝜖 = 1/𝜖 log(1/𝜖) iterations, we get

𝑢 𝑘 𝜖 0, 𝜖 -𝑢 * 0 ∞ ≤ 𝐶 𝜖 log 1 𝜖 . (3.25)
As for the optimal plan, Berman (2020) gives the estimate

𝜋 𝑘 𝜖 𝜖 ≤ 𝐵 𝜖 𝜇 0 ⊗ 𝜇 1 , (3.26)
where 𝜖 , in theory, has full support on the grid, it will be negligible, and even null in finite precision, for an increasing number of points as 𝜖 decreases. A detailed numerical description of these phenomena and a proposed heuristic multi-scale method in 𝜖 can be found in [START_REF] Schmitzer | Stabilized sparse scaling algorithms for entropy regularized transport problems[END_REF]. A GPU implementation can also be found in [START_REF] Feydy | Analyse de données géométriques, au delà des convolutions[END_REF]. They both observe numerically the desired 𝑂(𝑁 log(𝑁)) complexity arising from the sparsity of the transport plan. A convergence proof of the multiscale method using (3.26) is given in [START_REF] Benamou | Capacity constrained entropic optimal transport, Sinkhorn saturated domain out-summation and vanishing temperature[END_REF]. We are back to the performance of the Monge-Ampère and semi-discrete formulations but only when smooth Monge maps exist.

(𝑥 0 , 𝑥 1 ) → 𝐵 𝜖 (𝑥 0 , 𝑥 1 ) := 𝑝 𝜖 𝑝 exp - 𝑐(𝑥 1 , 𝑇(𝑥 0 )) 𝜖 𝑝 . ( 3

Soft displacement interpolation

Computing the Wasserstein interpolant should in principle be as simple as solving the original optimal transportation problem and using the Monge map (3.9) to push forward 𝜇 0 to 𝜌 𝜏 . As discussed in the introduction, Monge maps are concepts that unfortunately may be lost after discretization and will definitely be lost after the entropic regularization. So we instead resort to the simplified Kantorovich formulation (3.7)-(3.8) and will recover 𝜌 𝜏 as the intermediate marginal. The discretization of (1.2) presented in Section 1.3 was straightforward. The discretization of (3.7)-(3.8) is more subtle. We have seen that the optimal plan 𝜋 * 3 charges triplets (𝑥 0 , 𝑥 0 + 𝜏(𝑥 1 -𝑥 0 ), 𝑥 1 ) when the original optimal transportation solution 𝜋 * (𝑥 0 , 𝑥 1 ) also does. If 𝑋 𝑘 = {𝑥 𝑘,𝑖 } 𝑖=1, 𝑁 for 𝑘 = 0, 1, we need to cover all possible displacement interpolants for the discretization of the support 𝜌 𝜏 . This translates into {𝑥 0,𝑖 + 𝜏(𝑥 1, 𝑗 -𝑥 0,𝑖 )} (𝑖, 𝑗)∈ [1, 𝑁 ] 2 with possibly 𝑁 2 points. The 3-marginal transport plan will be a 𝑁 × 𝑁 2 × 𝑁 tensor. Solving the discretized (3.7)-(3.8) linear program is not possible for reasonable discretization sizes. One may, of course, use an arbitrary intermediate grid; the discretization of all three ground spaces may be, for instance, identical regular Cartesian grids of size 𝑁. In so doing, we add a constraint on the support of 𝜌 𝜏 . This often generates non-smooth oscillating results even for smooth 𝜇 0 and 𝜇 1 .

Entropic penalization offers a win-win solution to this problem, in terms of both computational efficiency and removing the grid pollution explained above. The suggested primal entropic version of (3.7) is inf

𝜋 3, 𝜖 ∈Π 2 (𝜇 0 , 𝜇 1 ) 𝑐 3 , 𝜋 3, 𝜖 𝑋 0 ×𝑋 ×𝑋 1 + 𝜖 KL(𝜋 3, 𝜖 |𝜇 0 ⊗ 1 𝑋 ⊗ 𝜇 1 ).
(3.28)

As 𝜋 3 and 𝜋 share the same marginal constraints (3.8), the dual problem is shown to be sup (𝑢 𝜖 ,0 ,𝑢 𝜖 ,1 )

𝑢 𝜖 ,0 , 𝜇 0 𝑋 0 + 𝑢 𝜖 ,1 , 𝜇 1 𝑋 1 -𝜖 e (𝑢 𝜖 ,0 ⊕𝑢 𝜖 ,1 -𝑐 3 )/𝜖 , 𝜇 0 ⊗ 𝜇 1 𝑋 0 ×𝑋 ×𝑋 1 . (3.29)
It differs from (3.11) only in the final term. Using 𝑐 3 (see 3.6), the primal-dual optimality conditions (3.22) are

𝜋 * 3, 𝜖 (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) = e 𝑢 * 0 (𝑥 0 )/𝜖 𝛾 𝜖 /(1-𝜏) (𝑥 𝜏 -𝑥 0,𝑖 ) 𝛾 𝜖 /𝜏 (𝑥 1 -𝑥 𝜏 ) e 𝑢 * 1 (𝑥 1 )/𝜖 . (3.30)
Remember that the non-entropic

𝑃 𝑋 0 ×𝑋 1 #𝜋 * 3 = 1, 𝜋 3 𝑋 is the Kantorovich plan 𝜋 *
0→1 solution of (1.2). Applying the same marginal integration inside (3.30) amounts to the convolution

{𝑥 𝜏 ↦ → 𝛾 𝜖 /(1-𝜏) (𝑥 𝜏 -𝑥 0 ) 𝛾 𝜖 /𝜏 (𝑥 1 -𝑥 𝜏 )}, 1 𝑋 𝑋 = 𝛾 𝜖 /(𝜏(1-𝜏)) (𝑥 1 -𝑥 0 ), (3.31)
where we are taking some liberties with the constants and used the convolution of Gaussian properties. We get

𝑃 𝑋 0 ×𝑋 1 #𝜋 * 3, 𝜖 (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) = e 𝑢 * 0 (𝑥 0 )/𝜖 𝛾 𝜖 /(𝜏(1-𝜏)) (𝑥 1 -𝑥 0 ) e 𝑢 * 1 (𝑥 1 )/𝜖 , (3.32)
and recognize 𝜋 * 𝜖 /(𝜏(1-𝜏)) the solution of OT 𝜖 (𝜇 0 , 𝜇 1 ), 𝜖 = 𝜖/(𝜏(1 -𝜏)) because it satisfies (3.24) and also belongs to Π(𝜇 0 , 𝜇 1 ) by construction. After discretization, we use (3.13) with 𝜖 to compute an approximation to (3.24) and of the potentials 𝑢 * 0,1 shared by (3.32) and (3.30). Then we discretize the formula (3.30) to smoothly interpolate

𝜌 𝜏 = 𝑃 𝑋 #𝜋 * 3, 𝜖 = 1, 𝜋 * 3, 𝜖 𝑋 0 ×𝑋 1 = 𝑥 𝜏 ↦ → (3.33) ∑ 𝑖 (𝑤 0,𝑖 e 𝑢 * 0,𝑖 /𝜖 )𝛾 𝜖 /(1-𝜏) (𝑥 𝜏 -𝑥 0,𝑖 ) ∑ 𝑗 𝛾 𝜖 /𝜏 (𝑥 1, 𝑗 -𝑥 𝜏 )(𝑤 1,𝑖 e 𝑢 *
1, 𝑗 /𝜖 ) .

30 The convolution of two Gaussians is a Gaussian. The expectation is the sum of the expectations and the variance is the product of the variances. This is the product of two discrete convolutions with Gaussians. If the discretization ℎ is small compared to the standard deviations 𝜖/(1 -𝜏) and 𝜖/(1 -𝜏), then the interpolation will be smooth. Note that the 𝜖 parameter needs to be small to get a reasonable approximation of the non-entropic interpolant. This bounds the 𝜏 parameter away from 0 and 1. Generalizing the interpolation, the computation of means or barycentres of images is a crucial step in the statistical analysis of medical images, for example, or in texture synthesis [START_REF] Bonneel | Wasserstein barycentric coordinates: Histogram regression using optimal transport[END_REF]. Formally, at least, Wasserstein barycentres are a direct generalization of (3.2). Given a set of 𝑃 probability measures

{𝜇 𝑝 } 𝑝 ∈ [1, 𝑃] , 𝜌 𝜆 := arg inf 𝜌∈P(𝑋 ) 𝑃 ∑ 𝑝=1 𝜆 𝑝 W 2 2 (𝜌, 𝜇 𝑝 ), (3.34)
where

{𝜆 𝑝 } 𝑝 ∈ [1, 𝑃] ∈ (R + )
𝑃 are given and 𝑃 𝑝=1 𝜆 𝑝 = 1. The Kantorovich formulation of this problem was worked out by [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF], who showed the equivalence of different formulation. Replacing W 2 with OT 𝜖 in (3.34) leads to a Sinkhorn-like algorithm (Benamou et al. 2015); see Figure 3.3 for an illustration.

For a meaningful central limit theorem see Eichinger and Carlier (2021), which also uses an entropic regularization but this time of 𝜌.

Dynamic optimal transportation

O , 33

Velocity discretization and Lagrangian computational fluid dynamics formulation

The 3-marginal formulation developed in Section 3.1 can be generalized to 𝑀 marginals. Let the sets {𝑋 𝑚 } be 𝑀 identical copies of a compact set 𝑋 ∈ 𝑋, and consider the problem

𝜋 * 𝑀 ∈ arg min 𝜋 𝑀 ∈Π 𝑀 (𝜌 0 ,𝜌 1 ) 𝑐 𝑀 , 𝜋 𝑀 𝑀 𝑚=0 𝑋 𝑚 , (4.1)
where

𝑐 𝑀 (𝑥 0 , 𝑥 1 , . . . , 𝑥 𝑀 ) = 1 2 𝑀 𝑀 -1 ∑ 𝑚=0 𝑥 𝑚+1 -𝑥 𝑚 2 (4.2) and Π 𝑀 (𝜌 0 , 𝜌 1 ) := 𝜋 𝑀 ∈ P 𝑀 𝑚=0 𝑋 𝑚 : 𝑃 𝑋 0 #𝜋 𝑀 = 𝜌 0 , 𝑃 𝑋 𝑀 #𝜋 𝑀 = 𝜌 𝑀 . (4.3)
What we have discussed in the 3-marginal case remains true for 𝑀 marginals, that is, the mass optimally transported from the optimal 𝑥 * 0 to 𝑥 * 𝑀 travels in straight lines:

𝑥 * 𝑚 = 𝑥 * 0 + 𝑚 𝑀 (𝑥 * 𝑀 -𝑥 * 0 ), 𝑚 = 0, 𝑀 .
On such an optimal chain {𝑥 * 𝑚 } 𝑚=0•••𝑀 the mass transported is constant and equal to 𝜌 0 (𝑥 * 0 ). This is because 𝜋 * 0→𝑀 = 𝑃 𝑋 0 ×𝑋 𝑀 #𝜋 * 𝑀 is also the solution of (1.2). This may all seem pointless as the underlying transport can be obtained via a static 𝜋 0→𝑀 computation, but the following Lagrangian interpretation will lead us to important variations of optimal transportation, in particular the computational fluid dynamics (CFD) formulation [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF].

The optimal chain {𝑥 * 𝑚 } 𝑚=1•••𝑀 can be interpreted as a straight continuous curve

X * : (𝜏, 𝑥 * 0 ) ∈ [0, 1] × 𝑋 0 ↦ → X * 𝜏 (𝑥 * 0 )
initially at 𝑥 * 0 and passing through all points 𝑥 * 𝑚 = X * 𝜏 𝑚 (𝑥 * 0 ) at times 𝜏 𝑚 = (1/𝑀) 0, 𝑀 . We note that X * 𝜏 (𝑥 * 0 ) = 𝑀(𝑥 * 𝑚+1 -𝑥 * 𝑚 ) is a constant for all 𝜏 ∈]𝑚/𝑀, (𝑚 + 1)/𝑀 [ and for all 𝑚. The map 𝑥 * 0 ↦ → X * 𝜏 (𝑥 * 0 ) is exactly the optimal transportation map between 𝜌 0 and 𝜌 𝜏 obtained via displacement interpolation (3.9). Hence

{𝑥 * 0 ↦ → X * 𝜏 (𝑥 * 0 )}#𝜌 0 = 𝜌 𝜏 . (4.4)
We have already explained in the discrete case (Section 1.3) that optimal transport trajectories cannot cross except perhaps at final time 1. This is also the case in the continuous case. Therefore there exists a vector field (a 'velocity' as d𝜏 :

= 1/𝑀 is a time step) (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → 𝑉 * 𝜏 (𝑥) ∈ R 𝑑 such that X * 𝜏 (𝑥 * 0 ) = V * 𝜏 (X * 𝜏 (𝑥 * 0 )) = (1/𝑀)(𝑥 * 𝑚+1 -𝑥 * 𝑚 )
for all 𝜏 ∈ (𝑚/𝑀, (𝑚 + 1)/𝑀) and for all 𝑚.

Wrapping all the information above in the cost function (4.1), we have

𝑐 𝑀 , 𝜋 * 𝑀 𝑀 𝑚=0 𝑋 𝑚 = 𝑀 2 𝑀 -1 ∑ 𝑚=0 Id 𝑋 𝑚+1 -Id 𝑋 𝑚 2 , 𝜋 * 𝑚→𝑚+1 𝑋 𝑚 ×𝑋 𝑚+1 = 1 2 ∫ 1 0 V * 𝜏 (X * ) 2 d𝜏, 𝜌 0 𝑋 0 , (4.5)
where 𝜋 * 𝑚→𝑚+1 := 𝑃 𝑋 𝑚 ×𝑋 𝑚+1 #𝜋 * 𝑀 is the (2-marginal) 𝑚 → 𝑚 +1 integration of 𝜋 * 𝑀 . It is legitimate to tighten the optimal transportation optimization to the class of Lagrangian curves and minimize their kinetic energy

W 2 2 (𝜌 0 , 𝜌 1 ) = inf (X ,V)∈CL(𝜌 0 ,𝜌 1 )
1 2

∫ 1 0 V 𝜏 (X 𝜏 ) 2 d𝜏, 𝜌 0 𝑋 0 , (4.6)
where

CL(𝜌 0 , 𝜌 1 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → V 𝜏 (𝑥) ∈ R 𝑑 and (𝜏, 𝑥 0 ) ∈ [0, 1] × 𝑋 0 ↦ → X 𝜏 (𝑥 0 ) ∈ 𝑋 such that X 𝜏 = V 𝜏 (X 𝜏 ), X 0 = Id 𝑥 0 , X 1 #𝜌 0 = 𝜌 1 . (4.7)
This is the Lagrangian CFD formulation. What are the computational benefits or pitfalls of this dynamic optimal transportation? Let us go back to (4.1) and consider an 𝑁-point Cartesian grid discretization, denoted 𝑋 𝑁 , of the identical sets {𝑋 𝑚 }. The marginals 𝜌 0 and 𝜌 1 are empirical measures supported on the grid and 𝑀 𝑚=0 𝑋 𝑁 is the fully discretized space-time cylinder. Solving (4.1) then involves 𝜋 𝑀 , a probability measure over a discrete set of size 𝑁 𝑀 . This is the set of all possible chains 𝑥 0 → 𝑥 1 → • • • → 𝑥 𝑀 lying on the grid and 𝜋 𝑀 mass charges these chains under the initial and final constraint (4.3). This problem is usually too large to be tackled numerically and can in any case be simplified to the simpler interpolation problem described in Section 3.1 corresponding to 𝑀 = 2. Using the same grid of points, the Lagrangian paths can be represented as 𝑁 piecewise affine chains {𝑥 𝑚 𝑖 } 𝑚= 1, 𝑀 , 𝑖 = 1, 𝑁 . The size of the problem is now 𝑂(𝑁 𝑀), and (4.6) becomes

1 2 𝑁 ∑ 𝑖=0 d𝜏 𝑀 -1 ∑ 𝑚=0 𝑥 𝑚+1,𝑖 -𝑥 𝑚,𝑖 2 𝜌 0 (𝑥 0,𝑖 ). (4.8)
32 However, we will see in Section 6.1 that entropic regularization again offers a tractable numerical method to compute the geodesic curve 𝜏 → 𝜌 𝜏 . We will also show in Sections 5 and 6.2 that the introduction of time allows us to enrich the problem with additional constraints on the 𝑀 marginals.

The mass-preserving constraint becomes

𝜌 1 (𝑥) = ∑ 𝑖 ∈ [0, 𝑁 ], s.t. 𝑥=𝑥 𝑀 ,𝑖
𝜌 0 (𝑥 0,𝑖 ) for all 𝑥 ∈ 𝑋 𝑁 .

(4.9)

The Lagrangian CFD formulation is an intermediate approach between the Monge-Ampère and semi-discrete methods and the Kantorovich problem (4.1) -when a Monge solution exists, of course. Instead of considering any set of paths between time 0 and 1, an Eulerian vector field governs the trajectories on the space-time cylinder [0, 1] × 𝑋, and no crossings are allowed. The constraint (4.9), however, makes the problem non-convex. The next section shows how the Eulerian formulation fixes the convexity.

Eulerian computational fluid dynamics and hidden convexity

To do so, again use the discrete time representation (4.5):

𝑐 𝑀 , 𝜋 * 𝑀 𝑀 𝑚=0 𝑋 𝑚 = 𝑀 2 𝑀 -1 ∑ 𝑚=0 Id 𝑋 𝑚+1 -Id 𝑋 𝑚 2 , 𝜋 * 𝑚→𝑚+1 𝑋 𝑚 ×𝑋 𝑚+1 = 1 2 𝑀 -1 ∑ 𝑚=0 d𝜏 1 d𝜏 (X * 𝜏 𝑚+1 -X * 𝜏 𝑚+1 ) 2 , 𝜌 0 𝑋 0 (4.10) = 1 2 𝑀 -1 ∑ 𝑚=0 d𝜏 1 d𝜏 (X * 𝜏 𝑚+1 • (X * 𝜏 𝑚 ) -1 -Id 𝑋 𝑚 ) 2 , 𝜌 𝑚 𝑋 𝑚 .
The last line is obtained via the change of variable 𝑥 𝑚 = X * 𝑚 (𝑥 0 ) and using (4.4). We have played with the time step d𝜏 = 1/𝑀 to produce a velocity at (𝜏 𝑚 , 𝑥 𝑚 = X * 𝜏 𝑚 ), and for small d𝜏 (we omit • * for clarity),

V 𝜏 𝑚 (𝑥 𝑚 ) 1 d𝜏 (X 𝜏 𝑚+1 (X -1 𝜏 𝑚 (𝑥 𝑚 )) -𝑥 𝑚 ). (4.11)
Neglecting the discretization error (which is nil for straight optimal curves), we complete the Lagrangian-to-Eulerian-coordinates transformation in (4.10) and obtain

𝑐 𝑀 , 𝜋 * 𝑀 𝑀 𝑚=0 𝑋 𝑚 1 2 𝑀 -1 ∑ 𝑚=0 d𝜏 V 𝜏 𝑚 2 , 𝜌 𝑚 𝑋 𝑚 .
(4.12)

We now turn to the measure-preserving constraint (4.4). On all time intervals [𝜏 𝑚 , 𝜏 𝑚+1 ] we can use the geodesic property of the curve 𝜏 → 𝜌 * 𝜏 (Section 3.1). For all 𝑚 = 0, 𝑀 -1, 𝑇 𝑚 = X * 𝜏 𝑚+1 • (X * 𝜏 𝑚 ) -1 is the optimal map from 𝜌 * 𝜏 𝑚 to 𝜌 * 𝜏 𝑚+1 and in particular 𝑇 𝑚 #𝜌 * 𝜏 𝑚 = 𝜌 * 𝜏 𝑚+1 . Using the approximation (4.11) and the Jacobian equation formulation (2), the constraint becomes

𝜌 𝜏 𝑚 = det(𝐷 𝑥 𝑇 𝑚 ) (𝜌 𝑚+1 • 𝑇 𝑚 ) det(1 𝑋 𝑚 + d𝜏 𝐷 𝑥 V 𝜏 𝑚 ) (𝜌 𝑚+1 • (Id 𝑋 𝑚 + d𝜏 V 𝜏 𝑚 )) (1 𝑋 𝑚 + d𝜏 Tr(𝐷 𝑥 V 𝜏 𝑚 )) (𝜌 𝑚+1 + d𝜏 𝐷 𝑥 𝜌 𝑚+1 • V 𝜏 𝑚 ). (4.13)
To first order in d𝜏, we end up with

1 d𝜏 (𝜌 𝑚+1 -𝜌 𝜏 𝑚 ) -𝜌 𝑚+1 • 𝐷V 𝜏 𝑚 -𝐷 𝑥 𝜌 𝑚+1 • V 𝜏 𝑚 = -div 𝑥 (𝜌 𝑚+1 V 𝜏 𝑚 ). (4.14)
The Eulerian CFD formulation is obtained by letting d𝜏 → 0 in (4.12) and (4.14):

W 2 2 (𝜌 0 , 𝜌 1 ) = inf (𝜌,V)∈CE(𝜌 0 ,𝜌 1 ) ∫ 1 0 1 2 V 𝜏 (𝑥) 2 , 𝜌 𝜏 𝑋 d𝜏, (4.15)
where

CE(𝜌 0 , 𝜌 1 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → (𝜌 𝜏 (𝑥), V 𝜏 (𝑥)) ∈ R + × R 𝑑 such that 𝜕 𝜏 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 V 𝜏 ) = 0 and 𝜌 𝜏=0,1 = 𝜇 0,1 . (4.16)
The partial differential equation 𝜕 𝜏 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 V 𝜏 ) = 0 is called the continuity equation in computational fluid dynamics. It must be understood as the measure-preserving constraint (4.4). Rigorous proofs of the Lagrangian/Eulerian formulation require careful mathematical analysis and can be found in [START_REF] Brenier | Examples of hidden convexity in nonlinear PDEs[END_REF] or Santambrogio (2015, §6.1).

The Eulerian formulation 'breaks' the Lagrangian path dynamics and replaces it with the continuity equation. The nonlinearity in the initial/final marginals (4.9) disappears. The continuity equation becomes linear under the simple change of variable, from velocity to momentum: (𝜌, V) → (𝜌, 𝑚 := 𝜌 V). 'Hidden convexity' is revealed when plugging it into the integrand of the cost function.

Indeed (𝜌, V) → 1 2 𝜌 V 2 becomes (𝜌, 𝑚) → 𝐽(𝜌, 𝑚) := 1 2 𝜌 𝑚 2 . (4.17)
Convex analysis tells us that, as a homogeneous function of degree one, this is a convex function; more precisely it is the Legendre-Fenchel dual of the characteristic function of a convex set. The change of variable holds at first glance only for 𝜌 > 0 but can be made legitimate by replacing 𝐽 with 𝐽 ★★ , its bi-dual Legendre-Fenchel transform (see footnote 13) extended to R × R 𝑑 . Keeping the notation 𝐽, we find that for all (𝑎, 𝐵) ∈ R × R 𝑑 , 𝐽 ★ (𝑎, 𝐵) := 𝜒 𝑎+ 1 2 𝐵 2 ≤0 , (4.18)

33 Tr is the trace of a square matrix, and we have used the classical identity det(Id + d𝜏 𝐴) 1 + d𝜏 Tr(𝐴) for small d𝜏.

and redefining 𝐽,

(𝐽 ★★ ) = 𝐽(𝜌, 𝑚) :=            1 2 𝜌 𝑚 2 if 𝜌 > 0 0 if 𝜌 = 0
𝐴 : 𝜎 := (𝜙, 𝜙 0 , 𝜙 1 ) ∈ C([0, 1], 𝑋) × C(𝑋 0 ) × C(𝑋 1 ) ↦ → (𝜕 𝜏 𝜙, 𝐷 𝑥 𝜙), 𝐴 : 𝑞 := (𝜌, 𝑚) ∈ (M([0, 1] × 𝑋)) 𝑑+1 ↦ → -(𝜕 𝜏 𝜌 𝜏 + div 𝑥 𝑚, 𝜌 0 , -𝜌 1 ), 𝐹 : (𝜙, 𝜙 0 , 𝜙 1 ) ↦ → 𝜙 0 , 𝜇 0 𝑋 0 + 𝜙 1 , 𝜇 1 𝑋 1 , 𝐹 ★ : (𝑎, 𝜌 0 , 𝜌 1 ) ↦ → 𝜒 𝑎=0,𝜌 0 =𝜇 0 ,𝜌 1 =𝜇 1 , 𝐺 : (𝑎, 𝑏) ↦ → 1 𝑎+ 1 2 𝑏 2 ≤0 , 𝐺 ★ : 𝑞 ↦ → ∫ 𝑋 ∫ 1 0 𝐽(𝑞) d𝑥 d𝜏 (4.20)
(we assume above that 𝜇 0,1 has compact support, to avoid discussion of boundary conditions (4.16)). With this particular choice, the primal-dual optimality conditions (3.20) take the form

𝜕 𝜏 𝜌 * 𝜏 + div 𝑥 (𝜌 * 𝜏 𝐷 𝑥 𝜙 * ) = 0, 𝜌 * 𝜏=0,1 = 𝜇 0,1 , 𝜕 𝜏 𝜙 * 𝜏 + 1 2 𝐷 𝑥 𝜙 * 𝜏 2 = 0, 𝜌 * 𝜏 a.e. on [0, 1] × 𝑋.
(4.21)

We have taken some liberties with the definition of the spaces above, but the general idea that this is a duality result between continuous functions and measures is valid.

The system (4.21) must be understood in a weak sense. For a rigorous treatment see Santambrogio (2015, §6.1).

The concave dual problem sup 𝜎 ∈𝐸 -𝐹(-𝜎) -𝐺(𝐴𝜎) is

(𝜙 * 0 , 𝜙 * 1 ) := arg inf { 𝜙 : 𝜕 𝜏 𝜙 𝜏 + 𝐷 𝑥 𝜙 𝜏 2 /2≤0} 𝜙 0 , 𝜇 0 𝑋 0 + 𝜙 1 , 𝜇 1 𝑋 1 . (4.22)
It is worth mentioning that the semi-dual formulation (e.g. Section 2.1 and equation (2.14)) can be recovered by substituting the explicit Lax-Oleinik formula solution of the Hamilton-Jacobi equation (4.21) 𝜙 1 (𝑥 1 ) = inf 𝑥 𝜙 0 (𝑥)+ 1 2 𝑥 1 -𝑥 2 in (4.22). The Lax-Oleinik formula also holds backwards, and setting

(𝑢 0 , 𝑢 1 ) = 1 2 Id 𝑋 0 2 -𝜙 0 , 1 2 Id 𝑋 1 2 -𝜙 1 (4.23) 34 Use inf(•) = -sup -(•).
35 A classical result in the calculus of variations: 𝜙 𝑡 (𝑥) = inf 𝑥 0 𝜙 0 (𝑥 0 ) + 𝑡 𝐻 ★ (𝑡 -1 (𝑥 -𝑥 0 )) is the solution of the Hamilton-Jacobi equation 𝜕 𝜏 𝜙 𝜏 + 𝐻(𝐷 𝑥 𝜙 𝜏 ) = 0 if 𝐻 is convex and superlinear.

we recover the maximization problem (2.1) and 𝑐-concave duality (2.5).

The non-smooth convex optimization problem (4.15) is amenable to first-order optimization methods and in particular a technique called proximal splitting. A notable advantage of these methods is to preserve the positivity of the density 𝜌 during the optimization and therefore its stability; this is detailed in Benamou and Carlier (2015). The discretization on regular grids is discussed in Papadakis, [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF]. For rigorous Galerkin discretizations in time and space of the CFD problem, see the recent papers by [START_REF] Lavenant | Unconditional convergence for discretizations of dynamical optimal transport[END_REF], [START_REF] Natale | A mixed finite element discretization of dynamical optimal transport[END_REF] and their references. See also [START_REF] Hug | On the convergence of augmented Lagrangian method for optimal transport between nonnegative densities[END_REF], [START_REF] Guittet | On the time-continuous mass transport problem and its approximation by augmented Lagrangian techniques[END_REF] and [START_REF] Andreev | Preconditioning the augmented Lagrangian method for instationary mean field games with diffusion[END_REF].

The recent book by Carlier (2021, §7.4) gives a comprehensive review of these methods and their convergence. The method advocated by [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] is a Douglas-Rachford solver called ALG2 (in [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems[END_REF] and belongs to this family of proximal splitting methods. Even though the convergence of such an optimization algorithm is slow, the CFD approach remained state-of-the-art for a decade, the hidden convexity and the reduced size of the problem 𝑀 𝑁 instead of 𝑁 2 giving an advantage over linear programming. At this time, I am not aware of any convincing way to apply second-order optimization methods such as the damped Newton's algorithm used for the Monge-Ampère and semi-discrete methods. The difficulty seems to be connected with the ability of the CFD formulation to handle locally vanishing Wasserstein geodesics.

Other instances of hidden convexity are given in [START_REF] Brenier | Examples of hidden convexity in nonlinear PDEs[END_REF] and in particular the theory of generalized incompressible flows [START_REF] Brenier | The least action principle and the related concept of generalized flows for incompressible perfect fluids[END_REF]), a precursor to modern optimal transportation and multi-marginal optimal transportation. See the next section.

Variational formulations for Euler equations

Euler geodesics

The general setting in this section, in particular the space-time cylinder, is identical to Section 4 above. A flow is a (possibly infinite) collection of Lagrangian particles {X 𝜏 (𝑥 0 )} for (𝜏, 𝑥 0 ) ∈ [0, 1] × 𝑋 0 . It is incompressible if the 'material density is constant within a fluid parcel'. This is mathematically translated into the Eulerian property 𝜌 𝜏 = 1 𝑋 ; for all 𝜏 ∈ [0, 1], 𝜌 𝜏 is the mean field density of particles at time 𝜏 when their number tends to ∞. The mass per volume is uniformly constant. As in Section 4.1, a curve of densities evolves in time via push-forwards by the Lagrangian map 𝑥 0 ↦ → X 𝜏 (𝑥 0 ). The Lagrangian flow interpretation of incompressibility can be written as

X 𝜏 #1 𝑋 0 = 1 𝑋 𝜏 for all 𝜏 ∈ [0, 1]. (5.1)
The density is constant, but an incompressible fluid does not necessarily remain still. 'Lagrangian stillness' would be X 𝜏 (𝑥 0 ) = 𝑥 0 for all (𝜏, 𝑥) ∈ [0, 1] × 𝑋 0 or V 𝜏 (𝑥) = 0 for all 𝜏 ∈ [0, 1] × 𝑋. For example, let us consider a fluid with non-zero initial velocity V 0 (𝑥). If we let the particles evolve in straight lines 𝑥 0 + 𝜏 V 0 (𝑥 0 ) then they may move closer to or further away from each other, hence violating the incompressibility, or exit the domain, or both. To maintain incompressibility and remain in the domain, there needs to be a change in velocity: an acceleration or deceleration. Under this condition a moving incompressible flow may still end up with a different final particle configuration at time 1, X 1 : 𝑥 0 → X 𝑓 (𝑥 0 ) ≠ 𝑥 0 , but the particle trajectories are not straight. The acceleration is provided by a mean pressure field generated by the particles themselves. Under the additional assumption that the flow is a diffeomorphism for all 𝜏, the Jacobian equation interpretation (see footnote 2) of incompressibility is det(𝐷 𝑥 0 X 𝜏 ) = 1. The set of such Lebesgue-preserving diffeomorphisms in 𝑋 is denoted SDiff(𝑋). As in the optimal transportation case (Section 4.1), the Eulerian vector field V 𝜏 can be defined and the familiar div 𝑥 V 𝜏 = 0 incompressibility condition can be derived from (4.14). This tells us, at the pointwise infinitesimal level, that dilations or compressions must be balanced coordinate-wise. Two centuries after Leonhard Euler, [START_REF] Arnold | Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF] gave a geodesic interpretation of the Euler equations. Replacing V 𝜏 with X 𝜏 , taking into account the final configuration constraint (X 𝑓 is given) and also the incompressibility, modify (4.6)-(4.7) to obtain

EG(X 𝑓 ) := inf (X 𝜏 ,V 𝜏 )∈𝐶 EG (X 𝑓 ) 1 2 ∫ 1 0 X 𝜏 2 d𝜏, 1 𝑋 0 𝑋 0 , (5.2)
where

𝐶 EG (X 𝑓 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → V 𝜏 (𝑥) ∈ R 𝑑 and (𝜏, 𝑥 0 ) ∈ [0, 1] × 𝑋 0 ↦ → X 𝜏 (𝑥 0 ) ∈ 𝑋 such that X 𝜏 = V 𝜏 (X 𝜏 ), X 0 = Id 𝑋 0 , X 1 = X 𝑓 , X 𝜏 ∈ SDiff(𝑋) for all 𝜏 ∈ [0, 1] . (5.3)
The cost EG(𝑋 𝑓 ) minimizes, with respect to the L 2 (𝑋, R 𝑑 ) metric, the length of the flow between Id 𝑋 0 and X 𝑓 , constrained onto the set of Lebesgue-preserving diffeomorphisms SDiff(𝑋). This is the definition of a geodesic for the optimal flow 𝜏 → X * 𝜏 . The particle acceleration must be tangent to SDiff(𝑋) and therefore the gradient of a 'pressure' potential 𝑝 * :

X * 𝜏 = 𝐷 𝑥 𝑝 * 𝜏 (X * 𝜏 ) for all 𝜏 ∈ [0, 1].
(5.4)

Together with constraints (5.3), this is the Lagrangian formulation of the Euler equation with initial and final time boundary conditions on the particle configuration. The equivalence with the Euler partial differential equations in (V 𝜏 , 𝑝 𝜏 )

𝜕 𝜏 V 𝜏 + V 𝜏 • 𝐷 𝑥 (V 𝜏 ) = 𝐷 𝑥 𝑝 𝜏 , div 𝑥 V 𝜏 = 0 (5.5)
is well known and can also be recovered (as in Section 4.2) using Fenchel-Rockafellar optimality conditions. This pure Eulerian formulation, however, cannot model the Lagrangian initial-final particle boundary conditions. The existence of global-in-time solutions to Euler equations remains a challenge. Arnold's approach offers a way to characterize solutions with the variational formulation (5.2)-( 5.3). However, SDiff is not closed in L 2 (𝑋, R 𝑑 ) and the existence of an Eulerian vector field is conditional on the flow remaining in SDiff. As such, (5.2)-( 5.3) can only be used with additional regularity assumptions on X 𝑓 ; see [START_REF] Ebin | Groups of diffeomorphisms and the motion of an incompressible fluid[END_REF].

Generalized geodesics

The idea in [START_REF] Brenier | The least action principle and the related concept of generalized flows for incompressible perfect fluids[END_REF] is to relax the constraint X 𝜏 ∈ SDiff(𝑋) in (5.3). Following (5.1), replace SDiff with S = {X ∈ L 2 (𝑋, 𝑋), X #1 𝑋 = 1 𝑋 }, i.e. still measure-preserving maps but not necessarily diffeomorphisms. Particles may cross, and the existence of an Eulerian velocity V is not guaranteed. The constraint X 𝜏 = V 𝜏 is relaxed and we optimize on X 𝜏 alone: (5.6) where

GEG(X 𝑓 ) := inf X 𝜏 ∈𝐶 GEG (X 𝑓 ) 1 2 ∫ 1 0 X 𝜏 2 d𝜏, 1 𝑋 0 ,
𝐶 GEG (X 𝑓 ) := (𝜏, 𝑥 0 ) ∈ [0, 1] × 𝑋 0 ↦ → X 𝜏 (𝑥 0 ) ∈ 𝑋 such that X 0 = Id 𝑋 0 , X 1 = X 𝑓 , X 𝜏 ∈ S(𝑋) for all 𝜏 ∈ [0, 1] . (5.7)
The set S is closed and is even the completion of SDiff for 𝑑 < 3. A numerical method based on a discretization of (5.3)-(5.7) has been proposed by [START_REF] Mérigot | Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport[END_REF]. It uses a finite number of particles following piecewise affine paths {𝑥 𝑚 𝑖 } 𝑚=0,..., 𝑀 , 𝑖 = 1, 𝑁 , as in the Lagrangian CFD optimal transportation formulation (4.8). Discrete incompressibility requires that particles, while moving, remain at every time 𝜏 𝑚 equi-distributed in space with respect to Lebesgue measure.

There are many such configurations corresponding to local minima of the nonconvex semi-discrete optimal transportation functional:

{𝑥 𝑚 𝑖 } 𝑖=1,..., 𝑁 → W 2 2 ∑ 𝑖=1,..., 𝑁
𝛿 𝑥 𝑚 𝑖 , 1 𝑋 .

(5.8) [START_REF] Mérigot | Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport[END_REF] propose penalizing (5.6) with (5.8) to approximately maintain the incompressibility instead of using the hard constraint in (5.7). The final configuration is also enforced through a penalization. The now unconstrained functional to minimize is inf

{𝑥 𝑚 𝑖 } 𝑚=0,..., 𝑀 ,𝑖=1,..., 𝑁 1 2 𝑁 ∑ 𝑖=0 d𝜏 𝑀 -1 ∑ 𝑚=0 𝑥 𝑚+1,𝑖 -𝑥 𝑚,𝑖 2 + 𝜆 1 𝑁 ∑ 𝑖=0 𝑥 𝑀 ,𝑖 -X 𝑓 (𝑥 0,𝑖 ) 2 + 𝜆 2 𝑀 -1 ∑ 𝑚=0 W 2 2 ∑ 𝑖=1,..., 𝑁
𝛿 𝑥 𝑚 𝑖 , 1 𝑋 .

(5.9)

A careful scaling of the parameters (𝑀, 𝑁, 𝜆 1 , 𝜆 2 ) is needed to ensure convergence to geodesics (5.2)-( 5.3) when they exist. The same technique has been adapted to the Euler-Dirichlet problem (Gallouët and Mérigot 2018).

The incompressible Beltrami flow on the unit square X = [0, 1] 2 is an analytical non-trivial global-in-time solution of Euler equations. It corresponds to the stationary velocity and pressure

V * (𝑥) = (-cos(𝜋 𝑎) sin(𝜋 𝑏), sin(𝜋 𝑎) cos(𝜋 𝑏)), 𝑝 * (𝑥) = 1 2 (sin(𝜋 𝑎) 2 + sin(𝜋 𝑏) 2 ), (5.10)
where 𝑥 = (𝑎, 𝑏) are the Cartesian coordinates. It can be used as a test case for Euler geodesics. Figure 5.1 is reproduced from [START_REF] Mérigot | Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport[END_REF]. Particles are labelled with just three different colours to provide for a general description of the flow. The first line shows the precomputed configuration X 𝑓 for different final times; this is the classical solution of X * 𝜏 = V * . The following lines show the Euler geodesic solution of (5.9). For small times and as predicted by the theory, the solution remains classical. For larger times a shorter geodesic is found in the larger set S by allowing particles close to the centre of rotation to cross and travel directly to the final configuration. The observed 'mixing' suggests that the solution is dependent on the discretization and points to a looser relaxation (see the next subsection).

Generalized incompressible flows

So far we have only dropped the diffeomorphism constraint. To go further, and following the Kantorovich paradigm, we now allow the unit of mass at all 𝑥 0 to split and be carried by more than one path. Let us consider the set Ω = C 0 ([0, 1], 𝑋) of all possible continuous paths in the space-time cylinder. The mass carried by every path 𝜏 ∈ [0, 1] → 𝜔(𝜏) will be modelled as 𝜋(𝜔), where 𝜋 ∈ P(Ω) now denotes probabilities on the space of curves Ω. The Eulerian density of paths 𝜌 𝜏 ∈ P(𝑋) (5.11)

The transport plan between times 𝜏 0 and 𝜏 1 is given by 𝜋 𝜏 0 →𝜏 1 := (𝑒 𝜏 0 , 𝑒 𝜏 1 )#𝜋.

(5.12)

If 𝜋 0→1 (𝑥 0 , X 𝑓 (𝑥 0 )) = 1 𝑋 0 , the unit of mass initially at 𝑥 0 ends up finally at X 𝑓 (𝑥 0 ) and only there. The final configuration Lebesgue-preserving map X 𝑓 : 𝑋 0 ↦ → 𝑋 1 therefore corresponds to the plan

𝜋 0→1 = (Id 𝑋 0 , X 𝑓 )#1 𝑋 0 .
(5.13)

A generalized incompressible flow will be the solution of

GIF(X 𝑓 ) := inf 𝜋 ∈𝐶 GIF (X 𝑓 ) 1 2 ∫ 1 0 𝜔(𝜏) 2 d𝜏, 𝜋 Ω , (5.14)
where

𝐶 GIF (X 𝑓 ) := 𝜋 ∈ P(Ω) : (𝑒 0 , 𝑒 1 )#𝜋 = (Id 𝑋 0 , X 𝑓 )#1 𝑋 0 𝑒 𝜏 #𝜋 = 1 𝑋 𝜏 for all 𝜏 ∈ [0, 1] .
(5.15)

In line with the Kantorovich formulation of optimal transportation, we are back to a 'simple' Eulerian linear programming problem but set on a 'large' infinitedimensional space P(Ω). The last line in (5.15) is Eulerian incompressibility and the middle line the initial/final boundary conditions. Note that even though (𝑒 0 , 𝑒 1 )#𝜋 is fully described by the map X 𝑓 , it does not imply that the generalized flow itself concentrates on a diffeomorphism in SDiff for all time. Mass can divide and flow from 𝑥 0 to X 𝑓 (𝑥 0 ) along an arbitrary number of paths, possibly crossing with different velocities. We follow the discretization in Section 4.1: let Ω 𝑁 , 𝑀 = 𝑀 𝑚=0 𝑋 𝑁 be the discrete (in space and time) cylinder, where 𝑋 𝑁 is a fixed 𝑁-point grid and {𝜏 𝑚 }, 𝑚 = 1, 𝑀 is a time discretization. The discretization of 𝜋 is a multimarginal plan 𝜋 𝑀 , 𝑁 ∈ P(Ω 𝑁 , 𝑀 ). 37 By definition (see Section 1) and for any measurable subset 𝐴 ⊆ 𝑋, 𝑒 𝜏 #𝜋(𝐴) = 𝜋({𝜔 ∈ Ω, 𝜔(𝜏) ∈

𝐴}).

38 For any pair of measurable subsets (𝐴 𝜏 0 , 𝐴 𝜏 1 ) ∈ 𝑋 𝜏 0 × 𝑋 𝜏 1 , we have (𝑒 𝜏 0 , 𝑒 𝜏 1 )#𝜋(𝐴 𝜏 0 , 𝐴 𝜏 1 ) = d𝜋({𝜔 ∈ Ω, (𝜔(𝜏 0 ), 𝜔(𝜏 1 )) ∈ 𝐴 𝜏 0 × 𝐴 𝜏 1 }). 39 Remember that 𝜌 1 = 𝑃 𝑋 1 #𝜋 0→1 is a probability measure and therefore necessarily 𝜋 0→1 (𝑥 0 , 𝑥 1 ) = 0 for all 𝑥 1 ≠ X 𝑓 (𝑥 0 ). We have already discussed (see (2.11)) how a transport map can be encoded into a plan.

Resorting once more to a piecewise affine approximation of paths (see (4.1)-(4.2)) and using an 𝑁-point grid 𝑋 𝑁 for the space discretization, we end up with inf (5.16) where Π 𝑀 , 𝑁 (X 𝑓 ) := 𝜋 𝑀 , 𝑁 ∈ P(Ω 𝑀 , 𝑁 ) :

𝜋 𝑀 , 𝑁 ∈Π 𝑀 , 𝑁 (X 𝑓 ) 𝑐 𝑀 , 𝜋 𝑀 , 𝑁 𝑀 𝑚=0 𝑋 𝑁 ,
𝑃 𝑋 0 ×𝑋 1 #𝜋 𝑀 , 𝑁 = (Id 𝑋 0 , X 𝑓 )#1 𝑋 𝑁 ,
𝑃 𝑋 𝑚 #𝜋 𝑀 , 𝑁 = 1 𝑋 𝑁 for all 𝑚 = 0 . . . ,𝑀 .

(5.17)

This problem minimizes the same action as the CFD optimal transportation problem (see (4.1)), but the incompressibility constraint is active on all marginals. We will see another example of multi-marginal constraints in Section 6.2. It is convenient to eliminate the 2-marginal 'plan' constraint (second line in (5.17)). It can be relaxed by adding a penalization to the displacement cost (𝜆 > 0 and 'large'):

𝑐 𝑀 ,X 𝑓 (𝑥 0 , , 𝑥 1 , . . . , 𝑥 𝑀 ) = 1 2 𝑀 𝑀 -1 ∑ 𝑚=0 𝑥 𝑚+1 -𝑥 𝑚 2 + 𝜆 𝑥 𝑀 -X 𝑓 (𝑥 0 ) 2 . (5.18)
The problem becomes inf (5.20) no longer depends on X 𝑓 . From the numerical point of view, this problem suffers from the same difficulties as the discrete Kantorovich problem (1.4) on an even larger scale (Section 1.3).

Entropic regularization for generalized incompressible flows

Following Benamou, Carlier and Nenna (2019a)

𝐴 : 𝜎 := (𝑢 0 , 𝑢 1 , . . . ,𝑢 𝑀 ) ∈ (R 𝑁 ) 𝑀 ↦ → 𝑀 ∑ 𝑚=0 𝑢 𝑚 ∈ R 𝑀 ×𝑁 , 𝐴 : 𝑞 := 𝜋 𝜖 ∈ R 𝑁 ×𝑀 ↦ → (𝑃 𝑋 0 #𝜋 𝜖 , 𝑃 𝑋 1 #𝜋 𝜖 , . . . , 𝑃 𝑋 𝑀 #𝜋 𝜖 ) ∈ (R 𝑁 ) 𝑀 , 𝐹 : (𝑢 0 , 𝑢 1 , . . . ,𝑢 𝑀 ) ∈ (R 𝑁 ) 𝑀 ↦ → 𝑀 ∑ 𝑚=0 𝑢 𝑚 , 1 𝑋 𝑁 , 𝐹 ★ : (𝜌 0 , 𝜌 1 , . . . , 𝜌 𝑀 ) ∈ (R 𝑁 ) 𝑀 ↦ → 𝑀 ∑ 𝑚=0 𝜒 𝜌 𝑚 =1 𝑋 𝑁 , 𝐺 : 𝑢 ∈ R 𝑀 ×𝑁 ↦ → 𝐺(𝑢) = 𝜖 e 𝑢/𝜖 , 𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 𝑀 𝑚=0 𝑋 𝑚 , 𝐺 ★ : 𝑞 ∈ R 𝑁 ×𝑀 ↦ → 𝐺 ★ (𝑞) = 𝜖 KL(𝜋 𝜖 |𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 ).
(5.24)

The primal-dual optimality conditions (3.20) take the form

𝜋 * 𝜖 = e ( 𝑀 𝑚=0 𝑢 * 𝑚 )/𝜖 𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 , 𝑃 𝑋 𝑚 #𝜋 * 𝜖 = 1 𝑋 for all 𝑚.
(5.25)

Likewise the coordinate ascent Sinkhorn algorithm (3.23) generalizes to

𝑢 𝑘 𝑚 = arg sup 𝑢 𝑚 -𝐹(-{𝑢 𝑘 0 , . . . , 𝑢 𝑘 𝑚-1 , 𝑢 𝑚 , 𝑢 𝑘-1 𝑚 , . . . , 𝑢 𝑘-1 𝑀 }) -𝐺(𝐴{𝑢 𝑘 0 , . . . , 𝑢 𝑘 𝑚-1 , 𝑢 𝑚 , 𝑢 𝑘-1 𝑚 , . . . , 𝑢 𝑘-1 𝑀 }) (5.26)
for all 𝑚 = 0, . . . , 𝑀 at each 𝑘 iteration. As for the 2-marginal Sinkhorn, these are strictly concave unconstrained problems, leading to the nonlinear set of equations:

𝜕 𝑢 𝑚 𝐹(-{𝑢 𝑘 0 , . . . , 𝑢 𝑘 𝑚-1 , 𝑢 𝑚 , 𝑢 𝑘-1 𝑚 , . . . , 𝑢 𝑘-1 𝑀 }) = 𝜕 𝑢 𝑚 𝐺(𝐴 {𝑢 𝑘 0 , . . . , 𝑢 𝑘 𝑚-1 , 𝑢 𝑚 , 𝑢 𝑘-1 𝑚 , . . . , 𝑢 𝑘-1 𝑀 }).
(5.27)

The left-hand side 𝜕 𝑢 𝑚 𝐹 = 1 𝑋 𝑁 is given. The right-hand side simplifies to

𝑢 𝑚 → 𝜕 𝑢 𝑚 𝐺(𝐴 {𝑢 0 , . . . , 𝑢 𝑚-1 , 𝑢 𝑚 , 𝑢 𝑚 , . . . , 𝑢 𝑀 }) = e ( 𝑀 𝑚 =0 𝑢 𝑚 )/𝜖 , 𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 𝑋 0 ו••×𝑋 𝑚-1 ×𝑋 𝑚+1 ו••𝑋 𝑀 = e 𝑢 𝑚 /𝜖 e ( 𝑚 ≠𝑚 𝑢 𝑚 )/𝜖 , 𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 𝑋 0 ו••×𝑋 𝑚-1 ×𝑋 𝑚+1 ו••×𝑋 𝑀 .
(5.28)

As in the classic Sinkhorn algorithm, the equation (5.27) is explicit in 𝑢 𝑚 . The numerical limits of this algorithm ultimately rest on the computational cost of the sum in (5.28). The product structure of (5.23) involves only two successive (in a circular way) marginal kernels, and the chain of summations in the expression above is 'broken' at

𝑚: 𝑚 + 1 → • • • → 𝑥 𝑀 → 𝑥 0 → • • • → 𝑥 𝑚-1 . Computing
(5.28) consists in 𝑀 independent matrix vector products

∑ 𝑗 𝛾 𝜖 /d𝜏 (𝑥 𝑝, 𝑗 -𝑥 𝑝+1,𝑖 ) e 𝑢 𝑝, 𝑗 /𝜖 ,
where 𝑝 runs over the chain and the indices (𝑝, 𝑗) and (𝑝 + 1, 𝑖) over the space discretization at times 𝑝 and 𝑝 + 1. We finally find that computing (5.28) for all {𝑚} requires 𝑂(𝑀 2 𝑁 2 ) operations which can be reduced to 𝑂(𝑀 2 𝑁 1+1/𝑑 ) using the separability of Euclidean quadratic cost (1.1). See Benamou et al. (2019a) for more on the convergence of the entropic regularization as 𝜖 → 0. For the convergence rate of multi-marginal Sinkhorn, see Di [START_REF] Marino | An optimal transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm[END_REF].

We applied this method to the numerical resolution of the Beltrami flow (5.10) with 𝑀 = 16 and 𝑁 = 64 2 . As in Section 5.2, we use three colours to label the mass initially. It is split into three non-overlapping subdomains of 𝑋 0 called 𝑅, 𝐺 and 𝐵 (red, green and blue). This is the first line in Figures 5.2 and 5.3. The first column shows the exact Beltrami flow for increasing times. The last three columns show where this mass has been sent by plotting

𝑥 𝑚 → 𝑃 𝑋 0 ×𝑋 𝑚 #𝜋 𝜖 (𝑅/𝐺/𝐵, 𝑥 𝑚 ),
(5.29) using a variable transparency depending on (5.29). The second column adds the three RGB channels. Figures 5.2 and 5.3 correspond to two different final configurations X 𝑓 given by the exact Beltrami flow at (small) time 𝑇 = 0.9 and (large) time 𝑇 = 𝜋. The diffusion induced by entropic regularization (𝜖 = 1𝐸 -04) smoothing can be observed on the short time geodesic (Figure 5.2). As in Figure 5.1, mass takes a shortest path in the large time simulation (Figure 5.3). There is also mass splitting and mixing, suggesting GIF is indeed a looser relaxation than generalized geodesics (Section 5.2).

41 𝜋 𝜖 is given by (5.25), and (5.29) is the measure of the mass transported from the grid points in regions 𝑅/𝐺/𝐵, respectively, to the grid point 𝑥 𝑚 at the discrete time 𝜏 𝑚 . 42 No mass has perfect transparency: one sees the white background, full mass 1 is the true opaque colour (blue, red, green), and if the mass is in between it is partially transparent.

The Schrödinger problem and transport by diffusion

The Schrödinger problem

At convergence, Sinkhorn equations (3.13) are

𝑢 * 0 = LSE 𝜖 𝜇 1 ,𝑋 1 (𝑢 * 1 ), 𝑢 * 1 = LSE 𝜖 𝜇 0 ,𝑋 0 (𝑢 * 0 ). (6.1)
This is a direct consequence of the optimality conditions (3.22). As mentioned in Section 3.2, the entropic regularization machinery holds in the continuous setting and this will make the exposition easier here. We will also assume that 𝜇 0,1 are smooth densities over R 𝑑 with finite second moments. The change of variable

( 𝑓 * 0 , 𝑔 * 1 ) = (2𝜋𝜖) 𝑑/2 (e 𝑢 * 0 /𝜖 𝜇 0 , e 𝑢 * 1 /𝜖 𝜇 1 ) (6.2)
in (6.1) gives (after a few lines)

𝑓 * 0 𝑔 0 = 𝜇 0 and 𝑓 1 𝑔 * 1 = 𝜇 1 , (6.3)
where ( 𝑓 1 , 𝑔 0 ) are obtained using forward and backward time integration of a heat flow, respectively. The heat flow is expressed using the Laplace operator Δ = div 𝑥 (𝐷 𝑥 ) (in our previous notation) as the generator of the semi-group, that is, (6.4) This reformulation relies on the interpretation of (3.15) as a heat kernel. It implies that the problem is now set on 𝑋 = R 𝑑 (the periodic torus with the corresponding periodic cost 𝑐 is also possible). The heat flows (6.4) have a Lagrangian stochastic differential equation interpretation, ( 𝑓 𝜏 , 𝑔 𝜏 ) are the probability laws of two independent random processes X ± 𝜏 following the standard Brownian motions B 𝜏 :

𝑔 𝜏 = e -(𝜖 𝜏/2) Δ 𝑔 * 1 , 𝜏 : 1 → 0 and 𝑓 𝜏 = e (𝜖 𝜏/2) Δ 𝑓 * 0 , 𝜏 : 0 → 1.
dX ± 𝜏 = ± √ 𝜖 dB 𝜏 , X + 0 ∼ 𝑓 * 0 and X - 1 ∼ 𝑔 * 1 . (6.5)
The Sinkhorn algorithm (3.12) may be interpreted as Picard-type iterations to determine the initial and final laws ( 𝑓 * 0 , 𝑔 * 1 ). We explain below how the density of transport in time is retrieved through the 'interference' product between the forward and backward probability laws (see also Figure 6.1). A similar approach was proposed by [START_REF] Guéant | Mean field games equations with quadratic Hamiltonian: A specific approach[END_REF] in the variational mean field games context (next section).

We know that ( 𝑓 * 0 , 𝑔 * 1 ) are the solutions of the dual problem (3.11). The Schrödinger problem is a probabilistic interpretation of the primal problem associated to (3.11) under the change of variable (6.2). It is a reformulation (again a few lines of calculations) of (3.16):

𝜋 * 𝜖 := arg inf 𝜋 𝜖 ∈Π(𝜇 0 , 𝜇 1 )
𝜖 KL(𝜋 𝜖 |𝜋 0 𝜖 1 𝑋 0 ×𝑋 1 ), (6.6) 43 ℎ 𝜏 (𝑥) = (e (𝜖 𝜏/2) Δ ℎ 0 ) (𝑥) := ∫ 𝛾 (𝜖 𝜏/2) (𝑥 -𝑥 ) ℎ 0 (𝑥 ) d𝑥 . 44 X 𝜏 ∼ 𝜌 𝜏 is to be understood as 'the probability law of X 𝜏 is 𝜌 𝜏 ': 𝑃(X 𝜏 ∈ 𝐴) = 𝜌 𝜏 (𝐴) for any measurable 𝐴 ⊂ 𝑋.

where the transport part of the cost is embedded in 𝜋 0 𝜖 (𝑥 0 , 𝑥 1 ) := 𝛾 𝜖 (𝑥 1 -𝑥 0 ). (6.7)

As in Section 5.3, let us use the time flow in 𝜏. Let 𝑅 𝜖 ∈ P(Ω) denote the Wiener measure. Then 𝜋 0 𝜖 is the density of 𝑅 0→1 𝜖 := (𝑒 0 , 𝑒 1 )#𝑅 𝜖 , the 0 → 1 transition probability measure associated to 𝑅 𝜖 . If, for example, 𝜇 0,1 = 1 𝑋 , i.e. Lebesgue measure, then the solution to (6.6) is simply 𝜋 * 𝜖 = 𝑅 0→1

𝜖

. More generally, 𝜋 * 𝜖 = (𝑒 0 , 𝑒 1 )#𝑄 * 𝜖 , where 𝑄 * 𝜖 solves a dynamic version of (6.6):

𝑄 * 𝜖 := arg inf 𝑄 𝜖 ∈P(Ω), (𝑒 0 ,𝑒 1 )#𝑄 𝜖 ∈Π(𝜇 0 , 𝜇 1 )
𝜖 KL(𝑄 𝜖 |𝑅 𝜖 ). (6.8)

This last step requires a careful analysis, which can be found in [START_REF] Léonard | A survey of the Schrödinger problem and some of its connections with optimal transport[END_REF].

The Lagrangian stochastic differential equation interpretation is interesting. The minimizer of (6.8) is Markovian and the law of a diffusion process with constrained initial and final time law:

dX * 𝜏 = -𝐷 𝑥 𝜙 * (X * 𝜏 ) d𝜏 + √ 𝜖 dB 𝜏 , X * 0,1 ∼ 𝜇 0,1 . (6.9)
The law of X * 𝜏 ∼ 𝑒 𝜏 #𝑄 * 𝜖 = 𝜌 * 𝜏 is the entropic optimal transportation interpolation. The drift 𝐷 𝑥 𝜙 * and volatility can be deduced from the PDE interpretation (6.12) below. As in the case of the generalized incompressible flows, the mass at 𝑥 0 can be split along 𝐶 0 -path solutions of (6.9). The measure 𝑄 * 𝜖 is closest to 𝑅 𝜖 in the sense of KL entropy (known in this context as the Boltzmann-Shannon entropy) and constrained to have 𝜇 0,1 as the initial and final densities in time, respectively. The value function can also be interpreted (using Sanov's theorem: see Léonard 2014) as (minus the log of) the event probability that, given 𝜇 0,1 , the Wiener measure satisfies (𝑒 0 , 𝑒 1 )#𝑅 𝜖 ∈ Π(𝜇 0 , 𝜇 1 ).

The Schrödinger interpretation of optimal transportation is somewhat technical for numerical analysts not well versed in probability theory. There is a more formal PDE interpretation that will, as in the non-entropic case, lead to a CFD-like model and link entropic optimal transportation with diffusion.

The optimality conditions (6.4) are rewritten as

(-𝜕 𝜏 -𝜖 Δ) 𝑔 𝜏 = 0 𝑔 1 = 𝑔 * 1 and (𝜕 𝜏 -𝜖 Δ) 𝑓 𝜏 = 0, 𝑓 0 = 𝑓 * 0 . (6.10)
We apply the 'Hopf-Cole-type' transformation

{ 𝑓 𝜏 , 𝑔 𝜏 } → 𝜌 𝜏 = 𝑓 𝜏 𝑔 𝜏 , 𝑔 𝜏 = exp - 1 2 𝜖 𝜙 𝜏 (6.11) 45 𝑅 𝜖 is a measure on Ω = C([0, 1], 𝑋), 𝑅 ∼ 1 (2𝜋𝜖) 𝑑/2 ∫ Law(𝑥 + √ 𝜖 B) d𝑥,
where B is the standard Brownian motion starting at 0, i.e. the Markov process whose generator is the operator above 𝜖Δ.

Figure 6.1. This is an illustration of the 'interference product' in (6.11) for two Dirac masses 𝜇 0,1 (one at the initial time and the other at the final time) and a diminishing entropy/diffusion parameter 𝜖. It also explains the numerical stability limit of the entropic approach. The normalization is achieved by the potentials; when the heat kernel is numerically 0 instead of extremely small, it is no longer possible.

to (6.10) and obtain the Eulerian version of (6.9):

(𝜕 𝜏 -𝜖 Δ) 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 𝐷 𝑥 𝜙 𝜏 ) = 0, (-𝜕 𝜏 -𝜖 Δ) 𝜙 𝜏 + 1 2 𝐷 𝑥 𝜙 𝜏 2 = 0. (6.12)
The second equation is a Hamilton-Jacobi-Bellman equation and the first equation a Fokker-Planck equation for the density 𝜌 𝜏 . The optimality condition (6.3) becomes the familiar marginal constraints 𝜌 0,1 = 𝜇 0,1 . (6.13)

The Fenchel-Rockafellar formalism (4.20) can be used simply by replacing 𝐴 with 𝐴 𝜖 : (𝜙, 𝜙 0 , 𝜙 1 ) ↦ → ((-𝜕 𝜏 -𝜖 Δ) 𝜙, 𝐷 𝑥 𝜙), the adjoint becoming (𝐴 𝜖 ) : (𝜌, 𝑚) ↦ → -((𝜕 𝜏 -𝜖 Δ) 𝜌 𝜏 + 𝐷 𝑥 𝑚, 𝜌 0 , -𝜌 1 ). The solution (𝜌 * , 𝜙 * ) of (6.12)-(6.13) is the (unique) minimizer of inf (𝜌,V)∈FP(𝜌 0 ,𝜌 1 )

∫ 1 0 1 2 V 𝜏 (𝑥) 2 , 𝜌 𝜏 𝑋 d𝜏, (6.14)
where

FP(𝜌 0 , 𝜌 1 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → (𝜌 𝜏 (𝑥), V 𝜏 (𝑥)) ∈ R + × R 𝑑 : (𝜕 𝜏 -𝜖 Δ) 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 V 𝜏 ) = 0 and 𝜌 𝜏=0,1 = 𝜇 0,1 . (6.15)
It is similar to the CFD formulation (4.15)-(6.15) but, as expected, with the entropic regularization the interpretation in terms of transport maps is completely lost, and the mass is transported along the stochastic paths (6.9).

Variational mean field games

Variational mean field games are a specific subclass of mean field games, introduced in Lasry and Lions ( 2007), which can be tackled with the numerical and theoretical tools used in dynamic optimal transportation. For a general introduction to mean field games and their link with Nash equilibria of multi-agent systems, see [START_REF] Achdou | Mean Field Games[END_REF]. We will restrict ourselves to the simplest variational mean field games generalization of (6.14)-(6.15):

inf (𝜌,V)∈FP(𝜌 0 ) ∫ 1 0 1 2 V 𝜏 (𝑥) 2 , 𝜌 𝜏 𝑋 d𝜏 + ∫ 1 0 𝐻(𝜌 𝜏 ) d𝜏 + 𝐻 1 (𝜌 1 ), (6.16)
where 𝐻 and 𝐻 1 are convex function in C(𝑃(𝑋), R) and

FP(𝜌 0 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → (𝜌 𝜏 (𝑥), V 𝜏 (𝑥)) ∈ R + × R 𝑑 : (𝜕 𝜏 -𝜖 Δ) 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 V 𝜏 ) = 0 and 𝜌 𝜏=0 = 𝜇 0 . (6.17)
Note that 𝜌 is now completely determined by the velocity V 𝜏 though the resolution of the initial value problem in (6.17). The final density 𝜌 1 appears in the cost 𝐻 1 (𝜌 1 ) in (6.16). This formulation belongs to the well-known class of optimal control of a system governed by partial differential equations [START_REF] Lions | Optimal Control of Systems Governed by Differential Equations, Die Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen[END_REF], here a Fokker-Planck equation representing (in the mean field games paradigm) the density of players whose trajectories are subject to white noise and optimized to achieve global minimum cost. Individual trajectories are not observable, only the mean field density. This formulation also holds when 𝜖 = 0, under the name of 'deterministic mean field games'. The Fokker-Planck equation becomes the continuity equation of the CFD formulation.

O , 53

With a slight variation of (4.20),

𝛿 𝜏 → (𝜕 𝜏 -𝜖 Δ) in 𝐴 , 𝐹 ★ → {(𝑎, 𝜌 0 ) ↦ → 𝜄 {𝑎=0,𝜌 0 =𝜇 0 } }, 𝐺 ★ → 𝑞 := (𝜌, 𝑃) ↦ → ∫ 𝑋 ∫ 1 0 𝐽(𝑞) d𝑥 d𝜏 + ∫ 1 0 𝐻(𝜌 𝜏 ) d𝜏 + 𝐻 1 (𝜌 1 ) ,
we can once more apply Fenchel-Rockafellar duality. The optimality system (4.21) becomes

(𝜕 𝜏 -𝜖 Δ)𝜌 * 𝜏 + div 𝑥 (𝜌 * 𝜏 𝐷 𝑥 𝜙 * ) = 0, 𝜌 * 0 = 𝜇 0 , (-𝜕 𝜏 -𝜖 Δ) 𝜙 * 𝜏 + 1 2 𝐷 𝑥 𝜙 * 𝜏 2 = 𝜕 𝜕 𝜌
𝐻(𝜌 𝜏 ), (6.18)

𝜙 * 1 = 𝜕 𝜕 𝜌 𝐻 1 (𝜌 1 ), 𝜌 * 𝜏 a.e. on [0, 1] × 𝑋.
This is the general form of our mean field games. Picking 𝐻 := 0 and the characteristic function 𝐻 1 := 𝜄 {𝜌 1 =𝜇 1 } , for example, brings us back to the Schrödinger problem (Section 6.1).

As in the entropic treatment of the generalized incompressible flows (Section 5.4), we are going to relax (6.16)-( 6.17) to measures on curves. By analogy with the PDE interpretation (4.15)-(6.15) of the Schrödinger problem in its dynamic form (6.8), (6.16)-( 6.17 Going back to the discretization-in-space-and-time technique developed in Section 5.4 and re-using all the notation of Section 5, we get from (5.22) inf

𝜋 𝜖 ∈R 𝑀×𝑁 𝜖 KL(𝜋 𝜖 |𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 ) + 𝑀 ∑ 𝑚=0 𝐻 𝑚 (𝑃 𝑋 𝑚 #𝜋 𝜖 ). (6.20)
This is a slightly generalized version as 𝐻 may now depend on time and 𝐻 𝑀 is the new notation for 𝐻 1 . The {𝐻 𝑚 } are assumed to be convex. The initial condition can be enforced with 𝐻 0 (𝜌) := 𝜄 {𝜌=𝜇 0 } . The gamma-convergence to variational mean field games (6.16)-(6.17) is established in Benamou, Carlier, Di Marino and Nenna (2019b). We again apply the Fenchel-Rockafellar duality (3.19):

𝐴 : 𝜎 := (𝑢 0 , 𝑢 1 , . . . , 𝑢 𝑀 ) ∈ (R 𝑁 ) 𝑀 ↦ → 𝑀 ∑ 𝑚=0 𝑢 𝑚 ∈ R 𝑀 ×𝑁 , 𝐴 : 𝑞 := 𝜋 𝜖 ∈ R 𝑁 ×𝑀 ↦ → (𝑃 𝑋 0 #𝜋 𝜖 , 𝑃 𝑋 1 #𝜋 𝜖 , . . . , 𝑃 𝑋 𝑀 #𝜋 𝜖 ) ∈ (R 𝑁 ) 𝑀 , 𝐹 : (𝑢 0 , 𝑢 1 , . . . , 𝑢 𝑀 ) ∈ (R 𝑁 ) 𝑀 ↦ → 𝑀 ∑ 𝑚=0 𝐻 ★ 𝑚 (𝑢 𝑚 ), 𝐹 ★ : (𝜌 0 , 𝜌 1 , . . . , 𝜌 𝑀 ) ∈ (R 𝑁 ) 𝑀 ↦ → 𝑀 ∑ 𝑚=0 𝐻 𝑚 (𝜌 𝑚 ), 𝐺 : 𝑢 ∈ R 𝑀 ×𝑁 ↦ → 𝐺(𝑢) = 𝜖 e 𝑢/𝜖 , 𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 𝑀 𝑚=0 𝑋 𝑚 , 𝐺 ★ : 𝑞 ∈ R 𝑁 ×𝑀 ↦ → 𝐺 ★ (𝑞) = 𝜖 KL(𝜋 𝜖 |𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 ). (6.21)
We recognize (5.24) with an abstract version of 𝐹. The primal-dual optimality condition (5.25) and Sinkhorn algorithm (5.27) are unchanged but for an important point. The chain of dependence of marginals in the cost 𝑥 0 → 𝑥 1 → • • • → 𝑥 𝑀 is no longer circular. If one is willing to pay the memory cost, the operation cost of one sweep of the Sinkhorn iteration can be taken down to 𝑂(𝑀 𝑁 2 ), and to 𝑂(𝑀 𝑁 1+1/𝑑 ) when using the separability of the Euclidean quadratic cost (1.1). Of course one still needs to iterate to reach convergence.

We illustrate this numerical approach with a simulation from Benamou et al. (2019b). Initial and final densities are prescribed with 𝐻 0, 𝑀 (𝜌) = 𝜄 {𝜌=𝜇 0,1 } . The 'agents' must avoid multiple obstacles moving in time and this is modelled letting 𝐻 𝑚 (𝜌) = 𝜄 𝑋 \𝑂 𝑚 , 𝜌 𝑋 , 𝑂 𝑚 be the moving sets, and we are paying an infinite price if some mass (and therefore agents) is present. In this case 𝐹 is linear, but entropic regularization still yields a strictly convex minimization problem. The boundaries of the obstacles are the white circles in the snapshots displayed in Figures 6.2-6.4. They correspond to different levels of diffusions: 𝜖 = 1, 10 -1 , 10 -2 . The discretization is 𝑀 = 32 and 𝑁 = 128 2 .

Martingale optimal transportation and transport by diffusion

Motivated by applications in finance (see [START_REF] Beiglböck | On a problem of optimal transport under marginal martingale constraints[END_REF][START_REF] Beiglböck | Monotone martingale transport plans and Skorokhod embedding[END_REF][START_REF] Ghoussoub | Structure of optimal martingale transport plans in general dimensions[END_REF] and the references therein), martingale optimal transportation is a recent branch of optimal transportation where the transport plan is constrained to satisfy a martingale constraint. As with 'standard' optimal transportation, it can be formulated as a dynamic problem, but we will start with the static version to simplify the exposition. The vector space 𝑋 ⊂ R 𝑑 describes all possible prices for 𝑑 underlying assets in a portfolio. The marginals 𝜇 0,1 ∈ P(𝑋 0,1 ) represent the state of the market at times 0 and 1, i.e. the distribution over the set of asset prices of a portfolio. The set of transport plans 𝜋 ∈ P(𝑋 0 × 𝑋 1 ) describes all market changes for the distribution of the portfolio between times 0 and 1. The displacement cost 𝑐 is now interpreted as the pay-off (a fixed gain or loss depending on the prices 𝑥 0 ∈ 𝑋 0 ⊂ 𝑋 and 𝑥 1 ∈ 𝑋 1 ⊂ 𝑋). Assuming the market changes are given by a fixed transport plan 𝜋, also called a model, the Monge-Kantorovich cost 𝑐, 𝜋 𝑋 0 ×𝑋 1 is the yield of the option (i.e. the right to buy the portfolio 𝜇 1 at time 1 knowing 𝜇 0 ). The buying (6.22) We recognize the familiar Kantorovich problem (1.2). This approach is called 'model-free hedging' because the sets of admissible changes on the market are really 'free': there are no constraints other than the marginals. The real world is not that simple, and the optimal price 𝑐, 𝜋 * 𝑋 0 ×𝑋 1 may be greatly underestimated. It seems generally agreed that, at least, a 'no-arbitrage' constraint must be added. Using the notation in this paper and assuming we are dealing with densities, it takes (6.23) This is the so-called martingale constraint. It states that considering the initial wealth 𝑥 0 𝜇 0 (𝑥 0 ) at 𝑥 0 , markets cannot 'arbitrage' a strategy to generate more (or less) wealth at time 1 from the portfolio components previously at 𝑥 0 . Wrapping up, martingale optimal transportation is

𝜋 * := arg inf 𝜋 ∈Π Mrt (𝜇 0 , 𝜇 1 )
𝑐, 𝜋 𝑋 0 ×𝑋 1 , (6.24)

where Π Mrt (𝜇 0 , 𝜇 1 ) := {𝜋 ∈ Π(𝜇 0 , 𝜇 1 ), Id 𝑋 1 , 𝜋 𝑋 1 = Id 𝑋 0 𝜇 0 }. (6.25) Because of the additional martingale constraints, Π Mrt (𝜇 0 , 𝜇 1 ) may be empty. Characterizing the class of 𝜇 0,1 such that there exists a martingale 𝜋 ∈ Π Mrt (𝜇 0 , 𝜇 1 ) is the subject of a classical result called Strassen's theorem.

In the mathematical finance literature, martingale optimal transportation is written in probabilistic notation. Let X 0 ∼ 𝜇 0 and X 1 ∼ 𝜇 1 be two random variables on 𝑋. Transport plans 𝜋 are the set of joint laws of (X 0 , X 1 ). Then (6.22) is reformulated as the mathematical expectation inf X 0,1 ∼𝜇 0,1 E(𝑐(X 0 , X 1 )).

(6.26)

The martingale constraint (6.23) is given using the conditional expectation E(X 1 |X 0 ) = X 0 , 𝜇 0 a.e. (6.27)

Departing from finance, let us look at the martingale optimal transportation problem when the displacement cost is the Euclidean quadratic cost (1.1). In the probabilistic formalism (6.26) and using (6.27), we get inf (6.28) the minimization of the conditional variance of X 1 knowing X 0 . This is in line with Monge optimization of mass transport: if mass follows a stochastic path X 0 → X 1 , and (X 0 , X 1 ) is a martingale, minimizing the transport work may only be achieved by acting on the variance or equivalently the volatility of the stochastic process.

X 0,1 𝜇 0,1 E( X 1 -E(X 1 |X 0 ) 2 |X 0 ),
In a striking parallel with Sections 4.1 and 6.1, Huesmann and Trevisan (2019) have provided a dynamic interpretation that sheds new light on the martingale constraint. All martingales (X 0 , X 1 ) ∼ 𝜋 ∈ Π Mrt (𝜇 0 , 𝜇 1 ) can be represented as a joint probability 𝜋 = (𝑒 0 , 𝑒 1 )#𝑄, where 𝑄 is the law of a diffusion process

dX 𝜏 = √ 𝛼 𝜏 dB 𝜏 , X 0,1 ∼ 𝜇 0,1 . (6.29)
Martingale optimal transportation is controlled by the volatility parameter 𝛼 𝜏 . The law of X 𝜏 , again called 𝜌 𝜏 , therefore satisfies the diffusion equation and initial/final boundary conditions:

(𝜕 𝜏 -𝛼 𝜏 Δ) 𝜌 𝜏 = 0, 𝜌 0,1 = 𝜇 0,1 . (6.30)

A dynamic generalization of martingale optimal transportation (6.26) is proposed in Huesmann and Trevisan 2019 related to the CFD formulation. Its simplest version is (the cost 𝑐 needs to be reinterpreted): inf (𝛼 𝜏 ,𝜌 𝜏 ) satisfies (6.30) ∫ 1 0 𝑐(𝛼 𝜏 ), 𝜌 𝜏 𝑋 d𝜏.

(6.31)

Needless to say, the mathematical tool is again convex duality, and the full optimal transportation machinery described in Section 6 applies. At this time, the interplay between the entropic-regularization-induced diffusion, if one uses the Sinkhorn algorithm, and the controlling diffusion is still unclear.

Classical dynamic optimal transport optimizes the velocities, while transport by diffusion optimizes the volatilities. Both can be combined into semi-martingale optimal transportation, where stochastic paths are also controlled using a deterministic drift {𝜏 ↦ → 𝛽 𝜏 }: (6.32) If 𝛼 𝜏 := 0, we are back to standard optimal transportation. If 𝛼 𝜏 is fixed and equal to 𝜖, this is entropic optimal transportation, and if 𝛽 𝜏 := 0, this is the martingale optimal transportation above. For these generalizations, see [START_REF] Léonard | A survey of the Schrödinger problem and some of its connections with optimal transport[END_REF] and [START_REF] Guo | Path dependent optimal transport and model calibration on exotic derivatives[END_REF] and the references therein.

dX 𝜏 = 𝛽 𝜏 d𝜏 + √ 𝛼 𝜏 dB 𝜏 , X 0,1 ∼ 𝜇 0,1 .
O , 59

7. Transport distances as loss/fidelity

Variational Russian dolls

So far we have described static and dynamic optimal transportation problems as linear or convex optimization problems set on the configuration space 𝑋 and the rich links existing between various deterministic or stochastic models, in particular in the dynamic case (Sections 4, 5 and 6). The optimizers of these variational problems were characterized as measures on 𝑋 for the primals and continuous functions on 𝑋 for the duals, respectively. This is our smallest Russian doll.

The intermediate Russian doll is at the P(𝑋) level. In Section 3 we saw that freezing one marginal (say the first 𝜌 0 ) in the W 2 cost gives a mathematical metrization on the space of probability measures P(𝑋). The distance-to-𝜌 0 (𝜌 0 ∈ P(𝑋)) information given by the Wasserstein loss functional

𝜌 ∈ P(𝑋) ↦ → W 2 2 (𝜌 0 , 𝜌) (7.1)
was sufficient to define meaningful (and useful) notions such as interpolation and barycentres in P(𝑋). A comprehensive review of transport-based statistical applications on this space can be found in [START_REF] Kolouri | Transport-based analysis, modeling, and learning from signal and data distributions[END_REF]. Statistical (or machine) learning can be seen as a powerful inference tool in M + (𝑋): it naturally explains the explosion of work connecting optimal transportation to this domain.

Finally, the theory of Wasserstein gradient flows, presented below in Section 7.2, may now be seen as an autonomous research topic and is the subject of several surveys; see [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF] or Santambrogio (2015, §8).

The last (and largest) fascinating Russian doll is the Gromov-Wasserstein distance on the space of metric-measured spaces introduced by [START_REF] Mémoli | Gromov-Wasserstein distances and the metric approach to object matching[END_REF] and [START_REF] Sturm | The space of spaces: Curvature bounds and gradient flows on the space of metric measure spaces[END_REF]. Using our notation, a metric-measured space is a triplet (𝑋, 𝜇, 𝑐), where 𝜇 is a reference measure on 𝑋 and the displacement cost 𝑐 a distance. Gromov-Wasserstein is a generalization of the optimal transportation distance between two such metric-measured spaces (𝑋 0 , 𝜇 0 , 𝑐 0 ) and (𝑋 1 , 𝜇 1 , 𝑐 1 ). It is important for applications where the spaces supporting the distributions 𝜇 0,1 cannot be embedded into a common 𝑋. The problem is known to be NP-hard; see [START_REF] Peyré | Gromov-Wasserstein averaging of kernel and distance matrices[END_REF], who propose an approximate entropic regularization method.

Below we discuss the middle Russian doll.

W 2 gradient flows

We are concerned with functionals defined on P(𝑋). Let us first recap the properties of (7.1). Using the dual formulation (2.1), this is equivalent to

𝜌 ∈ P(𝑋) ↦ → sup (𝑢 0 ,𝑢 1 )∈𝐶 𝐷 𝑢 0 , 𝜌 0 𝑋 0 + 𝑢 1 , 𝜌 𝑋 . (7.2)
The Wasserstein loss functional is the upper envelope of linear functionals in 𝜌, 𝜌 ↦ → 𝑢 0 , 𝜌 0 𝑋 0 + 𝑢 1 , 𝜌 𝑋 , and it is therefore convex. Formally, again, Danskin's theorem (see footnote 15) implies that (7.2) is differentiable and its functional gradient (we are in the intermediate Russian doll environment) is simply given by the Kantorovich optimal potential 𝑢 * 1 :

𝜕{𝜌 ↦ → W 2 2 (𝜌 0 , 𝜌)} 𝜕 𝜌 = 𝑢 * 1 . (7.3)
Looking back at the economics interpretation of the dual Kantorovich problem (Section 2.1), it makes sense. If the producer wants to increase revenue by lobbying the government into modifying the repartition of one of the goods, they just need to at the local prices.

A gradient flows in 𝑋 aims to find a minimizer of some convex differentiable 'energy' 𝑥 ∈ 𝑋 ↦ → 𝐹(𝑥) ∈ R. This is achieved by solving, from some initial X 0 = 𝑥 0 , the ordinary differential equation X 𝜏 = -𝐷 𝑥 𝐹(X 𝜏 ), until it becomes stationary and therefore reaches a minimum. Using an implicit discretization in time, we get

X 𝜏+ d𝜏 -X 𝜏 = -d𝜏 𝐷 𝑥 𝐹(X 𝜏+ d𝜏 ), (7.4)
which, if 𝐹 is strictly convex, admits the well-known variational formulation

X 𝜏+ d𝜏 := arg inf X ∈𝑋 1 2 X 𝜏 -X 2 + d𝜏 𝐹(𝑥). (7.5)
Starting with [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], this concept has been lifted to the space of probability measures P(𝑋) using the Wasserstein distance instead of the Euclidean quadratic cost (1.1):

𝜇 𝜏+ d𝜏 := arg inf 𝜇 ∈P(𝑋 ) 1 2 W 2 2 (𝜇 𝜏 , 𝜇) + d𝜏𝐹(𝜇). (7.6)
According to the dynamic optimal transportation formulation, the density 𝜇 minimizes its kinetic energy plus some penalization depending on the additional energy 𝐹. A sequence of probability measures {𝜌 𝜏 𝑚 } 𝑚=0,... (𝜏 𝑚 = 𝑚 d𝜏) minimizing the energy 𝐹 is built iterating (7.6). The minimization (7.6) is a particular instance of the mean field game (6.16)-(6.17) with 𝐻 := 0 and 𝐻 1 := d𝜏 𝐹, where we have interpreted the W 2 distance with its CFD formulation. Using the mean field game optimality condition (6.18), there is, on all intervals (𝜏 𝑚 , 𝜏 𝑚+1 ), a curve in time of (7.7)

In the limit d𝜏 → 0, it is possible to show rigorously (see Santambrogio 2015, §8) that the curve 𝜏 → 𝜌 𝜏 satisfies, in a weak sense, the initial value problem

𝜕 𝜏 𝜌 𝜏 -div 𝑥 𝜌 𝜏 𝐷 𝑥 𝜕𝐹 𝜕 𝜌 (𝜌 𝜏 ) = 0, 𝜌 0 given. (7.8)
This is a gradient descent of the energy 𝐹 with respect to the geometry of the support space P(𝑋) described by the Wasserstein distance. The simplest and famous example in [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], mathematically consistent with the second law of thermodynamics, is the negative Gibbs entropy 𝐹(𝜌) = KL(𝜌|1 𝑋 ), from which one recovers the dissipative heat equation; see also [START_REF] Gentil | The entropy, from Clausius to functional inequalities[END_REF] for a review of entropy and gradient flows. This approach to functional gradient flows is well documented and has been applied to many nonlinear dissipation/diffusion models, in particular to derive theoretical rates of convergence to equilibria. From the numerical point of view, all optimal transportation numerical techniques, with or without entropic regularization, can be applied and provide, naturally, mass conservation and non-negativity of the density; see Benamou, Carlier and Laborde (2016a), Cancès, Gallouët and Todeschi (2020), [START_REF] Matthes | Convergence of a variational Lagrangian scheme for a nonlinear drift diffusion equation[END_REF], Benamou, Carlier, Mérigot and Oudet (2016c) and [START_REF] Peyré | Entropic approximation of Wasserstein gradient flows[END_REF], among others.

The power of the Wasserstein variational approach is nicely illustrated in a model of crowd motion under congestion proposed in Maury, Roudneff-Chupin, [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF]. Individuals try to exit a room 𝑋. Let the door be a part of the boundary, and denote it by D. They follow the 'closest exit' direction signs posted everywhere in 𝑋. This is given as the gradient of the eikonal 𝐷 𝑥 𝐸, where

𝐸 := 𝑥 ↦ → inf 𝑥 0 ∈D 1 2 𝑥 -𝑥 0 . (7.9)
At the 'microscopic' level, people are modelled by hard spheres of positive radius 𝑅, which cannot overlap. At the macroscopic level, one asks that the mean field density does not exceed a fixed threshold, say 1. The density of the crowd of people may vary but it has a compressibility hard limit. The proposed energy for the Wasserstein gradient flow is While moving, the density strives to minimize the distance to the exit but cannot exceed the compressibility limit. The energy is convex but obviously not differentiable. It is nevertheless shown in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] that the curve-of-densities solution of (7.6) converges to a (weak) solution of

𝐹 := 𝜌 ∈ P(𝑋)) ↦ → 𝐸, 𝜌 𝑋 if 𝜌 ≤ 1 a.e. in
𝜕 𝜏 𝜌 𝜏 -div 𝑥 (𝜌 𝜏 (𝐷 𝑥 𝐸 -𝐷 𝑥 𝑝)) = 0, 𝜌 0 given, 0 ≤ 𝜌 ≤ 1, 𝑝 ≥ 0, 𝑝 (1 -𝜌) = 0, a.e. in 𝑋. (7.11)
The velocity 𝐷 𝑥 𝐸 points to the exit but, because of the congestion constraint and reminiscent of the Euler problem (Section 5), a pressure 𝑝 kicks in locally when the density reaches 1, to correct the trajectories. Figure 7.1 shows a convincing comparison with a microscopic simulation.

Inverse problems and unbalanced optimal transportation

The classical least-squares approach to the solution of ill-posed linear system 𝐴 𝑥 = 𝑏 * is the variational problem 𝑥 * := arg inf 𝑥 ∈𝑋 1 2 𝐴𝑥 -𝑏 * 2 . Likewise, defining an abstract model

Mod := (𝜃, 𝜌 0 ) ∈ Θ × P(X ) ↦ → 𝜌 1 = 𝑇 𝜃 #𝜌 0 ∈ P(𝑋) (7.12)
as a family of transport map 𝑇 𝜃 : P(𝑋) → P(𝑋) parametrized by a set of parameters 𝜃 ∈ Θ, the Wasserstein loss can be used in at least two ways. Assuming a computational definition of the model in terms of 𝜃 and 𝜌 0 (the forward map) is known, then, given some observed input and output of the model (𝜌 0 , 𝜌 1 ), the inverse problem

𝜃 * := arg inf 𝜃 ∈Θ 1 2 W 2 2 (𝑇 𝜃 #𝜌 0 , 𝜌 1 ) (7.13)
stands for the set of model parameters best approximating the observation density.

If, on the contrary, the model's dependence on 𝜃 is unknown or too complicated, another method, now called 'supervised learning', is to replace it with a surrogate convolutional neural network (CNN) model N 𝜃 (𝜌 0 ) Mod(𝜌 0 ) for all 𝜌 0 . The parameters, still denoted by 𝜃, are now independent of the model and characterize the CNN. Given a collection of observations (𝜌 0,𝑖 , 𝜌 1,𝑖 ),

𝜃 * := arg inf 𝜃 ∈Θ 1 2 ∑ 𝑖 W 2 2 (N 𝜃 (𝜌 0,𝑖 ), 𝜌 1,𝑖 ) (7.14)
is the optimal set of parameters for the CNN N 𝜃 * to approximate the model Mod.

The advantages of the transport distance have already been discussed in Section 3. First, as indicated by the name, it is a distance on the space of probability measures allowing a generalization to P(𝑋) of interpolation, barycentres, geodesics and more. It also 'metrizes weak convergence'. This property is important in supervised learning, for example, as the learning samples (𝜌 0,𝑖 , 𝜌 1,𝑖 ) are usually discrete observations: histograms or empirical measures. One expects the loss used in (7.14) to be at least continuous with respect to this sampling process.

There are limits to the use of standard transport distances for such inverse problems. On the curse of dimensionality suffered by the accuracy of the sampling process for W 2 , see Chizat et al. (2020) Vacher, Muzellec, Rudi, Bach and[START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF] and the references therein.

Another serious issue is linked to mass conservation, or rather the frequent lack of it for realistic models. Noise or discrete approximations may be a reason but the model itself may not be conservative: 1 𝑋 , Mod(𝜌 0 ) ≠ 1 𝑋 , 𝜌 0 . The simplest fix is to introduce an additional normalization, i.e. replace (7.12) with

Mod normalized := (𝜃, 𝜌 0 ) ↦ → 𝜌 1 = Mod(𝜌 0 ) 1 𝑋 , Mod(𝜌 0 ) , (7.15)
but it may seriously modify the mass distribution and changes the modelling.

A different approach was proposed in [START_REF] Benamou | Numerical resolution of an 'unbalanced' mass transport problem[END_REF], where the term 'unbalanced transport' was introduced. It corresponds to the mean field game (6.16) (Section 7.2), inf where 𝐻 := 0 and

(𝜌,V)∈FP(𝜌 0 ) ∫ 1 0 1 2 V 𝜏 (𝑥)
𝐻 1 : 𝜌 1 → 1 2 𝜌 1 -𝜇 1 2 L 2 (7.17)
is a relaxation of the target marginal condition. The Fokker-Planck equation constraint (6.17) corresponding to a noisy transport may be used and solved using the entropic optimal transportation machinery presented in Section 6. It is also possible to replace (𝜕 𝜏 -𝜖 Δ) with 𝜕 𝜏 , i.e. 𝜖 = 0, giving the deterministic continuity equation constraints (7.18) and use the proximal splitting methods mentioned in Section 4.2. This is the method followed in [START_REF] Benamou | Numerical resolution of an 'unbalanced' mass transport problem[END_REF]. In (7.16)-( 7.18) all the mass initially distributed as 𝜇 0 is transported. The relaxation (7.17) is the price to pay for the mass default. This part of the cost is static and allows for mass to be created or destroyed anywhere independently of the support 𝜌 0 . It is, however, known to fail to define a distance between non-negative Radon measures (𝜇 0 , 𝜇 1 ) ∈ M + × M + . The fix, independently proposed by Chizat, Peyré, Schmitzer and Vialard (2018b), Liero, Mielke and Savaré (2016) and [START_REF] Kondratyev | A new optimal transport distance on the space of finite Radon measures[END_REF], is based firstly on introducing a reaction term in the constraint allowing for mass creation/destruction:

CE(𝜌 0 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → (𝜌 𝜏 (𝑥), V 𝜏 (𝑥)) ∈ R + × R 𝑑 : 𝜕 𝜏 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 V 𝜏 ) = 0 and 𝜌 𝜏=0 = 𝜇 0 ,
𝑇 𝑅 (𝜌 0 ) := (𝜏, 𝑥) ∈ [0, 1] × 𝑋 ↦ → (𝜌 𝜏 (𝑥), V 𝜏 (𝑥), 𝑟 𝜏 (𝑥)) ∈ R + × R 𝑑 × R : 𝜕 𝜏 𝜌 𝜏 + div 𝑥 (𝜌 𝜏 V 𝜏 ) = 𝑟 𝜏 𝜌 𝜏 and 𝜌 𝜏=0,1 = 𝜇 0,1 . (7.19)
Note that, thanks to the reaction term, the equation is not conservative and we can impose the usual initial/final boundary conditions in time. Secondly, competition is established between transport and reaction in the cost function: Minimizing 𝐹 increases the proportion of dead cells modelled by 𝜌. It expands according to the 'geometry' of the reaction diffusion equation. in (7.10), there is a compressibility limit as the proportion of dead cells cannot exceed 1, resulting in the apparition of a pressure governing the flow of cells. The parameter 𝜆, called the 'homeostatic pressure', corresponds to the equilibrium between natural cell division and cell death. The reaction rate is positive only when the pressure is below 𝜆 and new cells appear; see Figure 7.2.

W 2 FR (𝜇 0 , 𝜇 1 ) := inf (𝜌,V ,𝑟 )∈𝑇 𝑅 (𝜌 0 ) ∫ 1 0 1 2 V 𝜏 (𝑥) 2 + 1 4 𝑟 𝜏 (𝑥)

Sinkhorn divergence

Use of the entropic OT 𝜖 (𝜇 0 , 𝜇 1 ) (defined in (3.11)) as a proxy for OT(𝜇 0 , 𝜇 1 ) is widespread, in particular for the applications presented in the above sections. Unfortunately OT 𝜖 is not a distance on P(𝑋), and in particular OT 𝜖 (𝜇 0 , 𝜇 0 ) > 0.

The identity Monge map cannot be represented by a diffuse entropic plan in the form (3.22): when 𝜖 → +∞, 𝜋 * 𝜖 will tend to the most diffuse admissible transport plan, 𝜇 0 ⊗ 𝜇 0 and OT 𝜖 (𝜇 0 , 𝜇 0 ) → 𝑐, 𝜇 0 ⊗ 𝜇 0 . Assuming, without loss of generality, 51 (•) + = max(0, .) is the positive part.

that the mean of 𝜇 0 is the origin and 𝑐 is the Euclidean quadratic cost (1.1), let us define the sequence of dilated distributions

𝜇 𝜆 := 𝑥 ↦ → 1 𝜆 𝑑 𝜇 0 𝑥 𝜆 for 𝜆 > 1.
Then 𝑐, 𝜇 0 ⊗ 𝜇 𝜆 = 1/2 + 1/𝜆 2 ) 𝑚 2 (𝜇) is strictly decreasing as 𝜆 → +∞. In particular, 𝑐, 𝜇 0 ⊗ 𝜇 0 > 𝑐, 𝜇 0 ⊗ 𝜇 𝜆 .

(7.23)

So OT 𝜖 obviously does not metrize weak convergence (see footnote 50), but also 𝜇 0 is not even guaranteed to be a minimizer of 𝜇 ↦ → OT 𝜖 (𝜇 0 , 𝜇) (depending of course 𝜖, 𝜆).

The tempting way to fix the problem is to decrease 𝜖 to reduce the entropic bias, but the stability and rate of convergence seriously deteriorate when 𝜖 → 0 (see Section 3.2). One remedy (see [START_REF] Feydy | Interpolating between optimal transport and MMD using Sinkhorn divergences[END_REF] and the references therein) is to replace OT 𝜖 with 𝑆 𝜖 (𝜇 0 , 𝜇 1 ) = OT 𝜖 (𝜇 0 , 𝜇 1 ) -1 2 (OT 𝜖 (𝜇 0 , 𝜇 0 ) + OT 𝜖 (𝜇 1 , 𝜇 1 )). (7.24)

We immediately see that, at least, the identity bias discussed previously is removed: 𝑆 𝜖 (𝜇, 𝜇) = 0. This loss is called Sinkhorn divergence, and [START_REF] Feydy | Interpolating between optimal transport and MMD using Sinkhorn divergences[END_REF] proved that it is symmetric in 𝜇 0 and 𝜇 1 , and remains positive and convex with respect to 𝜌 0 and 𝜌 1 . It also metrizes the weak convergence of measures. It yields a better optimal transportation cost proxy than OT 𝜖 at the same computational cost. Details, references and simulations can be found in [START_REF] Ramdas | On Wasserstein two sample testing and related families of nonparametric tests[END_REF] and the {𝜇 𝑚 } are given probability measures. Multi-marginal optimal transportation appeared in Section 3 and underlies all dynamic optimal transportation models (Sections 4, 5 and 6). The multi-marginal cost (4.2) there arises from the discretization of the kinetic energy ((4.15) reverting to (4.12)) and is closely linked to the Euclidean quadratic cost (1.1). et al. (2013) and [START_REF] Buttazzo | Optimal-transport formulation of electronic density-functional theory[END_REF]. All marginals 𝜇 𝑚 are a single 𝜌, the identical density probability distribution for the 𝑀 electrons in a given molecule. The optimal transportation value function is the Coulomb repulsive energy whose minimizers characterize the state of strictly correlated electrons. The dimensionality of this problem is daunting: it is naturally set in dimension 𝑑 = 3 but the plan itself is a probability measure over (R 𝑑 ) ×𝑀 . The sparsity of multi-marginal DFT plans is studied in Friesecke and Vögler (2017a) For a review of results on the sparsity -or lack -of multi-marginal optimal plans for general costs, see [START_REF] Pass | Multi-marginal optimal transport: Theory and applications[END_REF] and Di [START_REF] Marino | Optimal transportation theory with repulsive costs, in Topological Optimization and Optimal Transport[END_REF]. The global interactions between marginals prevents the tensorization gains used for the kernel (5.23); Altschuler and Boix-Adsera (2020) studied the complexity of available algorithms with respect to the structure of the multi-marginal cost. An entropic numerical approach was carried out by Benamou, Carlier and Nenna (2016b). Relaxations/simplifications are still under investigation; see [START_REF] Alfonsi | Approximation of optimal transport problems with marginal moments constraints[END_REF], Friesecke andVögler (2017b, 2018) and [START_REF] Cotar | Infinite-body optimal transport with Coulomb cost[END_REF], among others.

L 1 optimal transportation, the Beckman problem and optic flows

Contrary to the Euclidean quadratic cost (1.1), the L 1 optimal transportation cost (originally considered by Monge)

𝑐(𝑥 0 , 𝑥 1 ) = 𝑥 1 -𝑥 0 (8.4)

does not lead to a well-posed primal Kantorovich problem (1.2) and Monge map solutions (2.10). This is well documented in Santambrogio (2015, §3.1) and Peyré and Cuturi (2019, §6). The uniqueness issue is easily explained using the one-dimensional 'bookshelf' example. Books (all with the same weight) are arranged on a shelf and there is just one free spot at the furthest right. The librarian is a maniac and only admits free spots at the furthest left. You are in charge of fixing this problem but you are lazy. Will you pick the book on the furthest left and place it at the furthest right? Or will you shift all the books on the shelf one by one? The Euclidean quadratic cost (1.1) will tell you the second solution is optimal while the cost (8.4) says both strategies require the same amount of work. The only important choice is the direction of transport, and for 𝑑 = 1 this is an easy binary choice.

The optimality conditions (2.4) still hold. Using the distance property of 𝑐 and assuming 𝑋 0 = 𝑋 1 = 𝑋, the Kantorovich potentials can be shown to be 1-Lipschitz and complementary: 𝑢 * 0 = -𝑢 * 1 . The dual Kantorovich problem (2.1) may be tightened to in homogeneous space. Here 𝜙 * is the static phase and 𝑞 * = ∫ 1 0 𝐴 2 (𝜏, 𝑥) d𝜏 is the energy travelling through 𝑥. Problem (1.2) finds the optimal transportation kinematics 𝜙 * with prescribed initial/final amplitudes 𝐴 2 ({0, 1}, •) = 𝜇 0,1 .

A different high-frequency wave asymptotic analysis linking the CFD formulation to a paraxial approximation of the Helmholtz equation has been carried out by [START_REF] Rubinstein | A variational principle in optics[END_REF].

The link with optics is natural as the cost (8.4) is simply the length of rays (𝑥 0 → 𝑥 1 ), i.e. the travel time with index of refraction 𝑔 := 1. It can be replaced by the general Riemannian distance 𝑐 := 𝑑 𝑔 (𝑥 0 , 𝑥 1 ) = inf (8.9) yielding minimum travel times of rays from 𝑥 0 to 𝑥 1 . The geometry of the rays depends on 𝑔: the eikonal equation in (8.8) generalizes to 𝐷 𝑥 𝜙 * = 𝑔. See Figure 8.1 for an illustration in different configurations. More general metrics are possible, for example Finsler metrics [START_REF] Benamou | A numerical solution to Monge's problem with a Finsler distance as cost[END_REF] and nonlinear congestion problems [START_REF] Benmansour | Numerical approximation of continuous traffic congestion equilibria[END_REF], where 𝑔 depends on the energy 𝑞 * .

{X 𝜏 ∈𝑊 1,1 ([0,1],𝑋 ), X 0 =𝑥 0 , X 1 =𝑥 1 } ∫ 1 0 𝑔(X 𝜏 ) X 𝜏 d𝜏,

Signed measures and full seismic waveform inversion

Transport distances have been used recently to attack the problem of full seismic waveform inversion. This is a well-known nonlinear inverse problem: A simple computation shows that a shift in the distribution support 𝜏 ↦ → W 2 2 (𝜇 0 (•), 𝜇 0 (. + 𝜏)) = (𝜏 2 )/2 is convex and even quadratic for the W 2 distance. [START_REF] Yang | Analysis of optimal transport and related misfit functions in full-waveform inversion[END_REF] therefore proposed replacing L 2 with W 2 (8.10) to convexify the problem and eliminate cycle skipping. Time signals are, however, not probability measures and some data transformation is needed. There is little understanding of the interpretation of time oscillatory signals in terms of mass or probability measures. In particular, they are not positive. Just considering the energy of the signal to fix the problem (but losing the phase information), the source/target may also not have the same total mass because of acquisition noise or dissipation in the forward model Mod 𝜃 .

𝜃 * := arg inf 𝜃 1 2 𝑆 𝑡 -
The problem of defining an optimal transport distance for general signed measures is discussed in [START_REF] Mainini | A description of transport cost for signed measures[END_REF]. The simplest idea is to split the positive and the negative part of the signals and consider the sum of the Wasserstein distances (𝜇 0 , 𝜇 1 ) ↦ → W 2 2 (𝜇 + 0 , 𝜇 + 1 ) + 𝑊 2 2 (𝜇 - 0 , 𝜇 - 1 ). (8.11)

If this expression makes sense, i.e. the positive and negative parts are 'balanced' (the same total mass), then this cost remains a distance on P(𝑋). If not, then one may try replacing W 2 with the unbalanced W FR distance (7.20) as in Li, Lamoureux and Liao ( 2020), but it does not define a proper distance. Several signal transformations and normalizations are investigated in [START_REF] Yang | Analysis of optimal transport and related misfit functions in full-waveform inversion[END_REF].

The second idea is to use the W 1 distance (8.5). It remains well-defined and a distance for signed measures. There are no positivity constraints but the (signed) total mass balance 1, 𝜇 -𝜇 0 = 0 (8.12) needs to be satisfied. Using positive and negative parts, a simple computation shows that W 1 (𝜇 0 , 𝜇 1 ) = W 1 (𝜇 + 0 + 𝜇 - 1 , 𝜇 - 0 + 𝜇 + 1 ). (8.13)

It is not clear how to interpret this formula, as the distance optimizes transport between composite distributions aggregating information from the source and the target data. It was used with convincing results in [START_REF] Métivier | An optimal transport approach for seismic tomography: Application to 3D full waveform inversion[END_REF], where the total mass balanced constraint is relaxed.

The geophysical community represents seismograms as collections of independent time 'lines' (the time recordings at each location). Instead of considering the total acquisition as a time × surface image, summing line-wise the optimal transportation distance between observed and simulated lines at the same receivers seems to perform well, and to decrease the dimension of the optimal transportation problems to be solved. Lines may also be interpreted as shapes in the time × signal amplitude 'graph space' [START_REF] Métivier | A graph space optimal transport distance as a generalization of 𝐿 𝑝 distances: Application to a seismic imaging inverse problem[END_REF]. After a discretization in time, they give empirical measures living in R time × R amplitude space.

  Figure 2.1. The simplest non-trivial optimal transport maps are affine and correspond to data in the form of a translation and a dilation 𝑎𝜇 1 (𝑎 𝑥 0 + 𝑏) = 𝜇 0 (𝑥 0 ) (here 𝑎 = 0.5 and 𝑏 = -0.05); 𝜋 is null outside the graph {(𝑥 0 , 𝑥 1 = 𝑇(𝑥 0 ) := 𝑎 𝑥 0 + 𝑏)), 𝑥 0 ∈ 𝑋 0 }.

  Figure 2.2. (a) Different 𝜇 1 distributions for a fixed source 𝜇 0 . (b) Ray-tracing resimulation via Monge-Ampère computation of the reflector. Figure reproduced from Brix, Hafizogullari and Platen (2015) with permission. Copyright © 2015 World Scientific.

  Figure 2.3. The target distribution 𝜇 1 is a picture of Monge. Calculations were done with 𝑁 = 1000 paraboloids for the first row and 𝑁 = 15 000 paraboloids for the second row. (a) Paraboloid intersection diagram (the Laguerre cells) for an initial (𝑢 0,𝑛 ) 1≤𝑛≤ 𝑁 . (b,d) Final intersection diagram after optimization. (e) Reflector surface defined by the intersection of paraboloids. (c,f) Simulation of the illumination at infinity from a point light source illuminating the northern hemisphere (S 2 ) + uniformly, using L R , a physically accurate ray-tracing engine. Figure reproduced from de Castro, Mérigot and Thibert (2016) with permission. Copyright © 2016 Springer.

  Figure 3.1. Blue, 𝜇 0 ; red, 𝜇 1 ; black, 𝜌 𝜏 . The interpolation for (a) the W 2 distance and (b) the L 2 distance.

  Figure 3.2. Blue, 𝜇 0 ; red, 𝜇 1 . The colour map corresponds to 𝜋 *𝜖 for different values of 𝜖. The black line is the optimal transport map 𝑇 (𝜖 = 0). Observe how the entropic plan 'turns', 'aligns' and 'shrinks' from the product measure 𝜇 0 ⊗ 𝜇 1 to the graph of the transport map as 𝜖 decreases to 0.

Figure 3

 3 Figure 3.3. OT 𝜖 barycentres (the interior) and interpolations (the edges) of four densities (the corners) for different sets of (𝜆 𝑝 ) 𝑝=1..4 . Figure reproduced from Peyré and Cuturi (2019) with permission.

  Figure 5.1. (a-e) Beltrami flow in the unit square at various time steps, a classical solution to Euler's equation. The colour of the particles depends on their initial position. (f-j, k-o, p-t, u-y) Generalized fluid flows that are reconstructed for different final times 𝑡 max , i.e. using boundary conditions displayed in the first and last columns. When 𝜏 := 𝑡 max < 1, we recover the classical flow, while for 𝑡 max ≥ 1 the solution is no longer classical and includes mixing.Figure reproduced from Mérigot and Mirebeau (2016) with permission. Copyright © 2016 Society for Industrial and Applied Mathematics. All rights reserved.

  Figure 5.1. (a-e) Beltrami flow in the unit square at various time steps, a classical solution to Euler's equation. The colour of the particles depends on their initial position. (f-j, k-o, p-t, u-y) Generalized fluid flows that are reconstructed for different final times 𝑡 max , i.e. using boundary conditions displayed in the first and last columns. When 𝜏 := 𝑡 max < 1, we recover the classical flow, while for 𝑡 max ≥ 1 the solution is no longer classical and includes mixing.Figure reproduced from Mérigot and Mirebeau (2016) with permission. Copyright © 2016 Society for Industrial and Applied Mathematics. All rights reserved.

  Figure 5.2. Final time, 𝑇 = 0.9. Columns, classical colour tracking of the Lagrangian solution with no mixing in the first column, 𝑃 𝑅 + 𝑃 𝐺 + 𝑃 𝐵 in the second column and 𝑃 𝑅 /𝑃 𝐺 /𝑃 𝐵 in the remaining three columns. Rows (a-e), time evolution. The final Lagrangian configuration at the bottom left is the final datum 𝑋 𝑇 in 𝜋 0,𝑇 = (Id, 𝑋 𝑇 ) # 1. Figure reproduced from Benamou et al. (2019a) with permission. Copyright © 2019 Springer.

  Figure 5.3. Final time, 𝑇 = 𝜋. Columns, classical colour tracking of the Lagrangian solution with no mixing in the first column, 𝑃 𝑅 + 𝑃 𝐺 + 𝑃 𝐵 in the second column and 𝑃 𝑅 /𝑃 𝐺 /𝑃 𝐵 in the remaining three columns. Rows (a-e), time evolution. The final Lagrangian configuration at the bottom left is the final datum 𝑋 𝑇 in 𝜋 0,𝑇 = (Id, 𝑋 𝑇 ) # 1. Figure reproduced from Benamou et al. (2019a) with permission. Copyright © 2019 Springer.

  ) has its own dynamic Schrödinger version: inf 𝑄 𝜖 ∈P(Ω), 𝑒 0 #𝑄 𝜖 =𝜇 0 𝜖 KL(𝑄 𝜖 |𝑅 𝜖 ) + ∫ 1 0 𝐻(𝜌 𝜏 ) d𝜏 + 𝐻 1 (𝜌 1 ). (6.19)

  Figure 6.2. Planning mean field games on the torus with moving obstacles and densities at different time steps; 𝜖 = 1 and 32 time steps. Figure reproduced from Benamou et al. (2019b) with permission. Copyright © 2019 World Scientific.

  Figure 6.3. Planning mean field games on the torus with moving obstacles, fixed initial and final densities similar to Figure 6.2 and densities at different time steps; 𝜖 = 0.1 and 32 time steps. Figure reproduced from Benamou et al. (2019b) with permission. Copyright © 2019 World Scientific.

  Figure 6.4. Planning mean field games on the torus with moving obstacles and densities at different time steps; 𝜖 = 0.01 and 32 time steps. Figure reproduced from Benamou et al. (2019b) with permission. Copyright © 2019 World Scientific.

  Figure 7.1. (a) 'Microscopic' agent simulation at a fixed time. The colour map indicates the pressure. (b) The 'macroscopic' gradient flow at the same time. The density saturates to 1 in the darkest regions. Figure reproduced from Maury et al. (2011) with permission. Copyright © 2011 American Institute of Mathematical Sciences.

  Figure 7.2. Evolution on a non-convex two-dimensional domain ('bone' obstacles in black). (a) Evolution of the density. The colour map is linear from white to blue as the density goes from 0 to 1. (b) Pressure represented by a striped colour map to make the level sets apparent. The white area corresponds to 𝑝 = 𝜌 = 0. Figure reproduced from Di Marino and Chizat (2020) with permission. Copyright © 2020 Société de Mathématiques Appliquées et Industrielles.

  Figure 7.3 shows the comparison of a gradient descent 𝜇 𝜏+ d𝜏 := 𝜇 𝜏+ d𝜏 -d𝜏 𝜕 𝜕𝜇 Loss(𝜇, 𝜇 1 ) (7.25)for different losses and 𝜇 0 ≠ 𝜇 1 . It tests the ability of the Loss to 'drive' the distribution 𝜇 0 to the 'minimizer' 𝜇 1 and whether it metrizes the weak convergence. As expected, for the same regularization 𝜖 = 0.1, 𝑆 𝜖 performs much better than OT 𝜖 .8.A few missing topics amongst many . . .8.1. Multi-marginal optimal transportation and DFTThe prototype 𝑀-marginal problem is inf𝜋 𝑀 ∈Π 𝑀 𝑐 𝑀 , 𝜋 𝑀 𝑀 𝑚=0 𝑋 𝑚 , (8.1)52 𝑚 2 (𝜇) = Id 2 , 𝜇 is the second moment of 𝜇.

Figure 7 . 3 .

 73 Figure 7.3. Gradient descent (7.25) for different losses; 𝜇 0 is in green, 𝜇 1 in blue and 𝜇 𝜏 in red. Figure reproduced from Feydy et al. (2019) with permission.

W 1 (

 1 𝜇 0 , 𝜇 1 ) := sup {𝑢 0 , 𝐷 𝑥 𝑢 0 ≤1} 𝑢 0 , 𝜇 1 -𝜇 0 .(8.5)This new formulation has a well-posed primal called the Beckman problem. We can again use the Fenchel-Rockafellar formalism (3.19) and adapt dynamic CFD version (4.20) to the 'static' version:𝐴 : 𝜎 := 𝜙 ∈ C(𝑋) ↦ → -𝐷 𝑥 𝜙, 𝐴 : 𝑞 ∈ 𝑋(𝑋) 𝑑 ↦ → div 𝑥 𝑞, 𝐹 : 𝜙 ↦ → 𝜙, 𝜇 1 -𝜇 0 𝑋 , 𝐹 ★ : 𝑎, ↦ → 𝜒 𝑎=𝜇 1 -𝜇 0 , 𝐺 : 𝑏 ↦ → 1 𝑏 ≤1 ,𝐺 ★ : 𝑞 ↦ → choice the primal-dual optimality conditions (3.20) take the form div 𝑥 (𝑞 * ) = 𝜇 1 -𝜇 0 , on 𝑋, 𝑞 * = 𝑞 * 𝐷 𝑥 𝜙 * , 𝑞 a.e. on 𝑋.(8.7)For a rigorous derivation of these equations and the space settings, see Santambrogio (2015, §4.3); boundary issues in particular are delicate. The system (8.7) can be interpreted as the eikonal and transport equations of geometric optics𝐷 𝑥 𝜙 * = 1, div 𝑥 ( 𝑞 * 𝐷 𝑥 𝜙 * ) = 𝜇 1 -𝜇 0 (8.8)arising from the high-frequency wave equation asymptotic ansatz 𝑓 (𝜏, 𝑥) 𝐴(𝜏, 𝑥)𝑆(𝜏 -𝜙(𝑥))

  Figure 8.1. Arrow correspond to (1/𝑔) 𝐷 𝑥 𝜙 * and their length is proportional to 𝑞 * . Level curves correspond to the right-hand side density term of the divergence 𝜇 0,1 source target data to be transported (Gaussian densities except in (b), where the target is a collection of three Dirac masses). (a,b) 𝑔 := 1; (c,d) 𝑔 is a lens (c) and a two-layer (d) medium. Figure reproduced from Benamou and Carlier (2015) with permission. Copyright © 2015 Springer Nature.

  .27) See (2.10) for the formulation of 𝑇 as a function of 𝑐 and 𝑢 * 0 . The approximation 𝜋 𝑘 𝜖 𝜖 is to be understood as (3.22) for the running (3.13) iterates. The parameter 𝑝 is positive and depends on the marginals but is unfortunately not explicit. Estimate (3.27) is nevertheless important for two reasons. First it shows that the entropic bias decreases exponentially fast around the continuous transport map (see Figure 3.2). So even though 𝜋 𝑘 𝜖

  𝑃 𝑋𝑚 # 𝜋 𝜖 =1 𝑋 𝑁 + 𝜖 KL(𝜋 𝜖 |𝜋 0 𝜖 1 Ω 𝑀 , 𝑁 ), , 𝑥 1 , . . . , 𝑥 𝑀 ) := 𝛾 𝜖 /𝜆 (𝑥 𝑀 -X 𝑓 (𝑥 0 )) , 𝑥 1 , . . . , 𝑥 𝑀 ) ∈ Ω 𝑁 , 𝑀 and d𝜏 = 1/(2𝑀) a time step. Note that this is now a discrete optimization problem over 𝑀 × 𝑁 real tensors.The generalization of Fenchel-Rockafellar duality (3.19) (see also (3.21) on the larger multi-marginal space) is relatively straightforward:

	can be rewritten as
	inf			(5.22)
	𝜋 𝜖 ∈R 𝑀×𝑁
	where (𝛾 𝜖 is a Gaussian kernel: see (3.15); note that we have normalized the
	constant to 1) and
	𝜋 0 𝜖 (𝑥 0 𝑀 -1 𝑚=0	𝛾 𝜖 /d𝜏 (𝑥 𝑚+1 -𝑥 𝑚 )	(5.23)
	for (𝑥 0
			, we propose to use the entropic
	regularization (as presented in Section 3.3) and adapt the Sinkhorn algorithm (3.12)
	to solve (5.19)-(5.20). Dropping some of the indices for clarity,
	inf 𝜋 𝜖 ∈Π 𝑀 , 𝑁	𝑐 𝑀 ,X 𝑓 , 𝜋 𝜖 𝑀 𝑚=0 𝑋 𝑁 + 𝜖 KL(𝜋 𝜖 |1 Ω 𝑀 , 𝑁 )	(5.21)

40

𝑃 𝑋 𝑚 #𝜋 𝑀 , 𝑁 is the discrete analogue of 𝑒 𝜏 𝑚 #𝜋. 𝑀 ∑ 𝑚=0 𝜒

  48 The Fréchet derivative of 𝜌 ∈ P(X ) ↦ → 𝐹(𝜌), if it exists, is defined as, for any variation 𝜉 such that 𝜌 + 𝜉 ∈ P(𝑋), Reparametrizing time, it satisfies𝜕 𝜏 𝜌 * 𝜏 + div 𝑥 𝜌 * 𝜏 𝐷 𝑥 𝜙 * d𝜏 = 0, 𝜌 * 𝜏=𝜏 𝑚 , 𝜏 𝑚+1 = (𝜇 𝑚 , 𝜇 𝑚+1 ),

	1 d𝜏	𝜙 * 𝜏 𝑚+1 =	𝜕 𝜕 𝜌	𝐹(𝜇 𝜏 𝑚+1 ).
				lim 𝜉 →0	𝜉 1	𝐹(𝜌 + 𝜉) -𝐹(𝜌) =	𝜕 𝜌 𝜕𝐹	(𝜌), 𝜉

𝑋

.

probability measures {𝜇 * 𝜏 }.

  FR (𝜇 0 , 𝜇 1 ) is a well-defined distance on M + (𝑋). It shares a lot of the theory: a static Kantorovich formulation, entropic version, geodesics and barycentres. A gradient flow illustration follows.Di[START_REF] Marino | A tumor growth model of Hele-Shaw type as a gradient flow[END_REF] derived the free boundary 'Hele-Shaw' tumour growth model proposed in[START_REF] Mellet | A Hele-Shaw problem for tumor growth[END_REF] as the W FR gradient flow (replace W 2 with W FR in (7.6)) for the energy 𝐹 := 𝜌 ∈ M + (𝑋) ↦ →

		2 , 𝜌 𝜏	d𝜏. (7.20)
		𝑋
	The problem (7.20)-(7.19) remains convex and the Wasserstein-Fisher-Rao dis-
	tance W -𝜆 𝜌(𝑋) if 𝜌 +∞ else.	1 𝑋 ,	(7.21)

  The Coulomb 'repulsive cost'𝑐 𝑀 (𝑥 1 , . . . , 𝑥 𝑀 ) = |𝑥 𝑖 -𝑥 𝑗 | (8.3)appears in density functional theory, a branch quantum chemistry; see Friesecke

	𝑀 ∑	𝑀 ∑	1
	𝑖=1	𝑗>𝑖

  𝜇 1 2 L 2 , 𝑆 𝑡 = Mod 𝜃 (𝜇 0 ) (8.10) where one tries to recover the underground parameters 𝜃 * of a 'forward' wave propagation model Mod 𝜃 * mapping a wave source 𝑆 𝑡 (a surface explosion at a fixed source) to a 'seismogram' 𝜇 1 (time × surface signal recordings at fixed receivers). The classical least-squares loss in (8.10) is known to lead to an illposed minimization problem with many local minima. A pathology is seen, called 'cycle skipping', linked to phase shifts in the 𝜃 * -observed (or 𝜃-guessed) oscillatory signals.

(•) • (•) denotes the standard scalar product in R 𝑑 .

{𝑥 ∈ 𝑋 ↦ → 𝑓 (𝑥) ∈ 𝑌 }, a function with domain 𝑋 and valued in 𝑌 .

𝐷 𝑥 0 𝑐 is the gradient with respect to the first variable 𝑥 0 , and T 𝑥 0 𝑋 0 is a tangent vector to 𝑋 0 at 𝑥 0 .

𝐷 𝑥 𝑓 (𝑥) is the gradient of 𝑓 at 𝑥.

A two-dimensional code is available at https://gforge.inria.fr/scm/browser.php?group_id=9995.

Available software: https://www.cgal.org/, http://alice.loria.fr/index.php/software/4-library/75geogram.html, https://github.com/mrgt/MongeAmpere.

Written for vectors in the ambient space R 𝑑 .

A parabola reflects all incoming rays in the direction of its axis. See https://en.wikipedia.org/wiki/Parabola#Proof_of_the_reflective_property.

The marginal distributions are defined using the canonical projections 𝑃 𝑋 0 : (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) → 𝑥 0 and 𝑃 𝑋 0 #𝜋 3 = 1, 𝜋 3 𝑀 ×𝑋 1 and likewise 𝑃 𝑋 1 #𝜋 3 = 1, 𝜋 3 𝑋 0 ×𝑋 , 𝑃 𝑋 #𝜋 3 = 1, 𝜋 3 𝑋 0 ×𝑋 1 .

𝑃 is again a canonical projection: 𝑃 𝑋 0 ,𝑋 1 : (𝑥 0 , 𝑥 𝜏 , 𝑥 1 ) → (𝑥 0 , 𝑥 1 ).

An efficient off-the-shelf GPU implementation is available at https://www.kernel-operations.io/geomloss/.

𝜄 {𝑣 ≤0} := 0 if 𝑣 ≤ 0 and +∞.

lim 𝜖 →0 𝜖 e 𝑣/𝜖 = 𝜄 {𝑣 <0} .

This is a good exercise.

LSE is also known as 'soft-min': lim 𝜖 →0 LSE 𝜖 1 , 𝑋 (𝑔)(•) = inf 𝑥 ∈𝑋 𝑐(•, 𝑥) -𝑔(𝑥).

Diam 𝜌 is the diameter of supp 𝜌.

Also known as the relative entropy between 𝜈 and 𝜈 0 . It is strictly convex, vanishes and takes its minimum at 𝜈 0 , and has an infinite slope at 0. Its Fréchet derivative is formally given by 𝛿 KL(𝜈|𝜈 0 ), 𝛿𝜈 = log(𝜈/𝜈 0 ), 𝛿𝜈 .

The precise result involves interpolation of the discrete potential with the 𝑐-transform formula (2.10).

X 𝜏 is the time 𝜏 derivative along the curve 𝜏 → X 𝜏 (𝑥 0 ).

This is a classical result obtained formally by integration by parts (brackets denote the L 2 scalar product and I have skipped discussion of the boundary conditions), that is,V, 𝐷 𝑥 𝑝 =div 𝑥 V, 𝑝 L 2 (𝑋 ) = 0.

The characteristic function 𝜄 𝐵 (𝑥) = 0 if 𝑥 ∈ 𝐵 and +∞ else. 𝑋 \ 𝐵 is the complement of the set 𝐵 in 𝑋.

Formally, Π Mrt (𝜇 0 , 𝜇 1 ) ≠ ∅ if 𝜇 1 dominates 𝜇 0 in the convex order, meaning 𝜙, 𝜌 1 𝑋 1 ≥ 𝜙, 𝜌 0 𝑋 0 for all convex functions 𝜙.

A system particle evolves towards thermodynamic equilibrium by maximizing Gibbs entropy.

The calculation can be done using high-frequency asymptotics (see Symes 1998, §5, for example).

The wave speed is constant equal to 1, and 𝜏 ↦ → 𝑆(𝜏) is a wavelet in time triggered at a point source 𝑥 𝑆 . The acoustic wave equation is (𝜕 2 𝜏 -Δ) 𝑓 = 𝑆(𝜏) 𝛿 𝑥 𝑆 .

I am still wondering if it may be linked to high-frequency asymptotics of a nonlinear auto-(de)focusing optic model.

56 𝜇 + = max {𝜇, 0} and 𝜇 -= max {-𝜇, 0}.