A Convergent Quadrature Based Method For The Monge-Ampère Equation

by   Jake Brusca, et al.
New Jersey Institute of Technology

We introduce an integral representation of the Monge-Ampère equation, which leads to a new finite difference method based upon numerical quadrature. The resulting scheme is monotone and fits immediately into existing convergence proofs for the Monge-Ampère equation with either Dirichlet or optimal transport boundary conditions. The use of higher-order quadrature schemes allows for substantial reduction in the component of the error that depends on the angular resolution of the finite difference stencil. This, in turn, allows for significant improvements in both stencil width and formal truncation error. We present two different implementations of this method. The first exploits the spectral accuracy of the trapezoid rule on uniform angular discretizations to allow for computation on a nearest-neighbors finite difference stencil over a large range of grid refinements. The second uses higher-order quadrature to produce superlinear convergence while simultaneously utilizing narrower stencils than other monotone methods. Computational results are presented in two dimensions for problems of various regularity.


page 1

page 2

page 3

page 4


A C^0 finite element method for the biharmonic problem with Dirichlet boundary conditions in a polygonal domain

In this paper, we study the biharmonic equation with Dirichlet boundary ...

Error estimation of the Relaxation Finite Difference Scheme for the nonlinear Schrödinger Equation

We consider an initial- and boundary- value problem for the nonlinear Sc...

Singular integration towards a spectrally accurate finite difference operator

It is an established fact that a finite difference operator approximates...

A Neural Network Based Method to Solve Boundary Value Problems

A Neural Network (NN) based numerical method is formulated and implement...

cuSten -- CUDA Finite Difference and Stencil Library

In this paper we present cuSten, a new library of functions to handle th...

On the Reduction in Accuracy of Finite Difference Schemes on Manifolds without Boundary

We investigate error bounds for numerical solutions of divergence struct...

1. Introduction

In this article we introduce an integral representation of the Monge-Ampère equation


where is convex and the right-hand side is non-negative. This allows us to produce new monotone approximation schemes via quadrature. Because these schemes are monotone, they fit within several recently developed numerical convergence frameworks [4, 5, 14, 20, 21, 27]. Moreover, these new schemes offer significant advantages over existing monotone methods in terms of both accuracy and efficiency.

Recent years have seen a growing interest in Monge-Ampère type equations in the context of a diverse range of problems including design of optical systems [31], geophysics [12], mesh generation [7], medical image processing [17], meteorology [10]

, and data science 

[29]. This has encouraged the design of many new methods for the Monge-Ampère equation including [3, 6, 11, 13, 30].

The development of numerical methods that are guaranteed to converge to the correct solution, particularly in the absence of classical solutions, has proven to be more challenging. An early method [28] used a geometric interpretation of weak solutions to design a convergent, but computationally expensive, method for the 2D Monge-Ampère equation. Recently, convergence frameworks have been established for the Monge-Ampère equation with either Dirichlet boundary conditions [14, 20, 26, 27]:


or the second type boundary condition arising in optimal transport [4, 5, 21]:


These convergence proofs can be viewed as extensions of the powerful Barles and Souganidis convergence framework [1]

, which is valid for weak (viscosity) solutions of fully nonlinear partial differential equations. Critically, they are only valid for approximation schemes that are monotone. Construction of monotone schemes for degenerate elliptic PDE operators is not trivial: in fact, given any fixed finite difference stencil, it is possible to find linear elliptic operators for which no consistent, monotone approximation is possible on the given stencil 

[23, 25]. Circumventing this challenge requires the use of finite difference stencils that are allowed to grow wider as the grid is refined. Several monotone finite difference schemes are now available for the Monge-Ampère equation [5, 16, 14, 24, 27]. Because of the wide-stencil nature of these methods, the methods are computationally expensive and typically have low (sub-linear) accuracy. These challenges are magnified in three dimensions, where even evaluating the finite difference approximations (without attempting to solve the resulting nonlinear system) can be prohibitively expensive [18].

In this article, we propose to express the Monge-Ampère operator in terms of a Gaussian integral. This allows us to utilize higher-order quadrature schemes in order to simultaneously achieve improved consistency error and narrower finite difference stencils. The schemes are nevertheless monotone, and fit neatly within the existing proofs of convergence of numerical methods to the weak (viscosity) solution of the Monge-Ampère equation. We describe two different implementations of this approach in two dimensions and validate the performance using a range of standard benchmark problems for the Dirichlet problem. This new formulation of the Monge-Ampère equation holds particular promise for the development of computationally practical methods in three dimensions, as it provides a dimension-reduction as compared with a typical variational formulation of the 3D Monge-Ampère equation. It also extends naturally to more general Monge-Ampère type equations in optimal transport, including equations that are posed on the sphere [19].

