arXiv:1706.04702v1 [math.NA] 15 Jun 2017

Deep learning-based numerical methods for
high-dimensional parabolic partial differential equations
and backward stochastic differential equations
Weinan E1 , Jiequn Han2 , and Arnulf Jentzen3
1 Beijing

Institute of Big Data Research (China), Princeton University (USA),

and Peking University (China), e-mail: weinan (at) math.princeton.edu
2 Princeton
3 ETH

University (USA), e-mail: jiequnh (at) princeton.edu

Zurich (Switzerland), e-mail: arnulf.jentzen (at) sam.math.ethz.ch

June 16, 2017
Abstract
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 100-dimensional nonlinear PDEs
from physics and finance such as the Allen-Cahn equation, the Hamilton-JacobiBellman equation, and a nonlinear pricing model for financial derivatives.

Contents
1 Introduction

2

2 Main ideas of the algorithm
3
2.1 An example: a semilinear heat partial differential equation (PDE) . . . . .
4
2.2 Formulation of the PDE as a suitable stochastic control problem . . . . . .
4
2.3 The nonlinear Feynman-Kac formula . . . . . . . . . . . . . . . . . . . . .
5
2.4 Forward discretization of the backward stochastic differential equation (BSDE) 5
1

2.5
2.6

Deep learning-based approximations . . . . . . . . . . . . . . . . . . . . . .
Stochastic optimization algorithms . . . . . . . . . . . . . . . . . . . . . .

6
6

3 Details of the algorithm
7
3.1 Formulation of the proposed algorithm in the case of semilinear heat equations 7
3.2 Formulation of the proposed algorithm in the general case . . . . . . . . .
8
3.3 Comments on the proposed algorithm . . . . . . . . . . . . . . . . . . . . .
9
4 Examples for nonlinear partial differential equations (PDEs) and nonlinear backward stochastic differential equations (BSDEs)
4.1 Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Allen-Cahn equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 A Hamilton-Jacobi-Bellman (HJB) equation . . . . . . . . . . . . . . . . .
4.4 Pricing of European financial derivatives with different interest rates for
borrowing and lending . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Multidimensional Burgers-type PDEs with explicit solutions . . . . . . . .
4.6 An example PDE with quadratically growing derivatives and an explicit
solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7 Time-dependent reaction-diffusion-type example PDEs with oscillating explicit solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Appendix A: Special cases of the proposed algorithm
5.1 Stochastic gradient descent (SGD) . . . . . . . . . . . . .
5.2 Adaptive Moment Estimation (Adam) with mini-batches
5.3 Geometric Brownian motion . . . . . . . . . . . . . . . .
5.4 Euler-Maruyama scheme . . . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

6 Appendix B: Python and Matlab source codes
6.1 Python source code for an implementation of the deep BSDE solver in the
case of the Allen-Cahn PDE (35) in Subsection 4.2 . . . . . . . . . . . . .
6.2 Matlab source code for the Branching diffusion method used in Subsection 4.2
6.3 Matlab source code for the classical Monte Carlo method used in Subsection 4.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

9
10
12
13
18
20
24
24
27
27
27
28
29
29
29
34
36

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” [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
2

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 HamiltonJacobi equations, Darbon & Osher have recently developed an algorithm which performs
numerically well in the case of such high dimensional inviscid Hamilton-Jacobi equations;
see [9]. 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 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 O(dε−4 ) for semilinear heat equations, where d 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 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 [26, 28, 25]) since BSDEs
share a lot of common features with stochastic control problems.

2

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:
(i) Through the so-called nonlinear Feynman-Kac formula, we can formulate the PDEs
equivalently as BSDEs.
(ii) 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.
(iii) The (high dimensional) policy function can then be approximated by a deep neural
network, as has been done in deep reinforcement learning.
3

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.

2.1

An example: a semilinear heat partial differential equation
(PDE)

Let T ∈ (0, ∞), d ∈ N, ξ ∈ Rd , let f : R × Rd → R and g : Rd → R be continuous
functions, and let u = (u(t, x))t∈[0,T ],x∈Rd ∈ C 1,2 ([0, T ] × Rd , R) satisfy for all t ∈ [0, T ],
x ∈ Rd that u(T, x) = g(x) and

1
∂u
(t, x) + (∆x u)(t, x) + f u(t, x), (∇x u)(t, x) = 0.
∂t
2

(PDE)

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

2.2

Formulation of the PDE as a suitable stochastic control problem

More specifically, let (Ω, F, P) be a probability space, let W : [0, T ] × Ω → Rd be
a d-dimensional standard Brownian motion on (Ω, F, P), let F = (Ft )t∈[0,T ] be the normal filtration on (Ω, F, P) generated by W , let A be the set of all F-adapted Rd -valued
stochastic processes with continuous sample paths, and for every y ∈ R and every Z ∈ A
let Y y,Z : [0, T ] × Ω → R be an F-adapted stochastic process with continuous sample paths
which satisfies that for all t ∈ [0, T ] it holds P-a.s. that
Z t
Z t

y,Z
y,Z
f Ys , Zs ds +
hZs , dWs iRd .
Yt = y −
(1)
0

0

We now view the solution u ∈ C 1,2 ([0, T ]×Rd , R) 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 f it holds that the pair consisting of u(0, ξ) ∈ R
and ((∇x u)(t, ξ + Wt ))t∈[0,T ] ∈ A is the (up to indistinguishability) unique global minimum
of the function


(2)
R × A 3 (y, Z) 7→ E |YTy,Z − g(ξ + WT )|2 ∈ [0, ∞].
4

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 ∈ A using feedforward neural networks (see (11) and Section 4 below
for further details). The process u(t, ξ + Wt ), t ∈ [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).

2.3

The nonlinear Feynman-Kac formula

Let Y : [0, T ] × Ω → R and Z : [0, T ] × Ω → Rd be F-adapted stochastic processes with
continuous sample paths which satisfy that for all t ∈ [0, T ] it holds P-a.s. that
Z T
Z T
f (Ys , Zs ) ds −
hZs , dWs iRd .
(BSDE)
Yt = g(ξ + WT ) +
t

t

Under suitable additional regularity assumptions on the nonlinearity f we have that the
nonlinear parabolic PDE (PDE) is related to the BSDE (BSDE) in the sense that for all
t ∈ [0, T ] it holds P-a.s. that
Yt = u(t, ξ + Wt ) ∈ R

Zt = (∇x u)(t, ξ + Wt ) ∈ Rd

and

(3)

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

2.4

Forward discretization of the backward stochastic differential
equation (BSDE)

To derive the deep BSDE solver, we first plug the second identity in (3) into (BSDE)
to obtain that for all t ∈ [0, T ] it holds P-a.s. that
Z T
Z T

Yt = g(ξ + WT ) +
f Ys , (∇x u)(s, ξ + Ws ) ds −
h(∇x u)(s, ξ + Ws ), dWs iRd . (4)
t

t

In particular, we obtain that for all t1 , t2 ∈ [0, T ] with t1 ≤ t2 it holds P-a.s. that
Z t2
Z t2

Yt2 = Yt1 −
f Ys , (∇x u)(s, ξ + Ws ) ds +
h(∇x u)(s, ξ + Ws ), dWs iRd .
t1

(5)

t1

Next we apply a time discretization to (5). More specifically, let N ∈ N and let t0 , t1 , . . . , tN ∈
[0, T ] be real numbers which satisfy
0 = t0 < t1 < . . . < tN = T
5

(6)

and observe that (5) suggests for N ∈ N sufficiently large that
Ytn+1

(7)


≈ Ytn − f Ytn , (∇x u)(tn , ξ + Wtn ) (tn+1 − tn ) + (∇x u)(tn , ξ + Wtn ), Wtn+1 − Wtn

2.5

Rd

.

Deep learning-based approximations

In the next step we employ a deep learning approximation for
(∇x u)(tn , x) ∈ Rd ,

x ∈ Rd ,

n ∈ {0, 1, . . . , N },

(8)

but not for u(tn , x) ∈ R, x ∈ Rd , n ∈ {0, 1, . . . , N }. Approximations for u(tn , x) ∈ R,
x ∈ Rd , n ∈ {0, 1, . . . , N }, in turn, can be computed recursively by using (7) together with
deep learning approximations for (8). More specifically, let ρ ∈ N, let U θ ∈ R, θ ∈ Rρ , be
real numbers, let Vnθ : Rd → Rd , n ∈ {0, 1, . . . , N − 1}, θ ∈ Rρ , be continuous functions,
and let Y θ : {0, 1, . . . , N } × Ω → R, θ ∈ Rρ , be stochastic processes which satisfy for all
θ ∈ Rρ , n ∈ {0, 1, . . . , N − 1} that Y0θ = U θ and

θ
Yn+1
= Ynθ − f Ynθ , Vnθ (ξ + Wtn ) (tn+1 − tn ) + Vnθ (ξ + Wtn ), Wtn+1 − Wtn Rd . (9)
We think of ρ ∈ N as the number of parameters in the neural network, for all appropriate
θ ∈ Rρ we think of U θ ∈ R as suitable approximations
U θ ≈ u(0, ξ)

(10)

of u(0, ξ), and for all appropriate θ ∈ Rρ , x ∈ Rd , n ∈ {0, 1, . . . , N − 1} we think of
Vnθ (x) ∈ R1×d as suitable approximations
Vnθ (x) ≈ (∇x u)(tn , x)

(11)

of (∇x u)(tn , x).

2.6

Stochastic optimization algorithms

The “appropriate” θ ∈ Rρ can be obtained by minimizing the expected loss function
through stochastic gradient descent-type algorithms. For the loss function we pick the
squared approximation error associated to the terminal
 θ condition2 of
 the BSDE (BSDE).
ρ
More precisely, assume that the function R 3 θ 7→ E |YN −g(XN )| ∈ [0, ∞] has a unique
global minimum and let Λ ∈ Rρ be the real vector for which the function


Rρ 3 θ 7→ E |YNθ − g(XN )|2 ∈ [0, ∞]
(12)
is minimal. Minimizing the function (12) is inspired by the fact that


E |YT − g(XT )|2 = 0
6

(13)

