Hugo Lavenant

Sebastian Claici

Dynamical Optimal Transport on Discrete Surfaces

Keywords: CCS Concepts:, Computing methodologies → Shape analysis, • Mathematics of computing → Interpolation, Convex optimization, Partial differential equations, Optimal Transport, Wasserstein Distance, Discrete Differential Geometry

Fig. 1. Given two probability distributions over a discrete surface (left and right), our algorithm generates an interpolation that takes the geometric structure of the surface into account.

, adapted to discrete surfaces. Our structure-preserving construction yields a Riemannian metric on the (finitedimensional) space of probability distributions on a discrete surface, which translates the so-called Otto calculus to discrete language. From a practical perspective, our technique provides a smooth interpolation between distributions on discrete surfaces with less diffusion than state-of-the-art algorithms involving entropic regularization. Beyond interpolation, we show how our discrete notion of optimal transport extends to other tasks, such as distribution-valued Dirichlet problems and time integration of gradient flows.

INTRODUCTION

Probability distributions are key objects in geometry processing that can encode a variety of quantities, including uncertain feature locations on a surface, color histograms, and physical measurements like the density of a fluid. A central problem related to distributions is that of interpolation: Given two probability distributions over a fixed domain, how can one transition smoothly from the first to the second?

Optimal transport gives one potential solution. This theory lifts the geometric structure of a surface to a Riemannian structure on the space of probability distributions over the surface, the latter being endowed with the so-called Wasserstein metric; the set of distributions equipped with this metric is sometimes called Wasserstein space. To interpolate between two probability distributions, one computes a geodesic in Wasserstein space between the two. This definition is sometimes referred to as McCann's displacement interpolation [1997], applied to graphics e.g. in [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF].

Even though optimal transport theory is now well-understood [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF][START_REF] Villani | Topics in Optimal Transportation[END_REF][START_REF] Villani | Optimal Transport: Old and New[END_REF], the interpolation problem remains challenging numerically. Related problems, like the computation of Wasserstein distances or barycenters in Wasserstein space, can be tackled by fast and scalable algorithms like entropic regularization or semi-discrete methods, developed only a few years ago. Most of these methods, however, fail to reproduce the Riemannian structure of Wasserstein space and/or are prone to diffusion: The interpolation between two peaked probability distributions is more diffuse in the midpoint than optimal transport theory suggests. This drawback can inhibit application of transport in computer graphics practice, in which blurry interpolants are often undesirable.

As an alternative, we define a Riemannian structure on the space of probability distributions over a discrete surface, designed to mimic that of the Wasserstein distance between distributions over a smooth manifold. Our construction is inspired by the Benamou-Brenier formula [2000], previously discretized only on flat grids without structure preservation. This Riemannian structure automatically defines geodesics and distances between probability distributions. In particular, the geodesic problem can be recast as a convex problem and be tackled by iterative methods phrased using local operators familiar in geometry processing and finite elements (gradients, divergence and Laplacian on the surface). Our method does not require precomputation of pairwise distances between points on the surface.

Compared to other methods, our interpolation can be rephrased as a geodesic problem and numerically exhibits less diffusion when interpolating between peaked distributions. In cases where the sharpness captured by our method and predicted by optimal transport theory is undesirable visually, we provide a quadratic regularizer that controllably reduces congestion of the computed interpolant; unlike entropically-regularized transport, however, our optimization problem does not degenerate or become harder to solve when the regularization term vanishes. Although the computation of interpolants remains quite slow for meshes with more than a few thousand vertices and improving the scalability of numerical routines used to optimize our convex objective remains a challenging task for future work, we demonstrate application to tasks derived from transport, e.g. computation of harmonic mappings into Wasserstein space and integration of gradient flows.

In addition to our algorithmic contributions, we regard our work as a key theoretical step toward making optimal transport compatible with the language of discrete differential geometry (DDG). Our Riemannian metric induces a true geodesic distance-with a triangle inequality-on the space of distributions over a triangulated surface expressed using one value per vertex. Inspired by an analogous construction on graphs [START_REF] Maas | Gradient flows of the entropy for finite Markov chains[END_REF]], we leverage a non-obvious observation that a strong contender for structure-preserving discrete transport on meshes actually involves a real-valued external time variable, rather than discretizing transport as a linear program as in most previous work. The resulting geodesic problem naturally preserves convexity and other key properties from the theoretical case while suggesting an effective computational technique.

RELATED WORK 2.1 Linear programming and regularization

Landmark work by [START_REF] Kantorovich | On the translocation of masses[END_REF] showed that optimal transport can be phrased as a linear programming problem. If both probability distributions have finite support, we end up with a finitedimensional linear program solvable using standard convex programming techniques. A variety of solvers has been designed to tackle this linear program, which exploit the particular structure of the objective functional [START_REF] Edmonds | Theoretical improvements in algorithmic efficiency for network flow problems[END_REF][START_REF] Klein | A primal method for minimal cost flows with applications to the assignment and transportation problems[END_REF][START_REF] James B Orlin | A polynomial time primal network simplex algorithm for minimum cost flows[END_REF]]. These methods, however, usually require as input the pairwise distance matrix, a dense matrix that scales quadratically in the size of the support and is difficult to evaluate if the points are on a curved space.

A landmark paper by [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] reinvigorated interest in numerical transport by proposing adding an entropic regularizer to the problem, leading to the efficient Sinkhorn (or matrix rebalancing) algorithm. This algorithm, which involves iteratively rescaling the rows and columns of a kernel in the cost matrix, is highly parallelizable and well-suited to GPU architectures. When the cost matrix involves squared geodesic distances along a discrete surface, [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF] showed that Sinkhorn iterations can be written in terms of heat diffusion operators, eliminating the need to store the cost matrix explicitly. While they are efficient, these entropically-regularized techniques suffer from diffusion, making them less relevant to problems in which measures are sharp or peaked. They also do not define true distances on the space of distributions over mesh vertices.

When the transport cost is equal to geodesic distance, i.e. the 1-Wasserstein distance, optimal transport is equivalent to the Beckmann problem [Santambrogio 2015, Chapter 4], for which specific and efficient algorithms can be designed [START_REF] Li | A parallel method for earth mover's distance[END_REF][START_REF] Solomon | Earth mover's distances on discrete surfaces[END_REF]. These methods cannot be applied to the quadratic Wasserstein distance, which is needed to make transport-based interpolation nontrivial, namely to recover McCann's displacement interpolation [1997]. In particular, the optimal transport problem defining the 1-Wasserstein distance does not come with a time dependency allowing to define a smooth interpolation and suffers from nonuniqueness coming from the lack of strict convexity.

Semi-discrete optimal transport

