DeepBSDE
Deep BSDE solver in TensorFlow
view repo
We propose a new algorithm for solving parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) in high dimension, by making an analogy between the BSDE and reinforcement learning with the gradient of the solution playing the role of the policy function, and the loss function given by the error between the prescribed terminal condition and the solution of the BSDE. The policy function is then approximated by a neural network, as is done in deep reinforcement learning. Numerical results using TensorFlow illustrate the efficiency and accuracy of the proposed algorithms for several 100dimensional nonlinear PDEs from physics and finance such as the AllenCahn equation, the HamiltonJacobiBellman equation, and a nonlinear pricing model for financial derivatives.
READ FULL TEXT VIEW PDFDeep BSDE solver in TensorFlow
Deep Learningbased Numerical Methods for PDEs
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 wellknown, the difficulty lies in the “curse of dimensionality”
[1], 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 FeynmanKac formula and Monte Carlo methods to develop efficient algorithms to evaluate solutions at any given spacetime locations. For a class of inviscid HamiltonJacobi equations, Darbon & Osher have recently developed an algorithm which performs numerically well in the case of such high dimensional inviscid HamiltonJacobi equations; see [9]. Darbon & Osher’s algorithm is based on results from compressed sensing and on the Hopf formulas for the HamiltonJacobi equations. A general algorithm for (nonlinear) parabolic PDEs based on the FeynmanKac and BismutElworthyLi formula and a multilevel decomposition of Picard iteration was developed in [11] 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 for semilinear heat equations, where is the dimensionality of the problem and 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.,
[23]), natural language processing (cf., e.g.,
[20]), time series analysis, etc. (cf., e.g., [15, 24]). 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 [16] where deep learningbased algorithms were developed for high dimensional stochastic control problems, we explore a connection between (nonlinear) parabolic PDEs and backward stochastic differential equations (BSDEs) (see [26, 28, 25]) since BSDEs share a lot of common features with stochastic control problems.
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 socalled nonlinear FeynmanKac 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 modelbased 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.
Let , , , let and be continuous functions, and let satisfy for all , that and
(PDE) 
A key idea of this work is to reformulate the PDE (PDE) as an appropriate stochastic control problem.
More specifically, let
be a probability space, let
be a dimensional standard Brownian motion on , let be the normal filtration on generated by , let be the set of all adapted valued stochastic processes with continuous sample paths, and for every and every let be an adapted stochastic process with continuous sample paths which satisfies that for all it holds a.s. that(1) 
We now view the solution of (PDE) and its spatial derivative as the solution of a stochastic control problem associated to (1). More formally, under suitable regularity hypotheses on the nonlinearity it holds that the pair consisting of and is the (up to indistinguishability) unique global minimum of the function
(2) 
One can also view the stochastic control problem (1)–(2) (with being the control) as a modelbased reinforcement learning problem. In that analogy, we view as the policy and we approximate using feedforward neural networks (see (11) and Section 4 below for further details). The process , , corresponds to the value function associated to the stochastic control problem and can be computed approximatively by employing the policy (see (9) below for details). The connection between the PDE (PDE) and the stochastic control problem (1)–(2) is based on the nonlinear FeynmanKac formula which links PDEs and BSDEs (see (BSDE) and (3) below).
Let and be adapted stochastic processes with continuous sample paths which satisfy that for all it holds a.s. that
(BSDE) 
Under suitable additional regularity assumptions on the nonlinearity we have that the nonlinear parabolic PDE (PDE) is related to the BSDE (BSDE) in the sense that for all it holds a.s. that
(3) 
(cf., e.g., [25, Section 3] and [27]). The first identity in (3) is sometimes referred to as nonlinear FeynmanKac formula in the literature.
To derive the deep BSDE solver, we first plug the second identity in (3) into (BSDE) to obtain that for all it holds a.s. that
(4) 
In particular, we obtain that for all with it holds a.s. that
(5) 
Next we apply a time discretization to (5). More specifically, let and let be real numbers which satisfy
(6) 
and observe that (5) suggests for sufficiently large that
(7)  
In the next step we employ a deep learning approximation for
(8) 
but not for , , . Approximations for , , , in turn, can be computed recursively by using (7) together with deep learning approximations for (8). More specifically, let , let , , be real numbers, let , , , be continuous functions, and let , , be stochastic processes which satisfy for all , that and
(9) 
We think of as the number of parameters in the neural network, for all appropriate we think of as suitable approximations
(10) 
of , and for all appropriate , , we think of as suitable approximations
(11) 
of .
The “appropriate”
can be obtained by minimizing the expected loss function through stochastic gradient descenttype algorithms. For the loss function we pick the squared approximation error associated to the terminal condition of the BSDE (
BSDE). More precisely, assume that the function has a unique global minimum and letbe the real vector for which the function
(12) 
is minimal. Minimizing the function (12) is inspired by the fact that
(13) 
according to (BSDE) above (cf. (2) above). Under suitable regularity assumptions, we approximate the vector through stochastic gradient descenttype approximation methods and thereby we obtain random approximations of . For sufficiently large
we then employ the random variable
as a suitable implementable approximation(14) 
of (cf. (10) above) and for sufficiently large and all , we use the random variable as a suitable implementable approximation
(15) 
of (cf. (11) above). In the next section the proposed approximation method is described in more detail.
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 [21]) is not employed, and where the plainvanilla stochastic gradient descent approximation method with a constant learning rate and without minibatches 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.
Let , , , let and be functions, let be a probability space, let , , be independent dimensional standard Brownian motions on , let be real numbers with
(16) 
for every let , for every , let be a function, for every , let be the stochastic process which satisfies for all that and
(17) 
for every let be the function which satisfies for all , that
(18) 
for every let be a function which satisfies for all , that
(19) 
and let be a stochastic process which satisfy for all that
(20) 
Let , , , let , , and be functions, let be a probability space, let , , be independent dimensional standard Brownian motions on , let be real numbers with
(23) 
for every let , for every , , , let be a function, for every let and , , , be stochastic processes which satisfy for all , , that
(24) 
(25) 
for every , let be the function which satisfies for all , that
(26) 
for every , let be a function which satisfies for all , that
(27) 
let be a function, for every let and be functions, and let , , and be stochastic processes which satisfy for all that
(28) 
(29) 
The dynamics in (24) associated to the stochastic processes for allows us to incorporate different algorithms for the discretization of the considered forward stochastic differential equation (SDE) into the deep BSDE solver in Subsection 3.2. The dynamics in (29) associated to the stochastic processes , , and , , allows us to incorporate different stochastic approximation algorithms such as
stochastic gradient descent with or without minibatches (see Subsection 5.1 below) as well as
adaptive moment estimation (Adam) with minibatches (see Kingma & Jimmy
[22] and Subsection 5.2 below) into the deep BSDE solver in Subsection 3.2.The dynamics in (28) associated to the stochastic process , , allows us to incorporate the standardization procedure in batch normalization (see Ioffe & Szegedy [21] and also Section 4 below) into the deep BSDE solver in Subsection 3.2. In that case we think of ,
, as approximatively calculated means and standard deviations.
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 [22]) with minibatches with samples in each iteration step (see Subsection 4.1 for a detailed description).
In our implementation we employ fullyconnected feedforward neural networks to represent for , , (cf. also Figure 1 below for a rough sketch of the architecture of the deep BSDE solver). Each of the neural networks consists of layers ( input layer [dimensional], hidden layers [both dimensional], and output layer [dimensional]). The number of hidden units in each hidden layer is equal to . We also adopt batch normalization (BN) (see Ioffe & Szegedy [21]) right after each matrix multiplication and before activation. We employ the rectifier function
as our activation function for the hidden variables. All the weights in the network are initialized using a normal or a uniform distribution without any pretraining. Each of the numerical experiments presented below is performed in
Python using TensorFlow on a Macbook Pro with a Gigahertz (GHz) Intel Core i5 micro processor and 16 gigabytes (GB) of 1867 Megahertz (MHz) double data rate type three synchronous dynamic randomaccess memory (DDR3SDRAM). We also refer to the Python code 1 in Subsection 6.1 below for an implementation of the deep BSDE solver in the case of the dimensional AllenCahn PDE (35).Assume the setting in Subsection 3.2, assume for all that , , , , , let and be functions, let be a continuous and at most polynomially growing function which satisfies for all that , , and
(30) 
let , , , , , let , , be the functions which satisfy for all , that
(31) 
and assume for all , , that
(32) 
and
(33) 
In this remark we illustrate the specific choice of the dimension of in the framework in Subsection 4.1 above.
The first component of is employed for approximating the real number .
The next components of are employed for approximating the components of the matrix .
In each of the employed neural networks we use components of
to describe the linear transformation from the
dimensional first layer (input layer) to the dimensional second layer (first hidden layer) (to uniquely describe a real matrix).In each of the employed neural networks we use components of to uniquely describe the linear transformation from the dimensional second layer (first hidden layer) to the dimensional third layer (second hidden layer) (to uniquely describe a real matrix).
In each of the employed neural networks we use components of to describe the linear transformation from the dimensional third layer (second hidden layer) to the dimensional fourth layer (output layer) (to uniquely describe a real matrix).
After each of the linear transformations in items (iii)–(v) above we employ a componentwise affine linear transformation (multiplication with a diagonal matrix and addition of a vector) within the batch normalization procedure, i.e., in each of the employed neural networks, we use components of for the componentwise affine linear transformation between the first linear transformation (see item (iii)) and the first application of the activation function, we use components of for the componentwise affine linear transformation between the second linear transformation (see item (iv)) and the second application of the activation function, and we use components of for the componentwise affine linear transformation after the third linear transformation (see item (v)).
(34) 
In this section we test the deep BSDE solver in the case of an dimensional AllenCahn PDE with a cubic nonlinearity (see (35) below).
More specifically, assume the setting in the Subsection 4.1 and assume for all , , , , that , , , , , , , , , and . Note that the solution of the PDE (30) then satisfies for all , that and
(35) 
In Table 1 we approximatively calculate the mean of , the standard deviation of , the relative approximatin error associated to , the standard deviation of the relative approximatin error associated to , and the runtime in seconds needed to calculate one realization of against based on independent realizations ( independent runs) (see also the Python code 1 below). Table 1 also depicts the mean of the loss function associated to and the standard deviation of the loss function associated to against based on Monte Carlo samples and independent realizations ( independent runs). In addition, the relative approximation error associated to against is pictured on the left hand side of Figure 2 based on independent realizations ( independent runs) and the mean of the loss function associated to against is pictured on the right hand side of Figure 2 based on Monte Carlo samples and independent realizations ( independent runs). In the approximative computations of the relative approximation errors in Table 1 and Figure 2 the value of the exact solution of the PDE (35) is replaced by the value which, in turn, is calculated by means of the Branching diffusion method (see the Matlab code 2 below and see, e.g., [17, 19, 18] for analytical and numerical results for the Branching diffusion method in the literature).
Number  Mean  Standard  Relative  Standard  Mean  Standard  Runtime 
of  of  deviation  appr.  deviation  of the  deviation  in sec. 
iteration  of  error  of the  loss  of the  for one  
steps  relative  function  loss  realization  
appr.  function  of  
error  
0  0.4740  0.0514  7.9775  0.9734  0.11630  0.02953 

