Monotone and Consistent discretizationof the Monge-Ampere operator

Monotone and Consistent discretization
of the Monge-Ampere operator

Jean-David Benamou, Francis Collino111Mokaplan, INRIA, Domaine de Voluceau BP 105 78153, Le Chesnay Cedex France , Jean-Marie Mirebeau222CNRS, University Paris Dauphine, UMR 7534, Laboratory CEREMADE, Paris, France.
E-mail: mirebeau<at>ceremade.dauphine.fr, jean-david.benamou<at>inria.fr
Key words: Monge Ampere PDE, Monotone Finite differences scheme, Lattice Basis reduction, Stern-Brocot tree.
Mathematics Subject Classification (2000): 35J96, 65N06
Abstract

We introduce a novel discretization of the Monge-Ampere operator, simultaneously consistent and degenerate elliptic, hence accurate and robust in applications. These properties are achieved by exploiting the arithmetic structure of the discrete domain, assumed to be a two dimensional cartesian grid. The construction of our scheme is simple, but its analysis relies on original tools seldom encountered in numerical analysis, such as the geometry of two dimensional lattices, and an arithmetic structure called the Stern-Brocot tree. Numerical experiments illustrate the method’s efficiency.

1 Introduction

We introduce a new discretization of the Monge-Ampere operator, on two dimensional cartesian grids, which is consistent and preserves at the discrete level a fundamental property of the continuous operator: degenerate ellipticity. Discrete degenerate ellipticity [Obe06] implies strong guarantees for the numerical scheme: a comparison principle, convergence of discrete solutions towards the continuous one in the setting of viscosity solutions, and convergence of Euler iterative solvers for the discrete system [Obe06]. Some Degenerate Elliptic (DE) schemes for the Monge-Ampere (MA) Partial Differential Equation (PDE) already exist [FO11, Obe06], but they suffer from several flaws: they are strongly non-local, and only approximately consistent. Consistent non DE schemes such as [LR05, BN12] offer better accuracy, but require the PDE solution to be sufficiently smooth and the discrete numerical solver to be well initialized. Filtered schemes [FO13] nonlinearly combine several existing schemes, in order to cumulate their advantages (here degenerate ellipticity and consistency), or mitigate their defects. Their definition and their analysis are however complex, and their application requires to adjust several parameters. For a recent overview of the numerical approaches to solving the Monge-Ampère equation, see Glowinski, Feng and Neilan [FGN13].

We introduce a new numerical scheme, Monge-Ampère using Lattice Basis Reduction (MA-LBR), which is both consistent333Assuming the solution hessian condition number is uniformly bounded and degenerate elliptic. Lattice Basis reduction is a tool from discrete geometry, which arises here due to the interaction of the cartesian discretization grid, with the anisotropic nature of the Monge-Ampere operator. This operator is indeed invariant under all linear changes of variables with unit determinant, unlike e.g. the Laplacian which is merely invariant under orthogonal transformations. The MA-LBR belongs is inspired by the Wide-Stencil [Obe06] family of schemes. Using another arithmetic tool, the Stern-Brocot tree, we solve a second issue plaguing these methods (in addition to consistency errors): our discretization stencil needs not be chosen a priori (which usually involves a difficult arbitrage between scheme locality, consistency error and available CPU time), but can be generated automatically in a guaranteed, parameter free and solution adapted manner. Numerical experiments §4 illustrate the MA-LBR accuracy and robustness.

We fix throughout this paper a convex open bounded domain . Given a density , and some Dirichlet data , we set the goal of approximating numerically the unique viscosity solution [CIL92, Gut01] of

(1)

Our framework admittedly does not encompass solutions of the weaker Alexandrov type, where is merely a non-negative measure. If is convex but not strictly convex, then the Dirichlet data is assumed to be convex on any segment of . Let us point out that optimal transport, from to another domain , equipped with densities , , admits a PDE formulation similar in spirit to (1): , , convex. The gradient non-linearity, and the second boundary condition, raise difficulties [Urb97, BFO14] that we choose not to address in the present paper, focusing instead on the Monge Ampere operator .