When one of the distributions has a density w.r.t. Lebesgue while the other one is discrete, the transport problem can be reduced to a finite-dimensional convex problem whose number of unknowns scales with the cardinality of the support of the discrete distribution. Leveraging tools from computational geometry, this semi-discrete problem can be solved efficiently up to fairly large scale when the cost is Euclidean [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Kitagawa | Convergence of a Newton algorithm for semi-discrete optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for L 2 semi-discrete optimal transport in 3D[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF].

Semi-discrete transport has been used to tackle problems for which the precise structure of the optimal transportation map is relevant, as in fluid dynamics [de Goes et al. 2015b;[START_REF] Thomas | A Lagrangian scheme à la Brenier for the incompressible Euler equations[END_REF][START_REF] Mérigot | Minimal geodesics along volumepreserving maps, through semidiscrete optimal transport[END_REF]. It also has been used for approximating barycenters in the stochastic case [START_REF] Claici | Stochastic Wasserstein barycenters[END_REF]] and as a measure of proximity for shape reconstruction [START_REF] Fernando De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF][START_REF] Digne | Feature-preserving surface reconstruction and simplification from defect-laden point sets[END_REF]]. Extensions of semi-discrete transport to curved spaces can be found in [de Goes et al. 2014;[START_REF] Mérigot | An algorithm for optimal transport between a simplex soup and a point cloud[END_REF]. Although they can be fast and give explicit transport maps, these methods are not suited for the application we have in mind: They rely on the computation of transport maps between two probability distributions that are not of the same nature (one is discrete, the other has a density) and hence cannot be used to implement a distance or interpolation cleanly.

Fluid dynamic formulations

By switching from Lagrangian to Eulerian descriptions of transport, [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] proved that optimal transport could be rephrased using fluid dynamics: Instead of computing a coupling, they show that transport with quadratic costs is equivalent to finding a time-varying sequence of distributions smoothly interpolating between the two measures. The problem that they obtain is convex and solved via the Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Proof of the convergence of ADMM in the infinite-dimensional setting (i.e. when neither time nor the geometric domain is discretized) is provided in [START_REF] Guittet | On the time-continuous mass transport problem and its approximation by augmented Lagrangian techniques[END_REF][START_REF] Hug | On the convergence of augmented Lagrangian method for optimal transport between nonnegative densities[END_REF]. [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF] reread the ADMM iterations as a proximal splitting scheme and show how one can build different algorithms to solve the convex problem. This fluid dynamic formulation also appears in mean field games [Benamou and Carlier 2015].

In all of the above work, however, the authors work in a flat space and use finite difference discretizations of the densities and velocity fields. Hence their work does not contain a clear indication about how to handle the problem on a discrete curved space, and theoretical properties of their models after discretization remain unverified.

The algorithm for approximating 1-Wasserstein distances presented by [START_REF] Solomon | Earth mover's distances on discrete surfaces[END_REF] achieves some of the objectives mentioned above. Their vector field formulation is in some sense dynamical, and their distance satisfies properties like the triangle inequality after discretization. As mentioned above, however, their optimization problem lacks strict convexity and is not suitable for interpolation.

Dynamical transport on graphs and meshes

Maas [2011] defines a Wasserstein distance between probability distributions over the vertices of a graph. The (finite-dimensional) space of distributions in this case inherits a Riemannian metric with some structure preserved from the infinite-dimensional definition; for instance, the gradient flow of entropy corresponds to a notion of heat flow along the graph. A similar structure is proposed by [START_REF] Chow | Entropy dissipation semi-discretization schemes for Fokker-Planck equations[END_REF], but they recover a different heat flow. [START_REF] Erbar | Computation of Optimal Transport on Discrete Metric Measure Spaces[END_REF] propose a numerical algorithm for approximating the discrete Wasserstein distance introduced by Maas, but the distributions they produce have a tendency to diffuse along the graph. This flaw is not related to their numerical method but rather comes from the very definition of their optimal transport distance. It is also not obvious what is the best way to adapt their construction to discrete surfaces rather than graphs.

Interpolation and geodesics

Optimal transport is not the only way to interpolate between probability distributions; for instance, [START_REF] Azencot | Advection-based function matching on surfaces[END_REF] use a timeindependent velocity field to advect functions and match them. Their method, however, cannot be understood as a geodesic curve in the space of distributions. In another direction, [START_REF] Heeren | Time-Discrete Geodesics in the Space of Shells[END_REF] have provided an efficient way to discretize in time geodesics in a highdimensional space of thin shells. Their formulation is not well-suited for optimal transport where direct discretization of the Benamou-Brenier formula is possible. Finally, methods like [START_REF] Panozzo | Weighted averages on surfaces[END_REF] provide a means of averaging points on discrete surfaces, x y

t = 0 t = 1 /4 t = 1 /2 t = 3 /4 t = 1
x ′ y ′ dπ (x, y) (1). The initial distribution μ0 is on the left, and the final distribution μ1 is on the right. The quantity dπ (x, y) represents the amount of mass that is transported from x to y. The coupling π is chosen in such a way that the total cost is minimal. Bottom row: dynamical formulation between the same distributions (computed with the algorithm in Section 4). To go from the top to the bottom row, once one has the optimal π , a proportion dπ (x, y) of particles follows the geodesic (in this case a straight line) between x and y with constant speed. The macroscopic result of all these motions is a time-varying sequence of distributions, displayed in blue.

t = 0 t = 1 /3 t = 2 /3 t = 1
although it is not clear how to extend them to the more general distribution case.

OPTIMAL TRANSPORT ON A DISCRETE SURFACE 3.1 Optimal transport on manifolds

We begin by introducing briefly optimal transport theory on a smooth space. Let M be a connected and compact Riemannian manifold with metric ⟨ , ⟩ and induced norm ∥ ∥; define d : M ×M → R + to be geodesic distance.

Denote by P(M) the space of probability measures on M. This space is endowed with the quadratic Wasserstein distance from optimal transport: If μ0 , μ1 ∈ P(M), then the distance W 2 ( μ0 , μ1 ) between them is defined as

W 2 2 ( μ0 , μ1 ) := min π ∬ M×M 1 2 d(x, y) 2 dπ (x, y), (1)
where the minimum is taken over all probability measures π on the product space M × M whose first (resp. second) marginal is μ0 (resp. μ1 ).

The problem (1) can be interpreted as follows: dπ (x, y) denotes the quantity of particles located at x that are sent to y, and the cost for such a displacement is d(x, y) 2 . The constraint on the marginals enforces that π describes a way of moving the distribution of mass μ0 onto μ1 . Thus, the variational problem (1) reads: Find the cheapest way π to send μ0 onto μ1 , and the result (i.e. the minimal cost) is defined as the squared Wasserstein distance between μ0 and μ1 . In some generic cases [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Gangbo | The geometry of optimal transportation[END_REF], the optimal π is located on the graph of a map T : M → M, which means that a particle x ∈ M is sent onto a unique location y = T (x) ∈ M.

The space (P(M),W 2 ) is a complete metric space [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF][START_REF] Villani | Topics in Optimal Transportation[END_REF], and-at least formally-it has the structure of an (infinite-dimensional) Riemannian manifold. Revealing this manifold structure requires some manipulation and rephrasing of the original problem (1), detailed below.

As first noticed by [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF], the Wasserstein distance between μ0 and μ1 can be obtained by solving an alternative, physically-motivated optimization problem:

W 2 2 ( μ0 , μ1 ) =          min µ,v ∫ 1 0 ∫ M 1 2 ∥v t ∥ 2 dµ t dt s.t. µ 0 = μ0 , µ 1 = μ1 , ∂ t µ + ∇ • (µv) = 0.
(2)

As we will have to deal with functions and vectors depending both on time and space, here and moving forward we adopt the following convention: Upper indices denote time, and lower indices denote space. Moreover, t ∈ [0, 1] will denote an instant in time, and f will later denote a generic triangle (f for face) in a triangulation. In (2), the minimum is taken over all curves µ : [0, 1] → P(M) and all time-dependent velocity fields v : [0, 1] × M → T M such that the continuity equation ∂ t µ + ∇ • (µv) = 0 is satisfied in the sense of distributions. The optimal curve µ is known as McCann's displacement interpolation [1997].

The physical interpretation of this problem is as follows. Imagine probability distributions as distributions of mass, e.g. the density of a fluid. The curve µ represents an assembly of particles in motion, distributed as μ0 at t = 0 and μ1 at t = 1. At time t, a particle located at x ∈ M moves with velocity v t

x . The continuity equation ∂ t µ + ∇ • (µv) = 0 simply expresses the conservation of mass. For a given time t, the cost ∫ M 1 2 ∥v t ∥ 2 dµ t is the total kinetic energy of all the particles. Hence, the cost minimized in (2), i.e. the integral w.r.t. time of the kinetic energy, is the action of the curve. As there is no congestion cost-that is, the particles do not interact with each other-(2) is the least-action principle for a pressureless gas.

Formulation (1) is static, since it directly determines the target for each particle at t = 1 given the arrangement at t = 0. On the other hand, (2) is dynamical, recovering a curve in P(M) interpolating smoothly between μ0 and μ1 . To convert from the static to the dynamical formulation, one takes an optimal transport plan π from (1) and an assembly of particles distributed according to μ0 . If a particle located at x ∈ M at time t = 0 and is supposed, according to π , to be sent to y ∈ M, then this particle follows a constant-speed geodesic along M from x to y. The optimal curve µ in (2) is exactly the resulting macroscopic motion of all the particles, illustrated in Figure 3.

Calling m = µv the momentum and using the change of variables

(µ, v) ↔ (µ, m), problem (2) becomes convex, because the mapping (µ, v) → 1 /2∥v∥ 2 µ is not jointly convex while (µ, m) → 1 /2∥m∥ 2 /µ is. Its dual reads W 2 2 ( μ0 , μ1 ) = max φ ∫ M φ 1 d μ1 - ∫ M φ 0 d μ0 s.t. ∂ t φ + 1 2 ∥∇φ ∥ 2 ⩽ 0 on [0, 1] × M, (3)
where the maximization is performed over real-valued functions [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF][START_REF] Villani | Optimal Transport: Old and New[END_REF]]. The relation v = ∇φ holds whenever v (resp. φ) is a minimizer (resp. maximizer) of the primal (resp. dual) problem. In particular, in (2), we can restrict ourselves to the set of v such that v t = ∇φ t for every t ∈ [0, 1]. Equation ( 2) defines a formal Riemannian structure on P(M) [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF]]. Given µ ∈ P(M) with a density bounded from below by a strictly positive constant, the tangent space T µ P(M) is identified as the set of functions δ µ : M → R with 0-mean: δ µ is the partial derivative w.r.t. time of a curve whose value at time 0 is µ. If δ µ ∈ T µ P(M), we can compute φ : M → R the solution (unique up to translation by constants) of the elliptic equation

