G Carlier
email: carlier@ceremade.dauphine.fr

A Galichon
email: alfred.galichon@polytechnique.edu

F Santambrogio

From Knothe's transport to Brenier's map and a continuation method for optimal transport

Keywords: optimal transport, rearrangement of vector-valued maps, Knothe-Rosenblatt transport, continuation methods

A simple procedure to map two probability measures in R d is the so-called Knothe-Rosenblatt rearrangement, which consists in rearranging monotonically the marginal distributions of the last coordinate, and then the conditional distributions, iteratively. We show that this mapping is the limit of solutions to a class of Monge-Kantorovich mass transportation problems with quadratic costs, with the weights of the coordinates asymptotically dominating one another. This enables us to design a continuation method for numerically solving the optimal transport problem.

Introduction

Given two Borel probability measures µ and ν on R d , a Borel map S : R d → R d is said to be a transport map between µ and ν if S µ = ν where S µ denotes the push-forward (or image measure) of µ through ν (i.e. S µ(B) = µ(S -1 (B)) for every Borel B). In the present article, we will focus on two particular transport maps: the Knothe-Rosenblatt transport and the Brenier's map.

The Knothe-Rosenblatt transport. The Knothe-Rosenblatt rearrangement was independently proposed by Rosenblatt [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF] for statistical purposes and by Knothe [4] in order to extend the Brunn-Minkowski inequalities. The principle is the following, as explained in Villani [START_REF] Villani | Optimal transport: Old and New, lecture notes, Ecole d'été de probabilités de Saint-Flour[END_REF]. Let µ and ν be two Borel probability measures on R d and assume for simplicity for the moment that µ is absolutely continuous with respect to the Lebesgue measure. Let us denote by µ d (respectively ν d ) the d-th marginal of µ (respectively ν) and µ d-1

x d , µ d-2 (x d ,x d-1 ) ,. . . , µ 1 (x d ,...,x 2 ) (respectively ν d-1 y d , ν d-2

(y d ,y d-1 ) ,. . . , ν 1 (y d ,...,y 2 ) ) the successive disintegrations (or conditional measures) of µ (respectively ν) given x d , (x d , x d-1 ),..., (x d , ..., x 2 ) (respectively given y d , (y d , y d-1 ), ..., (y d , ..., y 2 )). Now let T d = T d (x d ) be the monotone nondecreasing map transporting µ d to ν d , such a map is well-defined and unique as soon as µ d has no atoms and in this case, it is explicitly given by T d = G T d (x d ) . One repeats the construction (well-known by statisticians under the name of conditional quantile transforms) iteratively until we define T 1 (x 1 , x 2 , ..., x d ), which is monotone in x 1 and transports µ 1 (x d ,...,x 2 ) onto ν 1

T 2 (x 2 ,...,x d ) . Finally, the Knothe-Rosenblatt rearrangement T is defined by

T (x) = (T 1 (x 1 , x 2 , ..., x d ), ..., T d-1 (x d-1 , x d ), T d (x d ))
. Obviously, T is a transport map from µ to ν, i.e. T µ = ν. By construction, the Knothe transport T has a triangular Jacobian matrix with nonnegative entries on its diagonal. Note also that the computation of the Knothe transport only involves one-dimensional monotone rearrangements and that it is well defined as soon the measures one transports have no atoms. The precise assumption is the following.

Assumption (H-source): the measure µ d , as well as µ d -almost all the measures µ d-1

x d , and the measures µ d-2 x d ,x d-1 for µ d -a.e. x d and µ d-1

x d -a.e. x d-1 . . . up to almost all the measures µ 1

x d ,...,x 2 , which are all measures on the real line, must have no atoms.

Notice that (H-source) is satisfied as soon as µ is absolutely continuous with respect to the Lebesgue measure.

The Monge-Kantorovich problem and the Brenier's map. Optimal transportation theory provides an alternative way to transport µ to ν. We recall that in the case of the quadratic cost, the Monge-Kantorovich problem reads as

inf π∈Γ(µ,ν) R d ×R d |x -y| 2 dπ(x, y) (1.1)
where Γ(µ, ν) denotes the set of transport plans between µ and ν i.e. the set of probability measure on R d × R d with marginals µ and ν. We refer to the books of Villani [START_REF] Villani | Topics in Optimal Transportation[END_REF], [START_REF] Villani | Optimal transport: Old and New, lecture notes, Ecole d'été de probabilités de Saint-Flour[END_REF] for a modern account of optimal transportation theory. The linear problem (1.1) is a relaxation of the Monge problem inf

S : S µ=ν R d |x -S(x)| 2 dµ(x) (1.2)
When µ is absolutely continuous with respect to the Lebesgue measure, Brenier [START_REF] Brenier | Polar factorization and monotone rearrangement of vectorvalued functions[END_REF] proved that (1.2) has a unique solution which is characterized by the fact that it is the gradient of some convex function. More precisely, there exists a unique (up to constants and µ-a.e. equivalence) convex function V : R d → R such that ∇V µ = ν. Also (id × ∇V ) µ is the only solution of (1.1) and ∇V is called the Brenier's map between µ and ν. The Knothe-Rosenblatt as a limit of optimal transportation plans.

Let us slightly modify the quadratic cost in (1.1) and replace it with the weighted quadratic cost

c ε (x, y) := d i=1 λ i (ε)(x i -y i ) 2
where the λ i (ε)'s are positive scalars depending on a parameter ε > 0. If µ is absolutely continuous with respect to the Lebesgue measure, the corresponding optimal transportation problem admits a unique solution T ε . When in addition, for all k ∈ {1, ..., d -1}, λ k (ε)/λ k+1 (ε) → 0 as ε → 0, it is natural to expect the convergence of T ε to the Knothe transport T . We will show that this convergence holds provided ν satisfies some additional condition, and namely

Assumption (H-target): the measure ν d , as well as ν d -almost all the measures ν d-1 y d , and the measures ν d-2 y d ,y d-1 for ν d -a.e. y d and ν d-1 y d -a.e. y d-1 . . . up to almost all the measures ν 2 y d ,...,y 3 , which are all measures on the real line, must have no atoms.