We assume that the PDE domain is discretized on a cartesian grid: , where is the grid scale, is an arbitrary rotation, and is an offset. For notational simplicity, and up to a linear change of coordinates, we limit our attention to the canonical values of these parameters, so that the discrete domain is

Definition 1.1.

We denote by the collection of discrete maps . A (discrete) operator is a map . It associates to each a collection of values , .

The notations and refer to the same object, which is a map , and are used interchangeably with the aim of improving readability. In numerical experiments, the values of on are the unknowns, while the values on are the supplied boundary data: . For each we introduce a second order differences operator , built so that , where and , and where denotes the euclidean scalar product on . In the simplest case where , we set

(2)

When is close to , the points or may not belong to . Denoting by the only element of such that , we define

(3)

Let us again point out that if , then the value is the supplied boundary data . On the other hand if , then (2) and (3) coincide. No other consistent approximation of can be built using the values , and .

Discretizations of the Monge-Ampere operator are typically built upon the operators . Consider for instance the Finite Differences (FD) discretization [LR05]

(4)

Given such a discrete operator , the discrete analog of (1) takes the form:

Find , such that on , and . (5)

This discrete system lacks a counterpart of the constraint of convexity in (1) because (i) there is no unique notion of discrete convexity but several competing approaches, see for instance [Mir14a, Obe13], and (ii) some form of discrete convexity constraint can often be embedded in the equation , see §1.2. From a theoretical and a practical standpoint, choosing in (5) is a risky bet: second order convergence can often be observed in numerical experiments, see §4, but only on rather easy cases and with a good initialization for the numerical solver. Robustness results (existence, uniqueness, and algorithmic guarantees) are limited to discretizations obeying an additional property: a counterpart of the ellipticity of the (opposite of the) Monge-Ampere operator .

We use the notion of discrete degenerate ellipticity [Obe06], slightly specialized due to our focus on MA. Degenerate Elliptic Monge-Ampere numerical schemes cannot be strictly local, unlike (4), but instead need to take into account some long range second order differences, indexed by a possibly wide stencil.

Definition 1.2.

A stencil is a finite set which is symmetric with respect to the origin (i.e. for each ).

Definition 1.3.

(DE2 scheme) A numerical scheme is Degenerate Elliptic, with stencil , iff for each the quantity is a non-decreasing, locally Lipschitz function of the second order differences , .

Figure 1: Examples of stencils , containing , , and elements respectively. Wider stencils yield smaller consistency errors, see Figure 3.

Observing that the second order difference can be expressed as a non-negative weighted sum of first order differences (3), we immediately find that a DE2 scheme is degenerate elliptic in the sense of [Obe06]. DE2 schemes are also positive difference operators in the sense of [KT92] in this paper schemes are indeed built using directional second order finite differences. In particular, for any , the slightly perturbed operator , defined by , is proper degenerate elliptic [FO13]. This in turn implies that the discrete system (5) associated with has a unique solution, which can be computed with a geometric convergence rate using an iterative Euler scheme. We refer to [FO13] and references therein for these results and will say no more on this analytic machinery in the rest of the paper, focusing instead on the algebraic structure of discrete Monge-Ampere operators.

Froese and Oberman [FO11] numerically address the MA PDE using a DE2 operator, referred to as the Wide Stencil (WS) scheme. Given a stencil , and denoting :

(6)

The minimum is taken over all pairs of vectors which are orthogonal, in the sense that . For instance , or . We introduce a variant of this operator, which does not rely on pairs of orthogonal stencil vectors, but on superbases of the lattice .

Definition 1.4.

A basis of is a pair such that .
A superbase of is a triplet such that , and is a basis of .

Figure 2: Left: A Stencil , a superbase , an orthogonal pair . Right: -obtuse superbase, see Definition 1.6, and ellipse for some .

The MA-LBR operator, associated to a stencil , is defined by

(7)

where for we define

(8)

Remark 1.8 provides a geometric interpretation for these at first abstruse formulas. The operators (6) and (7) are DE2 since the product is non-decreasing in each variable, as well as the function , see Lemma 3.6.

Outline.