φ : [0, 1] × M → R [
∇ • (µ∇φ) = -δ µ. (4)
Then, the norm of δ µ is defined as

∥δ µ ∥ 2 T µ P(M) := 1 2 ∫ M ∥∇φ ∥ 2 dµ. ( 5
)
Endowed with this scalar product obtained from the polarization identity ⟨x, y⟩ = 1 4 (∥x + y ∥ 2 -∥x -y ∥ 2 ), one can check, and the derivation appears in the supplemental material, that the Wasserstein distanceW 2 can be interpreted as the geodesic distance induced by ( 4) and ( 5). This is precisely the content of the Benamou-Brenier formula (2).

One needs to assume µ ⩾ c > 0 on M for the elliptic equation ( 4) to be well-posed. Nevertheless, one can still give a meaning to this Riemannian structure using tools from analysis in metric spaces [START_REF] Ambrosio | Gradient Flows: In Metric Spaces and in the Space of Probability Measures[END_REF]].

Discrete surfaces

The previous subsection contains only well-understood results. Let us now start our contribution: to mimic these definitions and properties when the manifold is replaced by a triangulated surface.

Instead of a smooth manifold M, we consider the case where we only have access to a triangulated surface S = (V , E,T ), which consists of a set V ⊂ R 3 of vertices, a set E ⊆ V ×V of edges linking vertices, and a set T ⊆ V × V × V of triangles containing exactly 3 vertices linked by 3 edges. For a given face f ∈ T , we denote by

V f ⊂ V the set of vertices v such that v ∈ f ; for a given vertex v ∈ V , we denote by T v ⊂ T the set of faces f such that v ∈ f . The area of a triangle f ∈ T is denoted by | f |. Each vertex v is associated to a barycentric dual cell (see Figure 4) whose area, equal to 1 3 f ∈T v | f |, is denoted by |v |.
Following standard constructions from first-order finite elements (FEM), a scalar function on M will be seen as having one value per vertex, i.e. belonging to R |V | . A distribution µ ∈ M will be also discretized by one value per vertex representing the density w.r.t. the volume measure. In other words, the volume of the dual cell centered at v ∈ V , measured with µ, is |v |µ v . We denote by P(S) the set of probability distributions on the discrete surface:

P(S) := µ ∈ R |V | s.t. µ v ⩾ 0 for all v ∈ V and v ∈V |v |µ v = 1 . (6)
For instance, the volume measure is represented by the vector in P(S) parallel to (1, 1, . . . , 1) ⊤ .

The set V of vertices can be interpreted as a discrete metric space, either by using directly the Eulidean distance on R 3 or by some version of the discrete geodesic distance along S. Hence, a natural attempt to discretize the 2-Wasserstein distance would be to use (1) and replace d by the distance between vertices. As pointed out in [START_REF] Gigli | Gromov-Hausdorff convergence of discrete transportation metrics[END_REF][START_REF] Maas | Gradient flows of the entropy for finite Markov chains[END_REF]], however, this discretization leads to a space without a smooth structure. For instance, there do not exist non-constant smooth (e.g., Lipschitz) curves valued in such a space; whereas in a space with a smooth structure (e.g. a Riemannian manifold), one expects the existence of non-constant Lipschitz curves, namely the (constant-speed) geodesics.

Let us briefly recall the argument. We take the simplest example of a space consisting of two points. If X = {x 0 , x 1 } contains two points separated by a given distance ℓ, a probability distribution µ on X is characterized by a single number

µ x 0 ∈ [0, 1], as µ x 1 = 1 -µ x 0 . If µ t is a curve valued in P(X ), one can compute W 2 (µ t , µ s ) = ℓ |µ t x 0 -µ s x 0 |. In particular, if µ is Lipschitz with Lipschitz constant L, our expression for W 2 implies |µ t x 0 -µ s x 0 | ⩽ L 2 ℓ 2 |t -s | 2 .
There is an exponent 2 on the r.h.s., but only 1 on the l.h.s.: it is precisely this discrepancy which is an issue. Indeed, dividing by |t -s | on both side and letting s → t, one sees that t → µ t

x 0 is differentiable everywhere with derivative 0, i.e. is constant.

For this reason, we prefer to discretize the Benamou-Brenier formulation (2), as it will automatically give a Riemannian structure on the space P(S). In this sense, the basic inspiration for our technique is the same as that of [START_REF] Maas | Gradient flows of the entropy for finite Markov chains[END_REF], although on a triangulated surface we enjoy the added structure afforded by an embedded manifold approximation of the domain rather than an abstract graph.

As (2) involves velocity fields, a choice has to be made about their representation [de Goes et al. 2015a]. To take full advantage of the triangulation, we want to use triangles and not only edges to define our objective functional. The latter choice leads to formulas similar to [START_REF] Chow | Entropy dissipation semi-discretization schemes for Fokker-Planck equations[END_REF][START_REF] Maas | Gradient flows of the entropy for finite Markov chains[END_REF], which, as we say above, exhibit strongly diffuse geodesics. We prefer to represent vector fields on triangles. More precisely, a (piecewise-constant) velocity field v is represented as an element of (R 3 ) |T | , i.e. as one vector per triangle, with the constraint that v f , which is a vector of R 3 , is parallel to the plane spanned by f , which means that our velocity fields lie in a subspace of dimension 2|T |.

If φ ∈ R |V | represents a real-valued function, we compute its gradient along the mesh using the first-order (piecewise-linear) finite element method [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]: For each triangle f , we compute φ, the unique affine function defined on f coinciding with φ on the vertices of f . Then, the gradient of φ in f is simply defined as the gradient of φ at any point of f ; as the gradient is constant on each triangle, we need to store only one vector per triangle. Since this operator is linear, let us denote by G ∈ R 3|T |×|V | its matrix representation. In particular, the Dirichlet energy of φ ∈ R |V | is defined as

Dir(φ) := 1 2 f ∈T | f |∥(Gφ) f ∥ 2 . ( 7
)
The sum is weighted by the areas of the triangles to discretize a surface integral. The first variation of this Dirichlet energy can be expressed in matrix form as

(G ⊤ M T G)φ, where M T ∈ R 3|T |×3|T |
is a diagonal weight matrix whose elements are the areas of the triangles. The matrix G ⊤ M T G is the so-called cotangent Laplace matrix of a triangulated surface [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF].

Dual problem on meshes

Let us introduce our discrete Benamou-Brenier formula by starting from its dual formulation (3). Since the objective functional is linear, its discrete counterpart is straightforward as both µ and φ are defined on vertices. On the other hand, in the constraint ∂ t φ +1 2 ∥∇φ ∥ 2 ⩽ 0, we would like to replace ∇φ by Gφ but then the two terms of the sum do not live on the same space.

The constraint ∂ t φ + 1 2 ∥∇φ ∥ 2 ⩽ 0 is a priori not coercive. Suppose φ satisfies the constraint, and take another function ψ with the property that φ + sψ satisfies the constraint for arbitrarily large s ⩾ 0. Expanding the inequality ∂ t (φ + sψ ) + 1 2 ∥∇φ + s∇ψ ∥ 2 ⩽ 0 and taking the limit s → +∞ shows that ψ satisfies this property if and only if ∥∇ψ ∥ = 0 and ∂ t ψ ⩽ 0; these two conditions together imply that the objective functional in (3) is smaller when evaluated at φ + sψ rather than at φ. This is a property that we would like to keep at the discrete level. To do so, we enforce a discrete analogue of the constraint at each vertex of the mesh. To go from ∥Gφ ∥ 2 , which is defined on triangles, to something defined on vertices, we first take the squared norm and subsequently average in space: 1 Definition 3.1. Let μ0 , μ1 ∈ P(S). The discrete (quadratic) Wasserstein distance W d ( μ0 , μ1 ) is defined as the solution of the following convex problem:

