Deep backward schemes for high-dimensional nonlinear PDEs

Introduction

This paper is devoted to the resolution in high dimension of nonlinear parabolic partial differential equations (PDEs) of the form

Due to the so called “curse of dimensionality”, the resolution of nonlinear PDEs in high dimension has always been a challenge for scientists. Until recently, only the BSDE (Backward Stochastic Differential Equation) approach first developed in [PP90] was available to tackle this problem: using the time discretization scheme proposed in [BT04], some effective algorithms based on regressions manage to solve non linear PDEs in dimension above 4 (see [GLW05, LGW06]). However this approach is still not implementable in dimension above 6 or 7 : the number of basis functions used for the regression still explodes with the dimension.

Quite recently some new methods have been developed for this problem, and several methodologies have emerged:

Some are based on the Feyman-Kac representation of the PDE. Branching techniques [HL+16] have been studied and shown to be convergent but only for small maturities and some small nonlinearities. Some effective techniques based on nesting Monte Carlo have been studied in [War18a, War18]: the convergence is proved for semi-linear equations. Still based on this Feyman-Kac representation some machine learning techniques permitting to solve a fixed point problem have been used recently in [CWNMW19]: numerical results show that it is efficient and some partial demonstrations justify why it is effective.

Multilevel Picard methods have been developed in [E+18] and [Hut+18] with algorithms based on Picard iterations, multi-level techniques and automatic differentiation. These methods permit to handle some high dimensional PDEs with non linearity in $u$ and its gradient $D_{x}u$, with convergence results as well as numerous numerical examples showing their efficiency in high dimension.

Another class of methods is based on the BSDE approach and the curse of dimensionality issue is partially avoided by using some machine learning techniques. The pioneering papers [HJE18, EHJ17] propose a neural-networks based technique called Deep BSDE, which was the first serious attempt for using machine learning methods to solve high dimensional PDEs. Based on an Euler discretization of the forward underlying SDE ${\cal X}_{t}$, the idea is to view the BSDE as a forward SDE, and the algorithm tries to learn the values $u$ and $z$ $=$ $\sigma^{\scriptscriptstyle{\intercal}}Du$ at each time step of the Euler scheme by minimizing a global loss function between the forward simulation of $u$ till maturity $T$ and the target $g({\cal X}_{T})$. This deep learning approximation has been extended to the case of fully nonlinear PDE and second order BSDE in [BEJ19].

At last, using some machine learning representation of the solution, [SS18] proposes with the so-called Deep Galerkin Method to use the automatic numerical differentiation of the solution to solve the PDE on a finite domain. The authors prove the convergence of their method but without information on the rate of convergence.

Like the second methodology, our approach relies on BSDE representation of the PDE and deep learning approximations: we first discretize the BSDE associated to the PDE by an Euler scheme, but in contrast with [EHJ17], we adopt a classical backward resolution technique. On each time step, we propose to use some machine learning techniques to estimate simultaneously the solution and its gradient by minimizing a loss function defined recursively by backward induction, and solving this local problem by a stochastic gradient algorithm. Two different schemes are designed to deal with the local problems:

The first one tries the estimate the solution and its gradient by a neural network.

The second one tries only to approximate the solution by a neural network while its gradient is estimated directly with some numerical differentiation techniques.

The proposed methodology is then extended to solve some variational inequalities, i.e., free boundary problems related to optimal stopping problems. We mention that the related recent paper [BCJ19] also proposes deep learning method for solving optimal stopping problems, but differently from our method, it relies on the approximation of (randomised) stopping decisions with a sequence of multilayer feedforward neural networks.

Convergence analysis of the two schemes for PDEs and variational inequalities is provided and shows that the approximation error goes to zero as we increase the number of time steps and the number of neurons/layers whenever the gradient descent method used to solve the local problems is not trapped in a local minimum. Notice that similar convergence result for the deep BSDE method has been also obtained in [HL18] with a posteriori error estimation of the solution in terms of the universal approximation capability of global neural networks.

In the last part of the paper, we test our algorithms on different examples. When the solution is easy to represent by a neural network, we can solve the problem in quite high dimension (at least $50$ in our numerical tests). We show that the proposed methodology improves the algorithm proposed in [HJE18] that sometimes does not converge or is trapped in a local minimum far away from the true solution. We then show that when the solution has a very complex structure, we can still solve the problem but only in moderate dimension: the neural network used is not anymore able to represent the solution accurately in very high dimension. Finally, we illustrate numerically that the method is effective to solve some system of variational inequalities: we consider the problem of American options and show that it can be solved very accurately in high dimension (we tested until $40$).