Notice that (H-target) is not natural as (H-source) is. Yet, we will show a counter-example to the convergence result when it is not satisfied. (Htarget) as well is satisfied should ν be absolutely continuous (actually, this assumption is slightly weaker then (H-source), since the last disintegration measures are not concerned).

This convergence result was conjectured by Y. Brenier as a very natural one, and actually its proof is not hard. Yet, it was not known before that extra assumptions on ν were needed. This makes one of the point of interest of this paper.

The other point is what we investigate later in the paper, i.e. the other direction: from Knothe to Brenier. We will study the dependence ε → T ε by means of the evolution with respect to ε of the dual variables. This will enable us, to design a numericaly strategy to approximate all the optimal transports T ε taking as initial condition the (cheap to compute) Knothe transport T .

An example. To illustrate the problem in a particular case where explicit solutions are available, take d = 2, and µ and ν two Gaussian measures where µ = N (0, I 2 ) and

ν = N 0, a b b c . Take λ 1 (ε) = ε and λ 2 (ε) = 1.
Then it can be verified that T ε is linear, and that its matrix in the canonical basis of R 2 is

T ε = 1 aε 2 + c + 2ε √ ac -b 2 aε + √ ac -b 2 b bε c + ε √ ac -b 2 which converges as ε → 0 to T = a -b 2 /c b/ √ c 0 √ c
, which is precisely the matrix of the Knothe transport from µ to ν. Organization of the paper. In section 2, we show, under suitable assumptions, that the optimal transportation maps for the cost c ε converge to Knothe's transport map as the parameter ε goes to 0, we will also emphasize that some conditions are to be imposed on ν for the convergence to hold. In section 3, we show that the evolution of the dual variables in the optimal transportation problem for cost the c ε is given by a well-posed ordinary differential equation. Finally in section 4, we discretize this equation and give several numerical results.

Knothe transport as a limit of quadratic optimal transports

We directly state our first result, whose proof, in the spirit of Γ-convergence developments (see [START_REF] Braides | Γ-convergence for beginners[END_REF]), will require several steps.

Theorem 2.1. Let µ and ν be two probability measures on R d satisfying (Hsource) and (H-target), respectively, with finite second moments, and γ ε be an optimal transport plan for the costs

c ε (x, y) = d i=1 λ i (ε)(x i -y i ) 2 , for some weights λ k (ε) > 0. Suppose that for all k ∈ {1, ..., d -1}, λ k (ε)/λ k+1 (ε) → 0 as ε → 0.
Let T be the Knothe-Rosenblatt map between µ and ν and γ K ∈ P(R d × R d ) the associated transport plan (i.e. γ K := (id × T ) µ). Then γ ε γ K as ε → 0. Moreover, should the plans γ ε be induced by transport maps T ε , then these maps would converge to T in L 2 (µ) as ε → 0.

Proof. Take the plans γ ε that are optimal for the Brenier-like cost c ε given by

c ε (x, y) = d i=1 λ i (ε)(x i -y i ) 2
(we suppose for simplicity λ d (ε) = 1 and λ i (ε)/λ i+1 (ε) → 0). Suppose (which is possible, up to subsequences) γ ε γ. We want to prove γ = γ K . By comparing γ ε to γ K and using optimality we first get

c ε dγ ε ≤ c ε dγ K (2.1)
and, passing to the limit as ε → 0, since c ε converges locally uniformly to

c (d) (x, y) = (x d -y d ) 2 , we get c (d) dγ ≤ c (d) dγ K .
Yet, the function c (d) only depends on the variables x d and y d and this shows that the measure (π d ) γ gets a better result than (π d ) γ K with respect to the quadratic cost (π d being the projection onto the last coordinates, i.e.

π d (x, y) = (x d , y d )
). Yet, the measure γ K has been chosen on purpose to get optimality from µ d to ν d with respect to this cost, and the two measures share the same marginals. Moreover, thanks to the assumptions on µ d , this optimal transport plan (which is actually induced by a transport map) is unique. This implies (π d ) γ = (π d ) γ K . Let us call γ d this common measure. We go back to (2.1) and go on by noticing that all the measures γ ε have the same marginals as γ K and hence their (separate) projection onto x d and y d are µ d and ν d , respectively. This implies that (π d ) γ ε must realize a result which is worse than (π d ) γ K as far as the quadratic cost is concerned and consequently we have

|x d -y d | 2 d(π d ) γ K (x d , y d ) + d-1 i=1 λ i (ε) (x i -y i ) 2 dγ ε ≤ c ε dγ ε ≤ c ε dγ K = |x d -y d | 2 d(π d ) γ K (x d , y d ) + d-1 i=1 λ i (ε) (x i -y i ) 2 dγ K ,
which implies, by simplifying the common term in d(π d ) γ K , dividing by λ d-1 (ε) and passing to the limit, c (d-1) dγ ≤ c (d-1) dγ K (we use the general notation c k (x, y) = |x k -y k | 2 ). We can notice that both integrals depend on the variables x d-1 and y d-1 only. Anyway, we can project onto the variables (x d-1 , x d ) and (y d-1 , y d ) (obtaining measures (π d-1 ) γ and (π d-1 ) γ K ) so that we disintegrate with respect to the measure γ d . We have

dγ d (x d , y d ) |x d-1 -y d-1 | 2 dγ d-1 (x d ,y d ) (x d-1 , y d-1 ) ≤ dγ d (x d , y d ) |x d-1 -y d-1 | 2 dγ d-1 (x d ,y d ),K (x d-1 , y d-1 ). (2.2)
It is is sufficient to prove that the measures γ d-1 (x d ,y d ) share the same marginals on x d-1 and y d-1 as the corresponding γ d-1

(x d ,y d ),K to get that their quadratic performance should be worse than the corresponding performance of γ d-1