We discuss in §1.1 the consistency of the MA-LBR, and show in particular that a finite stencil is sufficient to achieve consistency for all quadratic functions of condition number below a given bound. For more simplicity and efficiency we introduce in §1.2 an automatic stencil construction for the MA-LBR, which is adaptive, local, anisotropic, parameter free, and has good consistency guarantees. The proofs of the results appearing in §1.1 and §1.2 are postponed to §2 and §3 respectively.

Notations.

For each we denote . If then . Given pairwise distinct , we denote by the segment of endpoints , and by the triangle of vertices .

1.1 Consistency

Figure 3: Relative consistency error for quadratic functions (9), with the schemes MA-LBR (top) and WS (bottom), using the stencils of Figure 1. See (10) for the parametrization of symmetric matrices, by their condition number and their orientation . Note that the MA-LBR consistency error vanishes for a large set of matrices.

The consistency analysis of the numerical schemes FD, FO and the MA-LBR reveals significant differences. We denote by the collection of symmetric matrices of size , and by those which are positive definite. For each we introduce a quadratic map , defined by , . Since the second order difference operator is consistent, for any , it is exact for . Summarizing one has

(9)
Definition 1.5.

The consistency set of an operator is the collection of matrices for which , identically on .

One easily checks that the consistency set of the finite differences discretization , see (4), is the whole . In fact the identity also holds for non-definite matrices , although they are irrelevant for our application. Since that scheme is not DE, this consistency does not imply convergence results. As illustrated in Figures 3 and 4, schemes WS and MA-LBR have in contrast non-trivial consistency sets, depending on the chosen stencil. Matrices are parameterized in these figures by their condition number , and the orientation of their first eigenvector :

(10)

The consistency analysis of is based on Hadamard’s theorem [FO11]: for all , and any pair of non-zero orthogonal vectors, one has , with equality iff and are eigenvectors of . As a result, scheme is only consistent on a negligible subset of : those matrices which eigenvectors lie in , see Figures 3 and 4. From a theoretical standpoint, convergence results are obtained in [FO11] by increasing the stencil size, up to infinity, as the discretization grid scale tends to zero. In practical cases, finding the optimal stencil size is non-trivial, see §4.

Figure 4: Left: element of largest euclidean norm of an -obtuse superbase. Center left: an eigenvector of . The consistency set of (resp. ) is the union of the regions (resp. the dashed lines) corresponding to all stencil elements . Right: triangulations mentioned in Remark 1.8, and some 3D visualisations of the maximal convex extension .

The key concept in the MA-LBR consistency analysis is the notion of -obtuse superbase, which originates from lattice geometry [CS92] (a lattice is a discrete subgroup of containing a basis, such as ). It was already applied to PDE discretizations in [Mir14b, FM13].

Definition 1.6.

Let . A superbase of is said -obtuse iff for all .

Theorem 1.7 (Consistency).

A matrix is in the consistency set of iff there exists which form an -obtuse superbase.

The following remark attempts to give a geometrical interpretation of the function (8) and of Theorem 1.7. The results of this section, Theorem 1.7, Remark 1.8 and Theorem 1.9, are established in §2.

Remark 1.8 (Geometric interpretation).

Let , and let be a superbase of . Let be the maximal convex map bounded above by at the points . Then (the Lebesgue measure of the subgradient of at , which is a natural relaxation of the Monge-Ampere operator [Gut01]). The map is polygonal, on one of the four triangulations illustrated Figure 4. The identity holds for the first triangulation only, which corresponds to an -obtuse superbase .

Strikingly, one cannot hope for a DE2 scheme more localized than the MA-LBR. Finding well localized numerical schemes, involving small stencils, is a natural objective [Koc95].

Theorem 1.9 (Minimality).

Let be a DE2 scheme with stencil . If the consistency set of contains the neighborhood of a matrix , then there exists which form an -obtuse superbase.

The following algorithm and proposition, dating back to Selling [Sel74, CS92], constructively shows the existence of an -obtuse superbase for each , without which Theorems 1.7 and 1.9 would be mostly vacuous. It is worth noting that this algorithm extends to dimension three [CS92]. Proposition 1.10 also immediately implies that all matrices with condition number are simultaneously in the consistency set of the MA-LBR operator with stencil

