# Generalization Error Analysis of Deep Backward Dynamic Programming for
Solving Nonlinear PDEs

Du Ouyang<EMAIL_ADDRESS>Jichang Xiao
<EMAIL_ADDRESS>Xiaoqun Wang<EMAIL_ADDRESS>
###### Abstract

We explore the application of the quasi-Monte Carlo (QMC) method in deep
backward dynamic programming (DBDP) [1] for numerically solving high-
dimensional nonlinear partial differential equations (PDEs). Our study focuses
on examining the generalization error as a component of the total error in the
DBDP framework, discovering that the rate of convergence for the
generalization error is influenced by the choice of sampling methods.
Specifically, for a given batch size $m$, the generalization error under QMC
methods exhibits a convergence rate of $O(m^{-1+\varepsilon})$, where
$\varepsilon$ can be made arbitrarily small. This rate is notably more
favorable than that of the traditional Monte Carlo (MC) methods, which is
$O(m^{-1/2+\varepsilon})$. Our theoretical analysis shows that the
generalization error under QMC methods achieves a higher order of convergence
than their MC counterparts. Numerical experiments demonstrate that QMC indeed
surpasses MC in delivering solutions that are both more precise and stable.

###### keywords:

Neural networks , Nonlinear PDEs , Quasi-Monte Carlo , Backward stochastic
differential equations

###### MSC:

[2020] 41A63 , 65M75 , 65D40

[label1]organization=Department of Mathematical Sciences, Tsinghua
University,city=Beijing, postcode=100084, country=People’s Republic of China

##  1 Introduction

Nonlinear parabolic partial differential equations (PDEs) represent a general
form of linear Kolmogorov PDEs [2, 3, 4, 5, 6], which find wide applications
in financial mathematics [7]. The nonlinear PDE has the following form

$\left\\{\begin{aligned} \frac{\partial u}{\partial
t}+\mathscr{A}u+f(t,x,u,\Sigma^{\top}D_{x}u)&=0,\;\;\;\;\;\;\text{ for
}(t,x)\in[0,T)\times\mathbb{R}^{d},\\\ u(T,x)&=g(x),\;\;\;\;\;\mbox{ for
}x\in\mathbb{R}^{d},\end{aligned}\right.$ (1)

where
$f:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
and $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is called the terminal function.
The second-order generator $\mathscr{A}$ has the form

$\mathscr{A}:=\frac{1}{2}\sum_{i,j=1}^{d}(\Sigma\Sigma^{\top})_{ij}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}+\sum_{i=1}^{d}M_{i}\frac{\partial}{\partial x_{i}},$ (2)

where $(M_{1}(t,x),\dots,M_{d}(t,x))\coloneqq M(t,x)$ is a function defined on
$[0,T]\times\mathbb{R}^{d}$ with values in $\mathbb{R}^{d}$, $\Sigma(t,x)$ is
a function defined on $[0,T]\times\mathbb{R}^{d}$ with values in
$\mathbb{M}^{d}$, the set of $d\times d$ matrices.

The objective of this paper is to obtain the solution of high-dimensional PDE
(1) at $t=0$ over a region $[a,b]^{d}$, i.e., we aim to obtain $u(0,x)$ for
$x\in[a,b]^{d}$. Thanks to the nonlinear Feynman-Kac formula obtained by
Pardoux and Peng [8], the solution of the PDE (1) is connected to the backward
stochastic differential equations (BSDEs). The BSDE was first proposed by
Pardoux and Peng [9], and has the following form

$Y_{t}=g(\mathcal{X}_{T})+\int_{t}^{T}f(s,\mathcal{X}_{s},Y_{s},Z_{s})\mathrm{d}s-\int_{t}^{T}Z_{s}^{\top}\mathrm{d}W_{s}\
,\ \ \ 0\leq t\leq T,$ (3)

where $f$ and $g$ are identical to those in PDE (1). The
$\\{\mathcal{X}_{t}\\}_{0\leq t\leq T}$ satisfies the following forward
stochastic differential equation (SDE)

$\mathcal{X}_{t}=\eta+\int_{0}^{t}M(s,\mathcal{X}_{s})\mathrm{d}s+\int_{0}^{t}\Sigma(s,\mathcal{X}_{s})\mathrm{d}W_{s},\;\;\;0\leq
t\leq T,$ (4)

where $W$ is a $d$-dimensional Brownian motion and $\eta$ is uniformly
distributed on $\left[a,b\right]^{d}$. Consider the probability space
$(\Omega,\mathcal{F},P)$ equipped with a filtration
$\\{\mathcal{F}_{t}\\}_{0\leq t\leq T}$ generated by the Brownian motion $W$
and the initial random variable $\eta$ satisfying the usual conditions. $M$
and $\Sigma$ are idential to those in (1).

Under appropriate conditions (see [8] and [7]), the solutions of (3) and (1)
are connected in the form of

$Y_{t}=u(t,\mathcal{X}_{t}),\ \ 0\leq t\leq T,$ (5)

and if there are more smooth conditions, then

$Z_{t}=\Sigma^{\top}(t,\mathcal{X}_{t})\partial_{x}u(t,\mathcal{X}_{t}),\ \
0\leq t\leq T.$ (6)

The greatest advantage of the Feynman-Kac formula is its ability to transform
the solution of PDEs into the solution of BSDEs. This is beneficial because
general methods for solving PDEs, such as finite element methods, depend on
the partition of the state space and are subject to the “curse of
dimensionality”. The schemes in [10, 11, 12], using the time partition and
regression, overcome the issue of dimensionality. These schemes can be applied
up to dimension 6 or 7. However, the efficacy of the regression methods is
constrained by the choice of the basis functions (e.g., polynomial functions).
In cases of higher dimensions (e.g., 50), the sheer volume of necessary basis
functions renders these approaches highly inefficient.

To address this issue, the pioneering deep BSDE method [13, 14] applies deep
neural networks to efficiently solve high-dimensional BSDEs (nonlinear PDEs).
This approach utilizes the required initial values as parameters and employs
feedforward neural networks to approximate each time step to solve a global
optimization problem. In contrast, the Deep Backward Dynamic Programming
(DBDP) [1] simultaneously estimates the solution and its gradient through the
minimization of sequential loss functions via backward induction, showing
superior stability and accuracy compared to [13], as evidenced by several test
cases. Based on [1], Germain et al. [15] proposed a multistep version of DBDP.
The deep splitting (DS) method [16] estimates the PDE solution via backward
explicit local optimization problems, leveraging a neural network regression
method to calculate conditional expectations.

In the aforementioned methods, the loss function is typically represented by
an expectation of the states, which cannot be precisely calculated. To address
this issue, empirical risk minimization (ERM) principle aims to apply the
Monte Carlo (MC) method to sample states and using the sample mean (also
called empirical risk) to estimate the loss function. Then the problem becomes
finding a parameter that minimizes the empirical risk. Due to the ERM, the
total error of the algorithm will be influenced by the generalization error
(see Section 3). The generalization error (also called stochastic error in
[17, 18]), which correlates with the sampling method and sample size,
signifies the capacity of ERM. As a component of the total error, it affects
the performance of algorithm during the learning phase. Therefore, studying
the convergence order of the generalization error, which is not involved in
[16, 15, 1], is also important. Other research concerning the generalization
error can be found in [19, 20, 4, 21, 17, 18, 22].

The quasi-Monte Carlo (QMC) method is an efficient numerical method widely
used in financial engineering (see [23, 24]). For $m$ quadrature points, the
QMC method has a convergence rate $O(m^{-1+\varepsilon})$ with an arbitrarily
small $\varepsilon>0$, which is asymptotically faster than the MC rate
$O(m^{-1/2})$. The applications of QMC methods in machine learning have also
made some progress. Dick and Feischl [25] applied QMC methods to reduce the
computational effort required for data compression in machine learning. Liu
and Owen [26] demonstrated that randomized quasi-Monte Carlo (RQMC) methods
combined with stochastic quasi-Newton methods can significantly accelerate the
optimization. Longo et al. [27] and Mishra et al. [28] utilized QMC methods to
generate training data for training neural networks. Their results indicate
that the QMC-based deep learning algorithm is more efficient than the MC-based
one in approximating data-to-observables maps. Thus it is natural to ask
whether the QMC methods can achieve higher efficiency compared to MC in
solving nonlinear PDEs. Xiao et al [29] studied the efficiency of RQMC methods
for solving linear Kolmogorov PDEs by using deep learning. However, to the
best of our knowledge, there is no related research for the nonlinear PDEs
cases.

To this end, we decompose the total error of the DBDP method for solving
nonlinear PDEs into three parts: the first part is the scheme error, which
relates to the time partition; the second part is the approximation error,
related to the structure of the network; and the third part is the
generalization error (ignored by [1]), which is related to the sampling method
and batch size. We study the convergence rate of the generalization error
under MC and QMC frameworks. Our theoretical findings demonstrate that under
QMC, the convergence rate of generalization error is $O(m^{-1+\varepsilon})$
with $\varepsilon>0$ being arbitrarily small, which is better than the MC case
$O(m^{-1/2+\varepsilon})$. The efficacy of the QMC method is validated through
a series of experiments. The numerical results reveal that the QMC method
yields training outcomes with significantly reduced relative error in
comparison to the MC method when evaluated against the analytical solution.
Moreover, the results obtained by the QMC method exhibit enhanced stability,
characterized by a notably lower variance.

The structure of this paper is as follows. Section 2 introduces the neural
network and the DBDP scheme. Section 3 decomposes the total error of the DBDP
scheme and derives the generalization error. Section 4 studies the convergence
rates of generalization error under the MC and QMC frameworks. Section 5
presents the numerical result of several examples. Section 6 concludes this
paper.

##  2 Solving nonlinear PDE by using deep learning

Consider the time grid
$\pi=\\{(t_{0},\dots,t_{N}):0=t_{0}<t_{1}<\dots<t_{N}=T\\}$ with
$|\pi|=\max_{0\leq i\leq N-1}|t_{i+1}-t_{i}|$. Recall the nonlinear Feynman-
Kac formula (5) and (6). The approximation of the solution $Y_{t_{i}}$ and
$Z_{t_{i}}$ boils down to two parts. First, we use $X_{t_{i}}$ (this can be
obtained by Euler-Maruyama or Milstein methods, see [30, 31]) to estimate
$\mathcal{X}_{t_{i}}$. Second, we approximate the function $u(t,\cdot)$ by
neural network $\mathcal{U}$.

###  2.1 Neural network

######  Definition 1 (Feedforward Neutral Network).

Let $d,d_{1}\in\mathbb{N}$, a feedforward neural network is a function of the
form

$x\in\mathbb{R}^{d}\mapsto A_{k}\circ\rho\circ
A_{k-1}\circ\cdots\circ\rho\circ A_{1}(x)\in\mathbb{R}^{d_{1}},$ (7)

where $\rho:\mathbb{R}\rightarrow\mathbb{R}$ called activation function is a
nonlinear continuous function and applied component-wise on the vectors, i.e.,
$\rho(x_{1},\dots,x_{m})=(\rho(x_{1}),\dots,\rho(x_{m}))$, and every
$A_{i}:\mathbb{R}^{l_{i}}\rightarrow\mathbb{R}^{l_{i+1}},i=1,\dots,L$, is an
affine transformation given by

$A_{i}(x)=Q_{i}x+b_{i}.$ (8)

In the above equation, $Q_{i}$ and $b_{i}$ are parameters called weight and
bias, respectively. Denote $\theta=(Q_{i},b_{i})_{i=1}^{k}$ be the collections
of weight and bias, we can represent the functions in (7) by
$\mathcal{U}(x;\theta)$.

In this paper, we choose the activation function $\rho$ to be the hyperbolic
tangent function, i.e.,

$\rho(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}.$ (9)

Let $\mathcal{S}=(l_{1},\cdots,l_{k+1})$ be the structure of the network with
$l_{1}=d$ and $l_{k+1}=d_{1}$. Denote $D(\mathcal{S})=k$ be the depth of the
network, $\left\|\mathcal{S}\right\|_{\infty}=\max_{1\leq i\leq k+1}l_{i}$ be
the width of the network and $|\mathcal{S}|=\sum_{i=1}^{k}(l_{i}\times
l_{i+1}+l_{i+1})$ be the total number of the parameters of the network. We
consider the parameter space with bound $R$ and denote
$\Theta:=\\{\theta\mid\|\theta\|_{\infty}\leq R\\}$. Denote the class of
neural networks (7) with the same structure $\mathcal{S}$ by
$\mathcal{N}_{\mathcal{S},R}=\\{\mathcal{U}(x;\theta)\mid\theta\in\Theta\\}$.

The next lemma states some basic properties of the neural networks in
Definition 1. Similar results can be found in [20, 17].

######  Lemma 1.

Under Definition 1, suppose that $l_{k+1}=1$, for every
$\theta_{1},\theta_{2}\in\Theta$ and every $x\in\mathbb{R}^{d}$, we have

$\left|\mathcal{U}(x;\theta_{1})-\mathcal{U}(x;\theta_{2})\right|\leq
C_{\mathcal{S},R}\max\\{\left\|x\right\|_{\infty},1\\}\left\|\theta_{1}-\theta_{2}\right\|_{\infty},$

and

$\left|\mathcal{U}(x;\theta_{1})\right|\leq B_{\mathcal{S},R},$

where
$C_{\mathcal{S},R}=|\mathcal{S}|R^{D(\mathcal{S})-1}\left\|\mathcal{S}\right\|_{\infty}^{D(\mathcal{S})-1}$
and $B_{\mathcal{S},R}=2R\left\|\mathcal{S}\right\|_{\infty}$.