(x d ,y d ),K
(this is because the Knothe measure has been chosen exactly with the intention of being quadratically optimal on (x d-1 , y d-1 ) once x d and y d are fixed). Yet, (2.2) shows that, on average, the result given by the those measures is not worse than the results of the optimal ones. Thus, the two results coincide for almost any pair (x d , y d ) and, by uniqueness of the optimal transports (this relies on the assumptions on the measures µ d-1

x d ), we get γ d-1

(x d ,y d ) = γ d-1 (x d ,y d ),K .
To let this proof work it is sufficient to prove that the projections of the two measures coincide for γ d -a.e. pair (x d , y d ). For fixed (x d , y d ) we would like to prove, for any φ φ(x d-1 )dγ d-1

(x d ,y d ) = φ(x d-1 )dγ d-1 (x d ,y d ),K
(and to prove an analogous equality for functions of y d-1 ). Since we accept to prove it for a.e. pair (x d , y d ), it is sufficient to prove this equality:

dγ d (x d , y d )ψ(x d , y d ) φ(x d-1 )dγ d-1 (x d ,y d ) = dγ d (x d , y d )ψ(x d , y d ) φ(x d-1 )dγ d-1 (x d ,y d ),K
for any φ and any ψ. This means proving

ψ(x d , y d )φ(x d-1 )dγ d-1 = ψ(x d , y d )φ(x d-1 )dγ d-1 K ,
which is not trivial since we only know that the two measures γ d-1 and γ d- ) is actually a function of (x d-1 , x d ) only, and that equality holds when passing from γ to γ K .The same can be performed on functions ψ(x d , y d )φ(y d-1 ) but we have in this case to ensure that we can replace x d with a function of y d , i.e. that we can invert T d . This is possible thanks to the assumption on ν d , since T d is the optimal transport from µ d to ν d , but an optimal transport exists in the other direction as well and it gives the same optimal plan (thanks to uniqueness). These facts prove that the measures γ d-1 (x d ,y d ) and γ d-1 (x d ,y d ),K have the same marginals and hence, since they are both optimal, they coincide for a.e. pair (x d , y d ). This implies

γ d-1 = γ d-1
K . Now, it is possible to go on by steps: once we have proven that γ h = γ h K , let us take (2.1) and estimate all the terms with |x i -y i | 2 and i ≥ h thanks to the optimality of γ K , thus getting

i≥h λ i (ε) |x i -y i | 2 dγ K + h-1 i=1 λ i (ε) (x i -y i ) 2 dγ ε ≤ c ε dγ ε ≤ c ε dγ K = i≥h λ i (ε) |x i -y i | 2 dγ K + h-1 i=1 λ i (ε) (x i -y i ) 2 dγ K ,
and consequently, by dividing by λ h-1 (ε) and passing to the limit,

c (h-1) dγ ≤ c (h-1) dγ K .
We disintegrate with respect to γ h and we act exacly as before: proving that the marginals of the disintegrations coincide is sufficient to prove equality of the measures. Here we will use test-functions fo the form

ψ(x h , x h+1 , . . . , x d , y h , y h+1 , . . . , y d )φ(x h-1 )
and ψ(x h , x h+1 , . . . , x d , y h , y h+1 , . . . , y d )φ(y h-1 ).

The same trick as before, i.e. replacing the variables y with functions of the variables x is again possible. To invert the trick and replace x with y one needs to invert part of Knothe's transport. This is possible since our assumptions imply that all the monotone transports we get are invertible.

In the end we get, as before, γ h-1 = γ h-1 K . This procedure may go on up to h = 2, thus arriving at γ = γ K .

We have now proven γ ε γ K . Yet, if all these transport plans come from transport maps, it is well known that (T ε × id) µ (T × id) µ implies T ε → T in L p (µ), for any p > 1, as far as T ε is bounded in L p (µ). Actually, weak convergence is a simple consequence of boundedness: to go on, we can look at Young's measures. The assumption (the limit is a transport map as well) exactly means that all the Young measures are dirac masses, which implies strong convergence. In particular we get L 2 (µ) convergence and µ-a.e. convergence on a subsequence.

Let us remark here that if instead of considering the quadratic cost c ε , one considers the more general separable cost

c ε (x, y) := d i=1 λ i (ε)c i (x i -y i )
where each c i is a smooth strictly convex function (with suitable growth), then the previous convergence proof carries over.

A counterexample when the measures have atoms We now show that interestingly, and perhaps counterintuitively, the hypothesis of absence of atoms in theorem 2.1 is necessary not only for µ, but also for ν. We propose a very simple example in R 2 where µ is absolutely continuous with respect to the Lebesgue measure but ν does not satisfy (H-target), and we show that the conclusion of theorem 2.1 fails to hold. On the square Ω := [-1, 1] × [-1, 1], define µ such that µ(dx) = 1 {x 1 x 2 <0} dx/2 so that the measure µ is uniformly spread on the upper left and the lower right quadrants, and

ν = H 1 |S /2, being S the segment [-1, 1] × {0}.
The Knothe-Rosenblatt map is easily computed as (y 1 , y 2 ) = T (x) := (2(x 1 + sgn(x 2 )), 0). The solution of any symmetric transportation problem with λ ε = (ε, 1) is (y 1 , y 2 ) = T 0 (x) := (x 1 , 0) (no transport may do better than this one, which projects on the support of ν). Therefore, in this example the optimal transportation maps fail to tend to the Knothe-Rosenblatt map. The reason is the atom in the measure ν 2 = δ 0 .

Convergence even with atoms The convergence result of theorem 2.1 requires the absence of atoms in the projections of ν. This is obviously not the case if ν itself is purely atomic! Yet, this will precisely be the case we will consider in the algorithm we propose in the sequel. The same proof may be extended to this case under the following assumption. Keep the same assumptions on µ but suppose that ν is concentrated on a set S with the property y, z ∈ S, y

= z ⇒ y d = z d .
This means that, if we restrict ourselves to S, then all the variables y i for i < d are actually a function of the last variable y d . This is particularly useful when ν is purely atomic, concentrated on a finite (or countable) set of points with different y d components.

Just come back to the proof. The equality