W 2 d ( μ0 , μ1 ) =              sup φ v ∈V |v |φ 1 v μ1 v -v ∈V |v |φ 0 v μ0 v s.t. ∂ t φ t v + 1 2 f ∈T v | f |∥(Gφ) t f ∥ 2 3|v | ⩽ 0 for all (t, v) ∈ [0, 1] × V , (8)
where the unknown is a function φ

: [0, 1] × V → R.
The denominator 3|v | is nothing else, by definition, than

f ∈T v | f |.
In particular, the value

f ∈T v | f |∥(Gφ) t f ∥ 2 (3|v |) -1 is the average,
weighted by the areas of the triangles, of ∥(Gφ) t f ∥ 2 for f ∈ T v . One can check that the same reasoning as above can be performed. Indeed, if φ : [0, 1] ×V → R satisfies the constraint in (8) and φ + sψ also satisfies it for arbitrarily large s ⩾ 0, it implies, taking s → +∞, that 1 2

f s.t. v ∈f | f |∥(Gψ ) t f ∥ 2 3|v | ⩽ 0. ( 9
)
This inequality must hold for all (t, v) ∈ [0, 1]×V . Thus, we conclude (and it is for this implication that it is important to average after taking squares) that Gψ is identically 0. In other words, for all t ∈ [0, 1], the function ψ t is constant over the discrete surface. Plugging this information back into the constraint in ( 8) and taking again s → +∞, we see that ∂ t ψ ⩽ 0. Hence, the value ψ 0 (which is constant over the surface) is larger than ψ 1 . With this information (Gψ = 0 and ∂ t ψ ⩽ 0), the value of the objective functional must be smaller for φ + sψ than for φ as soon as s ⩾ 0.

Riemannian structure of the space of probabilities on a discrete surface

To recover an equation which looks like the primal formulation of the Benamou-Brenier formula (2), it is enough to write the dual of the discrete formulation (8). The latter formulation, as explained above, was important to justify the choice of the way we average quantities that do not live on the same grid. We introduce additional notation to deal with the averaging of the density µ. If µ ∈ P(S), we denote by μ ∈ R |T | the vector given by, for any f ∈ T , μf =

1 3 v ∈V f µ v . (10)
This is a natural way to average µ from vertices to triangles, appearing in the dual formulation given below:

Proposition 3.2. The following identity holds:

W 2 d ( μ0 , μ1 ) =                min µ,v ∫ 1 0 f ∈T 1 2 ∥v t f ∥ 2 | f | μt f dt s.t. µ 0 = μ0 , µ 1 = μ1 ∂ t (M V µ s v ) + (-G ⊤ M T [ μt v t ]) v = 0 for all (t, v) ∈ [0, 1] × V . (11) Recall that M T ∈ R 3|T |×3|T | and M V ∈ R |V |× |V |
are diagonal matrices corresponding to multiplication by the area of the triangles and of the dual cells respectively. Then, -G ⊤ M T represents a discrete version of the (integrated) divergence operator, suggesting that the constraint can be interpreted as a discrete continuity equation. The derivation of this result, detailed in the supplemental material, relies on an inf-sup exchange, similar to the case of a smooth surface M.

Proposition 3.2, very similar to (2), shows that W d is the geodesic distance for a Riemannian structure on the space P(S), at least for non-vanishing densities. Let us detail the metric tensor for a density µ ∈ P(S) with min v µ v > 0. As the set P(S) is a codimension-1 subset of the linear space R |V | , the tangent space at µ is naturally {x ∈ R |V | s.t. v ∈V |v |x v = 0}. In analogy to (4), take δ µ ∈ T µ P(S). We call φ ∈ R |V | a solution of where M μ ∈ R 3|T |×3|T | is a diagonal matrix corresponding to multiplication on each triangle by μ. As μ > 0 everywhere on V , this equation is well-posed: The kernel of G ⊤ M T M μG is of dimension one (it consists only of the constant functions), and M V δ µ lies in the image of this operator. When the distribution µ is uniform, ( 12) boils down to a Poisson equation, as the operator -(G ⊤ M T M μG) is proportional to the cotangent Laplacian.

M V δ µ = -(G ⊤ M T M μG)φ, (12)
One can then define the norm of δ µ on the tangent space T µ P(S)

as ∥δ µ ∥ 2 T µ P(S ) := 1 2 f ∈T ∥(Gφ) f ∥ 2 | f | μf . ( 13
)
The function φ is unique up to an additive constant, which lies in the kernel of the matrix G, so this norm is well-defined.

To put everything in one formula, the scalar product ⟨δ µ, δν ⟩ T µ P(S ) between two elements of the tangent space at µ can be expressed as (δν ) ⊤ P µ (δ µ), where the matrix P µ is expressed as

P µ = 1 2 M ⊤ V G -⊤ (M μ M T ) -1 G -1 M V . ( 14
)
One can check, and the derivation is provided in the supplemental material, that W d is exactly the geodesic distance induced by this metric tensor.

Proposition 3.3. The function W d : P(S) × P(S) is a distance.

Proof. It is a general fact that the geodesic distance on a manifold (defined by minimization over all possible trajectories) is a distance, see for instance [Jost 2008, Section 1.4].

□

A natural question is whether the space (P(S),W d ) converges to (P(M),W 2 ) as S becomes a finer and finer discretization of a manifold M. For a discrete Wasserstein distance like the one of [START_REF] Maas | Gradient flows of the entropy for finite Markov chains[END_REF], based on the graph structure of S-which corresponds to the case where velocity fields are discretized by their values on edges and a particular choice of scalar product-the answer is known to be positive in the case where M is the flat torus [START_REF] Gigli | Gromov-Hausdorff convergence of discrete transportation metrics[END_REF][START_REF] Garcia | Gromov-Hausdorff limit of Wasserstein spaces on point clouds[END_REF] in the sense of Gromov-Hausdorff convergence of metric spaces, while a very recent work by Gladbach, Kopfer and Mass [2018] has refined the analysis and exhibits necessary conditions for such a convergence to hold. The high technicality of the proofs of these results, however, indicates that the question for our particular definition is likely to be challenging and out of the scope of the present article.

TIME DISCRETIZATION OF THE GEODESIC PROBLEM 4.1 Discrete geodesic

So far, we have defined a structure-preserving notion of optimal transport on a triangle mesh. While our model has many properties in common with the continuum version of transport, the resulting optimization problem is infinite-dimensional since the unknown µ t is indexed by a time variable t ∈ [0, 1]. Our next step is to derive a time discretization that approximates this interpolant in practice. Put simply, we want to solve the geodesic problem, i.e., given μ0 , μ1 ∈ P(S), we want to approximate the solution µ of (11). To this end, we discretize in time the dual problem (8).

Our infinite-dimensional problem can be classified as a secondorder cone program (SOCP) [Boyd and Vandenberghe 2004, Section 4.4.2]; we choose a discretization that preserves this structure. The main issue is that with a standard finite difference scheme, the derivative ∂ t φ ends up on a grid staggered w.r.t. the one on which φ is defined. Hence, we average to define the constraint on a compatible grid. We apply the same idea as before: With the term involving ∥Gφ ∥ 2 , we average after taking the square to avoid the introduction of any spurious null space.

Let N be the number of discretization points in time. We consider two grids: the staggered grid G st time := {k/N : k = 0, 1, . . . , N } and the centered grid G c time := {(k + 1/2)/N : k = 0, 1, . . . N -1}, see Figure 4 The staggered grid has N + 1 elements whereas the centered one has only N . We call τ := 1/N the time step. The linear operator D : R G st time → R G c time defined by

(Dφ) t := φ t +τ /2 -φ t -τ /2 τ , (15)
f ∈Tv |f | ∥(Gφ) t f ∥ 2 3|v |
, which is defined on G st time , will be also averaged in time. In other words, the fully discrete problem reads:

