A partial differential equation for the $$\epsilon$$ -uniformly quasiconvex envelope

A partial differential equation for the $$\epsilon$$ -uniformly quasiconvex envelope 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 $$v<u$$ in $$\varOmega$$). This contradicts the uniqueness of (Ob). It is shown in Barron et al. (2013) that continuous QC functions necessarily satisfy $$\lambda_{\rm QC}[u] \geq 0$$ in the viscosity sense. The converse, however, is not true: consider the function $$u(x)=-x^4$$, which is not QC yet satisfies $$\lambda_{\rm QC}[u]=0$$ at $$x=0$$. In fact, using this as a sufficient condition is only possible when $$u$$ has no local maxima (Barron et al., 2013). This operator can be used to completely characterize the set of continuous $$\epsilon$$-robustly QC functions, defined below. Definition 1.6 ($$\epsilon$$-robustly QC) The function $$u:\varOmega\to\mathbb{R}$$ is $$\epsilon$$-robustly QC if $$u(x) +y^\intercal x$$ is QC for every $${\vert{y}\vert}\leq\epsilon$$. In particular, $$\epsilon$$-robustly QC functions are functions whose quasiconvexity is maintained under small linear perturbations. Write $$\lambda_{\rm QC}^\epsilon[u](x) \equiv L_\epsilon({\nabla} u(x),D^2u(x))$$. Proposition 1.7 (Characterization of $$\epsilon$$-robustly QC functions; Barron et al., 2013). The upper semicontinuous function, $$u:\varOmega\to\mathbb{R}$$, is $$\epsilon$$-robustly QC if and only if $$u$$ is a viscosity subsolution of $$\lambda_{\rm QC}^\epsilon[u]=0$$. Next, we show that viscosity subsolutions (defined below) of our operator are robustly-QC (which also implies they are QC). This allows us to avoid a technical argument from Barron et al. (2013). We believe that stronger results hold: see the formal analysis in Section 3. For simplicity, we first define the following abbreviation. Definition 1.8 ($$\epsilon$$-uniformly QC). The function $$u\in {\rm USC}(\overline\varOmega)$$ is $$\epsilon$$-uniformly QC if it is a viscosity subsolution of $$\epsilon-F^\epsilon[u]=0$$. Proposition 1.9. Suppose $$u\in {\rm USC}(\overline\varOmega)$$ is $$\epsilon$$-uniformly QC. Then, $$u$$ is $$\epsilon^2$$-robustly QC. Proof. First observe that for any $$\phi\in C^2$$ we have the following inequality: \begin{align*} F^\epsilon[\phi](x) - \frac{\alpha}{\epsilon} &= \min_{{\vert{v}\vert}=1}\left\{\tfrac{1}{{\epsilon}}{\vert{\nabla \phi(x)\cdot v}\vert} +v^\intercal D^2\phi(x)v - \tfrac{\alpha}{\epsilon} \right\}\\ &\leq \min_{\{{\vert{v}\vert}=1,{\vert{\nabla \phi(x)\cdot v}\vert}\le \alpha \}}\left\{v^\intercal D^2\phi(x)v\right\} = \lambda_{QC}^{\alpha}[\phi](x). \end{align*} Choosing $$\alpha=\epsilon^2$$ we see that any viscosity subsolution of $$\epsilon-F^\epsilon[u]= 0$$ is also a viscosity subsolution of $$-\lambda_{\rm QC}^{\epsilon^2}[u]= 0$$. Thus, $$u$$ is $$\epsilon^2$$-robustly QC. □ 2. Properties of solutions In this section, we present technical arguments proving the uniqueness of solutions of our PDEs and discuss some relevant properties. 2.1 Comparison principle In this subsection, we will show that a weak comparison principle holds for the Dirichlet problem of $$h-F^\epsilon[u]=0$$ for $$h>0$$ 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)<g(x)\}$$. Then, $$h-F^\epsilon[v]\ge 0$$ in $$\varOmega_g$$ and $$v=g$$ on $$\partial\varOmega_g$$, that is, $$v$$ is a viscosity supersolution of $$h-F^\epsilon[v]=0$$ in $$\varOmega_g$$. Now, by the definition of viscosity subsolutions we have $$u\le g$$ and $$h-F^\epsilon[u]\le 0$$ in $$\varOmega$$ and thus also $$\varOmega_g$$. This allows us to conclude that $$u\le v$$ on $$\partial\varOmega_g$$. Therefore, by Proposition 2.4, $$u\le v$$ in $$\varOmega_g$$. Concluding, in $$\varOmega\setminus\varOmega_g$$ we necessarily have $$v\ge g\ge u$$. □ 2.2 The $$\epsilon$$-uniformly QCE We formulate the $$\epsilon$$-uniformly QCE of a function $$g$$ as the unique viscosity solution of the following obstacle problem: ($$\epsilon$$-Ob) \begin{equation} \begin{cases} \max\{u(x)-g(x),\epsilon-F^\epsilon[u](x)\} = 0, & x\in\varOmega,\\ u(x) = g(x), & x\in\partial\varOmega. \end{cases} \end{equation} Remark 2.6 (Convergence of approximate solutions). It is clear that as $$\epsilon\to0$$, the penalization term in $$F^\epsilon[u]$$ tends to infinity. The result is that $$F^\epsilon[u]\to\lambda_{\rm QC}[u]$$ as $$\epsilon\to0$$. From this observation, the standard stability result of viscosity solutions and the quasiconvexity of subsolutions of $$\epsilon-F^\epsilon[u]=0$$, one can then apply the same argument presented in the proof of Barron et al. (2013, Theorem 5.3) to conclude that the unique viscosity solutions of (ϵ-QCE), (Dir) converge to the QCE as $$\epsilon\to0$$. This result allows us to compute asymptotic approximations of the QCE of a given obstacle. 2.3 Solution formula in one dimension In one dimension, $$\epsilon-F^\epsilon[u]$$ is simply the Eikonal operator with a small diffusion term. When considering a solution $$u$$ of (ϵ-Ob), whenever $$u(x) < g(x)$$, we have $$F^\epsilon[u](x) = \epsilon$$, which gives \begin{equation} \epsilon u_{xx} + {\vert{u_x}\vert}=\epsilon^2. \end{equation} (2.3) Define $$\varOmega_g\equiv\{x:u(x) < g(x)\}$$. Note that $$\varOmega_g$$ need not be connected; it can be written as the union of finitely many intervals (refer to Fig. 3). However, we can solve the equation in each interval. Lemma 2.7 (One-dimensional solution). The viscosity solution of the one-dimensional PDE for the operator (2.1) and (2.3), along with boundary conditions $$u(0)=0, u(W) = H$$, is given by the following. Set $$S\equiv H/\epsilon^2W$$. Then, \[ u^\epsilon(x) = \begin{cases} \pm \epsilon^2 x, & S=\pm1,\\ \epsilon^2x + C^+(1-\exp(-x/\epsilon)), & S>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 $$x<x_0$$ and $$u_x>0$$ 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]<F_\rho[v]\implies u\le v\quad\text{in}\,\varOmega^h. \] Proof. Suppose for contradiction that \[ \max_{x\in{\it{\Omega}}^h}\{u(x) - v(x)\} = u(x_0) - v(x_0) > 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)<g(x)$$ for every $$x\in\overline\varOmega$$. Moreover, one can check that $$(g_1)_{vv}^h(x)=0$$. Then, we observe \begin{align*} G^\epsilon_{h,{\rm d}\theta}[g_1](x) &= \max\left\{g_1(x)-g(x),\epsilon-\min_{v\in V}\left\{ \frac{\max\{g_1(x\pm hv)\}-g_1(x)}{\epsilon h} \right\}\right\}\\ &= \max\left\{g_1(x)-g(x),\epsilon-2\epsilon\min_{v\in V}\left\vert\sum\nolimits_{i=1}^nv_i\right\vert\right\}\\ &= \max\left\{g_1(x)-g(x),-\epsilon\right\}, \end{align*} where the last equality follows from the general assumption that the stencil width is at least 1 with directional resolution $${\rm d}\theta\le\pi/4$$. Next, define $$g_2(x)\equiv\max_{x\in\overline{\it{\Omega}}}g(x)$$. Then, in summary we have $$G^\epsilon_{h,{\rm d}\theta}[g_1]<0$$ and $$G^\epsilon_{h,{\rm d}\theta}[g_2]\ge0$$. Therefore, by Proposition 5.10, \[ \min_{x\in\overline\varOmega} g_1(x)\le u(x)\le g_2 \quad\text{for}\ x\in\varOmega^h. \] Moreover, this bound is independent of $$h$$ and $${\rm d}\theta$$. □ Proposition 5.12. The solutions of (4.4) converge locally uniformly on $$\varOmega$$ to the solutions of (ϵ-Ob) as $$h,{\rm d}\theta\to0$$. Proof. By Lemmas 5.8, 5.9, 5.11 we see that (4.4) is a consistent and elliptic scheme which has stable solutions. Moreover, a weak comparison principle holds by Corollary 2.5. Therefore, by Proposition 5.6, the numerical solutions converge uniformly on compact subsets of $$\varOmega$$. □ 6. Numerical results In this section we present the numerical results. The tolerance for the fixed-point iterations was taken to be $$10^{-6}$$ and unless otherwise stated we set $$\epsilon=h/2$$, where $$h$$ is the spatial resolution of the grid. Technical conditions on $$g$$ given by Remark 2.1 are needed to ensure that the boundary conditions are held in the strong sense. In practice, we violated this condition, resulting in solutions that were discontinuous at the boundary. Two natural choices for initialization of the solution are (i) the obstacle function, $$g(x)$$ and (ii) a QC function below $$g$$. We found that the first choice leads to very slow convergence. In particular, the parabolic equation takes time $$\mathcal{O}(1/\epsilon)$$ to converge. See Example 6.3 below. On the other hand, the second choice results in faster convergence: the simplest choice is simply the constant function with value the minimum of $$g$$. After one step of the iteration the boundary values are attained. Moreover, starting with this choice allows us to use the iterative method with $$\epsilon = 0$$ to find the QCE. Remark 6.1 It is an open question whether solving the $$\epsilon=0$$ time-dependent PDE with quasiconvex initial data below $$g$$ will converge. To apply our convergence results, we take $$\epsilon = h/2$$. 6.1 Results in one dimension We present examples demonstrating the convergence of approximate QCEs to the true QCE as $$\epsilon\to0^+$$. We also compare the iteration count when starting from below and at the obstacle. Example 6.2 (Convergence as $$\epsilon\to0$$). We consider the convergence of numerical solutions of (4.4) as $$\epsilon\to 0$$ with the following obstacle function: \[ g(x)=\min\{{\vert{x-0.5}\vert},{\vert{x+0.5}\vert}-0.3\}. \] The results are displayed in Fig. 5. Indeed, as expected we see convergence to the true QCE, from below, as $$\epsilon\to0$$. Fig. 5. View largeDownload slide Obstacles (solid) and the numerical solutions (dashed). Left: Example 6.2. Convergence in $$\epsilon$$. Bottom to top: $$\epsilon = 0.2\,,0.1\,,0.05\,,0.001\,,0.0001$$. Right: Example 6.3. Fixed-point iterations when $$u^0 \equiv g_2$$. Fig. 5. View largeDownload slide Obstacles (solid) and the numerical solutions (dashed). Left: Example 6.2. Convergence in $$\epsilon$$. Bottom to top: $$\epsilon = 0.2\,,0.1\,,0.05\,,0.001\,,0.0001$$. Right: Example 6.3. Fixed-point iterations when $$u^0 \equiv g_2$$. Example 6.3 (Visualization of iterations). Next, we consider the obstacle \[ g(x)= -x^2+1, \] whose QCE is simply $$g_2^{\rm QCE}\equiv0$$. We demonstrate the evolution of the iterations when the initial data are taken to be the obstacle. In this case, the solution corresponded closely to $$u(x,t) = \min(g(x), c - \epsilon t)$$, where $$c = \max g$$. This illustrates the slow speed of convergence, and the fact that the equation degenerates to a trivial operator as $$\epsilon \to 0$$. Results are displayed in Fig. 5. 6.2 Numerical results in two dimensions All examples take the initial function in the iterative scheme to be the minimum value $$g_{\rm min} = \min g$$ of the obstacle function. In all contour plots, the solid line represents the level sets of the original function and the dashed line represents the same level sets of the numerical solution. Unless otherwise stated, the two-dimensional numerical solutions shown are computed on a $$64\times64$$ grid, for illustration purposes. Computations were performed on larger-sized grids. We performed the same numerical experiments in Abbasi & Oberman (2016), and achieved similar results. This is expected, since solving our PDE with small $$\epsilon$$ is consistent with the QCE. We do not reproduce those results here to save space. The examples we present focus on the difference between the two operators for values of $$e$$ close to 1. Example 6.4. (Uniform convexification of level sets). Let $$g(x)$$ be the signed distance function (negative on the inside) to the square, $$S = \left\{x \mid \max_{i}{\vert{x_i}\vert}=1/2 \right\}$$. Although $$g$$ is convex, its level sets are not uniformly convex. We apply the iterative procedure to $$g$$ with $$\epsilon=0.5$$ so that we uniformly convexify the level sets. The results demonstrating this are displayed in Fig. 6. Fig. 6. View largeDownload slide Left: Examples 6.4, square level sets become uniformly convex with $$\epsilon = 1/2$$. Middle: Example 6.5, the solution with $$\epsilon=1$$ is indicated by the dashed contours. Right: Example 6.6 with $$\epsilon$$ near 0. Fig. 6. View largeDownload slide Left: Examples 6.4, square level sets become uniformly convex with $$\epsilon = 1/2$$. Middle: Example 6.5, the solution with $$\epsilon=1$$ is indicated by the dashed contours. Right: Example 6.6 with $$\epsilon$$ near 0. Example 6.5 (Nonconvex signed distance function). In this example, $$g(x)$$ is the signed distance function to the curve $${\it{\Gamma}}$$, \begin{align*} {\it{\Gamma}}(t) &= \frac{1}{2}\begin{cases} (\cos(t-\frac{\pi}{2}),\sin(t-\frac{\pi}{2})), & t\in[0,3\pi/2],\\[9pt] \left(0,\frac{t-7\pi/4}{\pi/4}\right)\!, & t\in[3\pi/2,7\pi/4], \\[9pt] \left(\frac{7\pi/4-t}{\pi/4},0\right)\!, & t\in[7\pi/4,2\pi]. \end{cases} \end{align*} We compare the results of using different $$\epsilon$$. In particular, we choose $$\epsilon=h/2$$ and $$\epsilon=1$$. In the latter case, we see that the level sets are more curved. Results are displayed in Fig. 6. Example 6.6. We consider the obstacle where $$g$$ is a cone with circular portions removed from its level sets. Take \[ g(x) = \begin{cases} 1, &x\in A,\\ {\vert{x}\vert}-1/2, &x\in {\it{\Omega}}\setminus A, \end{cases} \] where $$A=\{(x_1\pm1/2)^2+x_2^2\le 1/16\}\bigcup \{x_1^2+(x_2\pm1/2)^2\le1/16\}$$. Results showing the $$g=0.7$$ level set are found in Fig. 6. Example 6.7 (Comparison to the $$\epsilon$$-robust QCE). We compare the solutions of our PDE to the solutions of the line solver presented in Abbasi & Oberman (2016), which returns the robust QCE. In particular, we use a stencil width $$W=5$$, we used the same directions in the line solver and we set $$\epsilon=0.02$$ for the $$\epsilon$$-uniformly QCE, and we also used a regularization of $$0.02$$ for the robust QC (in principle we should have used $$0.02^2$$ but the matching value was better for illustration purposes). The results are found in Fig. 7. Fig. 7. View largeDownload slide Example 6.7. Top: Surface plots (inverted) of the function, robust QCE, $$\epsilon$$-uniformly QCE. Bottom: Corresponding contour plots. Fig. 7. View largeDownload slide Example 6.7. Top: Surface plots (inverted) of the function, robust QCE, $$\epsilon$$-uniformly QCE. Bottom: Corresponding contour plots. 6.3 Accelerating iterations using the line solver We found an effective method to reduce the computational time to find the solution. We implement the line solver for the robust QCE proposed in Abbasi & Oberman (2016), alternating with the PDE iterations. In particular, on an $$n\times n$$ grid, after every $$2n$$ iterations of the PDE we apply one iteration of the line solver (with the commensurate value of $$\epsilon$$) in each direction. Table 1 shows the number of iterations required for convergence using the same obstacle from Abbasi & Oberman (2016, Example 6.1). Note the line solver is for a different regularization of the QC operator; however, for small values of $$\epsilon,$$ the solutions are close. For larger values of $$\epsilon,$$ the operator still accelerates the solver, but it does not approach the solution to within arbitrary precision. Table 1 Number of iterations required for convergence using different solution methods $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 Table 1 Number of iterations required for convergence using different solution methods $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 In Table 1, we compare the number of iterations required for each method. The results were comparable across different examples. The number of iterations of the PDE operator required for convergence was typically a small constant times $$n^2$$. The computational cost of a single line solve was on the order of a small constant (say 10) times the cost of a single PDE iteration. So the combined method, requires $$\mathcal{O}(n)$$ iterations (possibly with a $$\log n$$ prefactor, which we ignore, since it is not significant at these values of $$n$$). So the combined method required (roughly) $$n$$ iterations. Recall that there are $$N =n^2$$ grid points. Each iteration has a cost proportional to the number of directions used and to the number of grid points, so the total cost is $$\mathcal{O}\left((WN^2) \right)$$ for the iterative method and $$\mathcal{O}(WN)$$ for the combined method (with constants of roughly 1 and 10, respectively). So combining the iterative solver with the line solver results in significant improvements to computation time. References Abbasi B. & Oberman A. M. ( 2016 ) Computing the quasiconvex envelope using a nonlocal line solver. Available at https://arxiv.org/abs/1612.05584. Avriel M. , Diewert W. E. , Schaible S. , & Zang I. ( 2010 ) Generalized concavity . Society for Industrial and Applied Mathematics . Bardi M. & Mannucci P. ( 2006 ) On the Dirichlet problem for non-totally degenerate fully nonlinear elliptic equations. Commun. Pure Appl. Anal. , 5 , 709 – 731 . Google Scholar CrossRef Search ADS Bardi M. & Mannucci P. ( 2013 ) Comparison principles and Dirichlet problem for fully nonlinear degenerate equations of Monge–Ampère type. Forum Math , 25 , 1291 – 1330 . Google Scholar CrossRef Search ADS Barles G. & Souganidis P. E. ( 1991 ) Convergence of approximation schemes for fully nonlinear second order equations. Asymptotic Anal. , 4 , 271 – 283 . Barron E. N. , Goebel R. & Jensen R. R. ( 2012a ) Functions which are quasiconvex under linear perturbations. SIAM J. Optim. , 22 , 1089 – 1108 . Google Scholar CrossRef Search ADS Barron E. N. , Goebel R. & Jensen R. R. ( 2012b ) The quasiconvex envelope through first-order partial differential equations which characterize quasiconvexity of nonsmooth functions. Discrete Contin. Dyn. Syst. Ser. B , 17 , 1693 – 1706 . Google Scholar CrossRef Search ADS Barron E. , Goebel R. & Jensen R. ( 2013 ) Quasiconvex functions and nonlinear pdes. Trans. Am. Math. Soc. , 365 , 4229 – 4255 . Google Scholar CrossRef Search ADS Barron E. N. & Jensen R. R. ( 2014 ) A uniqueness result for the quasiconvex operator and first order PDEs for convex envelopes. In Annales de l’Institut Henri Poincare (C) Non Linear Analysis , vol. 31 . Elsevier Masson , pp. 203 – 215 . Google Scholar CrossRef Search ADS Boyd S. & Vandenberghe L. ( 2004 ) Convex Optimization . Cambridge : Cambridge University Press . Google Scholar CrossRef Search ADS Caffarelli L. A. & Spruck J. ( 1982 ) Convexity properties of solutions to some classical variational problems. Comm. Partial Differential Equations , 7 , 1337 – 1379 . Google Scholar CrossRef Search ADS Cannarsa P. & Sinestrari C. ( 2004 ) Semiconcave Functions, Hamilton-Jacobi Equations, and Optimal Control , vol. 58. Springer Science & Business Media . Carlier G. & Galichon A. ( 2012 ) Exponential convergence for a convexifying equation. ESAIM Control Optim. Calc. Var. , 18 , 611 – 620 . Google Scholar CrossRef Search ADS Colesanti A. & Salani P. ( 2003 ) Quasi-concave envelope of a function and convexity of level sets of solutions to elliptic equations. Math. Nachr. , 258 , 3 – 15 . Google Scholar CrossRef Search ADS Crandall M. G. , Ishii H. , & Lions P. L. ( 1992 ) User’s guide to viscosity solutions of second order partial differential equations. Bull. Amer. Math. Soc. , 27 , 1 – 67 . Google Scholar CrossRef Search ADS Froese B. ( 2016 ) Convergent approximation of surfaces of prescribed Gaussian curvature with weak Dirichlet conditions. Available at https://arxiv.org/abs/1601.06315. Kawohl B. ( 1985 ) Rearrangements and Convexity of Level Sets in PDE . Lecture Notes in Mathematics , vol. 1150 . Berlin, Heidelberg : Springer-Verlag , pp. 1 – 134 . Google Scholar CrossRef Search ADS Kohn R. & Serfaty S. ( 2007 ) Second-order PDE’s and deterministic games. In Proc. Internat. Congress Ind. Appl. Math. (ICIAM’07) , pp. 239 – 249 . Lions P.-L. ( 1998 ) Identification du cône dual des fonctions convexes et applications. C. R. Acad. Sci. Paris Sér. I Math. , 326 , 1385 – 1390 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2004 ) A convergent monotone difference scheme for motion of level sets by mean curvature. Numer. Math. , 99 , 365 – 379 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2006 ) Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton-Jacobi equations and free boundary problems. SIAM J. Numer. Anal. , 44 , 879 – 895 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2007 ) The convex envelope is the solution of a nonlinear obstacle problem. Proc. Amer. Math. Soc. , 135 , 1689 – 1694 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2008a ) Computing the convex envelope using a nonlinear partial differential equation. Math. Models Methods Appl. Sci. , 18 , 759 – 780 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2008b ) Wide stencil finite difference schemes for the elliptic Monge-Ampère equation and functions of the eigenvalues of the Hessian. Discrete Contin. Dyn. Syst. Ser. B , 10 , 221 – 238 . Google Scholar CrossRef Search ADS Sethian J. A. ( 1999 ) Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science , vol. 3 . Cambridge University Press . © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

