Giuseppe Buttazzo

AND Guillaume Carlier
email: carlier@ceremade.dauphine.fr

Serena Guarino
email: sguarino@mail.dm.unipi.it

Lo O Bianco

Optimal regions for congested transport

Keywords: Shape optimization, transport problems, congestion effects, optimal networks. 2010 Mathematics Subject Classification: 49Q20, 49Q10, 90B20

come   L'archive ouverte pluridisciplinaire

Introduction

As everybody has experienced while traveling in urban traffic, the planning of an efficient network of roads is an extremely complex problem, involving many parameters such as the distribution of residences and working places, the period of the day one travels, the attitude of drivers, . . . . The congestion effects are also to be taken into account, since they are responsible for the formation of traffic jams and and have a social cost in terms of waste of time.

In the present paper we consider a very simplified model in which the densities of residents and of working places are known, represented by two probability measures f + and f -. The congestion effects in mass transportation theory has been deeply studied in the literature; we refer for instance to [START_REF] Carlier | Optimal transportation with traffic congestion and Wardrop equilibria[END_REF][START_REF] Brasco | On certain anisotropic elliptic equation arising in congested optimal transport: local gradient bounds[END_REF] and references therein. Denoting by f the difference f = f + -f -and by σ the traffic flux, the model, in the stationary regime, reduces to a minimization problem of the form min Ω H(σ) dx : -div σ = f in Ω, σ • n = 0 on ∂Ω .

(1.1)

Here Ω is the urban region under consideration, a bounded Lipschitz domain of R d , the boundary conditions at ∂Ω are usually taken imposing zero normal flux σ • n = 0, and H : R d → [0, +∞] is the congestion function, a convex nonnegative function with lim |s|→+∞ H(s) = +∞. The first order PDE -div σ = f in Ω, σ • n = 0 on ∂Ω has to be intended in the weak sense σ, ∇φ = f, φ for every φ ∈ C ∞ (Ω)

and it captures the equilibrium between the traffic flux σ and the difference between supply and demand f . In the case H(s) = |s| no congestion effect occurs, and the transport problem reduces to the Monge's transport, where mass particles travel along geodesics (segments in the Euclidean case). As it is well known, in the Monge's case the integral cost above is finite for every choice of the probabilities f + and f -. On the contrary, when H is superlinear, that is lim |s|→+∞ H(s) |s| = +∞, congestion effects may occur and the mass particles trajectories follow more complicated paths. In this case the integral cost can be +∞ if the source and target measures f + and f -are singular. For instance, if the congestion function H has a quadratical growth, in order to have a finite cost it is necessary that the signed measure f = f + -f -be in the dual Sobolev space H -1 ; thus, if d > 1 and the measures f + or f -contain some Dirac mass, the minimization problem (1.1) is meaningless. In other words, superlinear congestion costs prevent too high concentrations.

In the present paper, we aim to address the efficient design of low-congestion regions; more precisely, two congestion functions H 1 and H 2 are given, with H 1 ≤ H 2 , and the goal is to find an optimal region C ⊂ Ω where the less congested traffic may travel. Since reducing the congestion in a region C is costly (because of roads improvement, traffic devices, . . . ), a term m(C) will be added, to describe the cost of improving the region C, then penalizing too large low-congestion regions. On the region Ω \ C we then have the normally congested traffic governed by the function H 2 , while on the low-congestion region C the traffic is governed by the function H 1 . Throughout the paper, we will assume that H 1 and H 2 are two continuous convex functions such that 0

≤ H 1 ≤ H 2 and lim |s|→+∞ H i (s) |s| = +∞, i = 1, 2.
The mathematical formulation of the problem then is as follows, for every region C we may consider the cost function

F (C) = min Ω\C H 2 (σ) dx + C H 1 (σ) dx : -div σ = f in Ω, σ • n = 0 on ∂Ω , (1.2)
so that the optimal design of the low-congestion region amounts to the minimization problem min

F (C) + m(C) : C ⊂ Ω . (1.3)
Several cases will be studied in the sequel, according to the various constraints on the lowcongestion region C and the corresponding penalization/cost m(C). The simplest case is when C is a d-dimensional subdomain of Ω and the penalization m(C) involves the perimeter of C: in this situation an optimal region C is shown to exist and a necessary optimal condition is established.

When m(C) is simply proportional to the Lebesgue measure of C (that we will denote by L d (C) or by |C|), on the contrary an optimal domain C may fail to exist and a relaxed formulation of the problem has to be considered; in this case the optimal choice for the planner is to have a low-congestion area C 0 , a normally congested area C 1 , together with an area Ω \ (C 0 ∪ C 1 ) with intermediate congestion (that is a mixing of the two congestion functions occurs). For this case, we also give some numerical simulations in dimension two.

Another class of problems arises when the admissible sets C are networks, that is closed connected one-dimensional sets. In this case the penalization cost m(C) is proportional to the total lenght of C (the 1-dimensional Hausdorff measure H 1 (C)). In this case the precise definition of the cost function F (C) in (1.2) has to be reformulated more carefully in terms of measures (see Section 4). This one-dimensional problem has been extensively studied in the extremal case where H 1 = 0 and H 2 (s) = |s| (see for instance [START_REF] Buttazzo | Stationary configurations for the average distance functional and related problems[END_REF][START_REF] Buttazzo | Optimal transportation problems with free Dirichlet regions[END_REF][START_REF] Buttazzo | Optimal urban networks via mass transportation[END_REF][START_REF] Buttazzo | Optimal transportation networks as free Dirichlet regions for the Monge-Kantorovich problem[END_REF][START_REF] Buttazzo | Calculus of Variations: topics from the mathematical heritage of E. De Giorgi[END_REF]) providing an interesting geometric problem called average distance problem; an extended review on it can be found in [START_REF] Lemenant | A presentation of the average distance minimizing problem[END_REF]. We also point out that a similar problem arises in some models for the mechanics of damage, see for instance [START_REF] Braides | Evolution of damage in composites: the one-dimensional case[END_REF].

Perimeter constraint on the low-congestion region

In this section we consider the minimum problem (1.3), where the cost F (C) is given by (1.2) and m(C) = k Per(C), being k > 0 and Per(C) the perimeter of the set C in the sense of De Giorgi (see for instance [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]). Thanks to the coercivity properties of the perimeter with respect to the L 1 convergence of the characteristic functions (that we still call L 1 convergence of sets), we have the following existence result. Theorem 2.1. Assume that the cost F (C) is finite for at least a subset C of Ω with finite perimeter and that m(C) = k Per(C) with k > 0. Then there exists at least an optimal set C opt for problem (1.3).

Proof. Let (C n ) n∈N be a minimizing sequence for the optimization problem (1.3); then sequence Per(C n ) is bounded. Thanks to the compactness of the embedding of BV into L 1 , we may extract a (not relabeled) subsequence converging in L 1 to a subset C of Ω. We claim that this set C is an optimal set for the problem (1.3). Indeed, for the properties of the perimeter we have

Per(C) ≤ lim inf n Per(C n ).
Moreover, if we denote by σ n an optimal (or asymptotically optimal) function for

F (C n ) = Ω\Cn H 2 (σ n ) dx + Cn H 1 (σ n ) dx with -div σ n = f in Ω, σ n
• n = 0 on ∂Ω, by the superlinearity assumption on the congestion functions H 1 and H 2 , and by the De La Vallée Poussin compactness theorem, we have that (σ n ) n∈N is compact for the weak L 1 convergence and so we may assume that σ n weakly converges in L 1 (Ω) to a suitable function σ. This function σ still verifies the condition

-div σ = f in Ω, σ • n = 0 on ∂Ω.
Thanks to the convexity of H 1 and H 2 and the strong-weak lower semicontinuity theorem for integral functionals (see for instance [START_REF] Buttazzo | Semicontinuity, Relaxation and Integral Representation in the Calculus of Variations[END_REF]), we have

F (C) ≤ lim inf n F (C n ).
Therefore the set C is optimal and the proof is concluded.

Our aim now is to establish optimality conditions not only on an optimal flow σ but also on the corresponding optimal low-congestion regions C. Optimality conditions for σ can be directly derived from the duality formula:

F (C) = inf σ∈Γ f C H 1 (σ) dx + Ω\C H 2 (σ) dx = -inf u C H * 1 (∇u) dx + Ω\C H * 2 (∇u) dx - Ω uf dx ,
from which one easily infers that

σ = σ int in C, σ ext in Ω \ C where σ int = ∇H * 1 (∇u int ) in C, σ ext = ∇H * 2 (∇u ext ),
in Ω \ C the minimizer u in the dual is then given by:

u = u int in C, u ext in Ω \ C .
We have used the notations σ int , σ ext , u int and u ext to emphasize the fact that σ and ∇u may have a discontinuity when crossing ∂C. It is reasonable (by elliptic regularity and assuming smoothness of C) to assume that σ and ∇u are Sobolev on C and Ω \ C separately but they are a priori no better than BV on the whole of Ω. The functions u int and u ext are then at least formally characterized by the Euler-Lagrange equations

-div ∇H * 1 (∇u int ) = f in C, -div ∇H * 2 (∇u ext ) = f, in Ω \ C together with ∇H * 1 (∇u int ) • n = 0, on ∂Ω ∩ C, ∇H * 2 (∇u ext ) • n = 0, on ∂Ω ∩ Ω \ C,
and (assuming that f does not give mass to ∂C) the continuity of the normal component of σ across ∂C:

∇H * 1 (∇u int ) -∇H * 2 (∇u ext ) • n C = 0, on ∂C ∩ Ω,
where n C denotes the exterior unit vector to C. Now, we wish to give an extra optimality condition on C itself assuming that is smooth. To do so, we take a smooth vector field V such that V • n = 0 on ∂Ω, and we set C t = ϕ t (C), where ϕ t denotes the flow of V (i.e. ϕ 0 = id, ∂ t ϕ t (x) = V (ϕ t (x))). For t > 0, we then have

0 ≤ 1 t [F (C t ) -F (C) + k Per(C t ) -k Per(C)]. (2.1)
As for the perimeter term, it is well-known (see for instance [START_REF] Henrot | Variation et Optimisation de Formes[END_REF]) that the first-variation of the perimeter involves the mean curvature H of ∂C, more precisely, we have:

d dt Per(C t ) t=0 = ∂C H V • n C dH d-1 . (2.2)
For the term involving H, we observe that