(11)
Initialize ,  ,  . (Or any other initial superbase.)
While the superbase is not -obtuse do
Find such that , and set .
Algorithm 1 Construction of an -obtuse superbase (Selling [CS92]).

In order to analyse this algorithm, we associate to each the norm

(12)
Proposition 1.10 (Existence of an -obtuse superbase, Selling 1874).

Algorithm 2 terminates, and the final state of is an -obtuse superbase. Furthermore for each .

Proof.

To each superbase associate the energy . One easily checks that , for any . Denoting by the successive superbases generated in Algorithm 1, we observe that the energies are strictly decreasing by construction. Noticing that there exists only a finite number of superbases with energy below a given bound, we find that the algorithm terminates. At termination, the continuation criterion “the superbase is not -obtuse” is false, which establishes the first point.

Let , and let be the number of loop iterations. Then

(13)

which immediately implies the announced bound on the obtuse superbase elements norm. ∎

The MA-LBR consistency error is typically smaller than with the WS scheme, for a given stencil , see Figure 4. Furthermore while the WS consistency is an asymptotic property, depending on the stencil angular resolution, the MA-LBR has in contrast a consistency set of non-empty interior, and its elements can be identified with a simple test, see Theorem 1.7. Unfortunately, choosing the MA-LBR effective stencil before a numerical simulation remains at this point a puzzle for the practitioner. The option (11) is not practical because: (a) Uniform bounds on the hessian matrix condition number of solutions to (1), are seldom available. (b) Even if such a bound is available, the set (11) can be quite large, with cardinality . This becomes an issue if the bound is pessimistic, or if the solution hessian does degenerate in some places, such as along the domain boundary . A third issue (c) is that there is no clear way to a-posteriori validate the choice of a given stencil: would the numerical solution be improved with a larger one ?

Selling’s algorithm, in contrast with the inefficiency of (11), adaptively produces an -obtuse superbase in only few iterations. We present in the next section an adaptive, anisotropic, parameter free and guaranteed stencil refinement algorithm, which eliminates the implementation difficulties (a), (b), (c) above. Under the hood, it amounts to an adaptation of Selling’s algorithm to non-quadratic functions, see §3.3.

1.2 Hierarchical stencil refinement

The previous section fully characterized the consistency set of the MA-LBR operator , associated to a stencil . Larger stencils provide consistency on larger collections of matrices, as established in Theorem 1.7, and illustrated on Figure 4. Excessively large stencils are however unpractical, since the CPU cost of evaluating the MA-LBR operator (7) is proportional to their cardinality. Adapting Selling’s obtuse superbase construction, Algorithm 1, we show that one can emulate an MA-LBR with extremely large stencils for a limited numerical cost.

Our adaptive variant of the MA-LBR operator is defined by Algorithm 2 below, which is lines long and only involves elementary operations. Its analysis (and Definition 1.16 of a mild structural constraint on stencils) relies on an arithmetic construction named the Stern-Brocot tree, already used in [BOZ04, Mir13] for the discretization of anisotropic PDEs. Definitions 1.12, 1.14 introduce this structure. Propositions 1.11 and 1.15 are variants of commonly known facts on the Stern-Brocot tree which proof is, for completeness, presented in the appendix.

Proposition 1.11.

The identity defines a one to one correspondance between:

  • Vectors , such that and .

  • Direct acute bases of (i.e. , and ).

Definition 1.12.

We emphasize the unique décomposition introduced in Proposition 1.11 by using the notation Whenever we write , we implicitly limit our attention to vectors satisfying the assumptions of Proposition 1.11.

For instance , , and . If , then is a superbase of ; all superbases happen to be of that form, up to a permutation of their elements, see Lemma A.3. The next proposition shows how to generate numerous decompositions of the form of Proposition 1.11.

Proposition 1.13.

If , then and .

Proof.

We check , and . Hence is a direct acute basis of . Likewise for . ∎

Definition 1.14.

We introduce a graph , with vertices , and edges and for each .

We say that an edge of leaves from and arrives at . We denote by the eight point stencil illustrated on Figure 1 (left).

(14)
Proposition 1.15.

