MONOTONE DISCRETIZATION OF THE MONGE–AMPÈRE EQUATION OF OPTIMAL TRANSPORT

We design a monotone finite difference discretization of the second boundary value problem for the Monge–Ampère equation, whose main application is optimal transport. We prove the existence of solutions to a class of monotone numerical schemes for degenerate elliptic equations whose sets of solutions are stable by addition of a constant, and we show that the scheme that we introduce for the Monge–Ampère equation belongs to this class. We prove the convergence of this scheme, although only in the setting of quadratic optimal transport. The scheme is based on a reformulation of the Monge–Ampère operator as a maximum of semilinear operators. In dimension two, we recommend to use Selling’s formula, a tool originating from low-dimensional lattice geometry, in order to choose the parameters of the discretization. We show that this approach yields a closed-form formula for the maximum that appears in the discretized operator, which allows the scheme to be solved particularly efficiently. We present some numerical results that we obtained by applying the scheme to quadratic optimal transport problems as well as to the far field refractor problem in nonimaging optics. Mathematics Subject Classification. 65N06, 65N12, 35J70, 35J96. Received June 10, 2021. Accepted March 14, 2022.

the considered optimal transport problem, as suggested in this paper and previously in [3,26]. One benefit of this last approach is that, as illustrated by our numerical experiments in the setting of nonimaging optics in Section 6.5, it can be applied to optimal transport problems with various cost functions, provided that one uses a numerical scheme capable of handling Monge-Ampère equations with arbitrary coefficients. Another benefit is that convergence results may be established using the theory of monotone schemes for degenerate elliptic partial differential equations [1]. Note however that establishing theoretical guarantees is more complicated when considering general cost functions, and in this paper our convergence proof only applies in the setting of a quadratic transport cost.
We design a monotone finite difference discretization of the Monge-Ampère equation where is an open bounded subset of R containing the origin and and are bounded functions on ×R , whose values are respectively symmetric matrices and nonnegative numbers, and 1/ being Lipschitz continuous with respect to their second variables uniformly with respect to their first variables, and being continuous with respect to both its variables. For any symmetric matrix of size , we denoted (We use the Loewner order on the space of symmetric matrices: 1 ⪰ 2 if 1 − 2 is positive semidefinite. From now on, we denote respectively by , + , and ++ the sets of symmetric, symmetric positive semidefinite, and symmetric positive definite matrices of size .) Since we consider Monge-Ampère equations which are related to the problem of optimal transport, see Since we consider Section 5.1 and Remark 5.1, we also have to discretize the relevant boundary condition, described in Section 1.2. We prove the existence of solutions, under suitable assumptions, to the proposed finite difference scheme. We also prove the convergence of solutions to this scheme, but only in the setting of quadratic optimal transport, where the function is identically zero and the function is separable in the form ( , ) = ( )/ ( ). The Monge-Ampère equation is degenerate elliptic, meaning that it may be written in the form MA ( , ( ), 2 ( )) = 0 in , (1.2) where the operator MA : × R × → R is degenerate elliptic, that is, nondecreasing with respect to its last variable: MA ( , , 1 ) ≤ MA ( , , 2 ) if 1 ⪰ 2 . The degenerate ellipticity property has a discrete counterpart which we call monotonicity, see Definition 2.5. Convergence of monotone schemes for degenerate elliptic equations may often be proved using a general argument, which was introduced in [1]. We use the fundamental part of this argument, see Theorem 2.7. As we discuss below Theorem 2.7, the full convergence result stated in [1] requires the approximated equation to satisfy a strong comparison principle which does not hold for the Monge-Ampère equation equipped with the boundary condition (1.24). Therefore, in order to prove Theorem 5.25, our convergence result in the setting of quadratic optimal transport, we need to establish an appropriate substitute to this comparison principle, in the form of Theorems 5.11 and 5.12.
One way to define the operator MA ( , , ) so that it is both degenerate elliptic and consistent with (1.1) would be as ( , ) − det + ( − ( , )). (1.3) This is not the definition we use, however. The reason is that there is no obvious way to build a monotone scheme by directly discretizing (1.3). Instead, we use strategies described in [35,36] ⟩︀ )︃ = 0 (1. 6) and alternatively, following [23], where for any symmetric matrices and and nonnegative number , Note that the maximum in (1.7) is attained, as the maximum over a compact set of the continuous function ↦ → ( , ) (this function is also concave, by the Minkowski determinant inequality). On the contrary, the parameter set of the infimum in (1.6) is not compact. Both reformulations enforce the admissibility constraint (1.4): for instance in (1.7), for any unit vector ∈ R , choosing = ⊗ in the maximum yields the inequality

Discretization of the Monge-Ampère equation
For any discretization step ℎ > 0, we discretize the operator MA on a grid ℎ ⊂ ∩ ℎZ . Denoting by H the Hausdorff distance between compact subsets of R , which we recall is defined by we will assume that lim 10) or equivalently that if ⊂ is compact, then for sufficiently small ℎ > 0 one has ∩ ℎZ 2 ⊂ ℎ . We will also need the technical assumption (3.1) of uniform connectedness of the grid ℎ .
Before introducing the discretization of MA , we need to define some finite difference operators. For any function : ℎ → R, point ∈ ℎ , and vector ∈ Z , we define The constant +∞ in the definition of ℎ is related to the way we recommend discretizing the optimal transport boundary condition, discussed in Section 1.2.
In the whole paper, we denote by ( 1 , . . . , ) the canonical basis of Z . For any function : ℎ → R and point ∈ ℎ , we define the Laplacian approximation and, whenever it makes sense, the centered gradient approximation (1.13) We use Lax-Friedrichs approximations of the gradient of in ( , ( )) and ( , ( )). To this end, we let min ≤ 0, LF ≥ 0, and LF ≥ 0 be three constants independent of ℎ. We will assume that for any ∈ and , ′ ∈ R , ( , ) ⪰ min , (1.14) (1. 16) For any function : ℎ → R, point ∈ ℎ , and vector ∈ Z , we define (1.18) (In the whole paper, we denote respectively by ∨ and ∧ the maximum and the minimum of two real numbers and .) For any family = ( ) 1≤ ≤ of vectors of Z and any ∈ R , we define Finally, for any function : ℎ → R, point ∈ ℎ , and family of vectors of Z , we define Equivalently, if ⊂ ++ is compact, then for sufficiently small ℎ > 0 each element of can be written as ( ) where ∈ ℎ and ∈ R ( +1)/2 + . We will also need to assume that lim ℎ→0 ℎ max ∈ ℎ max ∈ | | = 0, (1.21) and that for any ℎ > 0, (1.22) where we recall that 1 denotes the first vector of the canonical basis of R . We discretize MA by the operator Coefficients of are required to be nonnegative in order for the discretization to result in a numerical scheme which satisfies the monotonicity property (defined rigorously in Def. 2.12). Note that the constraint Tr( ( )) = 1 may be rewritten as ∑︀ ( +1)/2 =1 | | 2 = 1. In dimension = 2, we recommend choosing ℎ as a set of superbases of Z 2 : Note that in the definition above, the constraint det( 1 , 2 ) = ±1 is equivalent to det( 2 , 3 ) = ±1 or det( 1 , 3 ) = ±1. We explain in Appendix B how a set ℎ of superbases of Z 2 satisfying the above assumptions may be constructed, using tools from the fields of lattice geometry and arithmetic known as the Selling's decomposition [43] and the Stern-Brocot tree [7]. We prove in Section 4 that when choosing ℎ in this way, the second maximum in (1.23) admits a closed-form expression, at least when no infinite values are involved (infinite values may stem from the handling of the boundary condition, see (1.11), and a simple modification of the formula of Thm. 1.2 allows to compute the maximum in this case, by excluding finite differences whose value is infinite): Theorem 1.2. If = ( 1 , 2 ) is a basis of Z 2 , then for any ≥ 0 and ∈ R 2 , where for any ∈ , ( ) is an open bounded convex nonempty subset of R . We assume that ( ) depends continuously on , for the Hausdorff distance H over compact subsets of R whose definition we recalled in (1.9).
In the particular setting of quadratic optimal transport, in which we will prove convergence of the proposed numerical scheme, the set ( ) does not depend on the variable . Note that despite being called a boundary condition, the constraint (1.24) involves the whole domain . Some numerical approaches for solving the second boundary value problem, although not the one that we describe in this paper, rely on the fact that, in some cases, the constraint (1.24) can be reformulated in a way that only involves the boundary of the domain , see for instance [3]. For now, let us consider the class of numerical schemes for equations (1.1) and (1.24) that are defined, for any discretization step ℎ > 0, by an operator ℎ MABV2 : R ℎ → R ℎ , and may be written as One property of equations (1.1) and (1.24) is that their expressions depend only on derivatives of the function and not on itself, and therefore that the set of solutions is stable by addition of a constant. Accordingly, we say that the operator ℎ MABV2 and the scheme (1.25) are additively invariant if for any function : ℎ → R and real number , ℎ MABV2 ( + ) = ℎ MABV2 . We adapt the approach introduced in [26] to build an operator ℎ MABV2 suitable for (1.25). The idea is to build ℎ MABV2 as a maximum of ℎ MA and of a monotone discretization ℎ BV2 : R ℎ → R ℎ of the left-hand side in a degenerate elliptic formulation of (1.24).
We use the following formulation of (1.24), initially introduced in [4]: (1.27) (We denote by ( ) the support function of the convex set ( ): for any ∈ R , ( ) ( ) := sup ∈ ( ) ⟨ , ⟩. Formally, if belongs to the boundary ( ) of ( ), then the maximum in the definition of BV2 is attained when is the unit outer normal of ( ) at the point .) For any function : ℎ → R, point ∈ ℎ , and vector ∈ R , we define the upwind finite difference In this setting, the scheme (1.25) is additively invariant. Additively invariant schemes of the form (1.25) are not well-posed: their sets of solutions are stable by addition of a constant, thus not a singleton. Moreover they often have no solutions. One way to see this formally is that a well-posed scheme would need an additional equality to guarantee uniqueness of solutions, for instance [0] = 0, but that then there would be one more equality than unknowns in the scheme. In the continuous setting, equations whose sets of solutions are stable by addition of a constant often admit solutions if and only if their coefficients satisfy some nonlocal condition, such as the mass balance condition (5.1) in the case of the Monge-Ampère equation of optimal transport; however, there may be no obvious discrete counterpart to this condition. See Section 2 for further discussion of this issue.
In order to get around this difficulty, we solve an altered form of the scheme (1.25), following the approach used in the numerical experiments of [3] (note that we present a fully detailed mathematical analysis of this alteration, in contrast with [3] where it was introduced essentially as a numerical trick). We add an unknown to the scheme, which must be a real number. For fixed , we define the operators ℎ, MA : R ℎ → R ℎ and ℎ, (1.29) The scheme that we actually solve, with respect to the extended unknown ( , ), is ℎ, (1.30)

