<figure><img src="image_1.png" /></figure>

HAL
open science

# Numerical resolution of Euler equations, through semi-discrete optimal transport

Jean-Marie Mirebeau

► To cite this version:

Jean-Marie Mirebeau. Numerical resolution of Euler equations, through semi-discrete optimal transport. 2015. hal-01237356

HAL Id: hal-01237356

https://hal.science/hal-01237356

Preprint submitted on 3 Dec 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical resolution of Euler equations,
through semi-discrete optimal transport

Jean-Marie Mirebeau*

December 3, 2015

Abstract

Geodesics along the group of volume preserving diffeomorphisms are solutions to Euler equations of inviscid incompressible fluids, as observed by Arnold [Arn66]. On the other hand, the projection onto volume preserving maps amounts to an optimal transport problem, as follows from the generalized polar decomposition of Brenier [Bre91].

We present, in the first section, the framework of semi-discrete optimal transport, initially developed for the study of generalized solutions to optimal transport [Abg09] and now regarded as an efficient approach to computational optimal transport. In a second and largely independent section, we present numerical approaches for Euler equations seen as a boundary value problem [Bre08, BCC$^{+}$15, MM15]: knowing the initial and final positions of some fluid particles, reconstruct intermediate fluid states. Depending on the data, we either recover a classical solution to Euler equations, or a generalized flow [Bre93] for which the fluid particles motion is non-deterministic, as predicted by [Shn94].

# 1 Semi-discrete optimal transport

The optimal transport problem, formulated by Monge in 1781, has been the subject of continued
research ever since. The objective is find a map $T$ pushing a distribution of mass $\mu$ onto another
one $\nu$, and which minimizes the average transport cost, see [Vil08] for an overview. We present
in this section the framework of *semi-discrete optimal transport*, in which one of the measures
$\nu$ is *discrete*, consisting of a finite weighted sum of Dirac masses, and the second measure $\mu$ is
absolutely *continuous* with respect to the Lebesgue measure.

Semi-discrete optimal transport was introduced by Pogorelov and Alexandrov as a theoretical
tool for the study of generalized solutions to optimal transport [Abg09], which are not solutions
to the Monge-Ampere equation in either the classical or the viscosity sense. Its numerical imple-
mentation yields some of the most reliable, fast and scalable methods for computational optimal
transport [AHA98, Mer11, Lév14].

We give in §1.1 a brief overview of the alternative computational approaches to optimal
transport, and proceed in the following subsections §1.2, §1.3, to the description of the semi-
discrete method. A prerequisite is to define the mathematical counterpart of (possibly non-
optimal) mass transport.

<footnote>*Laboratoire de mathématiques d'Orsay, Université Paris-Sud, CNRS, Université Paris-Saclay, 91405 Orsay, France</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: A map $T$ sends domain $X$ (blue boundary), equipped with a density $\mu = \rho \text{Leb}_X$ (grayscale) onto a finite set $Y$ (red dots). The push-forward $\nu = T_{\#}\mu$ is defined by $\nu(\{y\}) = \mu(T^{-1}(\{y\}))$ (mass of the green region), for all $y \in Y$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: The three paradigms of numerical optimal transport. I: Fully discrete, II: Fully continuous, III: Semi-Discrete. All illustrations of this section courtesy of Q. Merigot.</figcaption></figure>