###### Proof.

This lemma is a natural extension of [17, Lemma 5.11]. ∎

###  2.2 DBDP schemes

We consider the estimate of $Y_{t_{i}}$ as follows

$\displaystyle Y_{t_{i}}\approx\mathcal{U}_{i}(X_{t_{i}}),$ (10)

where $\mathcal{U}_{i}$ is a neural network having the form (7) and
$X_{t_{i}}$ is the approximation of $\mathcal{X}_{t_{i}}$. If SDE (4) does not
have an analytical solution, then $X_{t_{i}}$ can be obtained by some
numerical schemes, such as the Euler-Maruyama method and the Milstein method
(see [30, 31]).

The deep backward dynamic programming (DBDP) scheme (the DBDP1 scheme in [1])
is training $\mathcal{U}_{i}$ backward in time and approximating $Z_{t_{i}}$
by another neural network $\mathcal{Z}_{i}(X_{t_{i}})$. Assume that

$F(t,x,y,z,h,\Delta):=y-f(t,x,y,z)h+z^{\top}\Delta.$ (11)

For functions $\widehat{\mathcal{U}}$, $\mathcal{U}$ and $\mathcal{Z}$, define

$\displaystyle
L_{i}(\widehat{\mathcal{U}},\mathcal{U},\mathcal{Z}):=\mathbb{E}\big{|}\widehat{\mathcal{U}}(X_{t_{i+1}})-F(t_{i},X_{t_{i}},\mathcal{U}(X_{t_{i}}),\mathcal{Z}(X_{t_{i}}),\Delta
t_{i},\Delta W_{t_{i}})\big{|}^{2}.$ (12)

Assume $\mathcal{S}_{1}=(d,l_{2},\dots,l_{k},1)$ and
$\mathcal{S}_{2}=(d,l_{2},\dots,l_{k},d)$. DBDP has the following form,

  * •

Initialize $\widehat{\mathcal{U}}_{N}=g$.

  * •

For $i=N-1,\dots,0$, given $\widehat{\mathcal{U}}_{i+1}$, use
$(\mathcal{U}_{i}(\cdot,\theta),\mathcal{Z}_{i}(\cdot,\theta))\in\mathcal{N}_{\mathcal{S}_{1},R}\times\mathcal{N}_{\mathcal{S}_{2},R}$
for the approximation of
$(u(t_{i},\cdot),\Sigma^{\intercal}(t_{i},\cdot)\partial_{x}u(t_{i},\cdot))$,
and find (for example, by SGD) the optimal $\theta$ to minimize the loss
function

$\mathcal{L}_{i}(\theta)=L_{i}(\widehat{\mathcal{U}}_{i+1},\mathcal{U}(\cdot;\theta),\mathcal{Z}(\cdot;\theta)).$
(13)

  * •

Find the optimal $\theta_{i}^{*}\in\arg\min_{\theta}\mathcal{L}_{i}(\theta)$
and update $\widehat{\mathcal{U}}_{i}=\mathcal{U}_{i}(\cdot;\theta_{i}^{*})$.

##  3 Total error of DBDP

In this section, we study the total error of the DBDP scheme. Recall that
$Y=(Y_{t_{0}},\dots,Y_{t_{N-1}})$ and $Z=(Z_{t_{0}},\dots,Z_{t_{N-1}})$ are
the solution of BSDE (3). Let
$\mathcal{U}=(\mathcal{U}_{0},\dots,\mathcal{U}_{N-1})$ and
$\mathcal{Z}=(\mathcal{Z}_{0},\dots,\mathcal{Z}_{N-1})$ be two lists of neural
networks. Define the error between $(\mathcal{U},\mathcal{Z})$ and $(Y,Z)$ as

$\displaystyle\mathcal{E}[(\mathcal{U},\mathcal{Z}),(Y,Z)]=$
$\displaystyle\max_{i=0,\dots,N-1}\mathbb{E}\left|Y_{t_{i}}-\mathcal{U}_{i}(X_{t_{i}})\right|^{2}$
(14)
$\displaystyle+\mathbb{E}\left[\sum_{i=0}^{N-1}\int_{t_{i}}^{t_{i+1}}\left|Z_{t}-\mathcal{Z}_{i}(X_{t_{i}})\right|^{2}\mathrm{d}t\right].$

We make the following usual assumptions on the BSDE (3).

######  Assumption 1.

Assume that $f,g,M,\Sigma$ satisfy the lipschitz condition, i.e., there are
constants $C_{f},C_{g},C>0$, such that for all $(t_{1},x_{1},y_{1},z_{1})$ and
$(t_{2},x_{2},y_{2},z_{2})\in[0,T]\times\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}$,

$\left\\{\begin{aligned}
|f(t_{1},x_{1},y_{1},z_{1})-&f(t_{2},x_{2},y_{2},z_{2})|\\\
\leq&C_{f}\left(|t_{1}-t_{2}|^{1/2}+\|x_{1}-x_{2}\|_{2}+|y_{1}-y_{2}|+\|z_{1}-z_{2}\|_{2}\right),\\\
|g(x_{1})-g(x_{2})|\leq&C_{g}\|x_{1}-x_{2}\|_{2},\end{aligned}\right.$

and

$\displaystyle|M(t_{1},x_{1})-M(t_{2},x_{2})|+|\Sigma(t_{1},x_{1})-\Sigma(t_{2},x_{2})|\leq
C\left(|t_{1}-t_{2}|^{1/2}+\|x_{1}-x_{2}\|_{2}\right).$

######  Assumption 2.

Assume that the approximation $\\{X_{t_{i}}\\}$ of $\\{\mathcal{X}_{t_{i}}\\}$
satisfies

$\max_{i=0,\dots,N-1}\mathbb{E}\left[|X_{t_{i+1}}-\mathcal{X}_{t_{i+1}}|^{2}\right]=O(|\pi|),$

where $|\pi|=\max_{i=0,\dots,N-1}|t_{i+1}-t_{i}|$.

The following theorem can be proved by using the method in Theorem 4.1 of Hure
et al. [1].

######  Theorem 1.

Under Assumptions 1 and 2, for functions
$\mathcal{U}=(\mathcal{U}_{0},\dots,\mathcal{U}_{N-1})$ and
$\mathcal{Z}=(\mathcal{Z}_{0},\dots,\mathcal{Z}_{N-1})$, there is a constant
$C$, depending only on the coefficients of BSDE (3), such that the error (14)
has the following upper bound

$\displaystyle\mathcal{E}\left[\left(\mathcal{U},\mathcal{Z}\right),\left(Y,Z\right)\right]\leq
C\bigg{(}|\pi|+N\sum_{i=0}^{N-1}\Big{(}L_{i}\left(\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i}\right)-L_{i}\left(\mathcal{U}_{i+1},\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*}\right)\Big{)}\bigg{)}.$
(15)

where $\mathcal{U}_{N}=g$ and
$\left(\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*}\right)=\arg\min_{(\mathcal{U}_{i},\mathcal{Z}_{i})}L_{i}\left(\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i}\right)$.

###### Proof.

See A. ∎

If we fix a time partition $\pi$, then the total error is determined by the
second term on the right hand side of (15). Our goal is to find $\mathcal{U}$
and $\mathcal{Z}$ such that each term
$L_{i}(\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i})-L_{i}(\mathcal{U}_{i+1},\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*})$
is as small as possible, enabling $(\mathcal{U},\mathcal{Z})$ to approximate
$(Y,Z)$ with the minimal error. Suppose that $\mathcal{U}_{i}$ and
$\mathcal{Z}_{i}$ are neural networks with parameter $\theta_{i}\in\Theta$. We
have the following corollary, which can be proved by the method in the A.

######  Corollary 1.

Under the conditions of Theorem 1, suppose that
$\mathcal{U}_{i}=\mathcal{U}_{i}(\cdot;\theta_{i})$ and
$\mathcal{Z}_{i}=\mathcal{Z}_{i}(\cdot;\theta_{i})$, then

$\displaystyle L_{i}($
$\displaystyle\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i})-L_{i}(\mathcal{U}_{i+1},\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*})$
(16) $\displaystyle\leq$
$\displaystyle\underbrace{\inf_{\theta\in\Theta}\mathbb{E}\left[\left|\mathcal{U}_{i}^{*}(X_{t_{i}})-\mathcal{U}_{i}(X_{t_{i}};\theta)\right|^{2}\right]+\inf_{\theta\in\Theta}\mathbb{E}\left[\left|\mathcal{Z}_{i}^{*}(X_{t_{i}})-\mathcal{Z}_{i}(X_{t_{i}};\theta)\right|^{2}\right]}_{\text{Approximation
Error}}$
$\displaystyle+\mathcal{L}_{i}(\theta_{i})-\mathcal{L}_{i}(\theta_{i}^{*}),$

where $\mathcal{U}_{i}^{*}$ and $\mathcal{Z}_{i}^{*}$ are two functions in
Theorem 1 and

$\mathcal{L}_{i}(\theta_{i}):=L_{i}(\mathcal{U}_{i+1}(\cdot;\theta_{i+1}),\mathcal{U}_{i}(\cdot;\theta_{i}),\mathcal{Z}_{i}(\cdot;\theta_{i})),\quad\theta^{*}_{i}\in\arg\min_{\theta_{i}\in\Theta}\mathcal{L}_{i}(\theta_{i}).$

The first two terms on the right hand side of (16) have the same form as in
Theorem 4.1 of Hure et al. [1]. We refer to the sum of these two terms as the
approximation error since it quantifies how accurately the functions can be
approximated by the neural networks. The approximation error can be
arbitrarily small if the network and $\Theta$ are large enough, see [32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42]. The last term quantifies how well we can
optimize the loss function to obtain the optimal parameters $\theta_{i}^{*}$.
This term vanishes if we choose $\theta_{i}$ to be $\theta^{*}_{i}$ for
$i=N-1,\dots,0$ (this is what Hure et al. [1] do).

However, we cannot obtain $\theta^{*}_{i}$ by minimizing
$\mathcal{L}_{i}(\theta_{i})$ directly, since the expectation in (12) is
usually intractable. There are various methods to optimize the loss functional
formulated as an expectation, including the Stochastic Gradient Descent (SDG)
method [43, 44, 45] and the Empirical risk minimization (ERM) method [46, 47,
48]. In those methods, the true risk $\mathcal{L}_{i}$ is usually approximated
by the empirical risk $\mathcal{L}_{i}^{(m)}$, which is exactly the sample
mean over $m$ quadrature points. Typically, we use the MC method to simulate
the points, but the QMC method (see subsection 4.3.1) is an alternative
quadrature method which achieves smaller errors when approximating
expectations. From this prospective, the optimal tractable parameter is the
minimizer of the empirical risk $\mathcal{L}_{i}^{(m)}$, which can be obtained
by ERM. Denote the minimizer of the empirical risk computed by MC and QMC
methods by $\theta_{i}^{(m)}$ and $\widehat{\theta}_{i}^{(m)}$, respectively.
To study the impact of the sampling methods, we investigate the error between
$\theta_{i}^{(m)}$ (or $\widehat{\theta}_{i}^{(m)}$) and the global optimal
parameter $\theta^{*}_{i}$, which is defined as the following (under MC, for
illustration),

$\mathcal{L}_{i}(\theta_{i}^{(m)})-\mathcal{L}_{i}(\theta_{i}^{*}).$ (17)

We follow [19] to refer to (17) as the generalization error. The
generalization error is introduced by substituting the true risk with the
empirical risk, and it is the only component of the total error which depends
on the sampling method. It is vital to derive the convergence rate of
generalization error with respect to the sample size $m$.

Compared to the classical MC method, the QMC method provides more accurate
approximations of the true risk as the sample size increases. Consequently, it
is important to study whether QMC is a better choice for our problem. To this
end, we focus on studying the convergence rates of generalization error (17)
under MC and QMC frameworks. We also numerically compare the efficiency of MC-
based and QMC-based DBDP methods. These results are presented in Sections 4
and 5.

##  4 The analysis of generalization error

###  4.1 Problem specification

We consider two basic types of nonlinear parabolic partial PDEs: the nonlinear
heat PDEs and nonlinear Black-Scholes PDEs. For heat equations, the SDE is

$\mathcal{X}_{t}=\eta+\overline{\mu}t+\overline{\sigma}W_{t},\;\;\;0\leq t\leq
T,$ (18)

where $\overline{\mu}\in\mathbb{R}^{d},\overline{\sigma}\in\mathbb{R}^{d\times
d}$. For Black-Scholes equations, the SDE is

$\mathcal{X}_{t}=\eta+\int_{0}^{t}\mathrm{diag}(\mathcal{X}_{s})\overline{\mu}\mathrm{d}s+\int_{0}^{t}\mathrm{diag}(\mathcal{X}_{s})\overline{\sigma}\mathrm{d}W_{s},\;\;\;0\leq
t\leq T,$ (19)

where $\mathrm{diag}(x)$ denotes the diagonal matrix with its diagonal
elements equal to $x$. The solution of (19) is given by

$\displaystyle\mathcal{X}_{t,i}=\eta^{(i)}\exp\left(\left(\mu_{i}-\frac{1}{2}\left\|\sigma_{i}\right\|_{2}^{2}\right)t+\left<\sigma_{i},W_{t}\right>\right),\
\ 1\leq i\leq d,$ (20)