Find φ ∈ R G st time × |V | maximizing                  v ∈V |v |φ 1 v μ1 v -v ∈V |v |φ 0 v μ0 v s.t. (Dφ) t v + 1 2 i∈{-1,1} 1 2 f ∈T v | f |∥(Gφ) t +iτ /2 f ∥ 2 3|v | ⩽ 0 for all (t, v) ∈ G c time × V , (16)
The constraint still stays quadratic, and hence the fully-discrete problem is still a SOCP.

Algorithm

To tackle (16) algorithmically, we follow [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] by building an augmented Lagrangian and using the Alternating Direction Method of Multipliers (ADMM). The main issue is that the constraint is nonlocal-since it involves discrete derivatives-and nonlinear. We construct a splitting of the problem that decouples these two effects.

To this end, we introduce two additional variables A and B. We enforce the constraint A = Dφ, and hence A is defined on the grid G c time × V . On the other hand, the variable B stores the values of Gφ but with some redundancy. Each (Gφ) t f appears in more than one inequality constraint in ( 16), and B is chosen so that each component of B appears in only one inequality constraint. In detail, B is defined on the grid G c time × {±1} × T × V with the constraint that (f , v) ∈ T ×V is such that v ∈ f . We will impose the constraint that B t,i f ,v = (Gφ)

t +iτ /2 f for all (t, i, f , v) ∈ G c time × {±1} × T × V .
We introduce the notation q = (A, B) and write q = Λφ if A, B satisfy the relations written above. Define

F (φ) = v ∈V |v |φ 1 v μ1 v - v ∈V |v |φ 0 v μ0 v , (17)
and C to be the function such that C(A, B) = C(q) = 0 if

∀(t, v) ∈ G c time × V , A t v + 1 2 i ∈ {-1,1} 1 2 f ∈T v | f |∥B t,i v, f ∥ 2 3|v | ⩽ 0 (18)
and -∞ otherwise. The discrete problem ( 16) can be written

max q=Λφ F (φ) + C(q). ( 19
)
The idea is to introduce a Lagrange multiplier σ = (µ, m) associated to the constraint q = Λφ and to build the augmented Lagrangian

L(φ, q, σ ) = F (φ) + C(q) + ⟨σ , q -Λφ⟩ - r 2 ∥q -Λφ ∥ 2 . ( 20
)
In this equation, ⟨σ , q -Λφ⟩ = ⟨µ, A -Dφ⟩ V + ⟨m, B -Gφ⟩ T , where the scalar product ⟨ , ⟩ V (resp. ⟨ , ⟩ T ) is weighted by the areas of the vertices (resp. the triangles) and the time step τ . The saddle points of the Lagrangian (20)-which do not depend on the parameter r -are precisely the solutions to the problem (16), and µ, the first component of σ associated to the constraint A = Dφ, is an approximation of the time-continuous geodesic (11). On the other hand, the second component m is an approximation of the momentum µv.

To compute a saddle point, we use ADMM, which consists in iterations of the following form [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]]:

(1) Given q and σ , find φ that maximizes L.

(2) Given φ and σ , find q that maximizes L.

(3) Do a gradient descent step (with step r ) to update σ .

The parameter r > 0 is arbitrary and tuned to speed up the convergence; see [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] for discussion. In our case, details of the iterations are briefly presented below and summarized in Algorithm 1.

Maximization w.r.t. φ. The Lagrangian L is simply a quadratic function of φ, so its maximization amounts to inverting the matrix Λ ⊤ Λ which, in our case, behaves like a space-time Laplacian.

More precisely, writing φ ∈ R

r D ⊤ M V Dφ + 3(E ⊤ E)φ(G ⊤ M T G) = N ( μ1 I t =1 -μ0 I t =0 ) -D ⊤ M V (µ -rA) -(m -r B)M T G⊤ . ( 21
)
Again recall that the unknown here is φ; the remaining symbols are fixed matrices. In this equation, E ∈ R G c time ×G st time stands for the averaging in time defined by (Eφ) t = φ t -τ /2 +φ t +τ /2 2 . The matrices I t =0 and I t =1 ∈ R G st time ×V stand for the indicator of t = 0 (resp. t = 1), namely they contain zeros except on the first (resp. last) row which is full of ones. The factor 3 comes from the fact that each value of (Gφ) f is duplicated 3 times in B, one for each vertex which belongs to f . The operator G is almost the same as G but takes in account the fact that the values of Gφ are duplicated in B (hence in m): G corresponds to the adjoint of the second component of the operator Λ.

(D ⊤ M V D) is the discrete Laplacian in time, and G ⊤ M T G is the discrete Laplacian on S. In fact, ( 21) is a Poisson equation with a space-time Laplacian. Equation ( 21) admits more than one solution but they only differ by a constant whose value does not modify the value of L.

The linear operator to invert is the same at each iteration, and hence standard precomputation techniques can be used to speed up the application of its inverse.

Maximization w.r.t. A, B. The Lagragian L is also quadratic w.r.t. q, but there is a quadratic constraint on these two variables due to the presence of C(q). Because of the redundancy in B, each component of A or B is subject to only one constraint. More precisely, we can check that one needs, for each

(t, v) ∈ G c time × V , to minimize |v | A t v -(Dφ) t v - 1 r µ t v 2 + | f | 2 i ∈ {±1} v ∈V f B t,i f ,v -(Gφ) t +iτ /2 f - 1 r m t,i f ,v 2 (22)
under the constraint (18). This minimization amounts to a Euclidean projection on the set of A, B satisfying (18), which can be carried out by solving a cubic equation in one variable, independently on each point of G c time × V . These equations are solved using Newton's method.

Dual update. This gradient descent corresponds to the following operations:

µ t v ← µ t v -r A t v -(Dφ) t v m t,i f ,v ← m t,i f ,v -r B t,i f ,v -(Gφ) t +iτ /2 f , (23)
for any (t, v) ∈ G c time ×V and any (t, i, f , v) ∈ G c time × {±1} ×T ×V .

EXPERIMENTS

Recall that our main practical contribution is to be able interpolate between probability distributions using an optimal transport model that preserves structure from the non-discretized case. We will illustrate the robustness of our method: It can handle peaked distributions, and it lifts the intrinsic geometry of the discrete surface while being insensitive to the choice of mesh topology. The typical computation is the following: We enter the data μ0 , μ1 and compute a solution of the discrete problem ( 16). Then, we plot the evolution over time of µ, which approximates the geodesic in the Riemannian metric described in Subsection 3.4. As a byproduct of the optimization process, we also obtain the optimal momentum m = µv, which can be also plotted, see Figure 2. The code used to conduct all our experiments is available at https://github.com/ HugoLav/DynamicalOTSurfaces.

As the color map is sometimes normalized independently for different time instants on the same interpolation curve, let us underscore this fact: For every example, we have checked numerically that the densities are always nonnegative and that mass is always preserved over time.

Convergence of the ADMM iterations

For fixed boundary data μ0 and μ1 , we plot in Figure 5 the primal and dual error defined by [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], as a function of the number of iterations of the ADMM scheme. We usually need on the order of a few thousand iterations to satisfy our convergence criteria, this number being dependent of the boundary data μ0 , μ1 (the more diffuse, the fewer iterations are needed).

Because our objective functional is scaled according to the geometry of the mesh (i.e. scalar products are weighted by the areas of the triangles and the number of time steps), the number of iterations needed does not depend on the size of the resolution of the mesh nor the number of discretization points in time, but the computation time needed per iteration does. Typical values of the timings are given in Table 1, they are of the order of 1 second per ADMM iterations for meshes with a few thousand vertices.

CVX implementation

