Convergence rates for discretized Monge-Amp\`ere equations and quantitative stability of Optimal Transport

In recent works - both experimental and theoretical - it has been shown how to use computational geometry to efficently construct approximations to the optimal transport map between two given probability measures on Euclidean space, by discretizing one of the measures. Here we provide a quantative convergence analysis for the solutions of the corresponding discretized Monge-Amp\`ere equations. This yields L^{2}-converge rates, in terms of the corresponding spatial resolution h, of the discrete approximations of the optimal transport map, when the source measure is discretized and the target measure has bounded convex support. Periodic variants of the results are also established. The proofs are based on quantitative stability results for optimal transport maps, shown using complex geometry.


Introduction
The theory of Optimal Transports [39], which was originally motivated by applications to logistics and economics, has generated a multitude of applications ranging from meteorology and cosmology to image processing and computer graphics in more recent years [37,38]. This has led to a rapidly expanding literature on numerical methods to construct optimal transport maps, using an appropriate discretization scheme. From the PDE point of view this amounts to studying discretizations of the second boundary value problem for the Monge-Ampère operator. The present paper is concerned with a particular discretization scheme, known as semi-discrete Optimal Transport in the Optimal Transport literature (see [6] and references therein for other discretization schemes, based on finite differences). This approach uses computational geometry to compute a solution to the corresponding discretized Monge-Ampère equation and exhibits remarkable numerical performance, using a damped Newton iteration [32,31]. The convergence of the iteration towards the discrete solution was recently settled in [28] and one of the main aims of the present paper is to establish quantitative convergence rates of the discrete solutions, as the spatial resolution h tends to zero.
1.1. Optimal transport maps and the second boundary value problem for the Monge-Ampère operator. We start with the general setup. Let X and Y be open bounded domains in R n with Y assumed convex, µ a probability measure supported on X and ν a probability measure with support Y, which will be assumed to have a uniformly positive density: We recall that in the case when µ also admits a density then a map T in L ∞ (X, Y ) is said to be a transport map (with source µ and target ν) if and T is said to be an optimal transport map (wrt the Euclidean cost function |x − y| 2 ) if it minimizes the L 2 −distance to the identity operator I over all transport maps. By Brenier's 1 theorem [13] there exists a unique optimal transport map T and it has the characteristic property of being a gradient map: T = T φ := ∇φ (in the almost everywhere sense) for a convex function φ on X, called the potential of T (see [39] for further background on Optimal Transport). The potential φ is the unique (modulo an additive constant) convex solution to the corresponding second boundary value problem for the Monge-Ampère operator, i.e. the sub-gradient image of φ is contained in Y, and φ solves the equation where the Monge-Ampère measure M A g is defined by when φ is C 2 −smooth and the general definition, due to Alexandrov, is recalled in Section 2.2. We will say that (X, Y, µ, ν) is regular if the corresponding solution φ is in C 2 (X). By Cafferelli's regularity results [16,17,18] this is the case if X and Y are assumed strictly convex with C 2 −boundary and the densities of µ and ν are Hölder continuous onX andȲ , respectively. Then φ defines a classical solution of the corresponding PDE and the corresponding optimal transport map T φ yields a diffeomorphism between the closures of X and Y. In fact, as is well-known, for any probability measure µ there exists a solution (in the weak sense of Alexandrov) of the corresponding second boundary value problem, which is uniquely determined up to normalization. In the sequel it will be convenient to use the normalization condition that the integral of a solution over (X, dx) vanishes. A time-honored approach for discretizing Monge-Ampère equations, which goes back to the classical work of Alexandrov on Minkowski type problems for convex bodies and polyhedra, amounts to replacing the given probability measure µ with a sequence of discrete measures converging weakly towards µ (see [1,Thm 7.3.2 and Section 7.6.2] and [5,Section 17]). A standard way to obtain such a sequence is to first discretize X by fixing a sequence of "point clouds" (x 1 , ..., x N ) ∈ X N and a dual tessellation of X with N cells (C i ) N i=1 , i.e. x i ∈ C i and the intersection of different cells have zero Lebesgue measure . For example, given a point cloud the corresponding Voronoi tessellation of X provides a canonical dual tessellation of X. The "spatial resolution" of the discretization is quantified by where diam(C i ), denotes the diameter of the cell C i . The corresponding discretization of the measure µ is then defined by setting (1.4) where we have used the subindex h to emphasize that we are focusing on the limit when h → 0 (see also Section 3.4 for other discretizations). This discretization scheme corresponds, from the point of view of Optimal Transport, to the notion of semi-discrete Optimal Transport (since it corresponds to Optimal Transport between a continuous and a discrete measure; see [28,31] and references therein). From the point of view of numerics this kind of discretization scheme was first introduced in the different setting of the Dirichlet problem in [35].
Given a point-cloud on X with spatial resolution h we denote by φ h the normalized convex solution to the corresponding Monge-Ampère equation 1.3 with right-hand side given by the discrete measure µ h . It follows from a standard general convexity argument that, in the limit when h → 0 (and hence N → ∞) the functions φ h converge uniformly towards the solution φ and the gradients of φ h converge point-wise almost everywhere to the gradient of φ. But the general argument gives no control on the rate of convergence (in terms of h or N −1 ) and the main purpose of the present note is to provide such a result: Theorem 1.1. (Regular case) Assume that (X, Y, µ, ν) is regular and let µ h be a discretization of µ. Denote by φ and φ h the solutions to the corresponding second boundary value problem for the Monge-Ampère operator. There exists a constant C (depending on (X, Y, µ, ν)) such that As a consequence, if C P denotes the constant in the L 2 −Poincaré inequality on X, For a general bounded convex domain Y and for any domain X the estimates above still hold if X replaced by any given compact subset (if the densities of µ and ν are strictly positive and Hölder continuous).
From the point of view of Optimal Transport theory the previous theorem says that the optimal transport map T φ (defining a diffeomorphism between the closures of X and Y ) may be quantitatively approximated by the L ∞ −maps T h defined by the gradients of φ h , in the sense that As explained in Section 5 the maps T h are piecewise constant on the convex hull of the corresponding point cloud. We will also establish the following universal bound which applies in the general setting and which, in particular, yields a quantitative approximation of the corresponding optimal transport map if µ has a density. Theorem 1.2. (General case) Let X and Y be bounded domains in R n with Y assumed convex and ν a probability measure of the form 1.1. Then, for any given probability measure µ on X, where the constant C only depends on upper bounds on the diameters of X of Y and the lower bound δ of ν/dy, if the volumes of X and Y have been normalized.
By the Poincaré inequality on convex domains [36] such an inequality also holds for the The key new input in the proofs of the theorems above are estimates wich can be viewed as quantitative stability results for optimal transport maps (see Section 4).
1.2. Computional geometry. An important motivation for the present work comes from the recent result in [28], showing that the vector φ h := (φ h (x 1 ), ...φ h (x N )) ∈ R N -which solves a discrete variant of the Monge-Ampère equation (with target Y ) -may be effectively computed using a damped Newton iteration on R N , which converges globally at a linear rate towards φ h (and locally at quadratic rate if g is Lipschitz continuous). The iteration is defined in terms of computational geometry and the restriction of the solution φ h to the convex hull 3 X h of the points {x 1 , ..., x N } can then be recovered as the piecewise affine function defined by convex hull of the discrete graph of φ h . From this computational point of view Theorems 1.1, 1.2 above yield quantitative convergence results for the corresponding discrete objects defined on X h in the "continuous limit" when h → 0. This is explained in Section 5.
1.3. The periodic setting. Now assume that µ is a given Z n −periodic measure on R n normalized so that its total mass on a (or equivalently any) fundamental region X is equal to one (e.g. X = [0, 1[ n ). We then consider the corresponding Monge-Ampère equation 1.3 for a convex function φ on R n subject to the condition that ∂φ be periodic (which replaces the second boundary condition 1.2). Such a convex function will be called quasi-periodic.
In this periodic setting Theorem 1.1 still applies (with X = [0, 1[ n ) (see Section 6). In terms of Optimal Transport the induced diffeomorphism T φ of the torus (R/Z) n transporting µ to the Riemannian volume form on the flat torus is optimal wrt the cost function d(x, y) 2 where d denotes the Riemannian distance function on the flat torus [15,27]. The interpretation 1.7 still applies in the periodic setting if d L 2 is taken as the induced L 2 −metric on L ∞ −self maps of the torus.
1.4. Comparison with previous results. There seem to be no prior results giving convergence rates (in any norms) for the functions φ h or the vectors φ h in the limit when h → 0 (even in the model case of a uniform grid on the torus). Another quantitative stability result has previously been obtained by Ambrosio (reported in [23]) and it can be used to obtain rates of convergence for the corresponding approximations of the inverse optimal transport map from ν to µ. The corresponding approximations are the optimal maps transporting ν to µ h . This corresponds to an H 1 −estimate for the difference of Legendre transforms φ * h − φ * (see Section 4). In the different setting of the Oliker-Prussner discretization of the Dirichlet problem for the Monge-Ampère equation [35] with g = 1 (where the second boundary condition 1.2 is replaced by the vanishing of the solution φ at the boundary of X) convergence rates for H 1 −norms (and more generally, W 2 p −norms) have recently been obtained in [33]. It is assumed in [33] that the point-cloud is a uniform grid and then an H 1 −estimate of the form 1.5 is established, but with Ch replaced by Ch log(h −1 ), where C depends on the higher derivative norm φ C 3,1 (X) (which if finite under the assumption that f ∈ C 1,1 (X)). The results in [33] build on convergence rates for L ∞ −norms established in [34], which are still missing in the present setting of the second boundary value problem (and in the periodic setting considered in Section 6).
1.5. Proof by complexification. The proof of the H 1 −estimates in theorems above use a complexification argument to deduce the estimates from well-known inequalities in Kähler geometry and pluripotential theory (due to Aubin [2] and Blocki [11]). The universal dependence on the diameters in Theorem 1.2 is obtained by exploiting that a convex body Y induces a canonical toric Kähler-Einstein metric, whose analytical properties are controlled by the diameter of Y (thanks to the estimates in [22,29]). The proofs are particularly simple in the periodic setting as there are no boundary terms that need to be taken into account. This approach is, in spirit, similar to Gromov's approach [25] to the Brunn-Minkowski and Alexander-Fenchel inequalities for convex bodies, exploiting the complex geometry of toric varieties.
1.6. Organization. We start in Section 2 with preliminaries from convex and complex analysis. Since the complex analytic side may not be familiar to some readers a rather thorough presentation is provided. In Section 3 the theorems stated in the introduction are proven, starting with the special case when the density g is constant on Y and finally reducing to the special case. In Section 4 the relations to quantitative stability of Optimal Transport maps are spelled out. The relations to computational geometry are explained in Section 5.4, based on Prop 5.1. In the final section the periodic setting is considered. 1.7. Acknowledgments. I am grateful to Quantin Mérigot for illuminating comments and, in particular, for drawing my attention to the paper [23]. This work was supported by grants from the ERC and the KAW foundation.

Preliminaries
2.1. Notation. Given a convex function φ on R n taking values in ] − ∞, ∞] we will denote by ∂φ its sub-gradient, by φ * its Legendre transform and by M A(φ) its Monge-Ampère measure (these notations are recalled below). Given a subset K ⋐ R n we will denote by χ K its indicator function defined as 0 on K and ∞ on its complement. We will also use some standard complex analytic notions recalled below.
2.2. Convex analytic notions. Given a convex function φ on R n taking values in ] − ∞, ∞] we denote by ∂φ its subgradient, i.e. the set-valued function on R n defined by y · x − φ(x), the following formula holds: i.e. M A(φ) is the measure obtained as the push-forward of the Lebesgue measure dy under the L ∞ loc −map ∇φ * (the formula follows from point 2 and 3 below). More generally, given g ∈ L 1 (R n ) the measure M A g (φ) is defined by replacing dy in formula 2.1 (or, equivalently, in formula 2.2) with gdy.
We recall the following basic properties which hold for a given lower semi-continuous (lsc) convex function φ : R n →] − ∞, ∞]: If y ∈ (∂φ)(R n ), i.e. y ∈ ∂φ(x y ) for some x y ∈ R n , then φ * (y) = sup The function φ is piecewise affine on R n iff it is the max of a finite number of affine functions, i.e. there exists y 1 , ..., y M ∈ R n (assumed distinct) and c 1 , ..., c M in R such that φ(x) = max Denote by F i the closure of a maximal open region where the piecewise affine function φ is differentiable (or equivalently, affine). Then we can label F i such that ∇φ(x) = y i on the interior of F i . Moreover, the covering and ∂φ(x i ) is the convex hull of the vectors y j associated to all facets F j containing x i . A reference for point 1-5 is [39] and 6 could be taken as the definition of a convex piecewise affine function and then point 7 follows readily. As for point 8 it can be shown using the following observation: x is not in a 0−dimensional cell iff the convex hull C x (= ∂φ(x)) of the vectors y i corresponding to the facets F i containing x has dimension p < n (as can be seen by identifying the F i s with intersecting pieces of hyperplanes in the graph of φ in R n+1 and C x with the corresponding normal cone at (x, φ(x)). As a consequence, if x is not in a 0−dimensional cell, then there exists a whole neighbourhood U of x having the latter property and hence (∂φ)(U ) is a null-set for Lebesgue measure, i.e. M A(φ) = 0 in U. This shows that the support of M A(φ) is contained in {x 1 , ..., x n }. Conversely, if x is contained in the latter set, then p = 0 and hence (∂φ)(x) has dimension n, i.e. M A(φ){x} = 0, 2.2.1. The class C Y (R n ) of convex functions associated to a bounded convex domain Y.. Given a bounded convex domain Y, normalized to have unit volume, we denote by C Y (R n ) the space of all convex functions on R n such that (∂φ)(R n ) ⊂ Y . A reference element in C Y (R n ) is provided by the support function of Y : A leading role in the present paper will be played by the subspace C Y (R n ) + of C Y (R n ) consisting of the convex functions φ on R n with "maximal growth" in the sense that the reversed version of the previous inequality also holds: is a probability measure (by point 3 in the previous section). The converse is not true in general, but we will have great use the fact that if M A(φ) is a probability measure and moreover M A(φ) has compact support, then φ is in C Y (R n ) + , as follows from the following lemma this follows from the Sobolev inequality on Y (see [7,Prop 2.2]). Since Y is convex the corresponding Sobolev constant is dominated by diam (Y ) q(ǫ) C n,ǫ , where C n,ǫ only depends on n and ǫ and q(ǫ) → n as ǫ → 0 (by [20, Lemma 1.7.3]). q.e.d.

Complex analytic notions.
For the benefit of the reader lacking background in complex analysis and geometry we provide a (hopefully user friendly) recap of some complex analytic notions (see the book [21] for further general background and [7] for the case of toric varieties). Setting z := x + iy ∈ C n the space Ω 1 (C n ) of all complex one-forms on C n decomposes as a sum of the two subspaces spanned by {dz i } and {dz i }, respectively. This induces a decomposition of the exterior algebra of all complex differential forms Ω · (C n ) into forms of bidegree (p, q), where p ≤ n and q ≤ n. Accordingly, the exterior derivative d decomposes as d = ∂ +∂, where and taking its complex conjugate defines the (0, 1)−form∂φ. In particular, defines a real (1, 1)−form (the normalization above turns out to be convenient). Such a smooth form ω is said to be positive (Kähler), written as ω ≥ 0 (ω > 0), if the corresponding Hermitian matrix is semi-positive (positive definite) at any point. Positivity can also be defined for general (p, p)−forms, but for the purpose of the present paper it is enough to know that where the last inequality holds in the sense of measures.
loc is said to be psh if it is strongly upper semi-continuous and ω φ ≥ 0 holds in the weak sense of currents. More generally, given a real (1, 1)−form ω 0 a function u said to be ω 0 −psh if is defined by the local pluripotential theory of Bedford-Taylor. The current does not charge pluripolar subsets (i.e sets locally contained in the −∞−locus of a psh function) and in particular not analytic subvarieties. Accordingly, the Monge-Ampère measure of a locally bounded psh function φ(z), is defined by the measure We also recall that, if φ 0 and φ 1 are as above then the positive current i∂(φ 0 − φ 1 ) ∧∂(φ 0 − φ 1 ) may be defined by the formula Metrics on line bundles and ω 0 −psh functions. The local complex analytic notions above naturally extend to the global setting where C n is replaced by a complex manifold (since the decomposition 2.4 is invariant under a holomorphic change of coordinates). However, if X is compact, then all psh functions φ on X are constant (by the maximum principle). Instead the role of a (say, smooth) psh function φ on C n is played by a positively curved metric on a line bundle L → X (using additive notations for metrics). To briefly explain this first recall that a line bundle L over a complex manifold X is, by definition, a complex manifold (called the total space of L) with a surjective holomorphic map π to X such that the fibers L x of π are one-dimensional complex vector spaces and such that π is locally trivial. In other words, any point x ∈ X admits a neighborhood U such that π : L → U is (equivariantly) isomorphic to the trivial projection U × C → U. Fixing such an isomorphism holomorphic sections of L → U may be identified with holomorphic functions on U. In particular, the function 1 over U corresponds to a non-vanishing holomorphic section s U of L → U. Now, a smooth (Hermitian) metric · on the line bundle L is, by definition, a smooth family of Hermitian metrics on the one-dimensional complex subspaces L x , i.e. a one-homogeneous function on the total space of the dual line bundle L * which vanishes precisely on the zero-section. Given a covering U i of X and trivializations of L → U i a metric · on L may be represented by the following family of local functions φ U i on U i : (accordingly a metric on L is often, in additive notation, denoted by the symbol φ). Even if the functions φ U i do not agree on overlaps, the (normalized) curvature form ω of the metric · is a well-defined closed two-form on X, locally defined by Singular metrics on L may be defined in a similar way. In particular, a singular metrics is said to have positive curvature if the local functions φ U i are psh, i.e. the corresponding curvature form ω defines a positive (1, 1)−current on X. The difference of two metrics, written as φ 1 − φ 2 in additive notation, is always a globally well-defined function on X (as a consequence, the curature currents of any two metrics on L are cohomologous and represent the first Chern class c 1 (L) ∈ H 2 (X, Z)). Fixing a reference metric φ 0 and setting u := φ − φ 0 this means that the space all metrics φ on L with positive curvature current may be identified with the space of all ω 0 −psh functions u on X.
comes with a tautological line bundle whose total space is C m+1 and the line over a point The dual of the tautological line bundle is called the hyperplane line bundle and is usually denoted by O(1) (the notation reflects the fact that the metrics on O(1) may be identified with 1−homogeneous functions on C m+1 ). The Euclidean metric on C m+1 induces a metric on the tautological line bundle and hence on its dual O(1), called the Fubini-Study metric. In the standard affine chart U := C n ⊂ P n defined by all points where z 0 = 0 and with the standard trivializing section s U := (1, z 1 , ..., z n ) of the tautological line bundle the Fubini-Study metric is represented by defining a smooth metric with strictly positive curvature form. Another T n −invariant metric on O(1) with positive curvature current, which is continuous, but not smooth, is the "max-metric" defined by the one-homogeneous psh function max{|Z 0 |, ..., |Z n |}, which may be 8 represented by φ max (z) = log max{1, |z 1 | 2 , ..., |z n | 2 }, Given any complex subvariety X ⋐ P m one obtains a line bundle L over X by restricting O(1) to X (and a smooth positively curved metric φ by restricting the Fubini-Study metric). If X is singular then, by Hironaka's theorem it admits a smooth resolution, i.e. a smooth compact manifold X ′ with surjective and generically one-to-one projection π ′ to X. Pulling back L by π ′ yields a line bundle L ′ over X ′ (endowed with a smooth metric φ ′ ).
the standard Kähler form on C n . The form ω 0 descends to the Abelian variety C n /(Z + iZ).

2.4.2.
The complex torus C * n . Let Log be the map from C * n to R n defined by The real torus T n acts transitively on the fibers of the map Log. Pulling back a convex function φ(x) on R n by Log yields a T n −invariant function on C * n that we will, abusing notation slightly, denote by φ(z). The function φ(x) is convex iff the corresponding function φ(z) on C * n is plurisubharmonic. This can be seen by identifying C * n with C n /iπZ n (= R n + iπT n ) using the multivalued logarithmic coordinate 2 log z, and proceeding as in Section 2.4.1. Applying formula 2.6 also gives As recalled in the following section a given rational polytope Y in R n (with non-empty interior) determines a toric variety compactifying the complex torus C * n .

2.4.3.
The toric variety associated to a moment polytope Y . Let Y be a bounded closed convex polytope with non-empty interior and rational vertices. It determines a compact toric complex analytic variety X C and an ample line bundle L over X C . More precisely, this is the case if the rational polytope Y is replaced by the integer polytope kY for a sufficiently large positive integer k. But since scalings of Y will be harmless we may as well assume that k = 1. Then X C may be defined as the closure of the image of the holomorphic (algebraic) embedding defined by By construction, the following holds: • X C may be embedded in a complex projective space P N in such a way that L coincides with the restriction to X C of the hyperplane line bundle O(1) on P N . • The standard action of the real n−dimensional torus T n on C * n extends to a holomorphic action of X C which lifts to the line bundle L.
By the first point above we can identify C * n with an open dense subset of X C , whose complement in X C is an analytic subvariety. Now, the key point is that the function φ(x) is in the class C Y (R n ) (Section 2.2.1) iff φ(x) extends to a positively curved (singular) metric on L → X C . More precisely, we have the following [7, Prop 3.2] Lemma 2.4. A convex function φ(x) such that C Y (R n ) + may be identified with a L ∞ −metric on the line bundle L → X C with positive curvature current ω φ .
Proof. Since the lemma will play a key role in the proof of the main results we recall, for the benefit of the reader, the simple proof. Without loss of generality we can assume that m 0 := 0 is in Y. By construction the "max metric" on O(1) → P M −1 restricts to a continuous (and in particularly bounded) positively curved metric on L → X C . The map 2.9 may be factored as C * n → C N → P M −1 . Hence, the standard trivialization of O(1) over the affine piece C N pulls back to give a trivialization of L over C * n , where the restricted max metric is represented by φ Y (x), when expressed in the logarithmic coordinates x on R n . Now, any other L ∞ −metric φ on L → X satisfies φ − φ Y ∈ L ∞ (X C ) and, as a consequence, restricting to C * n and switching to the real logarithmic coordinate x shows that φ( Since C * n is dense in X C and u ≤ C it admits a canonical upper semi-continuous (usc) extension to all of X C (namely, the smallest one). Since φ Y (z) extends to define a continuous metric on L → X C this means that φ extends from C * n to a canonical usc metric on L → X C (in additive notation) which has a positive curvature current on C * n ⋐ X C . But the complement X C − C * n is an analytic subvariety of X C and hence it follows from basic local extension properties of psh functions that the corresponding metric on L → X C has positive curvature current. q.e.d.
We note that for any φ as above (2.10) 1 n! X C (ω φ 1 ) n =: Indeed, since M A C (φ) does not charge analytic subvarieties the integral can be restricted to C * n and then formula 2.8 can be invoked.
3. Proof of Theorems 1.1, 1.2 We will reduce the proof of the theorems stated in the introduction to the case when g = 1. Given a function φ which is smooth and strictly convex we will denote by g φ the corresponding Hessian metric, i.e. the Riemannian metric defined by the symmetric Hessian matrix ∇ 2 φ. Lemma 3.1. Let X and Y be open convex domains in R n , with Y bounded and normalized to have unit volume. Given a probability measure µ with compact support contained in X, a solution φ X to the corresponding second boundary value problem 1.3, 1.2 exists and is uniquely determined (mod R). Moreover, φ X is equal to the restriction to X of φ R n (mod R).
Proof. This result goes back to Alexandrov's work on Monge-Ampère equations, but for completeness we show here how to deduce the result from Brenier's theorem [13]. In the proof we will refer to point 1 − 8 in Section 2.2. Set ψ X := (χ X φ X ) * (where we have used that φ X extends uniquely to a Lipschitz continuous function on the closure of X). Since µ and ν have the same mass the second boundary condition 1.2 implies that ∂φ X (R n ) =Ȳ. Hence ψ X is finite on Y (by point 3) and the MA-equation 1.3 is equivalent to By Brenier's uniqueness theorem [13] the latter equation determines the L ∞ −map ∇ψ X a.e ν.
Since the support Y of ν is connected the restriction ψ of ψ X to Y is thus uniquely determined (mod R). Now, by point 1 χ X φ X = ψ * X and since ∂(χ X φ X )(R n ) ⋐Ȳ we get χ X φ X = (χ Y ψ X ) * (by point 5). In particular, on X we have φ X = (χ Y ψ) * , where ψ, as explained above, is independent of X. Hence, replacing X with R n reveals that the corresponding two solutions coincide on X (mod R) as desired. Also note that Brenier's existence theorem says that, given µ and ν, there exists some convex function ψ on Y satisfying the equation 3.1. Hence defining φ X in terms of ψ as above yields the existence of a solution to the second boundary value problem. q.e.d.
Remark 3.2. If X is not convex, then a solution need not be unique (mod R), but a canonical solution (mod R) is always provided by φ R n .
The following result gives the key inequality that we shall use.
Given bounded open domains X and Y with Y assumed convex, assume that φ 0 and φ 1 are convex functions on X, such that the the closures of the sub-gradient images (∂φ i )(X) are equal to Y . If φ ∈ C 2 (X) and there exists a positive constant C such that Proof. It will be enough to prove the following inequality: where c n is a positive constant only depending on n. In essence, the idea of the proof is to use integration by parts. But to handle boundary terms we will use a three step process of extension, complexification and compactification. We will denote by φ 0 and φ 1 the canonical extensions to R n solving the corresponding second boundary problem on all of R n (as in the previous lemma). In general, if φ(x) is a convex function on R n we will denote by φ(z) the corresponding T n −invariant plurisubharmonic function on C * n , obtained by pulling back φ to C * n using the Log map (as in Section 2.4.2). The starting point of the proof is the following basic formula, where φ(x) is strictly convex and C 2 −smooth on an open subset X of R n : for any two convex smooth functions φ 0 and φ 1 on X, where c n is a positive constant only on the dimension n. This is proved using a point-wise calculation, as in the proof of formula 2.8. Hence, the inequality in the proposition would follow from the following identity: since all the integrands in the rhs above are non-negative. Note that as follows directly from the algebraic identity 11 (3.4) This means that if integration by parts is justified, then the desired identity 3.3 would follow. However, the non-compactness of C * n poses non-trivial difficulties, so instead of working directly on C * n we will use a compactification argument, which applies when Y is a rational polytope (the general case then follows by approximation).
Compactification when Y is a convex polytope with rational vertices Let X C be the compact toric complex analytic variety determined by the moment polytopē Y and denote by L the corresponding ample line bundle over X C (Section 2.4.3). By Lemma 2.12.4 combined with Lemma 2.4 the T n −invariant psh functions φ 0 (z) and φ 1 (z) on C * n extend to define L ∞ −metrics on the line bundle L → X C with positive curvature currents ω φ 0 and ω φ 1 . In particular, φ 0 − φ 1 ∈ L ∞ (X C ). We claim that q.e.d.
This is a well-known identity in Kähler geometry (in the smooth case) and global pluripotential theory (in the general singular L ∞ −case) and follows from the general integration by parts formula for psh functions in L ∞ loc in [12, Thm 1.14] (see [9, Formula 2.9]). But for the benefit of the reader we provide the following alternative proof. First, assume that the metrics φ i on L are smooth (i.e. the restrictions to L → X of smooth metrics on O(1) → P N ). Expanding point-wise, as in formula 3.4 and using that the form Θ := (ω φ 0 ) n−j ∧ (ω φ 1 ) j is closed formula 3.5 then follows from Stokes formula if X C is non-singular and from Stokes formula on a nonsingular resolution of X C if X C is singular. To handle the general case we invoke the general fact that any (possibly singular) metric on an ample line bundle L over a projective complex variety can be written as a decreasing limit of smooth metrics [14] (in fact, this can be shown by a simple direct argument in the present toric setting). Formula 3.5 then follows from the previous smooth case, combined with the continuity of expressions of the form appearing in formula 3.5 under decreasing limits (see [12,Prop 2.8] and its proof for much more general convergence results).

Conclusion of proof
Let Y be as in the previous step. Since the complex Monge-Ampère measures of locally bounded psh functions do not charge complex analytic subvarietes formula 3.5 on C * n follows from formula 3.5 on X C . Since each term in the right hand side above is non-negative we deduce that By 2.8 and 3.2 this proves formula 3.5 when Y is a convex polytope with rational vertices. In the general case, we can write Y is an increasing limit of rational convex polytopes Y j . Then the general case follows from the previous case and basic stability properties for the solution of the second boundary value problem for M A on X wrt variations of the target domain Y.

3.1.
Proof of Theorem 1.1 when g = 1. First observe that the previous proposition implies that where diam(Y ) denotes the diameter of Y and W 1 is the Wasserstein L 1 −distance on X. Indeed, this follows directly from the Kantorovich-Rubinstein formula where the sup runs over all Lipschitz continuous functions on X with Lipschitz constant one (by taking u = −φ 0 /diam(Y ) and u = φ 1 /diam(Y )). Applying the inequality 3.6 to φ 0 = φ and φ h all that remains is thus to verify that To this end we rewrite Since u ∈ Lip 1 (X)we trivially have where the last inequality follows directly from the very definition of h. Hence, W 1 (µ 0 , µ 1 ) is bounded from above by h times the total mass of µ, which concludes the proof of the H 1 −estimate in Theorem 1.1 (when g = 1). The L 2 −Poincaré inequality then implies the corresponding L 2 −estimate. The last statement in Theorem 1.1 is shown by replacing X with a compact subset and using that the interior regularity of φ holds without any convexity assumption on X or smoothness assumptions on X and Y [17].

3.2.
Proof of Theorem 1.2 when g = 1. As in the proof of Theorem 1.1 it will be enough to prove the following Proposition 3.4. Given bounded open domains X and Y with Y assumed convex, assume that φ 0 and φ 1 are convex functions on X, such that the the closures of the sub-gradient images (∂φ i )(X) are equal to Y . Then there exists a constant C only depending on X and Y such that Proof. It will bee enough to show that To this end we fix a sequence of rational convex polytopes Y j decreasing to Y and some (say uniformly on R n . Just as in the proof of Theorem 1.1 it will be enough to prove that the inequality in the proposition holds when Y is a rational convex polytope and φ Y is replaced by φ Y j , if the corresponding constants C Y j are uniformly controlled by Y. To this end, we identify, just as in the proof of Theorem 1.1 a convex function φ(x) on R n with a psh function φ(z) on C * n . Setting ω Y j := ω φ Y j it will thus be enough to show that We will deduce this inequality from the following inequality of Blocki [11]. Let (X C , ω) be a compact Kähler manifold and u 0 and u 1 ω−psh functions on X in L ∞ (X). Then there exists a constant A such that The constant A only depends (in a continuous fashion) on upper bounds on the L ∞ −norms of u 0 and u 1 and the volume of ω. More generally, exactly the same proof as in [11] shows that the inequality holds more generally when ω is a semi-positive form with an L ∞ −potential and positive volume. In the present setting we set u 0 : originally defined on C * n . Just as in the proof of Theorem 1.1 the functions φ 0 , φ 1 and φ Y j may be identified with L ∞ −metrics on the ample line bundle L → X j over the toric variety X j determined by Y j . Moreover, the corresponding currents extend to positive currents on X j . Passing to a smooth resolution we may as well assume that X j is smooth (and L semi-ample). Applying the inequality 3.8 on X j all that remains is to verify that the corresponding constants A j are uniformly bounded in j. To this end first note that we may assume that φ 0 and φ 1 are normalized so that the integral of their Legendre transforms vanishes over Y which, by assumption is contained in Y j . Applying the Sobolev inequality on Y (see Lemma 2.1 and its proof) and using that, by assumption, Y is contained in Y j , gives In particular, if µ 0 := M A(φ 0 ) is supported in X then Hence, by the triangle inequality where, by assumption, the norm in the rhs above converges towards as j → ∞, where we have used that u 0 is unaltered when (φ 0 , φ Y ) are replaced by (φ 0 +c, φ Y +c) (which can be done before normalizing φ 0 as above). The same inequality also holds, for the same reasons, for u 1 . Moreover, since the volume of ω Y j on X C is equal to n! times the volume of Y j (formula 2.10), which converges to Y, this proves the uniformly in question. All in all this proves the proposition. q.e.d.

3.2.1.
The dependence on Y and X. The constant C appearing in the previous proposition can be made to only depend on upper bounds on the diameters of X and Y, if the the volumes are normalized (the lower bounds are automatic, by the volume normalization). To see this first note that, by basic equivariance under translations, we may as well assume that 0 ∈ X and that 0 is the barycenter of 1 Y dy. It will then be enough to show that the function φ Y can be chosen so that the corresponding constant A 1 , A 2 and A 3 only depend on the outer radius R(X) and R(Y ) of X and Y (defined as in formula 3.9). To this end we take φ Y to be the unique smooth and strictly convex function is a Kähler potential for a Kähler-Einstein metric on C * n ). By [7] such a solution exists and is uniquely determined. Morever, by [29] the solution satisfies ∇ 2 φ Y ≤ 2R(Y ) 2 I. Thus all that remains is to verify the bounds on A 2 and A 3 . To this end we first recall the following a priori estimate in [22,Section 3.3 ] (note that in our setting the constant C defined by [22, formula 22] vanishes): for any p > 0 where C p,Y only depends on p and on R(Y ).
Step 2 in the proof of the previous estimate in [22, Section 3.3 ] establishes a bound of the form and, as a consequence, for any x in R n (as is seen by integrating ∇φ Y along the line segment between 0 and x). This shows that A 2 is bounded from below by e −C R(Y ) −R(X)R(Y ) . Finally, combining 3.10 with Lemma 2.1 shows that A 3 is also under control.

3.3.
Reduction to the case g = 1. Fix φ 0 ∈ C Y (R n ) + and set µ 0 := M A g (φ 0 ). Consider the following functionals on C Y (R n ) + : To simplify the notation we will write I 1 = I and J 1 = J. The reduction to the case g = 1 is accomplished by the following Proposition 3.5. Setting δ := inf Y g we have The rest of the section is devoted to the proof of the proposition. We start with a number of lemmas.
Lemma 3.6. The Gateaux differential of J g at φ is represented by −M A g (φ) + µ 0 , i.e. fixing φ 1 and φ 2 in C Y (R n ) + and setting φ t := φ 1 + t(φ 2 − φ 1 ) we have Proof. This is essentially well-known in the literature of optimal transport, where J g appears as the Kantorovich functional (up to a change of sign and and additive constant). A simple direct proof is given in [7]. q.e.d.
Remark 3.7. The formula for the differential of J g implies that when g = 1 and Y is a rational convex polytope, then J g coincides with Aubin's J−functional defined with respect to the curvature form ω 0 corresponding to φ 0 (see [2] for the smooth setting and [9] for the general singular setting).
Using the previous lemma we next establish the following Lemma 3.8. The following inequality holds: Proof. Fix φ in C Y (R n ) + and set φ t := φ 0 + t(φ − φ 0 ) and I g (t) := I g (φ t ) and J g (t) := J g (φ t ).
By the previous lemma we have, for any t > 0. Hence, using that J g (t) is convex, by the previous lemma. Since I g (0) = J g (0)(= 0) it follows that I g (t) ≥ J g (t) for all t ≥ 0 and in particular for t = 1, which concludes the proof. q.e.d.
We next show that J g is controlled from below by J : Lemma 3.9. Setting δ := inf Y g we have J g ≥ δJ ≥ 0 Proof. We continue with the notation in the proof of the previous lemma. By a standard approximation argument we may as well assume that φ and φ 0 are smooth and strictly convex. Moreover, as before it will be enough to consider the case when Y is a rational polytope. First observe that it will be enough to prove the following claim: Claim: Indeed, by Lemma 3.6 the derivatives of both J g (t) and J(t) vanish at t = 0. As a consequence, if the claim holds, then dt for all t and since J g (0) = J(0) = 0 it follows that J g (t) ≥ δJ(t) for all t and in particular for t = 1, proving the lemma. In order to prove the claim above it is enough to establish the following formula: To this end first consider a general curve φ t (x) in C Y (R n ) + such that Φ(x, t) := φ t (x) is smooth on R n+1 . We complexify R n by C * n (as in Section 2.4.2) and R by C (in the usual way). Denote by (z, τ ) the corresponding complex coordinates on C * n × C. Denoting by π the natural projection from C * n × C to C * n we have the following general formula on C. This formula is a special case of the more general formula [10, formula 2.11]. Next, we note that, if φ t is affine, then M A C (Φ) ≤ 0. Indeed, a direct calculation gives Using formula 3.2 this concludes the proof of formula 1.2. q.e.d.
3.3.1. Conclusion of proof of Prop 3.5. Combing the previous lemmas with Lemma 3.8 gives Hence the proof is concluded by noting that J ≥ I/(n+1). In the complex setting this is a well-known inequality in the Kähler geometry literature which goes back to Aubin [2]. For a simple proof see [9]. The present real setting then follows from a complexification and approximation argument, just as before.
3.4. Other discretization schemes when µ has a density. Assume for simplicity that the domain X has been normalized to have unit volume and fix a sequence of point clouds (x 1 , ..x N ) on X. In the case when µ has a density f it may, from a computational point of view, by more convenient to consider discretizations ofμ h of the form for appropriate weights w i (independent of µ). For example, a convenient choice is to take w i to be equal to the volume of the intersection with X of the Voronoi cell corresponding to x i . We claim that the result in Theorem 1.1 then still holds if f is Lipschitz continuous. Indeed, by the proof of Theorem 1.1, it is enough to show that which, in turn, will follow from W 1 (dx, λ N ) ≤ h. But the latter inequality follows from the identity 3.7 applied to µ = dx. Similarly, if f is merely assumed to be in the Hölder space C α (X) for α < 1 then the same argument reveals the an analog of Theorem 1.1 holds, with the rate Ch 1/2 replaced by Ch α/2 .
Finally, in the case when the point cloud {x i } is equal to the intersection of X with the grid (ZN −1/n ) n and the distance between the point cloud and the boundary of X is of the order O(N −1/n ), then a slight variant of the previous argument reveals that the same results as above hold, with h := N −1/n , when all weights w i in formula 3.12 are taken to be equal to 1/N.

Relations to quantitative stability of optimal transport maps
Denote by P ac (R n ) the space of all probability measures on R n , which are absolutely continuous with respect to Lebesgue measure. Given µ and ν in P ac (R n ) with compact support there exists optimal transport maps T ν µ and T µ ν in L ∞ (X, Y ) and L ∞ (Y, X) transporting µ to ν and ν to µ, respectively. The maps are uniquely determined and inverses of each others almost every where wrt µ and ν, respectively. Theorem 4.1. (Quantitative stability wrt the source) Assume that (X, Y, µ, ν) is regular. Let µ t be a curve in µ t ∈ P ac (X) converging weakly towards µ ∈ P ac (X), as t → 0, with µ strictly positive and Hölder continuous onX. Then there exists a constant C, depending on µ and ν such that for any t, where d L 2 denotes the Euclidean L 2 −metric on L ∞ (X, Y ) and W p denotes the Wasserstein L p −metric on P(X). As a consequence, for any p ∈ [1, ∞[ Moreover, if µ and µ t have densities f and f t in L 2 (X), then In the general case when X and Y are assumed to be bounded domains in R n with Y convex the first inequality above holds if W 1 (µ, µ t ) 1/2 is replaced by W 1 (µ, µ t ) 1/2 n (and similarly for W p (µ, µ t ) 1/2 ) and then the constant C only depends on upper bounds on the diameters of X of Y and the lower bound δ of ν/dy, if the volumes of X and Y have been normalized.
Proof. By Brenier's theorem we can write T ν µt = ∇φ t for a convex function φ t on X. The first inequality then follows directly from the inequality 3.6. To prove the last inequality first note that Proposition 3.3 implies that Applying, the Poincaré inequality to φ 0 − φ 1 (whose integral wrt (X, dx) may be assumed to vanish) then concludes the proof. q.e.d.
The previous result is reminiscent of a previous result due to Ambrosio (reported in [23]), which says that giving quantitative stability wrt variations of the target measure, rather then the source measure. In particular, We note that the latter estimate yields an analog of Theorem 1.1, formulated in terms of the corresponding Legendre transforms: Under the same assumptions as in Theorem 1.1 there exists a constant C such that Proof. First recall that if T ν µ = ∇φ then T µ ν = ∇φ * (as follows from point 2 in Section 2.2). Next, observe that, for a fixed t, the inequalities 4.1 and 4.2 actually hold more generally for any given µ t ∈ P(X). Indeed, let µ j t be a sequence in P(X) converging weakly towards µ t (and hence, also wrt the W p (X)−Wasserstein distance for any p). Then it follows from well-known (non-quantitative) stability arguments that T µ j ν (y) → T µt ν (y) a.e. on X. Hence, since X is compact, it follows from the dominated convergence theorem that d L 2 (T µ ν , T , which proves the generalizations in question. Finally, setting t = h and using Brenier's theorem we have T µt ν = ∇φ * t for all t and hence the proof of the theorem is concluded by invoking the general form of the inequality 4.2. q.e.d.

Formulation in terms of computational geometry
In this section we show how to use Theorem 1.1 to obtain a quantitative approximation of the solution φ to the second boundary problem for the Monge-Ampère operator, using computational geometry, as developed in [3,4,32,26,31,28].
Assume given a discrete probability measure µ h on R n and a convex bounded domain Y of unit volume. We write Let φ h be the (finite) convex function on R n (uniquely determined up to an additive constant) by which we identify with a function on R n which is equal to φ h on {x 1 , ..., x N } and ∞ elsewhere. We will denote by X h the bounded convex polytope defined as the convex hull of {x 1 , ..., x N } : We will refer to the following equation for a vector φ ∈ R N as the discrete Monge-Ampère equation (associated to the discrete measure µ h on R n and the measure ν on Y ) : are the cells in the partition of R n induced by the piecewise affine function φ * (see point 6 and 7 in Section 2.2).
By the following proposition the previous equation is solved by the vector φ h and moreover, the graph of φ h over X h is simply the convex hull of the discrete graph over {x 1 , .., x N } of φ h . Proposition 5.1. Denote by ψ h := φ * h the Legendre transform of φ. Then • The vector φ h solves the discrete Monge-Ampère equation 5.1 and on Y the convex function ψ h is piecewise affine and given by • On the convex hull X h of {x 1 , ..., x N } the convex function φ h is piecewise affine and given by where {y i } M i=1 runs over the vertices in the polyhedral decomposition of Y determined by the convex piecewise affine function ψ h on Y. Moreover, the graph of φ over X h is equal to the convex hull of the graph of the discrete function φ(x i ) over X h . • The following duality relations hold: (∂φ h )(X h ) = Y and the multi-valued maps ∂φ h and ∂ψ h , which are inverses to each-other, yield a bijection between the facets and the vertices of the induced tessellations of X h and Y .
R n with support compact X. Assume that φ is a function on X such that (∇(χ X φ) * ) * ν = µ. Then φ is uniquely determined (mod R) a.e. wrt µ (and moreover, equal to the restriction to X of a convex function on R n with subgradient image in Y ). Indeed, setting ψ := (χ Y φ) * the argument in the proof of Lemma 5.3 gives that the convex function P Y φ := (χ Y ψ) * (= (χ Y φ) * ) on R n is uniquely determined and satisfies M A g (P Y φ) = µ. But, by general principles, we have that showing that φ is uniquely determined (mod R) on X. The vanishing 5.5 can, for example, by shown by noting that P Y φ is the upper envelope of all functions ϕ in C Y (R) dominated by φ. Then 5.5 follows from a local argument using the maximum principle for the Monge-Ampère operator (just as in the more general case of the complex Monge-Ampère operator considered in [8, Prop 2.10]). q.e.d.
There has been extensive numerical work on constructing solutions to the discrete Monge-Ampère equation 5.1 on R N , based on computational geometry [3,4,32,31,28]. In particular, it was shown in [28] that the solution φ h ∈ R N is the limit point of a damped Newton iteration, which converges globally at a linear rate (and locally at quadratic rate if g is Lipschitz continuous). Performing one step in the iteration, φ (m) → φ (m+1) , amounts to computing the convex hull in R n+1 of the discrete graph over {x 1 , ...x N } defined by φ (m) (or equivalently, the corresponding power diagram in R n and its intersection with Y ; see Section 5.1 below).
Here we show that Theorems 1.1, 1.2 yield a quantitative rate of convergence, in the limit when the spatial resolution h → 0, for the solutions φ h and the corresponding piecewise constant maps T h from X h to R n defined as follows. Consider the tessellation of X h by the convex polytopes F * j obtained by projecting to R n the facets of the "upper boundary" of the convex hull in R n+1 of the discrete graph Γ h of the solution φ h ∈ R N . Then T h maps F * j to the dual vertex y j ∈ R n . In geometric terms, y j is the projection to R n of the corresponding normal of the graph Γ h in R n+1 (normalized so that its projection to R is −1).
Theorem 5.4. Assume that (X, Y, µ, ν) is regular and let φ h ∈ R N be the normalized solution to the discrete Monge-Ampère equation 5.1 and denote by T h the corresponding piecewise constant maps from X h to R n . Then there exists a constant C, depending on (X, Y, µ, ν), such that where T is the optimal diffeomorphism transporting µ to ν. If moreover, the sequence of pointclouds (x 1 , ..., x N ) is quasi-uniform, i.e. there exists a constant A such that In the general case when X and Y are bounded domains in R n with Y assumed convex and ν a probability measure of the form 1.1 the estimates above hold if h is replaced by h 1/2 (n−1) and then the constant C only depends on upper bounds on the diameters of X and Y and the lower bound δ (if the volumes of X and Y are normalized).
Proof. Combing Theorems 1.1, 1.2 with Prop 5.1 and Lemma 5.3 immediately gives the first inequality in the theorem. Then the last one follows from the corresponding L 2 −estimate 1.6 and the fact that v h := φ h − φ is Lipschitz continuous with a uniform Lipschitz constant L.
To see this, set v h := φ − φ h , which we can extend to a Lipschitz continuous function with the same Lipschitz constant L on a neighborhood ofX. Then, the L 2 −estimate 1.6 yields Vol(B δ/2 (x i )), which, concludes the proof, thanks to the quasi-uniformity assumption. q.e.d.
As recalled in the following section the pairs (F * j , y j ) M j=1 defining the map T h correspond, from the point of view of computational geometry, to the facets and the dual vertices of the weighted Delaunay and Voronoi tessellations of X h and R n , respectively, determined by the points {x 1 , .., x N } and the vector φ h ∈ R N .

5.1.
Relations to duality in computational geometry. Proposition 5.1 is closely related to the well-known duality in computational geometry between weighted Voronoi tessellations (aka Laguerre tessellations or power diagrams) of R n and weighted Delaunay tessellations [3] of convex polytopes (see also [26]). To briefly explain this we first note that, by definition, the facets F Completing the square this means that F φ * h i = y ∈ R n : |x i − y| 2 + w i ≤ |x j − y| 2 + w j ∀j w i := 2φ h (x i ) − |x i | 2 , which is the definition of the cells in the weighted Voronoi tessellation of R n associated to the weighted points (x i ; w i ) N i=1 ) [3]. The corresponding weighted Delaunay tesselation is defined as the projection to X h of the polyhedral cell-complex in R n+1 defined by the convex hull of the graph of φ h . The duality in Proposition 5.1 implies that the corresponding sub-gradient map ∂φ h maps the facets F * i of the weighed Delaunay tesselation of X h to the vertices y i of the weighted Voronoi tessellation of R n . More generally, p−dimensional faces correspond to (n − p)−dimensional faces [3].
Remark 5.5. By the duality above, computing a convex hull of a discrete graph over N points in R n is equivalent to computing the corresponding weighted Voronoi tessellation (as emphasized in [3]). For example, when n = 2 the time-complexity is O (N log N ). Moreover, if g is constant and Y is a polytope then the computation of the volumes of the corresponding cells has timecomplexity O(N ) (since the cells are convex polytopes). More generally, efficient computation of the volume can be done if the density g is piecewise affine. However, the Newton iteration referred to above (as opposed to a steepest descent iteration) also requires computing the volumes of the intersections of the cells [32,31]. This has worst-case time-complexity O(N 2 ), but the experimental findings in [32,31] indicate much better near linear time-complexity.

The periodic setting
Assume given a Z n −periodic measures µ and ν on R n normalized so that their total mass on a (or equivalently any) fundamental region is equal to one. We will assume that ν has a positive density g with positive lower bound δ > 0. We will identify µ and ν with probability measures on the torus M := ( R Z ) n . We endow M with its standard flat Riemannian metric induced from the Euclidean metric on R n and denote by dx the corresponding volume form on M. Setting u(x) := φ(x) − |x| 2 /2 gives a bijection between quasi-convex functions φ on R n (i.e. such that ∂φ is periodic) and functions u ∈ C 0 (M ) which are quasi-convex in the sense that ∇ 2 u + I ≥ 0 holds in the weak sense of currents.
When u ∈ C 2 (M ) and µ = f dx the Monge-Ampère 1.3 for a quasi-convex function φ on R n is equivalent to the following equation for u : g(I + ∇u(x)) det(I + ∇ 2 u)dx = f dx We will say that (µ, ν) is regular if such a solution exists. This is the case if both f and g are Hölder continuous and strictly positive [19]. In the case of general (µ, ν) there always exists a weak solution, in the sense of Alexandrov, which is unique modulo an additive constant [15].
Let now µ h be a family of discrete measures on M, defined as in formula 1.4 (with h denoting the corresponding mesh norm). Then the following analog of Theorem 1.1 holds: