Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations

Introduction

Developing efficient numerical algorithms for high dimensional (say, hundreds of dimensions) partial differential equations (PDEs) has been one of the most challenging tasks in applied mathematics. As is well-known, the difficulty lies in the “curse of dimensionality” , namely, as the dimensionality grows, the complexity of the algorithms grows exponentially. For this reason, there are only a limited number of cases where practical high dimensional algorithms have been developed. For linear parabolic PDEs, one can use the Feynman-Kac formula and Monte Carlo methods to develop efficient algorithms to evaluate solutions at any given space-time locations. For a class of inviscid Hamilton-Jacobi equations, Darbon & Osher have recently developed an algorithm which performs numerically well in the case of such high dimensional inviscid Hamilton-Jacobi equations; see . Darbon & Osher’s algorithm is based on results from compressed sensing and on the Hopf formulas for the Hamilton-Jacobi equations. A general algorithm for (nonlinear) parabolic PDEs based on the Feynman-Kac and Bismut-Elworthy-Li formula and a multi-level decomposition of Picard iteration was developed in and has been shown to be quite efficient on a number examples in finance and physics. The complexity of the algorithm is shown to be $O(d\varepsilon^{-4})$ for semilinear heat equations, where $d$ is the dimensionality of the problem and $\varepsilon$ is the required accuracy.

In recent years, a new class of techniques, called deep learning, have emerged in machine learning and have proven to be very effective in dealing with a large class of high dimensional problems in computer vision (cf., e.g., ), natural language processing (cf., e.g., ), time series analysis, etc. (cf., e.g., ). This success fuels in speculations that deep learning might hold the key to solve the curse of dimensionality problem. It should be emphasized that at the present time, there are no theoretical results that support such claims although the practical success of deep learning has been astonishing. However, this should not prevent us from trying to apply deep learning to other problems where the curse of dimensionality has been the issue.

In this paper, we explore the use of deep learning for solving general high dimensional PDEs. To this end, it is necessary to formulate the PDEs as a learning problem. Motivated by ideas in where deep learning-based algorithms were developed for high dimensional stochastic control problems, we explore a connection between (nonlinear) parabolic PDEs and backward stochastic differential equations (BSDEs) (see ) since BSDEs share a lot of common features with stochastic control problems.

Main ideas of the algorithm

We will consider a fairly general class of nonlinear parabolic PDEs (see (30) in Subsection 4.1 below). The proposed algorithm is based on the following set of ideas:

Through the so-called nonlinear Feynman-Kac formula, we can formulate the PDEs equivalently as BSDEs.

One can view the BSDE as a stochastic control problem with the gradient of the solution being the policy function. These stochastic control problems can then be viewed as model-based reinforcement learning problems.

The (high dimensional) policy function can then be approximated by a deep neural network, as has been done in deep reinforcement learning.

Instead of formulating initial value problems, as is commonly done in the PDE literature, we consider the set up with terminal conditions since this facilitates making connections with BSDEs. Terminal value problems can obviously be transformed to initial value problems and vice versa.

In the remainder of this section we present a rough sketch of the derivation of the proposed algorithm, which we refer to as deep BSDE solver. In this derivation we restrict ourself to a specific class of nonlinear PDEs, that is, we restrict ourself to semilinear heat equations (see (PDE) below) and refer to Subsections 3.2 and 4.1 below for the general introduction of the deep BSDE solver.

A key idea of this work is to reformulate the PDE (PDE) as an appropriate stochastic control problem.

2 Formulation of the PDE as a suitable stochastic control problem

One can also view the stochastic control problem (1)–(2) (with $Z$ being the control) as a model-based reinforcement learning problem. In that analogy, we view $Z$ as the policy and we approximate $Z\in\mathcal{A}$ using feedforward neural networks (see (11) and Section 4 below for further details). The process $u(t,\xi+W_{t})$, $t\in[0,T]$, corresponds to the value function associated to the stochastic control problem and can be computed approximatively by employing the policy $Z$ (see (9) below for details). The connection between the PDE (PDE) and the stochastic control problem (1)–(2) is based on the nonlinear Feynman-Kac formula which links PDEs and BSDEs (see (BSDE) and (3) below).

3 The nonlinear Feynman-Kac formula

(cf., e.g., [25, Section 3] and ). The first identity in (3) is sometimes referred to as nonlinear Feynman-Kac formula in the literature.

4 Forward discretization of the backward stochastic differential equation (BSDE)

5 Deep learning-based approximations

In the next step we employ a deep learning approximation for

6 Stochastic optimization algorithms

is minimal. Minimizing the function (12) is inspired by the fact that

of $(\nabla_{x}u)(t_{n},x)$ (cf. (11) above). In the next section the proposed approximation method is described in more detail.