Since the optimization problem in Equation ( 25) is a convex cone problem, we have also used a straightforward implementation in CVX [Grant andBoyd 2008, 2014], with Mosek as a solver [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF]. This approach is provided as a simpler alternative to the ADMM implementation, and has comparable performance on small meshes with standard precision settings (fewer than 1000 vertices). In general, it is difficult to compare the error thresholds across the two implementations due to algebraic rearrangements performed by CVX. See Table 1. 1. Timing data for various meshes and boundary data from the figures (numbers listed in the table). N denotes the number of time discretization points and α is the value of the congestion regularization parameter (see Section 5.4). For the ADMM method, the number of iterations and timing are given. Iterations were stopped once an error of 10 -4 was reached for the L 2 norm of both the primal and dual residual. One can see that the time per iteration depends on the size of the mesh and the temporal grid, but the number of iteration is quite insensitive to these parameters and rather depends on the boundary conditions and the regularization parameter. For the CVX implementation of the optimization problem, the solver time and the total time (includes CVX pre-processing) are given. Standard precision settings were used, but are hard to interpret absolutely due to unknown algebraic rearrangement of the problem. * denotes that CVX reported a failure in this case. Results obtained on an 8-core 3.60GHz Intel i7 processor with 32GB RAM. . Amplitude of the primal and dual residual [Boyd et al. 2011, Section 3.3] in L 2 norm. The distributions μ0 and μ1 are delta functions located on respectively the right and left hand of the armadillo. We also show the midpoint µ 1/2 for different numbers of iterations (10,50 and 5000). After a few hundred iterations, there is no visible difference in µ 1/2 . There is a jump in the value of the dual residual at around 4600 iterations. It is due to a change in the value of the parameter r , which is updated according to the heuristic rule presented in Section 3.4.1 of [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF].

Convergence with discretization in space and time

As indicated in Section 3, it is not known theoretically whether our discrete distance converges to the true Wasserstein distance when the mesh is refined. This is also the case as far as the time discretization is concerned; one could likely adapt the method of proof of [START_REF] Erbar | Computation of Optimal Transport on Discrete Metric Measure Spaces[END_REF], but doing so is out of the scope of this article.

In Figure 6, however, we present some experiments indicating that convergence under space and time refinement is likely to be true. Number of discretization points in space Error with ground truth in L 1 norm Fig. 6. Top row: the test case. The inputs, i.e. probability distributions at times t = 0 and t = 1, correspond to the same density translated two different ways. Optimal transport predicts that at time t = 1 /2 we should observe the same density again, but translated to the midpoint between the two inputs; this gives us ground truth we can use to verify our algorithm's output. Bottom row: convergence plots. On the left: error, measured in L 1 norm, where the mesh is fixed (regular triangle mesh with 100 points per side of the square) and the number N of discretization points in time varies. On the right: error, measured in L 1 norm, where the number of discretization points in time is fixed (127 points) and the mesh is a regular triangle mesh whose number of points per side varies and is plotted on the x -axis.

t = 0 t = 1 /2 t = 1
These were conducted in the simplest case: translation of a given density on a flat space. For this problem, the ground truth is known, and for a flat space it is clear what it means to refine the mesh: We have use a regular triangle mesh with an increasing number of points per side. The error was evaluated at time t = 1 /2 between the computed geodesic and the ground truth. As a measure of error, as the distributions are compactly supported, we use a total variation norm (in other words the L 1 norm between the densities) rather 250:10 • Hugo Lavenant, Sebastian Claici, Edward Chien, and Justin Solomon than the Kullback-Leibler divergence. As expected, we observe a decrease in error as the temporal and spatial meshes are refined.

Congestion and regularization

In optimal transport there is no price paid for highly-congested densities. Imagine the probability distributions as an assembly of particles moving along the surface. Along a geodesic in Wasserstein space, each particle evolves in time by following a geodesic on the surface-but does not feel the presence of its neighbors. Now imagine, due to the particular structure of the triangle mesh, there is a small shortcut in terms of geodesic distance through which all geodesics tend to concentrate. This is likely to appear near a hyperbolic vertex [START_REF] Polthier | Straightest geodesics on polyhedral surfaces[END_REF]. Then all the particles have the incentive to take this shortcut, resulting in densely-populated zones, as they are not prone to congestion. As an example, see the first row of Figure 7 in which, to go from the left to the right of the armadillo, all the particles go through only two paths, leaving the rest of the mesh without any mass.

This effect, although visually unpleasant, would be observed on a smooth surface M as soon as geodesics concentrate in some regions. A way to remove this artifact is to penalize congestion; we can do so with little modification to the algorithm.

We penalize the densities by their L 2 norms: The choice of the exponent 2 is important, as it preserves the quadratic structure of the optimization problem. Namely, we add to the Lagrangian (20) the term ατ

2 t ∈ G c time v ∈V |v ||µ t v | 2 = sup λ t ∈ G c time v ∈V τ |v | λ t v µ t v - 1 2α (λ t v ) 2 ,
(24) where the parameter α tunes the scale of the congestion effect and λ ∈ R G c time ×V corresponds to the dual variable associated to the congestion constraint.

Using the notation from Section 4.2, one can write the problem as maximizing max

Λ(φ, λ)=q F (φ, λ) + C(q), (25)
but this time

F (φ, λ) = (17) - 1 2α t ∈ G c time v ∈V τ |v |(λ t v ) 2 (26)
and Λ(φ, λ) = (-λ, 0) -Λ(φ). Then one runs exactly the same algorithm, with a straightforward adaptation of the update formulas. After regularization, the interpolation is no longer a geodesic. For instance, the interpolation between two instances of the same probability distribution is not constant in time, because the L 2 norm potentially can be reduced by diffusing outward in the intermediate time steps. On the other hand, undesirable sharp features and oscillations can be removed, as seen in Figure 7. Note that regardless of the level of regularization, the interpolating curves are still valued in P(S), i.e. mass is still preserved along the interpolation.

The tuning of the parameter α allows our method to be robust to noisy mesh inputs, as shown in Figure 8. Noisy meshes have more local variation in curvature, leading to a higher tendency for congested trajectories, but this can be tamed via greater regularization.

Recall that the dynamical formulation of optimal transport can be interpreted as the least-action principle for a pressureless gas. The effect of the penalization of congested densities can be seen, from the modeling point of view, as adding a pressure force: the trajectories of the moving particles are no longer geodesics, they are bent by the pressure forces. The congestion term can also be see as an instance of variational mean field games, for which the augmented Lagrangian approach has been applied for flat spaces with grid discretization [START_REF] Benamou | Variational mean field games[END_REF].

Rather than a drawback, we see the regularization as an added feature of our method. Without regularization, one has a faithful discrete Benamou-Brenier formula on discrete surfaces with Riemannian structure. For applications in physics or gradient flows (Subsection 6.2), this is likely the preferable formulation. For graphics, where blurriness might be sharpened a posteriori, penalizing concentration of mass is reasonable. Either can be achieved thanks to our regularization term without additional computational cost and only by adding a few lines of code.

Intrinsic geometry

To illustrate the fact that the discrete Wasserstein metric is really associated to the geometric structure of the mesh, we perform the following experiment. We design a mesh where the right part is much coarser than the left one, and we let the density evolve. As one can see in Figure 9, the jump in coarseness does not affect the density and does not produce any numerical artifact.

Arbitrary topologies

The discrete formulation that we have chosen applies without change to meshes with boundary and those of non-spherical topology. This is illustrated in Figure 10 with two meshes topologically equivalent to a disc and a torus.

In the first example, the interpolating distribution stays near the boundary, approximately following the geodesic between the means of the endpoint distributions. In the second example, one can see the initial distribution splitting to travel both ways to the other side of a handle, before merging again to achieve the final distribution. [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF] provide a convolutional method for approximating the Wasserstein geodesic between two distributions supported on triangle meshes. Their approach solves a regularized optimal transport barycenter problem using a modified Sinkhorn algorithm, with a heat kernel taking the place of explicitly-calculated pairwise distances between vertices. As a result, their method blurs the input distributions, and the interpolated distributions are typically of higher entropy than the endpoints. This is combated with a nonconvex projection method that attempts to lower the entropy of intermediate distributions to an approximated bound.

Comparison to convolutional method

In comparing our methods, we found that [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]] also tends to produce interpolating distributions that do not travel with constant speed. This effect can be seen in Figures 11 and12, where their interpolating distributions remain mostly stationary for times near t = 0 and t = 1, but move with high speed for times near t = 1 /2. Loosening the entropy bound in the nonconvex step helps somewhat, but the problem persists regardless. Most likely this effect is due to the fact that the entropy reduction step of their algorithm is not geometry-aware but rather simply sharpens the regularized interpolant. Our method does not suffer from this issue, and the spread of our interpolating distributions is comparable or better in both cases. Furthermore, unless the regularizer α is large, our interpolating distributions tend to diffuse only in the direction of the geodesics along which particles are traveling, which better mimics the behavior of Wasserstein geodesics; this diffusion is reduced by adding more time steps to our interpolation problem.

