Convergence of implicit schemes for Hamilton-Jacobi-Bellman quasi-variational inequalities

TL;DR

This paper provides rigorous convergence proofs for implicit numerical schemes solving Hamilton-Jacobi-Bellman quasi-variational inequalities, using a novel nonlocal consistency approach and Barles-Souganidis analysis.

Contribution

It introduces a new nonlocal consistency framework and proves convergence of implicit schemes for HJBQVIs, filling a gap in existing computational methods.

Findings

Established convergence of implicit schemes for HJBQVIs

Introduced the concept of nonlocal consistency

Proved results using comparison principles

Abstract

In [Azimzadeh, P., and P. A. Forsyth. "Weakly chained matrices, policy iteration, and impulse control." SIAM J. Num. Anal. 54.3 (2016): 1341-1364], we outlined the theory and implementation of computational methods for implicit schemes for Hamilton-Jacobi-Bellman quasi-variational inequalities (HJBQVIs). No convergence proofs were given therein. This work closes the gap by giving rigorous proofs of convergence. We do so by introducing the notion of nonlocal consistency and appealing to a Barles-Souganidis type analysis. Our results rely only on a well-known comparison principle and are independent of the specific form of the intervention operator.

Equations

$$\{\begin{array}{ll}\min \left(-\underset{b\in B}{\sup }\left\{\frac{\partial u}{\partial t}+L_{b}u+f(\cdot ,b)\right\},u-\mathcal{M}u\right)=0,& \text{on }[0,T)\times {\mathbb{R}}^d\\ u(T,\cdot )=\max \left(g,\mathcal{M}u(T,\cdot )\right),& \text{on }{\mathbb{R}}^d\end{array}$$

$$L_{b}u(t,x):=⟨\mu (x,b),D_{x}u(t,x)⟩+\frac{1}{2}\mathrm{trace}(\sigma {\sigma }^{⊺}(x,b)D_x^{2}u(t,x))$$

$$\mathcal{M}u(t,x):=\underset{z\in Z(t,x)}{\sup }\left\{u(t,x+\Gamma (t,x,z))+K(t,x,z)\right\}.$$

$$\{\begin{array}{ll}\min \left(-\underset{b\in B}{\sup }\left\{\frac{\partial u^k}{\partial t}+L_b{u}^k+f(\cdot ,b)\right\},u^k-\mathcal{M}{u}^{k-1}\right)=0,& \text{on }[0,T)\times {\mathbb{R}}^d\\ u^k(T,\cdot )=\max \left(g,\mathcal{M}{u}^{k-1}(T,\cdot )\right),& \text{on }{\mathbb{R}}^d.\end{array}$$

$$\mathcal{M}u(t,x)=\underset{z\in [0,x_2]}{\sup }\left\{u(t,\max \{x_1-z,0\},x_2-z)+\kappa z-c\right\}$$

$$\begin{gathered} F((t,x),r,(a,p),A,\ell ):= \\ \{\begin{array}{ll}\min (-{\sup }_{b\in B}\{a+⟨\mu (x,b),p⟩+\frac{1}{2}\mathrm{trace}(\sigma {\sigma }^{⊺}(x,b)A)+f(t,x,b)\},r-\ell ),& \text{if }t<T\\ \min (r-g(x),r-\ell ),& \text{if }t=T\end{array} \end{gathered}$$

$$F(y,u(y),Du(y),D_x^{2}u(y),\mathcal{M}u(y))=0,\text{for }y=(t,x)\in [0,T]\times {\mathbb{R}}^d.$$

$$F(y,u(y),D\phi (y),D_x^2\phi (y),\mathcal{M}u(y))=0.$$

$$F(y,\phi (y),D\phi (y),D_x^2\phi (y),\mathcal{M}\phi (y))=0.$$

$$F(x,u(x),Du(x),D^{2}u(x),\mathcal{M}u(x))=0\text{ for }x\in \overline{\Omega }$$

$$F(x,r,p,A,{\ell }_1)⩾F(x,r,p,B,{\ell }_2)\text{ for }A⪯B\text{ and }{\ell }_1⩽{\ell }_2$$

$$F_{\ast }(x,u^{\ast }(x),D\phi (x),D^2\phi (x),\mathcal{M}{u}^{\ast }(x))⩽0$$

$$F^{\ast }(x,u_{\ast }(x),D\phi (x),D^2\phi (x),\mathcal{M}{u}_{\ast }(x))⩾0\text{)}.$$

$$F_{\ast }(x,\phi (x),D\phi (x),D^2\phi (x),\mathcal{M}\phi (x))⩽0.$$

$$0=(u^{\ast }-\phi )(x)⩾(u^{\ast }-\phi )(y).$$

$$F_{\ast }(x,\phi (x),D\phi (x),D^2\phi (x),\mathcal{M}\phi (x))⩽F_{\ast }(x,u^{\ast }(x),D\phi (x),D^2\phi (x),\mathcal{M}{u}^{\ast }(x))⩽0.$$

$$\mathcal{M}(v+ϵ)⩽\mathcal{M}v+\omega (ϵ)\text{(resp. }\mathcal{M}(v-ϵ)⩾\mathcal{M}v-\omega (ϵ)\text{)}$$

$$\mathcal{M}u⩽\mathcal{M}\psi ⩽\mathcal{M}(u+ϵ)⩽\mathcal{M}u+\omega (ϵ)=:{\mathcal{M}}_{ϵ}u$$

$$F_{\ast }(x,u(x),D\phi (x),D^2\phi (x),{\mathcal{M}}_{ϵ}u(x))⩽F_{\ast }(x,\psi (x),D\psi (x),D^2\psi (x),\mathcal{M}\psi (x))⩽0.$$

$$\chi (\cdot )=\sum _n{c}_n{\chi }_n(\cdot )\text{where }{c}_n:=\sup u(U_n).$$

$$\chi (y)=\sum _n{c}_n{\chi }_n(y)⩾\left(\underset{n\in \mathcal{N}}{\min }{c}_n\right)\sum _n{\chi }_n(y)=\underset{n\in \mathcal{N}}{\min }{c}_n⩾u(y).$$

$$u(w)-u(y)=\left|u(w)-u(y)\right|⩽\left|u(w)-u(x_m)\right|+\left|u(x_m)-u(y)\right|⩽ϵ.$$

$$\psi (y):=\zeta (y)\phi (y)+\left(1-\zeta (y)\right)\chi (y)$$

$$\text{if }u,w\in U\text{ are a (def.}k\text{)-subsolution}\text{ and (def.}k\text{)-supersolution of }\text{(refeq:pde)}\text{, then }{u}^{\ast }⩽w_{\ast }.$$

$$\overline{u}(x):=\underset{\begin{array}{c}\Delta x\downarrow 0\\ y\to x\end{array}}{\mathrm{limsup}}{u}^{\Delta x}(y)\text{and}\underline{u}(x):=\underset{\begin{array}{c}\Delta x\downarrow 0\\ y\to x\end{array}}{\mathrm{liminf}}{u}^{\Delta x}(y)$$

$$S(\rho ,x,u^{\rho }(x),[u^{\rho }{]}_x)=0\text{ for }x\in \overline{\Omega }.$$

$$[w{]}_x(y):=\{\begin{array}{ll}w(y)& \text{if }y\ne x\\ 0& \text{otherwise}.\end{array}$$

$$F(x,u(x),Du(x),D^{2}u(x))=0,$$

$$\underset{\begin{array}{c}\rho \downarrow 0\\ y\to x\\ \xi \to 0\end{array}}{\mathrm{limsup}}S(\rho ,y,\phi (y)+\xi ,[\phi +\xi {]}_y)$$

$$\underset{\begin{array}{c}\rho \downarrow 0\\ y\to x\\ \xi \to 0\end{array}}{\mathrm{liminf}}S(\rho ,y,\phi (y)+\xi ,[\phi +\xi {]}_y)$$

$$S(\rho ,x,u^{\rho }(x),[u^{\rho }{]}_x,{\mathcal{N}}^{\rho }{u}^{\rho }(x))=0\text{ for }x\in \overline{\Omega }$$

Full text

Key words. implicit numerical scheme, Hamilton-Jacobi-Bellman quasi-variational inequality (HJBQVI), viscosity solution, impulse control

AMS subject classifications. 49L25, 49N25, 65M06, 65M12, 93E20

1   Introduction

We consider numerical methods for approximating the viscosity solution of the Hamilton-Jacobi-Bellman quasi-variational inequality (HJBQVI)

[TABLE]

where $L_{b}$ is the second order operator

[TABLE]

and $\mathcal{M}$ is the so-called “intervention” operator

[TABLE]

In (2), $D_{x}u$ and $D_{x}^{2}u$ are the gradient vector and second derivative matrix of $u$ with respect to the spatial variable $x$. In (3), $\cup_{t,x}Z(t,x)$ is a subset of some metric space $Y$ and $\Gamma$ and $K$ are maps from $[0,T]\times\mathbb{R}^{d}\times Y$ to $\mathbb{R}^{d}$ and $\mathbb{R}$, respectively.

It is well-known (see, e.g., [29, Chapter 7]) that the solution of the HJBQVI can be obtained as the limit of the sequence $(u^{k})_{k}$ where $u^{0}=-\infty$ and, for $k\geqslant 1$,

[TABLE]

The above is a variational (not quasi-variational) inequality since the obstacle $\mathcal{M}u^{k-1}$ does not depend on the solution $u^{k}$. This approach, referred to as iterated optimal stopping, suggests a simple algorithm for solving the HJBQVI numerically: fix mesh sizes $\Delta t$ and $\Delta x$ and compute a finite difference approximation to $u^{1}$, use that to compute an approximation to $u^{2}$, etc. until convergence up to a desired tolerance. Unfortunately, this approach suffers from a high space complexity [5].

In light of this, a previous work of the first author [4] considers various schemes which do not resort to a variational approximation as per the previous paragraph. In [4], the focus is on the theory and implementation of computational methods for solving these schemes; no convergence proofs are given. In this work, we close the gap by giving rigorous convergence results. The schemes in [4] are “implicit” in the sense that second derivative terms are approximated using information from the current timestep, ensuring unconditional stability. This is in contrast to so-called “explicit” schemes, which approximate second derivative terms using information from the previously computed timestep, requiring a harsh timestep-size restriction of the form $\Delta t=\operatorname{const.}(\Delta x)^{2}$ to ensure stability (see, e.g., [27, Chapter 4]).

Our closest related works [18, 14] introduce implicit (in the sense above) finite difference schemes to solve HJBQVIs arising from insurance applications. However,

• •

[18, 14] assume a strong form of comparison principle for the HJBQVI which, to the best of our knowledge, does not exist in the literature. In contrast, we do not make any such assumptions, relying only on a well-known comparison principle (see Appendix A).

• •

Moreover, [18, 14] assume a specific form for the intervention operator $\mathcal{M}$ and as such, the convergence proofs do not generalize. For example, [18] considers a two-dimensional problem ($x=(x_{1},x_{2})$) with

[TABLE]

where $\kappa$ and $c$ are constants. In contrast, we do not assume specific forms for the functions $Z$, $\Gamma$, and $K$ in (3) (see, in particular, Lemma 12).

A significant part of focus is on the first issue, which is subtle in nature. To briefly sketch the issue, we first rewrite the HJBQVI in a more abstract form. Namely, let

[TABLE]

where we have used the notation $(a,p)$ to distinguish between the time and spatial derivatives (i.e., $a=\partial u/\partial t(t,x)$ and $p=D_{x}u(t,x)$) and $A$ to refer to the second spatial derivative (i.e., $A=D_{x}^{2}u(t,x)$). Letting $D\coloneqq(\partial u/\partial t,D_{x})$, (1) becomes

[TABLE]

Now, there are a few reasonable notions of viscosity solution for (5). Roughly speaking, one notion is obtained by replacing derivatives of $u$ with those of a test function $\varphi$:

[TABLE]

Another is obtained by replacing instead all terms involving $u$:

[TABLE]

Since the operator $\mathcal{M}$ is nondecreasing (i.e., $\mathcal{M}u\leqslant\mathcal{M}w$ whenever $u\leqslant w$), an upper semicontinuous (USC) subsolution in the sense of (def.1) is also a USC subsolution in the sense of (def.2). The same is true for lower semicontinuous (LSC) supersolutions and as such, a comparison principle in the sense of (def.2) implies a comparison principle in the sense of (def.1) (see Section 2 for details).

As hinted at above, the results of [18, 14] require the stronger (def.2)-comparison principle which, to our knowledge, does not exist in the literature for the HJBQVI. To be able to rely only on a well-known (def.1)-comparison principle, we introduce the notion of nonlocal consistency. Analogously to the Barles-Souganidis framework (BSF) for local equations [11], we establish that a monotone, stable, and nonlocally consistent scheme whose limiting equation satisfies a (def.1)-comparison principle is convergent. We use this to prove the convergence of the implicit schemes from [4].

This work is organized as follows. Section 2 discusses the different notions of viscosity solutions and their relationship with the BSF. Section 3 introduces the notion of nonlocal consistency and gives the BSF-type convergence result described in the previous paragraph. Section 4 applies this result to prove the convergence of various schemes for the HJBQVI.

2   Viscosity solutions of nonlocal equations

In this and the following section, instead of considering the HJBQVI directly, we consider the more general PDE

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ is open, $\overline{\Omega}$ is its closure, $u$ is in $B_{\operatorname{loc}}(\overline{\Omega})$ (the set of locally bounded real-valued maps from $\overline{\Omega}$), and $Du$ and $D^{2}u$ are the formal gradient vector and second derivative matrix of $u$. Note that no generality is lost in using $x$ instead of $(t,x)$ in (6) since a “time” variable $t$ can be incorporated by increasing the dimensionality $d$ of the space. With this in mind, it is clear that the HJBQVI is a special case of (6) which is obtained when $F$ and $\mathcal{M}$ are defined as per (3) and (4).

In our analyses, we will always assume that $\mathcal{M}$ maps $B_{\operatorname{loc}}(\overline{\Omega})$ to some subset of itself and that $F$ is a locally bounded real-valued map from $\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d}\times\mathcal{S}^{d}\times\mathbb{R}$ that is elliptic in the sense of [9, Assumption (E)]:

[TABLE]

where $\preceq$ is the usual semidefinite order on $\mathcal{S}^{d}$ (the real $d\times d$ symmetric matrices). To stress the significance of the operator $\mathcal{M}$, we call (6) a nonlocal second order equation.

For a locally bounded real-valued function $z$ mapping from a metric space, we define its upper and lower semicontinuous envelopes by $z^{*}(x)\coloneqq\limsup_{y\rightarrow x}z(y)$ and $z_{*}\coloneqq-(-z)^{*}$ where $\limsup_{y\rightarrow z}z(y)\coloneqq\lim_{\epsilon\downarrow 0}(\sup\{z(y)\colon|y-x|<\epsilon\})$. We are now ready to state the definition of viscosity solution hinted at in (def.1).

Definition 1**.**

$u\in B_{\operatorname{loc}}(\overline{\Omega})$* is a (viscosity) subsolution (resp. supersolution) of (6) if for all $\varphi\in C^{2}(\overline{\Omega})$ and $x\in\overline{\Omega}$ such that $u^{*}-\varphi$ (resp. $u_{*}-\varphi$) has a local maximum (resp. minimum) at $x$, we have*

[TABLE]

$u$* is said to be a (viscosity) solution of (6) if it is both a sub and supersolution of (6).*

Our usage of $F_{*}$ and $F^{*}$ above is so that we may write both the partial differential equation and its boundary conditions as a single expression (for a detailed explanation, see [11, Pg. 274]). We now state the definition of viscosity solution hinted at in (def.2).

Definition 2**.**

$u\in B_{\operatorname{loc}}(\overline{\Omega})$* is a subsolution of (6) if for all $\varphi\in C^{2}(\overline{\Omega})$ and $x\in\overline{\Omega}$ such that $(u^{*}-\varphi)(x)=0$ is a global maximum of $u^{*}-\varphi$, we have*

[TABLE]

Supersolutions are defined symmetrically.

Throughout this article, it should be assumed that unless otherwise specified, the terms subsolution/supersolution/solution refer to the concepts in Definition 1. When we wish to be explicit, we will write (def.1)-subsolution/supersolution/solution or (def.2)-subsolution/supersolution/solution.

If the operator $\mathcal{M}$ is nondecreasing (i.e., $\mathcal{M}u\leqslant\mathcal{M}w$ whenever $u\leqslant w$), then any (def.1)-subsolution (resp. supersolution) is trivially a (def.2)-subsolution (resp. supersolution). To see this, let $u$ be a (def.1)-subsolution, $\varphi\in C^{2}(\overline{\Omega})$, and $x\in\overline{\Omega}$ be such that $(u^{*}-\varphi)(x)=0$ is a global maximum of $u^{*}-\varphi$. In this case, for any $y\in\overline{\Omega}$,

[TABLE]

Therefore, $\varphi\geqslant u^{*}$, from which we obtain $\mathcal{M}\varphi\geqslant\mathcal{M}u^{*}$. Using the ellipticity of $F$,

[TABLE]

As for the converse, a (def.2)-subsolution (resp. supersolution) need not be a (def.1)-subsolution (resp. supersolution), even if $\mathcal{M}$ is nondecreasing. This is because the term $\mathcal{M}\varphi(x)$ is not determined by values of $\varphi$ in a neighbourhood of $x$. In order to control this term, we need to impose some additional continuity:

Proposition 3**.**

Suppose the operator $\mathcal{M}$ is nondecreasing and that there exists a function $\omega:[0,\infty]\rightarrow[0,\infty]$ satisfying $\lim_{\epsilon\downarrow 0}\omega(\epsilon)=0$ and

[TABLE]

for all $v\in C(\overline{\Omega})$ and $\epsilon>0$. If $u$ is a continuous (def.2)-subsolution (resp. supersolution), it is also a (def.1)-subsolution (resp. supersolution).

The above generalizes [24, Theorem 3.1], [8, Proposition 1.2], and [20, Section 4]. We point out that the requirement involving the modulus of continuity $\omega$ is rather innocuous. For example, when $\mathcal{M}$ is the intervention operator from the introduction, the choice of $\omega(\epsilon)=\epsilon$ satisfies the requirement under mild conditions (e.g., $Z(t,x)$ is nonempty and compact for each $(t,x)$ and $\Gamma$ and $K$ are continuous; see also assumptions (H1) and (H2) of the sequel).

Proof.

Let $u$ be a continuous (def.2)-subsolution (we omit the supersolution case, which is handled similarly). Let $\varphi\in C^{2}(\overline{\Omega})$ and $x\in\overline{\Omega}$ be such that $u-\varphi$ has a local maximum at $x$. Without loss of generality, we may assume $u(x)=\varphi(x)$. Let $\epsilon>0$. We claim that, by the continuity of $u$, we can construct a smooth function $\psi$ such that $u\leqslant\psi\leqslant u+\epsilon$ everywhere and $\psi=\varphi$ on a neighbourhood of $x$. In this case,

[TABLE]

and hence

[TABLE]

Taking limit inferiors of this inequality with respect to $\epsilon$ establishes the desired result.

We now return to the claim above. To simplify notation, we can, without loss of generality, assume $u$ and $\varphi$ are maps from $\mathbb{R}^{d}$. By Lindelöf’s lemma, we can find a countable cover $\{U_{n}\}_{n}$ of $\mathbb{R}^{d}$ by open balls $U_{n}$ centred at points $x_{n}$ and satisfying $|u(x_{n})-u(\cdot)|<\epsilon/2$ on $U_{n}$. By virtue of this, we can also find a smooth partition of unity $\{\chi_{n}\}_{n}\subset C^{\infty}(\mathbb{R}^{d})$ subordinate to $\{U_{n}\}_{n}$. Let

[TABLE]

Now, fix $y$ and let $\mathcal{N}\coloneqq\{n\colon\chi_{n}(y)\neq 0\}$. Then, since $\mathcal{N}$ is a finite set,

[TABLE]

Similarly, $\chi(y)\leqslant\max_{n\in\mathcal{N}}c_{n}=c_{m}$ for some $m\in\mathcal{N}$. By definition, $c_{m}=u(w)$ for some $w$ and hence

[TABLE]

Since $y$ was arbitrary, this establishes $u\leqslant\chi\leqslant u+\epsilon$.

Now, pick $r>0$ small enough so that $u\leqslant\varphi\leqslant u+\epsilon$ on $\overline{B(x,2r)}$, the closed ball centred at $x$ with radius $2r$. Let $\zeta$ be a smooth cutoff function satisfying $0\leqslant\zeta\leqslant 1$, $\zeta=1$ on $\overline{B(x,r)}$, and $\zeta=0$ outside of $B(x,2r)$. Then, the function $\psi$ defined by

[TABLE]

satisfies our requirements. ∎

We close this section by discussing the BSF. In order to do so, we first recall the notion of a comparison principle. Let $U$ be some subset of $B_{\operatorname{loc}}(\overline{\Omega})$. For $k=1,2$, we say (6) satisfies a (def.$k$)-comparison principle in $U$ whenever the following holds:

[TABLE]

The discussion preceding this paragraph implies that if $\mathcal{M}$ is nondecreasing, a (def.2)-comparison principle implies a (def.1)-comparison principle (but not vice versa).

The BSF provides a general approach for proving the convergence of finite difference schemes to the viscosity solution $u$ of a PDE [11]. The idea is as follows: letting $u^{\Delta x}:\overline{\Omega}\rightarrow\mathbb{R}$ denote the solution of a finite difference scheme (for a fixed grid size $\Delta x>0$), the BSF outlines three sufficient conditions (monotonicity, stability, and consistency) which guarantee the functions $\overline{u}$ and $\underline{u}$ defined by

[TABLE]

to be a subsolution/supersolution pair of the limiting equation. Applying a comparison principle, one concludes $\overline{u}=\underline{u}=u$. In our context, the subtlety here is that if $\overline{u}$ and $\underline{u}$ are a (def.2)-subsolution/supersolution pair, the above argument requires a (def.2)-comparison principle.

More generally, the BSF characterizes when a family of approximations $(u^{\rho})_{\rho>0}$ converges to $u$ as $\rho\downarrow 0$ (in the example above, $\rho=\Delta x$ was the grid size). The BSF represents an approximation $u^{\rho}$ as a solution of the equation

[TABLE]

The function $S$ is referred to as the approximation scheme. Given a function $w$, the symbol $[w]_{x}$ (cf. [10]) denotes a function which agrees with $w$ everywhere except possibly at $x$:

[TABLE]

While the BSF is posed in terms of local equations of the form

[TABLE]

the authors of [18, 14] (see, in particular, [18, Lemma 5.2] and [14, Definition 6.11]) noticed that the BSF could also be applied to nonlocal equations of the form (6) by changing the consistency condition of the BSF [11, Eq (2.4)] to

[TABLE]

The issue with this approach, however, is that it requires a (def.2)-comparison principle. In [18, 14], the authors avoid the issue by simply assuming that such a comparison principle holds (see, in particular, [18, Assumption 5.1] and [14, Assumption 6.3]), though it is not clear if such an assumption is justifiable for HJBQVIs.

We show in the next section that an alternative notion of consistency can be used so that only a (def.1)-comparison principle is required. Ultimately, we will use this to prove the convergence of various finite difference schemes for HJBQVIs, relying only on a well-known (def.1)-comparison principle given in Appendix A.

3   Convergence result

In our setting, approximation schemes take the form

[TABLE]

(compare with (7)). In particular, $S$ is a real-valued map from $(0,\infty)\times\overline{\Omega}\times\mathbb{R}\times B(\overline{\Omega})\times\mathbb{R}$, $\mathcal{N}^{\rho}$ maps $B(\overline{\Omega})$ to some subset of itself, and $B(\overline{\Omega})$ is the set of bounded real-valued maps from $\overline{\Omega}$. Intuitively, $\mathcal{N}^{\rho}$ will serve as an approximation of the operator $\mathcal{M}$. We refer to a function $u^{\rho}$ in $B(\overline{\Omega})$ as a solution of the scheme if it satisfies (9).

For brevity, let $\mathbb{R}_{+}\coloneqq(0,\infty)$. We will always assume the existence of functions $h_{1}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ and $h_{2}:\mathbb{R}_{+}\times\overline{\Omega}\rightarrow\mathbb{R}$ satisfying

[TABLE]

along with

[TABLE]

This assumption is a technical one required to allow solutions $u^{\rho}$ of the scheme to be discontinuous (see Remark 5 for an explanation). It is readily verified that all schemes in this work satisfy this assumption.

We will show that if a scheme $S$ is monotone, stable, and consistent, it converges to the unique solution of (6), provided that (6) satisfies a (def.1)-comparison principle. While our notions of monotonicity and stability are identical to those in the BSF, our notion of consistency will differ from (8). We now state precisely these notions.

A scheme $S$ is stable if (cf. [11, Eq. (2.3)])

[TABLE]

A scheme $S$ is monotone if (cf. [11, Eq. (2.2)])

[TABLE]

Note that the monotonicity requirement does not involve the operator $\mathcal{N}^{\rho}$ and as such, the reader may be tempted to guess that high order discretizations of the nonlocal operator $\mathcal{M}$ are possible. Unfortunately, this is not the case: high order discretizations are generally precluded by our notion of consistency (see Remark 13).

Before we introduce our notion of consistency, we recall half-relaxed limits. For a family $(z^{\rho})_{\rho>0}$ of real-valued maps from a metric space such that $(\rho,x)\mapsto z^{\rho}(x)$ is locally bounded, we define the upper and lower half-relaxed limits $\overline{z}$ and $\underline{z}$ by

[TABLE]

where

[TABLE]

and the limit inferior is defined similarly (cf. [7, Definition 1.4]). Now, a scheme $S$ is nonlocally consistent if for each family $(w^{\rho})_{\rho>0}$ of uniformly bounded real-valued maps from $\overline{\Omega}$, $\varphi\in C^{2}(\overline{\Omega})$, and $x\in\overline{\Omega}$, we have

[TABLE]

While (8) and (12) are aesthetically similar, the latter does not apply test functions to the operator $\mathcal{M}$. This is perhaps not too surprising if we think about (8) and (12) as notions of consistency relating to (def.2) and (def.1), respectively: in Definition 2, test functions are applied to the operator $\mathcal{M}$ while in Definition 1, test functions are not applied to $\mathcal{M}$. We are now ready to state the BSF-type convergence result.

Theorem 4**.**

Let $S$ be a monotone, stable, and nonlocally consistent scheme whose limiting equation (6) satisfies a (def.1)-comparison principle (in $B(\overline{\Omega})$). Then, as $\rho\downarrow 0$, the solution $u^{\rho}$ of (9) converges locally uniformly to the unique (def.1)-solution of (6) in $B(\overline{\Omega})$.

The proof is very similar to that of [11, Theorem 2.1]. Regardless, we give a detailed proof so that the reader can see the motivation behind nonlocal consistency.

Proof.

Let $\overline{u}\geqslant\underline{u}$ denote the half-relaxed limits of the family $(u^{\rho})_{\rho>0}$. We claim that $\overline{u}$ is a subsolution and $\underline{u}$ is a supersolution of (6). In this case, the comparison principle yields $\overline{u}=\overline{u}^{*}\leqslant\underline{u}_{*}=\underline{u}$, from which the desired result follows.

We now prove that $\overline{u}$ is a subsolution (that $\underline{u}$ is a supersolution is proved similarly). To this end, let $x\in\overline{\Omega}$ be a local maximum point of $\overline{u}-\varphi$ for some $\varphi\in C^{2}(\overline{\Omega})$. By definition, we can find a neighbourhood (relative to $\overline{\Omega}$) $U$ whose closure is compact and on which $x$ is a global maximum point of $\overline{u}-\varphi$. Without loss of generality, we may assume the maximum is strict, $\overline{u}(x)=\varphi(x)$, and $\varphi\geqslant\sup_{\rho}\|u^{\rho}\|_{\infty}$ outside $U$. By the definition of $\overline{u}$, we can find a sequence $(\rho_{n},x_{n})_{n}$ such that $\rho_{n}\downarrow 0$, $x_{n}\rightarrow x$, and $u^{\rho_{n}}(x_{n})\rightarrow\overline{u}(x)$. For each $n$, pick $y_{n}\in U$ such that

[TABLE]

where $h_{1}$ is the function in (10). Now, pick a subsequence of $(\rho_{n},x_{n},y_{n})_{n}$ such that its last argument converges to some point $\hat{y}\in\overline{U}$. With a slight abuse of notation, relabel this subsequence $(\rho_{n},x_{n},y_{n})_{n}$. It follows that

[TABLE]

Because $x$ was assumed to be a strict maximum point, the above implies $\hat{y}=x$. Letting $\xi_{n}\coloneqq(u^{\rho_{n}}-\varphi)(y_{n})+h_{1}(\rho_{n})$, we have $\xi_{n}\rightarrow 0$ and $u^{\rho_{n}}\leqslant\varphi+\xi_{n}$.

The definition of $u^{\rho}$ and the monotonicity of $S$ yield

[TABLE]

Taking limit inferiors, we get

[TABLE]

Now, employing (10), (11), and nonlocal consistency,

[TABLE]

which is the desired inequality, since $\overline{u}(x)=\varphi(x)$. ∎

Remark 5**.**

(10)* and (11) were needed in the proof above since the supremum in (13) is not necessarily attained at a point in $\overline{U}$. If we restrict our attention to the case in which solutions of the scheme $u^{\rho}$ are continuous, then this problem disappears. However, our analysis requires discontinuous solutions (see, e.g., (21) in the sequel). This issue is discussed further in the first author’s note [2].*

Remark 6**.**

Recall that the functions $(w^{\rho})_{\rho>0}$ appearing in the nonlocal consistency inequalities (12) are an arbitrary family of uniformly bounded real-valued maps. However, a close inspection of the proof of Theorem 4 reveals that we only use these inequalities with $w^{\rho}=u^{\rho}$ where $u^{\rho}$ is a solution of the scheme. Therefore, in establishing nonlocal consistency we can, without loss of generality, assume that $w^{\rho}$ is a solution of the scheme.

As in [11], the result above extends to solutions that are not necessarily bounded. In particular, let $r:\mathbb{R}\rightarrow\mathbb{R}$ be a map defined by $r(x)\coloneqq\operatorname{const.}\,(1+x^{d})$ where $d$ is a positive integer and let $R(\overline{\Omega})$ be the set of maps $u:\overline{\Omega}\rightarrow\mathbb{R}$ satisfying

[TABLE]

Now, relax the stability condition to read

[TABLE]

and the consistency condition by requiring (12) to hold more generally for families $(w^{\rho})_{\rho>0}\subset R(\overline{\Omega})$ not necessarily uniformly bounded. Then, by replacing instances of $B(\overline{\Omega})$ by $R(\overline{\Omega})$ in Theorem 4, we obtain a straightforward relaxation of Theorem 4 which allows for solutions of polynomial growth. Theorem 4 can also be extended to the case in which the PDE involves a family of nonlocal operators (cf. [9, Eq. (2)]):

[TABLE]

4   Hamilton-Jacobi-Bellman quasi-variational inequality

Using the results of the previous section, we now prove the convergence of two implicit schemes for HJBQVIs. Both schemes appear in the first author’s work [4]. We make the following assumptions about the various quantities appearing in the HJBQVI (1):

1. (H1)

$f$, $g$ $\mu$, $\sigma$, $\Gamma$, and $K$ are continuous functions with $f$ and $g$ bounded and $\mathcal{M}g\leqslant g$ where $\mathcal{M}g(x)\coloneqq\sup_{z\in Z(T,x)}\{g(x+\Gamma(T,x,z))+K(T,x,z)\}$. 2. (H2)

$B$ is a nonempty compact metric space and $Z(t,x)\subset Y$ is a nonempty compact set for each $(t,x)$. 3. (H3)

$\sup_{t,x,z}K(t,x,z)<0$.

It is well-known (see, e.g., [29]) that the HJBQVI is the dynamic programming equation associated to the optimal control problem whose value function is given by

[TABLE]

where

[TABLE]

Here, $W$ is a $d$-dimensional Brownian motion with filtration $(\mathcal{F}_{t})_{t\geqslant 0}$, $\gamma=(a,b)$ is a combined stochastic and impulse control111The interested reader can review [15, 25, 33, 29, 32, 20, 16, 19, 12, 17, 23, 6, 13, 26, 1] for information and recent developments on impulse control. We do not claim this list to be exhaustive. with $a=(\tau_{i},z_{i})_{i\geqslant 1}$, $(\tau_{i})_{i\geqslant 1}$ being an increasing sequence of stopping times, each $z_{i}$ being $\mathcal{F}_{\tau_{i}}$-measurable and taking values in some compact $Z(\tau_{i},X_{\tau_{i}-}^{t,x,\gamma})$, and $b$ being a progressively measurable process taking values in some compact $B$.

In the context of the optimal control problem (14), the assumption $\mathcal{M}g\leqslant g$ implies that it is suboptimal to perform an impulse at the terminal time $T$. The assumption (H3), needed to ensure stability, can be interpreted as the controller paying a cost for the right to perform an impulse. Lastly, we mention that the choice to write the functions $\mu$ and $\sigma$ as independent of the time coordinate $t$ was made only to simplify notation, and not due to a shortcoming of the theory.

4.1   Penalty scheme

In this subsection, we prove the convergence of a penalty scheme from [4, Section 5.2] to the solution of the HJBQVI. We start with the one-dimensional case ($d=1$), deferring for now a discussion of the complications that arise in higher dimensions. In this case, the operator $L_{b}$ simplifies to

[TABLE]

To simplify presentation, we assume uniformly spaced grid points $\{n\Delta t\}_{n=0}^{N}$ and $\{j\Delta x\}_{j=-M}^{M}$ in time and space with $T=N\Delta t$ and $M>0$. For brevity, we let $x_{j}\coloneqq j\Delta x$. We assume $\Delta t$ and $\Delta x$ are related through a positive parameter $\rho$ as per

[TABLE]

For PDEs on unbounded domains, the usual approach is to truncate the domain to a bounded region (e.g., $[-Q,Q]$) and impose artificial boundary conditions. Since these arguments are somewhat standard, we will assume

[TABLE]

in order to ignore artificial boundary conditions in our convergence analysis and focus instead on the main difficulties specific to the HJBQVI. For example, to satisfy (15), we may take $M=\operatorname{const.}\rho^{-1-\alpha}$ for some $\alpha>0$ so that $Q=\operatorname{const.}\rho^{-\alpha}$. Ignoring artificial boundary conditions to focus on the main points is a common practice in the literature on elliptic equations (see, e.g., [21, Section 3], in which the authors take an infinite grid to numerically solve a Hamilton-Jacobi-Bellman-Isaacs (HJBI) equation).

Writing $u_{j}^{n}$ for the quantity $u(n\Delta t,x_{j})$, let $u^{n}\coloneqq(u_{j}^{n})_{j=-M}^{M}$ be the image of $u(n\Delta t,\cdot)$ on the spatial grid. Note that $u^{n}$ is a vector in $\mathbb{R}^{2M+1}$. In order to present the penalty scheme, our first task is to define approximations $\mathcal{D}$ and $\mathcal{D}^{2}$ of the first and second derivatives $u_{x}$ and $u_{xx}$ in the sense that

[TABLE]

For the second derivative, we use a standard three point stencil:

[TABLE]

For the first derivative, we use an upwind stencil:

[TABLE]

These stencils should be modified accordingly in the case of a non-uniform spatial grid. Since the expressions (16) and (17) are not well-defined at $j=\pm M$, we set $\mathcal{D}$ and $\mathcal{D}^{2}$ to zero there for convenience:

[TABLE]

Noting that the upwind direction in (17) depends on the choice of $b$, we will sometimes write $\mathcal{D}_{b}$ in lieu of $\mathcal{D}$ to make this dependence explicit.

Remark 7**.**

Assumption (15) is computationally expensive to implement in practice and as such, a practical implementation should instead take $Q$ to be constant independent of $\rho$. In this case, we can interpret

[TABLE]

as imposing the Neumann boundary condition $u_{x}(t,\pm Q)=u_{xx}(t,\pm Q)=0$ (other boundary conditions can be imposed by modifying the scheme appropriately). In the context of the optimal control problem (14), this corresponds to (formally) modifying the drift $\mu$ and diffusion $\sigma$ of $X^{t,x,\gamma}$ by setting them to zero outside of $(-Q,Q)$.

Since the control set $B$ appearing in (1) may be infinite, we replace it with a nonempty finite subset $B^{\rho}\subset B$. Similarly, since the set $Z(t,x)$ appearing in the intervention operator (3) may be infinite at each point $(t,x)$, we replace it with a nonempty finite subset $Z^{\rho}(t,x)\subset Z(t,x)$. To ensure consistency, we require:

1. (H4)

As $\rho\downarrow 0$, $B^{\rho}$ converges to $B$ (in the Hausdorff metric) and $Z^{\rho}$ converges locally uniformly to $Z$ (in the Hausdorff metric). Moreover, $(t,x)\mapsto Z^{\rho}(t,x)$ is continuous (in the Hausdorff metric) for each $\rho$.

Unless otherwise mentioned, we will always assume the Hausdorff metric when discussing compact sets. Note that the above implies that $Z$ is also continuous.

Next, note that the point $x+\Gamma(t,x,z)$ appearing in the intervention operator (3) is not necessarily a point on the grid. Therefore, a discretization of the intervention operator requires interpolation. We use $\operatorname{interp}(u^{n},x)$ to denote the value of $u(n\Delta t,x)$ as approximated by a standard monotone linear interpolant. Precisely, for a point $x$ within the grid (i.e., $x_{-M}<x<x_{M}$), we define

[TABLE]

where $k$ is the unique integer satisfying $x_{k}\leqslant x<x_{k+1}$ and $\alpha=(x-x_{k})/(x_{k+1}-x)$. If $x$ is not within the grid (i.e., $x\leqslant x_{-M}$ or $x\geqslant x_{M}$), we define

[TABLE]

We can now discretize the intervention operator according to

[TABLE]

Using the notation above, $\mathcal{M}u(n\Delta t,x_{j})\approx(\mathcal{M}^{\rho}u^{n})_{j}$.

In [15], the authors show that a solution of the HJBQVI (1) can be constructed as the limit of $u^{\epsilon}$ as $\epsilon\downarrow 0$ where $u^{\epsilon}$ solves the so-called “penalized” problem

[TABLE]

where $(a)^{+}=\max(a,0)$. The penalty scheme is just a discretization of (19):

[TABLE]

where (recall that $d=1$)

[TABLE]

and we have introduced the shorthand $f_{j}^{n}(b)\coloneqq f(n\Delta t,x_{j},b)$, $\mu_{j}(b)\coloneqq\mu(x_{j},b)$, and $\sigma_{j}(b)\coloneqq\sigma(x_{j},b)$. To ensure consistency, we will always assume

[TABLE]

Remark 8**.**

(20)* defines $(2M+1)(N+1)$ “discrete” equations. In [4, Section 5.2], it is proved that these equations admit a unique solution. Morever, the solution at the $n$-th timestep $u^{n}$ can be computed by policy iteration or by value iteration. This is the motivation for discretizing (19) instead of the HJBQVI (1) directly since in the latter case, the resulting discrete equations may not admit any solutions, and even when they do, their computation may be nontrivial [4, Section 5.1].*

Note that the unique solution of the discrete equations (20) is defined only at grid points $(n\Delta t,x_{j})$. Following [11, Pg. 281], we extend it to all points in $[0,T]\times\mathbb{R}$ by

[TABLE]

where $E(n,j)\coloneqq[(n-\nicefrac{{1}}{{2}})\Delta t,(n+\nicefrac{{1}}{{2}})\Delta t)\times[(j-\nicefrac{{1}}{{2}})\Delta x,(j+\nicefrac{{1}}{{2}})\Delta x)$ and $\boldsymbol{1}_{G}$ is the indicator function of the set $G$. Note that the extension $u^{\rho}$ is a (discontinuous) simple function (recall Remark 5). We can now state the convergence result:

Theorem 9**.**

Suppose (H1)–(H4) and that $\mu$ and $\sigma$ are Lipschitz in $x$ uniformly in $b$ so that the HJBQVI (1) satisfies a (def.1)-comparison principle. Suppose also that (15) holds. Then, as $\rho\downarrow 0$, the solution $u^{\rho}$ of the penalty scheme converges locally uniformly to the unique bounded solution of the HJBQVI.

For the reader’s convenience, the comparison principle alluded to above is stated in Appendix A. We will spend the remainder of this subsection proving the above by appealing to Theorem 4 and establishing the monotonicity, stability, and nonlocal consistency of the penalty scheme.

First, note that the discrete equations (20) are equivalent to, by some algebra,

[TABLE]

where

[TABLE]

Next, we write the penalty scheme in the form (9) by taking

[TABLE]

where $[u]_{j}^{n}\coloneqq[u]_{(n\Delta t,x_{j})}$ and

[TABLE]

While the above only define $S$ and $\mathcal{N}^{\rho}$ at grid points $(n\Delta t,x_{j})$, it is understood that they are extended to all points in $[0,T]\times\mathbb{R}$ via the piecewise requirement (21).

4.1.1   Monotonicity

Due to our choice of upwind discretization, the penalty scheme is monotone by construction. We establish this below.

Monotonicity proof.

Let $(n\Delta t,x_{j})$ be an arbitrary grid point and let $u$ and $w$ be functions in $B([0,T]\times\mathbb{R})$ satisfying $u\geqslant w$ and $u_{j}^{n}=w_{j}^{n}$ where we have, as usual, employed the shorthand $u_{j}^{n}=u(n\Delta t,x_{j})$ and $w_{j}^{n}=w(n\Delta t,x_{j})$. If $n=N$, then

[TABLE]

Otherwise, defining $v\coloneqq w-u$, we have

[TABLE]

since the terms $v_{j}^{n+1}$ and $(L_{b}^{\rho}v^{n})_{j}$ are nonpositive. ∎

4.1.2   Stability

We will prove that a solution of the penalty scheme is bounded by $\|g\|_{\infty}+\|f\|_{\infty}T$. In the context of the optimal control problem (14), this bound has the natural interpretation of being the sum of the maximum possible continuous reward and the maximum possible reward obtained at the final time:

[TABLE]

Stability proof.

Let $u$ be a solution of the penalty scheme. Fix a particular timestep $n<N$ and let $j-$ be such that $u_{j-}^{n}=\min_{j}u_{j}^{n}$. It is easily verified that $(L_{b}^{\rho}u^{n})_{j-}\geqslant 0$. Therefore,

[TABLE]

Now, let $j+$ be such that $u_{j+}^{n}=\max_{j}u_{j}^{n}$. It is easily verified that $(L_{b}^{\rho}u^{n})_{j+}\leqslant 0$ and $(\mathcal{M}^{\rho}u^{n})_{j+}\leqslant u_{j+}^{n}$ (since $K\leqslant 0$ by (H3)). Therefore,

[TABLE]

By (25) and (26),

[TABLE]

and hence $\|u^{n}\|_{\infty}\leqslant\|u^{n+1}\|_{\infty}+\|f\|_{\infty}\Delta t$. By induction on $n$,

[TABLE]

and the desired result follows since $\|u^{N}\|_{\infty}\leqslant\|g\|_{\infty}$. ∎

4.1.3   Nonlocal consistency

To establish nonlocal consistency, we require three lemmas. The first lemma is a standard result from analysis:

Lemma 10**.**

Let $(a_{n})_{n}$, $(b_{n})_{n}$, and $(c_{n})_{n}$ be real sequences with $c_{n}\geqslant\min(a_{n},b_{n})$ (resp. $c_{n}\leqslant\min(a_{n},b_{n})$) for each $n$. Then,

[TABLE]

Proof.

If $c_{n}\geqslant\min(a_{n},b_{n})$ for each $n$, then

[TABLE]

Taking limits in the above and using the continuity of $(x,y)\mapsto\min(x,y)$ establishes the result. The case of $c_{n}\leqslant\min(a_{n},b_{n})$ is handled similarly. ∎

The second lemma justifies approximating the term $\sup_{b\in B}\{\cdot\}$ in the HJBQVI (1) by $\max_{b\in B^{\rho}}\{\cdot\}$. Recall that, as per (H4), $B^{\rho}$ converges to $B$ in the Hausdorff metric.

Lemma 11**.**

Let $Y$ be a compact metric space and $\zeta:Y\times B\rightarrow\mathbb{R}$ be a continuous function. Let $(\rho_{m})_{m}$ be a sequence of positive numbers converging to zero and $(y_{m})_{m}$ be a sequence taking values in $Y$ and converging to some $\hat{y}$ in $Y$. Then,

[TABLE]

Proof.

Let $\zeta_{m}\coloneqq\max_{b\in B^{\rho_{m}}}\zeta(y_{m},b)$. Since the sequence $(\zeta_{m})_{m}$ is bounded by $\|\zeta\|_{\infty}<\infty$ ($Y\times B$ is compact and $\zeta$ is continuous), it is sufficient to show that every convergent subsequence of $(\zeta_{m})_{m}$ converges to $\max_{b\in B}\zeta(\hat{y},b)$.

Therefore, consider an arbitrary convergent subsequence of $(\zeta_{m})_{m}$, and relabel it, with a slight abuse of notation, $(\zeta_{m})_{m}$. Choose $\hat{b}$ such that $\zeta(\hat{y},\hat{b})=\max_{b\in B}\zeta(\hat{y},b)$. Since $B^{\rho_{m}}\rightarrow B$, we can find a sequence $(b_{m})_{m}$ such that $b_{m}\in B^{\rho_{m}}$ for each $m$ and $b_{m}\rightarrow\hat{b}\in B$. Therefore,

[TABLE]

Furthermore, note that

[TABLE]

where we have obtained the last inequality by the continuity of $y\mapsto\max_{b\in B}\zeta(y,b)$.

To see that $y\mapsto\max_{b\in B}\zeta(y,b)$ is continuous, let $y_{0}\in Y$ and $\epsilon>0$ be arbitrary. Since $Y\times B$ is compact, $\zeta$ is uniformly continuous. Therefore, we can pick $\delta>0$ such that for all $b\in B$ and $y_{1}\in Y$ with $|y_{0}-y_{1}|<\delta$, $|\zeta(y_{0},b)-\zeta(y_{1},b)|<\epsilon$ and hence

[TABLE]

The last lemma characterizes $\mathcal{M}^{\rho}$ as a legitimate approximation of $\mathcal{M}$ as $\rho\downarrow 0$. To state it, we first need to introduce some notation. Recall that $\Delta t$ and $\Delta x$ are functions of the parameter $\rho$ (we write $\Delta t^{\rho}$ and $\Delta x^{\rho}$ to make this dependence explicit). Now, for a function $(\rho,n,j)\mapsto h_{j}^{n}(\rho)$ of the parameter $\rho$ and grid indices $n$ and $j$, we define

[TABLE]

The limit inferior is defined by $\liminf h_{j}^{n}(\rho)\coloneqq-\limsup(-h_{j}^{n}(\rho))$.

Lemma 12**.**

Let $(w^{\rho})_{\rho>0}$ be a family of uniformly bounded real-valued maps from $[0,T]\times\mathbb{R}$ with half-relaxed limits $\overline{w}$ and $\underline{w}$. Denote by $w^{n,\rho}\coloneqq(w^{\rho}(n\Delta t^{\rho},j\Delta x^{\rho}))_{j=-M}^{M}$ the image of $w^{\rho}(n\Delta t^{\rho},\cdot)$ on the spatial grid. Then, for any $(t,x)\in[0,T]\times\mathbb{R}$,

[TABLE]

Proof.

We first prove the leftmost inequality. Let $(\rho_{m},n_{m},j_{m})_{m}$ be an arbitrary sequence satisfying $\rho_{m}\downarrow 0$ and $(n_{m}\Delta t^{\rho_{m}},j_{m}\Delta x^{\rho_{m}})\rightarrow(t,x)$. Define $s_{m}\coloneqq n_{m}\Delta t^{\rho_{m}}$ and $y_{m}\coloneqq j_{m}\Delta x^{\rho_{m}}$ for brevity. Now, let $\delta>0$ and choose $z^{\delta}\in Z(t,x)$ such that

[TABLE]

By (H4), we can pick a sequence $(z_{m})_{m}$ such that $z_{m}\rightarrow z^{\delta}$ and $z_{m}\in Z^{\rho_{m}}(s_{m},y_{m})$ for each $m$. For brevity, we write $w_{j}^{m}$ for the quantity $w^{\rho_{m}}(s_{m},j\Delta x^{\rho_{m}})$. Noting that $w_{j}^{m}$ is the $j$-th entry of the vector $w^{n_{m},\rho_{m}}$ defined in the lemma statement,

[TABLE]

where $0\leqslant\alpha_{m}\leqslant 1$ and, for $m$ sufficiently large,

[TABLE]

(recall (18)). Therefore, by Lemma 10,

[TABLE]

Now, by (27), $(s_{m},k_{m}\Delta x^{\rho_{m}})\rightarrow(t,x+\Gamma(t,x,z^{\delta}))$ as $m\rightarrow\infty$. Therefore, by the definition of the half-relaxed limit $\underline{w}$,

[TABLE]

and

[TABLE]

Applying (29) and (30) to (28),

[TABLE]

Since $\delta$ is arbitrary, the desired result follows.

We now handle the rightmost inequality. Let $(\rho_{m},n_{m},j_{m})_{m}$ be a sequence as in the previous paragraph with $s_{m}$ and $y_{m}$ defined as before. Since $Z^{\rho_{m}}(s_{m},y_{m})$ is by definition a finite set (and hence compact), for each $m$, there exists a $z_{m}\in Z^{\rho_{m}}(s_{m},y_{m})$ such that

[TABLE]

where $0\leqslant\alpha_{m}\leqslant 1$ and $k_{m}$ satisfies (27). Now, consider a subsequence of $(z_{m})_{m}$ along which the limit superior

[TABLE]

is attained. With a slight abuse of notation, relabel this subsequence $(z_{m})_{m}$. By (H4), $(z_{m})_{m}$ is contained in a compact set and hence we can assume this subsequence was chosen to be convergent with limit $\hat{z}$. By (H4), the right hand side of the inequality

[TABLE]

approaches zero as $m\rightarrow\infty$ and hence we can conclude $\hat{z}\in Z(t,x)$. Therefore, similarly to how we established (28) and (31) in the previous paragraph,

[TABLE]

Remark 13**.**

The proof above relies heavily on the fact that $\operatorname{interp}$ is a linear interpolant and hence the “interpolation weights” $\alpha$ and $(1-\alpha)$ in (18) are nonnegative. A quadratic interpolant, for example, takes the form

[TABLE]

where it is not necessarily the case that the coefficients $\alpha$, $(1-\alpha-\beta)$, and $\beta$ are nonnegative. This suggests that higher order discretizations of $\mathcal{M}$ are generally precluded by the nonlocal consistency requirement.

We are now ready to establish the nonlocal consistency of the penalty scheme.

Nonlocal consistency proof.

Let $\varphi\in C^{1,2}([0,T]\times\mathbb{R})$. For brevity, we define $\varphi_{j}^{n}\coloneqq\varphi(n\Delta t,j\Delta x)$ and $\varphi^{n}=(\varphi_{0}^{n},\ldots,\varphi_{M}^{n})^{\intercal}$. Let $(w^{\rho})_{\rho>0}$ be a family of uniformly bounded real-valued maps from $[0,T]\times\mathbb{R}$ with half-relaxed limits $\overline{w}$ and $\underline{w}$. Let $(\rho_{m},s_{m},y_{m},\xi_{m})_{m}$ be an arbitrary sequence satisfying

[TABLE]

Without loss of generality, we will assume that the sequence is chosen such that $s_{m}=n_{m}\Delta t$ and $y_{m}=j_{m}\Delta x$ are grid points (we have omitted the dependence of $\Delta t$ and $\Delta x$ on $\rho_{m}$ in the notation). For brevity, we define $\alpha\coloneqq(t,x)$ and $\alpha_{m}\coloneqq(s_{m},y_{m})$.

By (22) and (23),

[TABLE]

where

[TABLE]

In the above, we have used the notation $w^{n,\rho}$ of Lemma 12. Now, by Lemma 11,

[TABLE]

Moreover, by Lemma 12,

[TABLE]

Suppose now that $t<T$. Since $s_{m}\rightarrow t$, we may assume $s_{m}<T$ (or, equivalently, $n_{m}<N$) for each $m$. In this case, taking limit inferiors of both sides of (33) and applying Lemma 10 yields

[TABLE]

Applying (35) and (36) to the above,

[TABLE]

where $D^{2}\varphi\equiv\varphi_{xx}$, $D\varphi\equiv(\varphi_{t},\varphi_{x})$, and $F$ is given by (4). In establishing the last equality in the above, we have used the fact that $F$ is continuous away from $t=T$ and hence $F=F_{*}$ there. Now, since $(\rho_{m},s_{m},y_{m},\xi_{m})_{m}$ is an arbitrary sequence satisfying (32), (37) implies

[TABLE]

which is exactly the nonlocal consistency inequality (12b) in the time-dependent case ($\alpha=(t,x)$). Symmetrically, we can establish the inequality

[TABLE]

which corresponds to the nonlocal consistency inequality (12a).

Suppose now that $t=T$. Since $s_{m}\rightarrow t$, it is possible that $s_{m}=T$ (or, equivalently, $n_{m}=N$) for one or more indices $m$ in the sequence. Therefore, by (33)

[TABLE]

An immediate consequence of the definition of $S_{m}^{(3)}$ in (34) is that

[TABLE]

Taking limit inferiors in (38) and applying Lemma 10, (35), (36), and (39),

[TABLE]

As in the previous paragraph, this implies the nonlocal consistency inequality (12b).

It remains to establish (12a) in the case of $t=T$. Since $\mathcal{M}g\leqslant g$ by (H1), it follows that $g(y_{m})=\max(g(y_{m}),\mathcal{M}g(y_{m}))$ for each $m$. By Remark 6, we can, without loss of generality, assume that $w^{\rho}$ is a solution of the scheme so that $w^{\rho}(N\Delta t,j\Delta x)=g(j\Delta x)$ for all $j$, corresponding to the terminal condition. Therefore,

[TABLE]

where in the last equality we have employed the fact that there is $O((\Delta x)^{2})$ error in approximating the intervention operator $\mathcal{M}$ by the discretized intervention operator $\mathcal{M}^{\rho}$ due to the linear interpolant. It follows that, letting

[TABLE]

we have $S_{m}^{(3)}=S_{m}^{(4)}$ whenever $n_{m}=N$. Therefore, by (33),

[TABLE]

Moreover, by Lemma 10 and Lemma 12,

[TABLE]

Taking limit superiors in (40) and applying Lemma 10, (35), (36), and (41),

[TABLE]

which establishes (12a), as desired. ∎

4.1.4   Higher dimensions

In higher dimensions, the operator $L_{b}$ may contain cross-derivatives which, in the context of the optimal control problem (14), corresponds to correlations between the components of $X^{t,x,\gamma}$. A naïve discretization of these cross-derivatives (e.g., using a standard five-point stencil) results in a nonmonotone scheme. Unfortunately, nonmonotone schemes are well-known to be capable of failing to converge to the viscosity solution [28].

It is possible to resolve these issues by discretizing $L_{b}$ using wide-stencils similarly to the implicit scheme described in [21, Corollary 5.1]. This results in a scheme that is first order accurate so long as the stencil is of length $O(\Delta x^{1/2})$. Since wide-stencils are well-understood, we simply mention that it is routine to extend the convergence analysis of this section to a wide-stencil version of the penalty scheme. In particular, it is trivial to extend Lemma 12 to higher dimensions by interpreting $j\equiv(j_{1},\ldots,j_{d})$ as a multi-index and $\operatorname{interp}(\cdot)$ as a monotone multi-dimensional linear interpolation (in this case, $x_{j}\equiv(x_{j_{1}},\ldots,x_{j_{d}})$ would be a point on a rectilinear grid, e.g., $(\{k\Delta x\}_{j=-M}^{M})^{d}$).

4.1.5   Infinite-horizon (steady state) case

If we assume the functions $f$, $Z$, $\Gamma$, and $K$ are independent of time (i.e., $f(t,x,b)=f(x,b)$, etc.), the infinite-horizon analogue of (1) is given by

[TABLE]

where $\beta>0$ is a positive discount factor. Note that in the above, $u$, $L_{b}u$, and $\mathcal{M}u$ are no longer functions of time and space but rather functions of space alone. Specifically, we interpret $\mathcal{M}$ above as $\mathcal{M}u(x)\coloneqq\sup_{z\in Z(x)}\{u(x+\Gamma(x,z))+K(x,z)\}$.

As with (1), (42) can be interpreted as the dynamic programming equation associated to the optimal control problem whose value function is given by (compare with (14))

[TABLE]

where

[TABLE]

We can adapt the scheme (20) to this setting by considering the discrete equations

[TABLE]

for $-M\leqslant j\leqslant M$, where we have written $\vec{u}\coloneqq(u_{-M},\ldots,u_{M})^{\intercal}$ to stress that the numerical solution is a vector in $\mathbb{R}^{2M+1}$. We interpret $\mathcal{M}^{\rho}$ in the above as $(\mathcal{M}^{\rho}\vec{u})_{j}\coloneqq\sup_{z\in Z^{\rho}(x_{j})}\{\operatorname{interp}(\vec{u},x_{j}+\Gamma(x_{j},z))+K(x_{j},z)\}$, with $Z^{\rho}$ satisfying a time-independent version of (H4). The analysis of this scheme is identical to that of (20), save that we obtain the stability bound $\|u\|_{\infty}\leqslant\beta^{-1}\|f\|_{\infty}$.

4.2   Semi-Lagrangian scheme

We now turn to a semi-Lagrangian scheme for the HJBQVI introduced in [4, Section 5.3], which can be used if the diffusion coefficient $\sigma$ does not depend on the control (i.e., $\sigma(x,b)=\sigma(x)$). As we will review shortly, the advantage of this scheme is that it only requires a single solution of a linear system per timestep (compare this with the penalty scheme, which requires an expensive iterative method; see Remark 8). As usual, we focus only on the one-dimensional case ($d=1$) since an extension to higher dimensions can be handled by wide-stencils (see Section 4.1.4). Roughly speaking, the semi-Lagrangian scheme

discretizes first order terms by tracing the path of a deterministic particle as determined by the Lagrangian derivative $\mathbb{D}$ (defined below) and 2. 2.

discretizes the intervention operator using information from the $(n+1)$-th timestep, adding a $O(\Delta t)$ term in order to “linearize” the scheme.

We start by discussing the first point (from which the scheme’s name is derived). Let $X^{x,b}$ denote a particle whose position at time $s$ is given by

[TABLE]

(in the above, $b$ is a fixed element of $B$; not a process taking values in $B$). For a “smooth enough” function $u$ of time and space, we consider the Lagrangian derivative

[TABLE]

where the approximation above is justified for $\Delta t>0$ small enough since in this case, we expect $X_{\Delta t}^{x,b}\approx x+\mu(x,b)\Delta t$. Substituting $t=n\Delta t$ and $x=x_{j}$ into (43) and using interpolation to approximate the term $u(t+\Delta t,x+\mu(x,b)\Delta t)$ yields

[TABLE]

As for the second point, we take

[TABLE]

as a discretization of the intervention operator, where we have used the shorthand $\sigma_{j}\coloneqq\sigma(x_{j})$. As we will see shortly, the $O(\Delta t)$ term in (45) allows us to express the scheme as a linear system at each timestep (at the cost of $O(\Delta t)$ discretization error).

Finally, we substitute (44) and (45) into the HJBQVI (1) to arrive at the semi-Lagrangian scheme:

[TABLE]

By some simple algebra, the discrete equations above are equivalent to

[TABLE]

where $(Au^{n})_{j}\coloneqq u_{j}^{n}-\frac{1}{2}\sigma_{j}^{2}(\mathcal{D}^{2}u^{n})_{j}\Delta t$. Note that we can also represent $A$ as the matrix

[TABLE]

interpreting $(Au^{n})_{j}$ as the $j$-th entry of the matrix-vector product $Au^{n}$. In particular, the matrix $A$ is strictly diagonally dominant and hence nonsingular.

Theorem 14**.**

If the diffusion coefficient $\sigma$ does not depend on the control (i.e., $\sigma(x,b)=\sigma(x)$) and the assumptions of Theorem 9 are satisfied, the solution of the semi-Lagrangian scheme converges locally uniformly to the unique bounded solution of the HJBQVI.

As with the penalty scheme, we write the semi-Lagrangian scheme in the form (9) by taking

[TABLE]

and

[TABLE]

Below, we handle the stability, monotonicity, and nonlocal consistency of the semi-Lagrangian scheme to establish Theorem 14.

4.2.1   Monotonicity

Monotonicity proof.

The proof is identical to that of Section 4.1.1, save that (24) becomes

[TABLE]

In the above, we have used the fact that the terms

[TABLE]

and $(\mathcal{D}^{2}v^{n})_{j}$ are nonpositive. ∎

4.2.2   Stability

Similar to the penalty scheme, we obtain a bound of $\|g\|_{\infty}+\|f\|_{\infty}T$ on the solution.

Stability proof.

Let $u$ be the solution of the semi-Lagrangian scheme. Fix a particular timestep $n<N$. Note that for any $j$, $\min_{j}u_{j}^{n+1}\leqslant\operatorname{interp}(u^{n+1},\cdot)\leqslant\max_{j}u_{j}^{n+1}$ and $(\mathcal{M}^{\rho}u^{n+1})_{j}\leqslant\max_{j}u_{j}^{n+1}$ (since $K\leqslant 0$ by (H3)). Now, let $j-$ and $j+$ be defined as in Section 4.1.2. It is readily verified that $(\mathcal{D}^{2}u^{n})_{j-}\geqslant 0$ and $(\mathcal{D}^{2}u^{n})_{j+}\leqslant 0$. Therefore,

[TABLE]

and

[TABLE]

The remainder of the proof proceeds as in Section 4.1.2. ∎

4.2.3   Nonlocal consistency

Nonlocal consistency proof.

Let $\varphi\in C^{1,2}([0,T]\times\mathbb{R})$. For brevity, we define $\varphi_{j}^{n}\coloneqq\varphi(n\Delta t,j\Delta x)$ and $\varphi^{n}=(\varphi_{0}^{n},\ldots,\varphi_{M}^{n})^{\intercal}$. Let $(\rho_{m},s_{m},y_{m},\xi_{m})_{m}$ be an arbitrary sequence satisfying (32) chosen such that $s_{m}=n_{m}\Delta t$ and $y_{m}=j_{m}\Delta x$ are grid points. Then,

[TABLE]

by a Taylor expansion. In the above, we used the fact that since $Q\rightarrow\infty$ as $\rho\downarrow 0$, we can always find $m$ large enough such that $y_{m}+\mu_{j_{m}}(b)\Delta t$ is a point inside $[-Q,Q]$.

The remaining arguments are a slight modification of those for the penalty scheme. In particular, we need only modify the proof by setting $\epsilon=0$ in (33) and replacing the definitions of $S_{m}^{(1)}$ and $S_{m}^{(2)}$ in (34) by

[TABLE]

and

[TABLE]

where we have used the symbol $\hat{\cdot}$ to distinguish new definitions from the old. By (46),

[TABLE]

Moreover, by Lemma 12,

[TABLE]

(compare with (36)). ∎

Appendix A Comparison principle for the HJBQVI

We now give a comparison principle for the HJBQVI. Recall that the HJBQVI is given by (6) when $\overline{\Omega}=[0,T]\times\mathbb{R}^{d}$ and $\mathcal{M}$ and $F$ are defined in (3) and (4).

Proposition 15**.**

Suppose (H1)–(H4) and that $\mu$ and $\sigma$ are Lipschitz in $x$ uniformly in $b$ (i.e., $\left|\mu(x,b)-\mu(y,b)\right|+\left|\sigma(x,b)-\sigma(y,b)\right|\leqslant\operatorname{const.}\left|x-y\right|$) Then, the HJBQVI satisfies a (def.1)-comparison principle in $B([0,T]\times\mathbb{R}^{d})$ (for any $d\geqslant 1$).

A proof of the above is found in [3, Theorem 2.13]; see also [32, Theorem 5.11] for a related comparison principle for solutions of polynomial growth. However, in both these works, the notion of viscosity solution used is the following:

Definition 16**.**

$u\in B_{\operatorname{loc}}(\overline{\Omega})$* is a subsolution of (6) if for all $\varphi\in C^{2}(\overline{\Omega})$ and $x\in\overline{\Omega}$ such that $u^{*}-\varphi$ has a local maximum (resp. minimum) at $x$, we have*

[TABLE]

Supersolutions are defined symmetrically.

Unlike Definition 1, the above does not use the envelopes $F_{*}$ and $F^{*}$. Fortunately, for the HJBQVI, these solution concepts are equivalent:

Proposition 17**.**

For the HJBQVI, Definitions 1 and 16 are equivalent.

We prove only the nontrivial direction. The proof is similar to [22, Remark 3.2].

Proof.

Let $u$ be a subsolution in the sense of Definition 1. Without loss of generality, assume $u$ is upper semicontinuous. Let $\varphi\in C^{1,2}([0,T]\times\mathbb{R}^{d})$ and $(T,x)$ be a maximum point of $u-\varphi$. Define $\psi(t,\cdot)\coloneqq\varphi(t,\cdot)+c(T-t)$ where $c$ is some positive constant, to be chosen later. Since $(T,x)$ is also a maximum point of $u-\psi$,

[TABLE]

By choosing $c$ large enough, the above reduces to a minimum of two terms so that

[TABLE]

Therefore, $u$ is also a subsolution in the sense of Definition 16. The case of supersolutions is handled similarly. ∎

Appendix B Overstepping error

We close with a discussion which, while not related to the theoretical convergence of the semi-Lagrangian scheme (Theorem 14), is of great importance for a practical implementation.

As discussed in Remark 7, a practical implementation should take $Q=\operatorname{const.}$ In this case, the point $x_{j}+\mu_{j}(b)\Delta t$ may not lie in the numerical domain $[-Q,Q]$ regardless of how small we take $\rho$. Therefore, if $Q=\operatorname{const.}$, the equality (46) fails to hold, and the semi-Lagrangian approximation introduces $O(1)$ “overstepping error”. Semi-Lagrangian schemes are well-known to exhibit this phenomenon (see, e.g., [31, 30]).

One way to avoid overstepping error is to modify the spatial grid so that the distance between boundary nodes approaches zero sublinearly as $\rho\downarrow 0$. For example, we may use a spatial grid satisfying (using Bachmann-Landau notation)

[TABLE]

We sketch the idea below.

Let

[TABLE]

If $x_{j}$ is a grid point satisfying $x_{j}<Q$, then

[TABLE]

Since $x_{M}-x_{M-1}=\omega(\rho)$ by (47), the above implies $x_{j}+\mu_{j}(b)\Delta t\leqslant Q$ for $\rho$ sufficiently small. By a symmetric argument, we conclude that if $|x_{j}|<Q$, then $|x_{j}\pm\mu_{j}(b)\Delta t|\leqslant Q$ for $\rho$ sufficiently small (i.e., there is no overstepping).

The $o(\sqrt{\rho})$ requirement in (47) is used to ensure that the error from interpolation vanishes. Indeed, if spatial grid points are not uniformly spaced, given a smooth function $\varphi$ of time and space, a computation similar to (46) yields

[TABLE]

for $|x_{j}|<Q$. Now, by the $o(\sqrt{\rho})$ requirement,

[TABLE]

The above computations suggest that, unsurprisingly, the local order of convergence at the boundaries is adversely affected by our new choice of grid. Therefore, it is worthwhile to point out that in the case of $\mu(-Q,\cdot)\geqslant 0$ and $\mu(Q,\cdot)\leqslant 0$, we can use the original uniform grid since no overstepping can occur.