The outline of the paper is organized as follows. In Section 2, we give a brief and useful reminder for neural networks. We describe in Section 3 our two numerical schemes and compare with the algorithm in [HJE18]. Section 4 is devoted to the convergence analysis of our machine learning algorithms, and we present in Section 5 several numerical tests.

Neural networks as function approximators

Multilayer (also called deep) neural networks are designed to approximate unknown or large class of functions. In contrast to additive approximation theory with weighted sum over basis functions, e.g. polynomials, neural networks rely on the composition of simple functions, and appear to provide an efficient way to handle high-dimensional approximation problems, in particular thanks to the increase in computer power for finding the “optimal” parameters by (stochastic) gradient descent methods.

We denote by $\Phi_{{}_{m}}(.;\theta)$ the neural network function defined in (2.1), and by ${\cal N}{\cal N}_{d,d_{1},L,m}^{\varrho}(\Theta_{m})$ the set of all such neural networks $\Phi_{{}_{m}}(.;\theta)$ for $\theta$ $\in$ $\Theta_{m}$, and set

as the class of all neural networks within a fixed structure given by $d$, $d_{1}$, $L$ and $\varrho$.

The fundamental result of Hornick et al. [HSW89] justifies the use of neural networks as function approximators:

Moreover, we have a universal approximation result for the derivatives in the case of a single hidden layer, i.e. $L$ $=$ $2$, and when the activation function is a smooth function, see [HSW90].

Deep learning-based schemes for semi-linear PDEs

The DBSDE algorithm proposed in [HJE18, EHJ17] starts from the BSDE representation (3.1) of the solution to (1.1), but rewritten in forward form as:

The forward process ${\cal X}$ in equation (1.3), when it is not simulatable, is numerically approximated by an Euler scheme $X$ $=$ $X^{\pi}$ on a time grid: $\pi$ $=$ $\{t_{0}=0<t_{1}<\ldots<t_{N}=T\}$, with modulus $|\pi|$ $=$ $\max_{i=0,\ldots,N-1}\Delta t_{i}$, $\Delta t_{i}$ $:=$ $t_{i+1}-t_{i}$, and defined as

where we set $\Delta W_{t_{i}}$ $:=$ $W_{t_{i+1}}-W_{t_{i}}$. To alleviate notations, we omit the dependence of $X$ $=$ $X^{\pi}$ on the time grid $\pi$ as there is no ambiguity (recall that we use the notation ${\cal X}$ for the forward diffusion process). The approximation of equation (1.1) is then given formally from the Euler scheme associated to the forward representation (3.2) by

one computes estimations ${\cal U}_{i}$ of $u(t_{i};X_{t_{i}})$ by forward induction via:

for $i$ $=$ $0,\ldots,N-1$. This algorithm forms a global deep neural network composed of the neural networks (3.7) of each period, by taking as input data (in machine learning language) the paths of $(X_{t_{i}})_{i=0,\ldots,N}$ and $(W_{t_{i}})_{i=0,\ldots,N}$, and giving as output ${\cal U}_{N}$ $=$ ${\cal U}_{N}(\theta)$, which is a function of the input and of the total set of parameters $\theta$ $=$ $({\cal U}_{0},\theta_{0},\ldots,\theta_{N-1})$. The output aims to match the terminal condition $g(X_{t_{N}})$ of the BSDE, and one then optimizes over the parameter $\theta$ the expected square loss function:

This is obtained by stochastic gradient descent-type (SGD) algorithms relying on training input data.

2 New schemes: DBDP1 and DBDP2

The proposed scheme is defined from a backward dynamic programming type relation, and has two versions:

Initialize from an estimation $\widehat{\cal U}_{N}^{(1)}$ of $u(t_{N},.)$ with $\widehat{\cal U}_{N}^{(1)}$ $=$ $g$

Then, update: $\widehat{\cal U}_{i}^{(1)}$ $=$ ${\cal U}_{i}(.;\theta_{i}^{*})$, and set $\widehat{\cal Z}_{i}^{(1)}$ $=$ ${\cal Z}_{i}(.;\theta_{i}^{*})$.

Initialize with $\widehat{\cal U}_{N}^{(2)}$ $=$ $g$

For $i$ $=$ $N-1,\ldots,0$, given $\widehat{\cal U}_{i+1}^{(2)}$, use a deep neural network ${\cal U}_{i}(.;\theta)$ $\in$ ${\cal N}{\cal N}_{d,1,L,m}^{\varrho}(\Theta_{m})$, and compute (by SGD) the minimizer of the expected quadratic loss function

where $\hat{D}_{x}{\cal U}_{i}(.;\theta)$ is the numerical differentiation of ${\cal U}_{i}(.;\theta)$. Then, update: $\widehat{\cal U}_{i}^{(2)}$ $=$ ${\cal U}_{i}(.;\theta_{i}^{*})$, and set $\widehat{\cal Z}_{i}^{(2)}$ $=$ $\sigma^{\scriptscriptstyle{\intercal}}(t_{i},.)\hat{D}_{x}{\cal U}_{i}(.;\theta_{i}^{*})$.

For the first version of the scheme, one can use independent neural networks, respectively for the approximation of $u(t_{i},.)$ and for the approximation of $\sigma^{\scriptscriptstyle{\intercal}}(t_{i},.)D_{x}u(t_{i},.)$. In other words, the parameters are divided into a pair $\theta$ $=$ $(\xi,\eta)$ and we consider neural networks ${\cal U}_{i}(.;\xi)$ and ${\cal Z}_{i}(.;\eta)$. $\Box$

In the sequel, we refer to the first and second version of the new scheme above as DBDP1 and DBDP2, where the acronym DBDP stands for deep learning backward dynamic programming.

The intuition behind DBDP1 and DBDP2 is the following. For simplicity, take $f$ $=$ , so that $F(t,x,y,z,h,\Delta)$ $=$ $y+z^{\scriptscriptstyle{\intercal}}\Delta$. The solution $u$ to the PDE (1.1) should then approximately satisfy (see (3.5))

Consider the first scheme DBDP1, and suppose that at time $i+1$, $\widehat{\cal U}_{i+1}^{(1)}$ is an estimation of $u(t_{i+1,.})$. The quadratic loss function at time $i$ is then approximately equal to

Therefore, by minimizing over $\theta$ this quadratic loss function, via SGD based on simulations of $(X_{t_{i}},X_{t_{i+1}},\Delta W_{t_{i}})$ (called training data in the machine learning language), one expects the neural networks ${\cal U}_{i}$ and ${\cal Z}_{i}$ to learn/approximate better and better the functions $u(t_{i},.)$ and $\sigma^{\scriptscriptstyle{\intercal}}(t_{i},)D_{x}u(t_{i},)$ in view of the universal approximation theorem [HSW90]. Similarly, the second scheme DPDP2, which uses only neural network on the value functions, learns $u(t_{i},.)$ by means of the neural network ${\cal U}_{i}$, and $\sigma^{\scriptscriptstyle{\intercal}}(t_{i},)D_{x}u(t_{i},)$ via $\sigma^{\scriptscriptstyle{\intercal}}(t_{i},)\hat{D}_{x}{\cal U}_{i}$. The rigorous arguments for the convergence of these schemes will be derived in the next section.

The advantages of our two schemes, compared to the Deep BSDE algorithm, are the following:

by decomposing the global problem into smaller ones, we may expect to help the gradient descent method to provide estimations closer to the real solution. The memory needed in [HJE18] can be a problem when taking too many time steps.

at each time step, we initialize the weights and bias of the neural network to the weights and bias of the previous time step treated : this trick is commonly used in iterative solvers of PDE, and allows us to start with a value close to the solution, hence avoiding local minima which are too far away from the true solution. Besides the number of gradient iterations to achieve is rather small after the first resolution step.

The small disadvantage is due to the Tensorflow structure. As it is done in python, the global graph creation takes much time as it is repeated for each time step and the global resolution is a little bit time consuming : as the dimension of the problem increases, the time difference decreases and it becomes hard to compare the computational time for a given accuracy when the dimension is above 5.

3 Extension to variational inequalities: scheme RDBDP

Let us consider a variational inequality in the form

which arises, e.g., in optimal stopping problem and American option pricing in finance. It is known, see e.g. [EK+97], that such variational inequality is related to reflected BSDE of the form

where $K$ is an adapted non-decreasing process satisfying

The extension of our DBDP1 scheme for such variational inequality, and refereed to as RDBDP scheme, becomes

Initialize $\widehat{\cal U}_{N}$ $=$ $g$

Then, update: $\widehat{\cal U}_{i}$ $=$ \max\big{[}{\cal U}_{i}(.;\theta_{i}^{*}),g], and set $\hat{\cal Z}_{i}$ $=$ ${\cal Z}(.;\theta_{i}^{*})$.

Convergence analysis

The main goal of this section is to prove convergence of the DBDP schemes towards the solution $(Y,Z)$ to the BSDE (3.1) (or reflected BSDE (3.11) for variational inequalities), and to provide a rate of convergence that depends on the approximation errors by neural networks.

We assume the standard Lipschitz conditions on $\mu$ and $\sigma$, which ensures the existence and uniqueness of an adapted solution ${\cal X}$ to the forward SDE (1.3) satisfying for any $p$ $>$ $1$,