α = 0 α = 10 -2 α = 10 -1 α = 1 t = 0 t = 1 /4 t = 1 /2 t = 3 /4 t = 1
Our formulation also has comparable runtimes to the convolutional method of [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]. For instance, the implementation of the convolutional method provided by the authors of that paper took 57 and 141 seconds to converge, on the punctured sphere (1020 vertices) and teapot (3900 vertices), respectively, for 13 time steps. This is to be put in comparision with the timings provided in Subsection 5.2.

The comparisons in this section were computed on a 3.60GHz Intel i7-7700 processor with 32GB of RAM. For the convolutional method, the heat kernel was used to diffuse to t = 0.0015 with 10 implicit Euler steps.

APPLICATIONS AND EXTENSIONS 6.1 Harmonic mappings

As P(M) can be viewed as a Riemannian manifold of infinite dimension, one can consider not only geodesics valued in this space, but also harmonic mappings. That is, we consider a domain Ω and a function µ : x ∈ Ω → µ x which takes fixed values on ∂Ω the boundary of Ω and minimizes the Dirichlet energy

Dir(µ) := 1 2 ∫ Ω ∥∇ Ω µ x ∥ 2 T µ x P(M) dx, ( 27
)
where the norm ∥ ∥ µ x is defined in (5). Such harmonic mappings have been introduced under the name soft maps by [START_REF] Solomon | Dirichlet energy for analysis and synthesis of soft maps[END_REF]2012] for the purpose of surface mapping; one can also find a formal definition and theoretical analysis in [START_REF] Lavenant | Harmonic mappings valued in the Wasserstein space[END_REF][START_REF] Lu | Properties of Soft Maps on Riemannian Manifolds[END_REF]].

As explained by [START_REF] Lavenant | Harmonic mappings valued in the Wasserstein space[END_REF], if some boundary conditions μ : ∂Ω → P(M) are given, the Dirichlet problem consists in minimizing the Dirichlet energy ( 27 on ∂Ω. Specifically, it is a convex problem whose dual reads

max φ ∫ ∂Ω ∫ M φ(x, •) • n Ω (x)d μx dx s.t. ∇ Ω • φ + 1 2 ∥∇ M φ ∥ 2 ⩽ 0 on Ω × M, (28)
where φ is defined on Ω × M and valued in T Ω (i.e. for a point (x, y) ∈ Ω ×M, one has φ(x, y) ∈ T x Ω), n Ω (x) is the outward normal to Ω, and dx is the integration on ∂Ω w.r.t. the surface measure. In the case where Ω is a segment, the dual problem (3) for the geodesics is recovered.

To discretize (28), we use the same strategy as for the geodesic problem. We assume that we have S Ω = (V Ω , E Ω ,T Ω ) a triangulation of the surface Ω. The discrete unknown φ maps every element of T Ω ×V onto a vector in R 3 (thought as the tangent space of S Ω ). The divergence ∇ Ω • φ is replaced by its discrete counterpart which lives on V Ω × V . On the other hand, ∇ M φ is naturally seen as a vector in R 3 for each pair of triangles in T Ω × T . We apply the same idea: For the term ∥∇ M φ ∥ 2 , first square and then average (weighting by the area of the triangles) to put it on the grid V Ω × V . Once we have a fully-discrete problem, we build an augmented Lagragian and use ADMM: The solution µ is the Lagrange multiplier associated to the constraint A = ∇ Ω • φ.

In Figure 13, we show an example where S Ω is a triangulation of an equilateral triangle and S the triangulation of a flat square. On the corners of the triangle we put some distributions, and on the side, as part of the boundary data, we have chosen to put the geodesic in the Wasserstein space between the distributions on the corners. This choice is arbitrary, we could have chosen other configurations on the edges of the boundary of the triangle Ω. This picture resembles the barycentric interpolation [Benamou et al. 2015;[START_REF] Cuturi | Fast computation of Wasserstein barycenters[END_REF][START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF], though no theoretical evidence indicates that harmonic and barycentric interpolation coincide.

Gradient flows in the Wasserstein space

The seminal work of [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], demonstrates that by considering the gradient flow

∂ t µ t = -∇ P(M ) F (µ t ), (29)
where ∇ P(M ) is the gradient w.r.t. the scalar product defined by ( 4) and ( 5), one can recover several well-known PDE (Fokker-Planck, porous medium equations, aggregation diffusion equations) by an appropriate choice of the functional F : P(M) → R [START_REF] Ambrosio | Gradient Flows: In Metric Spaces and in the Space of Probability Measures[END_REF][START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]]. Moreover, a natural implicit discretization comes with this equation: The idea is to take a time step s > 0 and

t = 0 mesh t = 1 t = 1 /4 t = 1 /2 t = 3 /4
Fig. 9. Top row: mesh and initial/final probability distributions. Notice the difference of coarseness in the mesh. Bottom row: interpolation shown at different times where no effect of the difference in coarseness is seen. We have used the regularization described in Subsection 5.4 with α = 0.1.

t = 0 t = 1 /3 t = 2 /3 t = 1
Fig. 10. Our formulation easily handles non-spherical topologies. In the top row is a punctured sphere, and in the bottom row is a genus-1 teapot mesh. These interpolations were generated with α = 0.02 and α = 0.2, respectively.

to define µ ks for k = 0, 1, . . . recursively in the following way:

µ (k+1)s := arg min µ W 2 2 (µ, µ ks ) 2s + F (µ) . (30)
The scheme above, defined in arbitrary metric spaces, is refered to as minimizing movement scheme, in the framework of optimal transport it is sometimes known as a JKO integrator, after the work of [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF].

In the case of a discrete surface S with the Riemannian structure on P(S), (29) makes sense and can be written

∂ t µ t = -(P µ t ) -1 ∇F (µ t ) (31)
where P µ is the metric tensor defined in ( 14) and ∇F is the usual gradient of F as a function from R |V | to R. We can still define an iterative scheme to compute solutions of (31) by solving (30), where W 2 has been replaced by the discrete Wasserstein distance W d . As long as F : P(S) → R is convex, the discrete version of (30) can be tackled by the same augmented Lagrangian and ADMM iterations at the price of introducing an additional variable [START_REF] Benamou | An augmented Lagrangian approach to Wasserstein gradient flows and applications[END_REF]].

We cannot recover cotangent Laplacian heat flow as a gradient flow for the Wasserstein distance W d ; to get such a result, one would need to define μ by some nonlinear averaging process rather than (10). Such a choice would increase the complexity of computing geodesics, and it would be likely to introduce more diffusion as in [START_REF] Erbar | Computation of Optimal Transport on Discrete Metric Measure Spaces[END_REF][START_REF] Maas | Gradient flows of the entropy for finite Markov chains[END_REF].

The advantages of our numerical method are that positivity is automatic and that mass is preserved. Moreover, as we expect the difference between two solutions µ ks and µ (k+1)s to be very small, we do not need a large number of discretization points in time N (we chose N = 5 in practice) for the computation of the discrete Wasserstein distance.

We apply our model to two different cases, illustrated in Figure 14. The first corresponds to

F (µ) = v ∈V |v |W v µ v if µ v ⩽ µ ⋆ ∀v ∈ V , +∞ else, ( 32
)
where W contains one value per vertex. This choice of F yields a crowd motion model [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Santambrogio | Crowd motion and evolution PDEs under density constraints[END_REF]]:

The probability distribution wants to flow to the areas where W is low, but at the same time its density is constrained to stay below a threshold µ ⋆ . This model can only be formulated in terms of a gradient flow in Wasserstein space and not as an evolution PDE, justifying the use of a JKO integrator. The second corresponds to

F (µ) = 1 m -1 v ∈V |v |(µ v ) m . ( 33
)
In the case of a continuous manifold M, this choice of F yields the porous medium equation ∂ t µ = ∆(µ m ) [Vázquez 2007]. With the scheme (30), we do not recover a discrete cotangent Laplacian, but the computed solution still exhibits features of the porous medium, like a finite speed of propagation and a convergence to a uniform probability distribution.

DISCUSSION AND CONCLUSION