A partial differential equation for the $$\epsilon$$ -uniformly quasiconvex envelope

Loading next page...
 
/lp/ou_press/a-partial-differential-equation-for-the-epsilon-uniformly-quasiconvex-FD0Rd99iZq
Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx068
Publisher site
See Article on Publisher Site

Abstract

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 $$v<u$$ in $$\varOmega$$). This contradicts the uniqueness of (Ob). It is shown in Barron et al. (2013) that continuous QC functions necessarily satisfy $$\lambda_{\rm QC}[u] \geq 0$$ in the viscosity sense. The converse, however, is not true: consider the function $$u(x)=-x^4$$, which is not QC yet satisfies $$\lambda_{\rm QC}[u]=0$$ at $$x=0$$. In fact, using this as a sufficient condition is only possible when $$u$$ has no local maxima (Barron et al., 2013). This operator can be used to completely characterize the set of continuous $$\epsilon$$-robustly QC functions, defined below. Definition 1.6 ($$\epsilon$$-robustly QC) The function $$u:\varOmega\to\mathbb{R}$$ is $$\epsilon$$-robustly QC if $$u(x) +y^\intercal x$$ is QC for every $${\vert{y}\vert}\leq\epsilon$$. In particular, $$\epsilon$$-robustly QC functions are functions whose quasiconvexity is maintained under small linear perturbations. Write $$\lambda_{\rm QC}^\epsilon[u](x) \equiv L_\epsilon({\nabla} u(x),D^2u(x))$$. Proposition 1.7 (Characterization of $$\epsilon$$-robustly QC functions; Barron et al., 2013). The upper semicontinuous function, $$u:\varOmega\to\mathbb{R}$$, is $$\epsilon$$-robustly QC if and only if $$u$$ is a viscosity subsolution of $$\lambda_{\rm QC}^\epsilon[u]=0$$. Next, we show that viscosity subsolutions (defined below) of our operator are robustly-QC (which also implies they are QC). This allows us to avoid a technical argument from Barron et al. (2013). We believe that stronger results hold: see the formal analysis in Section 3. For simplicity, we first define the following abbreviation. Definition 1.8 ($$\epsilon$$-uniformly QC). The function $$u\in {\rm USC}(\overline\varOmega)$$ is $$\epsilon$$-uniformly QC if it is a viscosity subsolution of $$\epsilon-F^\epsilon[u]=0$$. Proposition 1.9. Suppose $$u\in {\rm USC}(\overline\varOmega)$$ is $$\epsilon$$-uniformly QC. Then, $$u$$ is $$\epsilon^2$$-robustly QC. Proof. First observe that for any $$\phi\in C^2$$ we have the following inequality: \begin{align*} F^\epsilon[\phi](x) - \frac{\alpha}{\epsilon} &= \min_{{\vert{v}\vert}=1}\left\{\tfrac{1}{{\epsilon}}{\vert{\nabla \phi(x)\cdot v}\vert} +v^\intercal D^2\phi(x)v - \tfrac{\alpha}{\epsilon} \right\}\\ &\leq \min_{\{{\vert{v}\vert}=1,{\vert{\nabla \phi(x)\cdot v}\vert}\le \alpha \}}\left\{v^\intercal D^2\phi(x)v\right\} = \lambda_{QC}^{\alpha}[\phi](x). \end{align*} Choosing $$\alpha=\epsilon^2$$ we see that any viscosity subsolution of $$\epsilon-F^\epsilon[u]= 0$$ is also a viscosity subsolution of $$-\lambda_{\rm QC}^{\epsilon^2}[u]= 0$$. Thus, $$u$$ is $$\epsilon^2$$-robustly QC. □ 2. Properties of solutions In this section, we present technical arguments proving the uniqueness of solutions of our PDEs and discuss some relevant properties. 2.1 Comparison principle In this subsection, we will show that a weak comparison principle holds for the Dirichlet problem of $$h-F^\epsilon[u]=0$$ for $$h>0$$ 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)<g(x)\}$$. Then, $$h-F^\epsilon[v]\ge 0$$ in $$\varOmega_g$$ and $$v=g$$ on $$\partial\varOmega_g$$, that is, $$v$$ is a viscosity supersolution of $$h-F^\epsilon[v]=0$$ in $$\varOmega_g$$. Now, by the definition of viscosity subsolutions we have $$u\le g$$ and $$h-F^\epsilon[u]\le 0$$ in $$\varOmega$$ and thus also $$\varOmega_g$$. This allows us to conclude that $$u\le v$$ on $$\partial\varOmega_g$$. Therefore, by Proposition 2.4, $$u\le v$$ in $$\varOmega_g$$. Concluding, in $$\varOmega\setminus\varOmega_g$$ we necessarily have $$v\ge g\ge u$$. □ 2.2 The $$\epsilon$$-uniformly QCE We formulate the $$\epsilon$$-uniformly QCE of a function $$g$$ as the unique viscosity solution of the following obstacle problem: ($$\epsilon$$-Ob) \begin{equation} \begin{cases} \max\{u(x)-g(x),\epsilon-F^\epsilon[u](x)\} = 0, & x\in\varOmega,\\ u(x) = g(x), & x\in\partial\varOmega. \end{cases} \end{equation} Remark 2.6 (Convergence of approximate solutions). It is clear that as $$\epsilon\to0$$, the penalization term in $$F^\epsilon[u]$$ tends to infinity. The result is that $$F^\epsilon[u]\to\lambda_{\rm QC}[u]$$ as $$\epsilon\to0$$. From this observation, the standard stability result of viscosity solutions and the quasiconvexity of subsolutions of $$\epsilon-F^\epsilon[u]=0$$, one can then apply the same argument presented in the proof of Barron et al. (2013, Theorem 5.3) to conclude that the unique viscosity solutions of (ϵ-QCE), (Dir) converge to the QCE as $$\epsilon\to0$$. This result allows us to compute asymptotic approximations of the QCE of a given obstacle. 2.3 Solution formula in one dimension In one dimension, $$\epsilon-F^\epsilon[u]$$ is simply the Eikonal operator with a small diffusion term. When considering a solution $$u$$ of (ϵ-Ob), whenever $$u(x) < g(x)$$, we have $$F^\epsilon[u](x) = \epsilon$$, which gives \begin{equation} \epsilon u_{xx} + {\vert{u_x}\vert}=\epsilon^2. \end{equation} (2.3) Define $$\varOmega_g\equiv\{x:u(x) < g(x)\}$$. Note that $$\varOmega_g$$ need not be connected; it can be written as the union of finitely many intervals (refer to Fig. 3). However, we can solve the equation in each interval. Lemma 2.7 (One-dimensional solution). The viscosity solution of the one-dimensional PDE for the operator (2.1) and (2.3), along with boundary conditions $$u(0)=0, u(W) = H$$, is given by the following. Set $$S\equiv H/\epsilon^2W$$. Then, \[ u^\epsilon(x) = \begin{cases} \pm \epsilon^2 x, & S=\pm1,\\ \epsilon^2x + C^+(1-\exp(-x/\epsilon)), & S>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 $$x<x_0$$ and $$u_x>0$$ 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]<F_\rho[v]\implies u\le v\quad\text{in}\,\varOmega^h. \] Proof. Suppose for contradiction that \[ \max_{x\in{\it{\Omega}}^h}\{u(x) - v(x)\} = u(x_0) - v(x_0) > 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)<g(x)$$ for every $$x\in\overline\varOmega$$. Moreover, one can check that $$(g_1)_{vv}^h(x)=0$$. Then, we observe \begin{align*} G^\epsilon_{h,{\rm d}\theta}[g_1](x) &= \max\left\{g_1(x)-g(x),\epsilon-\min_{v\in V}\left\{ \frac{\max\{g_1(x\pm hv)\}-g_1(x)}{\epsilon h} \right\}\right\}\\ &= \max\left\{g_1(x)-g(x),\epsilon-2\epsilon\min_{v\in V}\left\vert\sum\nolimits_{i=1}^nv_i\right\vert\right\}\\ &= \max\left\{g_1(x)-g(x),-\epsilon\right\}, \end{align*} where the last equality follows from the general assumption that the stencil width is at least 1 with directional resolution $${\rm d}\theta\le\pi/4$$. Next, define $$g_2(x)\equiv\max_{x\in\overline{\it{\Omega}}}g(x)$$. Then, in summary we have $$G^\epsilon_{h,{\rm d}\theta}[g_1]<0$$ and $$G^\epsilon_{h,{\rm d}\theta}[g_2]\ge0$$. Therefore, by Proposition 5.10, \[ \min_{x\in\overline\varOmega} g_1(x)\le u(x)\le g_2 \quad\text{for}\ x\in\varOmega^h. \] Moreover, this bound is independent of $$h$$ and $${\rm d}\theta$$. □ Proposition 5.12. The solutions of (4.4) converge locally uniformly on $$\varOmega$$ to the solutions of (ϵ-Ob) as $$h,{\rm d}\theta\to0$$. Proof. By Lemmas 5.8, 5.9, 5.11 we see that (4.4) is a consistent and elliptic scheme which has stable solutions. Moreover, a weak comparison principle holds by Corollary 2.5. Therefore, by Proposition 5.6, the numerical solutions converge uniformly on compact subsets of $$\varOmega$$. □ 6. Numerical results In this section we present the numerical results. The tolerance for the fixed-point iterations was taken to be $$10^{-6}$$ and unless otherwise stated we set $$\epsilon=h/2$$, where $$h$$ is the spatial resolution of the grid. Technical conditions on $$g$$ given by Remark 2.1 are needed to ensure that the boundary conditions are held in the strong sense. In practice, we violated this condition, resulting in solutions that were discontinuous at the boundary. Two natural choices for initialization of the solution are (i) the obstacle function, $$g(x)$$ and (ii) a QC function below $$g$$. We found that the first choice leads to very slow convergence. In particular, the parabolic equation takes time $$\mathcal{O}(1/\epsilon)$$ to converge. See Example 6.3 below. On the other hand, the second choice results in faster convergence: the simplest choice is simply the constant function with value the minimum of $$g$$. After one step of the iteration the boundary values are attained. Moreover, starting with this choice allows us to use the iterative method with $$\epsilon = 0$$ to find the QCE. Remark 6.1 It is an open question whether solving the $$\epsilon=0$$ time-dependent PDE with quasiconvex initial data below $$g$$ will converge. To apply our convergence results, we take $$\epsilon = h/2$$. 6.1 Results in one dimension We present examples demonstrating the convergence of approximate QCEs to the true QCE as $$\epsilon\to0^+$$. We also compare the iteration count when starting from below and at the obstacle. Example 6.2 (Convergence as $$\epsilon\to0$$). We consider the convergence of numerical solutions of (4.4) as $$\epsilon\to 0$$ with the following obstacle function: \[ g(x)=\min\{{\vert{x-0.5}\vert},{\vert{x+0.5}\vert}-0.3\}. \] The results are displayed in Fig. 5. Indeed, as expected we see convergence to the true QCE, from below, as $$\epsilon\to0$$. Fig. 5. View largeDownload slide Obstacles (solid) and the numerical solutions (dashed). Left: Example 6.2. Convergence in $$\epsilon$$. Bottom to top: $$\epsilon = 0.2\,,0.1\,,0.05\,,0.001\,,0.0001$$. Right: Example 6.3. Fixed-point iterations when $$u^0 \equiv g_2$$. Fig. 5. View largeDownload slide Obstacles (solid) and the numerical solutions (dashed). Left: Example 6.2. Convergence in $$\epsilon$$. Bottom to top: $$\epsilon = 0.2\,,0.1\,,0.05\,,0.001\,,0.0001$$. Right: Example 6.3. Fixed-point iterations when $$u^0 \equiv g_2$$. Example 6.3 (Visualization of iterations). Next, we consider the obstacle \[ g(x)= -x^2+1, \] whose QCE is simply $$g_2^{\rm QCE}\equiv0$$. We demonstrate the evolution of the iterations when the initial data are taken to be the obstacle. In this case, the solution corresponded closely to $$u(x,t) = \min(g(x), c - \epsilon t)$$, where $$c = \max g$$. This illustrates the slow speed of convergence, and the fact that the equation degenerates to a trivial operator as $$\epsilon \to 0$$. Results are displayed in Fig. 5. 6.2 Numerical results in two dimensions All examples take the initial function in the iterative scheme to be the minimum value $$g_{\rm min} = \min g$$ of the obstacle function. In all contour plots, the solid line represents the level sets of the original function and the dashed line represents the same level sets of the numerical solution. Unless otherwise stated, the two-dimensional numerical solutions shown are computed on a $$64\times64$$ grid, for illustration purposes. Computations were performed on larger-sized grids. We performed the same numerical experiments in Abbasi & Oberman (2016), and achieved similar results. This is expected, since solving our PDE with small $$\epsilon$$ is consistent with the QCE. We do not reproduce those results here to save space. The examples we present focus on the difference between the two operators for values of $$e$$ close to 1. Example 6.4. (Uniform convexification of level sets). Let $$g(x)$$ be the signed distance function (negative on the inside) to the square, $$S = \left\{x \mid \max_{i}{\vert{x_i}\vert}=1/2 \right\}$$. Although $$g$$ is convex, its level sets are not uniformly convex. We apply the iterative procedure to $$g$$ with $$\epsilon=0.5$$ so that we uniformly convexify the level sets. The results demonstrating this are displayed in Fig. 6. Fig. 6. View largeDownload slide Left: Examples 6.4, square level sets become uniformly convex with $$\epsilon = 1/2$$. Middle: Example 6.5, the solution with $$\epsilon=1$$ is indicated by the dashed contours. Right: Example 6.6 with $$\epsilon$$ near 0. Fig. 6. View largeDownload slide Left: Examples 6.4, square level sets become uniformly convex with $$\epsilon = 1/2$$. Middle: Example 6.5, the solution with $$\epsilon=1$$ is indicated by the dashed contours. Right: Example 6.6 with $$\epsilon$$ near 0. Example 6.5 (Nonconvex signed distance function). In this example, $$g(x)$$ is the signed distance function to the curve $${\it{\Gamma}}$$, \begin{align*} {\it{\Gamma}}(t) &= \frac{1}{2}\begin{cases} (\cos(t-\frac{\pi}{2}),\sin(t-\frac{\pi}{2})), & t\in[0,3\pi/2],\\[9pt] \left(0,\frac{t-7\pi/4}{\pi/4}\right)\!, & t\in[3\pi/2,7\pi/4], \\[9pt] \left(\frac{7\pi/4-t}{\pi/4},0\right)\!, & t\in[7\pi/4,2\pi]. \end{cases} \end{align*} We compare the results of using different $$\epsilon$$. In particular, we choose $$\epsilon=h/2$$ and $$\epsilon=1$$. In the latter case, we see that the level sets are more curved. Results are displayed in Fig. 6. Example 6.6. We consider the obstacle where $$g$$ is a cone with circular portions removed from its level sets. Take \[ g(x) = \begin{cases} 1, &x\in A,\\ {\vert{x}\vert}-1/2, &x\in {\it{\Omega}}\setminus A, \end{cases} \] where $$A=\{(x_1\pm1/2)^2+x_2^2\le 1/16\}\bigcup \{x_1^2+(x_2\pm1/2)^2\le1/16\}$$. Results showing the $$g=0.7$$ level set are found in Fig. 6. Example 6.7 (Comparison to the $$\epsilon$$-robust QCE). We compare the solutions of our PDE to the solutions of the line solver presented in Abbasi & Oberman (2016), which returns the robust QCE. In particular, we use a stencil width $$W=5$$, we used the same directions in the line solver and we set $$\epsilon=0.02$$ for the $$\epsilon$$-uniformly QCE, and we also used a regularization of $$0.02$$ for the robust QC (in principle we should have used $$0.02^2$$ but the matching value was better for illustration purposes). The results are found in Fig. 7. Fig. 7. View largeDownload slide Example 6.7. Top: Surface plots (inverted) of the function, robust QCE, $$\epsilon$$-uniformly QCE. Bottom: Corresponding contour plots. Fig. 7. View largeDownload slide Example 6.7. Top: Surface plots (inverted) of the function, robust QCE, $$\epsilon$$-uniformly QCE. Bottom: Corresponding contour plots. 6.3 Accelerating iterations using the line solver We found an effective method to reduce the computational time to find the solution. We implement the line solver for the robust QCE proposed in Abbasi & Oberman (2016), alternating with the PDE iterations. In particular, on an $$n\times n$$ grid, after every $$2n$$ iterations of the PDE we apply one iteration of the line solver (with the commensurate value of $$\epsilon$$) in each direction. Table 1 shows the number of iterations required for convergence using the same obstacle from Abbasi & Oberman (2016, Example 6.1). Note the line solver is for a different regularization of the QC operator; however, for small values of $$\epsilon,$$ the solutions are close. For larger values of $$\epsilon,$$ the operator still accelerates the solver, but it does not approach the solution to within arbitrary precision. Table 1 Number of iterations required for convergence using different solution methods $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 Table 1 Number of iterations required for convergence using different solution methods $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 $$n$$ 32 64 128 256 2n PDE iterations 50 102 203 398 2n PDE It + 1 line solve 6 8 10 7 In Table 1, we compare the number of iterations required for each method. The results were comparable across different examples. The number of iterations of the PDE operator required for convergence was typically a small constant times $$n^2$$. The computational cost of a single line solve was on the order of a small constant (say 10) times the cost of a single PDE iteration. So the combined method, requires $$\mathcal{O}(n)$$ iterations (possibly with a $$\log n$$ prefactor, which we ignore, since it is not significant at these values of $$n$$). So the combined method required (roughly) $$n$$ iterations. Recall that there are $$N =n^2$$ grid points. Each iteration has a cost proportional to the number of directions used and to the number of grid points, so the total cost is $$\mathcal{O}\left((WN^2) \right)$$ for the iterative method and $$\mathcal{O}(WN)$$ for the combined method (with constants of roughly 1 and 10, respectively). So combining the iterative solver with the line solver results in significant improvements to computation time. References Abbasi B. & Oberman A. M. ( 2016 ) Computing the quasiconvex envelope using a nonlocal line solver. Available at https://arxiv.org/abs/1612.05584. Avriel M. , Diewert W. E. , Schaible S. , & Zang I. ( 2010 ) Generalized concavity . Society for Industrial and Applied Mathematics . Bardi M. & Mannucci P. ( 2006 ) On the Dirichlet problem for non-totally degenerate fully nonlinear elliptic equations. Commun. Pure Appl. Anal. , 5 , 709 – 731 . Google Scholar CrossRef Search ADS Bardi M. & Mannucci P. ( 2013 ) Comparison principles and Dirichlet problem for fully nonlinear degenerate equations of Monge–Ampère type. Forum Math , 25 , 1291 – 1330 . Google Scholar CrossRef Search ADS Barles G. & Souganidis P. E. ( 1991 ) Convergence of approximation schemes for fully nonlinear second order equations. Asymptotic Anal. , 4 , 271 – 283 . Barron E. N. , Goebel R. & Jensen R. R. ( 2012a ) Functions which are quasiconvex under linear perturbations. SIAM J. Optim. , 22 , 1089 – 1108 . Google Scholar CrossRef Search ADS Barron E. N. , Goebel R. & Jensen R. R. ( 2012b ) The quasiconvex envelope through first-order partial differential equations which characterize quasiconvexity of nonsmooth functions. Discrete Contin. Dyn. Syst. Ser. B , 17 , 1693 – 1706 . Google Scholar CrossRef Search ADS Barron E. , Goebel R. & Jensen R. ( 2013 ) Quasiconvex functions and nonlinear pdes. Trans. Am. Math. Soc. , 365 , 4229 – 4255 . Google Scholar CrossRef Search ADS Barron E. N. & Jensen R. R. ( 2014 ) A uniqueness result for the quasiconvex operator and first order PDEs for convex envelopes. In Annales de l’Institut Henri Poincare (C) Non Linear Analysis , vol. 31 . Elsevier Masson , pp. 203 – 215 . Google Scholar CrossRef Search ADS Boyd S. & Vandenberghe L. ( 2004 ) Convex Optimization . Cambridge : Cambridge University Press . Google Scholar CrossRef Search ADS Caffarelli L. A. & Spruck J. ( 1982 ) Convexity properties of solutions to some classical variational problems. Comm. Partial Differential Equations , 7 , 1337 – 1379 . Google Scholar CrossRef Search ADS Cannarsa P. & Sinestrari C. ( 2004 ) Semiconcave Functions, Hamilton-Jacobi Equations, and Optimal Control , vol. 58. Springer Science & Business Media . Carlier G. & Galichon A. ( 2012 ) Exponential convergence for a convexifying equation. ESAIM Control Optim. Calc. Var. , 18 , 611 – 620 . Google Scholar CrossRef Search ADS Colesanti A. & Salani P. ( 2003 ) Quasi-concave envelope of a function and convexity of level sets of solutions to elliptic equations. Math. Nachr. , 258 , 3 – 15 . Google Scholar CrossRef Search ADS Crandall M. G. , Ishii H. , & Lions P. L. ( 1992 ) User’s guide to viscosity solutions of second order partial differential equations. Bull. Amer. Math. Soc. , 27 , 1 – 67 . Google Scholar CrossRef Search ADS Froese B. ( 2016 ) Convergent approximation of surfaces of prescribed Gaussian curvature with weak Dirichlet conditions. Available at https://arxiv.org/abs/1601.06315. Kawohl B. ( 1985 ) Rearrangements and Convexity of Level Sets in PDE . Lecture Notes in Mathematics , vol. 1150 . Berlin, Heidelberg : Springer-Verlag , pp. 1 – 134 . Google Scholar CrossRef Search ADS Kohn R. & Serfaty S. ( 2007 ) Second-order PDE’s and deterministic games. In Proc. Internat. Congress Ind. Appl. Math. (ICIAM’07) , pp. 239 – 249 . Lions P.-L. ( 1998 ) Identification du cône dual des fonctions convexes et applications. C. R. Acad. Sci. Paris Sér. I Math. , 326 , 1385 – 1390 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2004 ) A convergent monotone difference scheme for motion of level sets by mean curvature. Numer. Math. , 99 , 365 – 379 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2006 ) Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton-Jacobi equations and free boundary problems. SIAM J. Numer. Anal. , 44 , 879 – 895 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2007 ) The convex envelope is the solution of a nonlinear obstacle problem. Proc. Amer. Math. Soc. , 135 , 1689 – 1694 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2008a ) Computing the convex envelope using a nonlinear partial differential equation. Math. Models Methods Appl. Sci. , 18 , 759 – 780 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2008b ) Wide stencil finite difference schemes for the elliptic Monge-Ampère equation and functions of the eigenvalues of the Hessian. Discrete Contin. Dyn. Syst. Ser. B , 10 , 221 – 238 . Google Scholar CrossRef Search ADS Sethian J. A. ( 1999 ) Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science , vol. 3 . Cambridge University Press . © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Nov 18, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off