1000  0.1446  0.0340  1.7384  0.6436  0.00550  0.00344  201 
2000  0.0598  0.0058  0.1318  0.1103  0.00029  0.00006  348 
3000  0.0530  0.0002  0.0050  0.0041  0.00023  0.00001  500 
4000  0.0528  0.0002  0.0030  0.0022  0.00020  0.00001  647 
In this subsection we apply the deep BSDE solver in Subsection 3.2 to a HamiltonJacobiBellman (HJB) equation which admits an explicit solution that can be obtained through the ColeHopf transformation (cf., e.g., Chassagneux & Richou [7, Section 4.2] and Debnath [10, Section 8.4]).
Assume the setting in the Subsection 4.1 and assume for all , , , , that , , , , , , , , , and . Note that the solution of the PDE (30) then satisfies for all , that and
(36) 
In Table 2 we approximatively calculate the mean of , the standard deviation of , the relative approximatin error associated to , the standard deviation of the relative approximatin error associated to , and the runtime in seconds needed to calculate one realization of against based on independent realizations ( independent runs). Table 2 also depicts the mean of the loss function associated to and the standard deviation of the loss function associated to against based on Monte Carlo samples and independent realizations ( independent runs). In addition, the relative approximation error associated to against is pictured on the left hand side of Figure 3 based on independent realizations ( independent runs) and the mean of the loss function associated to against is pictured on the right hand side of Figure 3 based on Monte Carlo samples and independent realizations ( independent runs). In the approximative computations of the relative approximation errors in Table 2 and Figure 3 the value of the exact solution of the PDE (35) is replaced by the value which, in turn, is calculated by means of Lemma 4.2 below (with , , , , in the notation of Lemma 4.2) and a classical Monte Carlo method (see the Matlab code 3 below).
Number  Mean  Standard  Relative  Standard  Mean  Standard  Runtime 
of  of  deviation  appr.  deviation  of the  deviation  in sec. 
iteration  of  error  of the  loss  of the  for one  
steps  relative  function  loss  realization  
appr.  function  of  
error  
0  0.3167  0.3059  0.9310  0.0666  18.4052  2.5090 

500  2.2785  0.3521  0.5036  0.0767  2.1789  0.3848  116 
1000  3.9229  0.3183  0.1454  0.0693  0.5226  0.2859  182 
1500  4.5921  0.0063  0.0013  0.006  0.0239  0.0024  248 
2000  4.5977  0.0019  0.0017  0.0004  0.0231  0.0026  330 