for some constant $C_{p}$ depending only on $p$, $b$, $\sigma$ and $T$. Moreover, we have the well-known error estimate with the Euler scheme $X$ $=$ $X^{\pi}$ defined in (3.4) with a time grid $\pi$ $=$ $\{t_{0}=0<t_{1}<\ldots<t_{N}=T\}$, with modulus $|\pi|$ s.t. $N|\pi|$ is bounded by a constant depending only on $T$ (hence independent of $N$):

Here, the standard notation $O(|\pi|)$ means that $\limsup_{|\pi|\rightarrow 0}\;|\pi|^{-1}O(|\pi|)$ $<$ $\infty$.

We shall make the standing usual assumptions on the driver $f$ and the terminal data $g$.

(H1) (i) There exists a constant $[f]_{{}_{L}}>0$ such that the driver $f$ satisfies:

(ii) The function $g$ satisfies a linear growth condition.

Recall that Assumption (H1) ensures the existence and uniqueness of an adapted solution $(Y,Z)$ to (3.1) satisfying

From the linear growth condition on $f$ in (H1), and (4.1), we also see that

Moreover, we have the standard $L^{2}$-regularity result on $Y$:

Let us also introduce the $L^{2}$-regularity of $Z$:

Let us first investigate the convergence of the scheme DBDP1 in (3.8), and define (implicitly)

for $i$ $=$ $0,\ldots,N-1$. Notice that $\widehat{\cal V}_{t_{i}}$ is well-defined for $|\pi|$ small enough (recall that $f$ is Lipschitz) by a fixed point argument. By the Markov property of the discretized forward process $(X_{t_{i}})_{i=0,\ldots,N}$, we note that there exists some deterministic functions $\hat{v}_{i}$ and $\overline{{\hat{z}_{i}}}$ s.t.

Let us now define a measure of the (squared) error for the DBDP1 scheme by

Our first main result gives an error estimate of the DBDP1 scheme in terms of the $L^{2}$-approximation errors of $\hat{v}_{i}$ and $\overline{{\hat{z}_{i}}}$ by neural networks ${\cal U}_{i}$ and ${\cal Z}_{i}$, $i=0,\ldots,N-1$, assumed to be independent (see Remark 3.1), and defined as

(Consistency of DBDP1) Under (H1), there exists a constant $C>0$, independent of $\pi$, such that

The error contributions for the DBDP1 scheme in the r.h.s. of estimation (4.16) consists of four terms. The first three terms correspond to the time discretization of BSDE, similarly as in [BT04], [GLW05], namely (i) the strong approximation of the terminal condition (depending on the forward scheme and the terminal data $g$), and converging to zero, as $|\pi|$ goes to zero, with a rate $|\pi|$ when $g$ is Lipschitz by (4.2) (see [Avi09] for irregular $g$), (ii) the strong approximation of the forward Euler scheme, and the $L^{2}$-regularity of $Y$, which gives a convergence of order $|\pi|$, (iii) the $L^{2}$-regularity of $Z$, which converges to zero, as $|\pi|$ goes to zero, with a rate $|\pi|$ when $g$ is Lipschitz. Finally, the better the neural networks are able to approximate/learn the functions $\hat{v}_{i}$ and $\overline{{\hat{z}_{i}}}$ at each time $i$ $=$ $0,\ldots,N-1$, the smaller is the last term in the error estimation. Moreover, given a prescribed accuracy for the neural network approximation error, the number of parameters of the employed deep neural networks grows at most polynomially in the PDE dimension, as recently proved in [Hut+19] in the case of semi-linear heat equations. $\Box$

In the following, $C$ will denote a positive generic constant independent of $\pi$, and that may take different values from line to line.

Step 1. Fix $i$ $\in$ $\{0,\ldots,N-1\}$, and observe by (3.1), (4.9) that

By using Young inequality: $(a+b)^{2}$ $\leq$ $(1+\gamma\Delta t_{i})a^{2}$ $+$ $(1+\frac{1}{\gamma\Delta t_{i}})b^{2}$ for some $\gamma$ $>$ to be chosen later, Cauchy-Schwarz inequality, the Lipschitz condition on $f$ in (H1), and the estimation (4.2) on the forward process, we then have

where we use in the last inequality the $L^{2}$-regularity (4.6) of $Y$.

Recalling the definition of $\bar{Z}$ as a $L^{2}$-projection of $Z$, we observe that

By multiplying equation (3.1) between $t_{i}$ and $t_{i+1}$ by $\Delta W_{t_{i}}$, and using Itô isometry, we have together with (4.9)

By Cauchy-Schwarz inequality, and law of iterated conditional expectations, this implies

Then, by plugging (4.23) and (4.28) into (4.21), and choosing $\gamma$ $=$ $8d[f]^{2}_{{}_{L}}$, we have