The set has one element in each connected component of . The four points of the form or are isolated. The four other points are the root of complete infinite binary trees, each one entirely contained in a quadrant of the plane.

Figure 5: Stern Brocot binary tree (left), and its local structure (center left). Right: some convex domains, and the associated discretization grids. Plain black arrows: set , for some . Light blue arrows: vectors such that , as in property (Reachability) of stencils.

The Stern-Brocot tree is the subgraph corresponding to the first quadrant, with vertices , see Figure 5. This complete infinite binary tree originates from arithmetic, and in the literature a vertex is often identified with the positive irreducible fraction .

The MA-LBR adaptive variant, presented below, requires stencils with a special structure. For each we introduce the set

(15)

where “” stands for “ and ”. Note that the continuous domain could in (15) be replaced with the discrete one : indeed for any , one has iff . Any subset of is regarded as a subgraph of , equipped with all edges of having their endpoints in .

Definition 1.16.

A family of stencils is the data, for each , of a stencil satisfying , and the additional structural properties:

  • (Hierarchy) The set has an element in each connected component of .

  • (Reachability) contains each such that .

Practical recommendations regarding the construction of stencils are discussed after Theorem 1.19. These structural constraints are in practice not hard to satisfy. Condition (Hierarchy) is natural in view of the Stern-Brocot tree structure, and is satisfied by all stencils illustrated on Figure 1. Condition (Reachability) ensures that the stencils can be extended, in the sense of Proposition 1.18 below. It is vacuous for all at distance from in general, and entirely vacuous in the case of a box domain, see Proposition 3.4 and Corollary 3.5.

The sets , , are small or empty when is close to , but typically huge when is far from , see Figure 5. They do not constitute a family of stencils, but they can be used to extend an existing family of stencils, as in the next definition.

Definition 1.17.

To each family of stencils , we associate the family of sets defined by , .

Proposition 1.18 (Extension of stencils).

If is a family of stencils, then also is.

We next introduce the MA-LBR operator associated to a family of stencils, as well as a hierarchical variant . (The expected superscript is omitted for readability.)

(16)
Initialize a variable , and list . Set also
While is non-empty do
bla Denote by the first element of , and set .
bla If , or [ and ] (Refinement test)
blabla then prepend to , and set
blabla else remove from and set .
Algorithm 2 Hierarchical operator (the final value of ).

Our main result Theorem 1.19 states that the MA-LBR operator associated to the large stencils coincides in all cases of interest with the hierarchical, adaptive variant . Algorithm 2 amounts to a depth-first transversal of a finite subtree of the Stern-Brocot tree, see §3.2 and [Mir13] where a similar approach is used for the discretization of Hamilton-Jacobi PDEs. This subtree is characterized by the stopping criterion (Refinement test), allowing to reject useless branches of where the minimum (16) defining cannot be attained. In the case of a quadratic map , , Algorithm 2 explores a single branch of the Stern-Brocot tree, just as Selling’s algorithm, see §3.3.

We say that a property holds “on ” iff it holds at each point of .

Theorem 1.19 (Adaptive pruning equals extensive sweeping).

Let be a family of stencils, and let . If on , or on , then we have on .

In our experiments §4 with , we use reasonably large stencils on a layer of a few grid points along , where is small or empty, see Figure 5. We use in contrast the minimal stencils elsewhere since they are adaptively completed by Algorithm 2.

The identity may break down when these two operators vanish at some points of . This is fortunately not an issue since (i) the problem of interest (5) has by assumption a positive right hand side, and (ii), the positivity of the MA-LBR operator is equivalent to the positivity of second order differences, see Proposition 1.20 below, which is a natural discrete counterpart of the convexity constraint present in the original Monge-Ampere problem (1).

Consider a smooth function , and a point such that . Then is convex (or concave) on a neighborhood of . The next proposition establishes a discrete analog of this property. Consider and a family of stencils. The discrete counterpart of is , while the counterpart of the convexity of locally around is the positivity of the second order differences centered at : , .

Proposition 1.20 (Discrete convexity).

Let , let be a family of stencils on , and let . Then