F (C t ) -F (C) ≤ Ct H 1 (σ) dx - C H 1 (σ) dx + Ω\Ct H 2 (σ) dx - Ω\C H 2 (σ) dx.
At this point, we have to be a little bit careful because of the discontinuity of σ at ∂C, but distinguishing the part of ∂C on which V • n C > 0 that is moved outside C by the flow, and that on which V • n C < 0 that is moved inside C by the flow, and arguing as in Theorem 5.2.2 of [START_REF] Henrot | Variation et Optimisation de Formes[END_REF], we arrive at:

lim sup t→0 F (C t ) -F (C) t ≤ ∂C H 1 (σ ext ) -H 2 (σ ext ) (V • n C ) + + H 2 (σ int ) -H 1 (σ int ) (V • n C ) -dH d-1 .
(

Combining (2.1), (2.2) and (2.3), we obtain

0 ≤ ∂C H 1 (σ ext ) -H 2 (σ ext ) + kH (V • n C ) + + H 2 (σ int ) -H 1 (σ int ) -kH (V • n C ) -dH d-1 .
But since V is arbitrary, we obtain the extra optimality conditions:

H 2 (σ int ) -H 1 (σ int ) ≥ kH ≥ H 2 (σ ext ) -H 1 (σ ext ) on ∂C ∩ Ω which, since H 2 ≥ H 1 , in particular implies that ∂C has nonnegative mean curvature.
The regularity of ∂C is an interesting open question. Note that when d = 2 and Ω is convex, replacing C by its convex hull diminishes the perimeter and also the congestion cost, so that optimal regions C are convex, this is a first step towards regularity, note also that convexity of optimal regions is consistent with the curvature inequality above.

Let us illustrate the previous conditions on the simple quadratic case where

H 1 (σ) = a 2 |σ| 2 , H 2 (σ) = b 2 |σ| 2 with 0 < a < b.
The optimality conditions for the pair u, σ then read as

-a∆u int = f in C -b∆u ext = f in Ω \ C, ∂u ∂n = 0 on ∂Ω, σ int = ∇u int a σ ext = ∇uext b , together with ∇u int a - ∇u ext b • n C = 0 on ∂C ∩ Ω and b -a 2 |σ int | 2 = b -a 2a 2 |∇u int | 2 ≥ kH ≥ b -a 2 |σ ext | 2 = b -a 2b 2 |∇u ext | 2 on ∂C ∩ Ω
where H again denotes the mean curvature of ∂C.

Relaxed formulation for the measure penalization

In this section we consider the case when the penalization on the low-congestion region is proportional to the Lebesgue measure, that is m(C) = k|C| with k > 0. The minimization problem we are dealing with then becomes min

σ,C C H 1 (σ) dx + Ω\C H 2 (σ) dx + k|C| : σ ∈ Γ f (3.1)
where

Γ f = σ ∈ L 1 (Ω; R d ) : div σ = f in Ω, σ • n = 0 on ∂Ω .
Passing from sets C to density functions θ with 0 ≤ θ(x) ≤ 1 we obtain the relaxed formulation of (3.1)

min σ,θ Ω θH 1 (σ) dx + Ω (1 -θ)H 2 (σ) dx + k Ω θ dx : σ ∈ Γ f . (3.2)
Writing the quantity to be minimized as

Ω H 2 (σ) + θ H 1 (σ) + k -H 2 (σ) dx the minimization with respect to θ is straightforward; in fact, if H 1 (σ) + k > H 2 (σ) we take θ = 0, while if H 1 (σ) + k < H 2 (σ) we take θ = 1.
In the region where H 1 (σ) + k = H 2 (σ) the choice of θ is irrelevant. In other words, we have taken

θ = 1 {H 1 (σ)+k<H 2 (σ)} ,
which gives

H 2 + θ H 1 + k -H 2 = H 2 -H 1 + k -H 2 -= H 2 ∧ H 1 + k .
Therefore, in the relaxed problem (3.2) the variable θ can be eliminated and the problem reduces to min

Ω H 2 (σ) ∧ H 1 (σ) + k dx : σ ∈ Γ f . (3.3)
Clearly the infimum in (3.3) coincides with that of (3.1) but since the new integrand H 2 ∧ H 1 + k is not convex, a further relaxation with respect to σ is necessary. This relaxation issue with a divergence constraint has been studied in [START_REF] Braides | Relaxation of functionals with constraint on the divergence[END_REF], where it is shown that the relaxation procedure amounts to convexify the integrand. We then end up with the minimum problem min

Ω H 2 (σ) ∧ H 1 (σ) + k * * dx : σ ∈ Γ f (3.4)
where * * indicates the convexification operation. Recalling that H 1 and H 2 are superlinear, we have that:

-in the region where

H 2 ∧ H 1 + k * * (σ) = H 2 (σ)
we take θ = 0. In other words, in this region, it is better not to spend resources for improving the traffic congestion; -in the region where

H 2 ∧ H 1 + k * * (σ) = H 1 (σ) + k
we take θ = 1. In other words, in this region, it is necessary to spend a lot of resources for improving the traffic congestion; -in the region where

H 2 ∧ H 1 + k * * (σ) < H 2 ∧ H 1 + k (σ)
we have 0 < θ(x) < 1 so that there is some mixing between the low and the high congestion functions. In other words, in this region the resources that are spent for improving the traffic congestion are proportional to θ.

The previous situation is better illustrated in the case where both functions H 1 and H 2 depend on |σ| and H 2 -H 1 increases with |σ|. In this case, we denote by r 1 the maximum number such that

H 2 ∧ H 1 + k * * (r) = H 2 (r)
and by r 2 the minimum number such that

H 2 ∧ H 1 + k * * (r) = H 1 (r) + k, then we have θ(x) = |σ| -r 1 r 2 -r 1 whenever r 1 < |σ| < r 2 .
In this case, for small values of the traffic flow (|σ| ≤ r 1 ), it is optimal not to spend any resource to diminish congestion, on the contrary when traffic becomes large (|σ| ≥ r 2 ), it becomes optimal to reduce the congestion to H 1 . Finally, for intermediate values of the traffic, mixing occurs with the coefficient θ above as a result of the relaxation procedure. Also, problem (3.4) is of type (1.1) and it is well-known, by convex analysis, that we have the dual formulation min

Ω H(σ) dx : σ ∈ Γ f = sup Ω u df - Ω H * (∇u) dx = -inf Ω H * (∇u) dx - Ω u df , (3.5)
where H(σ) = (H 2 (σ) ∧ (H 1 (σ) + k)) * * . Notice that the Euler-Lagrange equation of problem (3.5) is formally written as

-div ∇H * (∇u) = f in Ω ∇H * (∇u) • ν = 0 on ∂Ω. (3.6)
Moreover, the link between the flux σ and the dual variable u is σ = ∇H * (∇u).

In our case, the Fenchel tranform is easy computed and we have:

H * (ξ) = H * 2 (ξ) ∨ (H * 1 (ξ) -k).
As a conclusion of this paragraph, we observe that the treatment above is similar to the analysis of two-phase optimization problems. This consists in finding an optimal design for a domain that is occupied by two constituent media with constant conductivities α and β with 0 < α < β < +∞, under an objective function and a state equation that have a form similar to (3.5) and (3.6). We refer to [START_REF] Bucur | Variational Methods in Shape Optimization Problems[END_REF] (and references therein) for a general presentation of shape optimization problems and to [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] for a complete analysis of two-phase optimization problems together with numerical methods to treat them.

low-congestion transportation networks

In this section, our main unknown is a one-dimensional subset Σ of Ω; we consider a fixed number r > 0 and the low-congestion regions of the form C Σ,r = x ∈ Ω : dist(x, Σ) ≤ r = Σ r ∩ Ω, where Σ r := Σ + B r (0). and Σ is required to be a closed subset of Ω such that H 1 (Σ) < +∞. The penalization term m(C Σ,r ) is taken proportional to the Lebesgue measure of C Σ,r , so that our optimization problem becomes min

σ,Σ C Σ,r H 1 (σ) dx + Ω\C Σ,r H 2 (σ) dx + k|C Σ,r | : σ ∈ Γ f (4.1)
with k > 0. A key point in the existence proof below consists in remarking that the perimeter of an r-enlarged set Σ r can be controlled by its measure (see Appendix A). It also worth remarking that Σ r has the uniform interior ball of radius r property; for every x ∈ Σ r there exists y ∈ R d such that |x -y| ≤ r and B r (y) ⊂ Σ r . Clearly, r-enlarged sets have the uniform interior ball of radius r property and sets with this property are r-enlarged sets (i.e. can be written as the sum of a closed set and B r (0)), we refer to [START_REF] Alvarez | Existence and uniqueness for dislocation dynamics with nonnegative velocity[END_REF] for more on sets with the uniform interior ball property, and in particular estimates on their perimeter (which we could have used instead of the more elementary Lemma in Appendix A).

Proposition 4.1. Ler r > 0 be fixed, d = 2 and assume that F (C Σ,r ) < +∞ for some closed one-dimensional subset Σ of Ω. Then the optimization problem (4.1) admits a solution.

Proof. The sets C Σ,r satisfy the inequality (see for instance Appendix A)

Per(C Σ,r ) ≤ K r |C Σ,r |
for a suitable constant K depending only on the dimension d. Therefore, for a minimizing sequence (Σ n ) n∈N , the sets C n := C Σn,r = Σ r n ∩ Ω are compact in the strong L 1 convergence, we can thus extract a (not relabeled) subsequence such that C n converges strongly in L 1 (and a.e.) to some C. One can then repeat the proof of Theorem 2.1, to obtain

F (C) + k|C| ≤ inf (4.1).
It only remains to show that C can be obtained as C = C Σ,r (up to a negligible set) for some closed subset of Ω, Σ such that H 1 (Σ) < +∞. Up to an extraction, one can assume that Σ r n converges for the Hausdorff distance to some compact set E (which also satisfies the uniform interior ball property of radius r). Let us first check that C = E ∩ Ω (up to a negligible set), the inclusion C ⊂ E ∩ Ω is standard (see for instance [START_REF] Henrot | Variation et Optimisation de Formes[END_REF]). To prove the converse inclusion, it is enough to show that |C| = |E ∩ Ω| i.e. |C n | → |E ∩ Ω| as n → ∞. For this, we observe that

|C n | -|E ∩ Ω| ≤ |Σ r n \ E| + |E \ Σ r n | The convergence of |Σ r
n \ E| to 0 easily follows from the Hausdorff convergence of Σ r n to E and the fact that E is closed (see [START_REF] Henrot | Variation et Optimisation de Formes[END_REF] for details). As for the convergence of |E \ Σ r n | to 0, we proceed as follows: let ε > 0 and n be large enough so that E ⊂ Σ r n + B ε (0) = Σ r+ε n , using the coarea formula, as in Appendix A, and the fact that the sets Σ r+ε n have bounded perimeter (see again Appendix A or [START_REF] Alvarez | Existence and uniqueness for dislocation dynamics with nonnegative velocity[END_REF]), we get for n large enough

|E \ Σ r n | ≤ |Σ r n + B ε (0) \ Σ r n | ≤ Cε.
We thus have proved that C = E ∩ Ω (up to a negligible set). Let us finally denote by d the distance to R 2 \ E and set Σ := ∪ L l=1 d -1 ({lr}) where L is the integer part of r -1 max d. It is then not difficult to check that H 1 (Σ) < +∞ and Σ r = E because E satisfies the uniform interior ball property of radius r so that C = C Σ,r , which ends the proof. Remark 4.2. We have used the assumption that d = 2 only in the last step that is to prove that C = C Σ,r for some one-dimensional Σ. In higher dimensions, the same proof works if one requires H d-1 (Σ) < +∞ (however we believe the result remains true for one-dimensional sets in any dimension).

Remark 4.3. If the admissible sets Σ are supposed connected (in this case we call them networks), or with an a priori bounded number of connected components, then the penalization term |C Σ,r | can be replaced by the one-dimensional Hausdorff measure H 1 (Σ). In fact, for such sets we have

|C Σ,r | ≤ M 1 + H 1 (Σ)
where the constant M depends on the dimension d, on r, and on the number of connected components of Σ. Therefore the argument of Proposition 4.1 applies, providing the existence of an optimal solution.

We deal now with the case when the low-congestion region is a one-dimensional set Σ. We assume Σ connected (or with an a priori bounded number of connected components) and we take m(Σ) proportional to the one-dimensional Hausdorff measure H 1 (Σ). The integral on the low-congestion region has to be modified accordingly and we have to consider the problem formally written as

min σ,Σ Σ H 1 (σ) dH 1 + Ω H 2 (σ) dx + kH 1 (Σ) : σ ∈ Γ f (4.2)
with k > 0. Notice that, in view of the superlinearity assumption on the congestion functions H 1 and H 2 , the admissible fluxes u have to be assumed absolutely continuous measures with respect to L d ⌊Ω + H 1 ⌊Σ. Subsequently, the integral terms in the cost expression have to be intended as:

Σ H 1 dσ dH 1 dH 1 + Ω H 2 dσ dL d dx.
By an abuse of notation, when no confusion may arise, we continue to write the terms above as

Σ H 1 (σ) dH 1 + Ω H 2 (σ) dx.
In general, the optimization problem (4.2) does not admit a solution Σ opt , because the limits of minimizing sequences Σ n may develop multiplicities, providing as an optimum a relaxed solution made by a one-dimensional set Σ opt and function a ∈ L 1 (Σ opt ) with a(x) ≥ 1. The relaxed version of problem (4.2), taking into account these multiplicities, becomes

min σ,Σ,a Σ H 1 (σ/a)a dH 1 + Ω H 2 (σ) dx + k Σ a dH 1 : σ ∈ Γ f . (4.3)
The optimization with respect to a is easy: consider for simplicity the case

H 1 (σ) = α|σ| p with α > 0, p > 1;
then we have min

a≥1 ka + α |σ| p a p-1 = H(σ) =      α|σ| p + k if |σ| p ≤ k α(p -1) |σ|α 1/p p k p -1 1-1/p if |σ| p ≥ k α(p -1)
.

Therefore the relaxed problem (4.3) can be rewritten as

min σ,Σ Σ H(σ) dH 1 + Ω H 2 (σ) dx : σ ∈ Γ f
and the multiplicity density a(x) on Σ (that can be interpreted as the width of the road Σ at the point x) is given by

a(x) = 1 ∨ |σ(x)| α(p -1) k 1/p .

Numerical simulations

Here we wish to give a numerical example which clarifies and confirms what we expected from the analysis done in Section 3. In our examples, we mainly focus on the problem in the form (3.5):

min

Ω H * (∇u) dx - Ω f u dx .
The numerical simulation is based on a very simple situation that however seems quite reasonable. The two congestion function considered are both quadratic but with a different coefficient, say H 1 (σ) = a|σ| 2 and H 2 (σ) = b|σ| 2 with a < b. Then, in this case, the function H * involved in (3.5) is easy to compute:

H * (ξ) = ξ 2 4b ∨ ξ 2 4a -k
Before we start illustrating the numerical result, it is useful to do some considerations that justify the choice of some parameters in the following. The dual variable u has to be thought as a price system for a company handling the transport in a congested situation. An optimizer u then gives the price system which maximizes the profit of the company. When you take into account a congested transport between sources (here called f + and f -), the total mass plays an important role: as observed in [START_REF] Brasco | On certain anisotropic elliptic equation arising in congested optimal transport: local gradient bounds[END_REF], in the case of a small mass, hence of a large Lagrange multiplier k, the congestion effects are negligible, so one can expect in this case a distribution of the low-congestion region around the source distribution. On the contrary, for a large mass, hence for a small Lagrange multiplier k, we may expect a distribution of the low-congestion region also between the sources f + and f -.

In the following examples, we will consider as sources f + and f -two Gaussian distributions with variance λ, centered at two points x 0 and x 1

f + (x) = 1 √ 2πλ e -|x-x 0 | 2 /(2λ) , f -(x) = 1 √ 2πλ e -|x-x 1 | 2 /(2λ) .
Of course, a large value of λ means less concentration (and, on the contrary, a small λ captures more concentration); analogously, a large value of the penalization parameter k corresponds to a small quantity of available resources. Ending this consideration on parameters involved, we note that the traffic congestion parameters a, b and the "construction cost" parameter k are linked: we will change value of k according to a suitable choice of ratio a b , for fixed λ. Now, concerning the choice of the coefficients a, b we take a = 1 and b = 4, which means that the velocity in the low-congestion region is, at equal traffic density, four times the one in the region with normal congestion.

Using the equivalent dual formulation (3.5) of problem (3.2), we find numerically the solution u, hence the flux σ and the optimal density θ. Now, using the dual formulation of the problem, we find numerically the solution u of (3.5) and we obtain the flux σ as explained in Section 3. The numerical procedure to find u uses a Quasi-Newton method that updates an approximation of the Hessian matrix at each iteration (see [START_REF] Kelley | Iterative methods for optimization[END_REF] and reference therein). First we generate a finite element space with respect to a square grid. Then we implement the BFGS method, using a routine included in the packages of software FreeFem3D (available at http://www.freefem.org/ff3d ) that has the follow structure: BFGS(J,dJ,u,eps=1.e-6,nbiter=20)

The routine above means: find the optimal "u" for the functional J. The necessary parameters are the functional J, the gradient dJ and the u variable. Other parameters are optional and with clear means.

Example 5.1. The common setting of the simulation is a transportation domain Ω = [0, 1] 2 with a 30× 30 grid; we consider as initial and final distribution of resources two Gaussian approximations (with common variance λ) of Dirac delta function f -and f + respectively centered at x 0 = (0.3, 0.3) and x 1 = (0.7, 0.7). In the examples below we take different values of the parameters k and λ according to the considerations above, to show how the optimal distributions of the low-congestion regions may vary. Using the same notation as in Section 3, there are black and white region (respectively θ = 1 and θ = 0), passing through grey levels for the intermediate congestion.

In Figure 1 we take the variance parameter λ = 0.02, which provides the initial and final mass distributions not too concentrated, as depicted in Figure 1 (a). In Figure 1 (b) and (c) we take the penalization parameter k = 0.06 and k = 0.4 respectively; we see that in these cases, due to the low concentration of the initial and final mass distributions, the optimal density θ is higher in the region between x 0 and x 1 .

In Figure 2 we take the variance parameter λ = 0.001, which provides the initial and final mass distributions rather concentrated, as depicted in Figure 2 (a). In Figure 2 (b) and (c) we take the penalization parameter k = 0.01 and k = 0.05 respectively; we see that in these cases, due to the high concentration of the initial and final mass distributions, the optimal density θ is higher in the region around x 0 and x 1 . Proof. The inequality above can be deduced from the results in the appendix of [START_REF] Buttazzo | Asymptotic optimal location of facilities in a competition between population and industries[END_REF]; the present proof was obtained during a discussion with Giovanni Alberti, that we thank for his help. Since the set E r only depends on the closure of E, we may assume that E is closed; moreover, approximating E by smooth sets (for instance by the sets E s with s → 0), we may also assume that E is smooth.

4 Figure 1 (Figure 2 Appendix

 412 Figure 1

  Consider now the function f (r) = d|E r | -r Per(E r ); proving (A.1) amounts to show that f (r) ≥ 0 for every r > 0. Since E is assumed smooth, we have limr→0 |E r | = |E|, lim r→0 Per(E r ) = Per(E), so that lim r→0 f (r) = d|E| ≥ 0.By the coarea formula we have for all r < s|E s | -|E r | = Es\Er |∇ dist(x, E)| dx = s r Per(E t ) dt so that d dr |E r | = Per(E r ).Denoting by h(x) the mean curvature of ∂E r at x, and taking into account the definition of E r , we have h(x) ≤ (d -1)/r, so thatd dr Per(E r ) = ∂Er h(x) dH d-1 ≤ d -1 r Per(E r ).Therefore,f ′ (r) = d d dr |E r | -Per(E r ) -r d drPer(E r ) ≥ 0, which implies that f (r) ≥ 0 for every r > 0.

Acknowledgements.The work of the first and third authors is part of the project 2010A2TFX2 "Calcolo delle Variazioni" funded by the Italian Ministry of Research and University. The second author gratefully acknowledges the support of INRIA and the ANR through the Projects ISOTACE (ANR-12-MONU-0013) and OPTIFORM (ANR-12-BS01-0007).


G Carlier
email: carlier@ceremade.dauphine.fr

A Oberman
email: adam.oberman@mcgill.ca

E Oudet

Numerical methods for matching for teams and Wasserstein barycenters

Keywords: matching for teams, Wasserstein barycenters, duality, linear programming, numerical methods for nonsmooth convex minimization

Equilibrium multi-population matching (matching for teams) is a problem from mathematical economics which is related to multi-marginal optimal transport. A special but important case is the Wasserstein barycenter problem, which has applications in image processing and statistics. Two algorithms are presented: a linear programming algorithm and an efficient nonsmooth optimization algorithm, which applies in the case of the Wasserstein barycenters. The measures are approximated by discrete measures: convergence of the approximation is proved. Numerical results are presented which illustrate the efficiency of the algorithms.

Introduction

Optimal transport theory has received a lot of attention in the last decades and is now recognized as a powerful tool in PDEs, geometry, and functional inequalities (for which we refer to the monographs of Villani [START_REF] Villani | Topics in optimal transportation[END_REF]- [START_REF] Villani | Optimal transport: old and new[END_REF]). Given two Borel probability measures µ 1 , µ 2 , on metric spaces X 1 and X 2 , respectively, and a cost function c ∈ C(X 1 × X 2 , R), the Monge-Kantorovich optimal transport problem consists in finding the cheapest way to transport µ 1 to µ 2 for the cost c:

W c (µ 1 , µ 2 ) := inf γ∈Π(µ1,µ2) X×X c(x 1 , x 2 )γ(dx 1 , dx 2 ) (MK)
where Π(µ 1 , µ 2 ) denotes the set of transport plans between µ 1 and µ 2 , i.e. the set of probability measures on X 1 × X 2 having µ 1 , µ 2 , respectively as marginals.

Since this problem is of linear programming type, under very mild assumptions (e.g. when X 1 and X 2 are compact), the least transport cost W c (µ 1 , µ 2 ) admits a dual expression given by the Kantorovich duality formula:

W c (µ 1 , µ 2 ) := sup ϕ∈C(X1) X1 ϕ(x 1 )µ 1 (dx 1 ) + X2 ϕ c (x 2 )µ 2 (dx 2 ) (1.1)
where ϕ c denotes the c-transform of ϕ:

ϕ c (x 2 ) := min x1∈X1 {c(x 1 , x 2 ) -ϕ(x 1 )}.
A particularly important example is the quadratic case where X 1 = X 2 = R d , µ 1 and µ 2 have finite second moments, and c(x 1 , x 2 ) = |x 1x 2 | 2 . This case was first solved by Brenier [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF], who proved that whenever µ 1 is absolutely continuous, there is a unique optimal transport plan that is given by the gradient of a convex potential. This important result relates optimal transport to Monge-Ampère equations. We refer to [BFO14] and the references therein for numerical methods for optimal transport based on the Monge-Ampère equation. More generally, costs given by distances or convex power of distances are important because they lead to the so-called Wasserstein distances. More precisely, when X 1 = X 2 (a metric space with distance d) and c(x 1 , x 2 ) = d(x 1 , x 2 ) p for some p ≥ 1, the value W c (µ 1 , µ 2 ) in (MK) is the p-power of the so-called p-Wasserstein distance W p (µ 1 , µ 2 ) between µ 1 and µ 2 : W p (µ 1 , µ 2 ) := inf γ∈Π(µ1,µ2) X×X d(x 1 , x 2 ) p γ(dx 1 , dx 2 )

1/p .
In the present article, we are interested in solving numerically the following variant of the optimal transport problem which allows for more than two marginals. Given (compact metric, say) spaces X 1 , . . . , X I , equipped with Borel probability measures (µ 1 , . . . µ I ) ∈ P(X 1 ) × . . . × P(X I ), a (compact metric) space Z, and cost functions c i ∈ C(X i ×Z, R), we look for a probability measure ν on Z solving:

inf ν∈P(Z) J(ν) := I i=1 W ci (µ i , ν).
(1.2)

This problem was introduced in Carlier and Ekeland [START_REF] Carlier | Matching for teams[END_REF] in the framework of multi-population matching equlibrium; we will shortly recall in section 2 the economic interpretation of (1.2). Problem (1.2) is also a special case of multimarginal optimal transport (a variant of (MK) where more than two marginals are prescribed). Multi-marginal optimal transport is currently an active research field: compared to (two marginals) optimal transport, there are fewer theoretical results, and the complexity of general multi-marginal optimal transport problems typically increases exponentially in the number of marginals. Regarding the rapidly developing theory of multi-marginal optimal transport, we refer the reader to the recent papers by Pass [START_REF] Pass | Multi-marginal optimal transport and multi-agent matching problems: uniqueness and structure of solutions[END_REF], [START_REF] Pass | On the local structure of optimal measures in the multi-marginal optimal transportation problem[END_REF], by Ghoussoub and coauthors [START_REF] Ghoussoub | A self-dual polar factorization for vector fields[END_REF], [START_REF] Maurey | Remarks on multimarginals symmetric monge-kantorovich problems[END_REF] and the references therein for costs with special symmetry properties, motivated in particular by challenging computational issues in density functional theory in quantum physics.

We now discuss a special, but important case of (1.2) which has a clear geometric interpretation. Let all the X i 's and Z coincide with R d , the measures µ i have finite second moments, and the costs be quadratic (i.e. c i (x i , z) := λ i |x i -z| 2 for some weights λ i > 0, summing to 1 without loss of generality). In this case (1.2) takes the form:

inf ν∈P(R d ) J(ν) := I i=1 λ i W 2 2 (µ i , ν) (1.3)
where W 2 denotes the 2-Wasserstein distance. In analogy with the Euclidean case, a solution to (1.3) will be called a Wasserstein barycenter of the measures µ i with weights λ i . Properties of Wasserstein barycenters were studied by Agueh and Carlier [AC11]. Wasserstein barycenters interpolate between the measures µ i ; the idea of interpolating between points of a metric space by minimizing some weighted sum of squared distances goes back to the notion of Fréchet mean. The case I = 2 is well-known. Letting the weights (λ, 1λ) vary, one obtains the classical notion of McCann's interpolation [START_REF] Mccann | A convexity principle for interacting gases[END_REF] between two probability measures. This interpolating curve is also a geodesic for W 2 , and in their seminal paper [BB00] on the dynamic formulation of optimal transport, Benamou and Brenier gave a numerical scheme to compute this geodesic. Finding barycenters between more than three measures is more complicated (barycenters are not associative as soon as d ≥ 2). From a Partial Differential Equations problem, this problem requires to solve a system of Monge-Ampère equations, see (4.9)-(4.8) below. Interestingly, the Wasserstein barycenter problem recently found natural applications in image processing, see Peyré et al. [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] and statistics, see Bigot and Klein [START_REF] Bigot | Consistent estimation of a population barycenter in the wasserstein space[END_REF]. Of course, there are lots of variants of the interpolating scheme given by the W 2 -barycenter problem (1.3) and in particular one can replace W 2 by W p for some p ≥ 1 or even mix different powers of the distance. Slightly abusing the terminology, we will sometimes refer to barycenters even for these variants and even for the general form of the problem (1.2).

In the discrete setting, the transportation problem is classical. In fact, this problem motivated the historical development of optimization, by Kantorovich in 1939, working on Soviet railway transportation, and in the 1940's by Hitchcock [Sch03, Ch 21]. The "assignment problem" arises in case of integer values weights, it is a standard combinatorial optimization problem which can be solved by the Hungarian algorithm) [START_REF] Schrijver | Combinatorial optimization: polyhedra and efficiency[END_REF]Ch 21]. More generally, the "transportation problem" is a linear programming problem which arises when the weights are real-valued, it can be solved by the Hitchcock algorithm, [Sch03, Ch 21], or by modern commerical general linear programming software. Returning to the problem with continuous measures, it is natural to approxiate the measure by weighted sums of delta measures. In theory, the resulting problem can be solved using linear programming. However, the number of variables in the linear programming problem is quadratic in the number of variables used to represent the measures. In the discrete setting, current optimization algorithms are limited to approximately several thousands of variables for each of the measures. This problem size corresponds to a fairly coarse approximation of a two dimensional continuous measure. In special cases, or using specific approximations, improvements are available, see [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF] for quadratic costs, and for more references. For example, if each measure is represented by, for example, 40 2 = 1600 variables, the linear program has 40 4 = 2 560 000 which is near the limitations of linear programming algorithms (we performed experiments using CVX [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF] and calling several academic and commercial optimization packages). Enlarging the resolution of the measures quickly overwhelms the capabilities of the algorithms.

The problem (1.2) is even more challenging, since it involves multiple marginals and an additional unknown measure. A reasonable goal is to resolve each of the measures in two dimensions, with (on the order of, say) 1600 variables. Since we wish to allow for measures which are separated, this may require a full grid of size up to 200 2 . Resolving the barycenter measure on the full grid generally leads to an intractable problem. Our main contributions regarding numerical schemes for the general problem (1.2) or the particular case of Wasserstein barycenters (1.3) are as follows:

• We give a simple linear programming reformulation of (1.2) in subsection 2.3 whose size is proportional to the number of marginals. Together with a localization result that bounds the support of the unknown barycenter in subsection 2.2, one then obtains a tractable problem. We discretize the problem to arrive at a finite dimensional linear programming problem in subsection 2.4. We prove convergence, in the sense of weak convergence of measures, in subsection 2.5.

• Numerical results are presented in section 3. These illustrate the validity of this linear programming approach. Barycenter problems with different costs are solved, as well as a matching for teams problem.

• The second algorithm which is specialized to the case of Wasserstein barycenter measures (1.3), is described and implemented in section 4. This problem uses the dual formulation of the problem explained in section 4, and special features of the quadratic cost. The efficient nonsmooth optimization algorithm is described in subsection 4.3. Large size computational examples are presented (on grids of size 200 2 , and for measures resolved with 1500 points). The examples include barycenter measures using up to five measures, and a an example in texture synthesis in subsection 4.4.

Matching for teams and approximation

In this section, following [START_REF] Carlier | Matching for teams[END_REF], we first derive the generalized barycenter problem (1.2) as an equivalent reformulation of an equilibrium problem for multipopulation matching arising in economics. Next, we study localization of the barycenter measure. Then, we present an infinite dimensional linear programming reformulation of (2.5). This is followed by a discretization of the measures, which results in a finite dimensional linear programming problem that is tractable for moderate problem sizes. Finally, we address stability issues (in the sense of weak convergence of measures) when one approximates the measures µ i by some (discrete) measures.

Variational characterization of matching equilibria

The model of Carlier and Ekeland [START_REF] Carlier | Matching for teams[END_REF] deals with the equilibrium of a market for a quality good (e.g. house, school, hospitals, . . . ). Producing the good requires assembling a team consisting of a buyer and a set of producers. For instance, in the case of houses, the producers could be a plumber, an electrician and a mason. The quality good has a range of feasible qualities (location, surface, number of rooms, facilities etc...), denoted by Z which we assume to be a compact metric space. Each of the different populations (buyers, plumbers, electricians, masons...) is indexed by i ∈ {1, . . . , I}. The agents in each population are hetererogeneous, characterized by a certain type which affects their (quality dependent) cost function. For example, some masons are used to working with lower quality bricks, while other work with luxury stone, some electricians leave quite far from the location of the house they work on, consumers differ in their tastes... To be precise, for each population i, we are given a compact metric space of types, X i , and a continuous cost function c i ∈ C(X i × Z, R); with the interpretation that c i (x i , z) is the cost for an agent of population i with type x i to work in a team that produces good z. The distribution of type x i in population i is known and given by some Borel probability measure µ i ∈ P(X i ).

The goal is to find an equilibrium production line ν ∈ P(Z) (together with a price system) which clears both the quality good and the labor market. The equilibrium is described below, and as we shall see, it corresponds to the solution of the (generalized) barycenter measure problem (2.5). In this setting, one looks in particular for an equilibrium system of monetary transfers (paid by the buyer to the producers). A system of transfers is a collection of continuous functions ϕ 1 , . . . ϕ I : Z → R where ϕ i (z) is the amount paid to i by the other members of the team for producing z. An obvious equilibrium requirement is that teams are self-financed i.e. (2.1)

Given transfers, ϕ 1 , . . . ϕ I , an agent from population i with type x i ∈ X i , gets a net minimal cost given by the so-called c i -transform of ϕ i :

ϕ ci i (x i ) := min z∈Z {c i (x i , z) -ϕ i (z)}. (2.2) By construction, ϕ ci i (x i ) + ϕ i (z) ≤ c i (x i , z)
, and since agents are rational, they choose cost minimizing qualities, i.e. a z ∈ Z such that

ϕ ci i (x i ) + ϕ i (z) = c i (x i , z). (2.3)
The final unknown is a collection of plans, γ i ∈ P(X i × Z), such that γ i (A i × A) represents the probability that an agent in population i has a type in A i , and belongs to a team that produces a quality in A. At equilibrium, the first marginal of γ i should be µ i (this is equilibrium on the i-th labor market) and the second marginal of γ i should not depend on i (this is equilibrium on the quality good market), this common marginal represents the equilibrium quality line. An equilibrium can then be formally defined. It consists of a transfer system (ϕ 1 , . . . ϕ I ) ∈ C(Z, R) I , probability measures γ i ∈ P(X i × Z), and a probability measure ν ∈ P(Z), such that

• teams are self-financed i.e. (2.1) holds,

• γ i ∈ Π(µ i , ν) for i = 1, . . . , I (equilibrium on the labor markets and on the good market),

• (2.3) holds on the support of γ i for i = 1, . . . , I, (i.e. agents choose cost minimizing qualities).

If an equilibrium quality line, ν, was known, then clearly the last two conditions above would imply that the plan γ i should be optimal for the Monge-Kantorovich problem:

W ci (µ i , ν) := inf γ∈Π(µi,ν) Xi×Z c i (x i , z)γ(dx i , dz).
(2.4)

In fact, it was proved in [START_REF] Carlier | Matching for teams[END_REF] that there is a purely variational characterization of equilibria, which is tightly related to the following convex problem inf ν∈P(Z)

J(ν) := I i=1 W ci (µ i , ν) (2.5)
and its dual (concave maximization) formulation (see [START_REF] Carlier | Matching for teams[END_REF] or section 4 for details on this duality) sup

I i=1 Xi ϕ ci i (x i )µ i (dx i ) : I i=1 ϕ i = 0 . (2.6) Theorem 2.1. (ϕ i , γ i , ν
) is an equilibrium if and only if:

• ν solves (2.5),
• the transfers (ϕ 1 , . . . ϕ I ) solve (2.6),

• for i = 1, . . . , I, γ i solves the Monge-Kantorovich problem W ci (µ i , ν).

Localization

As noted in [START_REF] Carlier | Matching for teams[END_REF], the minimization problem (2.5), which characterizes equilibrium quality lines, can be reformulated as an optimal transport problem with multi-marginal constraints, as follows. First define the cost

c(x) := min z∈Z I i=1 c i (x i , z), (2.7)
where x = (x 1 , . . . , x I ). Let T (x) ∈ Z be a measurable selection of the solution of the above minimization, meaning that T (x) ∈ Z satisfies

I i=1 c i (x i , T (x)) = c(x).
Then consider the multi-marginal problem

inf γ∈Π(µ1,...,µ I ) X1×...×X I c(x)γ(dx), (2.8)
where Π(µ 1 , . . . , µ I ) denotes the set of probability measures on X 1 × . . . × X I having (µ 1 , . . . , µ I ) as marginals. It is not difficult to see that if γ solves (2.8) then ν := T # γ solves (2.5) (where as usual T # γ denotes the push forward of γ through T , i.e. T # γ(B) := γ(T -1 (B)) for every Borel B).

Conversely, one can relate the minimizers of (2.5) to those of (2.8). Indeed, let ν solve (2.5) and let γ i ∈ Π(µ i , ν) be an optimal plan for W ci (µ i , ν). Disintegrating γ i with respect to ν i.e. writing γ i = γ z i ⊗ ν and defining γ ∈ P(X 1 , × . . . × X I × Z) by:

γ := ⊗ I i=1 γ z i ⊗ ν
and η as the marginal of γ on the variables (x 1 , . . . , x I ), one easily checks that

• η ∈ Π(µ 1 , . . . , µ I ) solves (2.8),
• on the support of γ, spt(γ), one has

I i=1 c i (x i , z) = c(x),
• the previous relation, together with the fact that ν is the Z-marginal of γ and µ i its X i -marginal then imply a useful localization property: the support of the barycenter measure, spt(ν), is contained in the set of minimizers of the following problem.

min z∈Z I i=1 c i (x i , z) for some x i ∈ spt µ i , i = 1, . . . , I.
(2.9)

Since the support of ν is unknown, which causes difficulties in practice, the localization property (2.9) gives a practical method for bounding the unknown support of the barycenter measure. The condition above results in a reduction of the dimensionality of discretized problems since it gives an a priori information on the support of the unknown measure, at the expense of solving an optimization problem. However this optimization problem is decoupled on the domain Z: each point (or neighborhood) can be tested by looping through points (or small neighborhoods) in the domain Z and choices of points in the support sets spt(µ i ).

In the case where X i and Z coincide with some ball of R d , and the costs are powers of distance, c i (x i , z) = λ i |x i -z| p (with λ i > 0 and λ i = 1, say) for some p ≥ 1, one can easily derive an information on the unknown support. Indeed, using the optimality condition for the minimization problem (2.9), one deduces that spt(ν) is included in the convex hull of the supports of the µ i 's. If we particularize further to the the Wasserstein barycenter case, i.e. to the case p = 2 , the solution of (2.9) is explictly given by the barycenter z = I i=1 λ i x i so that the localization property (2.9) gives the following estimate on the barycenter measure ν:

spt(ν) ⊂ I i=1 λ i spt(µ i ).
(2.10)

Linear programming formulation

Multi-marginals optimal transport problems such as (2.8) are linear programs. For discrete marginals, such problems can in principle be solved exactly by the simplex method. In practice however, the number of variables explodes with the number of marginals, which makes the problem quickly untractable. We shall see below that one may take advantage of the fact that c is not any cost function but has the special structure (2.7). Interestingly, it was already proved by Pass [START_REF] Pass | Multi-marginal optimal transport and multi-agent matching problems: uniqueness and structure of solutions[END_REF] in the context of multi-marginal optimal transport that such costs are much more well-behaved than arbitrary costs of I variables.

To find a more tractable linear programming reformulation of the matching for teams problem, it is better to go back to the very definition of an equilibrium in terms of couplings and to reformulate problem (2.5) as inf (γ1,...,γ I )∈Π

I i=1 Xi×Z c i (x i , z)γ i (dx i , dz) (2.11)
where Π consists of all measures (γ 1 , . . . , γ I ) ∈ P(X 1 × Z) × . . . × P(X I × Z) such that

• the marginal of γ i on the

x i variable is µ i i.e. Xi×Z ψ(x i )γ(dx i , dz) = Xi ψ(x i )µ i (dx i ), ∀ψ ∈ C(X i ),
(2.12)

• the marginal of γ i on the z variable does not depend on i:

Z ϕ(z)γ 1 (dx 1 , dz) = . . . = X I ϕ(z)γ I (dx I , dz), ∀ϕ ∈ C(Z).
(2.13) Clearly, if the γ i 's solve (2.11) then their common marginal ν ∈ P(Z) solves (2.5) and the γ i 's are optimal for the optimal transport problem W ci (µ i , ν). In other words, the γ i 's are equilibrium couplings for the matching for teams problem.

The constraints above being linear, Π is a convex and weakly star compact subset of P(X 1 ×Z)×. . .×P(X I ×Z) so that (2.11) admits solutions. Moreover in the case of discrete µ i 's and ν supported by N points, the number of variables in the linear program (2.11) is linear (and not exponential as in the case of the multi-marginal optimal transport problem) in the number of marginals.

Discretization

The (a priori) infinite dimensional linear programming problem (2.11) of subsection 2.3 can be discretized as follows. Let {S i j } Ni j=1 be a partition of spt(µ i ) and let {S 0 k } N0 k=1 be a partition of Z (or better, of the support set estimated by the method of subsection 2.2). Approximate the measures by weighted sums of atoms

µ A i = j µ j i δ x i j , for i = 1, . . . , I, with µ j i = µ i (S i j ) ν A = k ν k δ z 0 k , with ν k = ν(S 0 k )
where x i j and z 0 k are representative points in the regions S i j , S 0 k , respectively. It is well-known, that µ A i converges weakly * to µ i as the diameter of the partition {S i j } Ni j=1 tends to 0. More precisely, denoting by W 1 the 1-Wasserstein distance (which metrizes the weak * topology on probability measures):

W 1 (µ A i , µ i ) ≤ max j=1,...,Ni diam(S i j ).
(2.14)

Inserting the approximation defined above into the linear programming problem (2.11, 2.12, 2.13) results in the following finite dimensional linear programming problem minimize

I i=1 j,k c i (x i j , z 0 k )γ i j,k subject to: k γ i j,k = µ j i ,
for all k = 1, . . . , N i , and for i = 1, . . . , I

j γ 1 j,k = . . . = j γ I j,k , for all k = 1, . . . , N 0 ., (2.15)
along with the non-negativity constraints γ i j,k ≥ 0. The (approximated) barycenter is then ν A = k ν k δ z 0 k where the weight ν k is given by the common value,

ν k = j γ i j,k ,
for any i = 1, . . . , I.

The linear programming problem above can be implemented in standard software packages. The size of the problem above is as follows. The number of variables is

N 0 × (N 1 + • • • + N I ) (or IN 2 if each N 0 = N 1 = • • • = N I = N ). The number of constraints is (N 0 +1)×(N 1 +• • •+N I )+IN 0 (that is I(N 2 +2N ) when N 0 = N 1 = • • • = N I = N ).
The size of this linear programming problem thus scales linearly with the number of marginals, for a given, fixed value of N (contrary to the multi-marginal formulation (2.8)).

Approximation and convergence

Since in practice, one considers approximation by discrete measures just as in subsection 2.4, we wish now to address the stability of the following convex problem when one replaces the measures µ i by some discrete approximation inf ν∈P(Z)

J(ν) := I i=1 W ci (µ i , ν). (2.16)
To do so, one has to control the dependence of W c (µ, ν) in its three arguments (c, µ, ν) ∈ C(X × Z) × P(X) × P(Z). We shall denote by d X and d Z the distances on X and Z, take now (c, µ, ν) ∈ C(X × Z) × P(X) × P(Z) and ( c, µ, ν) ∈ C(X × Z) × P(X) × P(Z) and let ω X and ω Z be respectively a modulus of continuity of c and c with respect to x uniform in z and a modulus of continuity of of c and c with respect to z uniform in x, that is

max(|c(x, z)-c(x ′ , z)|, | c(x, z)-c(x ′ , z)|) ≤ ω X (d X (x, x ′ )), ∀(x, x ′ , z) ∈ X×X×Z and max(|c(x, z)-c(x, z ′ )|, | c(x, z)-c(x, z ′ )|) ≤ ω Z (d Z (z, z ′ )), ∀(x, z, z ′ ) ∈ X ×Z ×Z.
Obviously, one has

|W c (µ, ν) -W c (µ, ν)| ≤ c -c ∞ .
(2.17)

Let ϕ ∈ C(Z) be a solution in the Kantorovich dual of W c (µ, ν), that is W c (µ, ν) = X ϕ c dµ + Z ϕdν
by the Kantorovich duality formula, we have

W c ( µ, ν) ≥ X ϕ c d µ + Z ϕdν.
Hence, for every θ ∈ Π(µ, µ), we have

W c (µ, ν) -W c ( µ, ν) ≤ X ϕ c d(µ -µ) = X×X (ϕ c (x) -ϕ c (x ′ ))θ(dx, dx ′ ).
we then observe that

ϕ c (x) -ϕ c (x ′ ) ≤ ω X (d X (x, x ′ )) so that W c (µ, ν) -W c ( µ, ν) ≤ W ω X (µ, µ) := inf θ∈Π(µ, µ) X×X ω X (d X (x, x ′ ))θ(dx, dx ′ ).
(2.18) Similarly

W c ( µ, ν) -W c ( µ, ν) ≤ W ω Z (ν, ν) := inf η∈Π(ν, ν) Z×Z ω Z (d Z (z, z ′ ))η(dz, dz ′ ).
(2.19) Putting everything together, we get

|W c (µ, ν) -W c ( µ, ν)| ≤ c -c ∞ + W ω X (µ, µ) + W ω Z (ν, ν).
(2.20)

We then observe that if µ n weakly * converges to µ then W ω X (µ, µ n ) → 0. Indeed, it is known to imply that the 1-Wasserstein distance (corresponding to W ω X for ω X (t) = t) between µ n and µ converges to 0, so that there is some θ n ∈ Π(µ, µ n ) which (up to a non relabeled subsequence) weakly * converges to some θ supported on the diagonal of X × X, hence

W ω X (µ, µ n ) ≤ X×X ω X (d(x, x ′ ))θ n (dx, dx ′ ) → 0.
Getting back to the approximation of (2.5), take sequences c n i ∈ C(X i × Z), µ n i ∈ P(X i ), and c i ∈ C(X i × Z), µ i ∈ P(X i ), such that

c n i -c i ∞ → 0, µ n i ⇀ * µ i (2.21)
and set:

J(ν) := I i=1 W ci (µ i , ν), J n (ν) := I i=1 W c n i (µ n i , ν), ∀ν ∈ P(Z). (2.22)
Denoting by ω i Xi and ω i Z common continuity modulus of the c n i (the first one in x i uniformly in z and the second in z, uniformly in x i just as above) we then have: Proposition 2.2. For every (ν, ν n ) ∈ P(Z) × P(Z) :

|J(ν) -J n (ν n )| ≤ I i=1 c i -c n i ∞ + I i=1 [W ω i X i (µ n i , µ i ) + W ω i Z (ν n , ν)] (2.23)
this implies in particular

• J n (ν n ) → J(ν) whenever ν n ⇀ * ν,
• a quantitative estimate for the stability of values:

| inf P(Z) J -inf P(Z) J n | ≤ I i=1 c i -c n i ∞ + I i=1 W ω i X i (µ n i , µ i ) (2.24)
• if ν n minimizes J n then, up to a subsequence, it weakly * converges to a minimizer of J.

Proof. The statements directly follow from estimate (2.20) and the already observed fact that the right-hand side of (2.23) converges to 0 as soon as ν n ⇀ * ν.

In the case where the cost functions c i = c n i are Lipschitz and the approximated measures µ n i satisfy (for the usual 1-Wasserstein distance) W 1 (µ n i , µ i ) ≤

C n , (2.23) above just takes the form

|J(ν) -J n (ν n )| ≤ I i=1 Lip(c i ) C n + W 1 (ν n , ν) .
3 Numerical simulations: Linear Programming

In this section, we present various numerical simulations using the Linear Programming approach of sections 2.3 and 2.4. The localization method of 2.2 is used to approximate the support of the barycenter measure. An alternate approach to approximating the support of the barycenter measure, which can be combined with localization is a a two stage solution approach: the first stage, using a coarse grid, gives an approximation of the support of the barycenter, the second stage gives a more accurate representation of the barycenter using information on the support obtained in the first stage. All computations in this section were performed in MATLAB on a Mid 2011 MacBook Air laptop. To solve (2.15) we use the software package CVX [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF] [GB08] which is callable from MATLAB. The CVX language allows for a very concise description of the convex optimization problem, and allows for the use of multiple solver libraries (e.g. MOSEK, Gurobi). The numerical solution obtained is the correct up to tolerances near numerical precision.

Geodesic paths between measures in the plane

We considered two measures in the plane, and by varying the weights in the quadratic cost function, we computed three points on the geodesic path (or McCann interpolant) represented in figure 1. The computational time was less Right: solution on the implied support, but more resolved. The size of the problems is the same, but the resolution is increased by a factor of three.

Comparing different cost functions

For the next set of examples, we took two uniform measures, the first corresponding to a vertically oriented rectangle, and the second corresponding to an horizonal rectangle. These measures are shown in Figure 3. First, we compared the convergence of the solutions for different grid sizes in Figure 4. Notice that the general support of the computed measures seems stable, but there are oscillations in the density, for different resolutions.

Next, we computed the barycenter with various power cost functions C(x, y) = |x-y| p , for p = 1, 2, 3, 4. The solutions we computed use grids of size 50 2 . Computational time was close to two minutes. A second run using grid size of 100 2 and a localization of the support, took 30 to 45 minutes. See figure 5. The densities are plotted using a grayscale which corresponds to the relative values, however the grayscale is different for each figure. We also include another view of the density for p = 1 in Figure 6.

The solutions have a complicated geometry. For the cost with p = 1, the support of the barycenter is the entire convex hull of the supports of the measures, although the density is highly concentrated at the intersection of the measures. The density ranges from about 0.01 at the edges to 0.12 in the center.

For the case p = 2, the density is supported on a square, but wider than the width of the rectangle. The density has some oscillations, but it strictly positive (taking values in the range [.005,.004]). For the case p = 3, the density is supported on a small octagonal shape, with zero density in the middle, and with larger oscillations. For the case p = 4, a much larger octagonal shape appears with a large zero density hole in the middle. The supports of the barycenter measures are close to the ones estimated by localization. Finally, we computed the barycenters using all four measures from figure 3, see figure 7.

Matching for teams

We considered the matching for teams problem and used measures and costs as follows, also see Figure 8. Set the quality domain Z = [0, 1] 2 and write z = (z 1 , z 2 ) for points in Z. Set M 0 , M 1 , M 2 to be measures which have constant density on their support, and let their supports be [1, 2] 2 , [1.25, 1.75] × [1, 2],

c 0 (x, z) = -5.5(x 1 z 1 + x 2 z 2 ) c 1 (x, z) = c 2 (x, z) = (x 1 + z 1 ) 2 + (x 2 + z 2 ) 2 .
The solution concentrates mass at the boundary, and especially at the corners of the domain.

Dual formulation

Duality and optimality conditions

Let us now explain why the variational problem (2.6) can be naturally be seen as a dual formulation of (2.11) (see [START_REF] Carlier | Matching for teams[END_REF] for more details on this duality). To that end, let us observe that (γ 1 , . . . , γ I ) ∈ Π if and only if (2.12) holds for every i (these are the fixed µ i marginals constraints) and

X I ×Z I i=1 ϕ i (z)γ i (dx i , dz) = 0, as soon as I i=1 ϕ i (z)dz = 0, ∀z ∈ Z. (4.1)
Indeed, clearly if the γ i 's have the same marginal on Z then (4.1) holds. Conversely assume (4.1), let i = j and ϕ ∈ C(Z) then applying (4.1) to the potentials ϕ i = ϕ, ϕ j = -ϕ and ϕ k = 0 for k / ∈ {i, j} we see that Xi×Z ϕ(z)γ i (dx i , dz) = Xj ×Z ϕ(z)γ j (dx j , dz). This proves that (4.1) characterizes the fact that the γ i 's share the same marginal on Z. This enables us to rewrite (2.11) in inf-sup form:

inf γi≥0 sup L((γ i ) i , (ψ i ) i , (ϕ i ) i ) : ψ i ∈ C(X i ), ϕ i ∈ C(Z) : I i=1 ϕ i = 0 (4.2)
where the Lagrangian L is given by

L((γ i ) i , (ψ i ) i , (ϕ i ) i ) := I i=1 Xi×Z (c i (x i , z) -ψ i (x i ) -ϕ i (z))γ i (dx i , dz) + I i=1 Xi ψ i (x i )µ i (dx i ).
To obtain the desired dual formulation, we formally switch the inf and the sup (again, we refer to [START_REF] Carlier | Matching for teams[END_REF] for a rigorous derivation):

sup

(ψi,ϕi), ϕi=0 inf γi≥0 L((γ i ) i , (ψ i ) i , (ϕ i ) i ).
We next observe that

inf γi≥0 L((γ i ) i , (ψ i ) i , (ϕ i ) i ) = I i=1 Xi ψ i (x i )µ i (dx i ) + I i=1 inf γi≥0 Xi×Z (c i (x i , z) -ψ i (x i ) -ϕ i (z))γ i (dx i , dz)
and the latter infimum is 0 when

c i (x i , z) ≥ ψ i (x i ) + ϕ i (z), ∀(x i , z) ∈ X i × Z (4.3)
and -∞ otherwise. The dual of (2.11) therefore consists in maximizing

I i=1 Xi ψ i (x i )µ i (dx i )
subject to the constraints (4.3) and I i=1 ϕ i = 0. For fixed ϕ i , the maximal ψ i that satisfies (4.3) being ψ i = ϕ ci i , we see that that the dual can be equivalently formulated as sup

I i=1 Xi ϕ ci i (x i )µ i (dx i ) : I i=1 ϕ i = 0 (4.4)
which is exactly (2.6). For the existence of solutions and the equality between the infimum in (2.11) and the supremum in (4.4) (which is obtained by a slightly different argument), we again refer to [START_REF] Carlier | Matching for teams[END_REF]. Now the optimality conditions for (2.11) and (4.4) are summarized by the equivalence between the following assertions:

• (γ i ) i ∈ Π solves (2.11) and (ϕ i ) i such that I i=1 ϕ i = 0 solves (4.4), • ((γ i ) i , (ϕ ci i ) i , (ϕ i ) i ) is a saddle point of L,
• for every i, one has

ϕ ci i (x i ) + ϕ i (z) = c i (x i , z) (4.5)
γ i -almost everywhere on X i × Z or, equivalently, by continuity, on the support of γ i .

The case of Wasserstein barycenters

From now on, we restrict ourselves to the quadratic case where all the X i 's and Z are some ball B (say) of R d and the costs c i are quadratic:

c i (x i , z) := λ i 2 |x i -z| 2
where the λ i 's are positive coefficients which we normalize in such a way that they sum to 1. In this case, (2.5) corresponds to inf

ν∈P(B) I i=1 λ i W 2 2 (µ i , ν) (4.6)
where W 2 2 stands for the squared 2-Wasserstein distance. This problem has been studied in details in [AC11] where uniqueness (under the assumption that one of the measures does not give mass to small sets), characterization, L p or L ∞ regularity results are established for Wasserstein barycenters as well as a close connection with the quadratic multimarginal optimal transport problem of Gangbo and Świȩch [START_REF] Gangbo | Optimal maps for the multidimensional monge-kantorovich problem[END_REF]. Since Wasserstein barycenters constitute a natural way to interpolate between an arbitrary number of measures, they therefore also find applications in image processing [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] and statistics [START_REF] Bigot | Consistent estimation of a population barycenter in the wasserstein space[END_REF].

Let us now informally give the optimality conditions for the barycenter using once again the dual formulation (4.4) (see [AC11] for details). In the present quadratic setting, the formula for the c i -transform takes the form

ϕ ci i (x i ) = inf z λ i 2 |x i -z| 2 -ϕ i (z) ,
which, defining

u i (x i ) := 1 2 |x i | 2 - ϕ ci i (x i ) λ i
can conveniently be rewritten as

u i = 1 2 |.| 2 - ϕ i λ i *
(where * denotes the standard Legendre transform). In particular, u i is convex (hence differentiable outside of a small set). Moreover the optimal coupling γ i is concentrated on the set where equality (4.5) holds which setting

v i = u * i is equivalent to the relation u i (x i ) + v i (z) = x i • z i.e. z = ∇u i (x i ) (provided u i is differentiable at x i
which is the case µ i a.e. as soon as µ i does not charge small sets...). This implies that the barycenter which is also the marginal ν that is common to all the γ i 's is actually given by ν = ∇u i# µ i for every i and ∇u i is the optimal transport between µ i and ν for W 2 2 . We then have

1 2 |.| 2 - ϕ i λ i ≥ 1 2 |.| 2 - ϕ i λ i * * = u * i = v i . (4.7)
But one can rewrite the fact that γ i -almost everywhere equality (4.5) holds as the fact that for γ i -almost every (x i , z), one has

1 2 |z| 2 - ϕ i (z) λ i = x i • z -u i (x i ) ≤ u * i (z) = v i (z).
With (4.7), we then deduce that for ν-a.e. z, one has

1 2 |z| 2 - ϕ i (z) λ i = v i (z).
Recalling that the ϕ i have to sum to 0, we deduce that

I i=1 λ i v i (z) = |z| 2 2 (4.8)
on the support of ν. The optimality conditions for the barycenter ν = ∇u i# µ i = ∇v i * # µ i therefore, at least formally take the form of the system of Monge-Ampère equations

ν = µ i (∇v i ) det(D 2 v i ), i = 1, . . . , I (4.9)
which is supplemented with equation (4.8) on the support of ν. We shall see in the next paragraph how to compute numerically in an efficient way the potentials ϕ i .

An efficient algorithm for Wasserstein barycenters

Discretization of the dual problem. We assume in all this section that the sets X i 's and Z are subsets of R d for some d = 1, 2 or 3. As described in the previous section, the computation of one Wasserstein quadratic barycenter is equivalent in its dual form to the maximization of

I i=1 Xi ϕ ci i (x i )µ i (dx i ) (4.10)
where

ϕ ci i (x i ) = inf z λ i 2 |x i -z| 2 -ϕ i (z) ,
under the linear constraint I i=1 ϕ i (z) = 0 for all z ∈ Z. This formulation leads to the following natural discretization of Wasserstein quadratic barycenter computation. Suppose (y j i , ν j i ) j=1,...,Ni ⊂ X i × R + is a convergent quantization of the measures µ i . More explicitly, we assume that for all i = 1, . . . , I

lim Ni→∞ c(N i ) Ni j=1 ν j i δ y j i = µ i
in the sense of the weak convergence of measures. In order to only consider probability measures, we set c

(N i ) = ( Ni j=1 ν j i ) -1
. Additionally, we suppose that (z k ) is a dense countable family of points of Z. Based on (4.10) and previous quantizations, for a given

N k ∈ N, our discrete optimization problem of I × N k variables reads Φ((ϕ k 1 ), . . . , (ϕ k I )) = I i=1 c(N i ) Ni j=1 ν j i min k=1,...,N k λ i 2 |y j i -z k | 2 -ϕ k i (4.11)
under the N k pointwise linear constraints:

i ϕ k i = 0, ∀k = 1, . . . , N k . (4.12)
This optimization problem in its dual form can be seen as a large scale nonsmooth concave maximization problem. We discuss in the next paragraph alternatives that have been developed to solve numerically this type of problems.

Non smooth algorithms. Many different approaches have been introduced in the last decades to approximate optimal solution of non-smooth concave (or convex) problems, e.g. gradient sampling methods [START_REF] Burke | A robust gradient sampling algorithm for nonsmooth, nonconvex optimization[END_REF] and bundle methods [START_REF] Lemaréchal | International Institute for Applied Systems Analysis[END_REF]. These algorithms make use of a partial or complete description of superdifferentials in order to identify ascent directions (see next paragraph). Even though Proposition 4.1 describes explicitly the whole superdifferential, finding an effective ascent direction in practice is made difficult by the large dimension of some superdifferentials. Additionally, those approaches are essentially of order one and follow the singular parts of the graph of the cost function. These two facts could explain a slow rate of convergence when starting from an initial point far away from any optimal vector. One surprisingly efficient alternative for minimizing non-smooth functions is the use of quasi-Newton methods. It is known [START_REF] Powell | Some global convergence properties of a variable metric algorithm for minimization without exact line searches[END_REF] that if the maximized function, Φ, is twice continuously differentiable and the suplevel set Φ ≥ Φ(x 0 ) is bounded, then the sequence of function values generated by the BFGS method with inexact Armijo-Wolfe line search, starting from x 0 converge to the maximal value of Φ. More recently, it has been pointed out by different authors [HUL96, LO09, LO08] that variable metric algorithms may produce in some cases sequences which converge to an optimal point in the sense of Clarke. The mathematical analysis of this good behavior has just been initiated in recent papers of Overton [START_REF] Lewis | Nonsmooth optimization via BFGS[END_REF][START_REF] Lewis | Behavior of BFGS with an exact line search on nonsmooth examples[END_REF]. This efficiency could be explained heuristically by the fact that the approximated inverse of the Hessian matrix has a spectral decomposition in two subspaces which describe the two different local behaviors of the cost Φ: a subspace associated to the regular directions of the cost function Φ, and the subspace of eigenvectors whose eigenvalues are small in absolute value, which correspond to the singular directions of Φ.

It has been observed in simple situations that L-BFGS algorithms are sometimes able to converge to an optimal point. In more standard examples, where concentration can occur for instance, L-BFGS approach fails to converge. This expected behavior for strongly non-smooth functions illustrates the need for using more specific non-smooth techniques close from the optimal point. The costly, but reliable, bundle type algorithms have demonstrated their efficiency in this context.

We will not give here a detailed study of quasi Newton methods applied to optimal transportation which would be out of the scope of this paper. We only point out that the L-BFGS (low memory version of Broyden-Fletcher-Goldfarb-Shanno algorithm) algorithm combined with a bundle approach gives a rather efficient way to solve this type of problem. We refer to [START_REF] Haarala | Globally convergent limited memory bundle method for large-scale nonsmooth optimization[END_REF] for a careful study and an efficient implementation of this kind of hybrid algorithm.

Gradient computation. The previous approach relies on the capability of providing at every iteration one supergradient vector of the current iterate. It is straightforward to obtain the following characterization of the supergradient of the discrete dual cost Φ: if and only if it is a convex combination of the finite set of extremal vectors defined in the following way. Let ϕ k(i,j) * i be any selection of minimizing values involved in the dual cost. That is ∀i, j ϕ

k(i,j) * i ∈ argmin k=1,...,N k λ i 2 |y j i -z k | 2 -ϕ k i . (4.13)
Then, the set of extremal vectors is the finite collection of all vectors

(e k * ) = ( j k -c(N i )ν j i δ k(i,j) * (k)) (4.14)
for any selection (k(i, j) * ).

A crucial observation that has been raised in [START_REF] Mérigot | Discrete optimal transport: complexity, geometry and applications[END_REF] is the fact that the computation of a vector of the superdifferential does not require generically an order of I × N k × N i operations. Actually the following formulation makes it possible to use a special data structure called kd-tree which in most cases reduces the complexity of finding one vector of ∂Φ((ϕ k 1 ), . . . , (ϕ k I )).

Let i, j be given and assume we want to evaluate the minimal value

M = min k λ i 2 |y j i -z k | 2 -ϕ k i .
Let us then define c i = max k ϕ k i and

M = -c i + min k λ i 2 |y j i -z k | 2 + c i -ϕ k i .
Since the latter term is positive:

M = -c i + min k ||P j i -Q k i || 2
where P j i = (

λ i 2 y j i , 0) is a fixed vector of R d+1 , Q k i = ( λ i 2 z k , c i -ϕ k i
) and ||.|| stands for the standard euclidean norm of R d+1 . Thus our supergradient request reduces to identify one closest point of P j i among points of (Q k i ) k . Observe that the family (Q k i ) k does not depend of the parameter j. This task is a standard operation in computational geometry which can be performed efficiently with kd-tree structures. By using such tools, we can reduce the complexity of the supergradient request in the generic case to an order of I × log N k × N i operations. Observe that if the N i 's and N k are of same order N , one request is generically of complexity of order I 2 N log N .

Numerical quantization and localization. We suppose that all the measures µ i are compactly supported. We use the discretization of Section 2.4. The support of the unknown barycenter measure is bounded using the results of Section 2.2, in particular, (2.10) Reconstruction of the barycenter density. One additional difficulty associated to the dual formulation is the fact that optimal dual vectors only give an implicit description of the barycentric measure. In order to recover the support and the density of the barycentric measure, we introduce the following least square procedure.

Every optimal dual vector (ϕ k i ) k must be associated to an optimal transport from µ i to the barycentric measure. A crucial observation is the fact that every associated map transports the µ i to the same measure. We exploit this optimality condition to recover the barycentric measure through the coefficients f j,k i which describe the mass transported from y j i to z k . By optimality, some mass can be transported from y j i to z k if and only if

G i,j (ϕ k i ) = min k ′ =1,...,N k G i,j (ϕ k ′ i ). (4.15) where G i,j (ϕ k ′ i ) = λi 2 |y j i -z k | 2 -ϕ k ′ i .
Let us fix some parameter ε > 0. Based on the previous observation, we only consider the unknown coefficients f j,k i for Figure 9: Classical McCann interpolation between translated measures which G i,j (ϕ k i ) is less than the optimal value (4.15) plus ε. In order to recover the barycentric measure, we look for the set of coefficients which generate the same measures in an optimal least square sense. More precisely, we solve the sparse least square problem: min

f j,k i l,m k ( j f j,k l ν j l - j ′ f j ′ ,k m ν j ′ m ) 2 (4.16)
among non negative coefficients less than one which satisfy the linear constraints ∀i, j, k f j,k i = 1.

Numerical results

As detailed in the previous paragraphs, our approach relies first on a non smooth optimization step using an hybrid LBFGS/Bundle algorithm and a fast computation of supergradient vectors. In a second phase, a sparse least square problem is solved in order to recover an approximation of barycentric density.

Let us point out that in all the following examples, the first optimization step was stopped after one hour of computation on a standard PC. This costly step could be dramatically sped up in using a straightforward parallelized cost evaluation. We validate our approach by considering different test cases for which analytic descriptions of barycenters are available. The simplest situation is the case of the barycenters of a measure of density ρ(.) and a translated measure of density ρ(.+V ) where V is some fixed vector. In this trivial case, the isobarycentric measure is of course the measure of density ρ(. + V /2). We display in figure 9, the barycentric measure obtained by our located approximation scheme applied to ρ = χ c where χ c is the characteristic function of a unit square. In this experiment, we used a grid of size 200 × 200 and a recovering parameter ε = 10 -5 . In the least square optimization problem (4.16), we obtain an error of order 10 -4 for every quadratic term.

Next we applied our approach to the case of Gaussian measures. A complete description of barycenters of Gaussian measures has been given in [AC11]: consider a family of Gaussian measures µ i (m i , S i ) of means (m i ) i and covariance matrices (S i ) i . Then, the barycentric measure associated to the non-negative weights (λ i ) i is a Gaussian measure of mean the barycenter of the (m i , λ i ) i . Moreover, its covariance matrix is the only definite positive matrix S solution of the equation i λ i (S 1/2 S i S 1/2 ) 1/2 = S.

(4.17)

Figure 10: Isobarycenter computation of three gaussian measures by a global (first row) and a localized approach (second row).

We denote by N (m i , σ i ) a Gaussian of mean m i and of covariance matrix equal to σ 2 i Id. We considered two different test cases and applied for both our global and localized approaches. In all the experiments the number of sampling points of the given measures and of the barycentric measure have been fixed for the global approach to ∀i, N i = N z = 15 × 10 3 . For the localized approach by Minkowski sum, we imposed ∀i, N i = N z = 5 × 10 3 . The first test case consists in approximating the isobarycenter of the three Gaussian measures of random standard variations N ((0.1, 0.8), 1/49.75), N ((-0.9, -1), 1/35.89) and N ((1, -0.9), 1/74.63). The resulting barycenter and the given measures are presented in figure 10. Our second test case is related to the approximation of the barycenter of the five gaussian measures (N (m i , σ i ), λ i ) i=1,...,5 where the m i are the vertices of a regular pentagon with σ i = 1/50 and λ i = 1/7 if i is odd and σ i = 1/100 and λ i = 2/7 otherwise. The resulting barycenter and the given measures are drawn in figure 11. Observe that in the localized illustrations, the support of the unknown measure is not anymore centered due to the loss of symmetry in Minkowski's sum. To conclude the study of those test cases, we provide in table 4.4, the errors between the theoretical and computed means and covariance matrices. As expected, even if the number of degree of freedom is smaller, the results obtained by the located approach are significantly better than the ones obtained by the first algorithm.

To conclude our numerical experiments, we provide large scale examples in which we interpolate three textures of images of size 150 × 150. This type of applications have been first studied in the framework of optimal transportation in [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF]. The texture mixing problem consists in synthesizing a texture from a family of given textures. The interest of using Wasserstein barycenters in this context is due to the spatial nature of Wasserstein distance which provides a more natural interpolation process than the naive pointwise means (see figures 12 and 13). We carried out similar experiments as the one depicted in [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF]. Our contribution with respect to [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] lies in the fact that we do not replace the quadratic Wasserstein distance by the easier to handle so-called sliced Wasserstein distance (which is an average over directions of one dimensional Wasserstein distances). We obtained by our method an approximation of the original model up to an error of 10 -3 for every quadratic term of (4.16).

Acknowledgements: The authors are grateful to the hospitality of BIRS where the present work was initiated at the occasion on a focused group meeting on Numerical methods for optimal transport. They are happy to thank Martial Agueh, Jean-David Benamou and Brendan Pass for many fruitful conversations. The first author gratefully acknowledges the support of the ANR, through the projects ISOTACE and OPTIFORM and INRIA through the "action exploratoire" MOKAPLAN.

  z) = 0, ∀z ∈ Z.

Figure 1 :Figure 2 :

 12 Figure 1: Given measures in black. Figure left, centre, right: solution of the geodesic problem with weights .25, .5, .75, respectively.

Figure 3 :Figure 4 :

 34 Figure 3: The measures m 1 , m 2 , m 3 , m 4 used in the examples which follow.

Figure 5 :Figure 6 :

 56 Figure 5: Barycenters of two rectangles m 1 , m 2 , with cost C(x, y) = |x -y| p for p = 1, 2, 3. Using grid size 100 2 .

Figure 7 :

 7 Figure 7: Barycenters of the measures m 1 , m 2 , m 3 , m 4 (four rotated rectangles) with cost C(x, y) = |x -y| p for p = 1, 2, 3. Using grid size 50 2 .

Figure 8 :

 8 Figure 8: Solution of the matching for teams problem. Top: the three measures, and the solution. Bottom: surface plot of the solution. The solution concentrates mass mostly on the corners with some mass on the edges.

  Proposition 4.1. Let (ϕ k 1 ), . . . , (ϕ k I ) ∈ R N k ×I . Then ((v k 1 ), . . . , (v k I )) ∈ ∂Φ((ϕ k 1 ), . . . , (ϕ k I ))

Figure 11 :

 11 Figure 11: Non uniform barycenter computation of five gaussian measures by a global (first row) and a localized approach (second row).

Figure 12 :

 12 Figure 12: Isobarycenter of the three textures of the first row. The pointwise mean of the three textures corresponds to the left picture of the second row. Wasserstein barycenter is presented on the right picture of the second row.

Table 1 :

 1 David Benamou and Yann Brenier. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math., 84(3):375-393, 2000. Upper bounds of the gap between theoretical and computed means and covariance coefficients

	[AC11]	Martial Agueh and Guillaume Carlier. Barycenters in the wasser-
		stein space. SIAM Journal on Mathematical Analysis, 43(2):904-
		924, 2011.
	[BB00]	Jean-||m th i -m i || |σ th i -σ i |
		First test case	0.003	0.005
		First test case localized	0.0002	0.0014
		Second test case	0.002	0.07
		Second test case localized	0.0002	0.0013