dγ d (x d , y d )ψ(x d , y d ) φ(x d-1 )dγ d-1 (x d ,y d ) = dγ d (x d , y d )ψ(x d , y d ) φ(x d-1 )dγ d-1 (x d ,y d ),K
only relied on y d being a function of x d , which is still true. The other equality, namely

dγ d (x d , y d )ψ(x d , y d ) φ(y d-1 )dγ d-1 (x d ,y d ) = dγ d (x d , y d )ψ(x d , y d ) φ(y d-1 )dγ d-1 (x d ,y d ),K
gives some extra troubles. It is not any more true that x d is a function of y d . Yet, it is true that y d-1 is a function of y d and this allows us to reduce the expression to functions of (x d , y d ) only, which is sufficient to get equality. The same procedure may be performed at subsequent steps as well.

An ODE for the dual variables

In this section, we consider for simplicity the case d = 2 (although our method extends to higher dimensions), µ uniform on some convex polyhedron Ω (for the sake of simplicity we will assume |Ω| = 1) and ν = 1 N N i=1 δ y i where all the points y i ∈ Ω have a different second coordinate y

(2) i . For every ε ≥ 0, let A ε be the diagonal 2 × 2 matrix with diagonal entries (ε, 1) and let c ε be the quadratic cost defined by c ε (x, y) = A ε (x -y)(x -y). We are interested in solving the family of optimal transportation problems inf

π∈Γ(µ,ν) R d ×R d c ε (x, y)dπ(x, y) (3.1)
for all values of the parameter ε ∈ [0, 1]. It is well-known, that (3.1) can be conveniently solved by the dual problem formulated in terms of prices:

sup p Φ(p, ε) := 1 N N i=1 p i + Ω p * ε (x)dx, (3.2)
where p * ε (x) = min i {c ε (x, y i ) -p i } and we impose as a normalization p 1 = 0. For each ε, there is a unique maximizer p(ε). For each (p, ε) we define

C(p, ε) i = {x ∈ Ω : inf j c ε (x, y j ) -p j = c ε (x, y i ) -p i }. It is easy to check that Φ ε := Φ(., ε) is concave differentiable and that the gradient of Φ ε is given by ∂Φ ε ∂p i (p) = 1 N -|C(p, ε) i |.
By concavity of Φ ε , the solution p(ε) of (3.2) is characterized by the equation ∇Φ ε (p(ε)) = 0. The optimal transportation between µ and ν for the cost c ε is then the piecewise map taking value y i in the cell C(p(ε), ε)) i . Our aim is to charcacterize the evolution of p(ε) as ε varies. Formally, differentiating the optimality condition ∇Φ(p(ε), ε) = 0, we obtain a differential equation for the evolution of p(ε):

∂ ∂ε ∇ p Φ(p(ε), ε) + D 2 p,p Φ(p(ε), ε) • dp dε (ε) = 0. (3.3)
Our aim now is to show that the equation (3.3) is well-posed starting with the initial condition p(0) (corresponding to horizontal cells of area 1/N ); this will involve computing the second derivatives of Φ in (3.3), proving their Lipschitz behavior as well as obtaining a negative bound on the larger eigenvalue of the negative semidefinite matrix D 2 p,p Φ . The price vector p(ε), along the evolution, will always be such that all the areas |C(p, ε) i | are equal (and are equal to 1/N ). Yet, we need to prove that the differential equation is well posed and we will set it in an open set,

O = (p, ε) : 1 2N < |C(p, ε) i | < 2 N for every i . (3.4)
The initial datum of the equation will be such that |C(p(0), 0) i | = 1/N and we will look at the solution only up to the first moment where it exits O.

Yet inside the set it will be well-posed and it will imply conservation of the areas. Hence we will never exit O.

All the quantities we are interested in depend on the position of the vertices of the cells C(p, ε) i , which are all polygons. Let us call x(p, ε) ± i,j the two extremal points of the common boundary between C(p, ε) i and C(p, ε) j , that we call D(ε, p) i,j (if such a common boundary exists; should it be a single point we consider the two points as coinciding). Each one of this points is obtained as the intersection of at least this common boundary with another one, or with the boundary of Ω (which is supposed to be a polygon as well, so that in the neighbourhood of almost any point locally the boundary is a line). We want to investigate the dependence of these points with respect to (p, ε) and prove that this dependence is Lipschitz. Notice that each point x(p, ε) ± i,j is not defined for any value of (p, ε) but only on a certain (closed) subset of the space R N × (0, 1).

Lemma 3.1. The positions of the vertices x(p, ε) ± i,j depend in a Lipschitz way on p and ε.

Proof. Locally it is true that the same point x(p, ε) ± i,j is defined either by a system of equations

A ε (x -y i )(x -y i ) -p i = A ε (x -y j )(x -y j ) -p j , A ε (x -y i )(x -y i ) -p i = A ε (x -y h )(x -y h ) -p h (3.5)
in the case of a point on the intersection of two common boundaries, or by a system

A ε (x -y i )(x -y i ) -p i = A ε (x -y j )(x -y j ) -p j , Lx = l 0 (3.6)
in the case of intersection with the boundary ∂Ω (locally being given by the equation Lx = l 0 ) . The first system, after simplifying, reads as

2A ε (y j -y i )(x) = A ε (y j )(y j ) -A ε (y i )(y i ) -p j + p i , 2A ε (y h -y i )(x) = A ε (y h )(y h ) -A ε (y i )(y i ) -p h + p i , (3.7)
and the second as well may be simplified the same way. This means that they are of the form

M (ε)x = V (ε, p)
for a matrix M (ε) which reads, in usual coordinates,

M (ε) = ε(y (1)
j -y

i ) y

(2)

j -y (2) i ε(y (1)
h -y

i ) y

(2)

h -y (2) i
in the first case and

M (ε) = ε(y (1)
j -y (1)
i ) y

(2)

j -y (2) i l 1 l 2 in the second.
The vector V (ε, p) is obtained regarding the right hand sides of the system, as in (3.7). Both M and V depend Lipschitzly on (ε, p), with uniformly bounded Lipschitz constants. Hence, to check that the dependence of x on (ε, p) is Lipschitz, we only need to bound (away from zero) det M (ε). The determinant of a 2 × 2 matrix is given by the product of the modulus of the two vector composing its lines, times the sinus of the angle between them.