where $\mathcal{X}_{t,i}$, $\mu_{i}$ and $\eta^{(i)}$ denote the $i$-th
components of $\mathcal{X}_{t}$, $\overline{\mu}$ and $\eta$, respectively,
and $\sigma_{i}$ denotes the $i$-th row of $\overline{\sigma}$. For
simplicity, we take $\overline{\mu}=\mu\bm{1}_{d}$ and
$\overline{\sigma}=\sigma\bm{I}_{d}$, where $\mu,\sigma\in\mathbb{R}$,
$\bm{1}_{d}$ is the d-dimensional vector where all elements equal 1 and
$\bm{I}_{d}$ is the d-dimensional identity matrix. All results in Section 4
can be extended to more general $\overline{\mu}$ and $\overline{\sigma}$.

For both nonlinear heat PDEs and nonlinear Black-Scholes PDEs, we can simulate
$\mathcal{X}_{t_{i}}$ without using Euler-Maruyama method. More precisely, for
the nonlinear heat PDEs, we use

$\displaystyle X_{t_{i+1}}$ $\displaystyle=X_{t_{i}}+\mu\Delta
t_{i}\bm{1}_{d}+\sigma\Delta W_{t_{i}}$ (21) $\displaystyle X_{t_{i}}$
$\displaystyle=\eta+\mu t_{i}\bm{1}_{d}+\sigma W_{t_{i}},$

and for nonlinear Black-Scholes PDEs, we use

$\displaystyle X_{t_{i+1},k}$
$\displaystyle=X_{t_{i},k}\exp\\{(\mu-\sigma^{2}/2)\Delta t_{i}+\sigma\Delta
W_{t_{i},k}\\}$ (22) $\displaystyle X_{t_{i},k}$
$\displaystyle=\eta_{k}\exp\\{(\mu-\sigma^{2}/2)t_{i}+\sigma W_{t_{i},k}\\},$

where $X_{t_{i},k}$, $\eta_{k}$ and $\Delta W_{t_{i},k}$ denote the $k$-th
components of $X_{t_{i}}$, $\eta$ and $\Delta W_{t_{i}}$, respectively. Note
that $X_{t_{i}}=\mathcal{X}_{t_{i}}$. Therefore, Assumption 2 holds. Since
$\eta$, $W_{t_{i}}$ and $\Delta W_{t_{i}}$ are independent random variables,
we can generate $(\eta$, $W_{t_{i}},\Delta W_{t_{i}})$ using a
$3d$-dimensional standard normal random variable $W$, i.e.,

$\displaystyle\eta=(b-a)\Phi(W_{1:d})+a\ \
,W_{t_{i}}=\sqrt{t_{i}}W_{(d+1):2d},\ \ \Delta W_{t_{i}}=\sqrt{\Delta
t_{i}}W_{(2d+1):3d},$ (23)

where $\Phi(\cdot)$ is the cumulative distribution function of the the
standard Gaussian distribution, acting on each component of $W_{1:d}$, and
$x_{i:j}$ denotes $(x_{i},\cdots,x_{j})$. The loss function (13) can be
written as

$\mathcal{L}_{i}(\theta)=\mathbb{E}\left[G_{i}(W;\theta)\right],$ (24)

where $G_{i}(\cdot;\theta)$ can be derived from (13) and (23). The following
lemma gives the upper bound and the lipschitz property of $G_{i}(W;\theta)$.

######  Lemma 2.

Denote $\mathcal{S}_{1}=(d,l_{2},\cdots,l_{L},1)$ and
$\mathcal{S}_{2}=(d,l_{2},\cdots,l_{L},d)$. Suppose that the neural networks
$\mathcal{U}_{i}\in\mathcal{N}_{\mathcal{S}_{1},R}$ and
$\mathcal{Z}_{i}\in\mathcal{N}_{\mathcal{S}_{2},R}$, $0\leq i\leq N-1$. For
nonlinear heat PDEs and nonlinear Black-Scholes PDEs, there exist constants
$C_{1},C_{2}>0$ that depend only on $T,\mu,\sigma$, such that for every
$i,\theta,\theta_{1},\theta_{2}$, we have

$\displaystyle\left|G_{i}(W;\theta)\right|$ $\displaystyle\leq
C_{1}C_{fg}d^{2}B^{2}_{\mathcal{S}_{1},R}\exp\\{C_{2}\left\|W\right\|_{\infty}\\},$
$\displaystyle\left|G_{i}(W;\theta_{1})-G_{i}(W;\theta_{2})\right|$
$\displaystyle\leq
C_{1}C_{fg}d^{2}B_{\mathcal{S}_{1},R}C_{\mathcal{S}_{1},R}\left\|\theta_{1}-\theta_{2}\right\|_{\infty}\exp\\{C_{2}\left\|W\right\|_{\infty}\\},$

where $C_{fg}=\max\\{C_{f}^{2},C_{g}^{2},f^{2}(0),g^{2}(0)\\}$,
$C_{\mathcal{S}_{1},R}=|\mathcal{S}_{1}|R^{D(\mathcal{S}_{1})-1}\left\|\mathcal{S}_{1}\right\|_{\infty}^{D(\mathcal{S}_{1})-1}$
and $B_{\mathcal{S}_{1},R}=2R\left\|\mathcal{S}_{1}\right\|_{\infty}$ (see
Definition 1 and Lemma 1).

###### Proof.

Denote

$H(X_{t_{i}},X_{t_{i+1}},\Delta
W_{t_{i}};\theta)=\widehat{\mathcal{U}}_{i+1}(X_{t_{i+1}})-F(t_{i},X_{t_{i}},\mathcal{U}_{i}(X_{t_{i}};\theta),\mathcal{Z}_{i}(X_{t_{i}};\theta),\Delta
t_{i},\Delta W_{t_{i}}).$

From (11), the lipschitz condition of functions $f,g$ and Lemma 1, we have

$\displaystyle\left|F(t,x,y,z,h,\Delta)\right|\leq\left|y\right|+\left|h\right|\left|f(0)\right|+\left|h\right|C_{f}\left(\sqrt{t}+\left|y\right|+\left\|x\right\|_{2}+\left\|z\right\|_{2}\right)+\left\|z\right\|_{\infty}\left\|\Delta\right\|_{1},$
$\displaystyle\left|F(t,x,y_{1},z_{1},h,\Delta)-F(t,x,y_{2},z_{2},h,\Delta)\right|\leq\left|y_{1}-y_{2}\right|+\left|h\right|C_{f}\left|y_{1}-y_{2}\right|$
$\displaystyle\qquad\qquad\qquad\qquad\quad+\left|h\right|C_{f}\sqrt{d}\left\|z_{1}-z_{2}\right\|_{\infty}+\left\|z_{1}-z_{2}\right\|_{\infty}\left\|\Delta\right\|_{1}.$

By applying Lemma 1, we can prove that

$\displaystyle\big{|}H(X_{t_{i}},X_{t_{i+1}},\Delta
W_{t_{i}};\theta)\big{|}\leq\left|\widehat{\mathcal{U}}_{i+1}(X_{t_{i+1}})\right|+C_{f}\left|\pi\right|\left\|X_{t_{i}}\right\|_{2}+B_{\mathcal{S}_{1},R}\left\|\Delta
W_{t_{i}}\right\|_{1}$ (25)
$\displaystyle\qquad\qquad\qquad\qquad\qquad+(1+C_{f}\left|\pi\right|+C_{f}\left|\pi\right|\sqrt{d})B_{\mathcal{S}_{1},R}+\left|f(0)\right|\left|\pi\right|+\sqrt{T}C_{f}\left|\pi\right|,$
$\displaystyle\big{|}H(X_{t_{i}},X_{t_{i+1}},\Delta
W_{t_{i}};\theta_{1})-H(X_{t_{i}},X_{t_{i+1}},\Delta
W_{t_{i}};\theta_{2})\big{|}$
$\displaystyle\leq(1+C_{f}\left|\pi\right|+C_{f}\left|\pi\right|\sqrt{d}+\left\|\Delta
W_{t_{i}}\right\|_{1})C_{\mathcal{S}_{1},R}\max\\{\left\|X_{t_{i}}\right\|_{\infty},1\\}\left\|\theta_{1}-\theta_{2}\right\|_{\infty}.$

For each $i$, we have

$\displaystyle\left|\widehat{\mathcal{U}}_{i+1}(X_{t_{i+1}})\right|$
$\displaystyle\leq\max\\{B_{\mathcal{S}_{1},R},\left|g(X_{t_{i+1}})\right|\\}$
(26)
$\displaystyle\leq\max\\{B_{\mathcal{S}_{1},R},\left|g(0)\right|+C_{g}\left\|X_{t_{i+1}}\right\|_{2}\\}.$

For the simulation scheme (21) and (23), we can easily prove that

$\displaystyle\left\|X_{t_{i}}\right\|_{\infty}$ $\displaystyle\leq
1+\left|\mu\right|T+\sigma\sqrt{T}\left\|W\right\|_{\infty},$ (27)
$\displaystyle\left\|X_{t_{i+1}}\right\|_{\infty}$ $\displaystyle\leq
1+2\left|\mu\right|T+2\sigma\sqrt{T}\left\|W\right\|_{\infty},$
$\displaystyle\left\|\Delta W_{t_{i}}\right\|_{\infty}$
$\displaystyle\leq\left|\pi\right|\left\|W\right\|_{\infty}.$

For the simulation scheme (22) and (23), we also have

$\displaystyle\left\|X_{t_{i}}\right\|_{\infty}$
$\displaystyle\leq\exp\\{\left|\mu-\sigma^{2}/2\right|T+\sigma
T\left\|W\right\|_{\infty}\\},$ (28)
$\displaystyle\left\|X_{t_{i+1}}\right\|_{\infty}$
$\displaystyle\leq\exp\\{2\left|\mu-\sigma^{2}/2\right|T+2\sigma
T\left\|W\right\|_{\infty}\\},$ $\displaystyle\left\|\Delta
W_{t_{i}}\right\|_{\infty}$
$\displaystyle\leq\left|\pi\right|\left\|W\right\|_{\infty}.$

Combining (25), (26), (27), (28) and the fact that

$G_{i}(W;\theta)=H^{2}(X_{t_{i}},X_{t_{i+1}},\Delta W_{t_{i}};\theta),$

the proof is completed after straightforward calculation. ∎

###  4.2 Generalization error for Monte Carlo methods

As shown in Section 3, the total error is bounded by the sum of three parts,
namely the scheme error, the approximation error and the generalization error.
In this section, we study the generalization error under the framework of MC.
For simplicity, we drop the subscript $i$ in (24).

Assume that there exists the optimal parameter $\theta^{*}$ such that

$\theta^{*}=\arg\min_{\theta\in\Theta}\mathcal{L}(\theta).$

Since the expectation is usually intractable, we apply the ERM method to
estimate $\theta^{*}$ by the minimizer $\theta^{(m)}$ of the empircal risk

$\mathcal{L}^{(m)}\left(\theta\right)=\frac{1}{m}\sum_{i=1}^{m}G(\bm{w}_{i};\theta)$
(29)

where $\left\\{\bm{w}_{i}:1\leq i\leq m\right\\}$ are random samples
satisfying $\bm{w}_{i}\sim W$.

We are interested in the generalization error defined by (see Theorem 1 and
Corollary 1 for details)

$\mathcal{L}(\theta^{(m)})-\mathcal{L}(\theta^{*}).$

The next lemma establishes the upper bound of the generalization error.

######  Lemma 3.

For every $m\in\mathbb{N}$, we have

$\displaystyle 0\leq\mathcal{L}(\theta^{(m)})-\mathcal{L}(\theta^{*})$
$\displaystyle\leq$ $\displaystyle
2\underset{\theta\in\Theta}{\mathrm{sup}}\left|\mathcal{L}(\theta)-\mathcal{L}^{(m)}(\theta)\right|$

###### Proof.

The lower bound follows by
$\theta^{*}=\mathrm{argmin}_{\theta\in\Theta}\mathcal{L}(\theta)$. For the
upper bounds, we have

$\displaystyle\mathcal{L}(\theta^{(m)})-\mathcal{L}(\theta^{*})$
$\displaystyle=$
$\displaystyle\mathcal{L}(\theta^{(m)})-\mathcal{L}^{(m)}(\theta^{(m)})+\mathcal{L}^{(m)}(\theta^{(m)})-\mathcal{L}^{(m)}(\theta^{*})$
$\displaystyle+\mathcal{L}^{(m)}(\theta^{*})-\mathcal{L}(\theta^{*})$
$\displaystyle\leq$
$\displaystyle\mathcal{L}(\theta^{(m)})-\mathcal{L}^{(m)}(\theta^{(m)})+\mathcal{L}^{(m)}(\theta^{(m)})-\mathcal{L}^{(m)}(\theta^{*})$
$\displaystyle\leq$ $\displaystyle
2\underset{\theta\in\Theta}{\mathrm{sup}}\left|\mathcal{L}(\theta)-\mathcal{L}^{(m)}(\theta)\right|,$

where the second inequality follows by
$\mathcal{L}^{(m)}(\theta^{(m)})\leq\mathcal{L}^{(m)}(\theta^{*})$. ∎

####  4.2.1 The convergence rate of mean error for MC

Under the MC framework, $(\bm{w}_{i})_{i=1}^{m}$ in (29) are identical and
independent standard Gaussian random variables. By Lemma 3, the mean
generalization error satisfies