Main contributions and relation to previous works
We introduce the numerical scheme (1.30) for the Monge-Ampère equation (1.1), equipped with the boundary condition (1.24). We prove the existence of solutions to a class of monotone additively invariant numerical schemes featuring an additional unknown ∈ R as in (1.30), see Section 2, and we show, in Section 3, that the scheme (1.30) belongs to this class. This scheme is based on a discretization of the reformulation (1.7) of the Monge-Ampère equation. We prove in Section 4 that this discretization admits a closed-form expression, as stated in Theorem 1.2. We prove convergence of the scheme in the setting of quadratic optimal transport, see Section 5; convergence in the setting of more general optimal transport problems remains an open problem. We present in Section 6 some numerical experiments, including an application to the far field refractor problem in nonimaging optics.
The closed-form expression obtained in Theorem 1.2 makes the implementation of the scheme particularly efficient, since no discretization of the parameter set of the maximum in (1.7) is needed. While to our knowledge the proposed discretization is the first one to admit such a closed-form expression among those that are based on the reformulation (1.7) of the Monge-Ampère equation, it is to be related to the MA-LBR scheme, introduced in [5] in the setting of the Dirichlet problem for the Monge-Ampère equation when the function is identically zero, and to the scheme we introduced in [9] for the Pucci equation. Both of the above-mentioned schemes involve the notion of superbases of Z 2 . We prove in Appendix A that the MA-LBR scheme is a discretization of (1.6), although it was not introduced as such in [5].
As opposed to (1.6), the reformulation (1.7) has the benefit that its left-hand side remains finite even when (1.4) is not satisfied. Thus schemes based on it are more stable numerically than those based on (1.6), and can handle the degenerate case of functions : × R → R + which are not everywhere positive, in which case the solutions to the Monge-Ampère problem typically satisfy (1.4) but not its strict variant. On the contrary, the MA-LBR scheme only applies in the case > 0, and in addition solving it using the damped Newton requires using extremely small steps so that the constraint (1.4) remains satisfied along the iterations. This behavior of the Newton method is illustrated numerically in Section 6.4, and does not occur with the scheme introduced in this paper.
Numerical schemes based on (1.7) were previously introduced in [23], and then in [14], although only in the setting of the Dirichlet problem for the Monge-Ampère equation when = 0. In those papers, no counterpart of Theorem 1.2 was proved, hence the parameter set of the maximum in (1.7) had to be discretized.
Convergence of schemes for the second boundary value problem was previously studied in [3,26] in the setting of the quadratic optimal transport problem. Schemes considered in those two papers were based on the MA-LBR scheme introduced in [5], and adapted in order to discretize the boundary condition (1.24).
In [3], convergence of a scheme of the form (1.25) was proved, but existence of solutions to this scheme was not. It turns out that solutions typically do not exist, due to the scheme being additively invariant. The approach used to solve the scheme in the numerical experiments was equivalent to adding an unknown ∈ R as in (1.30), but the proof of convergence was not extended to this setting. Remark 1.3 (Applicability of Thm. 2.15 to the scheme in [3]). The work [3] establishes the convergence of the solutions to a discretization of the optimal transport problem, under the assumption that they exist. The latter point is dubious, as discussed above and acknowledged by the authors of [3], and for that reason an altered variant of the scheme is considered in the numerical experiments section, featuring an additional unknown which is analogous to the parameter ∈ R in (1.29); the existence of solutions to this variant is observed numerically in [3], but left as an open problem from a theoretical standpoint.
While the detailed analysis of the scheme in [3] is out of the scope of this paper, let us discuss the applicability to this scheme of the assumptions of our existence result, Theorem 2.15. Those assumptions are continuity, monotonicity and stability, in the sense of Definition 2.12. The continuity of the scheme in [3] is easy to prove, since no infinite values are involved in the definition of the scheme operator, contrary to the scheme that we introduce in this paper. The monotonicity property is not satisfied by the scheme recommended by default in [3] due to the centered discretization of the gradient of the unknown ; however, a monotone Lax-Friedrichs discretization was described as an alternative in Section 4.4 of [3]. The remaining open question is the stability of the scheme in [3], again in the sense of Definition 2.12; we could not see an immediate proof of this property, but one could expect to develop one based on the sketches of the proofs of Proposition 3.6 and also of Proposition 4.3, item (5) from [3], which is a result about the Lipschitz continuity of the solutions to this scheme.
Another scheme of the form (1.25) is studied in [26]. In that work, a Dirichlet boundary condition is enforced on , which in our setting would amount to replacing +∞ with some fixed constant ∈ R in (1.11). Therefore the scheme considered in [26] is not additively invariant. The Dirichlet boundary condition is to be understood in a weak sense (the one of viscosity solutions, see Def. 2.3). It may formally be simplified to ( ) ≤ on , with equality at some point * ∈ . An important assumption in [26] is that the scheme satisfies a property of underestimation, discussed in Remark 1.4. Under this assumption, the existence and convergence of solutions is proved. The property of underestimation is satisfied in the case of quadratic optimal transport at the cost of a careful handling of the constraint (1.24), but it does not seem obvious that it is satisfied for similar schemes in the case of more general optimal transport problems, with ̸ = 0 in (1.1). No numerical experiments were performed in [26]. In our experience, the scheme introduced in that paper has the drawback that the numerical error of its solutions tends to be unevenly distributed. This effect is related to the particular role played in the discretization by the point * ∈ where the Dirichlet condition is satisfied in the classical sense, which leads to numerical artifacts and tends to decrease the accuracy of the scheme.
In our proof of convergence of the scheme (1.30), we use the arguments introduced in [26] when appropriate. However, the property of underestimation is not required in our setting. Remark 1.4 (Role of the property of underestimation in [3,26], and substitutes used in this paper). In each of the papers [3,26], the theoretical analysis of the considered scheme uses the fact that this scheme satisfies some property of underestimation. The schemes considered in those papers are both based on the MA-LBR discretization [5], represented in Appendix A by the operator Λ ℎ in (A.1), of the Monge-Ampère operator ↦ → det 2 (·). The property of underestimation used in [3] is formulated as the fact that the operator Λ ℎ overestimates the Lebesgue measure of the subgradient, in the sense that Λ ℎ [ ] ≥ ℎ − |˜( )| for any function : ℎ → R and for any suitable point ∈ ℎ , where˜: R → R denotes the convex envelope of . This property is described in more detail in Lemma 4.2 of [3]. Whether this property could be extended to the setting of the scheme (1.30) considered in this paper is not clear, since this scheme is based on a discretization of the reformulation (1.2) of the Monge-Ampère equation, which does not feature directly the Monge-Ampère operator ↦ → det 2 (·). The arguments in the convergence analysis in [3] that use the property of underestimation are based on the construction of solutions to semi-discrete optimal transport problems, and are completely different from the arguments in this paper, which are based on the theory of convergence of monotone schemes to viscosity solutions to degenerate elliptic equations.
In [26], two distinct definitions of the property of underestimation are given. The first one ( [26], Def. 3.8) is similar to the definition in [3]. The second one ( [26], Rem. 3.9) asks that ℎ MABV2 [ ] ≤ MABV2 ( ) in ℎ , for all smooth convex functions , where ℎ MABV2 is the discrete operator describing the whole scheme and MABV2 is its continuous counterpart. The second definition is claimed to correctly approximate the first one at small grid scales.
The property of underestimation is not only used in the convergence analysis in [26], but it is also crucial for the proof ( [26], Lem. 3.12), which guarantees the existence of a subsolution to the scheme in [26] and is an intermediary step for proving the existence of solutions. One technique that is often used to build a subsolution to a scheme is to consider a strict subsolution to the continuous problem, and to use the consistency of the scheme to show that it is also a subsolution to the scheme. However in the setting of [26] no strict subsolutions to the continuous problem may exist, as shown by some counterpart to Theorem 5.11 in this setting (with = 0, see [26], Thm. 2.1). Formally, the property of underestimation allows one to consider a solution to the continuous problem instead of a strict subsolution in the argument above (in [26], a slight variation of the solution, constructed by solving a semi-discrete optimal transport problem, is considered instead for technical reasons). In our setting, Theorem 5.11 does not prevent us to build a strict subsolution to the continuous problem for some value < 0, which is sufficient in order to apply our existence result Theorem 2.15, hence we have no need for the underestimation property at this stage.
Let us finally discuss why the property of underestimation is needed in the proof of the convergence result ( [26], Thm. 3.11). In the standard theory of convergence of monotone schemes [1], it is shown that some appropriately defined upper and lower limits and of sequences of solutions to the scheme are respectively a subsolution and supersolution to the continuous problem (see Thm. 2.7). In [26], the counterpart ( [26], Thm. 2.1) to Theorem 5.11 is used to deduce that the subsolution is actually a solution to the Monge-Ampère problem. From this point, it remains to prove that = . In the general setting of monotone schemes, this is often done by using a comparison principle for the continuous problem, but such a comparison principle is lacking in the setting of [26]. Another strategy is to prove that the solutions to the scheme are sufficiently regular so that the limits and coincide by definition, at least up to extraction of a subsequence: this is what we do in this paper, see our stability result Proposition 3.6. In [26], no such stability result is proved, and the arguments used instead to show that = rely on the assumption that the scheme satisfies the property of underestimation, indirectly through the use of the subsolutions to the scheme built in the previous paragraph.
Note that the scheme (1.30), and its continuous counterpart (3.2) below, which both feature an additional unknown or parameter ∈ R, fit in the framework of eigenvalue problems recently studied in [27]. Although our proof of convergence only applies to Monge-Ampère equation in the setting of quadratic optimal transport, our existence result, Theorem 2.15, is applicable to other such eigenvalue problems, as illustrated by the examples in Section 2.