In the first case the vectors are A ε (y j -y i ) and A ε (y h -y i ), while in the second they are A ε (y j -y i ) and l = (l 1 , l 2 ). In both cases they are the normal vectors to the sides of the cell we are considering. This implies, thanks to the lower bound on the areas of the cells, that the angles between these vectors may not be too small. Actually, since Ω is bounded and the cells are convex, the area of each cell is smaller than (diamΩ) 2 sin α/2, α being any angle between two neighbour sides. This implies a lower bound on sin α.

The lower bound on the moduli comes from the fact that the vectors l are chosen with modulus one and the modulus of A ε (y j -y i ) is always greater than its vertical component, which is y

(2) j -y (2)
i , which is supposed different from zero for any pair (i, j). Notice that in the first case the matrix has ε in the determinant, even if we proved a lower bound on such a determinant, independent on ε: this agrees with the fact that actually, this kind of crossing (between two common boundaries of two cells) will only happen for ε ≥ ε 0 (for ε < ε 0 we only have almost horizontal strips crossing non-horizontal sides of Ω).

Lemma 3.2. The function Φ admits pure second derivatives with respect to p and mixed second derivatives with respect to p and ε, and these second derivatives are Lipschitz continuous.

Proof. We have proven that the positions of the points x(p, ε) ± i,j depend Lipschitzly, with uniformly bounded Lipschitz constants, on p and ε. Since the volumes of the cells C(p, ε) i are Lipschitz functions of these points, this implies that ∇ p Φ is C 0,1 . Hence it admits derivatives almost everywhere and we can compute them in the following way.

The derivative of a volume of a polygonal cell is given, side by side, by the length of the side times the average of the components which are normal to such a side of the two derivatives ẋ(p, ε) + i,j and ẋ(p, ε) - i,j (the other terms -which are mainly the terms at the corners -are of higher order). The equation of the side D(ε, p) i,j is, as we know,

2A ε (y j -y i )(x) = A ε (y j )(y j ) -A ε (y i )(y i ) -p j + p i (3.8)
and the normal unit vector to the side is

n = A ε (y j -y i ) |A ε (y j -y i )| .
Let us start from the derivatives of the cell C(p, ε) i with respect to a variable p j with j = i. We differentiate (3.8) with respect to p j and we get

2A ε (y j -y i )( ẋ) = -1.
This formula only works for x = x(p, ε) + i,j and x = x(p, ε) - i,j . Obviously it ony works where they are differentiable, i.e. almost everywhere. Hence the derivative, by summing up and rescaling the normal vector, is given by

∂|C(p, ε) i | ∂p j = - l i,j 2|A ε (y j -y i )| , (3.9)
where l i,j is the length of D(ε, p) i,j .

As far as the derivative with respect to p i is concerned, it is not difficult to check that we have (by summing up the results on every side)

∂|C(p, ε) i | ∂p i = j l i,j 2|A ε (y j -y i )| , (3.10)
where the sum is performed on all the indices j such that the cell C(p, ε) j is in contact with the cell C(p, ε) i . Automatically, since these derivatives only depend on the values of l i,j , which depend in a Lipschitz manner on the positions of x = x(p, ε) ± i,j , they are Lipschitz as well. This proves that Φ ε is actually C 2,1 and that these derivatives are well defined and admit the previous expressions (3.9)-(3.10) everywhere.

The computation of the derivatives with respect to ε is a bit trickier. We derive again (3.8), but with respect to ε. Since dA ε /dε = B, we get 2A ε (y j -y i )( ẋ) = -2B(y j -y i )(x) + B(y j )(y j ) -B(y i )(y i ) = 2B(y j -y i )

y j + y i 2 -x .
Then we renormalize the normal vector, sum up the results for x = x(p, ε) + i,j and x = x(p, ε) - i,j , multiply by the lengths and sum up the results for all the sides, and get

∂|C(p, ε) i | ∂ε = j l i,j B(y j -y i )(y j + y i -x(p, ε) + i,j -x(p, ε) - i,j ) 2|A ε (y j -y i )| . ( 3

.11)

In this case as well the result is Lipschitz in (p, ε) and hence ∇ p Φ is differentiable everywhere with respect to ε, with Lipschitz derivative.

We can come now back to the evolution of p = p(ε) and consider again the differential equation (3.3). To solve this equation we need to prove that the matrix D 2 p,p Φ is actually invertible (for numerical purpose, we will also need to bound its eigenvalues away from zero). It is important to recall that we look at the evolution of the vector p = (p 2 , . . . , p N ), since we may assume p 1 (ε) = 0 for all ε. Hence, we will not look at the entries 1 in the vectors or the matrices. The matrix we consider is M := -(D 2 p,p Φ) i,j=2,...,N has the following properties:

• on each line, outside the diagonal we have negative terms M i,j = -

l i,j
2|Aε(y j -y i )| ;

• each element on the diagonal is the sum of minus all the others on the same line (hence it is positive), and possibly of the term which should be in the same line at the first column;

• an entry (i, j) of the matrix is non-zero if and only if the cells C(p, ε) i and C(p, ε) j share a common boundary with positive length;

• in particular, for any pair (i, j), even if the entry at place (i, j) is zero, it is possible to find a path i = i 0 , i 1 , i 2 , . . . , i k = j so that the matrix has non-zero values at all the positions (i h , i h+1 );

• the total of the entries of the first column (the one which is not present in the matrix) is strictly positive.

The invertibility of M is ensured by the following:

Lemma 3.3. Let the matrix M satisfy the following properties

(H1) for all i, M i,i ≥ j =i |M i,j |, (H2) there exists i such that M i,i > j =i |M i,j |, (
H3) for any pair (i, j) there is a sequence i 0 , i 1 , i 2 , . . . , i k with i 1 = i, i k = j, and M i h ,i h+1 = 0.

then M is invertible.