$\mathbb{E}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}(\theta^{*})\right]\leq
2\mathbb{E}\left[\sup_{\theta\in\Theta}\left|\frac{1}{m}\sum_{i=1}^{m}G(\bm{w}_{i};\theta)-\mathbb{E}\left[G(W;\theta)\right]\right|\right],$
(30)

where $W$ is a standard Gaussian random variable. As a result, it suffices to
obtain the upper bound of the right hand side of (30). To this end, we use the
Rademacher complexity technique.

######  Definition 2 (Rademacher complexity).

Suppose that $(\bm{x}_{i})_{i=1}^{m}$ are i.i.d. random variables on
$(\Omega,\mathcal{F},\mathbb{P})$ and take values in a metric space $S$.
Suppose $\mathcal{G}$ is a class of real-valued borel-measurable function over
$(S,\mathcal{B}(S))$ ($\mathcal{B}(S)$ is the Borel $\sigma$-field). The
Rademacher complexity of a function class $\mathcal{G}$ with respect to the
random sample $(\bm{x}_{i})_{i=1}^{m}$ is defined by

$R_{m}(\mathcal{\mathcal{G}};(\bm{x}_{i})_{i=1}^{m})=\mathbb{E}\left[{\underset{g\in\mathcal{G}}{\mathrm{sup}}\frac{1}{n}\sum_{i=1}^{m}\xi_{i}g(\bm{x}_{i})}\right],$

where $\left\\{\xi_{i}\right\\}_{i=1}^{m}$ are i.i.d. Rademacher variables
with discrete distribution
$\mathbb{P}\\{\xi_{i}=1\\}=\mathbb{P}\\{\xi_{i}=-1\\}=1/2$.

The next lemma provides an upper bound for the right hand side of (30) in
terms of the Rademacher complexity. We refer to [47, Lemma 26.2] and [48,
Theorem 3.3] for a detailed proof.

######  Lemma 4.

Under the setting of Definition 2, we have

$\mathbb{E}\left[{\underset{g\in\mathcal{G}}{\mathrm{sup}}\left|\frac{1}{m}\sum_{i=1}^{m}g(\bm{x}_{i})-\mathbb{E}[g(\bm{x})]\right|}\right]\leq
2R_{m}(\mathcal{\mathcal{G}};(\bm{x}_{i})_{i=1}^{m}).$

The Rademacher complexity has the following upper bound.

######  Lemma 5.

Under the setting of Definition 2, suppose $\mathcal{G}$ is a parameterized
class of functions
$\left\\{\mathcal{F}(\cdot;\theta):\mathcal{S}\to\mathbb{R}\mid\theta\in\Theta\right\\}$,
where $\Theta\subseteq\mathbb{R}^{p}$ is the parametric space satisfying
$\left\|\theta\right\|_{\infty}\leq R$. Suppose that there exist constants
$B_{1}$ and $B_{2}$ such that for all $\theta_{1},\theta_{2}\in\Theta$ it
holds that

$\displaystyle\underset{x\in
S}{\mathrm{sup}}\left|\mathcal{F}(x;\theta_{1})-\mathcal{F}(x;\theta_{2})\right|$
$\displaystyle\leq$ $\displaystyle
B_{1}\left\|\theta_{1}-\theta_{2}\right\|_{\infty},$
$\displaystyle\underset{x\in
S}{\mathrm{sup}}\left|\mathcal{F}(x;\theta_{1})\right|$ $\displaystyle\leq$
$\displaystyle B_{2}.$

Then we have

$\displaystyle R_{m}(\mathcal{G};(x_{i})_{i=1}^{m})$ $\displaystyle\leq$
$\displaystyle\frac{4}{\sqrt{m}}+\frac{6\sqrt{p}B_{2}}{\sqrt{m}}\sqrt{\log\left(2B_{1}R\sqrt{m}\right)}.$

For a detailed proof of Lemma 5, we refer to [17, Theorem 5.1]. By Lemma 4,
the right hand side of (30) is bounded by the Rademacher complexity of
function class $\\{G(\cdot;\theta)\\}_{\theta\in\Theta}$. However, we can not
directly use Lemma 5 to obtain the upper bound of this Rademacher complexity,
since $G(\cdot;\theta)$ is not bounded. To address this issue, we use the
truncation method to obtain the mean generalization error rate of MC.

######  Theorem 2.

In the cases of nonlinear heat PDEs and nonlinear Black-Sholes PDEs, suppose
that we use the MC method to obtain the samples $\left\\{\bm{w}_{i}:1\leq
i\leq m\right\\}$ in (29), then we have

$\displaystyle\mathbb{E}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}(\theta^{*})\right]=Cm^{-1/2}\log(m),$

where $C$ is a polynomial of
$(R,\|\mathcal{S}_{1}\|_{\infty},D(\mathcal{S}_{1}),C_{fg})$, independent of
$m$.

###### Proof.

For every $K>0$, denote $G^{(K)}(x;\theta):=\bm{1}_{\\{\|x\|_{\infty}\leq
K\\}}G(x;\theta)$ be the truncated version of $G(x;\theta)$. The truncated
function class $\mathcal{G}^{(K)}:=\\{G^{(K)}(x;\theta)\\}_{\theta\in\Theta}$
satisfies the conditions of Lemma 5. It follows from (30) that

$\displaystyle\mathbb{E}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}(\theta^{*})\right]\leq$
$\displaystyle
2\underset{\theta\in\Theta}{\mathrm{sup}}\left|\mathbb{E}\left[\bm{1}_{\\{\|W\|_{\infty}>K\\}}G(W;\theta)\right]\right|$
$\displaystyle+2\mathbb{E}\left[\underset{\theta\in\Theta}{\mathrm{sup}}\left|\bm{1}_{\\{\|W\|_{\infty}>K\\}}G(W;\theta)\right|\right]+4R_{m}(\mathcal{G}^{(K)};(\bm{w}_{i})_{i=1}^{m})$
$\displaystyle\leq$ $\displaystyle
4\mathbb{E}\left[\underset{\theta\in\Theta}{\mathrm{sup}}\left|\bm{1}_{\\{\|W\|_{\infty}>K\\}}G(W;\theta)\right|\right]+4R_{m}(\mathcal{G}^{(K)};(\bm{w}_{i})_{i=1}^{m}).$
(31)

Observed by Lemma 2, we have

$\displaystyle\left|G(W;\theta)\right|$ $\displaystyle\leq
C_{1}C_{fg}d^{2}B^{2}_{\mathcal{S}_{1},R}\exp\\{C_{2}\|W\|_{\infty}\\},$ (32)
$\displaystyle\left|G(W;\theta_{1})-G(W;\theta_{2})\right|$ $\displaystyle\leq
C_{1}C_{fg}d^{2}B_{\mathcal{S}_{1},R}L_{\mathcal{S}_{1},R}\exp\\{C_{2}\|W\|_{\infty}\\}\left\|\theta_{1}-\theta_{2}\right\|_{\infty}.$
(33)

For the first term on the right hand side of (31), it follows from (32) and
Hölder’s inequality that

$\displaystyle\mathbb{E}\left[\underset{\theta\in\Theta}{\mathrm{sup}}\left|1_{\\{\|W\|_{\infty}>K\\}}G(W;\theta)\right|\right]$
$\displaystyle\leq\sqrt{\mathbb{E}\left[\bm{1}_{\\{\|W\|_{\infty}>K\\}}\right]}\sqrt{\mathbb{E}\left[\sup_{\theta\in\Theta}|G(W;\theta)|^{2}\right]}$
(34) $\displaystyle\leq
C_{1}C_{fg}d^{3}B^{2}_{\mathcal{S}_{1},R}\sqrt{d\mathbb{E}\left[\exp\\{C_{2}|W_{1}|\\}\right]}\sqrt{\mathbb{P}\left[\|W\|_{\infty}>K\right]}$
$\displaystyle\leq\sqrt{2}C_{1}C_{fg}d^{4}B^{2}_{\mathcal{S}_{1},R}\sqrt{\mathbb{E}\left[\exp\\{C_{2}|W_{1}|\\}\right]}\exp\\{-K^{2}/4\\},$

where in the second inequality, we use the facts that

$\mathbb{E}\left[\exp\\{C_{2}\|W\|_{\infty}\\}\right]\leq
d\mathbb{E}\left[\exp\\{C_{2}|W_{1}|\\}\right],$

and

$\mathbb{P}\left[\left\|W\right\|_{\infty}>K\right]\leq\sum_{i=1}^{d}\mathbb{P}\left[\left|W_{i}\right|>K\right]\leq
2d\exp\\{-K^{2}/2\\}.$

For the second term on the right hand side of (31), it follows from (33) and
Lemma 5 that

$\displaystyle
R_{m}(\mathcal{G}^{(K)};(\bm{w}_{i})_{i=1}^{m})\leq\frac{4}{\sqrt{m}}+\frac{6\sqrt{|\mathcal{S}_{1}|+|\mathcal{S}_{2}|}B_{2}(K)}{\sqrt{m}}\sqrt{\log\left(2B_{1}(K)R\sqrt{m}\right)},$
(35)

where
$B_{1}(K)=C_{1}C_{fg}d^{2}B_{\mathcal{S}_{1},R}C_{\mathcal{S}_{1},R}\exp\\{C_{2}K\\}$,
$B_{2}(K)=C_{1}C_{fg}d^{2}B^{2}_{\mathcal{S}_{1},R}\exp\\{C_{2}K\\}.$

We take $K=C_{2}^{-1}\log(m)$. The upper bound of (34) becomes

$h_{1}\exp\left\\{-(C_{2}^{-1}\log(m))^{2}/4\right\\}\leq\exp\left\\{C_{2}^{2}\right\\}h_{1}m^{-1/2},$
(36)

where $h_{1}$ is a polynomial of $(d,R,\|\mathcal{S}_{1}\|_{\infty},C_{fg})$
and we use the fact that $(C_{2}^{-1}\log(m))^{2}+C_{2}^{2}\geq 2\log(m)$. The
upper bound of (35) becomes

$h_{2}m^{-1/2}\sqrt{\log(m\sqrt{m})}\leq h_{2}m^{-1/2}\log(m),$ (37)

where $h_{2}$ is a polynomial of
$(d,R,\|\mathcal{S}_{1}\|_{\infty},D(\mathcal{S}_{1}),C_{fg})$. Combining
(34), (36), (35) and (37), we prove the desired result. ∎

If we use the MC method to obtain the samples for computing the empirical
risk, Theorem 2 demonstrate that the generalization error converges at the
rate of $O(\log(m)m^{-1/2})$ with respect to the sample size $m$. This rate is
consistent with those found in classical learning theory [47, 48], where the
output is usually assumed to be bounded.

####  4.2.2 Break the curse of dimensionality

We denote
$\mathcal{L}_{i}^{*}=L_{i}(\mathcal{U}_{i+1},\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*})$.
For simplicity, we drop the subscript $i$ in this section. In Theorem 2, we
find that the implied constant is actually a polynomial of the constant
$C_{fg}$ (see Lemma 2) and the size of the neural networks, including the
parameter bound $R$, the width $\|\mathcal{S}_{1}\|_{\infty}$, the depth
$D(\mathcal{S}_{1})$. This allows us to prove that the ERM method breaks the
curse of dimensionality, provided that the approximation error in Corollary 1
satisfies certain assumptions.

######  Assumption 3.

Let $\mathcal{S}_{1}=(d,l_{2},\cdots,l_{L},1)$,
$\mathcal{S}_{2}=(d,l_{2},\cdots,l_{L},d)$. For each $i\in\\{1,2,\cdots,N\\}$,
there exists a poynomial $p$ such that for for every $\delta>0$, there exists
a neural network $\mathcal{U}_{\delta}\in\mathcal{N}_{\mathcal{S}_{1},R}$ and
$\mathcal{Z}_{\delta}\in\mathcal{N}_{\mathcal{S}_{1},R}$ such that

$\mathbb{E}\left[\left|\mathcal{U}_{i}^{*}(X_{t_{i}})-\mathcal{U}_{\delta}(X_{t_{i}})\right|^{2}\right]+\mathbb{E}\left[\left|\mathcal{Z}_{i}^{*}(X_{t_{i}})-\mathcal{Z}_{\delta}(X_{t_{i}})\right|^{2}\right]\leq\delta,$
(38)

and

$\max\\{R,D(\mathcal{S}_{1}),\left\|\mathcal{S}_{1}\right\|_{\infty}\\}\leq
p(d,\delta^{-1}).$ (39)

######  Theorem 3.

Under Assumption 3, suppose that the PDE (1) is nonlinear heat PDEs or
nonlinear Black-Sholes PDEs described in Section 4.1. Let
$C_{fg}=\max\\{C_{f}^{2},C_{g}^{2},f^{2}(0),g^{2}(0)\\}$
,$\mathcal{S}_{1}=(d,l_{2},\cdots,l_{L},1)$ and
$\mathcal{S}_{2}=(d,l_{2},\cdots,l_{L},d)$. Then there exists a polynomial $h$
such that for $m$ satisfying

$m\geq h\left(\delta^{-1},d,C_{fg}\right),$

we have

$\mathbb{E}[\mathcal{L}(\theta^{(m)})-\mathcal{L}^{*}]\leq\delta.$ (40)

Moreover, the coefficients of $h$ only depends on the constants $C_{1},C_{2}$
(see Lemma 2).

###### Proof.

This is a direct consequence of Corollary 1 and Theorem 2. ∎

Studying the probability inequality of the excess risk
$\mathcal{L}(\theta^{(m)})-\mathcal{L}^{*}$ is an important topic in the
analysis of the ERM method, see [47, 48, 20, 17, 49]. In some contexts, the
problem is for $\delta,\rho>0$, determining the $m$ such that

