Section: Research Program

Context

Optimal Mass Transportation is a mathematical research topic which started two centuries ago with Monge's work on “des remblais et déblais". This engineering problem consists in minimizing the transport cost between two given mass densities. In the 40's, Kantorovitch [54] solved the dual problem and interpreted it as an economic equilibrium. The Monge-Kantorovitch problem became a specialized research topic in optimization and Kantorovitch obtained the 1975 Nobel prize in economics for his contributions to resource allocations problems. Following the seminal discoveries of Brenier in the 90's [23] , Optimal Transportation has received renewed attention from mathematical analysts and the Fields Medal awarded in 2010 to C. Villani, who gave important contributions to Optimal Transportation and wrote the modern reference monograph [75] , arrived at a culminating moment for this theory. Optimal Mass Transportation is today a mature area of mathematical analysis with a constantly growing range of applications (see below).

In the modern Optimal Mass Transportation problem, two probability measures or "mass" densities : $d{\rho }_i\left(x_i\right)(={\rho }_i\left(x_i\right)d{x}_i),i=0,1$ such that ${\rho }_i\ge 0$, ${\int }_{X_0}{\rho }_0\left(x_0\right)d{x}_0={\int }_{X_1}{\rho }_1\left(x_1\right)d{x}_1=1$, $X_i\subset R^n$. They are often referred to, respectively, source and target densities, support or spaces. The problem is the minimization of a transportation cost, $ℐ\left(M\right)={\int }_{X_0}c(x,M\left(x\right)){\rho }_0\left(x\right)dx$ where $c$ is a displacement ground cost, over all volume preserving maps $M\in ℳ$ $ℳ=\{M:X_0\to X_1,M_{\#}d{\rho }_0=d{\rho }_1\}$. Assuming that $M$ is a diffeomorphism, this is equivalent to the Jacobian equation $det\left(DM\left(x\right)\right){\rho }_1\left(M\left(x\right)\right)={\rho }_0\left(x\right)$ . Most of the modern Optimal Mass Transportation theory has been developed for the Euclidean distance squared cost $c(x,y)={\parallel x-y)\parallel }^2$ while the historic monge cost was the simple distance $c(x,y)=\parallel x-y\parallel$.

In the Euclidean distance squared ground cost, the problem is well posed and in the seminal work of Brenier [24] , the optimal map is characterized as the gradient of a convex potential ${\phi }^{*}$ : $ℐ\left(\nabla {\phi }^{*}\left(x\right)\right)={min}_{M\in ℳ}ℐ\left(M\right)$. A formal substitution in the Jacobian equation gives the Monge-Ampère equation $det\left(D^2{\phi }^{*}\right){\rho }_1\left(\nabla {\phi }^{*}\left(x\right)\right)={\rho }_0\left(x\right)$ complemented by the second boundary value condition $\nabla {\phi }^{*}\left(X_0\right)\subset X_1$. Caffarelli [29] used this result to extend the regularity theory for the Monge-Ampère equation. He noticed in particular that Optimal Mass Transportation solutions, now called Brenier solutions, may have discontinuous gradients when the target density support $X_1$ is non convex and are therefore weaker that than the Monge-Ampère potentials associated to Alexandrov measures (see [50] for a review of the different notions of Monge-Ampère solutions). The value function $\sqrt{ℐ\left(\nabla {\phi }^{*}\right)}$ is also known to be the Wasserstein distance $W_2({\rho }_0,{\rho }_1)$ on the space of probability densities, see [75] . The Computational Fluid Dynamic formulation proposed by Brenier and Benamou in [2] introduces a time extension of the domain and leads to a convex but non smooth optimization problem : $ℐ\left(\nabla {\phi }^{*}\right)={min}_{(\rho ,V)\in 𝒞}{\displaystyle {\int }_0^1{\int }_X\frac{1}{2}\rho (t,x){\parallel V(t,x)\parallel }^{2}dxdt}.$ with constraints : $𝒞=\{(\rho ,V),\text{s.t}{\partial }_t\rho +div\left(\rho V\right)=0,\rho (\{0,1\},.)={\rho }_{\{0,1\}}(.)\}$. The time curves $t\to \rho (t,.)$ are geodesics between ${\rho }_0$ and ${\rho }_1$ for the Wassertein distance. This formulation is a limit case of Mean Fields games [55] , a large class of economic models introduced by Lasry and Lions. The Wasserstein distance and its connection to Optimal Mass Transportation also appears in the construction of semi-discrete Gradient Flows. This notion known as JKO gradient flows after its authors in [52] is a popular tool to study non-linear diffusion equations : the implicit Euler scheme ${\rho }_{k+1}^{dt}=argmi{n}_{\rho (.)}F\left(\rho (.)\right)+\frac{1}{2dt}{W}_2{(\rho (.),{\rho }_k^{dt})}^2$ can be shown to converge ${\rho }_k^{dt}(.))\to {\rho }^{*}(t,.)$ as $dt\to 0$ to the solution of the non linear continuity equation ${\partial }_t{\rho }^{*}+div\left({\rho }^{*}\nabla (-\frac{\partial F}{\partial \rho }\left({\rho }^{*}\right))\right)=0,{\rho }^{*}(0,.)={\rho }_0^{dt}(.)$. The prototypical example is given by $F\left(\rho \right)={\int }_X\rho \left(x\right)log\left(\rho \left(x\right)\right)+\rho \left(x\right)V\left(x\right)dx$ which corresponds to the classical Fokker-Planck equation. Extensions of the ground cost $c$ have been actively studied recently, some are mentioned in the application section. Technical results culminating with the Ma-Trudinger-Wang condition [58] which gives necessary condition on $c$ for the regularity of the solution of the Optimal Mass Transportation problem. More recently attention has risen on multi marginal Optimal Mass Transportation [49] and has been systematically studied in [67] [70] [68] [69] . The data consists in an arbitrary (and even infinite) number $N$ of densities (the marginals) and the ground cost is defined on a product space $c(x_0,x_1,....,x_{n-1})$ of the same dimension. Several interesting applications belong to this class of models (see below).

Our focus is on numerical method in Optimal Mass Transportation and applications. The simplest way to build a numerical method is to consider sum of dirac masses ${\rho }_0={\sum }_{i=1}^N{\delta }_{A_i}{\rho }_1={\sum }_{j=1}^N{\delta }_{B_j}$. In that case the Optimal Mass Transportation problem reduces to combinatorial optimisation assigment problem between the points $\left\{A_i\right\}$s and $\left\{B_i\right\}$s : ${min}_{\sigma \in Permut(1,N)}\frac{1}{N}{\sum }_{i=1}^N{C}_{i,\sigma \left(i\right)}{C}_{i,j}={\parallel A_i-B_j\parallel }^2$. The complexity of the best (Hungarian or Auction) algorithm, see [21] for example, is $O\left(N^{\frac{5}{2}}\right)$. An interesting variant is obtained when only the target measure is discrete. For instance $X_0=\{\parallel x\parallel <1\},{\rho }_0=\frac{1}{|X_0|}{\rho }_1=\frac{1}{N}{\sum }_{j=1}^N{\delta }_{y_j}.$ It corresponds to the notion of Pogorelov solutions of the Monge-Ampère equation [71] and is also linked to Minkowski problem [18] . The optimal map is piecewise constant and the slopes are known. More precisely there exists $N$ polygonal cells $C_j$ such that $X_0={\cup }_j{C}_j$, $|C_j|=\frac{1}{N}$ and $\nabla {\phi }^{*}{|}_{C_j}=y_j$. Pogorelov proposed a constructive algorithm to build these solutions which has been refined and extended in particular in [39] [66] [63] [62] . The complexity is still not linear : $O\left(N^{2}logN\right)$.

For general densities data, the original optimization problem is not tractable because of the volume preserving constraint on the map. Kantorovitch dual formulation is a linear program but with a large number of constraints set over the product of the source and target space $X_0\times X_1$. The CFD formulation [2] . preserves the convexity of the objective function and transforms the volume preserving constraint into a linear continuity equation (using a change of variable). We obtained a convex but non smooth optimization problem solved using an Augmented Lagrangian method [43] , as originally proposed in [2] . It has been reinterpreted recently in the framework of proximal algorithm [64] . This approach is robust and versatile and has been reimplemented many times. It remains a first order optimization method and converges slowly. The cost is also increased by the additional artificial time dimension. An empirical complexity is $O\left(N^{3}LogN\right)$ where $N$ is the space discretization of the density. Several variants and extension of these methods have been implemented, in particular in [27] [17] . It is the only provably convergent method to compute Brenier (non $C^1$) solutions.

When interested in slightly more regular solutions which correspond to the assumption that the target support is convex, the recent wide stencil monotone finite difference scheme for the Monge-Ampère equation [45] can be adapted to the Optimal Mass Transportation problem. This is the topic of [6] . This approach is extremely fast as a Newton algorithm can be used to solve the discrete system. Numerical studies confirm this with a linear empirical complexity.

For other costs, JKO schemes, multi marginal extensions, partial transport ... efficient numerical methods are to be invented.