Oberman [Obe13] numerically addressed variational problems posed on the cone of convex functions by imposing the positivity of second order differences, for all points and all vectors within some given stencil . It is also known, see Appendix A of [Mir14a], that any discrete map satisfying whenever , needs to coincide with a global convex function on the subsampled grid of points with even coordinates.

Remark 1.21.

Adaptivity in PDE discretizations often refers to the context where a sequence of discrete maps is generated along an iterative procedure, as well as a sequence of operators, and depends on . Our understanding in this paper is different: there is no underlying iteration, but a single operator which is evaluated in a subtle and cheap way as .

2 Proof of consistency and minimality

We establish the results announced in §1.1, and related to the MA-LBR consistency and optimal locality. Theorem 1.7 (Consistency) and Remark 1.8 are proved in §2.2, and Theorem 1.9 (Minimality) in §2.2.

2.1 Consistency

Our first result, Proposition 2.2 preceded with a technical lemma, shows that the MA-LBR operator (7) systematically overestimates the hessian determinant of quadratic functions. For any , defining as in (9), one has on . Equality holds iff contains an -obtuse superbase, which establishes the announced Theorem 1.7 (Consistency). We denote

(17)
(18)
Lemma 2.1.

Let . Then , with equality iff .

Proof.

If then by definition (8). Otherwise, we may assume without loss of generality that , so that and . ∎

Proposition 2.2.

Let , let be a superbase of , and let for . Then . Equality holds iff , equivalently iff is -obtuse.

Proof.

Given a permutation of we compute

Hence iff the superbase is -obtuse. We prove in the following that , which in view of Lemma 2.1 concludes the proof.

Special case of the superbase , , , with . We get , , . Inserting this into the expression (18), yields as announced .

General case. Let be a matrix such that and , so that by linearity . Note that . We obtain , for all , so that by the special case . ∎

The rest of this subsection is devoted to the proof of Remark 1.8. Let be a fixed superbase of . For each we introduce a polygon , defined by linear inequalities, and some of its edges ,

The area of is computed in Corollary 2.4, and this polygon (properly scaled and translated) is identified with a subgradient set in Proposition 2.5, concluding the proof of Remark 1.8. The proof unfortunately gives little geometric insight, hence it could be skipped at first reading. Given , and we use the notation .

Lemma 2.3.

Let . Then is a segment of length (i) if , (ii) if , (iii) if , or (iv) if (a case where is in fact empty).

Proof.

Let and let . One has iff

The equality is equivalent to . Recall that , and likewise . Hence the inequalities respectively hold iff belongs to the segment

Translating these two segments by yields , and . Finally the length of is

(19)

which coincides with the announced result. ∎

Corollary 2.4.

For any , one has .

Proof.

The triangle has area : half the height from the vertex at the origin, times the length of the opposite side. These three (possibly empty) triangles, with their opposites, partition . From this point the result follows from Lemma 2.3 and an easy calculation. ∎

Proposition 2.5.

Let , and let , . Let , and let be the maximal convex map bounded above by at the points and , . Then .

Proof.

For any , the following are equivalent:

  • , for all points of the hexagon of vertices , .

  • , for all , , .

In order to further simplify this expression, we write , , where , , and insert the expression (9) of . The following are then equivalent:

We recognize the inequalities defining , and the announced result follows. ∎

Proof of remark 1.8.

Proposition 2.5 and Corollary 2.4 imply as announced that . By Proposition 2.2 one has iff , which by Lemma 2.3 means that is an hexagon: each edge has a positive length (we exclude here for simplicity the limit case ).

The map is polygonal on a triangulation with vertices , , , which is symmetric with respect to . Only four such triangulations exist, as illustrated on Figure 4, and only the first one leads to an hexagonal subgradient , since the subgradient has one vertex for each triangle containing . This concludes the proof. ∎

2.2 Minimality

We prove Theorem 1.9 (Minimality), on the optimal locality of the MA-LBR. For that purpose we introduce some definitions, and establish in Proposition 2.8 a minimality property of obtuse superbases.

Definition 2.6.

We denote by , the closed convex cone spanned by two elements . We say that are trigonometrically consecutive elements of a set iff they are not collinear and no element of lies in the interior of .

Definition 2.7.

A matrix is said generic iff there exists no -orthogonal basis of . (i.e.