$\mathbb{P}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}^{*}\leq\delta\right]\geq
1-\rho.$ (41)

It follows from Markov inequality that for $\delta>0$,

$\mathbb{P}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}^{*}\leq\delta\right]\geq
1-\frac{\mathbb{E}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}^{*}\right]}{\delta}.$

Thus we have the natural extension of Theorem 3.

######  Corollary 2.

Under the setting of Theorem 3. There exists a polynomial $q$ such that for
$m$ satisfying

$m\geq q\left(\delta^{-1},\rho^{-1},d,C_{fg}\right),$

we have

$\mathbb{P}\left[\mathcal{L}(\theta^{(m)})-\mathcal{L}^{*}\leq\delta\right]\geq
1-\rho.$ (42)

Moreover, the coefficients of $q$ only depends on the constants $C_{1},C_{2}$
(see Lemma 2).

In fact, Theorem 3 and Corollary 2 show that, if the functions
$\mathcal{U}^{*}$ and $\mathcal{Z}^{*}$ can be approximated by the neural
networks without the curse of dimensionality, then the required sample size
$m$ to satisfy the probabilistic inequality grows polynomially, which means
the total error of the ERM principle over the neural networks overcomes the
curse of dimensionality as dimension $d$ increases.

There are numerous theoretical results on the approximation properties of the
neural networks. For universal approximation results, we refer the readers to
[36, 34, 35, 33, 32]. For the approximation of solutions to PDEs, various
studies have demonstrated that neural networks can overcome the curse of
dimensionality, see [37, 38, 39, 40, 41, 42]. Verifying the Assumption 3 is
challenging in general cases since the target functions are represented by
expectations. Additionally, it remains an open problem to prove that the
neural networks can overcome the curse of dimensionality when approximating
the solutions to the general forms of nonlinear heat equations and Black-
Scholes equations.

###  4.3 Generalization error for quasi-Monte Carlo

In this section, we consider the QMC method and study the corresponding
generalization error rate.

####  4.3.1 Quasi-Monte Carlo methods

QMC [50] is an efficient quadrature method to numerically solve the integral
problems on the unit cube $[0,1]^{d}$. While the MC method, relying on random
sampling, achieves a convergence rate of $O(m^{-1/2})$ with $m$ quadrature
points, QMC leverages low discrepancy sequences, such as the Halton sequence,
Faure sequence and Sobol’ sequence [51, 50, 52], which facilitates a faster
convergence rate of $O(m^{-1+\varepsilon})$ for an arbitrarily small
$\varepsilon>0$. This enhanced convergence rate of QMC is grounded in the
Koksma-Hlawka inequality [53], linking the error to both the variation of the
integrand in the sense of Hardy and Krause and the uniformity of sample points
as measured by star discrepancy. More precisely, the Koksma-Hlawka inequality
states that

$\left|\frac{1}{m}\sum_{i=1}^{m}f(\bm{y}_{i})-\int_{[0,1]^{d}}f(\bm{y})\mathrm{d}\bm{y}\right|\leq
V_{\mathrm{HK}}(f)D_{m}^{*}(\\{\bm{y}_{1},\dots,\bm{y}_{m}\\}),$ (43)

where $V_{\mathrm{HK}}(f)$ is the variation of $f$ in the sense of Hardy and
Krause and $D_{m}^{*}$ is the star discrepancy of point set
$\\{\bm{y}_{1},\dots,\bm{y}_{m}\\}$ with the magnitude $O(m^{-1}(\log
m)^{d-1})$ for low discrepancy point sets (see [50] for more details). The
convergence rate of QMC is validated if the function $f$ is of bounded
variation in the sense of Hardy and Krause (BVHK).

The randomized quasi-Monte Carlo (RQMC) method, which includes random shifts
[54] and scrambling [55], merges the benefits of randomization with low
discrepancy sequences to obtain unbiased estimations. If the integrand is
bounded and smooth, the nested scrambled method can achieve a faster
convergence rate of $O(n^{-3/2+\varepsilon})$, as detailed by [56]. Recently,
Ouyang et al. [57] have demonstrated that this high convergence rate can also
be maintained for unbounded and smooth integrands by employing appropriate
importance sampling methods.

####  4.3.2 The convergence rate of mean error for QMC

Combining QMC methods, the empirical risk is defined by

$\widehat{\mathcal{L}}^{(m)}(\theta):=\frac{1}{m}\sum_{i=1}^{m}G(\Phi^{-1}(\bm{y}_{i});\theta),$
(44)

where $\\{\bm{y}_{1},\dots,\bm{y}_{m}\\}$ is a low discrepancy point set and
$\Phi^{-1}$ is the inverse cumulative distribution function of standard
Gaussian distribution, acting on each component of $\bm{y}_{i}$. In this
paper, we consider the RQMC method with the scrambling approach (see Owen [55]
for more details). Define

$\widehat{\theta}^{(m)}:=\arg\min_{\theta\in\Theta}\widehat{\mathcal{L}}^{(m)}(\theta).$
(45)

By Lemma 3, the mean generalization error satisfies

$\displaystyle\mathbb{E}\left[\mathcal{L}(\widehat{\theta}^{(m)})-\mathcal{L}(\theta^{*})\right]$
$\displaystyle\leq
2\mathbb{E}\left[\sup_{\theta\in\Theta}\Big{|}\mathcal{L}(\theta)-\widehat{\mathcal{L}}_{m}(\theta)\Big{|}\right]$
(46)
$\displaystyle=2\mathbb{E}\left[\sup_{\theta\in\Theta}\Big{|}\mathbb{E}\left[G(W;\theta)\right]-\sum_{i=1}^{m}G(\Phi^{-1}(\bm{y}_{i});\theta)\Big{|}\right].$

Since $G(\Phi^{-1}(\cdot);\theta)$ is an unbounded function (it can not be of
BVHK), the Koksma-Hlawka inequality can not be directly used to provide a
convergence rate for the right hand side of (46). To address this issue, Xiao
et al. [29] applied the results in Owen [58] and obtained the convergence rate
for the case of linear Kolmogorov PDEs. Utilizing Owen [58] method would
become quite cumbersome for our case, given that the function $G$ is more
complex than those in Xiao et al [29]. Here, we recommend employing the
results in Ouyang et al. [57], where we only need to examine the growth
conditions of $G$. We introduce the following function class from [57].

######  Definition 3.

For $h:\mathbb{R}^{3d}\rightarrow\mathbb{R}$, define the exponential growth
function class as

$\mathcal{G}_{e}(A,B):=\left\\{h\in\mathcal{S}^{3d}(\mathbb{R}^{3d}):\sup_{\bm{u}\subset
1:3d}|\partial^{\bm{u}}h(\bm{x})|\leq Be^{A|\bm{x}|}\right\\},$

where $\bm{u}\subset 1:3d=\left\\{1,\dots,3d\right\\}$ and
$\partial^{\bm{u}}:=\prod_{i\in\bm{u}}\partial/\partial x_{i}$,
$\mathcal{S}^{3d}(\mathbb{R}^{3d})$ is the function class where the functions
have continuous mixed partial derivatives of order $1$.

The following lemma can be obtained from the proof of Corollary 4.10 in Ouyang
et al. [57], which gives the convergence rate of QMC methods regarding the
function class $\mathcal{G}_{e}(A,B)$.

######  Lemma 6.

Let $\left\\{\bm{y}_{1},\dots,\bm{y}_{m}\right\\}$ be a RQMC point set such
that each $\bm{y}_{j}$ is uniformly distributed on $[0,1]^{3d}$ and

$D^{*}_{m}(\left\\{\bm{y}_{1},\dots,\bm{y}_{m}\right\\})\leq C\frac{(\log
m)^{d-1}}{m}\quad a.s.,$ (47)

where $D^{*}_{m}$ is the star discrepancy and $C$ is a constant independent of
$m$. Then we have

$\mathbb{E}\left[\sup_{h\in\mathcal{G}_{e}(A,B)}\Big{|}\sum_{i=1}^{m}h(\Phi^{-1}(\bm{y}_{i}))-\mathbb{E}\left[h(W)\right]\Big{|}\right]=O(m^{-1+\varepsilon}),$
(48)

where $W$ is a standard $3d$-dimensional Gaussian random variable.

Noting that the nested scrambled Sobol’ point set satisfies the conditions in
Lemma 6, if we use this QMC point set $\\{\bm{y}_{i}\\}$ in (44), then by
Lemma 6 and (46), it suffices to verify that the function $G(\cdot;\theta)$ is
in the class $\mathcal{G}_{e}(A,B)$ for some constant $A,B$. The main result
of this section is the following theorem.

######  Theorem 4.

If we use nested scrambled Sobol’ point set in (44), then for the nonlinear
heat PDEs and the nonlinear Black-Scholes PDEs, we have

$\mathbb{E}\left[\mathcal{L}(\widehat{\theta}^{(m)})-\mathcal{L}(\theta^{*})\right]=O(m^{-1+\varepsilon})$
(49)

for every $\varepsilon>0$.

###### Proof.

By the same method used in Lemma 3, we have

$0\leq\mathcal{L}(\widehat{\theta}^{(m)})-\mathcal{L}(\theta^{*})\leq
2\sup_{\theta\in\Theta}\left|\mathcal{L}(\theta)-\widehat{\mathcal{L}}_{m}(\theta)\right|.$
(50)

Thus,

$\displaystyle\mathbb{E}\left[\mathcal{L}(\widehat{\theta}^{(m)})-\mathcal{L}(\theta^{*})\right]$
$\displaystyle\leq
2\mathbb{E}\left[\sup_{\theta\in\Theta}\Big{|}\mathcal{L}(\theta)-\widehat{\mathcal{L}}_{m}(\theta)\Big{|}\right]$
(51)
$\displaystyle=2\mathbb{E}\left[\sup_{\theta\in\Theta}\Big{|}\mathbb{E}\left[G(W;\theta)\right]-\sum_{i=1}^{m}G(\Phi^{-1}(\bm{y}_{i});\theta)\Big{|}\right].$

We consider the growth of $G(\cdot;\theta)$. Note that $G(\cdot;\theta)$ is
composed through addition, multiplication and composition of some base
functions (such as quadratic functions, linear functions, activation
functions), where the derivatives of any order of these functions are
dominated by a polynomial. By the results in Appendix of Ouyang et al. [57],
there are constants $A,B$, such that for any $\theta\in\Theta$,

$G(\cdot;\theta)\in\mathcal{G}_{e}(A,B).$ (52)

Note that the nested scrambled Sobol’ point set satisfies the conditions in
Lemma 6. By Lemma 6 and combining (51), we have

$\displaystyle\mathbb{E}\left[\mathcal{L}(\widehat{\theta}^{(m)})-\mathcal{L}(\theta^{*})\right]$
$\displaystyle\leq
2\mathbb{E}\left[\sup_{\theta\in\Theta}\Big{|}\mathbb{E}\left[G(W;\theta)\right]-\sum_{i=1}^{m}G(\Phi^{-1}(\bm{y}_{i});\theta)\Big{|}\right]$
$\displaystyle\leq
2\mathbb{E}\left[\sup_{h\in\mathcal{G}_{e}(A,B)}\Big{|}\sum_{i=1}^{m}h(\Phi^{-1}(\bm{y}_{i}))-\mathbb{E}\left[h(W)\right]\Big{|}\right]=O(m^{-1+\varepsilon}).$

This proves the desired result. ∎

Theorem 2 and 4 establish the convergence rates of the generalization error
under the MC and QMC methods. The results that the QMC method possesses a
higher convergence rate $O(m^{-1+\varepsilon})$, indicating that, with
sufficiently thorough training (achieving optimal parameters for each), the
QMC method can achieve a smaller generalization error compared to the MC
method. According to Theorem 1 and Corollary 1 (which deals with the
decomposition of the total error), the QMC method potentially enables the
final trained solution to have a smaller error. In the next section, we
validate our theoretical results with several numerical examples.

##  5 Numerical experiments

In this section, we conduct numerical experiments to compare the performance
of the DBDP scheme using different sample methods, specifically MC and RQMC
methods. To ensure a fair comparison, we use the same training settings for
different sampling methods, and to minimize the approximation error and
highlight the impact of the generalization error, we choose a network larger
than that in [15]. For each step, we use a feedforward neural network with a
depth of 4 and a width of $d+20$. The activation function is the Tanh function
(9), and the network is initialized by means of Xavier initialization [59].
The networks are trained using the Adam optimizer [44] with default weight
decay $0.01$. For the first step near the terminal condition, the initial
learning rate of is $0.01$ and is halved every $5000$ iterations, with a total
of $50000$ iterations. For the subsequent steps, the initial learning rate of
is $0.001$ and is halved every $500$ iterations, with a total of $5000$
iterations. We also add batch normalization layers [60] to enhance the
efficiency and robustness of the network training.

We measure the performance of the obtained solution by the relative $L_{2}$
error defined as

$\sqrt{\frac{\mathbb{E}\left[\big{|}\widehat{\mathcal{U}}(\eta)-u(0,\eta)\big{|}^{2}\right]}{\mathbb{E}[|u(0,\eta)|^{2}]}},$
(53)

where $\widehat{\mathcal{U}}$ is the trained neural network, $u(0,\cdot)$ is
the exact solution of the nonlinear parabolic PDEs (1), $\eta$ is uniformly
distributed over target domain $[a,b]^{d}$. Since the relative $L_{2}$ error
has no closed form, we approximated it by using the MC method as follows