2. Background

2.1. Elliptic equations

The Monge-Ampère equation is an example of a degenerate elliptic partial differential equation, which take the general form

Definition 1 (Degenerate Elliptic).

Let and denote by the set of symmetric matrices. The operator is said to be degenerate elliptic if

whenever and .

We note that the operator is defined on the closure of , and takes on the value of the relevant boundary conditions at . For the Dirichlet problem, which is the setting implemented in this article, the PDE operator at the boundary is defined as


The Monge-Ampère equation (1) does not immediately satisfy this definition of an elliptic equation; in fact, it holds only on the restricted class of convex functions. Going hand-in-hand with this difficulty is the fact that the solution of the Monge-Ampère equation is not expected to be unique; the additional constraint that is convex is needed in order to select a unique solution. A common remedy to these challenges is to define a globally elliptic extension of the Monge-Ampère equation that automatically enforces solution convexity [20]. This can be accomplished by considering the convexified Monge-Ampère operator


Here the modified determinant should agree with the usual determinant when operating on the Hessian of a convex function and should return a negative value otherwise. The particular choice utilized in this article is



are the eigenvalues of the symmetric matrix


In general, degenerate elliptic equations need not have classical solutions, and some notion of weak solution is required. The viscosity solution has proved to be particularly useful for this class of equations [9], and forms the foundation for most of the recently developed numerical convergence proofs for the Monge-Ampère equation. The idea of the viscosity solution is to use a maximum principle argument to pass derivatives onto smooth test functions that lie above or below the semi-continuous envelopes of the candidate weak solution.

Definition 2 (Semi-continuous envelopes).

Let be a bounded function. Then for , the upper and lower semi-continuous envelopes are defined, respectively, as

Definition 3 (Viscosity subsolutions (supersolutions)).

A bounded upper (lower) semi-continuous function is a viscosity subsolution (supersolution) of (4) if for every , that whenever has a local maximum (minimum) at , then

Definition 4 (Viscosity Solution).

A bounded function is a viscosity solution of (4) if is a viscosity subsolution and is a viscosity supersolution.

An important characteristic of many elliptic operators, which immediately yields solution uniqueness, is the comparison principle.

Definition 5 (Comparison Principle).

The operator (4) satisfies a strong comparison principle if whenever is an upper semi-continuous subsolution and a lower semi-continuous supersolution, then on

We remark that the Dirichlet problem for the Monge-Ampère equation (1),(5) does satisfy a comparison principle under reasonable assumptions on the data. However, this is no longer true when the right-hand side depends on the solution gradient or when the second type boundary condition (3) is considered [20, 21].

2.2. Convergence framework

A fruitful technique for numerically solving fully nonlinear elliptic equations involves finite difference schemes of the form


defined on a finite set of discretization points . Many key results on the convergence of finite difference methods to the viscosity solution of a degenerate elliptic PDE are based upon a set of criterion developed by Barles and Souganidis [1].

Definition 6 (Consistency).

The scheme (8) is consistent with (4) if, for any test function and , we have


To a consistent scheme we can also assign a local truncation error.

Definition 7 (Truncation error).

The truncation error of a scheme (8) on a set of admissible functions is a function such that for any there exists a constant such that

for every and sufficiently small .

Definition 8 (Monotonicity).

The scheme (8) is monotone if is a non decreasing function of its last two arguments.

Definition 9 (Stability).

The scheme (8) is stable if there exists some , independent of , such that every solution satisfies .

These simple concepts lead immediately to convergence of finite difference methods, provided the underlying PDE satisfies a strong comparison principle.

Theorem 10 (Convergence [1]).

Let be the unique viscosity solution of the PDE (4), where is a degenerate elliptic operator with a strong comparison principle. Let be any solution of (8) where is a consistent, monotone, stable approximation scheme. Then converges uniformly to as .

This result does apply to the Monge-Ampère equation (1) with Dirichlet boundary conditions (5) under mild assumptions on the data. However, many other Monge-Ampère equations of interest do not possess the strong comparison principle required by the theorem. In recent years, the convergence proof has been adapted to include discontinuous solutions of the non-classical Dirichlet problem [20] and the second boundary value problem [4, 5, 21].

2.3. Wide stencil methods

Several monotone finite difference approximations have been proposed for the Monge-Ampère operator [2, 5, 8, 14, 26]. These hinge upon different reformulations of the Monge-Ampère operator, which typically take a variational form


Here denotes the second directional derivative (SDD) in the direction , is some admissible set, and is a non-decreasing function. Generating a monotone approximation then requires (1) an appropriate discretization of the relevant SDDs and (2) an appropriate discretization of the admissible set.

On a structured grid, where aligned points , , and are available for some , a simple (negative) monotone approximation is