Step 2. By using Young inequality in the form: $(a+b)^{2}$ $\geq$ $(1-|\pi|)a^{2}$ $+$ $(1-\frac{1}{|\pi|})b^{2}$ $\geq$ $(1-|\pi|)a^{2}$ $-$ $\frac{1}{|\pi|}b^{2}$, we have

By plugging this last inequality into (4.32), we then get for $|\pi|$ small enough

From discrete Gronwall’s lemma (or by induction), and recalling the terminal condition $Y_{t_{N}}$ $=$ $g({\cal X}_{T})$, $\widehat{\cal U}_{i}^{(1)}(X_{t_{N}})$ $=$ $g(X_{T})$, the definition $\varepsilon^{Z}(\pi)$ of the $L^{2}$-regularity of $Z$, and (4.5), this yields

Step 3. Fix $i$ $\in$ $\{0,\ldots,N-1\}$. By using relation (4.11) in the expression of the expected quadratic loss function in (3.8), and recalling the definition of $\overline{{\widehat{Z}_{t_{i}}}}$ as a $L^{2}$-projection of $\widehat{Z}_{t}$, we have for all parameters $\theta$ $=$ $(\xi,\eta)$ of the neural networks ${\cal U}_{i}(.;\xi)$ and ${\cal Z}_{i}(.;\eta)$

By using Young inequality: $(a+b)^{2}$ $\leq$ $(1+\gamma\Delta t_{i})a^{2}$ $+$ $(1+\frac{1}{\gamma\Delta t_{i}})b^{2}$, together with the Lipschitz condition on $f$ in (H1), we clearly see that

On the other hand, using Young inequality in the form: $(a+b)^{2}$ $\geq$ $(1-\gamma\Delta t_{i})a^{2}$ $+$ $(1-\frac{1}{\gamma\Delta t_{i}})b^{2}$ $\geq$ $(1-\gamma\Delta t_{i})a^{2}$ $-$ $\frac{1}{\gamma\Delta t_{i}}b^{2}$, together with the Lipschitz condition on $f$, we have

By choosing $\gamma$ $=$ $4[f]_{{}_{L}}^{2}$, this yields

For $|\pi|$ small enough, and recalling (4.10), this implies

Plugging this last inequality into (4.39), we obtain

which proves the consistency of the $Y$-component in (4.16).

Step 5. Let us finally prove the consistency of the $Z$-component. From (4.23) and (4.28), we have for any $i$ $=$ $0,\ldots,N-1$:

By summing over $i$ $=$ $0,\ldots,N-1$, we get (recall (4.5))

where we change the indices in the last summation. Now, from (4.21), (4.34), we have

We now choose $\gamma$ $=$ $24d[f]^{2}_{{}_{L}}$ so that $\frac{8d[f]^{2}_{{}_{L}}}{\gamma}(1+\gamma|\pi|)/(1-|\pi|)$ $\leq$ $1/2$ for $|\pi|$ small enough, and by plugging into (4.52), we obtain (note also that \big{[}(1+\gamma|\pi|)/(1-|\pi|)-1\big{]} $=$ $O(|\pi|)$):

where we used (4.32) and (4.5) in the second inequality, and (4.47) and (4.49) in the last inequality.

and using (4.47), (4.62), we obtain after summation over $i$ $=$ $0,\ldots,N-1$, the required error estimate for the $Z$-component as in (4.49), and this ends the proof. $\Box$

2 Convergence of DBDP2

We shall consider neural networks with one hidden layer, $m$ neurons with total variation smaller than $\gamma_{m}$ (see Section 2), a $C^{3}$ activation function $\varrho$ with linear growth condition, and bounded derivatives, e.g., a sigmoid activation function, or a $\tanh$ function: this class of neural networks is then represented by the parametric set of functions

for some sequence $(\gamma_{m})_{m}$ converging to $\infty$, as $m$ goes to infinity, and such that

Notice that the neural networks in ${\cal N}{\cal N}_{d,1,2,m}^{\varrho}(\Theta_{m}^{\gamma})$ have their first, second and third derivatives uniformly bounded w.r.t. the state variable $x$. More precisely, there exists some constant $C$ depending only on $d$ and the derivatives of $\varrho$ s.t. for any ${\cal U}$ $\in$ ${\cal N}{\cal N}_{d,1,2,m}^{\varrho}(\Theta_{m}^{\gamma})$,

Let us investigate the convergence of the scheme DBDP2 in (3.9) with neural networks in ${\cal N}{\cal N}_{d,1,2,m}^{\varrho}(\Theta_{m}^{\gamma})$, and define for $i$ $=$ $0,\ldots,N-1$:

A measure of the (squared) error for the DBDP2 scheme is defined similarly as in DBDP1 scheme:

Our second main result gives an error estimate of the DBDP2 scheme in terms of the $L^{2}$-approximation errors of $\hat{v}_{i}$ and its derivative (which exists under assumption detailed below) by neural networks ${\cal U}_{i}$ $\in$ ${\cal N}{\cal N}_{d,1,2,m}^{\varrho}(\Theta_{m}^{\gamma})$, $i=0,\ldots,N-1$, and defined as

which are expected to be small in view of the universal approximation theorem (II), see discussion in Remark 4.2.

We also require the additional conditions on the coefficients:

(Consistency of DBDP2) Under (H1)-(H2), there exists a constant $C>0$, independent of $\pi$, such that

Proof. For simplicity of notations, we assume $d$ $=$ $1$, and only detail the arguments that differ from the proof of Theorem 4.16. From (4.67), and the Euler scheme (3.4), we have

where we use integration by parts in the second equality. Similarly, we have

Denoting by $\hat{f}_{i}(x)$ $=$ $f(t_{i},x,\hat{v}_{i}(x),\overline{{\hat{z}_{i}}}(x))$, it follows by the implicit function theorem, and for $|\pi|$ small enough, that $\hat{v}_{i}$ is $C^{1}$ with derivative given by

Under (H2), by the linear growth condition on $\sigma$, and using the bounds on the derivatives of the neural networks in ${\cal N}{\cal N}_{d,1,2,m}^{\varrho}(\Theta_{m}^{\gamma})$ in (4.66), we then have

Next, by the same arguments as in Steps 3 and 4 in the proof of Theorem 4.1 (see in particular (4.47)), we have for $|\pi|$ small enough,

for all $\theta$ $\in$ $\Theta^{N}$, and then with (4.78), and by definition of $\varepsilon_{i}^{NN,v,2}$:

On the other hand, by the same arguments as in Steps 1 and 2 in the proof of Theorem 4.1 (see in particular (4.39)), we have

Plugging (4.81) into this last inequality, together with (4.65), gives the required estimation (4.69) for the $Y$-component. Finally, by following the same arguments as in Step 5 in the proof of (4.1), we obtain the estimation (4.69) for the $Z$-component. $\Box$

The universal approximation theorem (II) [HSW90] is valid on compact sets, and one cannot conclude a priori that the error of network approximation $\varepsilon_{i}^{{\cal N},m}$ converge to zero as $m$ goes to infinity. Instead, we have to proceed into two steps:

where we set $\Delta_{i}(x;\theta)$ $:=$ $|\hat{v}_{i}(x)-{\cal U}_{i}(x;\theta)|^{2}$ + \Delta t_{i}\big{|}\sigma^{\scriptscriptstyle{\intercal}}(t_{i},x)\big{(}D_{x}\hat{v}_{i}(x)-D_{x}{\cal U}_{i}(x;\theta)\big{)}\big{|}^{2}.

Consider an increasing family of neural networks $\Theta_{m}^{\gamma^{N-1}}$ $\subset$ $\ldots$ $\subset$ $\Theta_{m}^{\gamma^{i}}$ $\subset$ $\ldots$ $\subset$ $\Theta_{m}^{\gamma^{0}}$ on which to minimize the approximation errors by backward induction at times $t_{i}$, $i$ $=$ $N-1,\ldots,0$, and where, $\gamma_{m}^{i}$ is defined by

The localized approximation error at time $t_{i}$, for $0\leq i\leq N-1$, should then be rewritten as

for some positive constant $C$ independent of $m,\pi$. We deduce by Cauchy-Schwarz and Chebyshev’s inequalities that for all $K$ $>$ , and $\theta$ $\in$ $\Theta_{m}^{\gamma^{i}}$, $i$ $=$ $0,\ldots,N-1$,

where we used (4.1) in the last inequality. This shows that

and thus, in theory, the error $\varepsilon_{i,N}^{{\cal N},m}$ can be made arbitrary small by suitable choices of large $m$ and $K$. $\Box$

3 Convergence of RDBDP

In this paragraph, we study the convergence of machine learning schemes for the variational inequality (3.10).

We first consider the case when $f$ does not depend on $z$, so that the component $Y_{t}$ $=$ $u(t,{\cal X}_{t})$ solution to the reflected BSDE (3.11) admits a Snell envelope representation, and we shall focus on the error on $Y$ by proposing an alternative to scheme (3.14), refereed to as RDBDPbis scheme, which only uses neural network for learning the function $u$:

Initialize $\widehat{\cal U}_{N}$ $=$ $g$

Then, update: $\widehat{\cal U}_{i}$ $=$ \max\big{[}{\cal U}_{i}(.;\theta_{i}^{*}),g].

Let us also define from the scheme (4.88)