$\sqrt{\frac{\sum_{i=1}^{m}[\widehat{\mathcal{U}}(\eta_{i})-u(0,\eta_{i})]^{2}}{\sum_{i=1}^{m}u^{2}(0,\eta_{i})}},$
(54)

where $(\eta_{i})_{i=1}^{m}$ are i.i.d. uniform samples on $[a,b]^{d}$. In the
following experiments, we take the sample size $m=2^{16}$.

Note that unlike the solution obtained by training only at a single point in
[15, 14, 1], we train solutions over an entire region. To demonstrate the
efficiency of the training results at single points, in addition to presenting
the $L_{2}$ error of the results, we also provide histograms of the single-
point errors.

###  5.1 Nonlinear heat equations

The first example is a nonlinear heat PDE from [15]. We take the drift
function $M(t,x)=\mu\bm{1}_{d}$, the diffusion function $\Sigma(t,x)=\sigma
I_{d}$, the terminal condition $g(x)=\cos(\overline{x})$ with
$\overline{x}=\sum_{i=1}^{d}x_{i}$, and the nonlinear term

$\displaystyle f(t,x,y,z)=$ $\displaystyle\left(\cos(\overline{x})+\mu
d\sin(\overline{x})\right)e^{\frac{T-t}{2}}-\frac{1}{2}\left(\sigma\sqrt{d}\sin(\overline{x})\cos(\overline{x})e^{T-t}\right)^{2}$
(55) $\displaystyle+\frac{1}{2d}(y(\bm{1}_{d}\cdot z))^{2}.$

The solution of this PDE is $u(t,x)=\cos(\overline{x})e^{\frac{T-t}{2}}$. In
our numerical experiment, we choose $\mu=0.2/d$, $\sigma=1/\sqrt{d}$ and the
region $[a,b]^{d}=[-0.5,0.5]^{d}$.

Table 1 presents the mean and standard deviation of the relative $L^{2}$ error
over 8 independent runs for solving the nonlinear heat equation (55), for
dimensions 20 and 50, respectively. The table demonstrates that the RQMC-based
algorithm outperforms the MC-based one. When the batch size is $2^{14}$, the
mean error of RQMC-based solution is only $1/2$ of the MC-based solution. The
standard deviations presented in Table 1 show that the variance of RQMC-based
solutions is much smaller than the MC ones with the variance reduction factor
(the ratio of MC variance to RQMC variance) varying from $4$ to $30$.

Figure 1 presents the density histogram of pointwise errors between the
trained networks and the exact solution across the sample
$(\eta_{i})_{i=1}^{65536}$. The batch size is set to $2^{14}$, and the neural
networks with the largest and smallest $L_{2}$ errors are selected from 8
independent runs. For the MC method, the smallest $L_{2}$ error is $1.65\times
10^{-2}$ and the largest $L_{2}$ error is $4.12\times 10^{-2}$. For the RQMC
method, the smallest $L_{2}$ error is $9.32\times 10^{-3}$ and the largest
$L_{2}$ error is $1.14\times 10^{-2}$. The figure demonstrates that the RQMC-
based solution achieves smaller pointwise errors across the majority of
samples. Figure 1 only displays the errors within the range of $[-0.1,0.1]$.
For the best trained networks using MC methods, 143 out of 65,536 errors fall
outside this range , while for the RQMC method, 0 out of 65,536 errors fall
outside this range. This indicates that the solution based on RQMC
approximates the exact solution with smaller error and smaller variance.

Dimension | Batch Size | Method | Mean Error | Standard Deviation
---|---|---|---|---
20 | $2^{12}$ | MC | $3.4\times 10^{-2}$ | $6.8\times 10^{-3}$
RQMC | $\bm{3.2\times 10^{-2}}$ | $\bm{8.7\times 10^{-4}}$
$2^{14}$ | MC | $1.5\times 10^{-2}$ | $6.0\times 10^{-3}$
RQMC | $\bm{6.9\times 10^{-3}}$ | $\bm{1.3\times 10^{-3}}$
$2^{16}$ | MC | $7.2\times 10^{-3}$ | $1.0\times 10^{-3}$
RQMC | $\bm{5.8\times 10^{-3}}$ | $\bm{3.6\times 10^{-4}}$
50 | $2^{12}$ | MC | $5.6\times 10^{-2}$ | $1.8\times 10^{-2}$
RQMC | $\bm{3.5\times 10^{-2}}$ | $\bm{3.7\times 10^{-3}}$
$2^{14}$ | MC | $2.7\times 10^{-2}$ | $7.3\times 10^{-3}$
RQMC | $\bm{1.0\times 10^{-2}}$ | $\bm{6.4\times 10^{-4}}$
$2^{16}$ | MC | $1.2\times 10^{-2}$ | $2.8\times 10^{-3}$
RQMC | $\bm{9.3\times 10^{-3}}$ | $\bm{2.6\times 10^{-3}}$
Table 1: The mean and standard deviation of the relative $L^{2}$ error over 8
independent runs for solving the nonlinear heat equation.

(a) Best trained neural networks

(b) Worst trained neural networks

Figure 1: The density histograms of the pointwise absolute errors for the best
and worst trained neural networks in solving the nonlinear heat equation. The
batch size is 16,384 and the dimension is 50.

###  5.2 Hamilton–Jacobi–Bellman (HJB) equations

The second example is an HJB equation from [16, 14]. This specific HJB
equation is given by

$\displaystyle\frac{\partial u}{\partial t}(t,x)+(\Delta_{x}u)(t,x)$
$\displaystyle=\left\|(\nabla_{x}u)(t,x)\right\|_{2}^{2},\;\;\;\;\;\;\text{
for }(t,x)\in[0,T)\times\mathbb{R}^{d},$ (56) $\displaystyle u(T,x)$
$\displaystyle=\left\|x\right\|_{2}^{1/2},\;\;\;\;\;\mbox{ for
}x\in\mathbb{R}^{d}.$

This HJB equation can be transformed into a linear kolmogorov PDE through the
Cole-Hopf transformation [61]. The solution can be expressed as

$u(t,x)=-\ln\left(\mathbb{E}\left[\exp\left(-\|x+\sqrt{2}W_{T-t}\|_{2}^{1/2}\right)\right]\right).$
(57)

Therefore, to approximate the relative $L_{2}$ error (54), we use MC method to
approximate the exact solution pointwise on samples $(\eta_{i})_{i=1}^{m}$. In
our numerical experiment, we choose the region $[a,b]^{d}=[0,1]^{d}$.

For the HJB equation (56), Table 2 presents the mean and standard deviation of
relative $L_{2}$ errors over 8 independent runs for dimensions 20 and 50,
respectively. The table reveals that RQMC results are superior to MC, with a
maximum error reduction to $1/4$ of MC error (e.g., when the dimension is 20
and the batch size is $2^{12}$). Figure 2 presents the density histogram of
pointwise errors between the trained networks and the exact solution. The
batch size is set to $2^{14}$, and the neural networks with the largest and
smallest $L_{2}$ errors are selected from 8 independent runs. For the MC
method, the smallest $L_{2}$ error is $4.47\times 10^{-4}$ and the largest
$L_{2}$ error is $9.89\times 10^{-4}$. For the RQMC method, the smallest
$L_{2}$ error is $1.05\times 10^{-4}$ and the largest $L_{2}$ error is
$6.54\times 10^{-4}$. Figure 2 indicates that the distributions of errors
roughly follow a normal distribution, with RQMC exhibiting smaller variance
and a mean closer to 0. This suggests that the solutions trained using RQMC
are closer to the true solutions with smaller errors, resulting in more stable
outcomes.

Dimension | Batch Size | Methods | Mean Error | Standard Deviation
---|---|---|---|---
20 | $2^{12}$ | MC | $1.7\times 10^{-3}$ | $8.5\times 10^{-4}$
RQMC | $\bm{4.4\times 10^{-4}}$ | $\bm{2.1\times 10^{-4}}$
$2^{14}$ | MC | $1.2\times 10^{-3}$ | $3.1\times 10^{-4}$
RQMC | $\bm{3.9\times 10^{-4}}$ | $\bm{2.0\times 10^{-4}}$
$2^{16}$ | MC | $5.4\times 10^{-4}$ | $3.3\times 10^{-4}$
RQMC | $\bm{1.9\times 10^{-4}}$ | $\bm{9.1\times 10^{-5}}$
50 | $2^{12}$ | MC | $1.3\times 10^{-3}$ | $4.4\times 10^{-4}$
RQMC | $\bm{2.8\times 10^{-4}}$ | $\bm{5.3\times 10^{-5}}$
$2^{14}$ | MC | $6.7\times 10^{-4}$ | $1.8\times 10^{-4}$
RQMC | $\bm{2.9\times 10^{-4}}$ | $\bm{1.5\times 10^{-4}}$
$2^{16}$ | MC | $3.6\times 10^{-4}$ | $1.3\times 10^{-4}$
RQMC | $\bm{1.7\times 10^{-4}}$ | $\bm{8.3\times 10^{-5}}$
Table 2: The mean and standard deviation of the relative $L^{2}$ error over 8
independent runs for solving the HJB equation.

(a) Best trained neural networks

(b) Worst trained neural networks

Figure 2: The density histograms of the pointwise absolute errors for the best
and worst trained neural networks in solving the HJB equation. The batch size
is 16,384 and the dimension is 50.

###  5.3 Nonlinear Black-Scholes equations

Based on the nonlinear heat equations in Section 5.1, we consider the
corresponding nonlinear Black-Scholes PDE. We take the drift function
$M(t,x)=\mu x$, the diffusion function $\Sigma(t,x)=\sigma\mathrm{diag}(x)$,
and the terminal condition $g(x)=\cos(\overline{x})$ with
$\overline{x}=\sum_{i=1}^{d}x_{i}$.

For the first case, we take

$\displaystyle f(t,x,y,z)=$
$\displaystyle\Big{(}\frac{1}{2}\cos(\overline{x})+\mu\overline{x}\sin(\overline{x})+\sigma^{2}\big{(}\sum_{i=1}^{d}x_{i}^{2}\big{)}\cos(\overline{x})\Big{)}e^{\frac{T-t}{2}}$
(58)
$\displaystyle-\frac{1}{2d}\left(\sigma\overline{x}\cos(\overline{x})\sin(\overline{x})e^{T-t}\right)^{2}+\frac{1}{2d}(y(\bm{1}_{d}\cdot
z))^{2}.$

For (58), the exact solution is $u(t,x)=\cos(\overline{x})e^{\frac{T-t}{2}}$.

For the second case, we take a different $f$ from (58), which is

$\displaystyle f(t,x,y,z)=$
$\displaystyle\left(\frac{1}{2}\Phi(\overline{x})-\mu\overline{x}\varphi(\overline{x})+\sigma^{2}\sum_{i=1}^{d}x_{i}^{2}\varphi(\overline{x})\right)e^{\frac{T-t}{2}}$
(59)
$\displaystyle-\frac{1}{2d}\left(\sigma\overline{x}\varphi(\overline{x})\Phi(\overline{x})e^{T-t}\right)^{2}+\frac{1}{2d}(y(\bm{1}_{d}\cdot
z))^{2},$

where $\Phi$ is the cumulative distribution function of the one-dimensional
standard Gaussian distribution and $\varphi$ is the corresponding probability
density function. We take the terminal function $g=\Phi(\overline{x})$. The
solution of this nonlinear Black-Scholes PDE is
$u(t,x)=\Phi(\overline{x})e^{\frac{T-t}{2}}$.

For the nonlinear term $\eqref{eq:f_BS_1}$, we take the region
$[a,b]^{d}=[-0.5,0.5]^{d}$. Table 3 presents the mean and standard deviation
of the relative $L^{2}$ error over 8 independent runs, for dimensions 20 and
50, respectively. It demonstrates that the RQMC-based algorithm outperforms
the MC-based one. Figure 2 presents the density histogram of pointwise errors
between the trained networks and the exact solution. The batch size is set to
$2^{14}$, and the neural networks with the largest and smallest $L_{2}$ errors
are selected from 8 independent runs. For the MC method, the smallest $L_{2}$
error is $1.74\times 10^{-2}$ and the largest $L_{2}$ error is $4.58\times
10^{-2}$. For the RQMC method, the smallest $L_{2}$ error is $9.01\times
10^{-3}$ and the largest $L_{2}$ error is $1.64\times 10^{-2}$. The figure
demonstrates that the RQMC-based algorithm achieves smaller errors across the
majority of samples. The histogram only displays the errors within the range
of $[-0.1,0.1]$. For the best trained networks using MC methods, 138 out of
16,384 errors fall outside this range, while for the RQMC method, only 2 out
of 16,384 errors fall outside this range.

For the nonlinear term $\eqref{eq:f_BS_2}$, we also take the region
$[a,b]^{d}=[-0.5,0.5]^{d}$. Table 4 presents the mean and standard deviation
of the relative $L^{2}$ error over 8 independent runs, for dimensions 20 and
50, respectively. Both tables demonstrate that the RQMC-based algorithm
outperforms the MC-based one. Figure 4 presents the density histogram of
pointwise errors between the trained networks and the exact solution. For the
MC method, the smallest $L_{2}$ error is $3.05\times 10^{-3}$ and the largest
$L_{2}$ error is $1.22\times 10^{-2}$. For the RQMC method, the smallest
$L_{2}$ error is $1.43\times 10^{-3}$ and the largest $L_{2}$ error is
$3.88\times 10^{-3}$. The figure demonstrates that the RQMC-based algorithm
achieves smaller errors across the majority of samples.