Proof. Let x ∈ Ker(M ) and let i be an index such that |x i | is maximal. We may suppose for simplicity that x i is positive. Then we have

0 = M i,i x i - j M i,j x j ≥ M i,i x i - j M i,j x i = x i M i,i - j M i,j ≥ 0.
This implies that all inequalities are equalities and in particular x j = x i whenever M i,j = 0. Hence, the entries of x on all the indices which are "neighbours" of i equal x i (and they are maximal as well). This allows to repeat the argument replacing i with another maximizing index j and so on... since any index is connected by a chain of neighbours to i, we get that all the entries are equal. But this implies that the vector in the kernel we selected must be a multiple of the vector (1, 1, . . . , 1). Yet, this vector is not in the kernel since the sum of the elements on each line is not zero for all lines, by assumption (H2). This proves that M is invertible.

Finding a lower bound for the modulus of the eigenvalues of M , i.e. quantifying its inversibility is not straightforward. Indeed, the properties (H1), (H2), (H3) are not sufficient to get this bound, even if we fix the norm of the remainding column as the following counter-example shows. The determinant of the matrices

M ε =   1 -ε 0 -ε 1 -(1 -ε) 0 -(1 -ε) 1   is 2ε(1 -ε) → 0,
which implies that some eigenvalue as well goes to 0. We will obtain a positive lower bound by a compactness argument, but we will use something stronger than simply assumptions (H1), (H2), (H3). The idea is that assumption (H3) is not closed, but it stays closed when we replace it with the stronger condition of the matrix being associated to a partition (as is the case for M = -D 2 p,p Φ). In this case if one connection degenerates (i.e. a common boundary reduces to a point), some other connections will play the role.

Lemma 3.4.

There is a positive uniform lower bound on the least eigenvalue of any matrix M associated to the cell partition corresponding to a pair (p, ε) ∈ O.

Proof. The proof will be obtained by contradiction. To this aim, take a sequence of partitions of Ω into sets (Ω n i ) i=1,...,N . We assume these sets to be convex polygons with a bounded number of sides. This is the case for the cells associated to pairs (p, ε). We also know that their areas are always bounded between 1/2N and 2/N ). These partitions give rise to a certain topology of connections between the cells. Up to subsequences, we may suppose that this topology is always the same on all the partitions of the sequence (since the number of possible topologies is finite). Up to subsequences, we also have convergence in the Hausdorff distance. This means that for any i we have Ω n i → Ω i , and this convergence, which is the same as the convergence of all the vertices, preserves the areas, the convexity, the upper bound on the number of sides, the fact of being a partition... The matrices associated to these partitions depend continuously on these sets (with respect to this convergence, since they actually depend on the positions of the vertices). Notice that it is possible that a side reduces its length along the sequence up to becoming a single point in the limit. Yet, two cells which share a boundary along the sequence will do it along the whole sequence (thanks to our choice of not changing the topology) and at the limit, either they share a side as well, or they share a point only, but in this case the terms l i,j converged to zero. Hence we can associate to all the partitions their matrices and this correspondence is continuous. We have a sequence of matrices M n → M and let us suppose that some eigenvalue λ (n) 1 goes to zero. This would imply that the matrix M is associated to a partition but has a zero eigenvalue. This is not possible, thanks to the proof of lemma 3.3 above. M is associated to a partition and hence satisfies assumptions (H1) and (H3). To check (H2) we observe that the column that we remove, the first one, is associated to the first cell and, up to some rescaling but bounded factor |A ε (y 1 -y j )|, its entries are the lengths of its sides. Yet, this cell conserves the area bounds we had on the sequence and its entry cannot be all zero.

From the previous results on the form and the regularity of the derivatives of ∇ p Φ, we deduce from the Cauchy-Lipschitz Theorem that the ODE (3.3) governing the evolution of the dual variables is well posed and actually characterizes the optimal prices: Theorem 3.5. Let p(ε) be the solution of the dual problem (3.2) (recall the normalization p 1 (ε) = 0), then it is the only solution of the ODE:

dp dε (ε) = -D 2 p,p Φ(p(ε), ε) -1 ∂ ∂ε ∇ p Φ(p(ε), ε) (3.12)
with initial condition p(0) such that all the horizontal strips C(p(0), 0) i have area 1/N .

4 Numerical results

Algorithm

The algorithm we propose consists simply in discretizing (3.12) (together with the initial condition p(0) determined as in Theorem 3.5) by an explicit Euler scheme. Let n be some positive integer and h := n -1 be the step size. Let us set p 0 = p(0) and define prices inductively as follows.

• While (p k , kh) belongs to the open set O defined by (3.4), compute:

A k := -D 2 p,p Φ(p k ; kh), δ k := ∂ ∂ε ∇ p Φ(p k , kh).
Note that computing A k and δ k by formulas (3.9)-(3.10) and (3.11), requires to construct the cells C(kh, p k ) i .

• Solve the linear system A k z = δ k ; by taking advantage of A k being positive definite, we use the conjugate gradient algorithm for minimizing

J k (z) = A k (z)(z) -2δ k • z to solve this system exactly in N -1 steps.
We denote by z k the solution.

• Update the prices by setting

p k+1 = p k + hz k .
Thanks to the Lipschitz properties established in section 3, it is easy to check that for h small enough, (p k , kh) always remain in O and then p k is well-defined for every k up to n. For such an h and since (3.12) is generated by a Lipschitz function on O, it is well-known that the convergence of the Euler scheme is linear (see for instance [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF]). Denoting by p h the piecewise constant function having values p k on intervals [kh, (k + 1)h), we thus get the following convergence: Theorem 4.1. For h small enough, the algorithm above is well-defined and the uniform error between p h and the optimal price p is O(h).

Numerical experiments

The construction of the cells at each step is achieved efficiently by an implementation in Matlab. In the setting described above where

