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 , geophysics , mesh generation , medical image processing , meteorology 
, and data science. 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  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]:
These convergence proofs can be viewed as extensions of the powerful Barles and Souganidis convergence framework 
, 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 .
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 .
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 . 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 , 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 .
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 ).
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  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 .
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 . A typical scaling for the maximal stencil width that optimizes truncation error is at least , 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).
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 quasi-uniformity constant of the angular discretization.
is a collection of non-negative quadrature weights summing to and satisfying
for some constant that depends only on the quasi-uniformity constant.
is a regularization parameter associated with the grid.
is the set of neighboring indices for such that for every we have .
described in (12) has the form
for every , where all .
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 Monge-Ampère operator at interior points .
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).
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).
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).
Suppose without loss of generality that the coordinates are chosen so that the eigenvector corresponding to the smallest eigenvalueis . 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).
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).
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)).
Lemma 20 (Truncation error (semi-degenerate)).
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.
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