Dimension | Batch Size | Method | Mean Error | Standard Deviation
---|---|---|---|---
20 | $2^{12}$ | MC | $4.6\times 10^{-2}$ | $8.8\times 10^{-3}$
RQMC | $\bm{3.4\times 10^{-2}}$ | $\bm{1.3\times 10^{-3}}$
$2^{14}$ | MC | $2.0\times 10^{-2}$ | $7.9\times 10^{-3}$
RQMC | $\bm{9.3\times 10^{-3}}$ | $\bm{1.3\times 10^{-3}}$
$2^{16}$ | MC | $9.8\times 10^{-3}$ | $2.5\times 10^{-3}$
RQMC | $\bm{5.1\times 10^{-3}}$ | $\bm{4.0\times 10^{-4}}$
50 | $2^{12}$ | MC | $7.2\times 10^{-2}$ | $2.1\times 10^{-2}$
RQMC | $\bm{4.0\times 10^{-2}}$ | $\bm{2.4\times 10^{-3}}$
$2^{14}$ | MC | $3.1\times 10^{-2}$ | $7.7\times 10^{-3}$
RQMC | $\bm{1.1\times 10^{-2}}$ | $\bm{2.3\times 10^{-3}}$
$2^{16}$ | MC | $1.4\times 10^{-2}$ | $4.2\times 10^{-3}$
RQMC | $\bm{6.2\times 10^{-3}}$ | $\bm{2.3\times 10^{-3}}$
Table 3: The mean and standard deviation of the relative $L_{2}$ errors over 8
independent runs for solving the nonlinear Black-Scholes equation with the
nonlinear term (58).

(a) Best trained neural networks

(b) Worst trained neural networks

Figure 3: The density histograms of the pointwise absolute errors for the best and worst trained neural networks in solving the nonlinear Black-Scholes equation with the nonlinear term (58). The batch size is 16,384, and the dimension is 50. Dimension | Batch Size | Method | Mean Error | Standard Deviation
---|---|---|---|---
20 | $2^{12}$ | MC | $1.5\times 10^{-2}$ | $5.7\times 10^{-3}$
RQMC | $\bm{7.5\times 10^{-3}}$ | $\bm{1.2\times 10^{-3}}$
$2^{14}$ | MC | $5.7\times 10^{-3}$ | $4.0\times 10^{-3}$
RQMC | $\bm{2.4\times 10^{-3}}$ | $\bm{2.5\times 10^{-4}}$
$2^{16}$ | MC | $3.6\times 10^{-3}$ | $1.1\times 10^{-3}$
RQMC | $\bm{2.2\times 10^{-3}}$ | $\bm{1.4\times 10^{-4}}$
50 | $2^{12}$ | MC | $1.8\times 10^{-2}$ | $8.2\times 10^{-3}$
RQMC | $\bm{8.1\times 10^{-3}}$ | $\bm{2.4\times 10^{-3}}$
$2^{14}$ | MC | $7.8\times 10^{-3}$ | $3.4\times 10^{-3}$
RQMC | $\bm{2.4\times 10^{-3}}$ | $\bm{8.3\times 10^{-4}}$
$2^{16}$ | MC | $3.6\times 10^{-3}$ | $1.2\times 10^{-3}$
RQMC | $\bm{1.9\times 10^{-3}}$ | $\bm{3.3\times 10^{-4}}$
Table 4: The mean and standard deviation of the relative $L^{2}$ errors over 8
independent runs for solving the nonlinear Black-Scholes equation with the
nonlinear term (59).

(a) Best trained neural networks

(b) Worst trained neural networks

Figure 4: The density histograms of the pointwise absolute errors for the best
and worst trained neural networks in solving the nonlinear Black-Scholes
equation with the nonlinear term (59). The batch size is 16,384, and the
dimension is 50.

##  6 Conclusion

Nonlinear PDEs (1) represent a general form of the linear Kolmogorov PDEs [2,
3, 4, 5, 6], related to BSDEs by nonlinear Feymann-Kac formula [9], which find
extensive applications in financial mathematics [7]. Therefore, efficiently
solving nonlinear PDEs constitutes a significant challenge. We investigated
the convergence orders of the generalization error within the DBDP [1] method
under both the MC and QMC approaches. Our theoretical findings indicate that
the generalization error under the QMC method exhibits a higher convergence
order compared to the MC method. Theoretically, this suggests that the use of
the QMC method for sampling during training yields superior outcomes. We
validated these theoretical insights through several numerical experiments.
The results obtained using the QMC method have smaller errors and improved
stability. It is an intriguing prospect to explore the trade-off among scheme
error, approximation error, and generalization error to achieve an optimal
total error with reduced computational costs. This is left for future work.

##  Appendix A Proof of Theorem 1

###### Proof.

In Theorem 4.1 of Hure et al. [1], they propose a mid-term for the proof

$\left\\{\begin{aligned}
\mathcal{V}_{i}(X_{t_{i}})&=\mathbb{E}_{i}\left[\mathcal{U}_{i+1}(X_{t_{i+1}})\right]+f(t_{i},X_{t_{i}},\mathcal{V}_{i}(X_{t_{i}}),\overline{\widehat{Z}}_{i}(X_{t_{i}}))\Delta
t_{i},\\\ \overline{\widehat{Z}}_{i}(X_{t_{i}})&=\frac{1}{\Delta
t_{i}}\mathbb{E}_{i}\left[\mathcal{U}_{i+1}(X_{t_{i+1}})\Delta
W_{t_{i}}\right],\end{aligned}\right.$ (60)

where
$\mathbb{E}_{i}\left[\cdot\right]:=\mathbb{E}\left[\cdot|X_{t_{i}}\right]$. By
the martingale representation theorem [62], there is a $\widehat{Z}_{s}$ such
that

$\mathcal{U}_{i+1}(X_{t_{i+1}})=\mathcal{V}_{i}(X_{t_{i}})-f(t_{i},X_{t_{i}},\mathcal{V}_{i}(X_{t_{i}}),\overline{\widehat{Z}}_{i}(X_{t_{i}}))\Delta
t_{i}+\int_{t_{i}}^{t_{i+1}}\widehat{Z}_{s}\mathrm{d}W_{s}.$ (61)

and by Itô isometry, we have

$\overline{\widehat{Z}}_{i}(X_{t_{i}})=\frac{1}{\Delta
t_{i}}\mathbb{E}_{i}\left[\int_{t_{i}}^{t_{i+1}}\widehat{Z}_{s}\mathrm{d}s\right].$
(62)

Hure et al. [1] obtain the following result

$\displaystyle\varepsilon\left[\left(\mathcal{U},\mathcal{Z}\right),\left(Y,Z\right)\right]\leq
C\bigg{(}$
$\displaystyle\mathbb{E}\left|g(\mathcal{X}_{T})-g(X_{N})\right|^{2}+|\pi|+\varepsilon^{Z}(\pi)$
(63)
$\displaystyle+\sum_{i=0}^{N-1}\Big{(}N\mathbb{E}\left|\mathcal{V}_{i}(X_{t_{i}})-\mathcal{U}_{i}(X_{t_{i}})\right|^{2}+\mathbb{E}\left|\overline{\widehat{Z}}_{i}(X_{t_{i}})-\mathcal{Z}_{i}(X_{t_{i}})\right|^{2}\Big{)}\bigg{)},$

where

$\varepsilon^{Z}(\pi)=\mathbb{E}\left[\sum_{i=0}^{N-1}\int_{t_{i}}^{t_{i+1}}\left|Z_{t}-\overline{Z}_{t_{i}}\right|^{2}\mathrm{d}t\right],\quad\overline{Z}_{t_{i}}=\frac{1}{\Delta
t_{i}}\mathbb{E}_{i}\left[\int_{t_{i}}^{t_{i+1}}Z_{t}\mathrm{d}t\right].$

Note that by substituting (61) into
$L_{i}(\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i})$, we can decompose
it as follows

$\displaystyle
L_{i}(\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i})=\tilde{L}_{i}(\mathcal{U}_{i},\mathcal{Z}_{i})+\mathbb{E}\left[\int_{t_{i}}^{t_{i+1}}\left|\widehat{Z}_{t}-\overline{\widehat{Z}}_{i}(X_{t_{i}})\right|^{2}\mathrm{d}t\right],$
(64)

where

$\displaystyle\tilde{L}_{i}(\mathcal{U}_{i},\mathcal{Z}_{i})=\Delta
t_{i}\mathbb{E}\left|\overline{\widehat{Z}}_{i}(X_{t_{i}})-\mathcal{Z}_{i}(X_{t_{i}})\right|^{2}$
(65) $\displaystyle+$
$\displaystyle\mathbb{E}\left|\mathcal{V}_{i}(X_{t_{i}})-\mathcal{U}_{i}(X_{t_{i}})+\left(f(t_{i},X_{t_{i}},\mathcal{U}_{i}(X_{t_{i}}),\mathcal{Z}_{i}(X_{t_{i}}))-f(t_{i},X_{t_{i}},\mathcal{V}_{i}(X_{t_{i}}),\overline{\widehat{Z}}_{i}(X_{t_{i}}))\right)\Delta
t_{i}\right|^{2}.$

Note that the second term of the right hand side of (64) is independent of
$(\mathcal{U}_{i},\mathcal{Z}_{i})$, thus it is the minimize of
$L_{i}(\mathcal{U},\mathcal{Z})$. By (60), $\mathcal{V}$ and
$\overline{\widehat{Z}}$ are functions of $X_{t_{i}}$. Therefore, the
minimizer $(\mathcal{U}^{*},\mathcal{Z}^{*})$ of (64) satisfies

$\mathcal{U}_{i}^{*}(X_{t_{i}})=\mathcal{V}_{t_{i}}(X_{t_{i}}),\quad\mathcal{Z}_{i}^{*}(X_{t_{i}})=\overline{\widehat{Z}}_{i}(X_{t_{i}}),$
(66)

Moreover,

$\tilde{L}_{i}(\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*})=0,\quad
L_{i}(\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*})=\mathbb{E}\left[\int_{t_{i}}^{t_{i+1}}\left|\widehat{Z}_{t}-\overline{\widehat{Z}}_{i}(X_{t_{i}})\right|^{2}\mathrm{d}t\right].$
(67)

By the lipschitz condition of $f$, there is a constant $C$, such that

$\left\\{\begin{aligned}
\tilde{L}_{i}\left(\mathcal{U}_{i},\mathcal{Z}_{i}\right)&\leq\left(1+C\Delta
t_{i}\right)\mathbb{E}\left|\mathcal{V}_{i}(X_{t_{i}})-\mathcal{U}_{i}(X_{t_{i}})\right|^{2}+C\Delta
t_{i}\mathbb{E}\left|\overline{\widehat{Z}}_{i}(X_{t_{i}})-\mathcal{Z}_{i}(X_{t_{i}})\right|^{2},\\\
\tilde{L}_{i}\left(\mathcal{U}_{i},\mathcal{Z}_{i}\right)&\geq\left(1-C\Delta
t_{i}\right)\mathbb{E}\left|\mathcal{V}_{i}(X_{t_{i}})-\mathcal{U}_{i}(X_{t_{i}})\right|^{2}+C\frac{\Delta
t_{i}}{2}\mathbb{E}\left|\overline{\widehat{Z}}_{i}(X_{t_{i}})-\mathcal{Z}_{i}(X_{t_{i}})\right|^{2}.\end{aligned}\right.$
(68)

It follows from (63), (64) and (68) that

$\displaystyle\varepsilon\left[\left(\mathcal{U},\mathcal{Z}\right),\left(Y,Z\right)\right]\leq
C\bigg{(}$
$\displaystyle\mathbb{E}\left|g(\mathcal{X}_{T})-g(X_{N})\right|^{2}+|\pi|+\varepsilon^{Z}(\pi)$
(69)
$\displaystyle+N\sum_{i=0}^{N-1}\left(L_{i}\left(\mathcal{U}_{i+1},\mathcal{U}_{i},\mathcal{Z}_{i}\right)-L_{i}\left(\mathcal{U}_{i+1},\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*}\right)\right)\bigg{)},$

where
$\left(\mathcal{U}_{i}^{*},\mathcal{Z}_{i}^{*}\right)=\arg\min_{(\mathcal{U},\mathcal{Z})}L_{i}\left(\mathcal{U}_{i+1},\mathcal{U},\mathcal{Z}\right)$,
which has the representation (66). Noting that $g$ has the lipschitz
condition, by Zhang [12], we have

$\varepsilon^{Z}(\pi)=O(|\pi|).$

Thus under the conditions of Theorem 1, (15) holds. ∎