according to (BSDE) above (cf. (2) above). Under suitable regularity assumptions, we
approximate the vector Λ ∈ Rρ through stochastic gradient descent-type approximation
methods and thereby we obtain random approximations Θ0 , Θ1 , Θ2 , . . . : Ω → Rρ of Λ ∈
Rρϑ . For sufficiently large N, ρ, m ∈ N we then employ the random variable U Θm : Ω → R
as a suitable implementable approximation
U Θm ≈ u(0, ξ)

(14)

of u(0, ξ) (cf. (10) above) and for sufficiently large N, ρ, m ∈ N and all x ∈ Rd , n ∈
{0, 1, . . . , N − 1} we use the random variable VnΘm (x) : Ω → R1×d as a suitable implementable approximation
(15)
VnΘm (x) ≈ (∇x u)(tn , x)
of (∇x u)(tn , 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.

3
3.1

Details of the algorithm
Formulation of the proposed algorithm in the case of semilinear heat equations

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 plain-vanilla stochastic gradient
descent approximation method with a constant learning rate γ ∈ (0, ∞) 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.
Framework 3.1 (Specific case). Let T, γ ∈ (0, ∞), d, ρ, N ∈ N, ξ ∈ Rd , let f : R×Rd → R
and g : Rd → R be functions, let (Ω, F, P) be a probability space, let W m : [0, T ] × Ω →
Rd , m ∈ N0 , be independent d-dimensional standard Brownian motions on (Ω, F, P), let
t0 , t1 , . . . , tN ∈ [0, T ] be real numbers with
0 = t0 < t1 < . . . < tN = T,

(16)

for every θ ∈ Rρ let U θ ∈ R, for every θ ∈ Rρ , n ∈ {0, 1, . . . , N − 1} let Vnθ : Rd → Rd be
a function, for every m ∈ N0 , θ ∈ Rρ let Y θ,m : {0, 1, . . . , N } × Ω → Rk be the stochastic
process which satisfies for all n ∈ {0, 1, . . . , N − 1} that Y0θ,m = U θ and

θ,m
Yn+1
= Ynθ,m − f Ynθ,m , Vnθ (ξ + Wtm
) (tn+1 − tn ) + Vnθ (ξ + Wtm
), Wtm
− Wtm
, (17)
n
n
n Rd
n+1
7

for every m ∈ N0 let φm : Rρ × Ω → R be the function which satisfies for all θ ∈ Rρ , ω ∈ Ω
that
2

φm (θ, ω) = YNθ,m (ω) − g(ξ + WTm (ω)) ,

(18)

for every m ∈ N0 let Φm : Rρ × Ω → Rρ be a function which satisfies for all ω ∈ Ω,
ρ
θ ∈ {v ∈ Rρ : (Rρ 3 w 7→ φm
s (w, ω) ∈ R is differentiable at v ∈ R )} that
Φm (θ, ω) = (∇θ φm )(θ, ω),

(19)

and let Θ : N0 × Ω → Rρ be a stochastic process which satisfy for all m ∈ N that
Θm = Θm−1 − γ · Φm (Θm−1 ).

(20)

Under suitable further hypotheses (cf. Sections 4 and 5 below), we think in the case
of sufficiently large ρ, N, m ∈ N and sufficiently small γ ∈ (0, ∞) of U Θm ∈ R as an
appropriate approximation
(21)
u(0, ξ) ≈ U Θm
of the solution u(t, x) ∈ R, (t, x) ∈ [0, T ] × Rd , of the PDE

∂u
1
(t, x) + (∆x u)(t, x) + f u(t, x), (∇x u)(t, x) = 0
∂t
2

(22)

for (t, x) ∈ [0, T ] × Rd .

3.2

Formulation of the proposed algorithm in the general case

Framework 3.2 (General case). Let T ∈ (0, ∞), d, k, ρ, %, N, ς ∈ N, ξ ∈ Rd , let f : [0, T ] ×
Rd ×Rk ×Rk×d → R, g : Rd → Rk , and Υ : [0, T ]2 ×Rd ×Rd → Rd be functions, let (Ω, F, P)
be a probability space, let W m,j : [0, T ] × Ω → Rd , m, j ∈ N0 , be independent d-dimensional
standard Brownian motions on (Ω, F, P), let t0 , t1 , . . . , tN ∈ [0, T ] be real numbers with
0 = t0 < t1 < . . . < tN = T,

(23)

for every θ ∈ Rρ let U θ ∈ Rk , for every θ ∈ Rρ , s ∈ Rς , n ∈ {0, 1, . . . , N − 1}, j ∈ N0 let
θ,s
Vn,j
: (Rd )N → Rk×d be a function, for every m, j ∈ N0 let X m,j : {0, 1, . . . , N } × Ω → Rd
and Y θ,s,m,j : {0, 1, . . . , N } × Ω → Rk , θ ∈ Rρ , s ∈ Rς , be stochastic processes which satisfy
for all θ ∈ Rρ , s ∈ Rς , n ∈ {0, 1, . . . , N − 1} that

m,j
Y0θ,s,m,j = U θ ,
Xn+1
= Υ tn , tn+1 , Xnm,j , Wtm,j
− Wtm,j
,
(24)
X0m,j = ξ,
n
n+1

θ,s
θ,s,m,j
Yn+1
= Ynθ,s,m,j − f tn , Xnm,j , Ynθ,s,m,j , Vn,j
({Xnm,i }i∈N ) (tn+1 − tn )
θ,s
+ Vn,j
({Xnm,i }i∈N ) (Wtm,j
− Wtm,j
), (25)
n
n+1

8

ρ
for every m, j ∈ N0 , s ∈ Rς let φm,j
s : R × Ω → R be the function which satisfies for all
ρ
θ ∈ R , ω ∈ Ω that
θ,s,m,j
φm,j
(ω) − g(XNm,j (ω))k2Rk ,
s (θ, ω) = kYN

(26)

ρ
ρ
for every m, j ∈ N0 , s ∈ Rς let Φm,j
s : R × Ω → R be a function which satisfies for all
ρ
ρ
m,j
ω ∈ Ω, θ ∈ {v ∈ R : (R 3 w 7→ φs (w, ω) ∈ R is differentiable at v ∈ Rρ )} that
m,j
Φm,j
s (θ, ω) = (∇θ φs )(θ, ω),

(27)

let S : Rς × Rρ × (Rd ){0,1,...,N −1}×N → Rς be a function, for every m ∈ N let ψm : R% → Rρ
and Ψm : R% × (Rρ )N → R% be functions, and let S : N0 × Ω → Rς , Ξ : N0 × Ω → R% , and
Θ : N0 × Ω → Rρ be stochastic processes which satisfy for all m ∈ N that

Sm = S Sm−1 , Θm−1 , {Xnm−1,i }(n,i)∈{0,1,...,N −1}×N ,
(28)

(Θm−1 )}j∈N ,
and
Θm = Θm−1 − ψm (Ξm ).
(29)
Ξm = Ψm Ξm−1 , {Φm−1,j
Sm

3.3

Comments on the proposed algorithm

The dynamics in (24) associated to the stochastic processes (Xnm,j )n∈{0,1,...,N } for m, j ∈
N0 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 Ξm , m ∈ N0 , and Θm , m ∈ N0 ,
allows us to incorporate different stochastic approximation algorithms such as
• 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 [22]
and Subsection 5.2 below) into the deep BSDE solver in Subsection 3.2.
The dynamics in (28) associated to the stochastic process Sm , m ∈ N0 , 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 Sm , m ∈ N0 , as approximatively calculated means and standard deviations.

4

Examples for nonlinear partial differential equations
(PDEs) and nonlinear backward stochastic differential equations (BSDEs)

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
9

method in Subsection 3.2 in conjunction with the Adam optimizer (cf. Example 5.2 below
and Kingma & Ba [22]) with mini-batches with 64 samples in each iteration step (see
Subsection 4.1 for a detailed description).
In our implementation we employ N − 1 fully-connected feedforward neural networks to
θ
for n ∈ {1, 2, . . . , N − 1}, j ∈ {1, 2, . . . , 64}, θ ∈ Rρ (cf. also Figure 1 below
represent Vn,j
for a rough sketch of the architecture of the deep BSDE solver). Each of the neural networks
consists of 4 layers (1 input layer [d-dimensional], 2 hidden layers [both d+10-dimensional],
and 1 output layer [d-dimensional]). The number of hidden units in each hidden layer is
equal to d + 10. We also adopt batch normalization (BN) (see Ioffe & Szegedy [21]) right
after each matrix multiplication and before activation. We employ the rectifier function
R 3 x 7→ max{0, x} ∈ [0, ∞) 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
pre-training. Each of the numerical experiments presented below is performed in Python
using TensorFlow on a Macbook Pro with a 2.90 Gigahertz (GHz) Intel Core i5
micro processor and 16 gigabytes (GB) of 1867 Megahertz (MHz) double data rate type
three synchronous dynamic random-access memory (DDR3-SDRAM). 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 100-dimensional Allen-Cahn PDE (35).

4.1

Setting

Assume the setting in Subsection 3.2, assume for all θ = (θ1 , . . . , θρ ) ∈ Rρ that k = 1,
ρ = d + 1 + (N − 1) (2d(d + 10) + (d + 10)2 + 4(d + 10) + 2d), % = 2ρ, U θ = θ1 , Ξ0 = 0,
let µ : [0, T ] × Rd → Rd and σ : [0, T ] × Rd → Rd×d be functions, let u : [0, T ] × Rd → R
be a continuous and at most polynomially growing function which satisfies for all (t, x) ∈
[0, T ) × Rd that u|[0,T )×Rd ∈ C 1,2 ([0, T ) × Rd , R), u(T, x) = g(x), and

1
∂u
(t, x) + Trace σ(t, x) [σ(t, x)]∗ (Hessx u)(t, x) + hµ(t, x), (∇x u)(t, x)i
∂t
2

+ f t, x, u(t, x), [(∇x u)(t, x)]∗ σ(t, x) = 0, (30)
999
9
let ε = 10−8 , X = 10
, Y = 1000
, J = 64, (γm )m∈N ⊆ (0, ∞), let Powr : Rρ → Rρ , r ∈ (0, ∞),
be the functions which satisfy for all r ∈ (0, ∞), x = (x1 , . . . , xρ ) ∈ Rρ that