To simplify the presentation we have restricted us in (PDE), (1), (2), (BSDE) above and Subsection 3.1 below to semilinear heat equations. We refer to Subsection 3.2 and Section 4 below for the general description of the deep BSDE solver.

Details of the algorithm

In this subsection we describe the algorithm proposed in this article in the specific situation where (PDE) is the PDE under consideration, where batch normalization (see Ioffe & Szegedy ) is not employed, and where the plain-vanilla stochastic gradient descent approximation method with a constant learning rate $\gamma\in(0,\infty)$ and without mini-batches is the employed stochastic algorithm. The general framework, which includes the setting in this subsection as a special case, can be found in Subsection 3.2 below.

2 Formulation of the proposed algorithm in the general case

3 Comments on the proposed algorithm

stochastic gradient descent with or without mini-batches (see Subsection 5.1 below) as well as

adaptive moment estimation (Adam) with mini-batches (see Kingma & Jimmy and Subsection 5.2 below) into the deep BSDE solver in Subsection 3.2.

In this section we illustrate the algorithm proposed in Subsection 3.2 using several concrete example PDEs. In the examples below we will employ the general approximation method in Subsection 3.2 in conjunction with the Adam optimizer (cf. Example 5.2 below and Kingma & Ba ) with mini-batches with $64$ samples in each iteration step (see Subsection 4.1 for a detailed description).

(cf. Example 5.2 below and Kingma & Ba ).

2 Allen-Cahn equation

In this section we test the deep BSDE solver in the case of an $100$-dimensional Allen-Cahn PDE with a cubic nonlinearity (see (35) below).

In Table 1 we approximatively calculate the mean of $\mathcal{U}^{\Theta_{m}}$, the standard deviation of $\mathcal{U}^{\Theta_{m}}$, the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$, the standard deviation of the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$, and the runtime in seconds needed to calculate one realization of $\mathcal{U}^{\Theta_{m}}$ against $m\in\{0,1000,2000,3000,4000\}$ based on $5$ independent realizations ($5$ independent runs) (see also the Python code 1 below). Table 1 also depicts the mean of the loss function associated to $\Theta_{m}$ and the standard deviation of the loss function associated to $\Theta_{m}$ against $m\in\{0,1000,2000,3000,4000\}$ based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs). In addition, the relative $L^{1}$-approximation error associated to $\mathcal{U}^{\Theta_{m}}$ against $m\in\{1,2,3,\dots,4000\}$ is pictured on the left hand side of Figure 2 based on $5$ independent realizations ($5$ independent runs) and the mean of the loss function associated to $\Theta_{m}$ against $m\in\{1,2,3,\dots,4000\}$ is pictured on the right hand side of Figure 2 based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs). In the approximative computations of the relative $L^{1}$-approximation errors in Table 1 and Figure 2 the value $u(0,\xi)=u(0,0,\dots,0)$ of the exact solution $u$ of the PDE (35) is replaced by the value $0.052802$ which, in turn, is calculated by means of the Branching diffusion method (see the Matlab code 2 below and see, e.g., for analytical and numerical results for the Branching diffusion method in the literature).

3 A Hamilton-Jacobi-Bellman (HJB) equation

In this subsection we apply the deep BSDE solver in Subsection 3.2 to a Hamilton-Jacobi-Bellman (HJB) equation which admits an explicit solution that can be obtained through the Cole-Hopf transformation (cf., e.g., Chassagneux & Richou [7, Section 4.2] and Debnath [10, Section 8.4]).

establish Item (iii). The proof of Lemma 4.2 is thus completed. ∎

4 Pricing of European financial derivatives with different interest rates for borrowing and lending

In this subsection we apply the deep BSDE solver to a pricing problem of an European financial derivative in a financial market where the risk free bank account used for the hedging of the financial derivative has different interest rates for borrowing and lending (see Bergman and, e.g., where this example has been used as a test example for numerical methods for BSDEs).

In Table 3 we approximatively calculate the mean of $\mathcal{U}^{\Theta_{m}}$, the standard deviation of $\mathcal{U}^{\Theta_{m}}$, the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$, the standard deviation of the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$, and the runtime in seconds needed to calculate one realization of $\mathcal{U}^{\Theta_{m}}$ against $m\in\{0,1000,2000,3000,4000\}$ based on $5$ independent realizations ($5$ independent runs). Table 3 also depicts the mean of the loss function associated to $\Theta_{m}$ and the standard deviation of the loss function associated to $\Theta_{m}$ against $m\in\{0,1000,2000,3000,4000\}$ based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs). In addition, the relative $L^{1}$-approximation error associated to $\mathcal{U}^{\Theta_{m}}$ against $m\in\{1,2,3,\dots,4000\}$ is pictured on the left hand side of Figure 4 based on $5$ independent realizations ($5$ independent runs) and the mean of the loss function associated to $\Theta_{m}$ against $m\in\{1,2,3,\dots,4000\}$ is pictured on the right hand side of Figure 4 based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs). In the approximative computations of the relative $L^{1}$-approximation errors in Table 3 and Figure 4 the value $u(0,\xi)=u(0,0,\dots,0)$ of the exact solution $u$ of the PDE (56) is replaced by the value $21.299$ which, in turn, is calculated by means of the multilevel-Picard approximation method in E et al. (see [11, $\rho=7$ in Table 6 in Section 4.3]).

