TY - JOUR AU1 - Abbasi,, Bilal AU2 - Oberman, Adam, M AB - Abstract Quasiconvex (QC) functions are functions whose level sets are convex. In a series of articles, Barron, Goebel and Jensen studied partial differential equations (PDEs) for QC functions. They introduced a regularization of the PDE that is more stable. In this article, we introduce a different regularization that is more amenable to numerical approximation. We prove existence and uniqueness of viscosity solutions of the PDE and demonstrate that level sets of solutions of the new PDE are uniformly convex. We build convergent finite difference approximations, comparing the QC envelope with the regularization. 1. Introduction In a series of articles from about four years ago, Barron, Goebel and Jensen introduced and studied partial differential equations (PDEs) for quasiconvexity (Barron et al., 2012a,b, 2013; Barron & Jensen, 2014). In this context, quasiconvexity means that the sublevel sets of a function are convex. The study of convexity of level sets for obstacle problems has a long history, which includes Caffarelli & Spruck (1982) and Kawohl (1985), as well as the more recent work by Colesanti & Salani (2003) and the references therein. Quasiconvex (QC) functions appear naturally in optimization, since they generalize convex functions. The property also appears in economics (Avriel et al., 2010). Earlier work by one of the authors studied a PDE for the convex envelope (Oberman, 2007), which led to a numerical method for convex envelopes (Oberman, 2008a,b). The notion of (scalar) quasiconvexity that we discuss here should not be confused with the well-known notion of quasiconvexity in vector variational problems. Our motivation for studying this PDE comes from two directions: geometric and algorithmic; it is natural to represent a set using a level set function (Sethian, 1999). The PDE allows us to find the convex hull of a set represented this way directly (in fact, every level set), rather than extracting the level set and then finding the convex hull, followed by generating a new level set function. In addition, we can use the parabolic version of the equation to progressively convexify a level set, or to make already convex level sets uniformly convex. Quasiconvexity is challenging because, unlike convexity, it is a nonlocal property (at least for functions that have flat parts). This means that, even using viscosity solutions, there is no local characterization for quasiconvexity. On the other hand, using the more stable notion of robust quasiconvexity, Barron, Goebel and Jensen showed that these functions are characterized in the viscosity sense (Crandall et al., 1992) by a PDE (Barron et al., 2013; Barron & Jensen, 2014). One motivation for this work was to build numerical solvers for the QC envelope PDE. However, we had difficulties with both the QC and the robust-QC operators: the former lacks uniqueness and the latter uses an operator defined over small slices of angles (see Fig. 3), which leads to poor accuracy when using wide-stencil finite difference schemes. An alternative, presented in Barron et al. (2013), was to use first-order nonlocal PDE solvers. In a companion article (Abbasi & Oberman, 2016), we built a nonlocal solver for the QC and robust-QC envelopes. The nonlocal PDEs can be solved explicitly and the problem implemented efficiently. By iteratively solving for the envelopes on lines in multiple directions, we approximated the solution of the problem in higher dimensions. However, we are still interested in the PDE approach, which has advantages arising from its nonlocal nature. In this article, we build on the results of Barron et al. (2012a,b) and Barron & Jensen (2014) to obtain a PDE for |$\epsilon$|-uniformly QC functions, which we define subsequently. Note that |$\epsilon$|-uniform quasiconvexity implies robust quasiconvexity. The choice of terminology lends itself to the fact that |$\epsilon$|-uniformly QC functions have uniformly convex level sets. Following the argument in Barron & Jensen (2014), we establish uniqueness of viscosity solutions for the PDE. Moreover, this operator is defined for all direction vectors, which makes it amenable to discretization using wide-stencil schemes. We consider the obstacle problem for the |$\epsilon$|-uniformly quasiconvex envelope (QCE). As is the case for robust-QC, we recover the QCE as the regularization parameter |$\epsilon \to 0$|⁠. While robust-QC functions can have corners in one dimension, strictly QC functions are smoother; see Fig. 3 below, and the explicit formula for the solution in one dimension in Section 2.3. We also build and implement convergent elliptic finite difference schemes (Barles & Souganidis, 1991; Oberman, 2006) for the envelopes. These are wide-stencil finite difference schemes, which can be developed using ideas similar to Oberman (2004, 2008a,b). Solutions to these PDEs can be found using an iterative method which is equivalent to the explicit Euler discretization of the parabolic equation (Oberman, 2006). However, the method has a nonlinear Courant-Friedrichs-Lewy (CFL) condition, which restricts the step size. We find that alternating the line solver with several iterations of the parabolic PDE solver significantly improves the speed of the solution. Numerical solutions show that we obtain very similar results to the line solver for the QCE with small |$\epsilon = h/2$|⁠, where |$h$| is the grid resolution. We also compare large |$\epsilon$| solutions with comparable robust QCE and find that solutions are smoother. Formally, we show that solutions are |$\epsilon$| uniformly convex. The QC operator in two dimensions recovers the level set curvature operator. We show that our discretization of our operator agrees with the Kohn–Serfaty (Kohn & Serfaty, 2007) first-order representation of the mean curvature operator in two dimensions. See Section 4.4. 1.1 Convexity of level sets of a function We give a brief informal derivation of the operator. Given a smooth function |$u: {\mathbb{R}^n} \to \mathbb{R}$|⁠, the direction of the gradient at |$x_0$|⁠, |$p = {\nabla} u(x_0)$| is the normal to the sublevel set |$\{ x \in {\mathbb{R}^n} \mid u(x) \leq u(x_0) \}$|⁠. The curvatures of the level set at |$x_0$| are proportional to the eigenvalues of the Hessian of |$u$|⁠, |$M = D^2u(x_0)$|⁠, projected onto the tangent hyperplane at |$x_0$|⁠, |$P = P_{x_0} = \{ v \in {\mathbb{R}^n} \mid v\cdot p = 0\}$|⁠. These curvatures are all positive if |$v^\intercal M v \geq 0$| for all |$v \in P$|⁠. Thus, formally, the condition for local convexity of the level set is non-negativity of the operator \begin{equation} L_0(p,M) \equiv \min_{{\vert{v}\vert} = 1} \left \{ v^\intercal M v \mid v\cdot p = 0 \right. \}. \end{equation} (1.1) This is the operator considered in Barron et al. (2013) to study QC functions. However, for technical reasons discussed below, they chose to relax the constraint |$v \cdot p =0$| to an inequality constraint |${\vert{v \cdot p}\vert} \leq \epsilon$|⁠, resulting in the operator \[ L_\epsilon(p,M) \equiv \min_{{\vert{v}\vert} = 1} \left \{ v^\intercal M v \mid {\vert{v \cdot p}\vert}\le\epsilon \right \}\!. \] Our operator is obtained by instead replacing the hard constraint |${\vert{v\cdot p}\vert} \leq \epsilon$| with a penalty in the objective function. So we define \[ F^\epsilon(p,M) \equiv \min_{{\vert{v}\vert} = 1}\left \{ v^\intercal M v + \tfrac{1}{{\epsilon}} {\vert{v^\intercal p}\vert} \right \}\!. \] This choice of penalty gives the operator for uniformly convex level sets, as we show subsequently. Refer to Figure 2 for a visualization of the feasible sets. Figure 3 illustrates the important differences in the notions of quasiconvexity discussed so far. Moreover, on a bounded domain, the |$\epsilon$|-robustly QCE can have an interval of global minima, whereas our solution has a unique (global) minimum. 1.2 Basic definitions We recall some basic definitions and establish our notation. For a reference, see Boyd & Vandenberghe (2004). The set |$S$| is convex if whenever |$x$| and |$y$| are in |$S$| then so is the line segment |$[x,y]\equiv\{tx + (1-t)y$|⁠: |$t \in [0,1]\}$|⁠. We say that the function |$u: {\mathbb{R}^n} \to \mathbb{R}$| is convex if \begin{equation} u(tx + (1-t)y) \le t u(x) + (1-t) u(y) \quad \text{ for every } x,y \in {\mathbb{R}^n}, \quad 0 \le t \le 1. \end{equation} (1.2) The convex envelope of a function |$g$|⁠, hereby denoted |$g^{\rm CE}$|⁠, is the largest convex function majorized by |$g$|⁠: \begin{equation} g^{\rm CE}(x) \equiv \sup \{ v(x)\mid v\,\text{is convex and}\,v\le g \}. \end{equation} (1.3) In terms of sets, if |$S$| is given by the sublevel set of a function |$u$|⁠, \begin{equation} S = S_\alpha(u)\equiv\{x \in {\mathbb{R}^n} \mid u(x)\leq\alpha\}, \end{equation} (1.4) then convexity of |$S$| is equivalent to the following condition: \[ u(x), u(y) \le \alpha \implies u(tx + (1-t)y) \le \alpha \quad \text{ for every } 0 \le t \le 1. \] This condition, when applied to every level set |$\alpha\in\mathbb{R}$|⁠, characterizes quasiconvexity of the function |$u$|⁠. That is, |$u$| is QC if every sublevel set |$S_\alpha(u)$| is convex. Equivalently, |$u$| is QC if \begin{equation} u(tx + (1-t)y) \le \max \{ u(x), u(y) \} \quad \text{ for every } x,y \in {\mathbb{R}^n},\quad 0 \le t \le 1. \end{equation} (1.5) Given |$g:{\mathbb{R}^n} \to \mathbb{R}$|⁠, the QCE of |$g$| is given similarly by \begin{equation} g^{\rm QCE}(x) \equiv \sup \{ v(x)\mid v\,\text{is QC and}\,v\le g\}. \end{equation} (1.6) Remark 1.1. Since the maximum of any two QC (convex) functions is QC (convex) as well, the suprema in (1.3) and (1.6) are well defined. Figure 1 provides a visual comparison between the convex envelope and the QCE in one dimension. Fig. 1. View largeDownload slide Comparison of a function (solid) and its envelopes (dashed). (a) Quasiconvex envelope. (b) Convex envelope. Fig. 1. View largeDownload slide Comparison of a function (solid) and its envelopes (dashed). (a) Quasiconvex envelope. (b) Convex envelope. Fig. 2. View largeDownload slide Constraint sets (dotted and solid) of different PDEs in two dimensions. Fig. 2. View largeDownload slide Constraint sets (dotted and solid) of different PDEs in two dimensions. Fig. 3. View largeDownload slide Comparison between solutions of different operators. Graphs arranged top to bottom: obstacle, QCE, |$\epsilon$|-robustly QCE, and |$\epsilon$|-uniformly QCE. Fig. 3. View largeDownload slide Comparison between solutions of different operators. Graphs arranged top to bottom: obstacle, QCE, |$\epsilon$|-robustly QCE, and |$\epsilon$|-uniformly QCE. 1.3 Viscosity solutions Suppose that |$\varOmega\subset{\mathbb{R}^n}$| is a domain. Let |$S^n$| be the set of real symmetric |$n\times n$| matrices and take |$N\le M$| to denote the usual partial ordering on |$S^n$|⁠, namely that |$N-M$| is negative semidefinite. Definition 1.2. The operator |$F(x,r,p,M):\varOmega\times\mathbb{R}\times{\mathbb{R}^n}\times S^{n} \to\mathbb{R}$| is degenerate elliptic if \[ F(x,r,p,M)\le F(x,s,p,N)\quad\text{whenever}\ r\le s\ \text{and}\ N\le M. \] Remark 1.3. For brevity, we use the notation |$F[u](x)\equiv F(x,u(x),\nabla u(x), D^2u(x))$|⁠. Definition 1.4 (Viscosity solutions). Suppose |$F$| is a degenerate elliptic operator, as defined above. We say that the upper semicontinuous (lower semicontinuous) function |$u:\varOmega\to\mathbb{R}$| is a viscosity subsolution (supersolution) of |$F[u]=0$| in |$\varOmega$| if for every |$\phi\in C^2(\varOmega)$|⁠, whenever |$u-\phi$| has a strict local maximum (minimum) at |$x\in\varOmega$|⁠, \[ F(x,u(x),\nabla\phi(x),D^2\phi(x))\le0\ (\ge0). \] Moreover, we say that |$u$| is a viscosity solution of |$F[u]=0$| if |$u$| is both a viscosity subsolution and a supersolution. 1.4 The QCE Use the notation |$\lambda_{\rm QC}[u](x) \equiv L_0({\nabla} u(x),D^2u(x))$|⁠. The obstacle problem for the QCE is given by \begin{equation} \begin{cases} \max\{u-g,-\lambda_{\rm QC}[u](x)\}=0, & x\in\varOmega,\\ u =g, & x\in\partial\varOmega. \end{cases} \end{equation} (Ob) This PDE can have multiple viscosity solutions for a given |$g$| (Barron et al., 2013). However, there is a unique QC solution which is the QCE of |$g$|⁠. Failure of uniqueness in general can be seen in the following counterexample, which is similar to Barron et al. (2013, Example 3.1). Example 1.5. Consider |$\varOmega=\{x^2+y^2<1\}, u(x,y) = -y^4, v(x,y) = -(1-x^2)^2$|⁠, so that |$u=v$| on |$\partial\varOmega$|⁠. Now consider |$g\equiv u$| on |$\varOmega$|⁠. We can calculate that |$\lambda_{\rm QC}[u]=\lambda_{\rm QC}[v]=0$|⁠. In summary, |$u$| and |$v$| are both solutions of (Ob), but |$u\neq v$| in |$\varOmega$| (in fact |$v0$| and |$\epsilon > 0$|⁠. Comparison also holds for the corresponding obstacle problem. The proof we present is based on the uniqueness proof of |$\lambda_{\rm QC}[u]=g$| for |$g>0$|⁠, presented in Barron & Jensen (2014). The result is simpler, because our operator is continuous as a function |$(p,M)\mapsto F^\epsilon(p,M)$| for |$\epsilon > 0$|⁠. Write \begin{equation} F^\epsilon[u](x)=F^\epsilon(\nabla u(x),D^2u(x)) \equiv \min_{{\vert{v}\vert}=1}\left\{\tfrac{1}{{\epsilon}}{\vert{\nabla u(x)\cdot v}\vert} + v^\intercal D^2u(x)v\right\}\!. \end{equation} (2.1) Note that |$-F^\epsilon(p,M)$| is elliptic by Definition 1.2. We consider the following PDEs: (⁠|$\epsilon$|-QC) \begin{gather} h(x)-F^\epsilon[u](x) = 0, \quad x\in\varOmega,\\ \end{gather} (⁠|$\epsilon$|-QCE) \begin{gather} \max\{u(x)-g(x),h(x)-F^\epsilon[u](x)\} = 0, \quad x\in\varOmega, \end{gather} where |$\varOmega\subset{\mathbb{R}^n}$| is an open, bounded and convex domain, and |$g:\overline\varOmega\to\mathbb{R}$| is continuous. In the latter equation, |$g$| is the obstacle. We also impose the following condition on |$h$|⁠: \begin{equation} h:\overline\varOmega\to\mathbb{R}\,\text{is continuous and positive.} \end{equation} (2.2) Enforce the following Dirichlet boundary data: \begin{equation} u(x) = g(x),\quad x\in\partial\varOmega. \end{equation} (Dir) Remark 2.1 (Continuity up to the boundary). In general, viscosity solutions of (ϵ-QCE) need not be continuous up to the boundary. To apply Barles & Souganidis (1991) for convergence of numerical schemes, we need a strong comparison principle, which requires that solutions be continuous up to the boundary. The following assumption ensures continuity up to the boundary There is a convex domain |$\varOmega_L\supset\varOmega$| and a continuous, QC function |$g_0: \varOmega_L \to \mathbb{R}$|⁠, with \[ g_0(x)\le g(x) \quad \text{ for } x \in \varOmega, \qquad g_0(x)=g(x) \quad \text{ for } x\in \varOmega_L\setminus\varOmega. \] Continuity follows because we have |$g_0 \leq u \leq g,$| which gives |$u = g$| on |$\partial \varOmega$|⁠. An alternative to this condition is to prove convergence in a neighbourhood of the boundary. The convergence proof in this setting can be found in Froese (2016). Next, we state a technical, but standard, viscosity solutions result, which gives the comparison principle in the case where we have strict subsolution and supersolution. Proposition 2.2 (Comparison principle for strict subsolutions; Crandall et al., 1992). Consider the Dirichlet problem for the degenerate elliptic operator |$F(p,M)$| on the bounded domain |$\varOmega$|⁠. Let |$u \in {\rm USC}(\bar \varOmega)$| be a viscosity subsolution and let |$v \in {\rm LSC}(\bar \varOmega)$| be a viscosity supersolution. Suppose further that for |$\sigma > 0$|⁠, \begin{align*} F[u] + \sigma &\le 0 \text{ in } \varOmega, \\ F[v] &\ge 0 \text{ in } \varOmega \end{align*} holds in the viscosity sense. Then, the comparison principle holds: \[ \text{ if }\,u\leq v\,\text{on}\,\,\partial\varOmega\,\,\text{then}\,u\le v\,\,\text{in}\,\varOmega. \] Remark 2.3. In Crandall et al. (1992, Section 5.C), it is explained how the main comparison theorem (Crandall et al., 1992, Theorem 3.3) can be applied when it is possible to perturb a subsolution to a strict subsolution. This version of the theorem is what we state in Proposition 2.2. This result was used in Bardi & Mannucci (2006, Theorem 3.1) and Bardi & Mannucci (2013) to prove a comparison principle. Using the the same perturbation technique from Barron & Jensen (2014), and consequently applying Proposition 2.2, we obtain the following comparison result. Proposition 2.4 (Comparison principle). Consider the Dirichlet problem given by (ϵ-QC) (Dir) and assume (2.2) holds. Let |$u\in {\rm USC}(\overline\varOmega)$| be a viscosity subsolution and |$v\in {\rm LSC}(\overline\varOmega)$| be a viscosity supersolution of (ϵ-QC). Then, the comparison principle holds: \[ \text{ if}\,\,u\leq v \,\, \text{on}\,\partial\varOmega \,\,\text{then}\,u\le v \,\,\text{in}\,\, \varOmega. \] Proof. We will show that we can perturb |$u$| to a function |$u^\sigma$| satisfying |$u^\sigma\le v$| on |$\partial\varOmega$| and \[ h-F^\epsilon[u^\sigma]<0 \quad\text{in}\, \varOmega \] in the viscosity sense. Applying Proposition 2.2, we will have |$u^\sigma\le v$| in |$\varOmega$|⁠. Taking |$\sigma\to0$| yields the desired result. Fix |$\sigma>0$| and define the following perturbation of |$u$|⁠, and notice that for |$x\in\partial\varOmega$| we have the following relation: \[ u^\sigma(x) \equiv u(x) - \sigma \left(\max_{y\in\partial\varOmega}u(y)-u(x)\right)\le u(x) - \sigma (u(x)-u(x)) = u(x)\le v(x); \] that is, |$u^\sigma\le v$| on |$\partial\varOmega$|⁠. Next, because |$u$| is a subsolution, we have \[ h-\frac{1}{1+\sigma}F^\epsilon[u^\sigma] \le 0, \] leading to \[ h - F^\epsilon[u^\sigma]\le -\sigma h \le -\sigma h_{\min} <0\quad\text{where}\,h_{\min}\equiv \min_{x\in\overline\varOmega}h(x). \] □ We use this result to prove a weak comparison principle for the obstacle problem given by (ϵ-QCE) and (Dir). Corollary 2.5 (Comparison principle for the obstacle problem). Consider the Dirichlet problem given by (ϵ-QCE), (Dir) and assume (2.2) holds. Let |$u\in {\rm USC}(\overline\varOmega)$| be a viscosity subsolution and |$v\in {\rm LSC}(\overline\varOmega)$| be a viscosity supersolution of (ϵ-QCE). Then, the comparison principle holds: \[ \text{ if}\,\, u\leq v\,\,\text{on}\,\,\partial\varOmega\,\,\text{then}\,\,u\le v\,\,\text{in}\,\,\varOmega. \] Proof. We begin by considering the domain |$\varOmega_g\equiv\{x\in\varOmega:v(x)1,\\ -\epsilon^2x + C^-(1-\exp(x/\epsilon)), & S<-1,\\ \epsilon^2{\vert{x-x^*}\vert} + \epsilon^3(\exp(-{\vert{x-x^*}\vert}/\epsilon)-1) + u_0, & {\vert{S}\vert}<1, \end{cases} \] where \[ C^\pm \equiv \frac{H\mp\epsilon^2 W}{\pm\epsilon(1-\exp(\mp W/\epsilon))} \] and |$u_0\in\mathbb{R}$| and |$x^*\in I = (0,W)$| are constants. Proof. If there exists |$x_0$| in the interval |$I$| such that |$u_x(x_0)=0$| then (2.3) implies |$u_{xx}(x_0)>0$|⁠, that is, if there exists an interior critical point, then |$u$| must be strictly convex in |$I$|⁠. This allows us to break down the analysis of (2.3), restricted to each interval, into several cases: (i) |$u_x<0$| in |$I$|⁠, (ii) |$u_x>0$| in |$I$| and (iii) |$u_x(x_0)=0$| for some |$x_0 \in I$|⁠. In each case, we can solve a linear second-order ordinary differential equation. The case where |$S = \pm 1$| is degenerate: the solution is linear. The case where |${\vert{S}\vert} > 1$| is not difficult. The final case, where |${\vert{S}\vert} < 1$| corresponds to (iii). In this case, |$u_x<0$| for |$x0$| for |$x>x_0$|⁠. Then, \begin{equation} u(x) = {\vert{x-x_0}\vert} + (\exp(-{\vert{x-x_0}\vert})-1) + u_0. \end{equation} (2.4) For some |$u_0\in\mathbb{R}$|⁠. Taylor expansion shows that |$u(x) = x^2/2 + \mathcal{O}(x^3)$| for |$x$| near |$x_0$|⁠. Finding |$x_0$| (or |$u_0$|⁠) analytically is infeasible. But we can argue that the solution is correct by a continuity argument. Given the function |$u$| defined by (2.4), define the continuous function |$S_1(y)\equiv\frac{u(y+W_1)-u(y)}{W_1}$|⁠. For small |$W_1$| and |$y\ll x_0$|⁠, we have that |$S_1(y)<-1$| (by the earlier discussion). Similarly for |$y\gg x_0$|⁠, we have |$S_1(y)>1$|⁠. Therefore, by the intermediate value theorem there exists |$y^*$| such that |$S_1(y^*)=S$|⁠. This leads to the correct choice of constants in (2.4) that achieve the boundary values. □ 3. Geometric properties of the operator To better understand the geometry of the solutions of |$\epsilon-F^\epsilon[u]=0$|⁠, we present, in this section, some results regarding the convexity of the level sets of subsolutions of related PDEs. We give a proof of uniform convexity of the level sets for |$C^2$| subsolutions of the equation. The general case is beyond the scope of this article. The proof of robust quasiconvexity for viscosity solutions is technical, involving the construction of special test functions (Barron et al., 2013). Analogous results can also be found in Lions (1998) for convexity, and Carlier & Galichon (2012) for semiconcavity of solutions of parabolic equations. Next, we quantify the degree to which a function can lack quasiconvexity, given that it is uniformly QC in only certain directions. This is relevant for numerical solutions that enforce the operator in only a finite number of grid directions. A related result on semiconcavity of solutions of a parabolic equation can be found in Carlier & Galichon (2012, Section 4). We first show that we can recover a convex function by enforcing strict convexity along a prescribed set of directions. This result is of independent interest for the related numerical method (Oberman, 2008a). Next, we prove an analogous result for QC functions. In what follows, it suffices to consider subsolutions of |$\epsilon - \lambda_{\rm QC}[u]=0$|⁠. This is because of the ordering of the operators: \begin{equation} -\lambda_{\rm QC}[u]\ge-\lambda_{\rm QC}^{\epsilon^2}[u]\ge \epsilon - F^\epsilon[u]. \end{equation} (3.1) The first inequality follows naturally from the restriction of the constraint set. The second inequality is evident from the proof of Proposition 1.9. 3.1 Uniform convexity of the level sets of solutions Definition 3.1 (Uniform convexity of sets). Let |$S\subset{\mathbb{R}^n}$| be a domain and suppose |$x,y\in\partial S$|⁠. Define the midpoint |$\bar x\equiv(x+y)/2$|⁠. Then we say |$S$| is uniformly convex if |${\text{dist}}(\bar x,\partial S) = \mathcal{O}(h^2)$|⁠. Proposition 3.2 (Uniformly convex level sets of subsolutions). Suppose |$u\in C^2(\overline\varOmega)$| is a subsolution of |$\epsilon-\lambda_{\rm QC}[u]=0$|⁠. Then, |$u$| has uniformly convex level sets. Proof. Suppose |$x,y\in\varOmega$| such that |$u(x)=u(y)=\alpha$|⁠. Let |$\bar x$| be the midpoint of the line segment joining |$x$| and |$y$|⁠, |$h = {\vert{\bar x - y}\vert}$| and |$d = (\bar x-y)/{\vert{\bar x -y}\vert}$|⁠. We see that |${\nabla} u(\bar x)^\intercal d= 0$|⁠, so that we have \begin{align*} \epsilon &\le \min_{{\vert{v}\vert}=1,{\nabla} u(\bar x)^\intercal v=0} v^\intercal D^2u(\bar x) v\\ &\le d^\intercal D^2u(\bar x) d\\ &= \frac{u(\bar x + hd) + u(\bar x -hd) - 2u(\bar x)}{h^2} +\mathcal{O}(h^2)\\ &= \frac{2\alpha - 2u(\bar x)}{h^2} + \mathcal{O}(h^2). \end{align*} Rearranging for |$u(\bar x)$| yields the following inequality: \[ u(\bar x) \le \alpha -\frac{\epsilon h^2}{2} + \mathcal{O}(h^4). \] □ 3.2 Directional quasiconvexity Definition 3.3. The function |$u:\mathbb{R}^n\rightarrow\mathbb{R}$| is QC along a direction |$v$| if for any |$x\in\mathbb{R}^n$|⁠, the function |$\tilde u(t)\equiv u(x+tv)$| is QC for |$t\in\mathbb{R}$|⁠. We also appeal to the following proposition, which is elementary from the definition of quasiconvexity. Proposition 3.4 The function |$u$| is QC if and only if |$u$| is QC in every direction |$v$|⁠. Proposition 3.4 provides a convenient characterization of QC functions. In practice, however, we are confined to a grid and thus cannot enforce quasiconvexity along every direction. Therefore, we relax the notion of directional quasiconvexity to only a finite set of directions. Doing so results in the notion of approximate quasiconvexity. We can quantify the degree to which a function might lack quasiconvexity, expressed in terms of the directional resolution which we define below. Definition 3.5 (Directional resolution). Let |$\mathcal{D}=\{d_1,\dots,d_n\}$| be a set of unit vectors. Then we define the directional resolution of |$\mathcal{D}$| as \begin{equation} {\rm d}\theta \equiv \max_{{\vert{w}\vert}=1}\min_{d\in \mathcal{D} }\cos^{-1}(w^\intercal d). \end{equation} (3.2) In two dimensions, |${\rm d}\theta$| is the largest angle an arbitrary unit vector can make with any vector in |$\mathcal{D}$|⁠. In two dimensions, |${\rm d}\theta$| is simply half the maximum angle between any two direction vectors. Recalling that a necessary condition for a function |$u$| to be QC is |$-\lambda_{\rm QC}[u]\le0$| in the viscosity sense (Barron et al., 2013), we have the following result. Proposition 3.6 (Approximate quasiconvexity). Suppose |$u\in C(\overline\varOmega)$| and |$\mathcal{D}$| is a set of directions with directional resolution |${\rm d}\theta\leq\frac{\pi}{4}$|⁠. Also, suppose |$u$| is QC along every |$d\in \mathcal{D}$|⁠. Then, we have the following approximate quasiconvexity estimate: \[ -\lambda_{\rm QC}[u]\leq\mathcal{O}({\rm d}\theta). \] Proof. Without loss of generality, assume |$x=0$| and |$u(0)=0$|⁠. Suppose |$u-\phi$| attains a global maximum at |$x=0$|⁠. We may assume |$\phi$| is locally quadratic, |$\phi(x)=x^\intercal A x+b^\intercal x$|⁠, for some real, symmetric matrix |$A$|⁠, and for |$b\in{\mathbb{R}^n}$|⁠. Next, suppose |$w$| is an arbitrary unit vector satisfying |$\nabla \phi(0)^\intercal w=b^\intercal w=0$|⁠. Decompose |$w$| as follows: \[ w = \cos(\theta)v + \sin(\theta) p, \] where |$v\equiv\text{argmin}_i(\cos^{-1}(w^\intercal v_i))$|⁠, |$\theta = \cos^{-1}(w^\intercal v)$| and |$p$| is some unit vector orthogonal to |$v$|⁠. By hypothesis, |$\theta\le\pi/4$|⁠. Taking |$\lambda_n$| to denote the largest eigenvalue of |$A$|⁠, we observe \begin{align*} 0=u(0) &= u\left(\tfrac{1}{2}(-\cos(\theta) v)+\tfrac{1}{2}(\cos(\theta) v)\right)\\ &\leq \max(u(-\cos(\theta) v),u(\cos(\theta) v))\quad\text{(by quasiconvexity)}\\ &\leq \max(\phi(-\cos(\theta) v),\phi(\cos(\theta) v))\\ &= \max(\phi(- w+\sin(\theta)p),\phi( w-\sin(\theta)p))\\ &=w^\intercal Aw+\sin^2(\theta)p^\intercal Ap-2\sin(\theta)w^\intercal Ap+\sin(\theta)|p^\intercal b|\\ &\leq w^\intercal Aw+ \sin^2(\theta)\lambda_n + 2\sin(\theta)\lambda_n + \sin(\theta){\vert{b}\vert}\\ &\leq w^\intercal Aw+ \sin^2({\rm d}\theta) \lambda_n + 2\sin({\rm d}\theta)\lambda_n+ \sin({\rm d}\theta){\vert{b}\vert}\\ &=w^\intercal Aw + \mathcal{O}({\rm d}\theta). \end{align*} In the above calculations, we used the fact that |$y^\intercal Ay\leq\lambda_n$| for every |$y$|⁠, |$\phi\ge u$| in |$\overline\varOmega$| and the fact that |$b^\intercal w=0$|⁠. Taking the minimum over all unit vectors |$w$| satisfying |$b^\intercal w=0$| yields the desired result. □ The result we present next is a follow-up to a formal result discussed in Oberman (2007) regarding the similar notion of a uniformly convex envelope. We show that by enforcing a degree of strict convexity along a prescribed set of directions, one can recover a convex function on the entire domain. The subsequent and analogous result for quasiconvexity is natural and requires only a slight modification of the proof. Before proceeding, we need to state the following equivalences, found in Cannarsa & Sinestrari (2004), which will provide a convenient framework for the proof. Lemma 3.7 (Strict convexity and semiconcavity; Cannarsa & Sinestrari, 2004). Suppose |$u\in C(\overline\varOmega)$| and |$C\ge0$|⁠. The following three conditions are equivalent for strict convexity: (i) |$D^2u(x) \ge CI$| in the viscosity sense; (ii) |$u - C\frac{|\cdot|^2}{2}$| is convex; (iii) |$u(x+h) - 2u(x) + u(x-h) \geq C\frac{h^2}2$| for all |$x\in\overline\varOmega$| and |$h>0$|⁠. The following three conditions are equivalent for semiconcavity: (i) |$D^2u(x) \le CI$| in the viscosity sense; (ii) |$u - C\frac{|\cdot|^2}{2}$| is concave; (ii) |$u(x+h) - 2u(x) + u(x-h) \leq C\frac{h^2}2$| for all |$x\in\overline\varOmega$| and |$h>0$|⁠. Proposition 3.8 (Recovering convexity). Suppose |$u\in C(\overline\varOmega)$| and |$\mathcal{D}$| is a set of directions with directional resolution |${\rm d}\theta\le\pi/4$|⁠. Suppose there exists |$C>0$| such that |$u$| satisfies the following in the viscosity sense for every |$x\in\overline\varOmega$|⁠: (i) |$D^2u(x)\le CI$|⁠; (ii) |$v^\intercal D^2u(x) v \ge C\,{\rm d}\theta^2$| for every |$v\in\mathcal{D}$|⁠. Then |$\lambda_1[u]\ge0$|⁠. Proof. Suppose |$y$| is an arbitrary unit vector and |$x\in\overline\varOmega$|⁠. Without loss of generality we can assume that |$x=0$|⁠. It suffices to show \[ \frac{u(hy)+u(-hy)}{2}-u(0) \ge 0\quad\text{for every} h>0. \] We may assume |$hy$| lies between some |$d_1,d_2$| such that |$d_i/{\vert{d_i}\vert}\in\mathcal{D}$| for each |$i$|⁠, as shown in Fig. 4. Then, by the first hypothesis and Lemma 3.7, we can write \begin{equation} \frac{u(d_i)+u(-d_i)}{2} - u(0) \ge \frac{C\,{\rm d}\theta^2{\vert{d_i}\vert}^2}{2},\quad i=1,2. \end{equation} (3.3) Next, define |$h_i\equiv {\vert{d_i-y}\vert}$| and |$e_i \equiv \frac{d_i-y}{h}$| for |$i=1,2$|⁠. Then we have the following inequalities by the second hypothesis and Lemma 3.7: \begin{align} \frac{h_1}{h_1+h_2}u(y+h_2e_2) + \frac{h_2}{h_1+h_2}u(y-h_1e_1) - u(hy) &\le C\frac{h_1h_2}{2}, \\[-2pt] \end{align} (3.4) \begin{align} \frac{h_1}{h_1+h_2}u(-y-h_2e_2) + \frac{h_2}{h_1+h_2}u(-y+h_1e_1) - u(-hy) &\le C\frac{h_1h_2}{2}. \end{align} (3.5) Averaging (3.4) and (3.5), and substituting |$d_i = y+he_i$| yields \begin{align} \frac{h_2}{h_1+h_2}\left(\frac{u(d_1)+u(-d_1)}{2}\right)+\frac{h_1}{h_1+h_2}\left(\frac{u(d_2)+u(-d_2)}{2}\right) -\frac{u(hy)+u(-hy)}{2} \le C\frac{h_1h_2}{2}. \end{align} (3.6) Now, weighting each equation in (3.3) by |$h_i/(h_1+h_2)$| (respectively) and averaging them yields \begin{align} & \frac{h_2}{h_1+h_2}\left(\frac{u(d_1)+u(-d_1)}{2}\right) + \frac{h_1}{h_1+h_2}\left(\frac{u(d_2)+u(-d_2)}{2}\right) - u(0)\nonumber\\[-2pt] &\quad{} \ge \frac{C\,{\rm d}\theta^2}{2}\left( \frac{h_1}{h_1+h_2}{\vert{d_1}\vert}^2+ \frac{h_2}{h_1+h_2}{\vert{d_2}\vert}^2\right)\!. \end{align} (3.7) Subtracting (3.6) from (3.7) recovers \[ \frac{u(hy)+u(-hy)}{2}-u(0) \ge \frac{C\,{\rm d}\theta^2}{2}\left( \frac{h_1}{h_1+h_2}{\vert{d_1}\vert}^2+ \frac{h_2}{h_1+h_2}{\vert{d_2}\vert}^2\right)-C\frac{h_1h_2}{2}. \] Defining |$\theta_i\equiv \cos^{-1}(y^\intercal d_i)$|⁠, the worst-case scenario occurs when |$\theta_1=\theta_2$|⁠, so that |${\vert{d_1}\vert}={\vert{d_2}\vert}\equiv{\vert{d}\vert}$| and |$h_1=h_2\equiv h_0$|⁠. Then the above inequality simplifies to \[ \frac{u(hy)+u(-hy)}{2}-u(0) \ge \frac{C}{2}\left({\vert{d}\vert}^2\,{\rm d}\theta^2 - h_0^2\right)\!. \] Now, referring to Fig. 4, we have |$h_0 = {\vert{d}\vert}\sin({\rm d}\theta)$|⁠, so that we get \[ \frac{u(hy)+u(-hy)}{2}-u(0) \ge \frac{C{\vert{d}\vert}^2}{2}\left({\rm d}\theta^2 - \sin({\rm d}\theta)^2\right)\!, \] yielding the desired result since |${\rm d}\theta\ge\sin({\rm d}\theta)$|⁠. □ Fig. 4. View largeDownload slide Illustration of the proof of recovering convexity. Fig. 4. View largeDownload slide Illustration of the proof of recovering convexity. Proposition 3.9 (Recovering quasiconvexity). Suppose |$u\in C(\overline\varOmega)$|⁠, |$\alpha\in\mathbb{R}$| and |$\mathcal{D}$| is a set of directions with directional resolution |${\rm d}\theta\le\pi/4$|⁠. Suppose there exists |$C>0$| such that |$u$| satisfies the following in the viscosity sense for every |$x\in\overline\varOmega$|⁠: (i) |$D^2u\le CI$|⁠; (ii) |${\vert{\!{\nabla} u}\vert}\le M$|⁠; (iii) |$v^\intercal D^2u v - \alpha \ge C\,{\rm d}\theta^2$| whenever |$v\in\mathcal{D}$| satisfies |${\vert{\!{\nabla} u^\intercal v}\vert} \le M\,{\rm d}\theta$|⁠. Then |$\lambda_{\text QC}[u]\ge\alpha$|⁠. Remark 3.10. Note that the third requirement is for |$u$| to be |$(M\,{\rm d}\theta)$| robustly QC along directions in |$\mathcal{D}$|⁠. Recalling the inequalities (3.1), this can be enforced numerically by solving for subsolutions of the following: \[ \alpha+C\,{\rm d}\theta^2+\sqrt{M\,{\rm d}\theta} - F^{\sqrt{M\,{\rm d}\theta}}[u] = 0. \] Proof. We begin by taking an arbitrary unit vector |$y$| satisfying |${\nabla} \!u^\intercal y=0$|⁠. As before, we would like to show \[ \frac{u(hy)+u(-hy)}{2}-u(0) \ge 0\quad\text{for every}\,h>0. \] The proof that follows is nearly identical to the proof of convexity. The only difference arises when trying to derive equations (3.3). In this case, we have to show |${\nabla} \!u^\intercal \left(\frac{d_i}{{\vert{d_i}\vert}}\right)\le M{\rm d}\theta$| to invoke the proper inequality. This is given by the following argument for each |$i$|⁠: \begin{equation} {\nabla}\!u^\intercal \left(\frac{d_i}{{\vert{d_i}\vert}}\right) = {\nabla}\!{u}^\intercal \left(\frac{y+h_ie_i}{{\vert{d_i}\vert}}\right) = \frac{h_i}{{\vert{d_i}\vert}}{\nabla}\! u ^\intercal e_i \le \frac{h_i}{{\vert{d_i}\vert}}{\vert{\!{\nabla}\! u}\vert}{\vert{e_i}\vert}\le \frac{h_i}{{\vert{d_i}\vert}}M, \end{equation} (3.8) but |$\frac{h_i}{{\vert{d_i}\vert}}=\sin(\theta_i)\le\sin({\rm d}\theta)\le {\rm d}\theta$| for small |${\rm d}\theta$|⁠. Therefore, we have \[ {\nabla}\!{u}^\intercal \left(\frac{d_i}{{\vert{d_i}\vert}}\right) \le M\,{\rm d}\theta. \] Then, by the third hypothesis, we are now in a position to write \[ \frac{u(d_i)+u(-d_i)}{2} - u(x) \ge \frac{C{\vert{d_i}\vert}^2\theta^2}{2}. \] Proceeding identically as in the proof of Proposition 3.8 yields the desired result. □ 4. Convergent finite difference schemes In what follows, we provide numerical schemes that discretize the above PDEs. As we will see, the schemes we present here fall into a general class of degenerate elliptic finite difference schemes for which there exists a convergence framework. Before we begin, we introduce some notation, which we will carry throughout the rest of the article. We will assume we are working on the hypercube |$[-1,1]^n\subset{\mathbb{R}^n}$|⁠. We write |$x=(x_1,\dots,x_n)\in[-1,1]^n$|⁠. For simplicity, we discretize the domain |$[-1,1]^n$| with a uniform grid, resulting in the following spatial resolution: \[ h\equiv\frac{2}{N-1}, \] where |$N$| is the number of grid points used to discretize |$[-1,1]$|⁠. Note that we use |$h$| here to denote the spatial resolution and it is not to be confused with the |$h$| in the formulations of (ϵ-QC) and (ϵ-QCE). We use the following notation for our computational domain: \[ \varOmega^h \equiv [-1,1]^n\cap h\mathbb{Z}^n. \] We define a grid vector, |$v$|⁠, as \[ v \equiv \frac{x-y}{h},\quad x,y\in\varOmega^h. \] 4.1 Finite difference equations and wide stencils Consider a degenerate elliptic PDE |$F[u]=0$|⁠, and its corresponding finite difference equation, \[ F_\rho[u]=0,\quad x\in\varOmega^h, \] where |$\rho$| is the discretization parameter. In our case, we take |$\rho=(h,\,{\rm d}\theta)$| for the schemes presented hereafter. In general, |$F_\rho[u]:C(\varOmega^h)\to C(\varOmega^h)$|⁠, where |$C(\varOmega^h)$| is the set of grid functions |$u:\varOmega^h\to\mathbb{R}$|⁠. We assume the following form: \[ F_\rho[u](x) = F_\rho(u(x),u(x)-u(\cdot))\quad\ \text{for}\ x\in\varOmega^h, \] where |$u(\cdot)$| corresponds to the value of |$u$| at points in |$\varOmega^h$|⁠. For a given set of grid vectors |$V$|⁠, we say |$F_\rho[u]$| has stencil width|$W$| if for any |$x\in\varOmega^h$|⁠, |$F_\rho[u](x)$| depends only on the values |$u(x + hv)$|⁠, where |$v\in V$| and |$\max_{v\in V}{\Vert{v}\Vert}_\infty\le W$|⁠. We assume stencils are symmetric about the reference point |$x$|⁠, and have at least width 1. Example 4.1. In two dimensions, the directional resolution can be written as |${\rm d}\theta=\arctan(1/W)/2$|⁠. A stencil of width |$W=1$| would correspond to the vectors |$V=\{(0,\pm1), (\pm1,0), (\pm1,\pm1), (\mp1,\pm1)\}$|⁠. The Dirichlet boundary conditions given by (Dir) are incorporated by setting \[ F_\rho[u](x) = u(x)-g(x)\quad\ \text{for}\ x\in\partial\varOmega^h. \] The PDEs we discuss involve a second-order directional derivative and a directional Eikonal operator. We use the following standard discretizations for our schemes: \begin{align*} {\vert{u_v^h(x)}\vert} &\equiv \max\left\{\frac{u(x \pm hv)-u(x)}{h}\right\}\!,\\ u_{vv}^h(x) &\equiv \frac{u(x+hv)+u(x-hv)-2u(x)}{h^2}, \end{align*} where |$u_v^h$| and |$u_{vv}^h$| denote the finite difference approximations for the first- and second-order directional derivatives in the grid direction |$v$|⁠. Recall that for any twice continuously differentiable function |$u$|⁠, the standard Taylor series computation yields the following consistency estimate: \begin{align} {\vert{u_v^h(x)}\vert} &= {\vert{\nabla u(x)^\intercal v}\vert} + \mathcal{O}(h), \\ \end{align} (4.1) \begin{align} u_{vv}^h(x) &= v^\intercal D^2u(x) v + \mathcal{O}(h^2). \end{align} (4.2) 4.2 A finite difference method for the PDE We begin by considering the full discretization of |$F^\epsilon[u]$|⁠, denoted |$F^\epsilon_{h,{\rm d}\theta}[u]$|⁠, where |$h$| is the spatial resolution and |${\rm d}\theta$| is the directional resolution of the set of grid vectors |$V$|⁠. We use the spatial discretizations given by (4.1) and (4.2): \begin{equation} F^\epsilon_{h,{\rm d}\theta}[u](x)\equiv \min_{v\in V}\left\{\tfrac{1}{\epsilon}{\vert{u_v^h(x)}\vert}+u_{vv}^h(x)\right\}\quad \text{for}\,x\in\varOmega^h. \end{equation} (4.3) We write the full discretization of (ϵ-Ob) as \begin{equation} \begin{cases} \max\{u(x)-g(x),\epsilon - F^\epsilon_{h,{\rm d}\theta}[u](x)\} = 0, & x\in\varOmega^h,\\ u(x) = g(x), & x\in\partial\varOmega^h, \end{cases} \end{equation} (4.4) where we make the following abbreviation: \[ G^\epsilon_{h,{\rm d}\theta}[u]\equiv\max\left\{u-g,\epsilon-F^\epsilon_{h,{\rm d}\theta}[u]\right.\}. \] 4.3 Iterative solution method We implement a fixed-point solver to recover the numerical solution of (4.4). The method is equivalent to a forward Euler time discretization of the parabolic version of (ϵ-QCE), \[ u_{t} + G^\epsilon[u]=0\quad \text{in}\,\varOmega, \] along with |$u = g$| on |$\partial\varOmega$| and |$u = u_0$| for |$t=0$|⁠. In particular, the following iterations are performed until a steady state is reached: \[ \begin{cases} u^{n+1} = u^{n} - \delta G^\epsilon_{h,d\theta}[u^{n}], & n\ge0,\\ u^{0} = u_0, \end{cases} \] where satisfies the CFL condition |$\delta\le 1/K^{h,\epsilon}$| (see Oberman, 2006) and |$K^{h,\epsilon}$| is the Lipschitz constant of the scheme. For our scheme, |$K$| is given by \[ K^{h,\epsilon} = \frac{1}{\epsilon h} + \frac{1}{h^2}, \] where the first and second terms come from the discretizations of the first- and second-order directional derivatives, respectively. 4.4 Using a first-order discretization In higher dimensions, we need the second-order term in (ϵ-QC) to guarantee quasiconvexity of the solutions. Interestingly, however, for a careful choice of |$\epsilon$|⁠, the discretization of only the first-order term is consistent with |$F^\epsilon[u]$|⁠. Indeed, observe that \begin{align*} -\min_{v\in V}\frac{{\vert{u_v^h}\vert}}{\epsilon}&=-\min_{v\in V}\left\{\frac {\max\{u(x\pm hv)-u(x)\}} {\epsilon h}\right\}\\ &= -\min_{v\in V}\left\{ \frac{{\vert{\!{\nabla}\! u(x)^\intercal v}\vert}}{\epsilon} + \frac{h}{2\epsilon}v^\intercal D^2u(x)v \right\}+ \mathcal{O}\left(\frac{h^2}{\epsilon}\right). \end{align*} Choosing |$\epsilon = h/2$| results in \[ -\min_{v\in V}\frac{{\vert{u_v^h}\vert}}{\epsilon} = -\min_{v\in V}\left\{ \frac{{\vert{\!{\nabla}\! u(x)^\intercal v}\vert}}{h/2} + v^\intercal D^2u(x)v \right\}+ \mathcal{O}(h). \] Thus, for |$\epsilon=h/2$|⁠, we see for fixed |$(h,{\rm d}\theta)$| that the discretization of the first-order term approximates |$F^\epsilon[u]$|⁠. Moreover, the time step in the above forward Euler discretization simplifies to |$\delta=\mathcal{O}(h^2)$|⁠. 5. Convergence of numerical solutions In this section we introduce the notion of degenerate elliptic schemes and show that the solutions of the proposed numerical schemes converge to the solutions of (ϵ-Ob) as the discretization parameters tend to 0. The standard framework used to establish convergence is that of Barles & Souganidis (1991), which we state below. In particular, it guarantees that the solutions of any monotone, consistent and stable scheme converge to the unique viscosity solution of the PDE. 5.1 Degenerate elliptic schemes Consider the Dirichlet problem for the degenerate elliptic PDE |$F[u]=0$|⁠, and recall its corresponding finite difference formulation: \[ \begin{cases} F_\rho(u(x),u(x)-u(\cdot)) = 0, & x\in\varOmega^h,\\ u(x)-g(x) = 0, & x\in\partial\varOmega^h, \end{cases} \] where |$\rho$| is the discretization parameter. Definition 5.1. The PDE |$F_\rho[u]$| is a degenerate elliptic scheme if it is nondecreasing in each of its arguments. Remark 5.2. Although the convergence theory is originally stated in terms of monotone approximation schemes (schemes with non-negative coefficients), ellipticity is an equivalent formulation for finite difference operators (Oberman, 2006). In fact, degenerate ellipticity implies monotonicity and stability in the nonexpansive norm. Definition 5.3. The finite difference operator |$F_\rho[u]$| is consistent with |$F[u]$| if for any smooth function |$\phi$| and |$x\in\varOmega$| we have \[ F_\rho[\phi](x)\to F[\phi](x)\quad \text{as}\,\rho\to0. \] Definition 5.4. The finite difference operator |$F_\rho[u]$| is stable if there exists |$M>0$| independent of |$\rho$| such that if |$F_\rho[u]=0$| then |$\| u \|_\infty\le M$|⁠. Remark 5.5. (Interpolating to the entire domain). The convergence theory assumes that the approximation scheme and the grid function are defined on all of |$\varOmega$|⁠. Although the finite difference operator acts only on functions defined on |$\varOmega^h$|⁠, we can extend such functions to |$\varOmega\setminus{\it{\Omega}}^h$| via piecewise interpolation. In particular, performing piecewise linear interpolation maintains the ellipticity of the scheme as well as all other relevant properties. Therefore, we can safely interchange |${\it{\Omega}}^h$| and |$\varOmega$| in the discussion of convergence without any loss of generality. 5.2 Convergence of numerical approximations Next we will state the theorem for convergence of approximation schemes, tailored to elliptic finite difference schemes, and demonstrate that the proposed schemes fit in the desired framework. In particular, we will show that the schemes are elliptic, consistent and have stable solutions. Proposition 5.6 (Convergence of approximation schemes; Barles & Souganidis, 1991). Consider the degenerate elliptic PDE |$F[u]=0$|⁠, with Dirichlet boundary conditions for which there exist a strong comparison principle. Let |$F_\rho[u]$| be a consistent and elliptic scheme. Furthermore, assume that the solutions of |$F_\rho[u]=0$| are bounded independently of |$\rho$|⁠. Then |$u^\rho\to u$| locally uniformly on |$\varOmega$| as |$\rho\to0$|⁠. Remark 5.7. Proposition 5.6 is an instance of the meta-theorem, which says that consistency and stability imply convergence. The proposition can also be seen as a special case of the fact that viscosity solutions are stable under perturbations of both the data and the operators, provided that the equation remains in the class of degenerate elliptic PDEs. The classical example of a perturbation is to add the Laplacian (or viscosity) to the operator: |$F^\epsilon = F + \epsilon {\it{\Delta}}$|⁠, which is the origin of the nomenclature ‘viscosity solutions’. The finite difference approximation of |$F$| can also be regarded as a perturbation of the operator. Provided that the approximation is (i) consistent and (ii) degenerate elliptic (Oberman, 2006) or monotone and stable (Barles & Souganidis, 1991) then the perturbation converges. Lemma 5.8. (Ellipticity). The scheme given by (4.4) is elliptic. Proof. The negatives of each of the spatial discretizations above are elliptic for fixed |$v$|⁠, so taking the minimum over all |$v$| maintains the ellipticity. Therefore, the discretization of |$F^\epsilon_{h,{\rm d}\theta}[u]$| is elliptic. Moreover, the obstacle term is trivially elliptic. Hence |$G^\epsilon_{h,{\rm d}\theta}[u]$|⁠, being the maximum of two elliptic schemes, is also elliptic. □ Lemma 5.9 (Consistency). Given a smooth function |$u$| and a set of grid vectors |$V$| with directional resolution |${\rm d}\theta\le\pi/4$|⁠, we have the following estimate: \begin{equation} G^\epsilon_{h,{\rm d}\theta}[u] - G^\epsilon[u] = \mathcal{O}(h+{\rm d}\theta). \end{equation} (5.1) In other words, the scheme (4.4) is consistent with (ϵ-Ob). Proof. It suffices to show that |$F^\epsilon_{h,{\rm d}\theta}[u] - F^\epsilon[u] = \mathcal{O}(h+{\rm d}\theta)$|⁠. Fix |$x\in\overline\varOmega$| and let |$u$| be a smooth function on |$\overline{\varOmega}$|⁠. Define \[ w\equiv{\text{argmin}}_{{\vert{v}\vert}=1}\left\{ \tfrac{1}{\epsilon}{\vert{\nabla u(x)^\intercal v}\vert} + v^\intercal D^2u(x) v\right\}. \] Let |$v\in V$| and make the following decomposition: \[ v = \cos(\theta)w + \sin(\theta)p, \] for some unit vector |$p$| orthogonal to |$w$|⁠. Performing a similar calculation to the proof in Proposition 3.6, we obtain the following consistency estimate from the directional resolution: \[ F^\epsilon[u] \le\min_{v\in V}\left\{\tfrac{1}{\epsilon}{\vert{\nabla u^\intercal v}\vert} + v^\intercal D^2u v\right\} \le F^\epsilon[u] +\mathcal{O}({\rm d}\theta). \] Finally, we conclude (5.1) by recalling the standard Taylor series consistency estimate from equations (4.1) and (4.2). □ Next, we establish stability of our numerical solutions by applying a discrete comparison principle for strict subsolution and supersolution. In particular, we use the following result. Proposition 5.10. (Discrete comparison principle for strict subsolutions). Let |$F_\rho[u]$| be a degenerate elliptic scheme defined on |$\varOmega^h$|⁠. Then for any grid functions |$u,v$| we have \[ F_\rho[u] 0. \] Then, \[ u(x_0) - u(x) > v(x_0) - v(x)\quad\text{for every}\,x\in\varOmega^h, \] and so by ellipticity (Lemma 5.8), \[ F_\rho[u](x_0)\ge F_\rho[v](x_0), \] which is the desired contradiction. □ Lemma 5.11 (Stability). The numerical solutions of (4.4) are stable. Proof. This a direct application of Proposition 5.10. Indeed, suppose |$u$| is a numerical solution such that |$G^\epsilon_{h,{\rm d}\theta}[u]=0$|⁠. Define |$g_1(x)\equiv C+2\epsilon^2 \sum_{i=1}^nx_i$|⁠, where |$C\in\mathbb{R}$| is chosen so that |$g_1(x)