Although techniques using entropic regularization or semi-discrete optimal transport can interpolate between distributions on a discrete surface, they do not provide a Riemannian structure and are subject to practical limitations that restrict the scenarios to which they can be applied. Using an intrinsic formulation of dynamical transport, we can realize the theoretical and practical potential of optimal transport on discrete domains enabled by the Riemannian structure on the space of probability distributions, the so-called Otto calculus.

Our technique can be phrased in familiar language from discrete differential geometry and is implementable using standard tools in that domain. The key ingredients, namely first-and second-order operators in geometry processing (gradient, divergence, Laplacian) as well as SOCP optimization, remain in the realm of what is already widely used.

We have demonstrated the power of our model by showing how it can handle a variety of geometries and peaked distributions, while introducing little diffusion. Mass may concentrate to yield a visually inelegant result, but this behavior is at the core of optimal transport theory and expected: No price is paid for mass congestion, and hence any concentration of geodesics will result in a concentration of mass. Nevertheless, as we have shown, one can easily modify the optimization problem to penalize congested densities, leading to smoother interpolants with a controllable level of diffusion. Unlike entropically-regularized transport, however, our optimization problem does not degenerate as the coefficient in front of the regularizer vanishes.

Beyond evaluation of transport distances, our framework extends to support other tasks involving transport terms. We can reliably compute harmonic mappings valued in this discrete Wasserstein space, and the JKO integrator based on our discrete Wasserstein distance exhibits expected qualitative behavior.

The main drawback of our approach remains its scalability. The bottleneck of the computations is the solution of a linear system t = 0 t = 1 t = 2 t = 3 t = 4 t = 5 t = 6 t = 7 whose number of unknowns is the product of the number of discretization points in time and the size of the mesh. This is an extremely structured linear system on a product manifold, for which specialized matrix inversion techniques may exist. In any event, with the current bottleneck our method can handle meshes with few thousand vertices but is not currently practical for larger meshes.

As one of the first structure-preserving discretizations of transport on meshes, our work also suggests several exciting avenues for future research. Many theoretical properties of our discrete Wasserstein distance remain to be explored. For instance, while we have shown that our formulation is a true Riemannian distance, one could verify the extent to which a wealth of other theoretical properties of transport are preserved. Convergence of our transport over meshes to the true transport in the limit of mesh refinement also remains an open problem for our techniques and others in a similar class. From a practical perspective, a natural next step is to accelerate the optimization procedure as much as possible; a faster solver for the convex optimization problem would clearly benefit our method.

Fig. 2 .

 2 Fig. 2. Top row: Interpolation of probability distributions. The left and right distributions are data and the middle ones are the output of our algorithm. Bottom row: Display of the momentum m = µv, where v is a time-dependent velocity-field advecting the left distrbution on the right one. We have used the regularization described in Subsection 5.4 with α = 0.1.

Fig. 3 .

 3 Fig. 3. Top row: schematic view of the static formulation of optimal transport

Fig. 4 .

 4 Fig. 4. Left: temporal grids G c time and G st time for N = 4. Right: a vertex (•) surrounded by 6 adjacent triangles, the dual barycentric cell is in gray.

  is a natural discretization of the time derivative. Next, we discretize φ ∈ R G st time × |V | a function depending both on space and time. The constraint ∂ t φ t v be imposed on the centered grid G c time . It is enough to replace ∂ t φ by Dφ. On the other hand, the term 1 2

  G st time ×V as a (N + 1) × |V | matrix (with rows indexed by time and columns by space), the equation 2018-09-17 17:48. Page 7 of 1-16. Algorithm 1 Geodesic Computation function Geodesic( μ0 , μ1 ) Initialize φ, A, B, µ, m ← 0 while PrimalResidual and DualResidual > ε do φ ← solution of (21) for s, v ∈ G c time × V do update A and B by solving (22) end for Update µ and m through (23) end while return µ end function satisfied by a maximizer of L over φ reads

  Fig.5. Amplitude of the primal and dual residual[Boyd et al. 2011, Section 3.3] in L 2 norm. The distributions μ0 and μ1 are delta functions located on respectively the right and left hand of the armadillo. We also show the midpoint µ 1/2 for different numbers of iterations (10,50 and 5000). After a few hundred iterations, there is no visible difference in µ 1/2 . There is a jump in the value of the dual residual at around 4600 iterations. It is due to a change in the value of the parameter r , which is updated according to the heuristic rule presented in Section 3.4.1 of[START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]].

Fig. 7 .

 7 Fig. 7. Effect of the regularizing parameter α penalizing congestion. On each row, the interpolation between the same boundary data (distributions located on the right and on the left of the armadillo) is shown. Different rows correspond to different to different values of α . The color in each image is normalized independently from the others, explaining the change in intensity. Mass is always preserved along the interpolation.

  ) of µ under the constraint that µ = μ 2018-09-17 17:48. Page 11 of 1-16. ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 250:12 • Hugo Lavenant, Sebastian Claici, Edward Chien, and Justin Solomon t = 0 t = 1 /5 t = 2 /5 t = 3 /5 t = 4 /5 t = 1

Fig. 8 .

 8 Fig. 8. Robustness to noisy meshes, after adjusting the parameter α . Top row: original mesh, α = 0.02; middle row: noisy mesh, α = 0.1; bottom row: very noisy mesh, α = 0.2. The bounding boxes of the meshes were of side length ~1.5. Noisy mesh vertices were obtained by uniformly random perturbation, in the normal direction, of magnitudes up to 0.02 and 0.04, for the middle and bottom row, respectively.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Constant-speed interpolation. Indicator distributions on handle ends of a pliers mesh are interpolated. Top row: our method, calculated with α = 0.001; middle row: method of Solomon et al. [2015], calculated with entropy bounded by that of the endpoint distributions; bottom row: method of Solomon et al. [2015], calculated with no entropy bound. As can be seen, the method of Solomon et al. [2015] stays mostly stationary except for the middle frames.

  Fig. 13. (a) Data: boundary conditions, i.e. value of the probability distributions at the boundary vertexes of a regular triangulation of an equilateral triangle. (b) Result: interpolation over the interior of the triangular mesh obtained by minimizing the Dirichlet energy with the boundary constraints. The mesh of the source domain Ω is a regular triangulation of a triangle with 15 vertexes per side but for clarity reasons we display the value of the probability distributions only a subset of the set of vertexes. The target domain M is a flat square triangulated with 30 vertexes by side.

Fig. 14 .

 14 Fig.14. First row: gradient flow in discrete Wasserstein space of the functional (32).The potential is chosen to be W (x, y, z) = z: the mass flows down and then saturates because of the density constraint. Congestion can be observed as the mass goes through the thinner part. Second row: gradient flow of the internal energy (33) starting from a delta function on the surface of the plier. The colors do not correspond to the same scale for different times; otherwise nothing would be seen for t > 0 thanks to the peaked distribution at t = 0. Notice the finite speed of propagation of the density.

  2018-09-17 17:48. Page 13 of 1-16. ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 250:14 • Hugo Lavenant, Sebastian Claici, Edward Chien, and Justin Solomon

	t = 0	t = 1 /7	t = 2 /7	t = 3 /7	t = 4 /7	t = 5 /7	t = 6 /7	t = 1

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 2 of 1-16.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 4 of 1-16.

If we do the opposite (averaging and then taking the square), there are spurious modes in the kernel of the quadratic part of the constraint, which leads to poor results when working with non-smooth densities.

2018-09-17 17:48. Page 5 of 1-16. ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 6 of 1-16.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 8 of 1-16.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 10 of 1-16.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 12 of 1-16.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 14 of 1-16.

ACM Trans. Graph., Vol. 37, No. 6, Article 250. Publication date: November 2018. 2018-09-17 17:48. Page 16 of 1-16.

ACKNOWLEDGEMENTS

J. Solomon acknowledges the generous support of Army Research Office grant W911NF-12-R-0011 ("Smooth Modeling of Flows on Graphs"), of National Science Foundation grant IIS-1838071 ("BIG-DATA:F: Statistical and Computational Optimal Transport for Geometric Data Analysis"), from the MIT Research Support Committee, from an Amazon Research Award, from the MIT-IBM Watson AI Laboratory, and from the Skoltech-MIT Next Generation Program. Most of this work was done during a visit of H. Lavenant to MIT; the hospitality of CSAIL and MIT is warmly acknowledged. The authors thank the reviewers for their helpful feedback.