Powr (x) = (|x1 |r , . . . , |xρ |r ),

(31)

and assume for all m ∈ N, x, y ∈ Rρ , (ϕj )j∈N ∈ (Rρ )N that

P
Ψm (x, y, (ϕj )j∈N ) = Xx + (1 − X) J1 Jj=1 ϕj , Yy + (1 − Y) Pow2
and


ψm (x, y) = ε + Pow1/2

y
(1 − Ym )
(cf. Example 5.2 below and Kingma & Ba [22]).
10

−1

γm x
.
(1 − Xm )

1
J

PJ

j=1

ϕj



(32)

(33)

Remark 4.1. In this remark we illustrate the specific choice of the dimension ρ ∈ N of
θ = (θ1 , . . . , θρ ) ∈ Rρ in the framework in Subsection 4.1 above.
(i) The first component of θ = (θ1 , . . . , θρ ) ∈ Rρ is employed for approximating the real
number u(0, ξ) ∈ R.
(ii) The next d-components of θ = (θ1 , . . . , θρ ) ∈ Rρ are employed for approximating the
∂
u)(0, ξ) σ(0, ξ) ∈ R1×d .
components of the (1 × d)-matrix ( ∂x
(iii) In each of the employed N − 1 neural networks we use d(d + 10) components of
θ = (θ1 , . . . , θρ ) ∈ Rρ to describe the linear transformation from the d-dimensional
first layer (input layer) to the (d + 10)-dimensional second layer (first hidden layer)
(to uniquely describe a real (d + 10) × d-matrix).
(iv) In each of the employed N − 1 neural networks we use (d + 10)2 components of
θ = (θ1 , . . . , θρ ) ∈ Rρ to uniquely describe the linear transformation from the (d + 10)dimensional second layer (first hidden layer) to the (d + 10)-dimensional third layer
(second hidden layer) (to uniquely describe a real (d + 10) × (d + 10)-matrix).
(v) In each of the employed N − 1 neural networks we use d(d + 10) components of θ =
(θ1 , . . . , θρ ) ∈ Rρ to describe the linear transformation from the (d + 10)-dimensional
third layer (second hidden layer) to the d-dimensional fourth layer (output layer) (to
uniquely describe a real d × (d + 10)-matrix).
(vi) 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 N −1 neural networks, we use 2(d+10) components of θ = (θ1 , . . . , θρ ) ∈ Rρ for
the componentwise affine linear transformation between the first linear transformation
(see item (iii)) and the first application of the activation function, we use 2(d + 10)
components of θ = (θ1 , . . . , θρ ) ∈ Rρ 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 2d components of θ = (θ1 , . . . , θρ ) ∈ Rρ for
the componentwise affine linear transformation after the third linear transformation
(see item (v)).
Summing (i)–(vi) results in
ρ=

+ d}
|1 {z

items (i)–(ii)


+ (N − 1) d(d + 10) + (d + 10)2 + d(d + 10)
{z
}
|
items (iii)–(v)

+ (N − 1) (2(d + 10) + 2(d + 10) + 2d)
|
{z
}
item (vi)


= d + 1 + (N − 1) 2d(d + 10) + (d + 10)2 + 4(d + 10) + 2d .

11

(34)

Figure 1: Rough sketch of the architecture of the deep BSDE solver.

4.2

Allen-Cahn equation

In this section we test the deep BSDE solver in the case of an 100-dimensional AllenCahn PDE with a cubic nonlinearity (see (35) below).
More specifically, assume the setting in the Subsection 4.1 and assume for all s, t ∈
3
,
[0, T ], x, w ∈ Rd , y ∈ R, z ∈ R1×d , m ∈ N that γm = 5 · 10−4 , d = 100, T = 10
√
√
d
N = 20, µ(t, x) = 0, σ(t, x)w = 2 w, ξ = (0, 0, . . . , 0) ∈ R , Υ(s, t, x, w) = x + 2 w,

−1
f (t, x, y, z) = y − y 3 , and g(x) = 2 + 25 kxk2Rd . Note that the solution u of the PDE
(30) then satisfies for all t ∈ [0, T ), x ∈ Rd that u(T, x) = g(x) and
∂u
(t, x) + u(t, x) − [u(t, x)]3 + (∆x u)(t, x) = 0.
∂t

(35)

In Table 1 we approximatively calculate the mean of U Θm , the standard deviation of U Θm ,
the relative L1 -approximatin error associated to U Θm , the standard deviation of the relative
L1 -approximatin error associated to U Θm , and the runtime in seconds needed to calculate
one realization of U Θm against m ∈ {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 Θm and the standard deviation of the loss function associated to Θm against m ∈ {0, 1000, 2000, 3000, 4000} based on 256 Monte Carlo
samples and 5 independent realizations (5 independent runs). In addition, the relative
L1 -approximation error associated to U Θm against m ∈ {1, 2, 3, . . . , 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 Θm against m ∈ {1, 2, 3, . . . , 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 L1 12

approximation errors in Table 1 and Figure 2 the value u(0, ξ) = u(0, 0, . . . , 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.,
[17, 19, 18] for analytical and numerical results for the Branching diffusion method in the
literature).
Number Mean
of
of U Θm
iteration
steps m

0
1000
2000
3000
4000

0.4740
0.1446
0.0598
0.0530
0.0528

Standard Relative Standard
Mean
Standard Runtime
deviation L1 -appr. deviation
of the
deviation
in sec.
Θm
of U
error
of the
loss
of the
for one
relative function
loss
realization
L1 -appr.
function
of U Θm
error
0.0514
7.9775
0.9734
0.11630
0.02953
0.0340
1.7384
0.6436
0.00550
0.00344
201
0.0058
0.1318
0.1103
0.00029
0.00006
348
0.0002
0.0050
0.0041
0.00023
0.00001
500
0.0002
0.0030
0.0022
0.00020
0.00001
647

Table 1: Numerical simulations for the deep BSDE solver in Subsection 3.2 in the case of
the PDE (35).

4.3

A Hamilton-Jacobi-Bellman (HJB) equation

In this subsection we apply the deep BSDE solver in Subsection 3.2 to a HamiltonJacobi-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]).
Assume the setting in the Subsection 4.1 and assume for all s, t ∈ [0, T ], x, w ∈ Rd ,
1
y ∈ R, z ∈ R1×d , m ∈ N that d = 100, T = 1, N = 20, γm = 100
, µ(t, x) = 0,
√
√
d
σ(t, x)w = 2 w, ξ = (0, 0, . . . , 0) ∈ R , Υ(s, t, x, w) = x + 2 w, f (t, x, y, z) = −kzk2R1×d ,
and g(x) = ln( 12 [1 + kxk2Rd ]). Note that the solution u of the PDE (30) then satisfies for
all t ∈ [0, T ), x ∈ Rd that u(T, x) = g(x) and
∂u
(t, x) + (∆x u)(t, x) = k(∇x u)(t, x)k2Rd .
∂t

(36)

In Table 2 we approximatively calculate the mean of U Θm , the standard deviation of U Θm ,
the relative L1 -approximatin error associated to U Θm , the standard deviation of the relative L1 -approximatin error associated to U Θm , and the runtime in seconds needed to
13

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 2: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 4000} in the case of the PDE (35). The deep BSDE approximation
U Θ4000 ≈ u(0, ξ) achieves a relative L1 -approximation error of size 0.0030 in a runtime of
595 seconds.
calculate one realization of U Θm against m ∈ {0, 500, 1000, 1500, 2000} based on 5 independent realizations (5 independent runs). Table 2 also depicts the mean of the loss function
associated to Θm and the standard deviation of the loss function associated to Θm against
m ∈ {0, 500, 1000, 1500, 2000} based on 256 Monte Carlo samples and 5 independent realizations (5 independent runs). In addition, the relative L1 -approximation error associated
to U Θm against m ∈ {1, 2, 3, . . . , 2000} is pictured on the left hand side of Figure 3 based
on 5 independent realizations (5 independent runs) and the mean of the loss function associated to Θm against m ∈ {1, 2, 3, . . . , 2000} is pictured on the right hand side of Figure 3
based on 256 Monte Carlo samples and 5 independent realizations (5 independent runs).
In the approximative computations of the relative L1 -approximation errors in Table 2 and
Figure 3 the value u(0, ξ) = u(0, 0, . . . , 0) of the exact solution u of the PDE (35) is replaced by the value 4.5901 which, in turn, is calculated by means of Lemma 4.2 below
(with d = 100, T = 1, α = 1, β = −1, g = Rd 3 x 7→ ln( 21 [1 + kxk2Rd ]) ∈ R in the notation
of Lemma 4.2) and a classical Monte Carlo method (see the Matlab code 3 below).
Lemma 4.2 (Cf., e.g., Section 4.2 in [7] and Section 8.4 in [10]). Let d ∈ N, T, α ∈
(0, ∞), β ∈ R\{0}, let (Ω, F, P) be a probability space, let W : [0, T ] × Ω → Rd be a ddimensional standard Brownian motion, let g ∈ C 2 (Rd , R) be a function which satisfies
supx∈Rd [βg(x)] < ∞, let f : [0, T ] × Rd × R × Rd → R be the function which satisfies for
all t ∈ [0, T ], x = (x1 , . . . , xd ), z = (z1 , . . . , zd ) ∈ Rd , y ∈ R that
P
f (t, x, y, z) = βkzk2Rd = β di=1 |zi |2 ,
(37)
14

Number Mean
of
of U Θm
iteration
steps m

0
500
1000
1500
2000

0.3167
2.2785
3.9229
4.5921
4.5977

Standard Relative Standard
Mean
Standard Runtime
1
deviation L -appr. deviation
of the
deviation
in sec.
Θm
of U
error
of the
loss
of the
for one
relative function
loss
realization
L1 -appr.
function
of U Θm
error
0.3059
0.9310
0.0666
18.4052
2.5090
0.3521
0.5036
0.0767
2.1789
0.3848
116
0.3183
0.1454
0.0693
0.5226
0.2859
182
0.0063
0.0013
0.006
0.0239
0.0024
248
0.0019
0.0017
0.0004
0.0231
0.0026
330

Table 2: Numerical simulations for the deep BSDE solver in Subsection 3.2 in the case of
the PDE (36).
and let u : [0, T ] × Rd → R be the function which satisfies for all (t, x) ∈ [0, T ] × Rd that
 h
 βg(x + W √2α) i
α
T −t
u(t, x) = ln E exp
.
(38)
β
α
Then
(i) it holds that u : [0, T ] × Rd → R is a continuous function,
(ii) it holds that u|[0,T )×Rd ∈ C 1,2 ([0, T ) × Rd , R), and
(iii) it holds for all (t, x) ∈ [0, T ) × Rd that u(T, x) = g(x) and

∂u
(t, x) + α(∆x u)(t, x) + f t, x, u(t, x), (∇x u)(t, x)
∂t
∂u
=
(t, x) + α(∆x u)(t, x) + βk(∇x u)(t, x)k2Rd
∂t
d
X
2
∂u
∂u
(t, x) + α(∆x u)(t, x) + β
(t, x) = 0.
=
∂t
∂xj
j=1

(39)

Proof of Lemma 4.2. Throughout this proof let c = αβ ∈ R\{0} and let V : Rd → (0, ∞)
and v : [0, T ] × Rd → (0, ∞) be the functions which satisfy for all t ∈ [0, T ], x ∈ Rd that
V(x) = exp

 g(x) 
c

 βg(x) 
= exp
α

and

15

√ 

v(t, x) = E V x + WT −t 2α .

(40)

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 3: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 2000}. The deep BSDE approximation U Θ2000 ≈ u(0, ξ) achieves a relative
L1 -approximation error of size 0.0017 in a runtime of 283 seconds.
Observe that the hypothesis that supx∈Rd [βg(x)] < ∞ ensures that for all ω ∈ Ω it holds
that
√ 
sup sup V x + WT −t (ω) 2α ≤ sup |V(x)| = sup V(x)
t∈[0,T ] x∈Rd