## References

  * Huré et al. [2020] authorC. Huré, authorH. Pham, authorX. Warin,  titleDeep backward schemes for high-dimensional nonlinear PDEs,  journalMathematics of Computation volume89 (year2020) pages1547–1579.
  * Beck et al. [2021] authorC. Beck, authorS. Becker, authorP. Grohs, authorN. Jaafari, authorA. Jentzen,  titleSolving the kolmogorov pde by means of deep learning,  journalJournal of Scientific Computing volume88 (year2021) pages1–28.
  * Berner et al. [2020] authorJ. Berner, authorM. Dablander, authorP. Grohs,  titleNumerically solving parametric families of high-dimensional kolmogorov partial differential equations via deep learning,  journalAdvances in neural information processing systems volume33 (year2020) pages16615–16627.
  * De Ryck and Mishra [2022] authorT. De Ryck, authorS. Mishra,  titleError analysis for physics-informed neural networks (pinns) approximating kolmogorov pdes,  journalAdvances in Computational Mathematics volume48 (year2022) pages79.
  * Glau and Wunderlich [2022] authorK. Glau, authorL. Wunderlich,  titleThe deep parametric PDE method and applications to option pricing,  journalApplied Mathematics and Computation volume432 (year2022) pages127355.
  * Richter and Berner [2022] authorL. Richter, authorJ. Berner,  titleRobust sde-based variational formulations for solving linear pdes via deep learning,  in: booktitleInternational Conference on Machine Learning, organizationPMLR, year2022, pp. pages18649–18666.
  * El Karoui et al. [1997] authorN. El Karoui, authorS. Peng, authorM. C. Quenez,  titleBackward stochastic differential equations in finance,  journalMathematical Finance volume7 (year1997) pages1–71.
  * Pardoux and Peng [1992] authorE. Pardoux, authorS. Peng,  titleBackward stochastic differential-equations and quasi-linear parabolic partial-differential equations,  journalLecture Notes in Control and Information Sciences volume176 (year1992) pages200–217.
  * Pardoux and Peng [1990] authorE. Pardoux, authorS. G. Peng,  titleAdapted solution of a backward stochastic differential-equation,  journalSystems & Control Letters volume14 (year1990) pages55–61.
  * Bouchard and Touzi [2004] authorB. Bouchard, authorN. Touzi,  titleDiscrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations,  journalStochastic Processes and their Applications volume111 (year2004) pages175–206.
  * Gobet et al. [2005] authorE. Gobet, authorJ.-P. Lemor, authorX. Warin,  titleA regression-based monte carlo method to solve backward stochastic differential equations,  journalThe Annals of Applied Probability volume15 (year2005).
  * Zhang [2004] authorJ. F. Zhang,  titleA numerical scheme for BSDEs,  journalAnnals of Applied Probability volume14 (year2004) pages459–488.
  * E et al. [2017] authorW. E, authorJ. Han, authorA. Jentzen,  titleDeep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations,  journalCommunications in Mathematics and Statistics volume5 (year2017) pages349–380.
  * Han et al. [2018] authorJ. Han, authorA. Jentzen, authorW. E,  titleSolving high-dimensional partial differential equations using deep learning,  journalProceedings of the National Academy of Sciences volume115 (year2018) pages8505–8510.
  * Germain et al. [2022] authorM. Germain, authorH. Pham, authorX. Warin,  titleApproximation error analysis of some deep backward schemes for nonlinear PDEs,  journalSIAM Journal on Scientific Computing volume44 (year2022) pagesA28–A56.
  * Beck et al. [2021] authorC. Beck, authorS. Becker, authorP. Cheridito, authorA. Jentzen, authorA. Neufeld,  titleDeep splitting method for parabolic PDEs,  journalSIAM Journal on Scientific Computing volume43 (year2021) pagesA3135–A3154.
  * Jiao et al. [2024] authorY. Jiao, authorY. Lai, authorY. Lo, authorY. Wang, authorY. Yang,  titleError analysis of deep ritz methods for elliptic equations,  journalAnalysis and Applications volume22 (year2024) pages57–87.
  * Jiao et al. [2023] authorY. Jiao, authorG. Shen, authorY. Lin, authorJ. Huang,  titleDeep nonparametric regression on approximate manifolds: Nonasymptotic error bounds with polynomial prefactors,  journalThe Annals of Statistics volume51 (year2023) pages691–716.
  * Beck et al. [2022] authorC. Beck, authorA. Jentzen, authorB. Kuckuck,  titleFull error analysis for the training of deep neural networks,  journalInfinite Dimensional Analysis, Quantum Probability and Related Topics volume25 (year2022) pages2150020.
  * Berner et al. [2020] authorJ. Berner, authorP. Grohs, authorA. Jentzen,  titleAnalysis of the generalization error: Empirical risk minimization over deep artificial neural networks overcomes the curse of dimensionality in the numerical approximation of Black–Scholes partial differential equations,  journalSIAM Journal on Mathematics of Data Science volume2 (year2020) pages631–657.
  * Jentzen and Welti [2023] authorA. Jentzen, authorT. Welti,  titleOverall error analysis for the training of deep neural networks via stochastic gradient descent with random initialisation,  journalApplied Mathematics and Computation volume455 (year2023) pages127907.
  * Mishra and Molinaro [2023] authorS. Mishra, authorR. Molinaro,  titleEstimates on the generalization error of physics-informed neural networks for approximating pdes,  journalIMA Journal of Numerical Analysis volume43 (year2023) pages1–43.
  * Glasserman [2004] authorP. Glasserman, titleMonte Carlo Methods in Financial Engineering, publisherSpringer, addressNew York, year2004.
  * Hull [2003] authorJ. C. Hull, titleOptions Futures and Other Derivatives, publisherPearson Education India, addressNew York, year2003.
  * Dick and Feischl [2021] authorJ. Dick, authorM. Feischl,  titleA quasi-Monte Carlo data compression algorithm for machine learning,  journalJournal of Complexity volume67 (year2021) pages101587.
  * Liu and Owen [2021] authorS. Liu, authorA. B. Owen,  titleQuasi-Monte Carlo quasi-newton in variational bayes,  journalJournal of Machine Learning Research volume22 (year2021) pages1–23.
  * Longo et al. [2021] authorM. Longo, authorS. Mishra, authorT. K. Rusch, authorC. Schwab,  titleHigher-order quasi-monte carlo training of deep neural networks,  journalSIAM Journal on Scientific Computing volume43 (year2021) pagesA3938–A3966.
  * Mishra and Rusch [2021] authorS. Mishra, authorT. K. Rusch,  titleEnhancing accuracy of deep learning algorithms by training with low-discrepancy sequences,  journalSIAM Journal on Numerical Analysis volume59 (year2021) pages1811–1834.
  * Xiao et al. [2024] authorJ. Xiao, authorF. Fu, authorX. Wang,  titleDeep learning based on randomized quasi-monte carlo method for solving linear kolmogorov partial differential equation,  journalJournal of Computational and Applied Mathematics (year2024) pages116088.
  * Kloeden and Platen [1992] authorP. E. Kloeden, authorE. Platen, titleNumerical Solution of Stochastic Differential Equations, publisherSpringer, addressBerlin, year1992.
  * Milstein and Tretyakov [2006] authorG. N. Milstein, authorM. V. Tretyakov,  titleNumerical algorithms for forward-backward stochastic differential equations,  journalSIAM Journal on Scientific Computing volume28 (year2006) pages561–582.
  * Hornik [1991] authorK. Hornik,  titleApproximation capabilities of multilayer feedforward networks,  journalNeural networks volume4 (year1991) pages251–257.
  * Hornik [1993] authorK. Hornik,  titleSome new results on neural network approximation,  journalNeural networks volume6 (year1993) pages1069–1072.
  * Gühring et al. [2020] authorI. Gühring, authorG. Kutyniok, authorP. Petersen,  titleError bounds for approximations with deep relu neural networks in $W^{s,p}$ norms,  journalAnalysis and Applications volume18 (year2020) pages803–859.
  * Siegel and Xu [2020] authorJ. W. Siegel, authorJ. Xu,  titleApproximation rates for neural networks with general activation functions,  journalNeural Networks volume128 (year2020) pages313–321.
  * Elbrächter et al. [2021] authorD. Elbrächter, authorD. Perekrestenko, authorP. Grohs, authorH. Bölcskei,  titleDeep neural network approximation theory,  journalIEEE Transactions on Information Theory volume67 (year2021) pages2581–2623.
  * Grohs et al. [2023] authorP. Grohs, authorF. Hornung, authorA. Jentzen, authorP. Von Wurstemberger, titleA proof that artificial neural networks overcome the curse of dimensionality in the numerical approximation of Black–Scholes partial differential equations, publisherAmerican Mathematical Society, year2023.
  * Jentzen et al. [2021] authorA. Jentzen, authorD. Salimova, authorT. Welti,  titleA proof that deep artificial neural networks overcome the curse of dimensionality in the numerical approximation of kolmogorov partial differential equations with constant diffusion and nonlinear drift coefficients,  journalCommunications in Mathematical Sciences volume19 (year2021) pages1167–1205.
  * Hutzenthaler et al. [2020a] authorM. Hutzenthaler, authorA. Jentzen, authorT. Kruse, authorT. A. Nguyen,  titleA proof that rectified deep neural networks overcome the curse of dimensionality in the numerical approximation of semilinear heat equations,  journalSN partial differential equations and applications volume1 (year2020a) pages10.
  * Hutzenthaler et al. [2020b] authorM. Hutzenthaler, authorA. Jentzen, authorT. Kruse, authorT. Anh Nguyen, authorP. von Wurstemberger,  titleOvercoming the curse of dimensionality in the numerical approximation of semilinear parabolic partial differential equations,  journalProceedings of the Royal Society A volume476 (year2020b) pages20190630.
  * Darbon and Osher [2016] authorJ. Darbon, authorS. Osher,  titleAlgorithms for overcoming the curse of dimensionality for certain Hamilton–Jacobi equations arising in control theory and elsewhere,  journalResearch in the Mathematical Sciences volume3 (year2016) pages19.
  * Hutzenthaler et al. [2023] authorM. Hutzenthaler, authorA. Jentzen, authorT. Kruse, authorT. Anh Nguyen,  titleOvercoming the curse of dimensionality in the numerical approximation of backward stochastic differential equations,  journalJournal of Numerical Mathematics volume31 (year2023) pages1–28.
  * Robbins and Monro [1951] authorH. Robbins, authorS. Monro,  titleA stochastic approximation method,  journalThe annals of mathematical statistics (year1951) pages400–407.
  * Kingma and Ba [2014] authorD. P. Kingma, authorJ. Ba,  titleAdam: A method for stochastic optimization,  journalarXiv preprint arXiv:1412.6980 (year2014).
  * Ruder [2016] authorS. Ruder,  titleAn overview of gradient descent optimization algorithms,  journalarXiv preprint arXiv:1609.04747 (year2016).
  * Vapnik [1991] authorV. Vapnik,  titlePrinciples of risk minimization for learning theory,  journalAdvances in neural information processing systems volume4 (year1991).
  * Shalev-Shwartz and Ben-David [2014] authorS. Shalev-Shwartz, authorS. Ben-David, titleUnderstanding machine learning: From theory to algorithms, publisherCambridge university press, year2014.
  * Mohri et al. [2018] authorM. Mohri, authorA. Rostamizadeh, authorA. Talwalkar, titleFoundations of machine learning, publisherMIT press, year2018.
  * Zhang [2004] authorT. Zhang,  titleStatistical behavior and consistency of classification methods based on convex risk minimization,  journalThe Annals of Statistics volume32 (year2004) pages56–85.
  * Niederreiter [1992] authorH. Niederreiter, titleRandom Number Generation and Quasi-Monte Carlo Methods, publisherSIAM, addressPhiladelphia, PA, year1992.
  * Caflisch et al. [1997] authorR. E. Caflisch, authorW. J. Morokoff, authorA. B. Owen,  titleValuation of mortgage-backed securities using Brownian bridges to reduce effective dimension,  journalJournal of Computational Finance volume1 (year1997) pages27–46.
  * Owen [2013] authorA. B. Owen, titleMonte Carlo Theory, Methods and Examples, publisherStanford, year2013.
  * Hlawka and Mück [1972] authorE. Hlawka, authorR. Mück,  titleÜber eine transformation von gleichverteilten folgen II,  journalComputing volume9 (year1972) pages127–138.
  * Kuo et al. [2006] authorF. Y. Kuo, authorG. W. Wasilkowski, authorB. J. Waterhouse,  titleRandomly shifted lattice rules for unbounded integrands,  journalJournal of Complexity volume22 (year2006) pages630–651.
  * Owen [1997] authorA. B. Owen,  titleScrambled net variance for integrals of smooth functions,  journalThe Annals of Statistics volume25 (year1997) pages1541–1562.
  * Owen [2008] authorA. B. Owen,  titleLocal antithetic sampling with scrambled nets,  journalThe Annals of Statistics volume36 (year2008) pages2319–2343.
  * Ouyang et al. [2023] authorD. Ouyang, authorX. Wang, authorZ. He,  titleQuasi-Monte Carlo for unbounded integrands with importance sampling,  journalarXiv preprint arXiv:2310.00650 (year2023).
  * Owen [2006] authorA. B. Owen,  titleHalton sequences avoid the origin,  journalSIAM Review volume48 (year2006) pages487–503.
  * Glorot and Bengio [2010] authorX. Glorot, authorY. Bengio,  titleUnderstanding the difficulty of training deep feedforward neural networks,  in: booktitleProceedings of the thirteenth international conference on artificial intelligence and statistics, organizationJMLR Workshop and Conference Proceedings, year2010, pp. pages249–256.
  * Ioffe and Szegedy [2015] authorS. Ioffe, authorC. Szegedy,  titleBatch normalization: Accelerating deep network training by reducing internal covariate shift,  in: booktitleInternational conference on machine learning, organizationpmlr, year2015, pp. pages448–456.
  * Cole [1951] authorJ. D. Cole,  titleOn a quasi-linear parabolic equation occurring in aerodynamics,  journalQuarterly of applied mathematics volume9 (year1951) pages225–236.
  * Karatzas and Shreve [2014] authorI. Karatzas, authorS. Shreve, titleBrownian motion and stochastic calculus, volume volume113, publisherspringer, year2014.