1. Introduction
In this article we introduce an integral representation of the MongeAmpère equation
(1) 
where is convex and the righthand side is nonnegative. 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 MongeAmpè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 MongeAmpè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 MongeAmpère equation. Recently, convergence frameworks have been established for the MongeAmpère equation with either Dirichlet boundary conditions [14, 20, 26, 27]:
(2) 
or the second type boundary condition arising in optimal transport [4, 5, 21]:
(3) 
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 MongeAmpère equation [5, 16, 14, 24, 27]. Because of the widestencil nature of these methods, the methods are computationally expensive and typically have low (sublinear) 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 MongeAmpère operator in terms of a Gaussian integral. This allows us to utilize higherorder 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 MongeAmpè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 MongeAmpère equation holds particular promise for the development of computationally practical methods in three dimensions, as it provides a dimensionreduction as compared with a typical variational formulation of the 3D MongeAmpère equation. It also extends naturally to more general MongeAmpère type equations in optimal transport, including equations that are posed on the sphere [19].
2. Background
2.1. Elliptic equations
The MongeAmpère equation is an example of a degenerate elliptic partial differential equation, which take the general form
(4) 
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
(5) 
The MongeAmpè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 handinhand with this difficulty is the fact that the solution of the MongeAmpè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 MongeAmpère equation that automatically enforces solution convexity [20]. This can be accomplished by considering the convexified MongeAmpère operator
(6) 
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
(7) 
where
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 MongeAmpè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 semicontinuous envelopes of the candidate weak solution.
Definition 2 (Semicontinuous envelopes).
Let be a bounded function. Then for , the upper and lower semicontinuous envelopes are defined, respectively, as
Definition 3 (Viscosity subsolutions (supersolutions)).
A bounded upper (lower) semicontinuous 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 semicontinuous subsolution and a lower semicontinuous supersolution, then on
We remark that the Dirichlet problem for the MongeAmpère equation (1),(5) does satisfy a comparison principle under reasonable assumptions on the data. However, this is no longer true when the righthand 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
(8) 
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).
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]).
This result does apply to the MongeAmpère equation (1) with Dirichlet boundary conditions (5) under mild assumptions on the data. However, many other MongeAmpè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 nonclassical 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 MongeAmpère operator [2, 5, 8, 14, 26]. These hinge upon different reformulations of the MongeAmpère operator, which typically take a variational form
(11) 
Here denotes the second directional derivative (SDD) in the direction , is some admissible set, and is a nondecreasing 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
(12) 
These approximations are typically allowed to have a widestencil 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 equispaced 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].
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 nontrivial 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:
(13) 
where is a symmetric positivedefinite matrix.
This provides an alternate characterization of the MongeAmpè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
(14) 
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 MongeAmpère operator in can be expressed as
(15) 
Integrating out the radius , we obtain
(16) 
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 twodimensional version of (16) in terms of polar coordinates.
Theorem 11 (Integral Representation).
Let be convex and be strictly convex. Then for every :
(17) 
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:
(18) 
This, in turn, is used to construct a relaxed version of the convexified MongeAmpère operator:
(19) 
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 twodimensional MongeAmpère equation. In section 5, we will describe two particular implementations.
4.1. Approximation scheme
We introduce the following notation.
Definition 12 (Notation).

is a bounded, open, convex domain with Lipschitz boundary .

is a finite set of discretization points , .

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

is a stencil width associated to the grid.

is a finite set of angles discretizing .

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

is the angular resolution of the discretization.

is the quasiuniformity constant of the angular discretization.

is a collection of nonnegative quadrature weights summing to and satisfying
for some constant that depends only on the quasiuniformity constant.

is a regularization parameter associated with the grid.

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

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

is the maximal truncation error of the finite difference approximations.

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 MongeAmpère operator at interior points .
(20) 
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.
Proof.
We note that the operator that appears in (20) can be written in the form
where and the are nonnegative by the (negative) monotonicity of the approximations .
Let have nonnegative components. We notice that
Since the max and min operators preserve monotonicity and the weights are nonnegative, we can immediately conclude that
Since has no dependence on its second argument, this completes the proof of monotonicity. ∎
Theorem 14 (Consistency).
We will break this result into three separate cases (Lemmas 1517), 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).
Proof.
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 MongeAmpè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).
Proof.
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 asLet 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).
Proof.
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 semidegenerate 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 (semidegenerate)).
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
(21) 
As required, the weights are all positive. As required by (N8) (Definition 12), they can also be bounded from below in terms of the quasiuniformity 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.
Higherorder 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
(22) 
where we identify and because of the periodicity of .
Rearranging, we find that the corresponding quadrature weights are
(23) 
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 quasiuniformity constant of the angular discretization is not too large. In particular, we note that is sufficient to guarantee that
Under the same assumption on quasiuniformity, we use the fact that
to verify that
as required by (N8) (Definition 12).
Similar results can be obtained using other higherorder quadrature schemes, which will place differing requirements on the quasiuniformity 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 equispaced 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