**Definition 1.1** (Push forward of a measure $\mu$ by a map $T$). Let $X, Y$ be metric spaces. Let $\mu$ be a Borel measure on $X$, and $T: X \to Y$ be a Borel measurable map. The push-forward of $\mu$ by $T$ is the Borel measure $\nu = T_{\#}\mu$ on $Y$ defined for any Borel subset $A \subseteq Y$ by

$$ \nu(A) := \mu(T^{-1}(A)). $$

Mass transport preserves total mass by construction: $\nu(Y) = \mu(T^{-1}(Y)) = \mu(X)$ if $\nu = T_{\#}\mu$, in particular $\mu \in \text{Prob}(X) \Rightarrow T_{\#}\mu \in \text{Prob}(Y)$ where $\text{Prob}$ denotes the collection of Borel probability measures. Note also that $T_{\#}\mu$ depends linearly on $\mu$, but non-linearly on $T$. In the following, $X, Y$ denote metric spaces which closed balls are compact. This allows to equip $\text{Prob}(X)$ with the weak* topology dual to continuous bounded functions $C_b^0(X)$. All subsets, maps and measures are implicitly assumed to be Borelian in the following.

The optimal transport problem, as originally stated by Monge, involves two probability measures $\mu \in \text{Prob}(X)$, $\nu \in \text{Prob}(Y)$, and a lower-semi-continuous transport cost $c: X \times Y \to \mathbb{R}_+ \cup \{+\infty\}$. The objective is to find a map $T: X \to Y$ pushing $\mu$ onto $\nu$ and minimizing the average transport cost:

$$ \min_{T_{\#}\mu=\nu} \int_X c(x, T(x)) d\mu(x). \qquad (1) $$

Under this form, the optimal transport problem may have no solution, or several solutions. Adequate problem relaxations, or structural assumptions on the cost and the measures, make the problem well posed. They also allow for numerical approaches, see the next subsection.

## 1.1 Computational optimal transport

The main inputs to the optimal transport problem are the measures $\mu \in \text{Prob}(X)$ and $\nu \in \text{Prob}(Y)$. We choose to categorize numerical approaches to optimal transport by the structure of this input data, which always belongs to one of the following two categories:

* Discrete measures, given as a finite sum $\mu = \sum_{x \in X} \mu_x \delta_x$ of Dirac masses.

* Absolutely continuous measures $\mu = \rho \text{Leb}_X$ with respect to Lebesgue on $X \subseteq \mathbb{R}^d$.
Discrete and absolutely continuous probabilities both form dense affine subspaces of Prob($X$), so there is in principle no loss of generality in assuming either form. Surprisingly enough, this input form has deep implications on the algorithm structure, the numerical cost, and the generality of the (weak) solutions to the optimal transport problem that can be extracted. A method is said fully discrete (resp. fully continuous) if both input probabilities are discrete (resp. continuous), and semi-discrete otherwise.

I Fully discrete optimal transport, via the linear relaxation.

Kantorowitch [Kan42] introduced a convex and linear relaxation of the optimal transport problem (1), for which the existence of a minimizer is guaranteed under the assumptions of §1 on the spaces $X, Y$ and the cost $c$. The relaxed transport cost from $\mu \in \text{Prob}(X)$ to $\nu \in \text{Prob}(Y)$ is

$$W(\mu, \nu) := \min\left\{\int_{X \times Y} c(x, y)d\pi(x, y); \pi \in \text{Prob}(X \times Y), P_{\#}^{X}\pi = \mu, P_{\#}^{Y}\pi = \nu\right\}, \quad (2)$$

where $P^X, P^Y$ denote projections onto first and second coordinates respectively. Monge's original formulation (1) corresponds to couplings of the special form $\pi = (\text{Id}, T)_{\#}\mu$, where $T$ is the unknown transport map. Kantorowitch's formulation is especially needed when both measures $\mu, \nu$ are discrete, due to the lack of measure preserving<fnref n="1" /> maps $T: (X, \mu) \to (Y, \nu)$.

When $\mu$ and $\nu$ both consist of $N$ weighted Dirac masses, Kantorowitch's formulation (2) amounts to minimizing a linear objective depending on $N^2$ non-negative variables (characterizing the coupling $\pi$), subject to $2N - 1$ linear constraints. General purpose linear program solvers such as the simplex method, or primal dual interior point methods, are directly applicable. Linear program solvers specialized in this specific instance also exist, such as the Hungarian algorithm, and Bertseka's auction's algorithm.

The performance of these methods is however disappointing on large problem instances $N \ge 5000$ Dirac masses due to the quadratic number of optimization variables. Since the support of the optimal $\pi$ is usually a small subset of $X \times Y$, concentrated along the graph of a deterministic map $T$, multiscale strategies have been proposed [OR15, Sch15] to eliminate inactive variables. Another approach [BCC$^+$15] consists in adding a small multiple of the entropy of $\pi$ with respect to $\mu \otimes \nu$ to the linear objective in (2), which makes the problem smooth and strictly convex.

II Fully continuous optimal transport, via the Monge-Ampere PDE formulation.

The PDE approach to optimal transport requires some structural assumptions: $X, Y$ are convex and compacts domains of $\mathbb{R}^d$, the cost $c(x, y) = \frac{1}{2}\|x - y\|^2$ is quadratic, and the measures $\mu = \rho \text{Leb}_X$, $\nu = \sigma \text{Leb}_Y$ have continuous densities bounded above and below. Monge's problem (1) then admits a unique solution, of the form $T = \nabla u$ where $u$ solves the following Monge-Ampere equation in the sense of viscosity solutions:

$$\left\{ \begin{array}{l} u: X \to \mathbb{R}, \text{ convex,} \\ \nabla u(X) = Y, \\ \det(\nabla^2 u(x)) = f(x)/g(\nabla u(x)), \quad \forall x \in \Omega. \end{array} \right. \qquad (3)$$

Each line of this formulation raises challenging discretization problems. First, the constraint of convexity can be discretized in various ways, none of which is particularly simple

<footnote n="1">¹Except in the special case where these measures are equidistributed on an equal number of points.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Left: Domain $X$ (blue boundary), population density $\rho$ (grayscale), and bakeries $Y$ (red dots). Center: Voronoi regions (red boundaries) attached to the bakeries. Right: Laguerre diagram, depending on both the bakeries positions $Y$ and their prices $\psi: Y \to \mathbb{R}$.</figcaption></figure>

or canonical [AM09, MO14, Mir15]. For the problem of interest, convexity can also be imposed through the discretization of the Monge-Ampere operator [FO11, BCM15]. Second, discretizations of the boundary condition $\nabla u(X) = Y$ have only been proposed recently [BFO14].

Finally, the Monge-Ampere operator $\det(\nabla^2 u)$ can be discretized through standard finite differences [LR05, BFO10]. However, it is desirable to preserve its monotony at the discrete level, in the spirit of [KT92], which requires more complex implementations [FO11, BCM15].

Other approaches to optimal transport have been proposed, such as the fluid dynamics formulation of Benamou and Brenier [BBG02]. The next section introduces the semi-discrete approach.

## 1.2 An economic metaphor

We illustrate semi-discrete optimal transport through an economic metaphor. The source domain $X \subseteq \mathbb{R}^d$, connected and bounded, is regarded as a large city on which population is distributed according to a density $\mu = \rho \text{Leb}_X$. The target domain $Y \subseteq \mathbb{R}^d$ is finite, and is regarded as the location of bakeries in the city, see Figure 3 (left). The function $c: X \times Y \to \mathbb{R}_+ \cup \{+\infty\}$ represents the cost for agent $x \in X$ to walk to bakery $y \in Y$.

It is natural for potential customers to walk to the closest bakery. The basin of attraction of bakery $y \in Y$ is called its Voronoi region, see Figure 3 (center):

$$ \text{Vor}(y) := \{x \in X; \forall z \in Y, c(x, y) \le c(x, z)\}. $$

However, if the bakeries are allowed to set different prices $\psi: Y \to \mathbb{R}$ for their products, then people $x \in X$ will make a compromise and minimize the sum $c(x, y) + \psi(y)$ of transport cost and purchase price. In this refined model, each bakery $y \in Y$ attracts a members of its Laguerre region, see Figure 3 (center)

$$ \text{Lag}_{\psi}(y) := \{x \in X; \forall z \in Y, c(x, y) + \psi(y) \le c(x, z) + \psi(z)\}. \quad (4) $$

When the cost is quadratic $c(x, y) := \frac{1}{2}|x - y|^2$, the Voronoi and Laguerre regions are convex polytopes. Their computation is a non-trivial yet classical problem of discrete geometry, for which mature software packages are available such as [cga].

The destination of each customer $x \in X$ is given by the price dependent transport map

$$ T_{\psi}(x) := \underset{y \in Y}{\operatorname{argmin}} c(x, y) + \psi(y). \quad (5) $$
Agents lying in the intersection of two distinct Laguerre regions $Lag_ψ(y) \cap Lag_ψ(z)$ do not have a well defined destination. We thus assume in the following, including Th 1.3, that these interesections are Lebesgue-negligible, hence also $\mu$-negligible, for any price function $\psi: Y \to \mathbb{R}$ and any distinct $y, z \in Y$. This property holds under mild non-degeneracy assumptions on the cost $c$, and is satisfied by the quadratic cost $c(x, y) := \frac{1}{2}|x - y|^2$.

The frequentation of the bakeries is encoded in the push-forward

$$T_{\psi\#\mu} = \sum_{y \in Y} \mu(Lag_{\psi}(y))\delta_y.$$

Customers $x \in X$, by their individual cost minimizations, construct an optimal transport plan.

**Lemma 1.2.** For any $\psi: Y \to \mathbb{R}$, the map $T_ψ$ is an optimal transport from $\mu$ to $T_{\psi\#\mu}$.

*Proof.* Let $T: X \to Y$ such that $T_{\#\mu} = T_{\psi\#\mu}$. The price dependent transport plan $T_ψ$ has the optimality property:

$$\forall x \in X, \quad c(x, T_ψ(x)) + \psi(T_ψ(x)) \le c(x, T(x)) + \psi(T(x)).$$

Integrating with respect to $\mu$ yields

$$\int_X c(x, T_ψ(x))d\mu(x) + \int_X \psi(T_ψ(x))d\mu(x) \le \int_X c(x, T(x))d\mu(x) + \int_X \psi(T(x))d\mu(x). \quad (6)$$

Denoting by $\nu \in \text{Prob}(Y)$ the common value of $T_{\#\mu} = T_{\psi\#\mu}$, we obtain

$$\int_X \psi(T_\psi(x)) d\mu(x) = \int_Y \psi d\nu = \int_X \psi(T) d\mu(x). \quad (7)$$

Subtracting from both sides of (6) the equal terms (7), namely the total revenue of the bakeries, yields as announced

$$\int_X c(x, T_\psi(x)) d\mu(x) \le \int_X c(x, T(x)) d\mu(x). \quad \square$$

## 1.3 The Kantorowitch functional

We continue with the economical metaphor of §1.2, and consider the situation where each bakery $y \in Y$ only has the capacity to cater for proportion $\nu_y$ of the population, and $\sum_{y \in Y} \nu_y = 1$. It turns out that it is possible to fix the prices $\psi$ so that the bakery frequentations $T_{\psi\#\mu}$ match their capacities $\nu$. Doing so amounts to an optimal transport problem.

**Theorem 1.3** (Aurenhammer, Hoffman, Aranov [AHA98]). Let $\mu = \rho \text{Leb}_X$ be a probability with density on the bounded domain $X \subseteq \mathbb{R}^d$, and let $\nu = \sum_{y \in Y} \nu_y \delta_y$ be a probability on the finite set $Y$. The following are equivalent:

(a) Finding an optimal transport from $\mu$ to $\nu$.

(b) Finding prices $\psi$ on $Y$ such that $T_{\psi\#\mu} = \nu$.

(c) Maximizing the concave function

$$\Phi(\psi) := \sum_{y \in Y} \int_{\text{Lag}_\psi(y)} (c(x, y) + \psi(y)) d\mu(x) - \sum_{y \in Y} \psi(y) \nu_y. \quad (8)$$
Interestingly, Monge's formulation (1) regards mass transport $T_{\#}\mu = \nu$ as a constraint and then optimizes the cost, whereas point (b) of Theorem 1.3 suggests exploring a family of cost-optimal mappings $T_{\psi}$ until the image measure matches $\nu$. The principles underlying these formulations of optimal transport seem opposite, and indeed Theorem 1.3 is a by-product of the dual formulation of the linear relaxation (2) of optimal transport, identified by Kantorowitch:

$$W(\mu, \nu) = \sup\{\int_X \varphi d\mu - \int_Y \psi d\nu; \varphi, \psi \text{ obey } \forall (x, y) \in X \times Y, \varphi(x) - \psi(y) \le c(x, y)\}, \quad (9)$$

where implicitly we also require $\varphi \in L^1(\mu)$ and $\psi \in L^1(\nu)$. If $Y$ is finite, then (9) amounts to maximize over the bakery prices $\psi: Y \to \mathbb{R}$ the following quantity, actually equal to (8)

$$\int_X \varphi_{\psi} d\mu - \sum_{y \in Y} \psi(y) \nu_y \quad \text{where } \forall x \in X, \varphi_{\psi}(x) := \min_{y \in Y} c(x, y) + \psi(y). \quad (10)$$

Note that $\varphi_{\psi}$ is the optimal choice of $\varphi$ in (9) for a given $\psi$.

The semi-discrete approach to numerical optimal transport is nothing else than the maximization of the Kantorowitch functional (8). Coordinate-wise increments were the first introduced optimization method [OP89], and have cubic complexity $O(N^3(\ln N)/\epsilon)$ where $N := \#(Y)$ and $\epsilon$ is the desired precision. Gradient based methods such as LBFGS are much more efficient in practice [Mer11, Lév14], and are themselves outperformed by Newton methods [dGBOD12]. Their implementation requires to identify the first and second derivatives of the Kantorowitch functional. Differentiation with respect to the price $\psi(y)$ at bakery $y \in Y$ is denoted by $\partial_y$ or $\partial/\partial_y$.

For any prices $\psi: Y \to \mathbb{R}$, and for each point $x \in X$ not at the intersection of two Laguerre cells, one has

$$\partial_y \varphi_{\psi}(x) = \begin{cases} 1 & \text{if } x \in \text{Lag}_{\psi}(y), \\ 0 & \text{otherwise,} \end{cases} \qquad \text{hence } \partial_y \Phi(\psi) = \mu(\text{Lag}_{\psi}(y)) - \nu_y. \quad (11)$$

Maximizing $\Phi$ amounts to solving $\nabla\Phi = 0$, which is thus equivalent to the discrete Monge-Ampere equation $\forall y \in Y, \mu(\text{Lag}_{\psi}(y)) = \nu_y$. The expressions of the second derivatives, first identified in [dGBOD12], involve integrating the density $\rho$ over Laguerre cell boundaries. Assuming a quadratic cost $c(x, y) = \frac{1}{2}|x - y|^2$, the Kantorowitch functional is twice differentiable, and for all distinct $y, z \in Y$

$$\frac{\partial^2 \Phi}{\partial y \partial z} = \int_{\text{Lag}_{\psi}(y,z)} \frac{\rho(s)}{\|y-z\|} ds, \qquad \frac{\partial^2 \Phi}{\partial y^2} = -\sum_{z \neq y} \frac{\partial^2 \Phi}{\partial y \partial z}, \quad (12)$$

where $\text{Lag}_{\psi}(y, z) := \text{Lag}_{\psi}(y) \cap \text{Lag}_{\psi}(z)$ and $ds$ denotes $(d-1)$-dimensional Lebesgue measure. As a result, $\nabla^2\Phi$ is the matrix of a Laplacian, on the graph dual to the Laguerre Diagram: for $v: Y \to \mathbb{R}$

$$v^{\mathrm{T}}(\nabla^{2}\Phi)v = -\sum_{y \neq z}(\partial_{yz}^{2}\Phi)(v(y) - v(z))^{2}. \quad (13)$$

This matrix is therefore negative semi-definite. In addition, $\mathrm{Ker}(\nabla^2\Phi)$ is the space of constant functions on $Y$, provided $\mu(\text{Lag}_{\psi}(y)) > 0$ for all $y \in Y$, and one has the continuity of $\rho$ and connectedness of $\{x \in X; \rho(x) > 0\}$. Our numerical implementation consists of maximizing the Kantorowitch functional using a damped Newton algorithm. The damping, a reduction in the step size, is chosen so as to enforce $\nu(\text{Lag}_{\psi}(y)) > 0$ at each iteration, hence also the invertibility of the Kantorowitch Hessian (13) up to constant functions.
<figure><img src="image_6.png" /><figcaption>Figure 4: Left: Laguerre cells and their dual graph. Center: Numerical test case setup. Right: Computation time versus number of Dirac masses.</figcaption></figure>

The performance of the method is evaluated by computing the optimal transport from the uniform density on the hollowed square $X = [0, 3]^2 \setminus [1, 2]^2$, to a cartesian grid $Y$ of Diracs approximating the the uniform density on the square $[0, 1]^2$. Note that the non-convexity of the domain $X$ would be a source of difficulty for PDE formulations of optimal transport (3).

Numerical experiments show that the computation time scales linearly in the cardinality of $Y$, see Figure 4, and that instances involving hundreds of thousands of Dirac masses can be solved in a few minutes on a modest 2.5Ghz Laptop using a single core. The semi-discrete approach to optimal transport has been extended to transport costs other than quadratic, in particular those subject to the Ma-Trudinger-Wang property [MTW05], and applied to inverse problems arising in geometric optics [DCMT14].

## 2 Numerical resolution of Euler equations

Euler discovered in 1755 the equations, which bear his name, governing the motion of incompressible fluid flows without viscosity. His mathematical notation, extremely modern at the time, closely resembles today's form: on a domain $X \subseteq \mathbb{R}^d$, and an interval $I$ of time

$$ \partial_t u + (u \cdot \nabla)u = -\nabla p, \qquad \mathrm{div} \, u = 0, \qquad (14) $$

where $p: I \times X \to \mathbb{R}$ denotes pressure and $u: I \times X \to \mathbb{R}^d$ is the velocity field, which is subject to a slipping condition $u \cdot n = 0$ on the boundary $\partial X$. Despite their age, the mathematical analysis of Euler equations is still largely open, see [FD12] for an overview.

Let us introduce the flow $s: I \times X \to X$ of the velocity field $u$, which describes the motion of fluid particles in Lagrangian coordinates. It is defined for any time $t \ge 0$ by an Ordinary Differential Equation (ODE), parametrized by the particle initial position $x \in X$

$$ s(0, x) = x, \qquad \partial_t s(t, x) = u(t, s(t, x)). $$

Depending on the occurrence, we often denote $s_t(x) := s(t, x)$. A fundamental observation of Arnold [Arn66] is that the Lagrangian coordinate system $(s_t)_{t \in I}$ formally obeys the equation of geodesics on the infinite dimensional space $\mathcal{SDiff}$ of diffeomorphisms of $X$ with unit jacobian, equipped with the $L^2(X)$ metric, see §2.1. This leads to two natural problems for Euler equations:

* The Cauchy problem, forward in time: given the initial position $s_0$ and the initial velocity $\dot{s}_0$ of the fluid particles, find their positions $(s_t)_{t \ge 0}$ at all positive times. This amounts to computing the exponential map on the space $\mathcal{SDiff}$.
<figure><img src="image_7.png" /><figcaption>Figure 5: The motion of inviscid incompressible fluids admits three formulations, either (I) Eulerian based on the local speed $v: [0, 1] \times X \to \mathbb{R}^d$, (II) Lagrangian based on diffeomorphisms $s(t, \cdot)$ which integrate the speed: $\partial_t s(t, x) = v(t, s(t, x))$, or (III) relaxed as a superposition of individual particles paths $\omega \in \Omega$, weighted by a measure $\mu$.</figcaption></figure>

• The boundary value problem: given the initial $s_0$ and final $s_1$ positions of the fluid particles, find their intermediate states $(s_t)_{0 \le t \le 1}$. This amounts to computing a minimizing geodesic on the space $\mathcal{SDiff}$.

This manuscript is devoted to the boundary value problem. Regarding the Cauchy problem, we refer for instance to the recent works of Székelyhidi et al [DLSJ09], on which a course was given at the Journées EDP 2015. Let us mention that the tools of semi-discrete optimal transport can also be used for the Cauchy problem [dGWH$^{+}$15].

We give in §2.1 more detail on Arnold’s formal interpretation of Euler equations as geodesics. Then we describe several relaxations in §2.2, 2.3, and the resulting discretizations including the one §2.4 proposed by Merigot and the author [MM15] and involving semi-discrete optimal transport.

## 2.1 Euler equations as a geodesic equation

In addition to incompressible fluid flows, Euler discovered in 1765 the equations governing the motion of a rigid body [Eul65]. Two centuries later, Arnold unified these two works [Arn66], identifying Eulerian motions of a fluid and of a rigid object as the geodesics on the group of its configurations, equipped with a suitable left-invariant metric. For a rigid object this is the group $SO_3$ of spatial rotations with respect to some reference position.

The configurations of a fluid on a domain $X$ are on the other hand described by maps $s: X \to X$, where $s(x) \in X$ is the present position of the fluid particle which reference position is $x \in X$. Smooth configurations are described by diffeomorphisms. Fluid incompressibility implies that the Jacobian $\det(\nabla s)$ is unit. In the following $X \subseteq \mathbb{R}^d$ is a connected bounded domain of unit area, and we let

$$ \mathcal{SDiff} := \{s \in C^{\infty}(X, X); \det(\nabla S) = 1 \text{ on } X\}. $$

Arnold [Arn66] embeds $\mathcal{SDiff}$ within the Hilbert space $\mathcal{M} := L^2(X, \mathbb{R}^d)$, which norm is denoted $\| \cdot \|$, and regards it as a submanifold. The energy and length of a path $s \in H^1([0, 1], \mathcal{M})$ are respectively

$$ \mathcal{E}(s) := \int_{0}^{1} \|\dot{s}_{t}\|^{2} dt, \quad \mathcal{L}(s) := \int_{0}^{1} \|\dot{s}_{t}\| dt, \quad \text{where } \|m\|^{2} := \int_{X} |m(x)|^{2} dx. \quad (15) $$
Here and below, $\dot{s}_t$ denotes the time derivative of $s_t$, hence $\|\dot{s}_t\|_{\mathbb{M}}^2 = \int_X |\dot{s}_t(x)|^2 dx$ is the kinetic energy of the fluid particles at time $t$. The geodesic distance between $s_*, s^* \in \mathbb{SDiff}$ is defined as

$$d(s_*, s^*)^2 := \inf\{\mathcal{E}(s); s \in H^1([0, 1], \mathbb{SDiff}), s_0 = s_*, s_1 = s^*\}. \quad (16)$$

By Cauchy-Schwartz's inequality, $\mathcal{E}(s) \ge \mathcal{L}(s)^2$, with equality if and only $s$ is parametrized at constant speed. The minimal geodesic for (16), if it exists, thus also minimizes path length, and has constant speed, which corresponds to the fluid kinetic energy conservation.

Following Arnold we derive below the equations of geodesics along $\mathbb{SDiff}$, only formally since $\mathbb{SDiff}$ is not closed in $L^2(X, \mathbb{R}^d)$. The tangent space to $\mathbb{SDiff}$ at the identity element $\mathrm{Id}$, consists of vector fields $v \in C^\infty(X, \mathbb{R}^d)$ such that $\mathrm{Id} + tv$ has unit jacobian up to second order terms as $t \to 0$.

$$\det(\mathrm{Id} + t\nabla v) = 1 + t \operatorname{Tr}(\nabla v) + \mathcal{O}(t^2) = 1 + t \operatorname{div} v + \mathcal{O}(t^2).$$

Hence $T_{\mathrm{Id}}\mathbb{SDiff} = \{v \in C^\infty(X, \mathbb{R}^d); \operatorname{div} v = 0\}$. From the group law we obtain for any $s \in \mathbb{SDiff}$

$$T_s\mathbb{SDiff} = s \cdot T_{\mathrm{Id}}\mathbb{SDiff} = \{v \circ s; v \in C^\infty(X, \mathbb{R}^d), \operatorname{div} v = 0\}.$$

Recall that a smooth path on a submanifold of a Riemannian manifold is a geodesic iff its acceleration is at each time normal to the submanifold tangent space. A geodesic path $(s_t)_{0 \le t \le 1}$ in $\mathbb{SDiff}$ thus satisfies for any divergence free vector field $v$ and any $t \in [0, 1]$

$$0 = \langle \ddot{s}_t, v \circ s_t \rangle = \langle \ddot{s}_t \circ s_t^{-1}, v \rangle,$$

where $\langle \cdot, \cdot \rangle$ denotes the $L^2$ scalar product. For the second identity, we used the fact that $s_t$ preserves the Lebesgue measure on $X$, hence that right composition with $s_t$ is an isometry of $\mathbb{M} := L^2(X, \mathbb{R}^d)$. On a simply connected domain such as $X$, a vector field is orthogonal to all divergence free vector fields if and only it is a gradient field. Hence there exists a pressure field $p_t \in C^\infty(X, \mathbb{R})$, for all $t \in [0, 1]$, such that

$$\ddot{s}_t = -\nabla p_t \circ s_t, \quad \text{i.e.} \quad \partial_t^2 s(t, x) = -\nabla p(t, s(t, x)) \quad (17)$$

Introducing the Eulerian velocity field, defined by $u(t, s(t, x)) = \partial_t s(t, x)$, we obtain

$$\begin{aligned} \partial_t^2 s(t, x) &= \frac{d}{dt}u(t, s(t, x)) = \partial_t u(t, s(t, x)) + \partial_t s(t, x)\nabla u(t, s(t, x)) \\ &= \partial_t u(t, s(t, x)) + \partial_t u(t, s(t, x))\nabla u(t, s(t, x)). \end{aligned}$$

Identifying these two expressions of the fluid particles acceleration $\partial_t^2 s(t, x)$, and composing with $s_t^{-1}$, we obtain Euler equations (12) as announced.

The main defect of this approach is that $\mathbb{SDiff}$ is not a submanifold of $\mathbb{M} := L^2(X, \mathbb{R}^d)$, in particular because it is not closed, hence geodesics need not exist. A possible remedy is to consider instead groups of diffeomorphisms with Sobolev regularity of index $s > d/2 + 1$, equipped with the associated Sobolev metric. The exponential map is then well defined and surjective, and minimizing geodesics exist between any pair of points [CP92, BV14]. Of course, geodesic equations with respect to these stronger metrics differ from Euler's equations, and from any known physical model. We present in the next section an opposite approach, which relaxes the set of solutions, rather than strengthening the metric.
## 2.2 The $L^2$-relaxation, and flows as permutations

Diffeomorphisms of unit jacobian admit, as a superset, the collection $\mathcal{S}$ of maps preserving the Lebesgue measure on $X$

$$ \mathcal{S} := \{s \in L^2(X, \mathbb{R}^d); s_{\#} \mathrm{Leb}_X = \mathrm{Leb}_X\}, $$

where $\mathrm{Leb}_X$ denotes the Lebesgue measure on the domain $X$. See Definition 1.1 for the push-forward of a measure by a map. The set $\mathcal{S}$ is closed in $L^2(X, \mathbb{R}^d)$, and is the completion of $\mathcal{SDiff}$ in dimension $d \ge 3$. In dimension $d \le 2$ however, rigidity phenomena prevent the density of $\mathcal{SDiff}$ in $\mathcal{S}$. For instance $\mathcal{SDiff} = \{\mathrm{Id}\}$ if $X = [0, 1]$, whereas $\mathcal{S}$ contains $\mathbb{1}$ the symmetric group of $\{1, \dots, N\}$ for all $N$, see below. In dimension $d \in \{1, 2\}$ the use of measure preserving maps can nevertheless be defended by arguing [Bre08] that fluid flows on $d$-dimensional domains $X$ generally arise by neglecting small dimensions in a three dimensional physical domain $X \times [0, \varepsilon]^{3-d}$.

**Permutations as measure preserving maps.** Assume that the domain $X$ can, for all $N \ge 1$, be partitioned into $N$ cubes, identical up to translation. For instance $[0, 1] = \cup_{k=0}^{N-1} [k/N, (k+1)/N]$. Denote by $X_N$ the collection of their centers, equipped with the uniform probability measure. Then the collection $\Sigma_N$ of permutations of $X_N$ isometrically embeds into $\mathcal{S}$, by permuting the partitioning cubes with their centers, and as $N \to \infty$ it becomes dense in $\mathcal{S}$.

A discretization of Euler's boundary value problem (16) based these permutations is proposed in [Bre89a]. Given a number $T$ of time steps, and $s_0, s_T \in \Sigma_N$ approximating some boundary data $s_*, s^* \in \mathcal{S}$, minimize over $(s_t)_{t=1}^{T-1}$ the path energy (15) counterpart:

$$ \frac{T}{N} \sum_{0 \le t < T} \sum_{x \in X_N} |s_{t+1}(x) - s_t(x)|^2. \qquad (18) $$

Figure 6 reproduces numerical experiments of [Bre08] for the one dimensional case $X = [0, 1]$, with boundary conditions $s_* = \mathrm{Id}$, $s^* = 1 - \mathrm{Id}$. The permutations $s_t$ associated to intermediate times $0 < t < T$ turn out to oscillate rapidly. They do not converge in $\mathcal{S}$, but do so in an even more relaxed setting allowing particles to split and be distributed over the domain $X$ at intermediate times, see the next subsection.

## 2.3 The linear relaxation, and the pressure field

Brenier [Bre89b], introduced a second relaxation of Euler equations, which is linear, and for which the existence of a minimizer is guaranteed. It involves the space of all possible particles paths

$$ \Omega := X^{[0,1]}, $$

which is compact when equipped with Tychonoff's product topology. Probability measures $\mu \in \mathrm{Prob}(\Omega)$ are called generalized flows. They are said incompressible if the particle positions equidistribute at all times on the domain: $\forall t \in [0, 1]$

$$ e_{t\#\mu} = \mathrm{Leb}_X, \qquad (19) $$

where $e_t: \Omega \to X: \omega \mapsto \omega(t)$ denotes the evaluation map at time $t$. The boundary conditions are also modeled by a marginal constraint, coupling the initial and end time:

$$ (e_0, e_1)_{\#\mu} = (s_*, s^*)_{\#\mathrm{Leb}_X}. \qquad (20) $$
<figure><img src="image_8.png" /><figcaption>Figure 6: Relaxed Euler boundary value problem on $X = [0, 1]$ with boundary conditions $s_* = \text{Id}$, $s^* = -\text{Id}$. The solution for this particular instance is unique and known [BFS09], it involves non-deterministic generalized flows see §2.3. Horizontal axis: initial particle position $x \in [0, 1]$. Vertical axis: particle position in $[0, 1]$ at some specified time $t$. Left: discretization based on permutations [Bre08] see (18). Right: convex relaxation as multi-marginal problem [BCC$^{+}$15], regularized via entropy penalization, see (23).</figcaption></figure>

Finally, the action is the mean with respect to $\mu$ of the kinetic energy of the paths

$$ \mathcal{E}(\mu) := \int_{\Omega} \int_{0}^{1} |\dot{\omega}(t)|^{2} dt d\mu(\omega). \qquad (21) $$

In this form, Euler's boundary value problem becomes the optimization of a lower-semi-continuous<fnref n="2" /> linear function $\mathcal{E}(\mu)$, subject to infinitely many marginal constraints(19), (20), which are linear and continuous. This optimization problem may have several solutions. Strikingly however, the dual linear problem has a unique<fnref n="3" /> solution, called the pressure field $p: [0, 1] \times X \to \mathbb{R}$, which accounts for the incompressibility constraint. For any minimizer $\mu \in \text{Prob}(\Omega)$ of the primal problem, $\mu$-almost every trajectory obeys $\ddot{\omega}(t) = -\nabla p(t, \omega(t))$, an equation which makes sense thanks to the regularity $p \in L_{\text{loc}}^{2}([0, 1[, \text{BV}(X))$ see [AF07].

**Discretization as a multi-marginal, entropy penalized problem.** Let $N, T$ be positive integers, let $X_N$ be a sampling of $N$ points equidistributed on the domain $X$, and let $s: X_N \to X_N$ be a permutation representing an incompressible map of interest. Benamou et al [BCC$^{+}$15] discretize the linear relaxation of Euler equations through the following multi-marginal optimization problem, which is strictly convex thanks to the added $\varepsilon$-entropy penalization: minimize over all generalized flows $\gamma \in \text{Prob}(X_N^T)$

$$ \sum_{x_1, \dots, x_T \in X_N} \gamma(x_1, \dots, x_T) c(x_1, \dots, x_T) + \varepsilon \gamma(x_1, \dots, x_T) (\ln \gamma(x_1, \dots, x_T) - 1), \\ \text{where } c(x_1, \dots, x_T) := n(|x_1 - x_2|^2 + \dots + |x_{T-1} - x_T|^2 + |x_T - s(x_1)|^2), \qquad (22) $$

subject to the marginal constraints $\sum_{x_1, \dots, x_T} \gamma(x_1, \dots, \hat{x}_i, \dots, x_T) = 1$ for any $1 \le i \le N$, and any fixed $\hat{x}_i \in X_N$, discretizing (20). The large number of variables $\#(X_N^T) = N^T$ makes a direct approach hopeless. However, the dual problem only involves $NT$ variables, instead of $N^T$,

<footnote n="2">²Because path energy $\omega \mapsto \int_0^1 ||\dot{\omega}(t)||^2 dt$ is lower-semi continuous on $\Omega$</footnote>

<footnote n="3">³Unique up to the trivial invariance by addition of a function depending only on time.</footnote>
associated to the incompressibility constraints and thus representing the pressure field: maximize over all $p_1, \cdots, p_T : X_N \to \mathbb{R}$

$$ \sum_{x \in X_N} \sum_{1 \le t \le T} p_t(x) - \sum_{x_1, \cdots, x_T} \varepsilon \exp((p_1(x_1) + \cdots + p_T(x_T) - c(x_1, \cdots, x_T))/\varepsilon). \quad (23) $$

This functional, which involves a sum over $N^T$ terms, can be evaluated using only $T$ products of $N \times N$ matrices, thanks to the additive structure of the cost (22) and the multiplicative properties of the exponential. Figure 6 reproduces numerical experiments of [BCC$^+$15]. The drawbacks of this approach are its numerical cost, as well as the blurring induced by the entropy penalization, which limits the visibility of the transition from classical solutions of (14) to non-deterministic generalized flows minimizing (21).

## 2.4 Polar factorization, and the distance to incompressible maps

We describe a discretization of Euler's boundary value problem, introduced by Merigot and the author in [MM15], which implements flow incompressibility by penalizing the distance to incompressible maps: for any $m \in \mathbb{M} := L^2(X, \mathbb{R}^d)$

$$ d^2(m, \mathbb{S}) := \inf_{s \in \mathbb{S}} \|m - s\|^2. \quad (24) $$

The distance from $m$ to $\mathbb{S}$ turns out to be the Wasserstein distance from the image measure of $m$ to the Lebesgue measure. This surprisingly links Euler equations of incompressible fluid flows, to optimal transport which describes a pressure-less gaz [BBG02].

**Theorem 2.1** (Brenier's polar factorization [Bre91]). For any $m \in L^2(X, \mathbb{R}^d)$, one has

$$ d(m, \mathbb{S}) = W_2(s_{\#} \mathrm{Leb}_X, \mathrm{Leb}_X), \quad (25) $$

where $W_2$ is the Wasserstein distance<fnref n="4" /> for the quadratic cost $c(x, y) := |x - y|^2$. Furthermore, if $m_{\#} \mathrm{Leb}_X$ is absolutely continuous with respect to the Lebesgue measure, then there exists a unique convex function $\varphi : X \to \mathbb{R}$, and a unique measure preserving map $s \in \mathbb{S}$, such that

$$ m = \nabla \varphi \circ s. \quad (26) $$

Brenier's polar factorization takes its name from the special instance where $X$ is the euclidean unit ball and $m : x \mapsto Mx$ is a linear function. Let $M = \Phi S$ be the polar decomposition of $M \in \mathrm{GL}_d \mathbb{R}$, into the product of a a positive definite matrix $\Phi \in S_d^+$, and an orthogonal transformation $S \in \mathcal{O}_n(\mathbb{R})$. Then $\varphi : x \mapsto \frac{1}{2} \langle x, \Phi x \rangle$ is convex, $s : x \mapsto Sx$ preserves the Lebesgue measure on the ball, and one indeed has (26).

Let $\mathcal{P}_N$ be a partition of $X$ into $N$ regions of area $1/N$ and diameter $\approx N^{-1/d}$, and let $\mathbb{M}_N \subseteq L^2(X, \mathbb{R}^d)$ be the subset of maps which are piecewise constant on the partition $\mathcal{P}_N$. The push-forward of the Lebesgue measure by some $m \in \mathbb{M}_N$ is a discrete probability measure equidistributed (if $m$ is injective) on the image of $m$

$$ m_{\#} \mathrm{Leb}_X = \frac{1}{N} \sum_{x \in \mathrm{Im}(m)} \delta_x. \quad (27) $$

Evaluating the distance (25), from $m \in \mathbb{M}_N$ to incompressible maps $\mathbb{S}$, thus amounts to a semi-discrete optimal transport problem. We regard this as a numerical commodity in view of recent progress in the scalability, speed and reliability of this method, see §1.3.

<footnote n="4">⁴The square root of the optimal transport cost.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 7: Left: partition $\mathcal{P}_N$ of the domain $X$ into $N$ regions. Right: illustration of the discretization (28). The maps $(m_t)_{t=0}^T$ belong to a linear subspace $M_N$ (black) of $L^2(X)$, whereas $\pi_S(m_t)$ denotes their projection onto the manifold $\mathbb{S}$ (blue).</figcaption></figure>

**Euler's boundary value problem via semi-discrete optimal transport.** Let $T$ and $N$ be positive integers, accounting respectively for the number of time steps and the number of regions in the partition $\mathcal{P}_N$ of $X$, and let $\lambda \ge 0$ be a penalization parameter. For all $m \in \mathbb{M}_N^{T+1}$ we define

$$ \mathcal{E}_{T,N,\lambda}(m) := T \sum_{0 \le i < T} \|m_i - m_{i+1}\|^2 + \lambda \left( \sum_{1 \le i < T} d^2(m_i, \mathbb{S}) + \|m_0 - s_*\|^2 + \|m_T - s^*\|^2 \right). \quad (28) $$

This functional depends on $N(T+1)$ variables. It implements the constraints of incompressibility, and the boundary values, by penalization as $\lambda \gg 1$. We show in [MM15], that minimizers of $\mathcal{E}(T, N_T, \lambda_T)$ converge in a weak sense to minimizers of the linear relaxation of Euler's equations §2.3, for suitable sequences $T, N_T, \lambda_T \to \infty$. The set $\mathbb{S}$ of measure preserving maps is non-convex, thus the distance $m \in L^2(X, \mathbb{R}^d) \mapsto d^2(m, \mathbb{S})$ is non convex, hence the optimized functional (28) is not convex either. Nevertheless a minimization using the LBFGS algorithm, and a multiscale (in time) initialization strategy, produced sensible results.

Figures 8 and 9 from [MM15] illustrate, for the first time, a sharp result of Brenier [Bre89b]: consider the linear relaxation of Euler equations on the time interval $[0, t_{\max}]$, and the unique associated pressure field $p$, see §2.3. If $t_{\max}^2 \nabla^2 p < \pi \text{Id}$, then the relaxation admits a unique minimizer, which is a classical solution to Euler equations. If this threshold is exceeded, then several minimizers may exist, and some or all of them may be non-classical generalized flows $\mu \in \text{Prob}(\Omega)$. In that case the support of the flow $\mu$ has dimension $> d$, so that fluid particles follow non-deterministic paths, forming complex and turbulent structures.

## References

[Abg09] R Abgrall. Construction of Simple, Stable, and Convergent High Order Schemes for Steady First Order Hamilton-Jacobi Equations. *SIAM Journal on Scientific Computing*, 31(4):2419-2446, January 2009.

[AF07] Luigi Ambrosio and Alessio Figalli. On the regularity of the pressure field of Brenier's weak solutions to incompressible Euler equations. *Calculus of Variations and Partial Differential Equations*, 31(4):497-509, August 2007.

[AHA98] F Aurenhammer, F Hoffmann, and B Aronov. Minkowski-Type Theorems and Least-Squares Clustering. *Algorithmica*, 20(1):61-76, January 1998.
<figure><img src="image_10.png" /><figcaption>Figure 8: (First row) The beltrami flow in the unit square, shown at various timesteps, is a classical solution to Euler equations. The particles color depends on their initial position. (Second to fifth row) Generalized fluid flows reconstructed by our algorithm, using the boundary conditions displayed in the first and last column. When $t_{max} < 1$ we recover the classical flow, while for $t_{max} \ge 1$ the solution is not classical anymore and includes some mixing.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Particles trajectories in the case $t_{max} = 1.5$, see Figure 8. Left: trajectories originating from inside this red circle seem strongly non-deterministic. Center: trajectories originating from this circle remain grouped. Right: particles around the domain center move little.</figcaption></figure>

[AM09] Nestor E Aguilera and Pedro Morin. On Convex Functions and the Finite Element Method. *SIAM Journal on Numerical Analysis*, 47(4):3139–3157, January 2009.

[Arn66] Vladimir Arnold. Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits. *Annales de l'institut Fourier*, 16(1):319–361, 1966.

[BBG02] J D Benamou, Y Brenier, and K Guittet. The Monge-Kantorovitch mass transfer and its computational fluid mechanics formulation. *International Journal for Numerical Methods in Fluids*, 40(1-2):21–30, 2002.

[BCC$^{+}$15] J D Benamou, G Carlier, M Cuturi, G Peyré, and L Nenna. *Iterative Bregman projections for regularized transportation problems*. Sci. Comp, 2015.

[BCM15] Jean-David Benamou, Francis Collino, and Jean-Marie Mirebeau. Monotone and Consistent discretization of the Monge-Ampère operator. *Mathematics of computation*, September 2015.

[BFO10] Jean-David Benamou, Brittany D Froese, and A M Oberman. Two numerical methods for the elliptic Monge-Ampère equation. *M2AN. Mathematical Modelling and Numerical Analysis*, 44(4):737–758, 2010.

[BFO14] Jean-David Benamou, Brittany D Froese, and A M Oberman. Numerical solution of the Optimal Transportation problem using the Monge–Ampère equation. *Journal of Computational Physics*, 260:107–126, March 2014.

[BFS09] Marc Bernot, Alessio Figalli, and Filippo Santambrogio. Generalized solutions for the Euler equations in one and two dimensions. *Journal de Mathématiques Pures et Appliquées*, 91(2):137–155, February 2009.

[Bre89a] Yann Brenier. A combinatorial algorithm for the Euler equations of incompressible flows. In *Proceedings of the Eighth International Conference on Computing Methods in Applied Sciences and Engineering (Versailles, 1987)*, pages 325–332, 1989.

[Bre89b] Yann Brenier. The least action principle and the related concept of generalized flows for incompressible perfect fluids. *Journal of the American Mathematical Society*, 2(2):225–255, 1989.
[Bre91] Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. *Communications on Pure and Applied Mathematics*, 44(4):375–417, June 1991.

[Bre93] Y Brenier. The dual least action principle for an ideal, incompressible fluid. *Archive for rational mechanics and analysis*, 1993.

[Bre08] Y Brenier. Generalized solutions and hydrostatic approximation of the Euler equations. *Physica D: Nonlinear Phenomena*, 2008.

[BV14] Martins Bruveris and François Xavier Vialard. On Completeness of Groups of Diffeomorphisms. *arXiv.org*, March 2014.

[cga] CGAL, Computational Geometry Algorithms Library. http://www.cgal.org.

[CP92] Andrew Carverhill and Franz J Pedit. Global solutions of the Navier-Stokes equation with strong viscosity. *Annals of Global Analysis and Geometry*, 10(3):255–261, 1992.

[DCMT14] Pedro Machado Manhães De Castro, Q Merigot, and Boris Thibert. Intersection of paraboloids and application to Minkowski-type problems. In *Computational geometry (SoCG'14)*, pages 308–317. ACM, New York, 2014.

[dGBOD12] Fernando de Goes, Katherine Breeden, Victor Ostromoukhov, and Mathieu Desbrun. Blue noise through optimal transport. *ACM Transactions on Graphics (TOG)*, 31(6):171, November 2012.

[dGWH$^{+15}$] F de Goes, C Wallez, J Huang, U Zhejiang, and D Pavlov. Power Particles: An incompressible fluid solver based on power diagrams. *geometry.caltech.edu*, 2015.

[DLSJ09] C De Lellis and L Székelyhidi Jr. The Euler equations as a differential inclusion. *Annals of mathematics*, 2009.

[Eul65] Leonhard Euler. Theoria Motus Corporum Solidorum Seu Rigidorum, 1765.

[FD12] Alessio Figalli and S Daneri. Variational models for the incompressible Euler equations. *HCDTE Lecture Notes, Part II*, 2012.

[FO11] Brittany D Froese and A M Oberman. Convergent Finite Difference Solvers for Viscosity Solutions of the Elliptic Monge–Ampère Equation in Dimensions Two and Higher. *SIAM Journal on Numerical Analysis*, 49(4):1692–1714, January 2011.

[Kan42] L V Kantorovich. *Kantorovich: On the transfer of masses*. Dokl. Akad. Nauk. SSSR, 1942.

[KT92] Hung-Ju Kuo and Neil S Trudinger. Discrete Methods for Fully Nonlinear Elliptic Equations. *SIAM Journal on Numerical Analysis*, 29(1):123–135, February 1992.

[Lév14] Bruno Lévy. A numerical algorithm for $L_2$ semi-discrete optimal transport in 3D. *arXiv.org*, September 2014.

[LR05] Grégoire Loeper and Francesca Rapetti. Numerical solution of the Monge-Ampère equation by a Newton’s algorithm. *Comptes Rendus Mathématique. Académie des Sciences. Paris*, 340(4):319–324, 2005.
[Mer11] Q Merigot. A Multiscale Approach to Optimal Transport. *Computer Graphics Forum*, 30(5):1583–1592, August 2011.

[Mir15] Jean-Marie Mirebeau. Adaptive, anisotropic and hierarchical cones of discrete convex functions. *Numerische Mathematik*, pages 1–47, 2015.

[MM15] Q Merigot and Jean-Marie Mirebeau. Minimal geodesics along volume preserving maps, through semi-discrete optimal transport. *arXiv.org*, May 2015.

[MO14] Q Merigot and Edouard Oudet. Handling convexity-like constraints in variational problems. *SIAM Journal on Numerical Analysis*, 52(5):2466–2487, 2014.

[MTW05] Xi-Nan Ma, Neil S Trudinger, and Xu-Jia Wang. Regularity of Potential Functions of the Optimal Transportation Problem. *Archive for rational mechanics and analysis*, 177(2):151–183, 2005.

[OP89] V I Oliker and L D Prussner. On the numerical solution of the equation $\frac{\partial^2 z}{\partial x^2} \frac{\partial^2 z}{\partial y^2} - \left(\frac{\partial^2 z}{\partial x \partial y}\right)^2 = f$ and its discretizations, I. *Numerische Mathematik*, 54(3):271–293, May 1989.

[OR15] A M Oberman and Yuanlong Ruan. An efficient linear programming method for Optimal Transportation. *arXiv.org*, September 2015.

[Sch15] Bernhard Schmitzer. A sparse algorithm for dense optimal transport. In *Scale Space and Variational Methods in Computer Vision*, pages 629–641. Springer International Publishing, Cham, May 2015.

[Shn94] A I Shnirelman. Generalized fluid flows, their approximation and applications. *Geometric and Functional Analysis*, 4(5):586–620, 1994.

[Vil08] C Villani. Optimal Transport: Old and New. 2008.