These approximations are typically allowed to have a wide-stencil flavor, with the spacing , being potentially larger than the characteristic spacing of grid points. See Figure 1. We also remark that in the special case of equi-spaced neighboring points (), such as on a uniform Cartesian grid, this reduces to the usual centered difference approximation with second order truncation error. Monotone approximations are also possible on unstructured grids, though they are typically less accurate [16].

Figure 1. Wide finite difference stencils.

The width of stencils required by the approximations of (11) is determined by the discretization of the admissible set . Optimal discretization of this set is itself a non-trivial problem and evaluating (11) may involve minimization over a prohibitively large set of candidate directions, particularly in three dimensions [18]. A typical scaling for the maximal stencil width that optimizes truncation error is at least  [16], though in some cases it is not clear what the optimal choice is.

3. Integral Formulation

In this section, we present an integral representation of which can be used to create a monotone discretizaton through the use of quadrature.

To motivate this, we recall a well know result about the integrals of multivariate Gaussians:


where is a symmetric positive-definite matrix.

This provides an alternate characterization of the Monge-Ampère operator if we let be the Hessian of the potential function . Provided is strictly convex, its Hessian is positive definite and we can write


To express this in a form that is easily discretized, we convert to spherical coordinates. Let and and denote by the second directional derivative of in the direction of . Then the Monge-Ampère operator in can be expressed as


Integrating out the radius , we obtain


For simplicity, the results in the remainder of this paper are presented in two dimensions. However, they are certainly generalizable to higher dimensions. We introduce the notation

and note that

Then we easily obtain the two-dimensional version of (16) in terms of polar coordinates.

Theorem 11 (Integral Representation).

Let be convex and be strictly convex. Then for every :


The characterization in theorem 11 only holds when . However, we are also interested in degenerate cases where . In these cases, we know that has at least eigenvalue equal to zero, and the integrand in (17) becomes singular.

We introduce the following relaxation to approximate the integral in these cases:


This, in turn, is used to construct a relaxed version of the convexified Monge-Ampère operator:


Here we have represented as , which is equivalent via the minimax principle.

4. Quadrature Scheme

In this section, we describe a very general framework for utilizing the integral formulation (15) to produce a consistent, monotone approximation of the two-dimensional Monge-Ampère equation. In section 5, we will describe two particular implementations.

4.1. Approximation scheme

We introduce the following notation.

Definition 12 (Notation).
  1. is a bounded, open, convex domain with Lipschitz boundary .

  2. is a finite set of discretization points , .

  3. is the spatial resolution of the grid. In particular, every ball of radius contained in contains at least one discretization point .

  4. is a stencil width associated to the grid.

  5. is a finite set of angles discretizing .

  6. is the local angular resolution of the discretization, where we define .

  7. is the angular resolution of the discretization.

  8. is the quasi-uniformity constant of the angular discretization.

  9. is a collection of non-negative quadrature weights summing to and satisfying

    for some constant that depends only on the quasi-uniformity constant.

  10. is a regularization parameter associated with the grid.

  11. is the set of neighboring indices for such that for every we have .

  12. described in (12) has the form

    for every , where all .

  13. is the truncation error of the finite difference scheme for approximating the second directional derivative on the admissible set .

  14. is the maximal truncation error of the finite difference approximations.

  15. is the truncation error of the quadrature scheme

    for approximating the integral on the admissible set .

Then we propose the following scheme for approximating the convexified Monge-Ampère operator at interior points .


4.2. Convergence

We now provide conditions under which the scheme (20) is consistent and monotone. As an immediate consequence, it fits directly into the convergence proofs developed in [4, 5, 14, 20, 21, 27].

Theorem 13 (Monotonicity).

The approximation scheme (20) is monotone.


We note that the operator that appears in (20) can be written in the form

where and the are non-negative by the (negative) monotonicity of the approximations .

Let have non-negative components. We notice that

Since the max and min operators preserve monotonicity and the weights are non-negative, we can immediately conclude that

Since has no dependence on its second argument, this completes the proof of monotonicity. ∎

Theorem 14 (Consistency).

Consider discretizations of such that the corresponding parameters

as and the corresponding quasi-uniformity constants are bounded uniformly. Then the approximation scheme (20) is consistent with the convexified Monge-Ampère operator (6).

We will break this result into three separate cases (Lemmas 15-17), depending on the sign of , the smallest eigenvalue of the Hessian. We note that the scheme appearing in (20) has no dependence on the first two arguments, which allows us to simplify slightly the verification of consistency.

Lemma 15 (Consistency with positive eigenvalues).

Under the assumptions of Theorem 14, let and consider such that . Then the scheme (20) satisfies


Since the smallest eigenvalue is strictly positive and , we are assured that