x∈Rd

x∈Rd



supx∈Rd [βg(x)]
= exp
< ∞.
α

(41)

Combining this with Lebesgue’s theorem of dominated convergence ensures that v : [0, T ]×
Rd → (0, ∞) is a continuous function. This and the fact that
∀ (t, x) ∈ [0, T ] × Rd : u(t, x) = c ln(v(t, x))

(42)

establish Item (i). Next note that the Feynman-Kac formula ensures that for all t ∈ [0, T ),
x ∈ Rd it holds that v|[0,T )×Rd ∈ C 1,2 ([0, T ) × Rd , (0, ∞)) and
∂v
(t, x) + α(∆x v)(t, x) = 0.
∂t

(43)

This and (42) demonstrate Item (ii). It thus remains to prove Item (iii). For this note that
the chain rule and (42) imply that for all t ∈ [0, T ), x = (x1 , . . . , xd ) ∈ Rd , i ∈ {1, 2, . . . , d}
it holds that
∂u
c
∂v
(t, x) =
·
(t, x)
∂t
v(t, x) ∂t

and

16

∂u
c
∂v
(t, x) =
·
(t, x).
∂xi
v(t, x) ∂xi

(44)

Again the chain rule and (42) hence ensure that for all t ∈ [0, T ), x = (x1 , . . . , xd ) ∈ Rd ,
i ∈ {1, 2, . . . , d} it holds that

2
∂u2
c
∂v 2
c
∂v
(t, x) =
·
(t, x) −
·
(t, x) .
(45)
∂x2i
v(t, x) ∂x2i
[v(t, x)]2 ∂xi
This assures that for all t ∈ [0, T ), x = (x1 , . . . , xd ) ∈ Rd , i ∈ {1, 2, . . . , d} it holds that
2
d 
X
αc
αc
∂v
α(∆x u)(t, x) =
· (∆x v)(t, x) −
(t, x)
·
v(t, x)
[v(t, x)]2 i=1 ∂xi
αc(∆x v)(t, x) αc k(∇x v)(t, x)k2Rd
=
.
−
v(t, x)
[v(t, x)]2

(46)

Combining this with (44) demonstrates that for all t ∈ [0, T ), x ∈ Rd it holds that
∂u
(t, x) + α(∆x u)(t, x) + β k(∇x u)(t, x)k2Rd
∂t
c
∂v
αc(∆x v)(t, x) αc k(∇x v)(t, x)k2Rd
=
+ β k(∇x u)(t, x)k2Rd .
·
(t, x) +
−
2
v(t, x) ∂t
v(t, x)
[v(t, x)]

(47)

Equation (43) hence shows that for all t ∈ [0, T ), x = (x1 , . . . , xd ) ∈ Rd it holds that
∂u
(t, x) + α(∆x u)(t, x) + β k(∇x u)(t, x)k2Rd
∂t
αc k(∇x v)(t, x)k2Rd
2
= β k(∇x u)(t, x)kRd −
[v(t, x)]2
" d
#
2
X ∂u
αc k(∇x v)(t, x)k2Rd
=β
(t, x)
.
−
2
∂x
[v(t,
x)]
i
i=1

(48)

This and (44) demonstrate that for all t ∈ [0, T ), x = (x1 , . . . , xd ) ∈ Rd it holds that
∂u
(t, x) + α(∆x u)(t, x) + β k(∇x u)(t, x)k2Rd
∂t "
#
d
2
X
αc k(∇x v)(t, x)k2Rd
c
∂v
=β
·
(t, x)
−
v(t, x) ∂xi
[v(t, x)]2
i=1
" d
#
2
X ∂v
αc k(∇x v)(t, x)k2Rd
c2 β
=
−
(t, x)
[v(t, x)]2 i=1 ∂xi
[v(t, x)]2
=

[c2 β − cα] k(∇x v)(t, x)k2Rd
= 0.
[v(t, x)]2
17

(49)

This and the fact that

 
g(x)
= g(x)
∀ x ∈ R : u(T, x) = c ln(v(T, x)) = c ln(V(x)) = c ln exp
c
d

(50)

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

4.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 [4] and, e.g., [12, 2, 3, 5, 8, 11] where this example has been used as a test
example for numerical methods for BSDEs).
2
4
6
6
, σ̄ = 10
, Rl = 100
, Rb = 100
, and
Assume the setting in Subsection 4.1, let µ̄ = 100
d
assume for all s, t ∈ [0, T ], x = (x1 , . . . , xd ), w = (w1 , . . . , wd ) ∈ R , y ∈ R, z ∈ R1×d ,
m ∈ N that d = 100, T = 1/2, N = 20, γm = 5 · 10−3 = 0.005, µ(t, x) = µ̄x, σ(t, x) =
σ̄ diagRd×d (x1 , . . . , xd ), ξ = (100, 100, . . . , 100) ∈ Rd , and






g(x) = max
max xi − 120, 0 − 2 max
max xi − 150, 0 ,
(51)
1≤i≤100

1≤i≤100


2





σ̄
2


(t − s) exp σ̄ diagRd×d (w1 , . . . , wd ) x,
( " d #
)
d
l X
X
(µ̄
−
R
)
1
f (t, x, y, z) = −Rl y −
zi + (Rb − Rl ) max 0,
zi − y .
σ̄
σ̄
i=1
i=1
Υ(s, t, x, w) = exp

µ̄ −

(52)
(53)

Note that the solution u of the PDE (30) then satisfies for all t ∈ [0, T ), x = (x1 , x2 , . . . , xd ) ∈
Rd that u(T, x) = g(x) and
d

X ∂u

∂u
(t, x) + f t, x, u(t, x), σ̄ diagRd×d (x1 , . . . , xd )(∇x u)(t, x) + µ̄
xi
(t, x)
∂t
∂xi
i=1
d
σ̄ 2 X
∂ 2u
+
|xi |2 2 (t, x) = 0. (54)
2 i=1
∂xi

Hence, we obtain for all t ∈ [0, T ), x = (x1 , x2 , . . . , xd ) ∈ Rd that u(T, x) = g(x) and
d
σ̄ 2 X
∂ 2u
∂u
(t, x) +
|xi |2 2 (t, x)
∂t
2 i=1
∂xi
n hP
i

hP

d
d
l
∂u
+ max Rb
x
(t,
x)
−
u(t,
x)
,
R
i=1 i ∂xi
i=1 xi

∂u
∂xi



i
o
(t, x) − u(t, x)
= 0.
(55)

18

This shows that for all t ∈ [0, T ), x = (x1 , x2 , . . . , xd ) ∈ Rd it holds that u(T, x) = g(x)
and
d