5 Multidimensional Burgers-type PDEs with explicit solutions

In this subsection we consider a high-dimensional version of the example analyzed numerically in Chassagneux [6, Example 4.6 in Subsection 4.2].

(cf. Lemma 4.3 below [with $\alpha=d^{2}$, $\kappa=\nicefrac{{1}}{{d}}$ in the notation of Lemma 4.3 below]). On the left hand side of Figure 5 we present approximatively the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$ against $m\in\{1,2,3,\dots,60\,000\}$ based on $5$ independent realizations ($5$ independent runs) in the case

On the right hand side of Figure 5 we present approximatively the mean of the loss function associated to $\Theta_{m}$ against $m\in\{1,2,3,\dots,60\,000\}$ based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs) in the case (59). On the left hand side of Figure 6 we present approximatively the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$ against $m\in\{1,2,3,\dots,30\,000\}$ based on $5$ independent realizations ($5$ independent runs) in the case

On the right hand side of Figure 6 we present approximatively the mean of the loss function associated to $\Theta_{m}$ against $m\in\{1,2,3,\dots,30\,000\}$ based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs) in the case (60).

Throughout this proof let $\beta,\gamma\in(0,\infty)$ be the real numbers given by

The proof of Lemma 4.3 is thus completed. ∎

6 An example PDE with quadratically growing derivatives and an explicit solution

In this subsection we consider a high-dimensional version of the example analyzed numerically in Gobet & Turkedjiev [13, Section 5]. More specifically, Gobet & Turkedjiev [13, Section 5] employ the PDE in (76) below as a numerical test example but with the time horizont $T=\nicefrac{{2}}{{10}}$ instead of $T=1$ in this article and with the dimension $d\in\{3,5,7\}$ instead of $d=100$ in this article.

On the left hand side of Figure 7 we present approximatively the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$ against $m\in\{1,2,3,\dots,4000\}$ based on $5$ independent realizations ($5$ independent runs). On the right hand side of Figure 7 we present approximatively the mean of the loss function associated to $\Theta_{m}$ against $m\in\{1,2,3,\dots,4000\}$ based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs).

7 Time-dependent reaction-diffusion-type example PDEs with oscillating explicit solutions

In this subsection we consider a high-dimensional version of the example PDE analyzed numerically in Gobet & Turkedjiev [14, Subsection 6.1]. More specifically, Gobet & Turkedjiev [14, Subsection 6.1] employ the PDE in (78) below as a numerical test example but in two space-dimensions ($d=2$) instead of in hundred space-dimensions ($d=100$) as in this article.

(cf. Lemma 4.4 below). On the left hand side of Figure 7 we present approximatively the relative $L^{1}$-approximatin error associated to $\mathcal{U}^{\Theta_{m}}$ against $m\in\{1,2,3,\dots,24000\}$ based on $5$ independent realizations ($5$ independent runs). On the right hand side of Figure 7 we present approximatively the mean of the loss function associated to $\Theta_{m}$ against $m\in\{1,2,3,\dots,24000\}$ based on $256$ Monte Carlo samples and $5$ independent realizations ($5$ independent runs).

The proof of Lemma 4.4 is thus completed. ∎

Appendix A: Special cases of the proposed algorithm

2 Adaptive Moment Estimation (Adam) with mini-batches

In this subsection we illustrate how the so-called Adam optimizer (see ) can be employed in conjunction with the deep BSDE solver in Subsection 3.2 (cf. also Subsection 4.1 above).

3 Geometric Brownian motion

4 Euler-Maruyama scheme

Appendix B: Python and Matlab source codes

2 Matlab source code for the Branching diffusion method used in Subsection 4.2

3 Matlab source code for the classical Monte Carlo method used in Subsection 4.3

Christian Beck and Sebastian Becker are gratefully acknowledged for useful suggestions regarding the implementation of the deep BSDE solver. This project has been partially supported through the Major Program of NNSFC under grant 91130005, the research grant ONR N00014-13-1-0338, and the research grant DOE DE-SC0009248.

References