(Case $f$ independent of $z$: Consistency of RDBDPbis) Let Assumption (H1) hold, with $g$ Lipschitz. Then, there exists a constant $C>0$, independent of $\pi$, such that

which is of the same order than the error estimate in Theorem 4.1 when $g$ is Lipschitz. $\Box$

Proof. Let us introduce the discrete-time approximation of the reflected BSDE

Fix $i$ $=$ $0,\ldots,N-1$. From (4.89), (4.93), we have

from the Lipschitz condition on $f$ in (H1), and then for $|\pi|$ small enough

By Minkowski inequality, this yields for all $\theta$

and the expected squared loss function of the DBDP3 scheme can be written as

From the Lipschitz condition on $f$, and by Minkowski inequality, we have for all $\theta$

We finally turn to the general case when $f$ may depend on $z$, and study the convergence of the RDBDP scheme (3.14) towards the variational inequality (3.10) related to the solution $(Y,Z)$ of the reflected BSDE (3.11) by showing an error estimate for

The result is obtained under one of the following additional assumptions

(H3) $g$ is $C^{1}$, and $g$, $D_{x}g$ are Lipschitz.

(H4) $\sigma$ is $C^{1}$, with $\sigma$, $D_{x}\sigma$ both Lipschitz, and $g$ is $C^{2}$, with $g$, $D_{x}g$, $D_{x}^{2}g$ all Lipschitz.

(Consistency of RDBDP) Let Assumption (H1) hold. There exists a constant $C>0$, independent of $\pi$, such that

with $\varepsilon(\pi)$ $=$ $O(|\pi|^{\frac{1}{2}})$ under (H3), and $\varepsilon(\pi)$ $=$ $O(|\pi|)$ under (H4).

Proof. Let us introduce the discrete-time approximation of the reflected BSDE

with $\varepsilon(\pi)$ $=$ $O(|\pi|^{\frac{1}{2}})$ under (H3), and $\varepsilon(\pi)$ $=$ $O(|\pi|)$ under (H4).

Fix $i$ $=$ $0,\ldots,N-1$. By writing that

and proceeding similarly as in Step 1 in the proof of Theorem 4.1, we have by Young inequality and Lipschitz condition on $f$

From (4.106), (4.109), Cauchy-Schwarz inequality, and law of iterated conditional expectations, we have similarly as in Step 1 in the proof of Theorem 4.1:

Then, by plugging into (4.113) and choosing $\gamma$ $=$ $4d[f]^{2}_{{}_{L}}$, we have for $|\pi|$ small enough:

Next, by using Young inequality as in Step 2 in the proof of Theorem 4.1, we obtain for all $\theta$ $=$ $(\xi,\zeta)$:

and the expected squared loss function of the RDBDP scheme can be written as

By the same arguments as in Step 3 in the proof of Theorem 4.1, using Lipschitz condition on $f$ and Young inequality, we show that for all $\theta$ $=$ $(\xi,\eta)$

Combining with (4.110), this proves the error estimate (4.108) for the $Y$-component. The error estimate (4.108) for the $Z$-component is proved along the same arguments as in Step 5 in the proof of Theorem 4.1, and is omitted here. $\Box$

Numerical results

In the first two subsections, we compare our schemes DBDP1 (3.8), DBDP2 (3.9) and the scheme proposed by [HJE18] on some examples of PDEs and BSDEs.

We first test our algorithms on some PDEs with bounded solutions and quite a simple structure (see section 5.1), and then try to solve some PDEs with unbounded solutions and more complex structures (see section 5.2). Our goal is to emphasize that solutions with simple structure easily represented by a neural network can be evaluated by our method even in very high-dimension, whereas the solution with complex structure can only be evaluated in moderate dimension.

Finally, we apply the scheme described in section 3.3 to an American option problem and show its accuracy in high dimension (see section 5.3).

If not specified, we use in the sequel a fully connected feedforward network with two hidden layers, and $d+10$ neurons on each hidden layer, to implement our schemes (3.8) and (3.9). We choose tanh as activation function for the hidden layers in order to avoid some explosion while calculating the numerical gradient $Z$ in scheme (3.9) and choose identity function as activation function for the output layer. We renormalize the data before entering the network. We use Adam Optimizer, implemented in TensorFlow and mini-batch with $1000$ trajectories for the stochastic gradient descent.

We begin with a simple example in dimension one. It is not hard to find test cases where the scheme proposed in [HJE18] fails even in dimension one. In fact the latter scheme works well for small maturities and with a starting point close to the solution.

It is always interesting to start by testing schemes in dimension one as one can easily compare graphically the numerical results to the theoretical solution. Then we take some examples in higher dimensions and show that our method seems to work well when the dimension increases higher.

We take the following parameters for the BSDE problem defined by (1.3) and (3.1):