d = 2, Ω = [0, 1] 2 , µ is the uniform distribution on Ω, ν = 1 N N k=1 δ y k , our algorithm computes the cells Ω ε k = {x ∈ Ω : T ε (x) = y k }
as well as the prices p ε k of the cell y k for ε = 0 to ε = 1. For ε = 0 which is the case where the transportation plan is the Knothe one, the computational task amounts to sorting the second component of the y k 's. Appropriate discretization steps are then chosen for the transition ε = 0 to ε = 1. At each step, the tesselation of Ω into the polyhedral cells Ω ε k is computed based on the prices p ε k . Adjacency information on these cells is computed, as well as the length of the facet between two cells and the coordinates of its extreme points. This information allows one to formulate a discretized version of ODE (3.12) using an Euler discretization scheme.

For geometric computations we use the Multi-Parametric Toolbox library, available online at http://control.ee.ethz.ch/ mpt. In particular, polytope computes the tesselation of Ω into the polyhedral cells Ω ε k . (A slighlty modified version of) mpt buildAdjacency extracts adjacency information on this tesselation, from which the vertices and the lengths of the sides between two cells can be deduced. The library also incorporates convenient graphical routines.

Error analysis Some numerical examples are presented below, for which relative errors in cell areas (i.e. deviation from the optimality conditions) are given as well as a comparison between the tesselations obtained with our method and the true solution. Another way to test our method is as follows. Let us consider the set

C := z(γ) := R 2 ×R 2 x 1 y 1 dγ, R 2 ×R 2 x 2 y 2 dγ , γ ∈ Γ(µ, ν) .
C is a closed convex subset of R 2 and it is strictly convex in the sense that its boundary contains no line segment. Denoting by γ ε the solution of (3.1), it is easy to check that z(ε) := z(γ ε ) is an extreme point of C and that εz 1 + z 2 = εz 1 (ε) + z 2 (ε) is the equation of a supporting line of C at z(ε) and this supporting line intersects C only at z(ε). If we consider the correlation curve ε ∈ (0, 1) → z(ε) it can be represented as the graph of a concave decreasing function whose slope (when it exists) at point z(ε) is -ε. In our numerical test, we will also present graphs comparing the true concave correlation curve to the one computed by our method.

We give three instances of executions of our algorithm, with samples of respectively 5, 10, and 15 points. Taking weights (ε, ε -1 ) with ε ∈ (0, +∞) (rather than (ε, 1) with ε ∈ (0, 1)) we get the full evolution of the optimal transports from one Knothe's transport (horizontal strips) to the other (vertical strips). Further examples as well as videos can be found at http://alfred.galichon.googlepages.com/anisotropic. It should also be pointed that our method does not approximate the solution of a single optimal transportation problem but a whole family of such problems (which actually explains relatively high running times).

Five sample points. We take as our sample set a sample of five points. We get the following errors:

# steps Relative errors in cell areas Time 100 -4.41% 2.66% 3.41% -2.46% 0.80% 66 s 500 -0.88% 0.54% 0.68% -0.49% 0.16% 349 s for which we draw in Figure 1 the partition obtained for ε = 1 using an exact method, as well as the true evolutions of the componentwise correlations of the x's and the y's from ε = 0 until ε = 1.

Running the continuation algorithm with 100 and 500 discretization steps we obtained partition sketched below, where the relative error on the cells area is inferior to 5%. The evolution of the componentwise correlations from ε = 0 until ε = 1 are also sketched and compared with their exact conterparts (the above curve). Running the algorithm with 500 discretization steps on a standard laptop took 349 seconds, and yields to the results below, where the maximal relative error on the cells area is 1%. The evolution of the componentwise correlations from ε = 0 until ε = 1 are also sketched and compared with their exact counterparts.

Ten sample points. Taking a sample of ten points, we obtain the following errors:

# steps Relative errors in cell areas Time 2,000 -1.37% -12% -2.44% -0.09% 2.27% ≤ 1 h 3,000 -0.92% -8.23% 1.63% 0.06% 1.52% ≤ 1,5 h (continued) 2,000 0.27% 5.71% 2.12% 7.36% -1.58% ≤ 1 h 3,000 0.18% 3.84% 1.42% 4.94% -1.05% ≤ 1,5 h

We draw in Figure 3 the partition obtained for ε = 1 using an exact method, as well as the true evolutions of the componentwise correlations of the x's and the y's from ε = 0 until ε = 1.

Fifteen sample points. With a sample of 15 points, we get the following errors:

#

  -1 d • F d (with F d (α) := µ((-∞, α]) and G d (α) := ν((-∞, α])). Then let T d-1 = T d-1 (x d-1 , x d ) be such that T d (., x d ) is monotone and maps µ d-1 x d to ν d-1

2 Figure 4 :

 24 Figure 4: Ten sample points: evolution of the tesselation for ε = 0 to ε = +∞ (from top left to bottom right).

  0