Degenerate elliptic additively invariant equations
In this section, we study numerical schemes for a general degenerate elliptic equation of the form Typically, is discontinuous and ( , , ) is defined differently depending on whether belongs to or to , in order to take into account the boundary condition in equation (2.1). The equation without the boundary condition would then be (︀ , Let us recall the definition of degenerate ellipticity: We say that equations (2.1) and (2.2) are additively invariant since, for reasonable notions of solutions, their sets of solutions are stable by addition of a constant, due to the fact that at any point , the left-hand sides of those equations depend only on the derivatives ( ) and 2 ( ) of the function , and not on its value ( ). This is not a standard property, and we will show that it is a source of difficulty in the design of monotone numerical schemes. Typically, an additively invariant equation only has solutions if its coefficients are well-chosen and satisfy a particular nonlocal property. where the degenerate elliptic operator ex : (The choice ex ( , , ) = ( ) − , rather than ex ( , , ) = − ( ), is required in order for the equation (2.3) to be degenerate elliptic. The choice ex (−1, , ) = − and ex (1, , ) = , rather than ex (−1, , ) = and ex (1, , ) = − , is not dictated by the degenerate ellipticity property, but is the standard formulation of Neumann boundary conditions for degenerate elliptic equations, and the one which allows a comparison principle to hold in the context of more favorable equations such as − ′′ + = 0, see Section 7.B of [16].) The equation (2.3) only has solutions (respectively subsolutions, supersolutions) if ∫︀ 1 −1 ( ) d = 0 (respectively ≤ 0, ≥ 0), which we assume. Notice the similarity with the mass balance condition (5.1) which occurs in the setting of optimal transport. An appropriate notion of solutions for degenerate elliptic equations, and for the study of discretizations of such equations, is the one of viscosity solutions. Before defining them, let us recall the definitions of the upper semicontinuous envelope * and lower semicontinuous envelope * of a function : → R, being a subset of R : for any ∈ ,

Definition 2.3 (Viscosity solution). A function :
→ R is a viscosity subsolution to (2.1) if (i) it is upper semicontinuous and (ii) for any function in 2 ( ) and local maximum of − in , It is a viscosity supersolution if (i) it is lower semicontinuous and (ii) for any function in 2 ( ) and local minimum of − in , It is a viscosity solution if it is both a viscosity subsolution and a viscosity supersolution. The same definitions, with replaced by , apply to equation (2.2).
Note that if a viscosity subsolution (respectively supersolution) to (2.1) is twice differentiable at some point , then is a classical subsolution (respectively supersolution) to (2.1) at the point .

Discretization
For any discretization step ℎ > 0, let ℎ be a finite nonempty subset of containing the origin. In the rest of this paper, it is required that ℎ be a subset of the Cartesian grid ∩ ℎZ ; however, this is not necessary in this section. What will be required in our definition of consistency is that Note that in the case that ℎ is included in ∩ ℎZ , then (2.4) is implied by (1.10). We represent discretizations of the operator by operators : R ℎ → R ℎ that are additively invariant, according to the following definition: For now, we let ℎ : R ℎ → R ℎ be an additively invariant operator, for any ℎ > 0, and we consider a numerical scheme of the form  A framework is outlined in [1] for the proof of convergence of monotone schemes. The following fundamental result follows directly from the proof of Theorem 2.1 from [1]: Theorem 2.7. Assume that there exist a sequence (ℎ ) ∈N of discretization steps ℎ > 0 converging to zero and a sequence ( ) ∈N of solutions : ℎ → R to (2.5) with ℎ = ℎ such that [ ] is bounded, uniformly over ∈ N and ∈ ℎ . If (2.5) is monotone and consistent with equation (2.1), then functions , : → R defined by are respectively a viscosity subsolution and supersolution to (2.1).
The definition of consistency in Definition 2.5 is slightly simpler than the one in [1], due to the assumption that operators ℎ are additively invariant. In the framework of [1], in which the left-hand side in (2.1) may also depend on ( ), a strong comparison principle, that is, a result stating that viscosity subsolutions to (2.1) are always less than viscosity supersolutions, is used after applying Theorem 2.7 to prove that ≤ , which allows to conclude that = , since ≥ by definition. Obviously, no strong comparison principle may hold if the set of viscosity solutions is nonempty and stable by addition of a constant. In our proof of convergence in the setting of quadratic optimal transport, we use Theorems 5.11 and 5.12 as a substitute to this comparison principle.
An important difficulty that we encounter is that numerical schemes of the form (2.5) typically have no solutions.
Then the scheme ℎ ex [ ] = 0 in ℎ is monotone and consistent with equation (2.3). Solving this scheme is equivalent to solving a square linear system, since the scheme operator ℎ ex : R ℎ → R ℎ is an affine map. However, this linear system is noninvertible, since all constant functions belong to the kernel of the associated linear operator.
To get around this difficulty, we add a parameter ∈ R to the equation (2.1), yielding a new equation where for any ∈ R, : × R × → R is a given operator, typically degenerate elliptic. The idea is to choose so that 0 = and (2.7) has no viscosity subsolutions when > 0 and no viscosity supersolutions when < 0. Example 2.9. For any ∈ R, we define ex : Then equation ex ( , ′ ( ), ′′ ( )) = 0 in coincides with (2.3) when = 0, and only has solutions (respectively subsolutions, supersolutions) if . Recall that we assumed that Accordingly, we add an unknown ∈ R to the numerical scheme. For any ℎ > 0 and ∈ R, we let ℎ, : R ℎ → R ℎ be an additively invariant operator, and we consider the scheme Example 2.10. In the setting of Example 2.8, for any ℎ > 0 and ∈ R, we define ℎ, [ ] = 0 in ℎ may easily be constructed explicitly.
The definition of solutions ( , ) ∈ R × R ℎ to (2.8) is obvious, but we will also need a notion of subsolutions (we could define supersolutions similarly, but this will not be needed): Since is an unknown of the scheme, and not simply a fixed parameter, Definition 2.5 needs to be adapted to this new setting. We also define some other properties that the scheme (2.8) may satisfy. Conceptually, the following definition is intended for schemes such that ℎ, [ ] is nondecreasing with respect to . -Monotone if for any ∈ R, the scheme (2.5) with ℎ = ℎ, is monotone in the sense of Definition 2.5.
-Consistent with the parametrized equation (2.7) if for any family of real numbers ( ℎ ) ℎ>0 converging to some ∈ R as ℎ approaches zero, the scheme (2.5) with ℎ = ℎ, ℎ is consistent with equation (2.7) in the sense of Definition 2.5.
takes finite values and is continuous.
(ii) There exists a nonincreasing function : R → R + such that for any small ℎ > 0, any subsolution ( , ) ∈ R × R ℎ to (2.8), and any 1 , 2 ∈ ℎ , one has (iii) There exists 0 ∈ R such that for any small ℎ > 0 and any subsolution ( , ) ∈ R × R ℎ to (2.8), one has ≤ 0 . (iv) There exists 1 ∈ R such that for any small ℎ > 0 and any solution ( , ) ∈ R × R ℎ to (2.8), one has ≥ 1 . -Equicontinuously stable if it satisfies all items in the definition of stability above, with (ii) replaced by the following: (ii') There exists a function : R × R + → R + , nonincreasing with respect to its first variable and satisfying lim →0 ( , ) = 0 for any ∈ R, such that for any small ℎ > 0, any subsolution ( , ) ∈ R × R ℎ to (2.8), and any 1 , 2 ∈ ℎ , one has Note that the value = 0 does not play a special role in Definition 2.12. The role of the functions in the definitions of stability and equicontinuous stability is to allow schemes to become unstable when → −∞.
Obviously, if (2.8) is equicontinuously stable, then it is stable. In the case of the scheme considered in this paper for the Monge-Ampère equation, subsolutions will be established to be uniformly Lipschitz continuous, which is stronger than equicontinuity, see the proof of Proposition 3.6. In particular, the boundary condition ( ) − ∞ = 0 on (to be understood in the viscosity sense, as mentioned in Sect. 1) does not induce a boundary layer.
Theorem 2.7 is easily adapted to the scheme (2.8): Corollary 2.13. Assume that there exist a sequence (ℎ ) ∈N of discretization steps ℎ > 0 converging to zero, a sequence ( ) ∈N of real numbers converging to some ∈ R, and a sequence ( ) ∈N of functions [ ] is bounded, uniformly over ∈ N and ∈ ℎ . If (2.8) is monotone and consistent with (2.7), then limits superior and inferior , : → R defined as in (2.6) are respectively a viscosity subsolution and supersolution to (2.7) in .
If (2.8) is equicontinuously stable, then Corollary 2.13 is simplified by the fact that, by the Arzelà-Ascoli theorem, sequences ( ) ∈N and ( ) ∈N converge uniformly, up to extracting a subsequence, to some ∈ R, and to some continuous function : → R, which coincides with the limits superior and inferior and for this subsequence.

Existence
Our main result in this section concerns existence of solutions to the scheme (2.8). The proof is an adaptation of discrete Perron's method to our setting. Remark 2.14 (Discrete Perron's method). In the context of this remark, let us consider a scheme of the form (2.5) where the operator ℎ : R ℎ → R ℎ is not necessarily additively invariant. Discrete Perron's method states that, if this scheme is monotone in the sense of Definition 2.5 and if, at some fixed step size ℎ > 0, the map ℎ : R ℎ → R ℎ takes finite values and is continuous, then the function˜: ℎ → R defined bỹ is a solution to the scheme, provided that this function takes finite values. While the discrete version of Perron's method is not as extensively described in the literature as its continuous part ( [16], Sect. 4), some of its variants are stated and proved in Theorem 2.3 of [40] and Theorem 3.5 of [42]. Proof. We define the set of subsolutions to (2.8). Since we assumed that (2.8) is stable, is nonempty and there exists ∈ R defined by (2.10) Let us show that there exists : ℎ → R such that ( , ) is a subsolution to (2.8). Let (( , )) ∈N be a maximizing sequence in the definition of , and let * := min ∈N . We may assume, up to adding a constant to , that [0] = 0 for any ∈ N. Then by stability, for any ∈ N and ∈ R ℎ . This means that the sequence ( ) ∈N is bounded in R ℎ and thus that it converges, up to extracting a subsequence, to some functionˆ: ℎ → R. By continuity of the scheme, ( ,ˆ), as the limit of subsolutions (( , )) ∈N , is a subsolution to (2.8).
Among all functions : ℎ → R such that ( , ) is a subsolution to (2.8), we choose one which maximizes the cardinality of the set * : . We will show that such a function may be turned into another function˜: ℎ → R such that ( ,˜) is a solution to (2.8), by a maximization procedure similar to the construction of the function˜in Remark 2.14.
One important difference between our proof and the classical proof of Perron's method is that in the classical proof, no specific subsolution has to be used as a basis for constructing the solution˜. In particular the assumption about the maximal cardinality of the set * is specific to our setting.
Note that * cannot be equal to ℎ , since in this case, by continuity of the scheme, there would exist ′ > such that ( ′ , ) ∈ (choose ′ close enough to ), and this would contradict (2.10).
Knowing that * ̸ = ℎ , and using stability, we may define, for small > 0, the function˜: To ensure that the supremum above is the one of a nonempty set, we choose small enough so that itself is suitable choice of function . In the following two paragraphs, we show, using arguments from the classical proof of Perron's method, that ( ,˜) is a subsolution to (2.8) and that ℎ, By continuity of the scheme, we may pass to the limit in maximizing sequences and deduce that for any ∈ ℎ , there exists : ].
Note that the right-hand side is the limit of a bounded nondecreasing sequence. By continuity, ℎ,˜[ ] = 0 in * and ( ,˜) is a subsolution to (2.8). Let us show that it is a solution. If it is not the case, then there exists * ∈ ℎ ∖ * such that ℎ,˜[ * ] < 0. By continuity, there exists > 0 such that ℎ,˜[ * ] < 0. Since ( ,˜) is a subsolution to (2.8) and ℎ,˜[ ] < 0 in * , this contradicts the assumption that * is of maximal cardinal. Thus ( ,˜) is necessarily a solution to (2.8).
Remark 2.16. Since ℎ > 0 is fixed in Theorem 2.15, the subsolution, the function , and the number 0 in (i), (ii), and (iii) in the definition of stability of the scheme (Def. 2.12) only need to exist for this fixed value of ℎ. Also, (iv) is not needed.

Properties of the proposed scheme
In this section, we show that the scheme (1.30) satisfies the properties we defined in Section 2. First note that for any ℎ > 0 and ∈ R, the operator ℎ, MABV2 : R ℎ → R ℎ is additively invariant. Proof. Let ℎ > 0, ∈ R, ∈ ℎ , and , : ℎ → R be such that [ ] = [ ] and ≥ in ℎ . We need to show that ℎ, By the definition (1.29) of the operator ℎ, MABV2 , it suffices to prove that both ℎ The second inequality follows directly from the definition (1.28) of ℎ BV2 , so let us prove the first one.
By the definition (1.23) of ℎ MA , it suffices to prove that for any family = ( 1 , . . . , ) of vectors of Z and any ∈ R + , by definition of ℎ , implies that ± ℎ ∈ ℎ for any ∈ {1, . . . , }), then, using (1.16) for the second inequality, and otherwise, using (1.15), From the grid ℎ , we may build a graph whose nodes are the points of ℎ and whose edges are pairs of points that are neighbors on the grid, that is, between whom the Euclidean distance is equal to ℎ. To prove other properties of the scheme, we need the technical assumption that the distance on this graph, multiplied by ℎ, is equivalent to the Euclidean distance, uniformly over small ℎ > 0. Equivalently, we require that there exists some positive constant , such that for any small ℎ > 0 and any function : ℎ → R, which concludes the proof.
Let us now study the consistency of the scheme (1.30) with the degenerate elliptic equation where for any ∈ R, ∈ , ∈ R , and ∈ , and MA ( , , ) and BV2 ( , ) are defined respectively in (1.8) and (1.27). We first prove a consistency property that is stronger to the one we introduced in Definition 2.12, and that will be useful in the study of stability of the scheme. . Let ∈ ∞ ( ) and ( ℎ ) ℎ>0 be a family of real numbers converging to some ∈ R as ℎ approaches zero. Then uniformly over ∈ ℎ and ∈ R. Moreover, for any compact subset of ,
We deduce (3.5) using (1.20) and that the affine map is continuous, uniformly over and belonging to compact sets. For the operator ℎ MA , we distinguish two cases: (General case) If there exist 1 > 0 and 2 ∈ (0, 1) such that the following refinements of (1.20) and (1.21) hold: then, refining the last equality in (3.7) and using that the map (3.8) is 1/ -Hölder continuous, one has, for small ℎ > 0 and uniformly over ∈ ∩ ℎ , In dimension = 2, when choosing ℎ as in Remark B.9, one has 1 = 2 and 2 = , hence ℎ MA is consistent with MA to the order 1 ∧ (2 − 2 ) ∧ , and the optimal choice for is = 2/3, yielding consistency to the order 2/3. (Smooth case) The consistency is improved if (1.2) admits a solution ∈ 2 ( ) such that, uniformly over , 2 ( )− ( , ( )) ∈ ++ has condition number less than some constant > 1. In this setting, the maximum )︀ −1 )︁ , which has condition number less than for all ∈ . We thus recommend choosing the set ℎ independently of ℎ, but such that any ∈ ++ with condition number less than is of the form = ( ) for some ∈ ℎ and ∈ R ( +1)/2 + (see Appendix B for a suitable construction of ℎ in dimension = 2). Then (1.20) is not satisfied, but in a neighborhood of the solution , the operator ℎ MA is still consistent with MA , to the order one, uniformly over ∈ . In practice, one may choose to implement the scheme with Lax-Friedrichs relaxation parameters LF = LF = 0, as we do in Section 6. The drawback of doing this is that (1.15) and (1.16), and thus Proposition 3.1, do not hold anymore unless ( , ) and ( , ) do not depend on . The benefit is that consistency is improved. In the setting of the smooth case described above, if LF = LF = 0, then, in a neighborhood of and uniformly over ∈ , ℎ MA is consistent with MA to the order two. Note that the order of consistency of the whole scheme (1.30) is the minimum of the ones of ℎ BV2 and ℎ MA , but for a fixed point , the order is the one of the operator for which the maximum is reached in (1.29), which in practice is ℎ, MA = ℎ MA + at most points of the grid. Proof. We have to show that if , ( ℎ ) ℎ>0 , and are as in Proposition 3.3, then for any ∈ , If ∈ , then (3.9) and (3.10) follow respectively from (3.3) and (3.4), taking first the limit over ℎ and then the limit over ′ . If ∈ , then (3.9) follows from (3.3) and (3.10) is always true, since Finally, we establish stability of the proposed scheme.
The existence of a suitable function in Proposition 3.6 is a natural assumption in the setting of optimal transport. We defer discussion of this assumption to Section 5.1, and in particular to Remark 5.1.