for which, the explicit analytic solution is equal to $u(t,x)=e^{\frac{T-t}{2}}\cos(x)$.

We want to estimate the solution $u$ and its gradient $D_{x}u$ from our schemes. This example is interesting, because with $T=1$, the method proposed in [HJE18], initializing $u(0,.)$ as the solution of the associated linear problem associated ($f=0$) and randomly initializing $D_{x}u(0,.)$ works very well. However, for $T=2$, the method in [HJE18] always fails on our test whatever the choice of the initialization: the algorithm is either trapped in a local minimum when the initial learning rate associated to the gradient method is too small or explodes when the learning rate is taken higher. This numerical failure is not dependent on the considered network: using some LSTM networks as in [CWNMW19] gives the same result.

Because of the high non-linearity, we discretize the BSDE using $N$ $=$ $240$ time steps, and implemented hidden layers with $d+10$ $=$ $11$ neurons. Figure 1 (resp. Figure 2) depicts the estimated functions $u(t,.)$ and $D_{x}u(t,.)$ estimated from DBDP1 (resp. DBDP2) scheme.

1.2 Increasing the dimension

We extend the example from the previous section to the following $d$-dimensional problem:

We take $N$ $=$ $120$ in the Euler scheme, and $d+10$ neurons for each hidden layer. We take $1000$ trajectories in mini batch, use data renormalization, and check the loss convergence every $50$ iterations. For this small maturity, the scheme [HJE18] generally converges, and we give the results obtained with the same network and initializing the scheme with the linear solution of the problem. Results in dimension 5 to 50 are given in Tables 2, 3, 4 and 5. Both schemes (3.8) and (3.9) work well with results very close to the solution and close to the results calculated by the scheme [HJE18]. As the dimension increases, scheme (3.8) seems to be the most accurate.

In dimension $50$, the initial learning rate in scheme [HJE18] is taken small in order to avoid a divergence of the method. In fact, running the test 3 times (with 10 runs each time), we observed convergence of the algorithm two times, and in the last test: one of the ten run exploded, and another one clearly converged to a wrong solution. $\Box$

2 PDEs with unbounded solution and more complex structure

In this section with take the following parameters

where the function $k$ is chosen such that the solution to the PDE is equal to

Notice that the structure of the solution is more complex than in the first example. We aim at evaluating the solution at $x=0.51{\rm I}_{d}$. We take $120$ time steps for the Euler time discretization and $d+10$ neurons in each hidden layers. As shown in Figures 3 and 4 as well as in Table 6, the three schemes provide accurate and stable results in dimension $d$ $=$ $1$.

In dimension 2, the three schemes provide very accurate and stable results, as shown in Figures 5 and 6, as well as in Table 7.

Above dimension 3, the scheme [HJE18] always explodes no matter the chosen initial learning rate and the activation function for the hidden layers (among the $\tanh$, ELU, ReLu and sigmoid ones). Besides, taking $3$ or $4$ hidden layers does not improve the results.

We reported the results obtained in dimension $d$ $=$ $5$ and $8$ in Table 8 and 9. Scheme (3.8) seems to work better than scheme (3.9) as the dimension increases. Note that the standard deviation increases with the dimension of the problem.

When $d\geq 10$, schemes (3.8) and (3.9) both fail at providing correct estimates of the solution, as shown in Table 10. Increasing the number of layers or neurons does not improve the result.

3 Application to American options

Consider the stock price $X_{t}$ $=$ $(X^{1}_{t},\dots,X^{d}_{t})$ of $d$ assets with the following dynamics under the risk neutral probability measure:

The value at time $t$ of an American option with payoff $g$ and maturity $T$ is given by:

where ${\cal T}_{t,T}$ is the set of stopping time with values in $[t,T]$, and is solution of the variational inequality

Let us define the change of function $v$ by: $u(t,x)$ $=$ $e^{rt}v(t,\log(x))$, (where $\log$ is applied component-wise), which is solution of the following variational inequality

In this section, we test the scheme described in section 3.3 on (5.7) in the special case of a geometrical put with strike $K=1$ , $T=1$, $r=0.05$, $X_{0}^{i}=1$, $\sigma_{i}=0.2$ for $i=1$ to $d$, and payoff $(K-\prod_{i=1}^{d}X^{i}_{t})_{+}$, as considered previously in [BW12]. In dimension $d$, the case boils down to the resolution of an American option in dimension $d$ $=$ $1$: indeed, the option payoff involving only the product of the asset values, it can be written as the payoff of a single asset with a trend equal to $dr$ and a volatility $\sigma_{1}\sqrt{d}$, so that it can be very accurately estimated e.g. with a tree-based method. Results given in Table 11 show that scheme (3.14) is very accurate for the pricing of American options.

References