for all sufficiently close to and sufficiently small . Then using the consistency error in the components of this scheme, we can compute

Since is continuous in

and bounded away from zero, the two sums in the last line can both be interpreted as consistent quadrature schemes. Thus we can further estimate

Recalling now the integral formulation of the Monge-Ampère operator (17), we conclude that

Since and all eigenvalues of the Hessian are strictly positive, we conclude that

Lemma 16 (Consistency with a negative eigenvalue).

Under the assumptions of Theorem 14, let and consider such that . Then the scheme (20) satisfies


Suppose without loss of generality that the coordinates are chosen so that the eigenvector corresponding to the smallest eigenvalue

is . Then the second directional derivative of in the direction can be expressed as

Let us consider in particular the first angle in the angular discretization. Since , the second directional derivative of in this direction satisfies

Since is strictly negative, it is certainly the case that for sufficiently close to and small enough :

Now we perform a crude estimate on the sum in (20) by considering only a single term:

Using the same estimates on the discrete second directional derivatives, we can also estimate the term

By combining these estimates and recalling that , we conclude that

Lemma 17 (Consistency with a vanishing eigenvalue).

Under the assumptions of Theorem 14, let and consider such that . Then the scheme (20) satisfies


By the regularity of , we know that

Suppose without loss of generality that the coordinates are chosen so that the eigenvector corresponding to the smallest eigenvalue is . Then the second directional derivative of in the direction can be expressed as

Now we are going to estimate the sum in (20) by considering only the angles that are close to zero, which corresponds to the direction of the eigenvector . To this end, we define

and let be the number of nodes in the interval .

We notice that for any we have

Then using the lower bound allows us to obtain the following bounds on the sum.

This allows us to obtain bounds on the following value appearing in the scheme:

Recalling from the definition of that and allows us to simplify this as follows:

Thus under the conditions of Theorem 14, we find that

We also observe that

which implies that

We conclude that

which coincides with the value of

As an immediate consequence of the consistency proof (in particular, Lemma 15), we obtain the formal truncation error of the scheme as points where the function is “locally” strictly convex. This will be used to inform and optimize the particular implementations of this method in section 5.

Corollary 18 (Truncation error).

Under the assumptions of Theorem 14, let and consider such that . Then there exists a constant such that for all sufficiently small ,

We are also interested in approximating functions that are convex, but not necessarily strictly convex. In this case, the integrand in (17) is singular and we cannot directly use the formal truncation error of the quadrature rule. However, we can easily bound the resulting sums directly in the case where at least one eigenvalue is known to vanish. We consider two separate cases: (1) the fully degenerate case () and (2) the semi-degenerate case ().

Lemma 19 (Truncation error (fully degenerate)).

Under the assumptions of Theorem 14, let and consider such that . Then there exists a constant such that for all sufficiently small ,

Lemma 20 (Truncation error (semi-degenerate)).

Under the assumptions of Theorem 14, let and consider such that . Then there exists a constant such that for all sufficiently small ,

4.3. Quadrature rules

In designing a scheme of the form (20), the choice of quadrature rule

is a key factor that will influence the overall cost and accuracy.

A simple choice is the trapezoid rule, which utilizes the weights


As required, the weights are all positive. As required by (N8) (Definition 12), they can also be bounded from below in terms of the quasi-uniformity constant via

In general, the truncation error of the trapezoid rule is . However, in the special case of a uniform angular discretization ( for all ), the trapezoid rule is spectrally accurate. In this case, the truncation error satisfies for every given a integrand.

Higher-order quadrature is also possible. In fact, as we will demonstrate in section 5, this can be exploited in order to design approximation schemes that simultaneously improve the formal consistency error and reduce the required stencil width.

As an example, we consider Simpson’s rule. Suppose that , the number of angles in the angular discretization, is even. Then Simpson’s rule takes the form


where we identify and because of the periodicity of .

Rearranging, we find that the corresponding quadrature weights are


The truncation error associated with Simpson’s rule is .

However, unlike with the trapezoid rule, these quadrature weights are not automatically positive. Instead, positivity is guaranteed only if the quasi-uniformity constant of the angular discretization is not too large. In particular, we note that is sufficient to guarantee that

Under the same assumption on quasi-uniformity, we use the fact that

to verify that

as required by (N8) (Definition 12).

Similar results can be obtained using other higher-order quadrature schemes, which will place differing requirements on the quasi-uniformity constant in order to ensure positivity of the weights.

5. Implementation

We now present two specific implementations of a quadrature scheme based upon the formulation of (20).

The first implementation relies on a hexagonal tiling of the domain. The underlying structure of the grid allows us to design an angular discretization of involving twelve equi-spaced angles. The use of the trapezoid rule then leads to a compact finite difference stencil that in practice achieves spectral convergence in the angular parameter