Closed-form formula in dimension two
This section is devoted to the proof of Theorem 1.2, whose motivation is to improve the numerical efficiency of the scheme. Recall that the scheme residues are defined as the value (1.23) of a maximization problem. In Remark 4.1, we contrast the numerical cost of computing this maximal value using the explicit formula of Theorem 1.2 with a more traditional approach based on a grid search in the parameter space. In practice, and in the numerical experiments Section 6, the objective is to solve the scheme using a Newton method, which requires the following additional ingredients: (i) generating the sparse Jacobian matrix of the scheme, (ii) solving the linearized scheme, and (iii) iterating the previous two steps until the residues fall below a given threshold. Point (i) is addressed using a custom automatic differentiation library 1 , combining sparse and dense forward differentiation, and which takes advantage of the envelope theorem ( [13], Sect. 6.1) so as to efficiently differentiate the maximal value (1.23). Point (ii) relies on the standard SuperLU sparse direct solver. Point (iii) usually terminates in less than a dozen steps in practice, and the proposed scheme compares favorably to alternatives in this regard, see Section 6.4. Eventually, the evaluation of the scheme residues nevertheless accounts for a substantial part of the complexity of the proposed numerical method, and is also its most specific ingredient.
Remark 4.1 (Numerical complexity of the scheme numerical evaluation). Consider a two-dimensional Cartesian grid ℎ with ( 2 ) points. Assume that at any point ∈ ℎ , one has to perform respectively MA and BV2 operations in order to compute ℎ MA [ ] and ℎ BV2 [ ]. Then the overall numerical complexity of the scheme on the grid ℎ is ( 2 ( MA + BV2 )).
When using Theorem 1.2 in the implementation of the scheme, MA is proportional to the number of superbases in the set ℎ . As in Remark 3.4, we distinguish between the smooth case and the general case. In the smooth case, ℎ does not depend on , hence MA = (1). In the general case, if ℎ is built as in Remark B.9, with = 2/3 as suggested by Remark 3.4, then by Proposition B.10, MA = ( 2/3 log ).
For comparison, one could choose to discretize the parameter set of the maximum in the definition (1.8) of the operator ℎ MA instead of using Theorem 1.2, and in this case MA would be proportional to the number of points in this discretization. Since the set of symmetric positive semidefinite matrices of size two and of unit trace has dimension two, in order to guarantee consistency of the scheme to some order > 0, one should choose at least MA = ( 2 ). This is more costly than using Theorem 1.2, both in the smooth case (in which the desired order, according to Remark 3.4, is = 1, or even = 2 if LF = LF = 0) and in the general case (in which the desired order is = 2/3).
There is also a maximum in the definition (1.28) of ℎ BV2 which, depending on the expression of the set-valued function in (1.24), either admits a closed-form formula or needs to be discretized. If it admits a closed-form formula, then BV2 does not depend on . If it needs to be discretized, then BV2 is proportional to the number of points in the discretization and, in order to guarantee consistency of the operator ℎ BV2 with BV2 at some order > 0, one should choose BV2 = ( ), since the parameter set is one-dimensional. The numerical cost of this discretization is negligible in the general case, but not in the smooth case. In practice, in many applications, the maximum in (1.29) is only attained by ℎ BV2 [ ] at points ∈ ℎ that are close to . A perspective for future research would be to prove that one may use a variant of the scheme (1.30) which would only require computing ℎ BV2 [ ] at such points, reducing the numerical cost of handling the boundary condition (1.24).
In dimension = 2, choosing ℎ as a family of superbases of Z 2 (see Def. 1.1) is motivated by Selling's formula [43]: for any family = ( 1 , 2 , 3 ) of vectors of Z 2 , recall that we defined : and let us also define : 2 → R 3 by where we consider the indices of the elements of modulo three, and where if = ( , ) ∈ R 2 , we denote ⊥ := (− , ).
Proof of Theorem 1.2. We prove separately the two statements of the theorem.
Case of bases. Let = ( 1 , 2 ) be a basis of Z 2 , ≥ 0, and = ( 1 , 2 ) ∈ R 2 . Note that is the segment of endpoints We compute that for any ∈ [−1, 1], det (︃(︃ 1 + using the definition of for the first equality, that det( ⊗ + ⊗ ) = det( , ) 2 for any , ∈ R 2 for the second equality, and that det( 1 , 2 ) = ±1 for the third equality. After defining (0) ∈ R and (1) , (2) ∈ R 2 by (0) := 1 This is the maximum of a concave function over [−1, 1]. Writing the first order optimality condition yields that the optimal must satisfy from which we deduce the expected formula Case of superbases. We use that in the space of symmetric matrices size two equipped with the Frobenius norm, the set of symmetric positive semidefinite matrices of unit trace is a disk. More precisely, let us define the affine map D : R 2 → 2 by Note that the above definition is closely related to Pauli matrices in quantum mechanics. It is easily proved that Moreover, for any ∈ R such that | | ≤ 1, ) be a superbase of Z 2 , ≥ 0, and ∈ R 3 . The Minkowski determinant inequality states, in any dimension ∈ N, the function det(·) 1/ is concave over + . Hence the function Recall that ∈ 3 and ∈ R 3 were defined in the statement of the theorem, and note that = ⊤ . It is easily computed that for any ∈ R 2 , (D( )) = − , and thus, using also (4.4), that Therefore, * ( , ) is the argmax of a concave function over the unit disk, and writing the first-order optimality condition yields which concludes the proof.

Application to quadratic optimal transport
We specialize in this section the proposed scheme to the quadratic optimal transport problem and provide a convergence proof, taking advantage of specific tools in this setting such as Aleksandrov solutions and the mass balance equation, in addition to the generic tools introduced in Section 3.

The quadratic optimal transport problem
Let be an open bounded convex nonempty subset of R and : → R + and : → R * + be two densities satisfying the mass balance condition being continuous almost everywhere and bounded and being Lipschitz continuous. For convenience, in this paper we extend the function to the whole domain R in such a manner that −1/ : R → R * + is bounded and Lipschitz continuous.
In the quadratic optimal transport problem between and , one aims to solve the minimization problem where the unknown is a Borel map : → and the constraint # = means that for any Borel subset of , ∫︁ In the literature, it is typically assumed that: the set is convex.
For simplicity, we will sometimes assume instead that: the set is strongly convex. (5.5) It was proved in [10] (see also [44], Thm. 2.12) that, under assumption (5.4), the optimal transport problem (5.2) admits a solution which is the gradient of a convex function : → R, called the potential function of the problem. Then, if is smooth enough, it may be deduced by performing the change of variables = ( ) in the right-hand side of (5. Additionally, the constraint that ( ) = ( ) ∈ , for any ∈ , may be written as (1.24), where for any ∈ , ( ) = . (5.7) Note that in this setting, a possible choice of function in Proposition 3.6 is given by ( ) := ⟨ , 0 ⟩, for some 0 ∈ . Remark 5.1 (General optimal transport). In the general optimal transport problem, a cost function ∈ 2 (R × R ) is given, and one aims to solve If is defined by ( , ) = | − | 2 , this problem reduces to (5.2). It is also equivalent to (5.2) when ( , ) = −⟨ , ⟩, as follows directly from the equality | − | 2 = | | 2 + | | 2 − 2⟨ , ⟩.
Under suitable assumptions (see [20,37]), there exists a solution : → to (5.8) of the form ( ) = -exp ( ( )), where for any ∈ and , ∈ R , the function -exp : R → R is such that and where the function (called the potential function) is -convex, in the sense that for any 0 ∈ , there exists 0 ∈ R and 0 ∈ R such that

Weak solutions to the Monge-Ampère equation
If the open set is convex, and if : → R is a convex function and is a subset of , then we denote by ( ) the union ⋃︀ ∈ ( ), where ( ) is the subgradient of at point : A notion of weak solutions to the Monge-Ampère equation that is directly related to the optimal transport problem (5.2) is the one of Brenier solutions.  If is only defined in (respectively ), we define in the same manner after having extended with value +∞ outside (respectively ). In addition to the convex envelope : R × R of , let us define the function : R → R by ( ) := sup ∈ (⟨ , ⟩ − ( )).
One motivation for the definition of (which is similar to the definition of the function˜in Sect. It is a minimal Aleksandrov solution to (1.1) and (1.24) if moreover ( ) ⊂ .
In our setting, Brenier and Aleksandrov solutions coincide, see for instance [24] (noting that the relevant part ( [24], Sect. 1) is not specific to the dimension two): This is related to the fact that is convex and is nonnegative in , and that this does not remain true in more general settings.
We will also need to Below is the adaptation of Theorem 1.6.2 from [29] to our setting. For simplicity, it is assumed that ( ) = 1 for any ∈ R , which turns (1.1) into the basic Monge-Ampère equation det + ( 2 ( )) = ( ). Note however that we only use Theorem 5.7 as an intermediary result and that our convergence result, Theorem 5.25, is not limited to the case ( ) = 1.

Reformulation of the Monge-Ampère equation
Let us now study the reformulation of the Monge-Ampère equation (1.1) in the form (1.2), in the setting of quadratic optimal transport. We sum up the idea of the reformulation in the following proposition: Proof. We refer to Lemma 3.2.2 of [35] for the proof of the equivalence max ∈ + Tr( )=1 ( , ) = 0 ⇐⇒ = det . (5.16) Also, the first equality in (1.5) is proved in Lemma 3.2.1 of [35] (it is related to the inequality of arithmetic and geometric means applied to eigenvalues of the product 1/2 1/2 ), while the second one follows from the identity From (1.5), we deduce that Then (5.14) follows from the continuity of ( , ) with respect to ∈ + , and (5.15) follows from (5.14) and (5.16). The proof is an adaptation of the one of Proposition 1.3.4 from [29]. It uses Lemma 1.4.1 of [29], which we recall below in our setting: Proof of Proposition 5.9. We adapt the proof of Proposition 1.3.4 from [29], which is a particular case of this proposition. First let us show that is a viscosity subsolution to (1.2). Let ∈ 2 ( ), and let 0 ∈ be a local maximum of − . Since is convex, 2 ( ) must be positive semidefinite. We may assume without loss of generality that is convex, that ( 0 ) = ( 0 ), and that 0 is a strict local maximum. For any small > 0, there exists an open set such that ⊂ , ≤ + in , = + on , and lim →0 H ( , { 0 }) = 0 (see [29] for detail). By Lemma 5.10, ( ) = ( + )( ) ⊂ ( ). Thus, since is an Aleksandrov solution, Passing to the limit in , we deduce that * ( 0 ) ≤ ( ( 0 )) det 2 ( 0 ). By Proposition 5.8, it follows that ( MA ) * ( 0 , ( 0 ), 2 ( 0 )) ≤ 0, and thus that is a viscosity subsolution to (1.2). Now let us show that is a viscosity supersolution to (1.2). Let ∈ 2 ( ), and let 0 ∈ be a local minimum of − . If there exists a unit vector ∈ R such that ⟨︀ , 2 ( 0 ) ⟩︀ ≤ 0, then choosing = ⊗ in the maximum in the definition (1.8) of the operator MA yields If on the contrary 2 ( 0 ) is positive definite, then we may assume without loss of generality that is convex, that ( 0 ) = ( 0 ), and that 0 is a strict local minimum. By the same reasoning as above, we prove that * ( 0 ) ≥ ( ( 0 )) det 2 ( 0 ), and we deduce using Proposition 5.8 that ( MA ) * (︀ 0 , ( 0 ), 2 ( 0 ) )︀ ≥ 0. Therefore is a viscosity supersolution to (1.2).
In order to prove convergence of a family of monotone numerical schemes for the Monge-Ampère equation, we need to study under which conditions viscosity solutions to (3.2)  Note also those two theorems are particularly strong results, since they apply to viscosity subsolutions and supersolutions and not only to viscosity solutions as one would have expected. In the particular case of viscosity solutions, one has the following immediate corollary, which is also particularly strong since it does not assume that = 0, but instead proves this equality: Proof. Since is a viscosity solution, it is both a viscosity subsolution and supersolution. By Theorem 5.12, if ≤ 0, then = 0. This means that in any case ≥ 0. Therefore Theorem 5.11 applies, and concludes the proof.
Corollary 5.13 is the main original argument that we use in the proof of our convergence result Theorem 5.25, in combination with standard arguments [1] about the convergence of monotone schemes for degenerate elliptic equations.
Note that in the proof of Corollary 5.13, we did not use the part of Theorem 5.12 about being a minimal Aleksandrov solution to (1.1) and (1.24). We mention this fact nevertheless in the statement of Theorem 5.12 since it is a direct consequence of our proof that = 0.
Remark 5.14 (Sketch of proof of Thms. 5.11 and 5.12). The rigorous proof of Theorems 5.11 and 5.12 are delayed to the end of Section 5.3, but let us first explain the main arguments that we use in those proofs. In this remark we will only discuss the proof that (respectively ) is a minimal Aleksandrov solution to (1.1) and (1.24); in order to prove that = 0 one just needs to sufficiently refine the arguments below.
Theorem 5.11 is very close to Theorem 2.1 of [26] (although the considered reformulation of the Monge-Ampère equation is not the same) and thus we follow the same sketch of proof. If is a viscosity subsolution to (3.2) with ≥ 0, then it is a viscosity subsolution to both (1.2) and (1.26). From the fact that is a subsolution to (1.2), we deduce that it is a convex function, see Lemma 5. 16. The fact that is a subsolution to (1.24) means that the optimal transport boundary condition ( ) ⊂ is satisfied, see Lemma 5.17. We deduce from the optimal transport boundary condition the inequality On the other hand, we are able to deduce from the fact that is a viscosity subsolution to the (reformulated) Monge-Ampère equation (1.2) that for any Borel set ⊂ , one has the inequality ∫︀ ( ) ( ) d ≥ ∫︀ ( ) d , which is the inequality variant of the equality in the definition of Aleksandrov solutions. Observe that the two previous inequalities are in competition with each other. Thus we are able to show that they are actually equalities and that is therefore a minimal Aleksandrov solution to (1.1) and (1.24).
Contrary to Theorem 5.11, no counterpart to Theorem 5.12 is established in [26]. The proof of Theorem 5.11 does not translate directly to the setting of Theorem 5.12. This is because, for an arbitrary viscosity supersolution to (3.2) with ≤ 0, on the one hand one cannot expect to be a viscosity supersolution to (1.2) (contrary to the case of subsolutions), and on the other hand is not even guaranteed to be convex, so for instance the optimal transport boundary condition ( ) ⊂ does not make sense. We get around those difficulties by considering the modified convex envelope instead of the function itself. By construction, is guaranteed to be convex and to satisfy the optimal transport boundary condition ( ) ⊂ . By analyzing the meaning of the Dirichlet boundary condition − ∞ ≥ 0 in the viscosity sense, we are able to prove the converse inclusion ⊂ ( ), see Lemma 5.19. We are also able to show that, contrary to , is guaranteed to be a viscosity supersolution to (1.2), see Lemma 5.21. From this point the proof of Theorem 5.12 is similar to the one of Theorem 5.11, although the rigorous proof of the inequality ∫︀ ( ) ( ) d ≤ ∫︀ ( ) d , for Borel sets ⊂ , is a bit more technical than its counterpart in the setting of Theorem 5.11 and involves Lemma 5.22 in addition to Propositions 5.8 and 5.24.
Let us now turn to the intermediary results needed in the proof of Theorems 5.11 and 5. 12 We will need the following comparison principle for equation (1.2). The assumptions that we make on the function are more restrictive than in the rest of the paper, but this does not affect the generality of our main results since we will only need to apply this comparison principle to the case of a constant function , see Lemma 5.22. Proof. Let 0 ∈ . For any > 0, let : → R be defined by so that ≤ ≤ on . Let 1 ∈ , ∈ 2 ( ), and := + ( /2)| · − 0 | 2 . Then 1 is a local maximum of − if and only if it is a local maximum of − . For some constant > 0 and for = 1/(2 ), using that | ( 1 ) − ( 1 )| ≤ and 2 ( 1 ) = 2 ( 1 ) + , it holds for any ∈ + satisfying Tr( ) = 1 that Thus if 1 is a local maximum of − , Then by Theorem 3.3 and Section 5.C of [16], ≤ in , and we conclude by letting approach zero.
Notice that we did not need to assume (5.6); however, if (5.6) holds, it may be shown that the assumption that diam( ) ≤ is not necessary, see Theorem V.2 of [32] for the argument.
We will also need the following lemmas.
Proof. Let ∈ 2 ( ) and 0 be a local maximum of − in . Then, using that is a viscosity subsolution and choosing = ⊗ in the maximum in the definition of MA , Thus is a viscosity subsolution to By Theorem 1 of [41], it follows that is convex. The proof of Lemma 5.17 is a direct transposition to our setting to the one of Lemma 2.5 from [26], so we do not reproduce it here.
Proof. This standard result follows directly from the facts that is not twice differentiable at points of this set (since { 1 , 2 } ⊂ ( )) and that , as a convex, hence locally Lipschitz function, is differentiable almost everywhere, by Rademacher's theorem ( [22], Thm. 3.1.2).
In the lemma below, the right-hand side in (5.17) is to be understood as the integral of function which coincides almost everywhere with ( (·)) det 2 (·). Indeed, the convex function is twice differentiable almost everywhere by Aleksandrov's theorem ( [22], Thm. 6.4.1). In particular, points where is not twice differentiable do not contribute to the integral in the right-hand side, while they do contribute to the one in the left-hand side. If moreover ( ′ ) has Lebesgue measure zero for any subset ′ of of Lebesgue measure zero, then the above inequality is an equality.

=1
and, for any ∈ N * ,˜: Using Lemma 5.23, is a singleton for almost every ), with equality if ( ∖˜) has Lebesgue measure zero (note that ∖˜always has Lebesgue measure zero).
By the change of variables formula ( [22], Thm. 3.3.2) which is a corollary of the area formula of geometric measure theory, for any ∈ N * ,

It follows that
which concludes the proof.
Let us now prove the main Theorems 5.11 and 5.12.
Proof of Theorem 5.11. If : → R is a viscosity subsolution to (3.2) with ≥ 0, it is both a viscosity subsolution to (1.2) and (1.26). Thus by Lemmas 5.16 and 5.17, it is convex in and ( ) ⊂ . By Aleksandrov's theorem ( [22], Thm. 6.4.1), is twice differentiable almost everywhere. Thus it is almost everywhere a classical subsolution to (3.2). It follows that for almost every ∈ , MA (︀ , ( ), 2 ( ) )︀ ≤ 0, with a strict inequality if > 0. Then, using Proposition 5.8, for any Borel subset of , with a strict inequality if > 0 and has positive Lebesgue measure. By Lemma 5.24, we deduce that with a strict inequality if > 0 and has positive Lebesgue measure. The same is true when replacing by ∖ , and by Lemma 5.23, ( ) ∩ ( ∖ ) has Lebesgue measure zero. Thus with a strict inequality if > 0, since at least one of the sets and ∖ has positive Lebesgue measure. On the other hand, since ( ) ⊂ , one has the converse inequality Therefore the case > 0 cannot happen, and moreover from which it follows that is a minimal Aleksandrov solution to (1.1) and (1.24).

Convergence
We are now able to prove convergence of a family of numerical schemes (which includes the scheme (1.30), see Sect. 3) for the Monge-Ampère equation, in the setting of quadratic optimal transport. Proof. Let (ℎ ) ∈N be a sequence of small discretization steps ℎ > 0 converging to zero. Since (2.8) is equicontinuously stable, the sequence ( ℎ ) ∈N is bounded, and ( ℎ ) ∈N is uniformly bounded and uniformly equicontinuous. Then by the Arzelà-Ascoli theorem, up to extracting a subsequence, ℎ converges to some ∈ R and ℎ converges uniformly to some continuous function : → R, satisfying (0) = 0. By Corollary 2.13, is a viscosity solution to (3.2). By Corollary 5.13, = 0 and is the minimal Aleksandrov solution to (1.1) and (1.24), which concludes the proof. Figure 1. Image of a Cartesian grid by the approximated optimal transport maps, for the quadratic optimal transport problem with some given source and target densities.

Numerical experiments
6.1. Approximation of optimal transport maps for some quadratic optimal transport problems We apply the scheme (1.30) to the numerical resolution of some problems of the form (5.2), see Figure 1. The problems considered are inspired by the numerical experiments in [3]. The source and target domains and are chosen as the unit disk 2 (0, 1). The source density : 2 (0, 1) → R + is chosen among the ones depicted in the top row of Figure 1 and may have a non-convex or non-connected support, while the target density : 2 (0, 1) → R * + (extended to all of R 2 for numerical purposes, as explained in Sect. 5.1) is either the uniform density : ↦ → 1/ or the following combination of Gaussian densities and of a small uniform density: In all our numerical experiments, we choose the discretization step ℎ > 0 as ℎ = 2/ , for some ∈ 2N * , so that the Cartesian grid [−1, 1] 2 ∩ ℎZ 2 contains exactly ( + 1) 2 points. We define the Cartesian grid ℎ := ∩ ℎZ 2 , as well as the smaller grid ℎ := { ∈ ℎ | ∀ ∈ {1, 2}, + ℎ ∈ ℎ and − ℎ ∈ ℎ }.
In Figure 1, we choose = 128. Following Appendix B and in particular Table B.1, we choose ℎ = with = 2 + √ 5 ≈ 4.236. We solve the numerical scheme (1.30) on the grid ℎ , in the setting described by (5.6) and (5.7), and we denote by ( ℎ , ℎ ) ∈ R × R ℎ the solution to the scheme. We approximate the optimal transport map : → by ℎ ℎ :˜ℎ → R , where ℎ is the centered finite difference operator defined in (1.13). The grids displayed in Figure 1 are the image of˜ℎ (coarsened for readability) by the approximate transport map ℎ ℎ in each of the settings considered.   Figure 1, we display the number of iterations required in order to achieve convergence of the Newton method for each of the problems considered, with initialization [ ] = | | 2 and with the stopping criterion We observe that more iterations seem to be required when the size of the support of the source density is small comparatively to the source domain .

Numerical convergence analysis
In Figures 2 and 3, we display the approximation errors obtained when using the recommended scheme in order to solve some Monge-Ampère problems whose solution : → R is known exactly. According to Remark 3.4, the discretization of the Monge-Ampère operator by the operator ℎ MA defined in (1.23) is expected to achieve second-order consistency in favorable cases. However, consistency to an order higher than one cannot be expected for the whole scheme (1.30) due to the fact that the discretization (1.28) of the optimal transport boundary condition is only first-order consistent. In order to study separately the errors stemming from the discretizations of the Monge-Ampère operator and of the optimal transport boundary condition, we consider both the second and the first boundary value problems for the Monge-Ampère equation. We use the scheme (1.30) in order to approximate the solution to the second boundary value problem. The first boundary value problem involves a Dirichlet boundary condition; the variant (C.1) of our numerical scheme devoted to this setting is described in Appendix C. : → R in the case of the Dirichlet boundary condition; the source density : → R + and the target domain ⊂ R 2 in (5.7) in the case of the optimal transport boundary condition) appropriately so that the function is the solution to the problem. For the quartic and 1 problem, the target domain is = 2 (0, 1). For the singular problem, we only consider the Dirichlet boundary condition, since with the optimal transport boundary condition the target domain would be unbounded and non-convex, which does not fit in our framework.
We define the ∞ approximation error between the exact solution and the numerical solution ℎ as We also display the error between and its approximation obtained by applying to ℎ the centered finite difference operator ℎ defined in (1.13).
It may be of practical interest to approximate by ℎ ℎ since, at least in the setting of the second boundary value problem, coincides with the optimal transport map for the associated optimal transport problem, see Section 5.1. Note however that theoretical guarantees on the convergence of ℎ ℎ towards are unreachable using the techniques developed in this paper.
In the case of the singular problem, convergence of ℎ ℎ towards is not observed in the ∞ norm, which is expected since is unbounded in . For this reason, we display instead the 1 error According to Appendix B and Table B.1, we choose the set of superbases ℎ in (1.23) as ℎ = for some > 1. In Figure 2, we use the value = 1 + √ 2. According to Proposition B.8 larger values of may need to be used for small discretization steps ℎ in order to observe convergence. This is illustrated in Figure 3, where the values = 2 + √ 5 and = 3 + √ 10 are also considered.

Effect of the set of superbases on the pointwise approximation error
In order to describe more visually the effect of the choice of the parameter on the numerical solution to the Monge-Ampère problem, we display in Figure 4  . We observe as expected that the solution is better reconstructed, especially near its singularity, for larger values of , which correspond to wider finite difference stencils.

Comparison between the recommended scheme and the MA-LBR scheme
We study the behavior of the Newton method applied to the resolution of the recommended scheme (C.1) and the MA-LBR scheme (C.2), in the setting of the Dirichlet problem  on domains = 2 (0, 1)∪] − 1, 1[ 2 and = 2 (0, 1) ∖ [−1, 1] 2 . While the second of those domains does not fit in standard theoretical frameworks for the Monge-Ampère equation due to being non-convex, this choice has to be considered as a stress test for the considered numerical methods.
Let us denote by 0 : ℎ → R our initial guess in the Newton method. The iterates of the Newton method are defined as := −1 + 2 − , where is the Newton descent direction and ∈ N is a damping parameter. In the case of the recommended scheme (C.1), no damping is required, so we always choose = 0. In the case of the MA-LBR scheme (C.2), one has to use damping so that the constraint (C.3) remains satisfied along the iterates, as discussed in Section 1.3. More precisely, let us introduce the following quantitative variant of the constraint (C. (in the setting of our experiments, one has˜ℎ [ ]/2 = 1/2). Following [5,39], and in the spirit of [33], we assume that (6.2) is satisfied for 0 and at each iteration we let be the smallest natural number such that (6.2) holds.
We use the grid size = 120 and the set of superbases ℎ = , = 2 + √ 5. We initialize the Newton method with 0 [ ] := | | 2 − 2 (since the simpler initialization 0 [ ] := | | 2 does not satisfy (6.2) close to , in view of the boundary condition = 0 on ) and we use the stopping criterion where either ℎ = ℎ MABV1 or ℎ = ℎ MA-LBR-BV1 as appropriate. On the domain = 2 (0, 1)∪] − 1, 1[ 2 , the Newton method for the recommended scheme (C.1) converges in 9 iterations without damping. The Newton method for the MA-LBR scheme converges in 47 iterations and the largest observed value for is = 4, corresponding to a damping step 2 − = 0.0625. On the domain = 2 (0, 1) ∖ [−1, 1] 2 , the Newton method for the recommended scheme converges in 7 iterations without damping. The Newton method for the MA-LBR scheme converges in 52 iterations and the largest observed value for is = 5, corresponding to a damping step 2 − = 0.03125.

Application of the scheme to the far field refractor problem in nonimaging optics
We apply the scheme (1.30) to the far field refractor problem [30] in nonimaging optics. Various numerical methods for solving this problem, and some of its variants such as the reflector problem, have been previously studied in the literature [6,11,18,19,28]. We illustrate the refractor problem in Figure 5. Light rays emanate from a point source of light located at the origin, in directions belonging to some subset of the upper hemisphere {︀ ∈ 2 | 3 > 0 }︀ of the unit sphere 2 . In our experiments we assume that the intensity of those light rays is constant with respect to their direction. The rays propagate in the ambient medium, whose index of refraction we denote by 1 > 0, until they hit a lens, which is located at distance ℓ 0 from the origin and whose index of refraction we denote by 2 > 0. The rays are refracted by the lens, then continue to propagate in the ambient medium until they hit a screen, represented by the plane R 2 × {ℓ}, where ℓ > ℓ 0 denotes the distance from the screen to the origin. The aim is to find a suitable shape for the lens that refracts the light rays so that a given target image is reconstructed on the screen.
The far field refractor problem is studied in the limit ℓ → ∞. In this limit, it has been shown [30] to be equivalent to an optimal transport problem with a specific, non-quadratic cost, and thus to reduce to the second boundary value problem for the Monge-Ampère equation (1.1) with coefficients of the form (5.10) to (5.12) (as opposed to (5.6) and (5.7)). In the above-mentioned optimal transport problem, the source and target densities, which may be exchanged up to an appropriate transformation of the transport map, are the density describing the image that has to be reconstructed in the refractor problem and the one describing the intensity of light rays depending on their initial direction of emission.  The difference between (6.3) and (1.1) is that in (6.3) the coefficients and depend on the values of the function and not only on its derivatives. The study of the extension of the scheme (1.30) to generated Jacobian equations is outside the scope of this paper and is an opportunity for future research.
We approximate the solution to the far field refractor problem by solving the scheme (1.30) on the unit disk = 2 (0, 1), choosing as the source density the one describing the target image in the problem. We use the grid size = 120 and the set of superbases ℎ = , = 2 + √ 5 (see Tab. B.1). We choose the indices of refraction 1 = 1 and 2 = 1.5, corresponding to a glass-air interface. Eleven iterations of the Newton method are needed in order to solve the scheme. Up to an appropriate postprocessing of the solution to the scheme (1.30), we obtain an approximation of the height map : 2 (0, ) → R + describing the upper interface {( , ( )) | ∈ 2 (0, )} of the lens, where > 0 denotes the radius of the lens.
We display our numerical results in Figure 6. On the representation of the approximation of the pointwise curvature of , we observe that the parts of the refractor corresponding to dark areas of the image have a small area, compared to the ones corresponding to bright areas. This is consistent with the fact that the total intensity of the light traversing them should be low, in order for the image to be properly reconstructed. In order to validate our results, we inject the shape that we obtain for the lens in a simulation of the scene that we perform using the appleseed R ○2 rendering engine.

Conclusion and perspectives
We were able to adapt Perron's method in order to prove the existence of solutions to a class of monotone numerical schemes whose sets of solutions are stable by addition of a constant. We designed a finite difference scheme for the Monge-Ampère equation that belongs to this class, and proved convergence of the scheme in the setting of quadratic optimal transport. We showed that in dimension two, the discretization of the Monge-Ampère operator admits a closed-form formulation, and thus yields a particularly efficient numerical method, when carefully choosing its parameters using Selling's formula. We validated the method by numerical experiments in the context of the far field refractor problem in nonimaging optics.
A natural perspective is the adaptation of the proof of convergence of the scheme to the setting of more general optimal transport problems. The extension of the scheme to generated Jacobian equations such as (6.3) could also be studied. This would require adapting both the proof of convergence and the one of existence of solutions to the scheme, since the invariance in the set of solutions would not be the same in this case.
Another perspective is the extension of this work to higher dimensions. While our existence and convergence results are valid in any dimension, the closed-form formula that we obtain for the maximum in the discretized Monge-Ampère operator is specific to the dimension two. In higher dimensions, this maximum could be approximated either by sampling the parameter set or by resorting to a numerical optimization procedure, since (1.23) is an instance of a semidefinite program. We expect the second approach to be more efficient, but developing such an optimization procedure is still an open research direction. The design of this procedure could be made easier by an appropriate choice of the set ℎ in (1.23). In dimension three, one could choose it as a set of superbases of Z 3 , benefiting from the fact that Selling's formula (described in Prop. 4.2 in dimension two) also holds in dimension three; in this case, which superbases exactly the set ℎ should contain is an open question, since the selection criterion based on the Stern-Brocot tree and presented in Appendix B is two-dimensional only. In dimensions four and higher, one could resort to Voronoi's first reduction of quadratic forms [15], which generalizes Selling's formula to those dimensions.
Contrary to the operator ℎ MA , the operator ℎ MA-LBR is only intended to be applied to functions : ℎ → R satisfying the constraint Λ ℎ [ ] > 0, ∀ ∈ ℎ , (A. 4) which is a discrete counterpart to the strict variant of the admissibility constraint (1.4). If this condition fails, then the Jacobian matrix of the scheme is not invertible [39], which breaks the Newton method relied upon. Recall that, for any superbase ∈ ℎ and any ∈ R 3 , one has ( ) := ∑︀ 3

=1
⊗ . The following proposition shows that the MA-LBR scheme may be seen as a discretization of the reformulation (1.6) of the Monge-Ampère equation: By Lemma A.1, one has ( ) = 0 ( ) and > 0 elementwise. Using that is a superbase of Z 2 , and that therefore det( , ) = ±1 for any 1 ≤ < ≤ 3, one can compute that for any ∈ R 3 + , .

Appendix B. Choosing the set of superbases in dimension two
In this appendix, we explain how one may choose, in dimension = 2 and for any ℎ > 0, a finite set ℎ of superbases of Z 2 satisfying (1.20)- (1.22). The motivation is to use this set ℎ in (1.23). The construction of ℎ is based on the Stern-Brocot tree of bases of Z 2 (see [7] for a similar approach in the setting of Hamilton-Jacobi-Bellman equations): Definition B.1. A pair ( , ) of vectors of Z 2 is a direct basis of Z 2 if det( , ) = 1.

Definition B.2. The Stern-Brocot tree
is the collection of direct bases of Z 2 defined inductively as follows: (i) the canonical basis belongs to , and (ii) for any ( , ) ∈ , one has ( , + ) ∈ and ( + , ) ∈ . Remark B.3. In classical descriptions of the Stern-Brocot tree, the vector = ( , ) is often identified with the ratio / , which is a nonnegative rational, or with +∞, and likewise for = ( , ) (note that and are nonnegative and coprime by construction).
For any > 1, we define := ⋃︁ )︀ Then the scheme that we use for the Dirichlet problem may be written as The complete study of the theoretical properties of this scheme, such as monotonicity and convergence, is outside the scope of this paper. In Section 6.4, we also apply the MA-LBR scheme, see Appendix A, to the Dirichlet problem. The scheme The MA-LBR scheme in this setting may be written as The operator ℎ MA-LBR-BV1 is intended to be applied to functions : ℎ → R satisfying the admissibility con-straintΛ which is the natural counterpart to (A.4) in the Dirichlet setting.