Figure 5 : 2 Figure 6 :

 526 Figure 5: Fifteen sample points. Top row: exact algorithm (gradient method). Middle row: continuation algorithm, 3000 steps. Bottom row: continuation algorithm, 10,000 steps. 25

  But here there is a function of the three variables (x d-1 , x d , y d ).Yet, we know that the measure γ d is concentrated on the set y d = T d (x d ) for a certain map T d , and this allows to replace the expression of y d , thus getting rid of one variable. This proves that the function ψ(x d , y d )φ(x d-1

	1
	K
	have the same marginals with respect to the pairs (x d-1 , x d ), (y d-1 , y d ) (since
	they have the same projections onto x and onto y) and (x d , y d ) (since we just
	proved it).

  Figure 2: Five sample points: evolution of the tesselation for ε = 0 to ε = +∞ (from top left to bottom right).

		1 1			1 1												1 1						0.335 0.355						1 1
					0.9 0.9												0.9 0.9												0.9 0.9
		0.9 0.9			0.8 0.8												0.8 0.8												0.8 0.8
																							0.354
					0.7 0.7												0.7 0.7												0.7 0.7
		0.8 0.8
					0.6 0.6												0.6 0.6												0.6 0.6
		0.7 0.7		x 2 x 2	0.5 0.5											x 2 x 2	0.5 0.5						0.353					x 2 x 2	0.5 0.5
					0.4 0.4												0.4 0.4						0.33						0.4 0.4
		0.6 0.6			0.3 0.3												0.3 0.3						0.352						0.3 0.3
					0.2 0.2												0.2 0.2												0.2 0.2
	x 2 x 2	0.5 0.5			0.1 0.1												0.1 0.1						0.351						0.1 0.1
		0.4 0.4			0 0	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 x 1 0.5 x 1	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1	0 0	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 x 1 0.5 x 1	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1	0 0	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 x 1 0.5 x 1	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1
					1 1												1 1						0.35						1 1
		0.3 0.3			0.9 0.9												0.9 0.9						0.325						0.9 0.9
		0.2 0.2			0.8 0.8												0.8 0.8						0.349						0.8 0.8
					0.7 0.7												0.7 0.7												0.7 0.7
		0.1 0.1			0.6 0.6												0.6 0.6						0.348						0.6 0.6
				x 2 x 2	0.5 0.5											x 2 x 2	0.5 0.5											x 2 x 2	0.5 0.5
		0 0	0 0		0.4 0.4	0.1 0.1		0.2 0.2		0.3 0.3		0.4 0.4	0.5 0.5		0.6 0.6	0.7 0.7	0.4 0.4		0.8 0.8	0.9 0.9		1 1							0.4 0.4
					0.3 0.3								x 1 x 1				0.3 0.3						0.15 0.32 0.315 0.347	0.155	0.16 0.32	0.165 0.3 0.325 0.3	0.17		0.175 0.33		0.18	0.185 0.335		0.19 0.34	0.195	0.2 0.345
		1 1			0.2 0.2												0.2 0.2						0.335 0.356						0.2 0.2
		0.9 0.9			0 0.1 0 0.1	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1	0 0.1 0 0.1	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5 0.355	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1	0 0.1 0 0.1	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1
		0.8 0.8			1						x 1 x 1						1						x 1 0.354 x 1						1						x 1 x 1
					0.9												0.9												0.9
		0.7 0.7
					0.8												0.8						0.33 0.353						0.8
		0.6 0.6			0.7												0.7												0.7
	x 2 x 2	0.5 0.5		x 2	0.5 0.6											x 2	0.5 0.6						0.352					x 2	0.5 0.6
		0.4 0.4			0.4												0.4						0.351						0.4
					0.3												0.3												0.3
		0.3 0.3			0.2 0.2												0.2						0.325 0.35						0.2
					0.1 0.1												0.1												0.1
		0.2 0.2			0 0	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1 1	0	0	0.1	0.2	0.3	0.4	0.5 0.349	0.6	0.7	0.8	0.9	1	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
											x 1												x 1												x 1
		0.1 0.1			1												1												1
					0.9												0.9						0.348						0.9
		0 0	0 0		0.8	0.1 0.1		0.2 0.2		0.3 0.3		0.4 0.4	0.5 0.5		0.6 0.6	0.7 0.7	0.8		0.8 0.8	0.9 0.9		1 1							0.8
					0.7								x 1 x 1				0.7						0.15 0.32 0.25 0.347	0.155 0.26	0.16 0.27	0.165 0.7 0.28	0.17 0.29		0.175 0.3		0.18 0.31	0.185 0.32		0.19 0.33	0.195 0.34	0.2 0.35
		1 1			0.6												0.6						0.335 0.356						0.6
				x 2	0.5											x 2	0.5											x 2	0.5
		0.9 0.9			0.4												0.4						0.355						0.4
		0.8 0.8			0.3												0.3												0.3
		0.7 0.7			0.1 0.2		steps Relative errors in cell areas 0.1 0.2 0.354	0.1 0.2								Time
	x 2 x 2	0.4 0.4 0.5 0.5 0.6 0.6		3,000 10,000 3,000 0.7 0.8 0.9 1 0 0.1 0.2 0	0.3	0.4	0.5 x 1	-2.24% -0.12% -0.58% 0.20% -3.20% 0.17% -2.66% 4.25% -2.43% 28.85% 0.7 0.7 0.8 0.8 0.351 (continued) 0.9 1 0.9 1 0.352 -0.67% -0.04% -0.17% 0.06% -0.97% ≤ 4 h 2 h 2 h 0.33 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 x 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x 1 0 0.353	1
		0.2 0.2 0.3 0.3		10,000 0.4 0.5 x 2 0.6				0.05% (continued) 0.4 0.5 x 2 -0.82% 1.28% -0.74% 9.42% ≤ 4 h 0.4 0.5 x 2 0.349 0.325 0.6 0.6 0.35
					0.3												0.3												0.3
		0.1 0.1		3,000 0.2				-30.11% 2.51% 1.82% 2.92% 0.63% 0.2 0.2 0.348			2 h
		0 0	0 0	10,000 0.1 0.2 0.1 0.1 0.2 0 0.1 0.2 0	0.3	0.3 0.3 0.4	0.5 x 1	0.4 0.4 -9.79% 0.76% 0.55% 0.89% 0.19% ≤ 4 h 0.5 0.6 0.7 0.8 0.9 1 0.1 0.1 0.5 0.6 0.7 0.8 0.9 1 x 1 0.15 0.155 0.16 0.165 0.17 0.175 0.18 0.185 0.19 0.195 0.32 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 x 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 x 1 x 1 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.347	0.2 0.35
	Figure 1: Five sample points. Top row: exact algorithm (gradient method). Figure 3: Ten sample points. Top row: exact algorithm (gradient method).
	Middle row: continuation algorithm, 100 steps. Bottom row: continuation Middle row: continuation algorithm, 2000 steps. Bottom row: continuation
	algorithm, 500 steps. algorithm, 3000 steps.									21 23

Acknowledgements G.C. and F.S. gratefully acknowledge the support of the Agence Nationale de la Recherche via the research project OTARIE. A.G. gratefully acknowledge the support of Chaire EDF-Calyon "Finance et Développement Durable," of Chaire AXA "Assurance et risques majeurs," and of Chaire Société Générale "Risques Financiers". The authors wish to warmly thank Yann Brenier and Alessio Figalli for stimulating discussions.