∂ 2u
∂u
σ̄ 2 X
(t, x) +
|xi |2 2 (t, x)
∂t
2 i=1
∂xi
 

 
d
d
X
X
∂u
∂u
l
b
xi
xi
(t, x) , R u(t, x) −
(t, x)
= 0. (56)
− min R u(t, x) −
∂xi
∂xi
i=1
i=1
In Table 3 we approximatively calculate the mean of U Θm , the standard deviation of U Θm ,
the relative L1 -approximatin error associated to U Θm , the standard deviation of the relative
L1 -approximatin error associated to U Θm , and the runtime in seconds needed to calculate
one realization of U Θm against m ∈ {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 Θm and the standard deviation of the loss function associated to Θm against
m ∈ {0, 1000, 2000, 3000, 4000} based on 256 Monte Carlo samples and 5 independent realizations (5 independent runs). In addition, the relative L1 -approximation error associated
to U Θm against m ∈ {1, 2, 3, . . . , 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 Θm against m ∈ {1, 2, 3, . . . , 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 L1 -approximation errors in Table 3 and
Figure 4 the value u(0, ξ) = u(0, 0, . . . , 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. [11] (see [11, ρ = 7 in Table 6 in Section 4.3]).
Number Mean
of
of U Θm
iteration
steps m

0
1000
2000
3000
4000

16.964
20.309
21.150
21.229
21.217

Standard Relative Standard
Mean
Standard Runtime
deviation L1 -appr. deviation
of the
deviation
in sec.
Θm
of U
error
of the
loss
of the
for one
relative function
loss
realization
L1 -appr.
function
of U Θm
error
0.882
0.204
0.041
53.666
8.957
0.524
0.046
0.025
30.886
3.076
194
0.098
0.007
0.005
29.197
3.160
331
0.034
0.003
0.002
29.070
3.246
470
0.043
0.004
0.002
29.029
3.236
617

Table 3: Numerical simulations for the deep BSDE solver in Subsection 3.2 in the case of
the PDE (56).

19

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 4: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 4000} in the case of the PDE (56). The deep BSDE approximation
U Θ4000 ≈ u(0, ξ) achieves a relative L1 -approximation error of size 0.0039 in a runtime of
566 seconds.

4.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].
More specifically, assume the setting in Subsection 4.1, and assume for all s, t ∈ [0, T ],
x = (x1 , . . . , xd ), w = (w1 , . . . , wd ) ∈ Rd , y ∈ R, z = (zi )i∈{1,2,...,d} ∈ R1×d that µ(t, x) = 0,
σ(t, x)w = √d2 w, ξ = (0, 0, . . . , 0) ∈ Rd , Υ(s, t, x, w) = x + √d2 w, and
g(x) =

exp(T +

1
d

1 + exp(T +

Pd

i=1 xi )
,
P
d
1
x
)
i
i=1
d


f (t, x, y, z) =

2+d
y−
2d



d
P


zi .

(57)

i=1

Note that the solution u of the PDE (30) then satisfies for all t ∈ [0, T ), x = (x1 , x2 , . . . , xd ) ∈
Rd that u(T, x) = g(x) and

 d

P ∂u
d2
2+d
∂u
(t, x) + (∆x u)(t, x) + u(t, x) −
d
(t, x) = 0
(58)
∂t
2
2d
i=1 ∂xi
(cf. Lemma 4.3 below [with α = d2 , κ = 1/d in the notation of Lemma 4.3 below]). On the
left hand side of Figure 5 we present approximatively the relative L1 -approximatin error
associated to U Θm against m ∈ {1, 2, 3, . . . , 60 000} based on 5 independent realizations (5
independent runs) in the case
T = 1,

d = 20,

N = 80,

∀ m ∈ N : γm = 10(1[1,30000] (m)+1[1,50000] (m)−4) .
20

(59)

On the right hand side of Figure 5 we present approximatively the mean of the loss function
associated to Θm against m ∈ {1, 2, 3, . . . , 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 L1 -approximatin error associated to U Θm
against m ∈ {1, 2, 3, . . . , 30 000} based on 5 independent realizations (5 independent runs)
in the case
2
d = 50,
N = 30,
∀ m ∈ N : γm = 10(1[1,15000] (m)+1[1,25000] (m)−4) . (60)
T = ,
10
On the right hand side of Figure 6 we present approximatively the mean of the loss function
associated to Θm against m ∈ {1, 2, 3, . . . , 30 000} based on 256 Monte Carlo samples and
5 independent realizations (5 independent runs) in the case (60).

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 5: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 60 000} in the case of the PDE (58) with (59). The deep BSDE approximation U Θ60 000 ≈ u(0, ξ) achieves a relative L1 -approximation error of size 0.0073 in a
runtime of 20 389 seconds.
Lemma 4.3 (Cf. Example 4.6 in Subsection 4.2 in [6]). Let α, κ, T ∈ (0, ∞), d ∈ N, let
u : [0, T ] × Rd → R be the function which satisfies for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd
that
P
exp(t + κ di=1 xi )
1
=
,
(61)
u(t, x) = 1 −
P
P
(1 + exp(t + κ di=1 xi ))
(1 + exp(t + κ di=1 xi ))
and let f : [0, T ] × Rd × R1+d → R be the function which satisfies for all t ∈ [0, T ], x ∈ Rd ,
y ∈ R, z = (z1 , . . . , zd ) ∈ Rd that

 d 
P
1
ακ
f (t, x, y, z) = ακy −
−
zi .
(62)
dκ
2
i=1
21

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 6: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 30 000} in the case of the PDE (58) with (60). The deep BSDE approximation U Θ30 000 ≈ u(0, ξ) achieves a relative L1 -approximation error of size 0.0035 in a
runtime of 4281 seconds.
Then it holds for all t ∈ [0, T ], x ∈ Rd that

∂u
α
(t, x) + (∆x u)(t, x) + f t, x, u(t, x), (∇x u)(t, x) = 0.
∂t
2

(63)

Proof of Lemma 4.3. Throughout this proof let β, γ ∈ (0, ∞) be the real numbers given
by
1
ακ
β = ακ
and
γ=
+
(64)
dκ
2
and let w : [0, T ] × Rd → (0, ∞) be the function which satisfies for all t ∈ [0, T ], x =
(x1 , . . . , xd ) ∈ Rd that
!
d
X
w(t, x) = exp t + κ
xi .
(65)
i=1

Observe that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd , i ∈ {1, 2, . . . , d} it holds that
u(t, x) = 1 − [1 + w(t, x)]−1 =

and

[1 + w(t, x)]
1
w(t, x)
−
=
,
[1 + w(t, x)] [1 + w(t, x)]
1 + w(t, x)

(66)

∂u
∂w
w(t, x)
(t, x) = [1 + w(t, x)]−2 ·
(t, x) =
,
∂t
∂t
[1 + w(t, x)]2

(67)

∂w
∂u
(t, x) = [1 + w(t, x)]−2 ·
(t, x) = κ w(t, x) [1 + w(t, x)]−2 .
∂xi
∂xi

(68)

22

Note that (66), (67), and (68) ensure that for all t ∈ [0, T ], x ∈ Rd it holds that

α
∂u
(t, x) + (∆x u)(t, x) + f t, x, u(t, x), (∇x u)(t, x)
∂t
2
d

P ∂u
∂u
α
=
(t, x) + (∆x u)(t, x) + (βu(t, x) − γ)
(t, x)
∂t
2
i=1 ∂xi

α
∂u
∂u
(t, x) + (∆x u)(t, x) + d
(t, x) βu(t, x) − γ
=
∂t
2
∂x1


w(t, x)
α
dκw(t, x)
βw(t, x)
=
+ (∆x u)(t, x) +
−γ .
2
[1 + w(t, x)]2
[1 + w(t, x)]2 [1 + w(t, x)]

(69)

Moreover, observe that (68) demonstrates that for all t ∈ [0, T ], x ∈ Rd it holds that
∂ 2u
∂w
∂w
(t, x) = κ
(t, x) [1 + w(t, x)]−2 − 2κw(t, x) [1 + w(t, x)]−3
(t, x)
2
∂xi
∂xi
∂xi


κ2 w(t, x)
2κ2 |w(t, x)|2
κ2 w(t, x)
2w(t, x)
=
−
=
1−
.
[1 + w(t, x)]
[1 + w(t, x)]2 [1 + w(t, x)]3
[1 + w(t, x)]2
Hence, we obtain that for all t ∈ [0, T ], x ∈ Rd it holds that


dα ∂ 2 u
dακ2 w(t, x)
α
2w(t, x)
(∆x u)(t, x) =
(t, x) =
.
1−
2
2 ∂x21
[1 + w(t, x)]
2 [1 + w(t, x)]2

(70)

(71)

Combining this with (69) implies that for all t ∈ [0, T ], x ∈ Rd it holds that

α
∂u
(t, x) + (∆x u)(t, x) + f t, x, u(t, x), (∇x u)(t, x)
∂t
2
w(t, x) [1 − dκγ] α
dβκ|w(t, x)|2
+
=
(∆
u)(t,
x)
+
x
2
[1 + w(t, x)]2
[1 + w(t, x)]3


dακ2 w(t, x)
w(t, x) [1 − dκγ]
2w(t, x)
dβκ|w(t, x)|2
=
+
1
−
+
[1 + w(t, x)]
[1 + w(t, x)]2
2 [1 + w(t, x)]2
[1 + w(t, x)]3


2
w(t, x) 1 − dκγ + dακ
dακ2 |w(t, x)|2 dβκ|w(t, x)|2
2
=
−
+
.
[1 + w(t, x)]2
[1 + w(t, x)]3
[1 + w(t, x)]3

(72)

The fact that ακ2 = βκ hence demonstrates that for all t ∈ [0, T ], x ∈ Rd it holds that

∂u
α
(t, x) + (∆x u)(t, x) + f t, x, u(t, x), (∇x u)(t, x)
∂t
2
2
w(t, x) 1 − dκγ + dακ
2
=
.
[1 + w(t, x)]2
23

(73)

This and the fact that 1 +

dακ2
2

= dκγ show that for all t ∈ [0, T ], x ∈ Rd it holds that


∂u
α
(t, x) + (∆x u)(t, x) + f t, x, u(t, x), (∇x u)(t, x) = 0.
∂t
2

(74)

The proof of Lemma 4.3 is thus completed.

4.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 = 2/10 instead of T = 1 in this article and with the dimension d ∈ {3, 5, 7}
instead of d = 100 in this article.
Assume the setting in Subsection 4.1, let α = 4/10, let ψ : [0, T ]×Rd → R be
 the function
d
2 α
which satisfies for all (t, x) ∈ [0, T ] × R that ψ(t, x) = sin [T − t + kxkRd ] , and assume
for all s ∈ [0, T ], t ∈ [0, T ), x, w ∈ Rd , y ∈ R, z ∈ R1×d , m ∈ N that T = 1, d = 100,
5
= 0.005, µ(t, x) = 0, σ(t, x)w = w, ξ = (0, 0, . . . , 0) ∈ Rd ,
N = 30, γm = 5 · 10−3 = 1000
Υ(t, s, x, w) = x + w, g(x) = sin(kxk2α
), and
Rd
f (t, x, y, z) = kzk2R1×d − k(∇x ψ)(t, x)k2Rd −

∂ψ
1
(t, x) − (∆x ψ)(t, x).
∂t
2

(75)

Note that the solution u of the PDE (30) then satisfies for all t ∈ [0, T ), x = (x1 , x2 , . . . , xd ) ∈
Rd that u(T, x) = g(x) and
∂ψ
∂u
1
1
(t, x) + k(∇x u)(t, x)k2Rd + (∆x u)(t, x) =
(t, x) + k(∇x ψ)(t, x)k2Rd + (∆x ψ)(t, x).
∂t
2
∂t
2
(76)
On the left hand side of Figure 7 we present approximatively the relative L1 -approximatin
error associated to U Θm against m ∈ {1, 2, 3, . . . , 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 Θm against m ∈ {1, 2, 3, . . . , 4000} based on 256
Monte Carlo samples and 5 independent realizations (5 independent runs).

4.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
24

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 7: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 4000} in the case of the PDE (76). The deep BSDE approximation
U Θ4000 ≈ u(0, ξ) achieves a relative L1 -approximation error of size 0.0009 in a runtime of
957 seconds.
but in two space-dimensions (d = 2) instead of in hundred space-dimensions (d = 100) as
in this article.
√
Assume the setting in Subsection 4.1, let κ = 6/10, λ = 1/ d, assume for all s, t ∈ [0, T ],
1
x = (x1 , . . . , xd ), w = (w1 , . . . , wd ) ∈ Rd , y ∈ R, z ∈ R1×d , m ∈ N that γm = 100
= 0.01,
d
T = 1, d = 100, N = 30, µ(t,
= 0, σ(t, x)w = w, ξ = (0, 0, . . . , 0) ∈ R , Υ(s, t, x, w) =
Px)
d
x + w, g(x) = 1 + κ + sin(λ i=1 xi ), and
n 
o

Pd
λ2 d(t−T ) 2
f (t, x, y, z) = min 1, y − κ − 1 − sin λ i=1 xi exp
.
(77)
2
Note that the solution u of the PDE (30) then satisfies for all t ∈ [0, T ), x = (x1 , x2 , . . . , xd ) ∈
Rd that u(T, x) = g(x) and
n 
o

P
2
∂u
) 2
(t, x)+min 1, u(t, x)−κ−1−sin λ di=1 xi exp λ d(t−T
+ 21 (∆x u)(t, x) = 0 (78)
2
∂t
(cf. Lemma 4.4 below). On the left hand side of Figure 7 we present approximatively the
relative L1 -approximatin error associated to U Θm against m ∈ {1, 2, 3, . . . , 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 Θm against m ∈
{1, 2, 3, . . . , 24000} based on 256 Monte Carlo samples and 5 independent realizations (5
independent runs).
Lemma 4.4 (Cf. Subsection 6.1 in [14]). Let T, κ, λ ∈ (0, ∞), d ∈ N and let u : [0, T ] ×
Rd → R be the function which satisfies for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd that

P
2
)
u(t, x) = 1 + κ + sin λ di=1 xi exp λ d(t−T
.
(79)
2
25

(a) Relative L1 -approximation error

(b) Mean of the loss function

Figure 8: Relative L1 -approximation error of U Θm and mean of the loss function against
m ∈ {1, 2, 3, . . . , 24000} in the case of the PDE (78). The deep BSDE approximation
U Θ24000 ≈ u(0, ξ) achieves a relative L1 -approximation error of size 0.0060 in a runtime of
4872 seconds.
Then it holds for all t P
∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd that u ∈ C 1,2 ([0, T ] × Rd , R),
u(T, x) = 1 + κ + sin(λ di=1 xi ), and
n 

P
∂u
(t, x)+min 1, u(t, x)−κ−1−sin λ di=1 xi exp
∂t

λ2 d(t−T ) 2
2

o

+ 12 (∆x u)(t, x) = 0. (80)

Proof of Lemma 4.4. Note that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd it holds that

P
λ2 d
∂u
(t, x) =
sin λ di=1 xi exp
∂t
2

λ2 d(t−T ) 
.
2

(81)

In addition, observe that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd , k ∈ {1, 2, . . . , d} it holds
that

P
2
∂u
)
(t, x) = λ cos λ di=1 xi exp λ d(t−T
.
(82)
2
∂xk
Hence, we obtain that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd , k ∈ {1, . . . , d} it holds that

P
∂ 2u
(t, x) = −λ2 sin λ di=1 xi exp
2
∂xk

λ2 d(t−T ) 
.
2

(83)

This ensures that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd it holds that
(∆x u)(t, x) = − d λ2 sin λ

26

Pd

i=1


xi exp

λ2 d(t−T ) 
.
2

(84)

Combining this with (81) proves that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd it holds that
1
∂u
(t, x) + (∆x u)(t, x) = 0.
∂t
2

(85)

This demonstrates that for all t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd it holds that
n 

P
∂v
(t, x) + min 1, v(t, x) − κ − 1 − sin λ di=1 xi exp
∂t
∂v
1
=
(t, x) + (∆x v)(t, x) = 0.
∂t
2

λ2 d(t−T )
2

2 o

+ 21 (∆x v)(t, x)

(86)

The proof of Lemma 4.4 is thus completed.

5

Appendix A: Special cases of the proposed algorithm

In this subsection we illustrate the general algorithm in Subsection 3.2 in several special
cases. More specifically, in Subsections 5.1 and 5.2 we provide special choices for the
functions ψm , m ∈ N, and Ψm , m ∈ N, employed in (29) and in Subsections 5.3 and 5.4
we provide special choices for the function Υ in (24).

5.1

Stochastic gradient descent (SGD)

Example 5.1. Assume the setting in Subsection 3.2, let (γm )m∈N ⊆ (0, ∞), and assume
for all m ∈ N, x ∈ R% , (ϕj )j∈N ∈ (Rρ )N that
% = ρ,

Ψm (x, (ϕj )j∈N ) = ϕ1 ,

and

ψm (x) = γm x.

(87)

Then it holds for all m ∈ N that
Θm = Θm−1 − γm Φm−1,1 (Θm−1 ).

5.2

(88)

Adaptive Moment Estimation (Adam) with mini-batches

In this subsection we illustrate how the so-called Adam optimizer (see [22]) can be employed in conjunction with the deep BSDE solver in Subsection 3.2 (cf. also Subsection 4.1
above).
Example 5.2. Assume the setting in Subsection 3.2, assume that % = 2ρ, let Powr : Rρ →
Rρ , r ∈ (0, ∞), be the functions which satisfy for all r ∈ (0, ∞), x = (x1 , . . . , xρ ) ∈ Rρ that
Powr (x) = (|x1 |r , . . . , |xρ |r ),
27

(89)

let ε ∈ (0, ∞), (γm )m∈N ⊆ (0, ∞), (Jm )m∈N0 ⊆ N, X, Y ∈ (0, 1), let m, M : N0 × Ω → Rρ
be the stochastic processes which satisfy for all m ∈ N0 that Ξm = (mm , Mm ), and assume
for all m ∈ N, x, y ∈ Rρ , (ϕj )j∈N ∈ (Rρ )N that


Pm
Pm
ϕj , Yy + (1 − Y) Pow2 J1m Jj=1
ϕj
(90)
Ψm (x, y, (ϕj )j∈N ) = Xx + (1 − X) J1m Jj=1
and


ψm (x, y) = ε + Pow1/2

y
(1 − Ym )

−1

γm x
.
(1 − Xm )

(91)

Then it holds for all m ∈ N that


Θm = Θm−1 − ε + Pow1/2
mm = X mm−1 +

(1 − X)
Jm

−1
Mm
γm mm
,
m
(1 − Y )
(1 − Xm )
!
Jm
X
Φm−1,j
(Θm−1 ) ,
Sm

Mm = Y Mm−1 + (1 − Y) Pow2

5.3

(92)

j=1

!
Jm
1 X
m−1,j
(Θm−1 ) .
Φ
Jm j=1 Sm

Geometric Brownian motion

Example 5.3. Assume the setting in Section 3.2, let µ̄, σ̄ ∈ R, and assume for all s, t ∈
[0, T ], x = (x1 , . . . , xd ), w = (w1 , . . . , wd ) ∈ Rd that



σ̄ 2
Υ(s, t, x, w) = exp µ̄ −
(t − s) exp(σ̄ diagRd×d (w1 , . . . , wd )) x.
(93)
2
Then it holds for all m, j ∈ N0 , n ∈ {0, 1, . . . , N } that



σ̄ 2
m,j
θ,m,j
Xn
= exp µ̄ −
tn IdRd +σ̄ diagRd×d (Wtn ) ξ.
2

(94)

In the setting of Example 5.3 we consider under suitable further hypotheses (cf. Subsection 4.4 above) for every sufficiently large m ∈ N0 the random variable U Θm as an
approximation of u(0, ξ) where u : [0, T ] × Rd → Rk is a suitable solution of the PDE
∂u
(t, x)
∂t

+

σ̄ 2
2

d
P
i=1

|xi |2

∂2u
∂x2i

d

P
(t, x) + µ̄ xi
i=1

∂u
∂xi


(t, x)


+ f t, x, u(t, x), σ̄ ( ∂u
)(t, x) diagRd×d (x1 , . . . , xd ) = 0 (95)
∂x
with u(T, x) = g(x) for t ∈ [0, T ], x = (x1 , . . . , xd ) ∈ Rd .
28

5.4

Euler-Maruyama scheme

Example 5.4. Assume the setting in Section 3.2, let µ : [0, T ] × Rd → Rd and σ : [0, T ] ×
Rd → Rd be functions, and assume for all s, t ∈ [0, T ], x, w ∈ Rd that
Υ(s, t, x, w) = x + µ(s, x) (t − s) + σ(s, x) w.

(96)

Then it holds for all m, j ∈ N0 , n ∈ {0, 1, . . . , N − 1} that
Xnm,j = Xnm,j + µ(tn , Xnm,j ) (tn+1 − tn ) + σ(tn , Xnm,j ) (Wtn+1 − Wtn ).

(97)

In the setting of Example 5.4 we consider under suitable further hypotheses for every
sufficiently large m ∈ N0 the random variable U Θm as an approximation of u(0, ξ) where
u : [0, T ] × Rd → Rk is a suitable solution of the PDE
∂u
(t, x)
∂t

+

1
2

d
P

2

(d)

(d) 

( ∂∂xu2 )(t, x) σ(t, x) ej , σ(t, x) ej

j=1

+ ( ∂u
)(t, x) µ(t, x)
∂x


+ f t, x, u(t, x), ( ∂u
)(t, x) σ(t, x) = 0 (98)
∂x
(d)

(d)

with u(T, x) = g(x), e1 = (1, 0, . . . , 0), . . . , ed = (0, . . . , 0, 1) ∈ Rd for t ∈ [0, T ],
x = (x1 , . . . , xd ) ∈ Rd (cf. (PDE) in Section 2 above).

6

Appendix B: Python and Matlab source codes

6.1

1
2
3
4
5
6
7
8
9
10
11

Python source code for an implementation of the deep BSDE
solver in the case of the Allen-Cahn PDE (35) in Subsection 4.2

import time
import math
import tensorflow as tf
import numpy as np
from tensorflow . python . training . moving_averages \
import a s s i g n _ m o v i n g _ a v e r a g e
from scipy . stats import m ul t i va r i at e _ no r m al as normal
from tensorflow . python . ops import control_flow_ops
from tensorflow import r a n d o m _ n o r m a l _ i n i t i a l i z e r as norm_init
from tensorflow import r a n d o m _ u n i f o r m _ i n i t i a l i z e r as unif_init
from tensorflow import c o n st a n t _ i n i t i a l i z e r as const_init

12
13
14
15

class SolveAllenCahn ( object ):
""" The fully - connected neural network model ."""
def __init__ ( self , sess ):

29

16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

self . sess = sess
# parameters for the PDE
self . d = 100
self . T = 0.3
# parameters for the algorithm
self . n_time = 20
self . n_layer = 4
self . n_neuron = [ self .d , self . d +10 , self . d +10 , self . d ]
self . batch_size = 64
self . valid_size = 256
self . n_maxstep = 4000
self . n_displaystep = 100
self . learning_rate = 5e -4
self . Yini = [0.3 , 0.6]
# some basic constants and variables
self . h = ( self . T +0.0)/ self . n_time
self . sqrth = math . sqrt ( self . h )
self . t_stamp = np . arange (0 , self . n_time )* self . h
self . _extra_train_ops = []

35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

def train ( self ):
start_time = time . time ()
# train operations
self . global_step = \
tf . get_variable ( ’ global_step ’ , [] ,
initializer = tf . c o n s t a n t _ i n i t i a l i z e r (1) ,
trainable = False , dtype = tf . int32 )
trainable_vars = tf . t ra i n ab l e _v a r ia b l es ()
grads = tf . gradients ( self . loss , trainable_vars )
optimizer = tf . train . AdamOptimizer ( self . learning_rate )
apply_op = \
optimizer . apply_gradients ( zip ( grads , trainable_vars ) ,
global_step = self . global_step )
train_ops = [ apply_op ] + self . _extra_train_ops
self . train_op = tf . group (* train_ops )
self . loss_history = []
self . init_history = []
# for validation
dW_valid , X_valid = self . sample_path ( self . valid_size )
feed_dict_valid = { self . dW : dW_valid ,
self . X : X_valid ,
self . is_training : False }
# initialization
step = 1
self . sess . run ( tf . g l o b a l _ v a r i a b l e s _ i n i t i a l i z e r ())
temp_loss = self . sess . run ( self . loss ,
feed_dict = feed_dict_valid )
temp_init = self . Y0 . eval ()[0]
self . loss_history . append ( temp_loss )

30

65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

self . init_history . append ( temp_init )
print " step : % 5u , loss : %.4 e , " % \
(0 , temp_loss ) + \
" Y0 : % .4 e , runtime : %4 u " % \
( temp_init , time . time () - start_time + self . t_bd )
# begin sgd iteration
for _ in range ( self . n_maxstep +1):
step = self . sess . run ( self . global_step )
dW_train , X_train = self . sample_path ( self . batch_size )
self . sess . run ( self . train_op ,
feed_dict ={ self . dW : dW_train ,
self . X : X_train ,
self . is_training : True })
if step % self . n_displaystep == 0:
temp_loss = self . sess . run ( self . loss ,
feed_dict = feed_dict_valid )
temp_init = self . Y0 . eval ()[0]
self . loss_history . append ( temp_loss )
self . init_history . append ( temp_init )
print " step : % 5u , loss : %.4 e , " % \
( step , temp_loss ) + \
" Y0 : % .4 e , runtime : %4 u " % \
( temp_init , time . time () - start_time + self . t_bd )
step += 1
end_time = time . time ()
print " running time : % .3 f s " % \
( end_time - start_time + self . t_bd )

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

def build ( self ):
start_time = time . time ()
# build the whole network by stacking subnetworks
self . dW = tf . placeholder ( tf . float64 ,
[ None , self .d , self . n_time ] ,
name = ’ dW ’)
self . X = tf . placeholder ( tf . float64 ,
[ None , self .d , self . n_time +1] ,
name = ’X ’)
self . is_training = tf . placeholder ( tf . bool )
self . Y0 = tf . Variable ( tf . random_uniform ([1] ,
minval = self . Yini [0] ,
maxval = self . Yini [1] ,
dtype = tf . float64 ));
self . Z0 = tf . Variable ( tf . random_uniform ([1 , self . d ] ,
minval = -.1 ,
maxval =.1 ,
dtype = tf . float64 ))
self . allones = \
tf . ones ( shape = tf . pack ([ tf . shape ( self . dW )[0] , 1]) ,
dtype = tf . float64 )

31

114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

Y = self . allones * self . Y0
Z = tf . matmul ( self . allones , self . Z0 )
with tf . variable_scope ( ’ forward ’ ):
for t in xrange (0 , self . n_time -1):
Y = Y - self . f_tf ( self . t_stamp [ t ] ,
self . X [: , : , t ] , Y , Z )* self . h
Y = Y + tf . reduce_sum ( Z * self . dW [: , : , t ] , 1 ,
keep_dims = True )
Z = self . _one_time_net ( self . X [: , : , t +1] ,
str ( t +1))/ self . d
# terminal time
Y = Y - self . f_tf ( self . t_stamp [ self . n_time -1] ,
self . X [: , : , self . n_time -1] ,
Y , Z )* self . h
Y = Y + tf . reduce_sum ( Z * self . dW [: , : , self . n_time -1] , 1 ,
keep_dims = True )
term_delta = Y - self . g_tf ( self .T ,
self . X [: , : , self . n_time ])
self . clipped_delta = \
tf . clip_by_value ( term_delta , -50.0 , 50.0)
self . loss = tf . reduce_mean ( self . clipped_delta **2)
self . t_bd = time . time () - start_time

136
137
138
139
140
141
142
143
144
145
146
147
148

def sample_path ( self , n_sample ):
dW_sample = np . zeros ([ n_sample , self .d , self . n_time ])
X_sample = np . zeros ([ n_sample , self .d , self . n_time +1])
for i in xrange ( self . n_time ):
dW_sample [: , : , i ] = \
np . reshape ( normal . rvs ( mean = np . zeros ( self . d ) ,
cov =1 ,
size = n_sample )* self . sqrth ,
( n_sample , self . d ))
X_sample [: , : , i +1] = X_sample [: , : , i ] + \
np . sqrt (2) * dW_sample [: , : , i ]
return dW_sample , X_sample

149
150
151
152

def f_tf ( self , t , X , Y , Z ):
# nonlinear term
return Y - tf . pow (Y , 3)

153
154
155
156

def g_tf ( self , t , X ):
# terminal conditions
return 0.5/(1 + 0.2* tf . reduce_sum ( X **2 , 1 , keep_dims = True ))

157
158
159
160
161
162

def _one_time_net ( self , x , name ):
with tf . variable_scope ( name ):
x_norm = self . _batch_norm (x , name = ’ layer0_normal ’)
layer1 = self . _one_layer ( x_norm , self . n_neuron [1] ,
name = ’ layer1 ’)

32

163
164
165
166
167

layer2 = self . _one_layer ( layer1 , self . n_neuron [2] ,
name = ’ layer2 ’)
z = self . _one_layer ( layer2 , self . n_neuron [3] ,
activation_fn = None , name = ’ final ’)
return z

168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183

def _one_layer ( self , input_ , out_sz ,
activation_fn = tf . nn . relu ,
std =5.0 , name = ’ linear ’ ):
with tf . variable_scope ( name ):
shape = input_ . get_shape (). as_list ()
w = tf . get_variable ( ’ Matrix ’ ,
[ shape [1] , out_sz ] , tf . float64 ,
norm_init ( stddev = \
std / np . sqrt ( shape [1]+ out_sz )))
hidden = tf . matmul ( input_ , w )
hidden_bn = self . _batch_norm ( hidden , name = ’ normal ’)
if activation_fn != None :
return activation_fn ( hidden_bn )
else :
return hidden_bn

184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211

def _batch_norm ( self , x , name ):
""" Batch normalization """
with tf . variable_scope ( name ):
params_shape = [ x . get_shape ()[ -1]]
beta = tf . get_variable ( ’ beta ’ , params_shape ,
tf . float64 ,
norm_init (0.0 , stddev =0.1 ,
dtype = tf . float64 ))
gamma = tf . get_variable ( ’ gamma ’ , params_shape ,
tf . float64 ,
unif_init (0.1 , 0.5 ,
dtype = tf . float64 ))
mv_mean = tf . get_variable ( ’ moving_mean ’ ,
params_shape ,
tf . float64 ,
const_init (0.0 , tf . float64 ) ,
trainable = False )
mv_var = tf . get_variable ( ’ moving_variance ’ ,
params_shape ,
tf . float64 ,
const_init (1.0 , tf . float64 ) ,
trainable = False )
# These ops will only be preformed when training
mean , variance = tf . nn . moments (x , [0] , name = ’ moments ’)
self . _extra_train_ops . append (\
a s s i g n _ m o v i n g _ a v e r a g e ( mv_mean , mean , 0.99))
self . _extra_train_ops . append (\

33

a s s i g n _ m o v i n g _ a v e r a g e ( mv_var , variance , 0.99))
mean , variance = \
control_flow_ops . cond ( self . is_training ,
lambda : ( mean , variance ) ,
lambda : ( mv_mean , mv_var ))
y = tf . nn . b a t ch _ n or m a li z a ti o n (x , mean , variance ,
beta , gamma , 1e -6)
y . set_shape ( x . get_shape ())
return y

212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241

def main ():
tf . re s e t_ d e fa u l t _g r a ph ()
with tf . Session () as sess :
tf . set_random_seed (1)
print " Begin to solve Allen - Cahn equation "
model = SolveAllenCahn ( sess )
model . build ()
model . train ()
output = np . zeros (( len ( model . init_history ) , 3))
output [: , 0] = np . arange ( len ( model . init_history )) \
* model . n_displaystep
output [: , 1] = model . loss_history
output [: , 2] = model . init_history
np . savetxt ("./ AllenCahn_d100 . csv " ,
output ,
fmt =[ ’% d ’ , ’ %.5 e ’ , ’ %.5 e ’] ,
delimiter =" ," ,
header =" step , loss function , " + \
" target value , runtime " ,
comments = ’ ’)

242
243
244
245

if __name__ == ’ __main__ ’:
np . random . seed (1)
main ()

Matlab code 1: A Python code for the deep BSDE solver in Subsection 3.2 in the case
of the PDE (35).

6.2

1
2
3
4
5

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

function Branching_Matlab ()
% Parameters for the model
T = 0.3; t0 = 0; x0 = 0; d = 100; m = d ;
mu = zeros (d ,1); sigma = eye ( d )* sqrt (2);
a = [0 2 0 -1] ’;

34

6

g = @ ( x ) 1./(1+0.2* norm ( x )^2)*1/2;

7
8
9

% Parameters for the algorithm
rng ( ’ default ’ ); M = 10^7; beta = 1; p = [0 0.5 0 0.5] ’;

10
11
12
13
14

% Branching method
tic ;
[ mn , sd ] = MC_BM ( mu , sigma , beta , p , a , t0 , x0 , T , g , M );
runtime = toc ;

15
16
17
18
19
20
21
22

% Output
disp ([ ’ Terminal condition : u (T , x0 ) = ’ num2str ( g ( x0 )) ’; ’ ]);
disp ([ ’ Branching method : u (0 , x0 ) ~ ’ num2str ( mn ) ’; ’ ]);
disp ([ ’ Estimated standard deviation : ’ num2str ( sd ) ’; ’ ]);
disp ([ ’ Estimated L2 - appr . error = ’ num2str ( sd / sqrt ( M )) ’; ’ ]);
disp ([ ’ Elapsed runtime = ’ num2str ( runtime ) ’; ’ ]);
end

23
24
25
26
27
28
29
30
31
32

function [ mn , sd ] = MC_BM ( mu , sigma , beta , p , a , t0 , x0 , T , g , M )
mn = 0; sd = 0;
for m =1: M
result = BM_Eval ( mu , sigma , beta , p , a , t0 , x0 , T , g );
mn = mn + result ;
sd = sd + result ^2;
end
mn = mn / M ; sd = sqrt ( ( sd - mn ^2/ M )/ M );
end

33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

function result = BM_Eval ( mu , sigma , beta , p , a , t0 , x0 , T , g )
bp = BP ( mu , sigma , beta , p , t0 , x0 , T );
result = 1;
for k =1: size ( bp {1} ,2)
result = result * g ( bp {1}(: , k ) );
end
if norm (a - p ) > 0
for k =1: length ( a )
if p ( k ) > 0
result = result * ( a ( k )/ p ( k ) )^( bp {2}( k ) );
elseif a ( k ) ~= 0
error ( ’a ( k ) zero but p ( k ) non - zero ’ );
end
end
end
end

50
51
52
53
54

function bp = BP ( mu , sigma , beta , p , t0 , x0 , T )
bp = cell (2 ,1);
bp {2} = p *0;
tau = exprnd (1/ beta );

35

55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

new_t0 = min ( tau + t0 , T );
delta_t = new_t0 - t0 ;
m = size ( sigma ,2);
new_x0 = x0 + mu * delta_t + sigma * sqrt ( delta_t )* randn (m ,1);
if tau >= T - t0
bp {1} = new_x0 ;
else
[ tmp , nonlinearity ] = max ( mnrnd (1 , p ));
bp {2}( nonlinearity ) = bp {2}( nonlinearity ) + 1;
for k =1: nonlinearity -1
tmp = BP ( mu , sigma , beta , p , new_t0 , new_x0 , T );
bp {1} = [ bp {1} tmp {1} ];
bp {2} = bp {2} + tmp {2};
end
end
end

Matlab code 2: A Matlab code for the Branching diffusion method in the case of the
PDE (35) based on M = 107 independent realizations.

6.3

1
2

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

function Mont eCarlo _Matla b ()
rng ( ’ default ’ );

3
4
5
6
7
8
9

% Parameters for the model
d = 100;
g = @ ( x ) log ( (1+ norm ( x )^2)/2 );
T = 1;
M = 10^7;
t = 0;

10
11
12
13
14
15
16
17
18
19

% Classical Monte Carlo
tic ;
MC = 0;
for m =1: M
dW = randn (1 , d )* sqrt (T - t );
MC = MC + exp ( - g ( dW * sqrt (2) ) );
end
MC = - log ( MC / M );
runtime = toc ;

20
21
22
23

% Output
disp ([ ’ Solution : u (T ,0) = ’ num2str ( g (0)) ’; ’ ]);
disp ([ ’ Solution : u (0 ,0) = ’ num2str ( MC ) ’; ’ ]);

36

24
25

disp ([ ’ Time = ’ num2str ( runtime ) ’; ’ ]);
end

Matlab code 3: A Matlab code for a Monte Carlo method related to the PDE (36)
based on M = 107 independent realizations.
Acknowledgements
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
[1] Bellman, R. Dynamic programming. Princeton Landmarks in Mathematics. Princeton University Press, Princeton, NJ, 2010. Reprint of the 1957 edition, With a new
introduction by Stuart Dreyfus.
[2] Bender, C., and Denk, R. A forward scheme for backward SDEs. Stochastic
Processes and their Applications 117, 12 (2007), 1793–1812.
[3] Bender, C., Schweizer, N., and Zhuo, J. A primal-dual algorithm for BSDEs.
arXiv:1310.3694 (2014), 36 pages.
[4] Bergman, Y. Z. Option pricing with differential interest rates. Review of Financial
Studies 8, 2 (1995), 475–500.
[5] Briand, P., and Labart, C. Simulation of BSDEs by Wiener chaos expansion.
Ann. Appl. Probab. 24, 3 (06 2014), 1129–1171.
[6] Chassagneux, J.-F. Linear multistep schemes for BSDEs. SIAM J. Numer. Anal.
52, 6 (2014), 2815–2836.
[7] Chassagneux, J.-F., and Richou, A. Numerical simulation of quadratic BSDEs.
Ann. Appl. Probab. 26, 1 (2016), 262–304.
[8] Crisan, D., and Manolarakis, K. Solving backward stochastic differential equations using the cubature method: Application to nonlinear pricing. SIAM Journal on
Financial Mathematics 3, 1 (2012), 534–571.

37

[9] Darbon, J., and Osher, S. Algorithms for overcoming the curse of dimensionality
for certain Hamilton-Jacobi equations arising in control theory and elsewhere. Res.
Math. Sci. 3 (2016), Paper No. 19, 26.
[10] Debnath, L. Nonlinear partial differential equations for scientists and engineers,
third ed. Birkhäuser/Springer, New York, 2012.
[11] E, W., Hutzenthaler, M., Jentzen, A., and Kruse, T. On full history
recursive multilevel Picard approximations and numerical approximations for highdimensional nonlinear parabolic partial differential equations and high-dimensional
nonlinear backward stochastic differential equations. arXiv:1607.03295 (2017), 46
pages.
[12] Gobet, E., Lemor, J.-P., and Warin, X. A regression-based Monte Carlo method
to solve backward stochastic differential equations. Ann. Appl. Probab. 15, 3 (2005),
2172–2202.
[13] Gobet, E., and Turkedjiev, P. Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions. Math. Comp. 85, 299
(2016), 1359–1391.
[14] Gobet, E., and Turkedjiev, P. Adaptive importance sampling in least-squares
Monte Carlo algorithms for backward stochastic differential equations. Stochastic
Process. Appl. 127, 4 (2017), 1171–1203.
[15] Goodfellow, I., Bengio, Y., and Courville, A. Deep Learning. MIT Press,
2016. http://www.deeplearningbook.org.
[16] Han, J., and E, W. Deep Learning Approximation for Stochastic Control Problems.
arXiv:1611.07422 (2016), 9 pages.
[17] Henry-Labordère, P. Counterparty risk valuation: a marked branching diffusion
approach. arXiv:1203.2369 (2012), 17 pages.
[18] Henry-Labordère, P., Oudjane, N., Tan, X., Touzi, N., and Warin, X.
Branching diffusion representation of semilinear PDEs and Monte Carlo approximation. arXiv:1603.01727 (2016), 30 pages.
[19] Henry-Labordère, P., Tan, X., and Touzi, N. A numerical algorithm for a
class of BSDEs via the branching process. Stochastic Process. Appl. 124, 2 (2014),
1112–1140.
[20] Hinton, G. E., Deng, L., Yu, D., Dahl, G., Mohamed, A., Jaitly, N.,
Senior, A., Vanhoucke, V., Nguyen, P., Sainath, T., and Kingsbury, B.
38

Deep neural networks for acoustic modeling in speech recognition. Signal Processing
Magazine 29 (2012), 82–97.
[21] Ioffe, S., and Szegedy, C. Batch normalization: accelerating deep network training by reducing internal covariate shift. Proceedings of The 32nd International Conference on Machine Learning (ICML), June 2015.
[22] Kingma, D., and Ba, J. Adam: a method for stochastic optimization. Proceedings
of the International Conference on Learning Representations (ICLR), May 2015.
[23] Krizhevsky, A., Sutskever, I., and Hinton, G. E. Imagenet classification
with deep convolutional neural networks. Advances in Neural Information Processing
Systems 25 (2012), 1097–1105.
[24] LeCun, Y., Bengio, Y., and Hinton, G. Deep learning. Nature 521 (2015),
436–444.
[25] Pardoux, É., and Peng, S. Backward stochastic differential equations and quasilinear parabolic partial differential equations. In Stochastic partial differential equations and their applications (Charlotte, NC, 1991), vol. 176 of Lecture Notes in Control
and Inform. Sci. Springer, Berlin, 1992, pp. 200–217.
[26] Pardoux, É., and Peng, S. G. Adapted solution of a backward stochastic differential equation. Systems Control Lett. 14, 1 (1990), 55–61.
[27] Pardoux, E., and Tang, S. Forward-backward stochastic differential equations and
quasilinear parabolic PDEs. Probab. Theory Related Fields 114, 2 (1999), 123–150.
[28] Peng, S. G. Probabilistic interpretation for systems of quasilinear parabolic partial
differential equations. Stochastics Stochastics Rep. 37, 1-2 (1991), 61–74.

39

