# Convergence Rate Analysis for Deep Ritz Method 111This work has been
presented at workshop on ”Data Driven Scientific Computing” held in Southern
University of Science and Technology, China, on November 21-22, 2020.

Chenguang Duan School of Mathematics and Statistics, Wuhan University, Wuhan
430072, P.R. China<EMAIL_ADDRESS>Yuling Jiao School of Mathematics
and Statistics, and Hubei Key Laboratory of Computational Science, Wuhan
University, Wuhan 430072, P.R. China<EMAIL_ADDRESS>Yanming Lai
School of Mathematics and Statistics, Wuhan University, Wuhan 430072, P.R.
China<EMAIL_ADDRESS>Xiliang Lu School of Mathematics and
Statistics, and Hubei Key Laboratory of Computational Science, Wuhan
University, Wuhan 430072, P.R. China<EMAIL_ADDRESS>Jerry Zhijian
Yang School of Mathematics and Statistics, and Hubei Key Laboratory of
Computational Science, Wuhan University, Wuhan 430072, P.R. China.
<EMAIL_ADDRESS>
###### Abstract

Using deep neural networks to solve PDEs has attracted a lot of attentions
recently. However, why the deep learning method works is falling far behind
its empirical success. In this paper, we provide a rigorous numerical analysis
on deep Ritz method (DRM) [48] for second order elliptic equations with
Neumann boundary conditions. We establish the first nonasymptotic convergence
rate in $H^{1}$ norm for DRM using deep networks with $\mathrm{ReLU}^{2}$
activation functions. In addition to providing a theoretical justification of
DRM, our study also shed light on how to set the hyper-parameter of depth and
width to achieve the desired convergence rate in terms of number of training
samples. Technically, we derive bound on the approximation error of deep
$\mathrm{ReLU}^{2}$ network in $H^{1}$ norm and bound on the Rademacher
complexity of the non-Lipschitz composition of gradient norm and
$\mathrm{ReLU}^{2}$ network, both of which are of independent interest.

##  1 Introduction

Partial differential equations (PDEs) have broad applications in physics,
chemistry, biology, geology and engineering. A great deal of efforts have been
devoted to studying numerical methods for solving PDEs [5, 7, 22, 44, 16].
However, it is still a challenging task to develop numerical scheme for
solving PDEs in high-dimension. Due to the success of deep learning for high-
dimensional data analysis in computer vision and natural language processing,
people have been paying more attention to using (deep) neural network to solve
PDEs in high dimension with may be complex domain, an idea that goes back to
1990’s [21, 19]. In the last few years, there are growing literatures on
neural network based numerical methods for PDEs These works can be roughly
classified into two categories.

In the first category, deep neural networks are used to improve classical
methods. [10] designs a neural network to estimate artificial viscosity in
discontinuous Galerkin schemes, see also [6]. [32] trains a neural network
serving as a troubled-cell indicator in high-resolution schemes for
conservation laws. [42] proposes a universal discontinuity detector using
convolution neural network and applies it in conjunction of solving nonlinear
conservation. [47] uses reinforcement learning to find new and potentially
better data-driven solvers for conservation laws.

In the second category, deep neural networks are utilized to approximate the
solution of the PDEs directly. Being benefit from the excellent approximation
power of deep neural networks and SGD training, these methods have been
successfully applied to solve PDEs in high-dimension. [3, 9] convert nonlinear
parabolic PDEs into backward stochastic differential equations and solve them
by deep neural networks, which can deal with high-dimensional problems.
Methods based on the strong form of PDEs [31, 40] are also proposed. In [31],
physics-informed neural networks (PINNs) use the squared residuals on the
domain as the loss function and treat boundary conditions as penalty term.
There are several extensions of PINNs for different types of PDEs, including
fractional PINNs [30], nonlocal PINNs [29], conservative PINNs [18], eXtended
PINNs [17], among others. A similar method presented in [25] proposes a
residual-based adaptive refinement method to improve the training efficiency.

In contrast to minimizing squared residuals of strong form, a natural
alternative approach to derive loss functions are based on the variational
form of PDEs [48, 51]. Inspired by Ritz method, [48] proposes deep Ritz method
(DRM) to solve variational problems arising from PDEs. The idea of Galerkin
method has also been used in [51], where, they propose a deep Galerkin method
(DGM) via reformulating the problem of finding the weak solution of PDEs into
an operator norm minimization problem induced by the weak formulation.

###  1.1 Related works and contributions

Although there are great empirical achievements in recent years as mentioned
above, a challenging and interesting question is that can we give rigorous
analysis to guarantee their performances as people has done in the classical
counterpart such as finite element method (FEM) [7] and finite difference
method [22] ? Several recent efforts have been devoted to making processes
along this line. [26] consider the optimization and generalization error of
second-order linear PDEs with two-layer neural networks in the scenario of
over-parametrization. [37, 27, 38] study the convergence of PINNs with deep
neural networks. When we were about to finish our draft, we aware that [24]
give an error analysis that focuses on analyzing one hidden layer shallow
networks with ReLU-Cosine activation functions to solve elliptic PDEs whose
solutions are restricted to spectral Barron space, see also [49] for handling
general equations with solutions living in spectral Barron space via two layer
$\mathrm{ReLU}^{k}$ networks.  Two important questions have not been addressed
in the above mentioned related study are those: what is the influence of the
topological structure of the networks, say the depth and width, in the
quantitative error analysis ? How to determine these hyper-parameters to
achieve a desired convergence rate ? In this paper, we give a firm answers on
these questions by studying convergence ratec of the deep Ritz method to solve
second order elliptic equations with Neumann boundary conditions by using
$\mathrm{ReLU}^{2}$ networks with arbitrary depth. As far as we know, we
establish the first nonasymptotic bound on DRM. The main contritions of this
paper are summarized as follows.

  * •

We derive a bound on the approximation error of deep $\mathrm{ReLU}^{2}$
network in $H^{1}$ norm, which is of independent interest, i.e., we prove that
for any $u^{*}\in H^{2}(\Omega)$, there exist a $\mathrm{ReLU}^{2}$ network
$\bar{u}_{\bar{\phi}}$ with depth $\mathcal{D}\leq\lceil\log_{2}d\rceil+3,$
width
$\mathcal{W}\leq\mathcal{O}(4d\left\lceil\frac{1}{\epsilon}-4\right\rceil^{d})$
such that

$\|u^{*}-\bar{u}_{\bar{\phi}}\|_{H^{1}\left(\Omega\right)}\leq\epsilon.$

  * •

We establish a bound on the statistical error in DRM with the tools of Pseudo
dimension, especially we give an bound on

$\mathbb{E}_{Z_{i},\sigma_{i},i=1,...,n}[\sup_{u_{\phi}\in\mathcal{N}^{2}}\frac{1}{n}|\sum_{i}\sigma_{i}\|\nabla
u_{\phi}(Z_{i})\|^{2}|],$

i.e., the Rademacher complexity of the non-Lipschitz composition of gradient
norm and $\mathrm{ReLU}^{2}$ network, via calculating the Pseudo dimension of
networks with both $\mathrm{ReLU}$ and $\mathrm{ReLU}^{2}$ activation
functions. We believe that the technique we used here is helpful for bounding
the statistical errors for other deep PDEs solvers where the Rademacher
complexity of non-Lipschitz composition is hard to handle.

  * •

Based on the above to error bounds we establish the first nonasymptotic
convergence rate of deep Ritz method. We prove that if we set the depth and
width in $\mathrm{ReLU}^{2}$ networks to be

$\mathcal{D}\leq\lceil\log_{2}d\rceil+3,\mathcal{W}\leq\mathcal{O}(4d\left\lceil
n^{\frac{1}{d+2+\nu}}-4\right\rceil^{d}),$

the $H^{1}$ norm error of DRM in expectation is

${\mathcal{O}}(n^{-1/(2d+4+\nu)}),$

where $n$ is the number of training samples on both the domain and the
boundary, $\nu$ is a positive number that can be an arbitrary small. Our
theory shed lights on choosing the topological structure of the employed
networks to achieve the desired convergence rate in terms of number of
training samples.

  * •

By comparing the known results in nonparametric regression, where the optimal
convergence rate in $H^{1}$ norm for estimating functions in $H^{2}$ with $n$
paired samples is $\mathcal{O}(n^{-\frac{1}{4+d}})$ [41], we conjecture that
the optimal convergence rate of DRM in $H^{1}$ norm is also
${\mathcal{O}}(n^{-1/(d+4)}).$

The rest of the paper are organized as follows. In Section 2, we give some
preliminaries. In Section 3, we present the detail analysis on the convergence
rate of DRM. We give conclusion and short discussion in Section 4.

##  2 Preliminaries

Consider the following elliptic equation with Neumann boundary conditions

$\left\\{\begin{aligned} -\triangle u+wu&=f\text{ in }\Omega\\\ \frac{\partial
u}{\partial n}&=g\text{ on }\partial\Omega,\\\ \end{aligned}\right.$ (1)

where, $\Omega$ is a bounded open subset of $\mathbb{R}^{d},d>1$, $f(x)\in
L^{2}(\Omega)$, $w(x)\in L^{\infty}(\Omega)$ satisfying $w(x)\geq c_{1}>0$
a.e., and $g(s)\in L^{2}(\partial\Omega)$. Without loss of generality we
assume $\Omega=(0,1)^{d}.$ Define

$\mathcal{L}(u)=\frac{1}{2}|u|_{H^{1}(\Omega)}^{2}+\frac{1}{2}\|u\|_{L^{2}(\Omega;{w})}^{2}-\langle
u,f\rangle_{L^{2}({\Omega})}-\langle{Tu},g\rangle_{L^{2}({\partial\Omega})},$
(2)

where $T$ is the trace operator.

######  Lemma 1.

The unique weak solution $u^{*}\in H^{1}(\Omega)$ of (1) is the unique
minimizer of $\mathcal{L}(u)$ over $H^{1}(\Omega)$. Moreover, $u^{*}\in
H^{2}(\Omega)$.

###### Proof.

Well known results, see for example [13]. ∎

A function $\mathbf{f}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{N_{L}}$
implemented by a neural network is defined by

$\begin{array}[]{l}\mathbf{f}_{0}(\mathbf{x})=\mathbf{x},\\\
\mathbf{f}_{\ell}(\mathbf{x})=\varrho_{\ell}\left(A_{\ell}\mathbf{f}_{\ell-1}+\mathbf{b}_{\ell}\right)\quad\text{
for }\ell=1,\ldots,L-1,\\\
\mathbf{f}=\mathbf{f}_{L}(\mathbf{x}):=A_{L}\mathbf{f}_{L-1}+\mathbf{b}_{L},\end{array}$

where $A_{\ell}\in\mathbb{R}^{N_{\ell}\times N_{\ell-1}}$,
$\mathbf{b}_{\ell}\in\mathbb{R}^{N_{\ell}}$ and the activation function
$\varrho_{\ell}$ is understood to act component-wise (note that here we allow
different activation functions in different layers). $L$ is called the depth
of the network and $\max\\{N_{\ell},\ell=0,\cdots,L\\}$ is called the width of
the network. We will use $\mathcal{L}$ and $\mathcal{W}$ to denote the depth
and width of neural networks $\mathbf{f}$, respectively.
$\sum_{\ell=1}^{L}N_{\ell}$ is called number of unites of $\mathbf{f}$ and
$\phi=\\{A_{\ell},\mathbf{b}_{\ell}\\}_{\ell}$ are called the weight
parameters. For simplicity we also use $\mathbf{f}_{\phi}$ to refer to the
network. We use $\mathcal{N}^{2}_{{\mathcal{D},\mathcal{W},\mathcal{B}}}$ to
denote the set of neural networks with depth $\mathcal{D}$, width
$\mathcal{W}$, output of the function values and their square norm of
gradients bounded by $\mathcal{B}$, activation function
$\mathrm{ReLU}^{2}(x)=\max\\{0,x^{2}\\}$. Denote
$\mathcal{N}^{1,2}_{{\mathcal{D},\mathcal{W},\mathcal{B}}}$ as the set of
neural networks with depth $\mathcal{D}$, width $\mathcal{W}$, output bounded
by $\mathcal{B}$, activation functions $\mathrm{ReLU}(x)=\max\\{0,x\\}$ and
$\mathrm{ReLU}^{2}(x)=\max\\{0,x^{2}\\}$.

Obviously,

$\displaystyle\mathcal{L}(u)$
$\displaystyle=|\Omega|\mathbb{E}_{{X\sim\mathrm{U}(\Omega)}}[\|\nabla
u(X)\|_{2}^{2}/2+w(X)u^{2}(X)/2-u(X)f(X)]$
$\displaystyle-|\partial{\Omega}|\mathbb{E}_{{Y\sim\mathrm{U}(\partial\Omega)}}[{Tu}(Y)g(Y)],$

where $\mathrm{U}(\Omega),\mathrm{U}(\partial\Omega)$ are the uniform
distributions on $\Omega$ and $\partial\Omega$, respectively. The main idea of
deep Ritz method (DRM) [48] is employing a
$u_{\phi}\in\mathcal{N}^{2}:=\mathcal{N}_{\mathcal{D},\mathcal{W},\mathcal{B}}^{2}$
to approximate the minimizer $u^{*}$ of $\mathcal{L}$, i.e., finding
$u_{\phi}$ such that $\mathcal{L}(u_{\phi})$ closes to $\mathcal{L}(u^{*})$.
To this end, by Lemma 1, one may consider the following empirical loss
minimization problem

$\widehat{u}_{\phi}\in\min_{u_{\phi}\in\mathcal{N}^{2}}\widehat{\mathcal{L}}(u_{\phi}),$
(3)

where,

$\displaystyle\widehat{\mathcal{L}}(u_{\phi})$
$\displaystyle=\frac{|\Omega|}{N}\sum_{i=1}^{N}[\frac{\|\nabla
u_{\phi}(X_{i})\|_{2}^{2}}{2}+\frac{w(X_{i})u_{\phi}^{2}(X_{i})}{2}-u_{\phi}(X_{i})f(X_{i})]$
$\displaystyle-\frac{\partial\Omega}{M}\sum_{j=1}^{M}[u_{\phi}(Y_{j})g(Y_{j})],$
(4)

is a discrete version of the functional $\mathcal{L}(u_{\phi})$ with
$\\{X_{i}\\}_{i=1}^{N}$ being identically and independently distributed
(i.i.d.) according to $\mathrm{U}(\Omega)$, $\\{Y_{j}\\}_{j=1}^{M}$ being
identically and independently drawn from $\mathrm{U}(\partial\Omega)$. Then,
we call a (random) solver $\mathcal{A}$, say SGD, to minimize (3) and denote
the output of $\mathcal{A}$, say $u_{\phi_{\mathcal{A}}}\in\mathcal{N}$, as
the final solution.

##  3 Error Analysis

In this section we prove the convergence rate analysis for DRM with deep
$\mathrm{ReLu}^{2}$ networks. The following Lemma play an important role by
decoupling the total errors into three types of errors.

######  Lemma 2.

$\displaystyle\|u_{\phi_{\mathcal{A}}}-u^{*}\|_{H^{1}(\Omega)}^{2}$
$\displaystyle\leq\frac{2}{c_{1}\wedge
1}[\underbrace{\frac{\|w\|_{L^{\infty}(\Omega)}\vee
1}{2}\cdot\inf_{\bar{u}\in\mathcal{N}^{2}}\|\bar{u}-u^{*}\|_{H^{1}(\Omega)}^{2}}_{\mathcal{E}_{app}}+\underbrace{2\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)|}_{\mathcal{E}_{sta}}+\underbrace{\widehat{\mathcal{L}}(u_{\phi_{\mathcal{A}}})-\widehat{\mathcal{L}}(\widehat{u}_{\phi})}_{\mathcal{E}_{opt}}].$

###### Proof.

For any $\bar{u}\in\mathcal{N}^{2}$, we have

$\begin{split}&\mathcal{L}\left(u_{\phi_{\mathcal{A}}}\right)-\mathcal{L}\left(u^{*}\right)\\\
&=\mathcal{L}\left(u_{\phi_{\mathcal{A}}}\right)-\widehat{\mathcal{L}}\left(u_{\phi_{\mathcal{A}}}\right)+\widehat{\mathcal{L}}\left(u_{\phi_{\mathcal{A}}}\right)-\widehat{\mathcal{L}}\left(\widehat{u}_{\phi}\right)+\widehat{\mathcal{L}}\left(\widehat{u}_{\phi}\right)-\widehat{\mathcal{L}}\left(\bar{u}\right)\\\
&\quad+\widehat{\mathcal{L}}\left(\bar{u}\right)-\mathcal{L}\left(\bar{u}\right)+\mathcal{L}\left(\bar{u}\right)-\mathcal{L}\left(u^{*}\right)\\\
&\leq\left[\mathcal{L}\left(\bar{u}\right)-\mathcal{L}\left(u^{*}\right)\right]+2\sup_{u\in\mathcal{N}^{2}}\left|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)\right|+\left[\widehat{\mathcal{L}}\left(u_{\phi_{\mathcal{A}}}\right)-\widehat{\mathcal{L}}\left(\widehat{u}_{\phi}\right)\right],\end{split}$

where the last step is due to the fact that
$\widehat{\mathcal{L}}\left(\widehat{u}_{\phi}\right)-\widehat{\mathcal{L}}\left(\bar{u}\right)\leq
0$. Since $\bar{u}$ can be any element in $\mathcal{N}^{2}$, we take the
infimum of $\bar{u}$ on both side of the above display,

$\displaystyle\mathcal{L}\left(u_{\phi_{\mathcal{A}}}\right)-\mathcal{L}\left(u^{*}\right)$
$\displaystyle\leq\inf_{\bar{u}\in\mathcal{N}^{2}}\left[\mathcal{L}\left(\bar{u}\right)-\mathcal{L}\left(u^{*}\right)\right]+2\sup_{u\in\mathcal{N}^{2}}\left|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)\right|$
$\displaystyle+\left[\widehat{\mathcal{L}}\left(u_{\phi_{\mathcal{A}}}\right)-\widehat{\mathcal{L}}\left(\widehat{u}_{\phi}\right)\right].$
(5)

Now for any $u\in\mathcal{N}^{2}$, set $v=u-u^{*}$, then

$\begin{split}&\mathcal{L}\left(u\right)=\mathcal{L}\left(u^{*}+v\right)\\\
&=\frac{1}{2}(\nabla(u^{*}+v),\nabla(u^{*}+v))_{L^{2}(\Omega)}+\frac{1}{2}(u^{*}+v,u^{*}+v)_{L^{2}(\Omega;w)}\\\
&\quad\ -\langle
u^{*}+v,f\rangle_{L^{2}({\Omega})}-\langle{Tu^{*}+Tv},g\rangle_{L^{2}({\partial\Omega})}\\\
&=\frac{1}{2}(\nabla u^{*},\nabla
u^{*})_{L^{2}(\Omega)}+\frac{1}{2}(u^{*},u^{*})_{L^{2}(\Omega;w)}-\langle
u^{*},f\rangle_{L^{2}({\Omega})}-\langle{Tu^{*}},g\rangle_{L^{2}({\partial\Omega})}\\\
&\quad\ +\frac{1}{2}(\nabla v,\nabla
v)_{L^{2}(\Omega)}+\frac{1}{2}(v,v)_{L^{2}(\Omega;w)}\\\ &\quad\
+\left[(\nabla u^{*},\nabla
v)_{L^{2}(\Omega)}+(u^{*},v)_{L^{2}(\Omega;w)}-\langle
v,f\rangle_{L^{2}({\Omega})}-\langle{Tv},g\rangle_{L^{2}({\partial\Omega})}\right]\\\
&=\mathcal{L}\left(u^{*}\right)+\frac{1}{2}(\nabla v,\nabla
v)_{L^{2}(\Omega)}+\frac{1}{2}(v,v)_{L^{2}(\Omega;w)},\end{split}$

where the last equality is due to the fact that $u^{*}$ is the weak solution
of equation $(\ref{mse})$. Hence

$\displaystyle\frac{c_{1}\wedge
1}{2}\|v\|_{H^{1}(\Omega)}^{2}\leq\mathcal{L}\left(u\right)-\mathcal{L}\left(u^{*}\right)$
$\displaystyle=\frac{1}{2}(\nabla v,\nabla
v)_{L^{2}(\Omega)}+\frac{1}{2}(v,v)_{L^{2}(\Omega;w)}$
$\displaystyle\leq\frac{\|w\|_{L^{\infty}(\Omega)}\vee
1}{2}\|v\|_{H^{1}(\Omega)}^{2},$

that is,

$\frac{c_{1}\wedge
1}{2}\|u-u^{*}\|_{H^{1}(\Omega)}^{2}\leq\mathcal{L}\left(u\right)-\mathcal{L}\left(u^{*}\right)\leq\frac{\|w\|_{L^{\infty}(\Omega)}\vee
1}{2}\|u-u^{*}\|_{H^{1}(\Omega)}^{2}.$ (6)

Combining $(\ref{errdec1})$ and $(\ref{errdec2})$ yields

$\displaystyle\|u_{\phi_{\mathcal{A}}}-u^{*}\|_{H^{1}(\Omega)}^{2}$
$\displaystyle\leq\frac{2}{c_{1}\wedge 1}[\frac{\|w\|_{L^{\infty}(\Omega)}\vee
1}{2}\cdot\inf_{\bar{u}\in\mathcal{N}^{2}}\|\bar{u}-u^{*}\|_{H^{1}(\Omega)}^{2}$
$\displaystyle+2\sup_{u\in\mathcal{N}^{2}}\left|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)\right|+\left[\widehat{\mathcal{L}}\left(u_{\phi_{\mathcal{A}}}\right)-\widehat{\mathcal{L}}\left(\widehat{u}_{\phi}\right)\right]].$

∎

The approximation error $\mathcal{E}_{app}$ describes the expressive power of
the $\mathrm{ReLU}^{2}$ networks $\mathcal{N}^{2}$ in $H^{1}$ norm, which
corresponds to the approximation error in FEM known as the C$\acute{e}$a’s
lemma [7]. The statistical error $\mathcal{E}_{sta}$ is caused by the monte
carlo discritization of $\mathcal{L}(\cdot)$ defined in (2) with
$\widehat{\mathcal{L}}(\cdot)$ in (2). While, the optimization error
$\mathcal{E}_{opt}$ indicates the perforce of the solver $\mathcal{A}$ we
utilized. In contrast, this error is is corresponding to the error of solving
the linear systems in FEM. In this paper we focus on the first tow errors,
i.e, considering the scenario of perfect training with $\mathcal{E}_{opt}=0$.

###  3.1 Approximation error

The current literature on network approximation theory are mainly focus on the
$L^{p},p\in[1,+\infty]$ norm for deep networks [50, 43, 36, 28, 33, 23, 39].
The approximation error of ReLU network in Sobolev norm are considered in [15,
12]. However, the ReLU network may not be suitable for solving PDEs since the
term $\nabla u_{\phi}$ in the loss function will become piece-wise constant
with respect to $\phi$, which will prohibit using SGD for training. In this
section we derive an upper bound on the approximation error
$\mathrm{ReLU}^{2}$ networks $\mathcal{N}^{2}$ in $H^{1}$ norm, which is of
independent interest.

######  Theorem 3.1.

Assume $\|u^{*}\|_{H^{2}(\Omega)}\leq c_{2}$, then there exist an
$\mathrm{ReLU}^{2}$ network $\bar{u}_{\bar{\phi}}\in\mathcal{N}^{2}$ with

$\mathcal{D}\leq\lceil\log_{2}d\rceil+3,\ \ \mathcal{W}\leq
4d\left\lceil\frac{Cc_{2}}{\epsilon}-4\right\rceil^{d}$

such that

$\mathcal{E}_{app}\leq\frac{\|w\|_{L^{\infty}(\Omega)}\vee
1}{2}\|u^{*}-\bar{u}_{\bar{\phi}}\|_{H^{1}\left(\Omega\right)}^{2}\leq\frac{\|w\|_{L^{\infty}(\Omega)}\vee
1}{2}\epsilon^{2},$

where $C$ is a constant depending only on $d$.

###### Proof.

Our proof is based on some classical approximation results of B-splines [35,
8]. Let us recall some notation and useful results. We denote by $\pi_{l}$ the
dyadic partition of $[0,1]$, i.e.,

$\pi_{l}:t_{0}^{(l)}=0<t_{1}^{(l)}<\cdots<t_{2^{l}-1}^{(l)}<t_{2^{l}}^{(l)}=1,$

where $t_{i}^{(l)}=i\cdot 2^{-l}(0\leq i\leq 2^{l})$. The cardinal B-spline of
order $3$ with respect to partition $\pi_{l}$ is defined by

$N_{l,i}^{(3)}(x)=(-1)^{k}\left[t_{i}^{(l)},\ldots,t_{i+3}^{(l)},(x-t)_{+}^{2}\right]\cdot\left(t_{i+3}^{(l)}-t_{i}^{(l)}\right),\quad
i=-2,\cdots,2^{l}-1$

which can be rewritten in the following equivalent form,

$N_{l,i}^{(3)}(x)=2^{2l-1}\sum_{j=0}^{3}(-1)^{j}\left(\begin{array}[]{l}3\\\
j\end{array}\right)(x-i2^{-l}-j2^{-l})_{+}^{2},\quad i=-2,\cdots,2^{l}-1.$ (7)

The multivariate cardinal B-spline of order $3$ is defined by the product of
univariate cardinal B-splines of order $3$, i.e.,

${N}_{l,\mathbf{i}}^{(3)}(\mathbf{x})=\prod_{j=1}^{d}N_{l,i_{j}}^{(3)}\left(x_{j}\right),\quad\mathbf{i}=\left(i_{1},\ldots,i_{d}\right),-3<i_{j}<2^{l}.$

Denote

$S_{l}^{(3)}([0,1]^{d})=\text{span}\\{N_{l,\mathbf{i}}^{(3)},-3<i_{j}<2^{l},j=1,2,\cdots,d\\}.$

Then, the element $f$ in $S_{l}^{(3)}([0,1]^{d})$ are piecewise polynomial
functions according to to partition $\pi_{l}^{d}$ with each piece being degree
$2$ and in $C^{1}([0,1]^{d})$. Since

$S_{1}^{(3)}\subset S_{2}^{(3)}\subset S_{3}^{(3)}\subset\cdots,$

We can further denote

$S^{(3)}([0,1]^{d})=\bigcup_{l=1}^{\infty}S_{l}^{(3)}([0,1]^{d}).$

The following approximation result of cardinal B-splines in Sobolev spaces
which is a direct consequence of theorem 3.4 in [34] play an important role in
the proof of this Theorem.

######  Lemma 3.

Assume $u^{*}\in H^{2}([0,1]^{d})$, there exists
$\\{c_{j}\\}_{j=1}^{(2^{l}-4)^{d}}\subset\mathbb{R}$ with $l>2$ such that

$\|u^{*}-\sum_{j=1}^{(2^{l}-4)^{d}}c_{j}{N}_{l,\mathbf{i}_{j}}^{(3)}\|_{H^{1}(\Omega)}\leq\frac{C}{2^{l}}\|u^{*}\|_{H^{1}(\Omega)},$

where $C$ is a constant only depend on $d$.

######  Lemma 4.

The multivariate B-spline ${N}_{l,\mathbf{i}}^{(3)}(\mathbf{x})$ can be
implemented exactly by a $\mathrm{ReLU}^{2}$ network with depth
$\lceil\log_{2}d\rceil+2$ and width $4d$.

###### Proof.

Denote

$\sigma(x)=\left\\{\begin{array}[]{ll}x^{2},&x\geq 0\\\
0,&\text{else}\end{array}\right.$

as the activation function in $\mathrm{ReLU}^{2}$ network. By definition of
$N_{l,i}^{(3)}(x)$ in (7), it’s clear that $N_{l,i}^{(3)}(x)$ can be
implemented by $\mathrm{ReLU}^{2}$ network without any error with depth $2$
and width $4$. On the other hand $\mathrm{ReLU}^{2}$ network can also realize
multiplication without any error. In fact, for any $x,y\in\mathbb{R}$,

$xy=\frac{1}{4}[(x+y)^{2}-(x-y)^{2}]=\frac{1}{4}[\sigma(x+y)+\sigma(-x-y)-\sigma(x-y)-\sigma(y-x)].$

Hence multivariate B-spline of order $3$ can be implemented by
$\mathrm{ReLU}^{2}$ network exactly with depth $\lceil\log_{2}d\rceil+2$ and
width $4d$. ∎

For any $\epsilon>0$, by Lemma 3 and 4 with
$\frac{1}{2^{l}}\leq\left\lceil\frac{C\|u^{*}\|_{H^{2}}}{\epsilon}\right\rceil$,
there exists $\bar{u}_{\bar{\phi}}\in\mathcal{N}^{2}$, such that

$\left\|u^{*}-\bar{u}_{\bar{\phi}}\right\|_{H^{1}(\Omega)}\leq\epsilon.$ (8)

The depth $\mathcal{D}$ and width $\mathcal{W}$ of $\bar{u}_{\bar{\phi}}$ are
satisfying $\mathcal{D}\leq\lceil\log_{2}d\rceil+3$ and $\mathcal{W}\leq
4dn=4d\left\lceil\frac{C\|u^{*}\|_{H^{2}}}{\epsilon}-4\right\rceil^{d}$,
respectively. ∎

###  3.2 Statistical error

In this section, we bound the statistical error

${\mathcal{E}_{sta}}=2\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)|.$

######  Lemma 5.

$\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)|\leq\sum_{j=1}^{4}\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}_{j}(u)-\widehat{\mathcal{L}_{j}}(u)|,$

where,

$\mathcal{L}_{1}(u)=|\Omega|\mathbb{E}_{{X\sim\mathrm{U}(\Omega)}}[w(X)u^{2}(X)/2],\
\
\widehat{\mathcal{L}_{1}}(u)=\frac{|\Omega|}{N}\sum_{i=1}^{N}[\frac{w(X_{i})u^{2}(X_{i})}{2}],$
$\mathcal{L}_{2}(u)=|\Omega|\mathbb{E}_{{X\sim\mathrm{U}(\Omega)}}[u(X)f(X)],\
\
\widehat{\mathcal{L}_{2}}(u)=\frac{|\Omega|}{N}\sum_{i=1}^{N}[u(X_{i})f(X_{i})],$
$\mathcal{L}_{3}(u)=|\partial{\Omega}|\mathbb{E}_{{Y\sim\mathrm{U}(\partial\Omega)}}[{Tu}(Y)g(Y)],\
\
\widehat{\mathcal{L}_{3}}(u)=\frac{\partial\Omega}{M}\sum_{j=1}^{M}[u(Y_{j})g(Y_{j})],$
$\mathcal{L}_{4}(u)=|\Omega|\mathbb{E}_{{X\sim\mathrm{U}(\Omega)}}[\|\nabla
u(X)\|_{2}^{2}/2],\ \
\widehat{\mathcal{L}_{4}}(u)=\frac{|\Omega|}{N}\sum_{i=1}^{N}[\frac{\|\nabla
u(X_{i})\|_{2}^{2}}{2}].$

###### Proof.

This Lemma holds by the direct consequence of triangle inequality. ∎

By Lemma 5, we have to bound the the maximum value of four random processes
indexed by $u\in\mathcal{N}^{2}$. To this end, we recall tools in empirical
process [45, 20]. Denote

$\mathcal{B}=\max\\{\|\bar{u}_{\bar{\phi}}\|_{L^{\infty}(\Omega)},\|\|\nabla\bar{u}_{\bar{\phi}}\|_{2}^{2}\|_{L^{\infty}(\Omega)}\\},$
(9)

where $\bar{u}_{\bar{\phi}}$ is the best approximation of $u^{*}$ in Theorem
3.1. Let $f,g,w$ be bounded, say by some constant $c_{3}$, i.e, we assume that

$\|f\|_{L^{\infty}(\Omega)}\vee\|w\|_{L^{\infty}(\Omega)}\vee\|g\|_{L^{\infty}(\partial\Omega)}\vee\mathcal{B}\leq
c_{3}<\infty.$

We use $\mu$ to denote $\mathrm{U}(\Omega)(\mathrm{U}(\partial\Omega))$. Given
$n=N(M)$ i.i.d samples $\mathbf{Z}_{n}=\\{Z_{i}\\}_{i=1}^{n}$ from $\mu$, with
$Z_{i}=X_{i}(Y_{i})\sim\mu$, we need the following Rademacher complexity to
measure the capacity of the given function class $\mathcal{N}$ restricted on
$n$ random samples $\mathbf{Z}_{n}$.

######  Definition 6.

The Rademacher complexity of a set $A\subseteq\mathrm{R}^{n}$ is defined as

$\mathfrak{R}(A)=\mathbb{E}_{\mathbf{Z}_{n},\Sigma_{n}}[\sup_{a\in
A}\frac{1}{n}|\sum_{i}\sigma_{i}a_{i}|],$

where, $\Sigma_{n}=\\{\sigma_{i}\\}_{i=1}^{n}$ are $n$ i.i.d Rademacher
variables with
$\mathbb{P}(\sigma_{i}=1)=\mathbb{P}(\sigma_{i}=-1)=\frac{1}{2}.$ The
Rademacher complexity of function class $\mathcal{N}$ associate with random
sample $\mathbf{Z}_{n}$ is defined as as

$\mathfrak{R}(\mathcal{N})=\mathbb{E}_{\mathbf{Z}_{n},\Sigma_{n}}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i}\sigma_{i}u(Z_{i})|].$

######  Lemma 7.

Let $\Psi_{1}(x,y)=\frac{w(x)y^{2}}{2}:\mathbb{R}^{d}\times\mathbb{R},|y|\leq
c_{3}$, $\Psi_{2}(x,y)=f(x)y:\mathbb{R}^{d}\times\mathbb{R},|y|\leq c_{3}$,
$\Psi_{3}(x,y)=g(x)y:\mathbb{R}^{d}\times\mathbb{R},|y|\leq c_{3}$. Then
$\Psi_{1}(x,y)$, $\Psi_{2}(x,y)$ and $\Psi_{3}(x,y)$ are $c_{3}^{2}$, $c_{3}$
and $c_{3}$-Lipschitz continuous on $y$ for all $x$ and
$\Psi_{j}(x,0)=0,j=1,2,3$.

###### Proof.

We give the proof for $\Psi_{1}$ and omit the details for $\Psi_{2},\Psi_{3}$
since they can be shown similarly. For arbitrary $y_{1},y_{2}$ with
$|y_{i}|\leq c_{3},i=1,2$,

$\displaystyle|\Psi_{1}(x,y_{1})-\Psi_{1}(x,y_{2})|$
$\displaystyle=|\frac{w(x)y_{1}^{2}}{2}-\frac{w(x)y_{2}^{2}}{2}|$
$\displaystyle=\frac{|w(x)(y_{1}+y_{2})|}{2}|y_{1}-y_{2}|\leq
c_{3}^{2}|y_{1}(x)-y_{2}|.$

∎

By Corollary 3.17 in [20] and Lemma 7, we have the following Lipschitz
contraction results on Rademacher complexity.

######  Lemma 8.

Let $\mathcal{N}=\\{u(x):\|u\|_{L^{\infty}(\Omega)}\leq c_{3}\\}$. Define,

$\Psi_{j}\circ\mathcal{N}=\\{\mathrm{composition\ \ of}\ \ \Psi_{j}\ \
\mathrm{and}\ \
\mathcal{N}:x\rightarrow\Psi_{j}(x,u(x)):u\in\mathcal{N}\\},j=1,2,3.$

Then, $\mathfrak{R}(\Psi_{1}\circ\mathcal{N})\leq
2c_{3}^{2}\mathfrak{R}(\mathcal{N}),$
$\mathfrak{R}(\Psi_{i}\circ\mathcal{N})\leq
2c_{3}\mathfrak{R}(\mathcal{N}),i=2,3$.

The following symmetrization result shows that the Rademacher complexity
$\mathfrak{R}(\Psi_{j}\circ\mathcal{N}^{2})$ gives upper bound on
$\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}_{j}(u)-\widehat{\mathcal{L}_{j}}(u)|,$
$j=1,...,3.$

######  Lemma 9.

$\mathbb{E}_{\mathbf{Z}_{n}}[\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}_{j}(u)-\widehat{\mathcal{L}_{j}}(u)|]\leq\mathfrak{R}(\Psi_{j}\circ\mathcal{N}^{2})$,
$j=1,...,3.$

###### Proof.

We give the proof for $j=1$ and omit the proof for $\Psi_{2}$ and $\Psi_{3}$
since they can be shown similarly. Let
$\tilde{\mathbf{Z}}_{n}=\\{\tilde{Z}_{i}\\}_{i=1}^{n}$ be an i.i.d ghost
sample from $\mu$ and $\tilde{\mathbf{Z}}_{n}$ is independent of
${\mathbf{Z}}_{n}$.

$\displaystyle\mathbb{E}_{\mathbf{Z}_{n}}[\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}_{1}(u)-\widehat{\mathcal{L}_{1}}(u)|]$
$\displaystyle=\frac{|\Omega|}{n}\mathbb{E}_{\mathbf{Z}_{n}}[\sup_{u\in\mathcal{N}^{2}}|\sum_{i=1}^{n}(\mathbb{E}_{\tilde{\mathbf{Z}}_{n}}[\frac{w(\tilde{Z}_{i})u^{2}(\tilde{Z}_{i})}{2}]-\frac{w(Z_{i})u^{2}(Z_{i})}{2})|]$
$\displaystyle\leq\frac{|\Omega|}{n}\mathbb{E}_{\mathbf{Z}_{n},\tilde{\mathbf{Z}}_{n}}[\sup_{u\in\mathcal{N}^{2}}|\sum_{i=1}^{n}(\frac{w(\tilde{Z}_{i})u^{2}(\tilde{Z}_{i})}{2}-\frac{w(Z_{i})u^{2}(Z_{i})}{2})|]$
$\displaystyle=\frac{|\Omega|}{n}\mathbb{E}_{\mathbf{Z}_{n},\tilde{\mathbf{Z}}_{n},\Sigma_{n}}[\sup_{u\in\mathcal{N}^{2}}|\sum_{i=1}^{n}\sigma_{i}(\frac{w(\tilde{Z}_{i})u^{2}(\tilde{Z}_{i})}{2}-\frac{w(Z_{i})u^{2}(Z_{i})}{2})|]$
$\displaystyle\leq\frac{|\Omega|}{n}\mathbb{E}_{\mathbf{Z}_{n},\Sigma}[\sup_{u\in\mathcal{N}^{2}}|\sum_{i=1}^{n}\frac{w({Z}_{i})u^{2}({Z}_{i})}{2}|]+\frac{|\Omega|}{n}\mathbb{E}_{\tilde{\mathbf{Z}}_{n},\Sigma}[\sup_{u\in\mathcal{N}^{2}}|\sum_{i=1}^{n}\frac{w(\tilde{Z}_{i})u^{2}(\tilde{Z}_{i})}{2}|]$
$\displaystyle=2|\Omega|\mathfrak{R}(\Psi_{1}\circ\mathcal{N}^{2})$

where, the first inequality follows from the Jensen’s inequality, and the
second equality holds since both
$\sigma_{i}(\frac{w(\tilde{Z}_{i})u^{2}(\tilde{Z}_{i})}{2}-\frac{w({Z}_{i})u^{2}({Z}_{i})}{2})$
and
$(\frac{w(\tilde{Z}_{i})u^{2}(\tilde{Z}_{i})}{2}-\frac{w({Z}_{i})u^{2}({Z}_{i})}{2})$
are governed by the same law, and the last equality holds since the
distribution of the two terms are the same. ∎

Next we give a upper bound of $\mathfrak{R}(\mathcal{N}^{2})$ in terms of the
covering number of $\mathcal{N}^{2}$ by using the Dudley’s entropy formula
[11].

######  Definition 10.

Suppose that $W\subset\mathbb{R}^{n}$. For any $\epsilon>0$, let
$V\subset\mathbb{R}^{n}$ be a $\epsilon$-cover of $W$ with respect to the
distance $d_{\infty}$, that is, for any $w\in W$, there exists a $v\in V$ such
that $d_{\infty}(u,v)<\epsilon$, where $d_{\infty}$ is defined by

$d_{\infty}(u,v):=\|u-v\|_{{\infty}}.$

The covering number $\mathcal{C}(\epsilon,W,d_{\infty})$ is defined to be the
minimum cardinality among all $\epsilon$-cover of $W$ with respect to the
distance $d_{\infty}$.

######  Definition 11.

Suppose that $\mathcal{N}$ is a class of functions from $\Omega$ to
$\mathbb{R}$. Given $n$ sample
$\mathbf{Z}_{n}=(Z_{1},Z_{2},\cdots,Z_{n})\in\Omega^{n}$,
$\mathcal{N}|_{\mathbf{Z}_{n}}\subset\mathbb{R}^{n}$ is defined by

$\mathcal{N}|_{\mathbf{Z}_{n}}=\\{(u(Z_{1}),u(Z_{2}),\cdots,u(Z_{n})):u\in\mathcal{N}\\}.$

The uniform covering number $\mathcal{C}_{\infty}(\epsilon,\mathcal{N},n)$ is
defined by

$\mathcal{C}_{\infty}(\epsilon,\mathcal{N},n)=\max_{\mathbf{Z}_{n}\in\Omega^{n}}\mathcal{C}(\epsilon,\mathcal{N}|_{\mathbf{Z}_{n}},d_{\infty})$

######  Lemma 12.

Assume $0\in\mathcal{N}$ and the diameter of $\mathcal{N}$ is less than
$\mathcal{B}$, i.e., $\|u\|_{L^{\infty}(\Omega)}\leq\mathcal{B},\forall
u\in\mathcal{N}$. Then

$\mathfrak{R}(\mathcal{N})\leq\inf_{0<\delta<\mathcal{B}}\left(4\delta+\frac{12}{\sqrt{n}}\int_{\delta}^{\mathcal{B}}\sqrt{\log(2\mathcal{C}\left(\epsilon,\mathcal{N},n\right))}\mathrm{d}\epsilon\right).$

###### Proof.

The proof follows from the chaining method [45]. We need the Massart’s finite
class Lemma [4].

######  Lemma 13.

For any finite set $V\in\mathbb{R}^{n}$ with diameter $D=\sup_{v\in
V}\|v\|_{2}$, then

$\mathbb{E}_{\Sigma_{n}}[\sup_{v\in
V}\frac{1}{n}|\sum_{i}\sigma_{i}v_{i}|]\leq\frac{D}{n}\sqrt{2\log(2|V|)}.$

By definition

$\mathfrak{R}(\mathcal{N})=\mathfrak{R}(\mathcal{N}|_{\mathbf{Z}_{n}})=\mathbb{E}_{\mathbb{Z}_{n}}[\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i}\sigma_{i}u(Z_{i})||\mathbf{Z}_{n}]].$

Thus, it suffice to show

$\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i}\sigma_{i}u(Z_{i})|]\leq\inf_{0<\delta<\mathcal{B}}\left(4\delta+\frac{12}{\sqrt{n}}\int_{\delta}^{\mathcal{B}}\sqrt{\log\mathcal{C}\left(\epsilon,\mathcal{N}^{2},n\right)}\mathrm{d}\epsilon\right)$

by conditioning on $\mathbf{Z}_{n}$. Given an positive integer $K$, let
$\epsilon_{k}=2^{-k+1}\mathcal{B}$, $k=1,...K$. Let $C_{k}$ be a cover of
$\mathcal{N}|_{\mathbf{Z}_{n}}\subseteq\mathbb{R}^{n}$ whose covering number
is denoted as
$\mathcal{C}(\epsilon_{k},\mathcal{N}|_{\mathbf{Z}_{n}},d_{\infty})$. Then, by
definition, $\forall u\in\mathcal{N},$ there $\exists$ $c^{k}\in C_{k}$ such
that

$d_{\infty}(u|_{\mathbf{Z}_{n}},c^{k})=\max\\{|u(Z_{i})-c^{k}_{i}|,i=1,...,n\\}\leq\epsilon_{k},k=1,...,K.$

Moreover, we denote the best approximate element of $u$ in $C_{k}$ with
respect to $d_{\infty}$ as $c^{k}(u)$. Then,

$\displaystyle\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}u(Z_{i})|]$
$\displaystyle=\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}(u(Z_{i})-c^{K}_{i}(u))+\sum_{j=1}^{K-1}\sum_{i=1}^{n}\sigma_{i}(c^{j}_{i}(u)-c^{j+1}_{i}(u))+\sum_{i=1}^{n}\sigma_{i}c^{1}_{i}(u)|]$
$\displaystyle\leq\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}(u(Z_{i})-c^{K}_{i}(u))|]+\sum_{j=1}^{K-1}\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}(c^{j}_{i}(u)-c^{j+1}_{i}(u))|]$
$\displaystyle+\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}c^{1}_{i}(u)|].$

Since $0\in\mathcal{N}$, and the diameter of $\mathcal{N}$ is smaller than
$\mathcal{B}$, we can choose $C_{1}=\\{0\\}$ such that the third term in the
above display vanishes. By Hölder’s inequality, we deduce that the first term
can be bounded by $\epsilon_{K}$ as follows.

$\displaystyle\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}(u(Z_{i})-c^{K}_{i}(u))|]$
$\displaystyle\leq\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}(\sum_{i=1}^{n}|\sigma_{i}|)(\sum_{i=1}^{n}\max_{i=1,...,n}\\{|u(Z_{i})-c^{K}_{i}(u)|\\})]$
$\displaystyle\leq\epsilon_{K}.$

Let $V_{j}=\\{c^{j}(u)-c^{j+1}(u):u\in\mathcal{N}\\}$. Then by definition, the
number of elements in $V_{j}$ and $C_{j}$ satisfying

$|V_{j}|\leq|C_{j}||C_{j+1}|\leq|C_{j+1}|^{2}.$

And the diameter of $V_{j}$ denoted as $D_{j}$ can be bounded as

$\displaystyle D_{j}=\sup_{v\in
V_{j}}\|v\|_{2}\leq\sqrt{n}\sup_{u\in\mathcal{N}}\|c^{j}(u)-c^{j+1}(u)\|_{\infty}$
$\displaystyle\leq\sqrt{n}\sup_{u\in\mathcal{N}}\|c^{j}(u)-u\|_{\infty}+\|u-c^{j+1}(u)\|_{\infty}$
$\displaystyle\leq\sqrt{n}(\epsilon_{j}+\epsilon_{j+1})$ $\displaystyle\leq
3\sqrt{n}\epsilon_{j+1}.$

Then,

$\displaystyle\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{j=1}^{K-1}\sum_{i=1}^{n}\sigma_{i}(c^{j}_{i}(u)-c^{j+1}_{i}(u))|]\leq\sum_{j=1}^{K-1}\mathbb{E}_{\Sigma}[\sup_{v\in
V_{j}}\frac{1}{n}|\sum_{i=1}^{n}\sigma_{i}v_{j}|]$
$\displaystyle\leq\sum_{j=1}^{K-1}\frac{D_{j}}{n}\sqrt{2\log(2|V_{j}|)}$
$\displaystyle\leq\sum_{j=1}^{K-1}\frac{6\epsilon_{j+1}}{\sqrt{n}}\sqrt{\log(2|C_{j+1}|)},$

where we use triangle inequality in the first inequality, and use Lemma 13 in
the second inequality. Putting all the above estimates together, we get

$\displaystyle\mathbb{E}_{\Sigma}[\sup_{u\in\mathcal{N}}\frac{1}{n}|\sum_{i}\sigma_{i}u(Z_{i})|]\leq\epsilon_{K}+\sum_{j=1}^{K-1}\frac{6\epsilon_{j+1}}{\sqrt{n}}\sqrt{\log(2|C_{j+1}|)}$
$\displaystyle\leq\epsilon_{K}+\sum_{j=1}^{K}\frac{12(\epsilon_{j}-\epsilon_{j+1})}{\sqrt{n}}\sqrt{\log(2\mathcal{C}\left(\epsilon_{j},\mathcal{N},n\right))}$
$\displaystyle\leq\epsilon_{K}+\frac{12}{\sqrt{n}}\int_{\epsilon_{K+1}}^{\mathcal{B}}\sqrt{\log(2\mathcal{C}\left(\epsilon,\mathcal{N},n\right))}\mathrm{d}\epsilon$
$\displaystyle\leq\inf_{0<\delta<\mathcal{B}}(4\delta+\frac{12}{\sqrt{n}}\int_{\delta}^{\mathcal{B}}\sqrt{\log(2\mathcal{C}\left(\epsilon,\mathcal{N},n\right))}\mathrm{d}\epsilon).$

where, last inequality holds since for $0<\delta<\mathcal{B}$, we can choose
$K$ to be the largest integer such that $\epsilon_{K+1}>\delta$, at this time
$\epsilon_{K}\leq 4\epsilon_{K+2}\leq 4\delta.$ ∎

Now we turn to handle the most difficult term
$\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}_{4}(u)-\widehat{\mathcal{L}_{4}}(u)|$,
where we need to bound the Rademacher complexity of the non-Lipschitz
composition of gradient norm and $\mathrm{ReLU}^{2}$ network. We believe that
the technique we used here is helpful for bounding the statistical errors for
other deep PDEs solvers where the main difficulties is bounding the Rademacher
complexity of non-Lipschitz composition induced by the gradient operator.

######  Lemma 14.

$\displaystyle\mathbb{E}_{\mathbf{Z}_{n}}[\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}_{4}(u)-\widehat{\mathcal{L}_{4}}(u)|]$
$\displaystyle\leq\mathbb{E}_{\mathbf{Z}_{n},\Sigma_{n}}[\sup_{u\in\mathcal{N}^{2}}\frac{1}{n}|\sum_{i}\sigma_{i}\|\nabla
u(Z_{i})\|^{2}|]$ (10)
$\displaystyle\leq\mathfrak{R}(\mathcal{N}^{1,2})=\mathbb{E}_{\mathbf{Z}_{n},\Sigma_{n}}[\sup_{u\in\mathcal{N}^{1,2}}\frac{1}{n}|\sum_{i}\sigma_{i}u(Z_{i})|]$
(11)
$\displaystyle\leq\inf_{0<\delta<\mathcal{B}}\left(4\delta+\frac{12}{\sqrt{n}}\int_{\delta}^{\mathcal{B}}\sqrt{\log(2\mathcal{C}\left(\epsilon,\mathcal{N}^{1,2},n\right))}\mathrm{d}\epsilon\right).$
(12)

###### Proof.

The proof of (10) is based on the symmetrization method used in the proof of
Lemma 9, we omit the detail here. The proof of (11) is a direct consequence of
the following claim. Claim: Let $u$ be a function implemented by a
$\mathrm{ReLU}^{2}$ network with depth $\mathcal{D}$ and width $\mathcal{W}$.
Then $\|\nabla u\|_{2}^{2}$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network with depth $\mathcal{D}+3$ and
width $d\left(\mathcal{D}+2\right)\mathcal{W}$.

Denote $\mathrm{ReLU}$ and $\mathrm{ReLU}^{2}$ as $\sigma_{1}$ and
$\sigma_{2}$, respectively. As long as we show that each partial derivative
$D_{i}u(i=1,2,\cdots,d)$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network respectively, we can easily obtain
the network we desire, since, $\|\nabla
u\|_{2}^{2}=\sum_{i=1}^{d}\left|D_{i}u\right|^{2}$ and the square function can
be implemented by $x^{2}=\sigma_{2}(x)+\sigma_{2}(-x)$.

Now we show that for any $i=1,2,\cdots,d$, $D_{i}u$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network. We deal with the first two layers
in details since there are a little bit difference for the first two layer and
apply induction for layers $k\geq 3$. For the first layer, since
$\sigma_{2}^{{}^{\prime}}(x)=2\sigma_{1}(x)$, we have for any
$q=1,2\cdots,n_{1}$

$D_{i}u_{q}^{(1)}=D_{i}\sigma_{2}\left(\sum_{j=1}^{d}a_{qj}^{(1)}x_{j}+b_{q}^{(1)}\right)=2\sigma_{1}\left(\sum_{j=1}^{d}a_{qj}^{(1)}x_{j}+b_{q}^{(1)}\right)\cdot
a_{qi}^{(1)}$

Hence $D_{i}u_{q}^{(1)}$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network with depth $2$ and width $1$. For
the second layer,

$D_{i}u_{q}^{(2)}=D_{i}\sigma_{2}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}u_{j}^{(1)}+b_{q}^{(2)}\right)=2\sigma_{1}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}u_{j}^{(1)}+b_{q}^{(2)}\right)\cdot\sum_{j=1}^{n_{1}}a_{qj}^{(2)}D_{i}u_{j}^{(1)}$

Since
$\sigma_{1}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}u_{j}^{(1)}+b_{q}^{(2)}\right)$
and $\sum_{j=1}^{n_{1}}a_{qj}^{(2)}D_{i}u_{j}^{(1)}$ can be implemented by two
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ `subnetworks`, respectively, and the
multiplication can also be implemented by

$\begin{split}x\cdot y&=\frac{1}{4}\left[(x+y)^{2}-(x-y)^{2}\right]\\\
&=\frac{1}{4}\left[\sigma_{2}(x+y)+\sigma_{2}(-x-y)-\sigma_{2}(x-y)-\sigma_{2}(-x+y)\right],\end{split}$

we conclude that $D_{i}u_{q}^{(2)}$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network. We have

$\mathcal{D}\left(\sigma_{1}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}u_{j}^{(1)}+b_{q}^{(2)}\right)\right)=3,\mathcal{W}\left(\sigma_{1}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}u_{j}^{(1)}+b_{q}^{(2)}\right)\right)\leq\mathcal{W}$

and

$\mathcal{D}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}D_{i}u_{j}^{(1)}\right)=2,\mathcal{W}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(2)}D_{i}u_{j}^{(1)}\right)\leq\mathcal{W}.$

Thus $\mathcal{D}\left(D_{i}u_{q}^{(2)}\right)=4,$
$\mathcal{W}\left(D_{i}u_{q}^{(2)}\right)\leq\max\\{2\mathcal{W},4\\}$.

Now we apply induction for layers $k\geq 3$. For the third layer,

$D_{i}u_{q}^{(3)}=D_{i}\sigma_{2}\left(\sum_{j=1}^{n_{2}}a_{qj}^{(3)}u_{j}^{(2)}+b_{q}^{(3)}\right)=2\sigma_{1}\left(\sum_{j=1}^{n_{2}}a_{qj}^{(3)}u_{j}^{(2)}+b_{q}^{(3)}\right)\cdot\sum_{j=1}^{n_{2}}a_{qj}^{(3)}D_{i}u_{j}^{(2)}.$

Since
$\mathcal{D}\left(\sigma_{1}\left(\sum_{j=1}^{n_{2}}a_{qj}^{(3)}u_{j}^{(2)}+b_{q}^{(3)}\right)\right)=4$,
$\mathcal{W}\left(\sigma_{1}\left(\sum_{j=1}^{n_{2}}a_{qj}^{(3)}u_{j}^{(2)}+b_{q}^{(3)}\right)\right)\leq\mathcal{W}$
and

$\mathcal{D}\left(\sum_{j=1}^{n_{2}}a_{qj}^{(3)}D_{i}u_{j}^{(2)}\right)=4,\mathcal{W}\left(\sum_{j=1}^{n_{1}}a_{qj}^{(3)}D_{i}u_{j}^{(2)}\right)\leq\max\\{2\mathcal{W},4\mathcal{W}\\}=4\mathcal{W},$

we conclude that $D_{i}u_{q}^{(3)}$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network and
$\mathcal{D}\left(D_{i}u_{q}^{(3)}\right)=5$,
$\mathcal{W}\left(D_{i}u_{q}^{(3)}\right)\leq\max\\{5\mathcal{W},4\\}=5\mathcal{W}$.

We assume that $D_{i}u_{q}^{(k)}(q=1,2,\cdots,n_{k})$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network and
$\mathcal{D}\left(D_{i}u_{q}^{(k)}\right)=k+2$,
$\mathcal{W}\left(D_{i}u_{q}^{(3)}\right)\leq(k+2)\mathcal{W}$. For the
$(k+1)-$th layer,

$\displaystyle
D_{i}u_{q}^{(k+1)}=D_{i}\sigma_{2}\left(\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}u_{j}^{(k)}+b_{q}^{(k+1)}\right)$
$\displaystyle=2\sigma_{1}\left(\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}u_{j}^{(k)}+b_{q}^{(k+1)}\right)\cdot\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}D_{i}u_{j}^{(k)}.$

Since
$\mathcal{D}\left(\sigma_{1}\left(\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}u_{j}^{(k)}+b_{q}^{(k+1)}\right)\right)=k+2$,
$\mathcal{W}\left(\sigma_{1}\left(\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}u_{j}^{(k)}+b_{q}^{(k+1)}\right)\right)\leq\mathcal{W}$
and
$\mathcal{D}\left(\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}D_{i}u_{j}^{(k)}\right)=k+2$,
$\mathcal{W}\left(\sum_{j=1}^{n_{k}}a_{qj}^{(k+1)}D_{i}u_{j}^{(k)}\right)\leq\max\\{(k+2)\mathcal{W},4\mathcal{W}\\}=(k+2)\mathcal{W}$,
we conclude that $D_{i}u_{q}^{(k+1)}$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network and
$\mathcal{D}\left(D_{i}u_{q}^{(k+1)}\right)=k+3$,
$\mathcal{W}\left(D_{i}u_{q}^{(k+1)}\right)\leq\max\\{(k+3)\mathcal{W},4\\}=(k+3)\mathcal{W}$.

Hence we derive that $D_{i}u=D_{i}u_{1}^{\mathcal{D}}$ can be implemented by a
$\mathrm{ReLU}$-$\mathrm{ReLU}^{2}$ network and
$\mathcal{D}\left(D_{i}u\right)=\mathcal{D}+2$,
$\mathcal{W}\left(D_{i}u\right)\leq\left(\mathcal{D}+2\right)\mathcal{W}$.
Finally we obtain that $\mathcal{D}\left(\|\nabla
u\|^{2}\right)=\mathcal{D}+3$, $\mathcal{W}\left(\|\nabla u\|^{2}\right)\leq
d\left(\mathcal{D}+2\right)\mathcal{W}$.

The proof of (12) follows from Lemma 12. ∎

By Lemma 9, Lemma 8, Lemma 12 and Lemma 14, we have to find upper bounds for
the converging numbers $\mathcal{C}\left(\epsilon,\mathcal{N}^{2},n\right)$
and $\mathcal{C}\left(\epsilon,\mathcal{N}^{1,2},n\right)$ used in the
Dudley’s entropy formula. To this end, we need the VC-dimension [46] and
Pseudo-dimension [1].

######  Definition 15.

Let $\mathcal{N}$ be a set of functions from $X=\Omega(\partial\Omega)$ to
$\\{0,1\\}$. Suppose that $S=\\{x_{1},x_{2},\cdots,x_{n}\\}\subset X$. We say
that $S$ is shattered by $\mathcal{N}$ if for any $b\in\\{0,1\\}^{n}$, there
exists a $u\in\mathcal{N}$ satisfying

$u(x_{i})=b_{i},\quad i=1,2,\dots,n$

######  Definition 16.

The VC-dimension of $\mathcal{N}$, denoted as $\mathrm{VCdim}(\mathcal{N})$,
is defined to be the maximum cardinality among all sets shattered by
$\mathcal{N}$.

VC-dimension reflects the capability of a class of functions to perform binary
classification of points. The larger VC-dimension is, the stronger the
capability to perform binary classification is. For more discussion of VC-
dimension, readers are referred to [1].

For real-valued functions, we can generalize the concept of VC-dimension into
pseudo-dimension [1].

######  Definition 17.

Let $\mathcal{N}$ be a set of functions from $X$ to $\mathbb{R}$. Suppose that
$S=\\{x_{1},x_{2},\cdots,x_{n}\\}\subset X$. We say that $S$ is pseudo-
shattered by $\mathcal{N}$ if there exists $y_{1},y_{2},\cdots,y_{n}$ such
that for any $b\in\\{0,1\\}^{n}$, there exists a $u\in\mathcal{N}$ satisfying

$\mathrm{sign}(u(x_{i})-y_{i})=b_{i},\quad i=1,2,\dots,n$

and we say that $\\{y_{i}\\}_{i=1}^{n}$ witnesses the shattering.

######  Definition 18.

The pseudo-dimension of $\mathcal{N}$, denoted as
$\mathrm{Pdim}(\mathcal{N})$, is defined to be the maximum cardinality among
all sets pseudo-shattered by $\mathcal{N}$.

The following proposition showing a relation between uniform covering number
and pseudo-dimension.

######  Proposition 19 (Theorem 12.2, [1]).

Let $\mathcal{N}$ be a set of real functions from a domain $X$ to the bounded
interval $[0,\mathcal{B}]$. Let $\epsilon>0$. Then

$\mathcal{C}_{\infty}(\epsilon,\mathcal{N},n)\leq\sum_{i=1}^{\mathrm{Pdim}(\mathcal{N})}\left(\begin{array}[]{c}n\\\
i\end{array}\right)\left(\frac{\mathcal{B}}{\epsilon}\right)^{i},$

which is less than
$\left(\frac{en\mathcal{B}}{\epsilon\cdot\mathrm{Pdim}(\mathcal{N})}\right)^{\mathrm{Pdim}(\mathcal{N})}$
for $n\geq\mathrm{Pdim}(\mathcal{N})$.

We now present the bound of pseudo-dimension for the $\mathcal{N}^{1,2}$, the
class of network functions with $\mathrm{ReLU}$ and $\mathrm{ReLU}^{2}$
activation functions. We first need a lemma stated below.

######  Lemma 20.

Let $p_{1},\cdots,p_{m}$ be polynomials with $n$ variables of degree at most
$d$. If $n\leq m$, then

$|\\{(\mathrm{sign}(p_{1}(x)),\cdots,\mathrm{sign}(p_{m}(x))):x\in\mathbb{R}^{n}\\}|\leq
2\left(\frac{2emd}{n}\right)^{n}$

###### Proof.

See Theorem 8.3 in [1]. ∎

######  Theorem 3.2.

Let

$\begin{split}\mathcal{N}:=\\{&u\in[0,1]^{d}:u\text{ can be implemented by a
neural network}\\\ &\text{with depth no more than }\mathcal{D}\text{ and width
no more than }\mathcal{W},\\\ &\text{and activation function in each unit be
the ReLU or the}\ \ \mathrm{ReLU}^{2}.\\}\end{split}$

Then

$\mathrm{Pdim}(\mathcal{N})=\mathcal{O}(\mathcal{D}^{2}\mathcal{W}^{2}(\mathcal{D}+\log\mathcal{W})).$

###### Proof.

The argument is follows from the proof of Theorem 6 in [2]. The result stated
here is somewhat stronger then Theorem 6 in [2] since
$\mathrm{VCdim}(\mathrm{sign}(\mathcal{N}))\leq\mathrm{Pdim}(\mathcal{N})$.

We consider a new set of functions:

$\mathcal{\widetilde{N}}=\\{\widetilde{u}(x,y)=\mathrm{sign}(u(x)-y):u\in\mathcal{H}\\}$

It is clear that
$\mathrm{Pdim}(\mathcal{N})\leq\mathrm{VCdim}(\mathcal{\widetilde{N}})$. We
now bound the VC-dimension of $\mathcal{\widetilde{N}}$. Denoting
$\mathcal{M}$ as the total number of parameters(weights and biases) in the
neural network implementing functions in $\mathcal{N}$, in our case we want to
derive the uniform bound for

$K_{\\{x_{i}\\},\\{y_{i}\\}}(m):=|\\{(\operatorname{sign}(f(x_{1},a)-y_{1}),\ldots,\operatorname{sign}(u(x_{m},a)-y_{m})):a\in\mathbb{R}^{\mathcal{M}}\\}|$

over all $\\{x_{i}\\}_{i=1}^{m}\subset X$ and
$\\{y_{i}\\}_{i=1}^{m}\subset\mathbb{R}$. Actually the maximum of
$K_{\\{x_{i}\\},\\{y_{i}\\}}(m)$ over all $\\{x_{i}\\}_{i=1}^{m}\subset X$ and
$\\{y_{i}\\}_{i=1}^{m}\subset\mathbb{R}$ is the growth function
$\mathcal{G}_{\mathcal{\widetilde{N}}}(m)$. In order to apply Lemma 20, we
partition the parameter space $\mathbb{R}^{\mathcal{M}}$ into several subsets
to ensure that in each subset $u(x_{i},a)-y_{i}$ is a polynomial with respcet
to $a$ without any breakpoints. In fact, our partition is exactly the same as
the partition in [2]. Denote the partition as $\\{P_{1},P_{2},\cdots,P_{N}\\}$
with some integer $N$ satisfying

$N\leq\prod_{i=1}^{\mathcal{D}-1}2\left(\frac{2emk_{i}(1+(i-1)2^{i-1})}{\mathcal{M}_{i}}\right)^{\mathcal{M}_{i}}$
(13)

where $k_{i}$ and $\mathcal{M}_{i}$ denotes the number of units at the $i$th
layer and the total number of parameters at the inputs to units in all the
layers up to layer $i$ of the neural network implementing functions in
$\mathcal{N}$, respectively. See [2] for the construction of the partition.
Obviously we have

$K_{\\{x_{i}\\},\\{y_{i}\\}}(m)\leq\sum_{i=1}^{N}|\\{(\operatorname{sign}(u(x_{1},a)-y_{1}),\cdots,\operatorname{sign}(u(x_{m},a)-y_{m})):a\in
P_{i}\\}|$ (14)

Note that $u(x_{i},a)-y_{i}$ is a polynomial with respect to $a$ with degree
the same as the degree of $u(x_{i},a)$, which is equal to
$1+(\mathcal{D}-1)2^{\mathcal{D}-1}$ as shown in [2]. Hence by Lemma 20, we
have

$\displaystyle|\\{(\operatorname{sign}(u(x_{1},a)-y_{1}),\cdots,\operatorname{sign}(u(x_{m},a)-y_{m})):a\in
P_{i}\\}|$ $\displaystyle\leq
2\left(\frac{2em(1+(\mathcal{D}-1)2^{\mathcal{D}-1})}{\mathcal{M}_{\mathcal{D}}}\right)^{\mathcal{M}_{\mathcal{D}}}.$
(15)

Combining $(\ref{pdimb1}),(\ref{pdimb2}),(\ref{pdimb3})$ yields

$K_{\\{x_{i}\\},\\{y_{i}\\}}(m)\leq\prod_{i=1}^{\mathcal{D}}2\left(\frac{2emk_{i}(1+(i-1)2^{i-1})}{\mathcal{M}_{i}}\right)^{\mathcal{M}_{i}}.$

We then have

$\mathcal{G}_{\mathcal{\widetilde{N}}}(m)\leq\prod_{i=1}^{\mathcal{D}}2\left(\frac{2emk_{i}(1+(i-1)2^{i-1})}{\mathcal{M}_{i}}\right)^{\mathcal{M}_{i}},$

since the maximum of $K_{\\{x_{i}\\},\\{y_{i}\\}}(m)$ over all
$\\{x_{i}\\}_{i=1}^{m}\subset X$ and $\\{y_{i}\\}_{i=1}^{m}\subset\mathbb{R}$
is the growth function $\mathcal{G}_{\mathcal{\widetilde{N}}}(m)$. Some
algebras as that of the proof of Theorem 6 in [2], we obtain

$\mathrm{Pdim}(\mathcal{N})\leq\mathcal{O}\left(\mathcal{D}^{2}\mathcal{W}^{2}\log\mathcal{U}+\mathcal{D}^{3}\mathcal{W}^{2}\right)=\mathcal{O}\left(\mathcal{D}^{2}\mathcal{W}^{2}\left(\mathcal{D}+\log\mathcal{W}\right)\right)$

where $\mathcal{U}$ refers to the number of units of the neural network
implementing functions in $\mathcal{N}$. ∎

###  3.3 Main results

With the above preparation we present the main results of this paper in the
scenario $\mathcal{E}_{opt}=0$.

######  Theorem 3.3.

Let $u^{*}$ is the solution of (1) with bounded $f,g,w$. $\widehat{u}_{\phi}$
is the minimizer of deep Ritz method defined in (3) with $n=N(M)$ random
samples. If we set the network parameters depth and width in the
$\mathrm{ReLU}^{2}$ network
$\mathcal{N}^{2}_{\mathcal{D},\mathcal{W},\mathcal{B}}$ as

$\mathcal{D}\leq\lceil\log_{2}d\rceil+3,\mathcal{W}\leq\mathcal{O}(4d\left\lceil
n^{\frac{1}{d+2+\nu}}-4\right\rceil^{d}).$

Let $\mathcal{B}$ be the constraint on the output of function value and
gradient norm of $u\in\mathcal{N}^{2}$ according to (9). Then,

$\mathbb{E}_{\mathbf{X},\mathbf{Y}}[\|\widehat{u}_{\phi}-u^{*}\|_{H^{1}(\Omega)}^{2}]\leq
C_{\mathcal{B},c_{1},c_{2},c_{3},d}{\mathcal{O}}(N^{-1/(d+2+\nu)}+M^{-1/(d+2+\nu)}),$

where $\nu>0$ but can be arbitrary small.

###### Proof.

In order to apply Lemma 12, we need to bound the term

$\begin{split}&\frac{1}{\sqrt{n}}\int_{\delta}^{\mathcal{B}}\sqrt{\log(2\mathcal{C}(\epsilon,\mathcal{N},n))}d\epsilon\\\
&\leq\frac{\mathcal{B}}{\sqrt{n}}+\frac{1}{\sqrt{n}}\int_{\delta}^{B}\sqrt{\log\left(\frac{en\mathcal{B}}{\epsilon\cdot\mathrm{Pdim}(\mathcal{N})}\right)^{\mathrm{Pdim}(\mathcal{N})}}d\epsilon\\\
&\leq\frac{\mathcal{B}}{\sqrt{n}}+\left(\frac{\mathrm{Pdim}(\mathcal{N})}{n}\right)^{1/2}\int_{\delta}^{B}\sqrt{\log\left(\frac{en\mathcal{B}}{\epsilon\cdot\mathrm{Pdim}(\mathcal{N})}\right)}d\epsilon,\end{split}$

where in the first inequality we use Proposition 19. Now we calculate the
integral. Set

$t=\sqrt{\log\left(\frac{en\mathcal{B}}{\epsilon\cdot\mathrm{Pdim}(\mathcal{N})}\right)}$

then $\epsilon=\frac{en\mathcal{B}}{\mathrm{Pdim}(\mathcal{N})}\cdot
e^{-t^{2}}$. Denote
$t_{1}=\sqrt{\log\left(\frac{en\mathcal{B}}{\mathcal{B}\cdot\mathrm{Pdim}(\mathcal{N})}\right)}$,
$t_{2}=\sqrt{\log\left(\frac{en\mathcal{B}}{\delta\cdot\mathrm{Pdim}(\mathcal{N})}\right)}$.
And

$\begin{split}&\int_{\delta}^{\mathcal{B}}\sqrt{\log\left(\frac{en\mathcal{B}}{\epsilon\cdot\mathrm{Pdim}(\mathcal{N})}\right)}d\epsilon\\\
&=\frac{2en\mathcal{B}}{\mathrm{Pdim}(\mathcal{N})}\int_{t_{1}}^{t_{2}}t^{2}e^{-t^{2}}dt\\\
&=\frac{2en\mathcal{B}}{\mathrm{Pdim}(\mathcal{N})}\int_{t_{1}}^{t_{2}}t\left(\frac{-e^{-t^{2}}}{2}\right)^{\prime}dt\\\
&=\frac{en\mathcal{B}}{\mathrm{Pdim}(\mathcal{N})}\left[t_{1}e^{-t_{1}^{2}}-t_{2}e^{-t_{2}^{2}}+\int_{t_{1}}^{t_{2}}e^{-t^{2}}dt\right]\\\
&\leq\frac{en\mathcal{B}}{\mathrm{Pdim}(\mathcal{N})}\left[t_{1}e^{-t_{1}^{2}}-t_{2}e^{-t_{2}^{2}}+(t_{2}-t_{1})e^{-t_{1}^{2}}\right]\\\
&\leq\frac{en\mathcal{B}}{\mathrm{Pdim}(\mathcal{N})}\cdot
t_{2}e^{-t_{1}^{2}}=\mathcal{B}\sqrt{\log\left(\frac{en\mathcal{B}}{\delta\cdot\mathrm{Pdim}(\mathcal{N})}\right)}.\end{split}$

Choosing
$\delta=\mathcal{B}\left(\frac{\mathrm{Pdim}(\mathcal{N})}{n}\right)^{1/2}\leq\mathcal{B}$,
by Lemma 12 and the above display, we get for both

$\mathcal{N}=\mathcal{N}^{2},\\\ \mathrm{and}\ \
\mathcal{N}=\mathcal{N}^{1,2}$

there holds

$\displaystyle\mathfrak{R}(\mathcal{N})$ $\displaystyle\leq
4\delta+\frac{12}{\sqrt{n}}\int_{\delta}^{\mathcal{B}}\sqrt{\log(2\mathcal{C}(\epsilon,\mathcal{N},n))}d\epsilon$
$\displaystyle\leq
4\delta+\frac{12\mathcal{B}}{\sqrt{n}}+12\mathcal{B}\left(\frac{\mathrm{Pdim}(\mathcal{N})}{n}\right)^{1/2}\sqrt{\log\left(\frac{en\mathcal{B}}{\delta\cdot\mathrm{Pdim}(\mathcal{N})}\right)}$
$\displaystyle\leq
28\sqrt{\frac{3}{2}}\mathcal{B}\left(\frac{\mathrm{Pdim}(\mathcal{N})}{n}\right)^{1/2}\sqrt{\log\left(\frac{en}{\mathrm{Pdim}(\mathcal{N})}\right)}.$
(16)

Then by Lemma 5, 8, 9, 12, 14 and equation (16), we have

$\displaystyle{\mathcal{E}_{sta}}=2\sup_{u\in\mathcal{N}^{2}}|\mathcal{L}(u)-\widehat{\mathcal{L}}(u)|$
$\displaystyle\leq
2(c_{3}+c_{3}+c_{3}^{2})\mathfrak{R}(\mathcal{N}^{2})+2\mathfrak{R}(\mathcal{N}^{1,2})$
$\displaystyle\leq
56\sqrt{\frac{3}{2}}(2c_{3}+c_{3}^{2})\mathcal{B}\left(\frac{\mathrm{Pdim}(\mathcal{N}^{2})}{n}\right)^{1/2}\sqrt{\log\left(\frac{en}{\mathrm{Pdim}(\mathcal{N}^{2})}\right)}$
$\displaystyle+56\sqrt{\frac{3}{2}}\mathcal{B}\left(\frac{\mathrm{Pdim}(\mathcal{N}^{1,2})}{n}\right)^{1/2}\sqrt{\log\left(\frac{en}{\mathrm{Pdim}(\mathcal{N}^{1,2})}\right)}.$

Since $\mathrm{Pdim}(\mathcal{N}^{2})\leq\mathrm{Pdim}(\mathcal{N}^{1,2})$,
plugging the upper bound of $\mathrm{Pdim}(\mathcal{N}^{1,2})$ derived in
Theorem 3.2 into the above display and using the relationship of depth and
width between $\mathcal{N}^{2}$ and $\mathcal{N}^{1,2}$ proved in Lemma 14, we
get

$\mathcal{E}_{sta}\leq
C_{\mathcal{B},c_{3}}\left[d(\mathcal{D}+3)(\mathcal{D}+2)\mathcal{W}\sqrt{\frac{\mathcal{D}+3+\log(d(\mathcal{D}+2)\mathcal{W})}{n}}\right]^{1-\nu},$
(17)

where $\nu>0$ can be arbitrarily small. Combing (17) with the approximation
error in Theorem 3.1 and taking
$\epsilon^{2}=C_{\mathcal{B},c_{1},c_{2},c_{3}}C\left(\frac{1}{n}\right)^{\frac{1}{d+2+\nu}}$,
we get that

$\displaystyle\mathbb{E}_{\mathbf{X},\mathbf{Y}}[\|\widehat{u}_{\phi}-u^{*}\|_{H^{1}(\Omega)}^{2}]$
$\displaystyle\leq\frac{2}{c_{1}\wedge
1}\left[C_{\mathcal{B},c_{3}}\left(d(\mathcal{D}+3)(\mathcal{D}+2)\mathcal{W}\sqrt{\frac{\mathcal{D}+3+\log(d(\mathcal{D}+2)\mathcal{W})}{n}}\right)^{1-\mu}+\frac{c_{3}+1}{2}\epsilon^{2}\right]$
$\displaystyle\leq\frac{2}{c_{1}\wedge
1}\left[C_{\mathcal{B},c_{3}}\left(4d^{2}(\lceil\log d\rceil+6)(\lceil\log
d\rceil+5)\left\lceil\frac{Cc_{2}}{\epsilon}-4\right\rceil^{d}\cdot\right.\right.$
$\displaystyle\quad\left.\left.\sqrt{\frac{\lceil\log
d\rceil+6+\log\left(d(\lceil\log d\rceil+5)\cdot
4d\left\lceil\frac{Cc_{2}}{\epsilon}-4\right\rceil^{d}\right)}{n}}\right)^{1-\mu}+\frac{c_{3}+1}{2}\epsilon^{2}\right]$
$\displaystyle\leq
C_{\mathcal{B},c_{1},c_{2},c_{3},d}{\mathcal{O}}\left(n^{-1/(d+2+\nu)}\right).$

∎

######  Remark 3.1.

Deep Ritz method is actually a kinds of deep nonparametric estimation method
where we estimate functions from random samples. The benefit is that we can
use DRM to handle PDEs in high dimension since only a small bach of samples is
needed during SGD training. In contrast, we have form the loading matrix and
vector explicitly in FEM. However, what we have to pay is that the convergence
rate as illustrate as follows. Let $N=N=\mathcal{O}(\frac{1}{h^{d}})$, where
$h$ is the size of the mesh in FEM. From Theorem 3.3, we get

$\mathbb{E}_{\mathbf{X},\mathbf{Y}}[\|\widehat{u}_{\phi}-u^{*}\|_{H^{1}(\Omega)}^{2}]=\mathcal{O}(h^{\frac{d}{d+2+\nu}}).$

In the case $d=2$, by Markov’s inequality and the above display, we get with
high probability,

$\|\widehat{u}_{\phi}-u^{*}\|_{H^{1}(\Omega)}\leq\mathcal{O}(h^{\frac{1}{4+\nu}}).$

Comparing the well known results of FEM where the convergence rate in $H^{1}$
norm is $\mathcal{O}(h)$, the rate proved here for the DRM method is far from
satisfactory.

In the literature of nonparametric regression, where functions living in
certain function classes are estimated from $n$ paired random samples, the
best convergence rate in $H^{1}$ norm for estimating functions in $H^{2}$ is
$\mathcal{O}(n^{-\frac{1}{4+d}})$ [41]. Obviously, the nonparametric learning
task in DRM is not easier than nonparametric regressions. What we proved here
is $\mathcal{O}(n^{-\frac{1}{4+2d+\nu}})$ in $H^{1}$ norm with $\nu$ can be
arbitrary small. We conjecture that the best convergence rate in $H^{1}$ norm
of DRM with both $n$ training samples in both the domain and boundary is also
$\mathcal{O}(n^{-\frac{1}{4+d}})$. Upon this conjecture, we get the optimal
convergence rate of DRM for $d=2$ is
$\|\widehat{u}_{\phi}-u^{*}\|_{H^{1}(\Omega)}\leq\mathcal{O}(h^{\frac{1}{3}})$
with high probability.

The convergence rate of DRM proved here suffers the curse of dimensionality.
One possible direction to improve the convergence rate and reduce the curse of
dimensionality is considering solutions of PDE (1) with higher regularity. For
example, if we assume $f\in H^{s}(\Omega),s\geq 1$, deep Ritz method using
deep neural network with $\mathrm{ReLU}^{s+2}$ activation functions will
reduce the curse since the higher regularity assumption can improve both the
approximation and statistical error. We leave the detail of this idea in a
following up work.

##  4 Conclusion and extension

In this paper, we provide a rigorous numerical analysis on deep Ritz method
(DRM) [48] for second order elliptic equations with Neumann boundary
conditions. We establish the first nonasymptotic convergence rate in $H^{1}$
norm on DRM for general deep networks with $\mathrm{ReLU}^{2}$ activation
functions. In addition to provide theoretical justification of DRM, our study
also provides guidance on how to set the hyper-parameter of depth and width to
achieve the desired convergence rate in terms of number of training samples.
Technically, we derive bounds on the approximation error of deep
$\mathrm{ReLU}^{2}$ network in $H^{1}$ norm and on the Rademacher complexity
of the non-Lipschitz composition of gradient norm and $\mathrm{ReLU}^{2}$
network.

There are several directions for our future exploration. First, it is easy to
extend the current analysis for general second order elliptic equations with
variational form under Dirichlet or Robin boundary conditions. Second, the
approximation and statistical error bounds deriving here can be used for
studying the nonasymptotic convergence rate for residual based method (PINNS).
Studying deep DGM by combing current analysis with the tools for analyzing GAN
[14] is also of immense interest.

## Acknowledgements

Y. Jiao was supported in part by the National Science Foundation of China
under Grant 11871474 and by the research fund of KLATASDSMOE. X. Lu is
partially supported by the National Science Foundation of China (No.
11871385), the National Key Research and Development Program of China
(No.2018YFC1314600) and the Natural Science Foundation of Hubei Province (No.
2019CFA007). J. Yang is supported by National Science Foundation of China (No.
12071362 and 11671312), the National Key Research and Development Program of
China (No. 2020YFA0714200), the Natural Science Foundation of Hubei Province
(No. 2019CFA007). The numerical simulations in this work have been done on the
supercomputing system in the Supercomputing Center of Wuhan University.

## References

  * [1] Martin Anthony and Peter L Bartlett, Neural network learning: Theoretical foundations, cambridge university press, 2009.
  * [2] Peter L Bartlett, Nick Harvey, Christopher Liaw, and Abbas Mehrabian, Nearly-tight vc-dimension and pseudodimension bounds for piecewise linear neural networks., J. Mach. Learn. Res., 20 (2019), pp. 1–17.
  * [3] Christian Beck, Weinan E, and Arnulf Jentzen, Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations, Journal of Nonlinear Science, 29 (2019), pp. 1563–1619.
  * [4] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart, Concentration inequalities: A nonasymptotic theory of independence, Oxford university press, 2013.
  * [5] Susanne Brenner and Ridgway Scott, The mathematical theory of finite element methods, vol. 15, Springer Science & Business Media, 2007.
  * [6] Fan Chen, Jianguo Huang, Chunmei Wang, and Haizhao Yang, Friedrichs learning: Weak solutions of partial differential equations via deep learning, arXiv preprint arXiv:2012.08023, (2020).
  * [7] Philippe G Ciarlet, The finite element method for elliptic problems, SIAM, 2002.
  * [8] Carl De Boor and Carl De Boor, A practical guide to splines, vol. 27, springer-verlag New York, 1978.
  * [9] A. Debussche, M. Fuhrman, and G. Tessitore, Optimal control of a stochastic heat equation with boundary-noise and boundary-control, ESAIM: Control, Optimisation and Calculus of Variations, 13 (2007), pp. 178–205.
  * [10] Niccol Discacciati, Jan S. Hesthaven, and Deep Ray, Controlling oscillations in high-order discontinuous galerkin schemes using artificial viscosity tuned by neural networks, Journal of Computational Physics, 409 (2020), p. 109304.
  * [11] R.M Dudley, The sizes of compact subsets of hilbert space and continuity of gaussian processes, Journal of Functional Analysis, 1 (1967), pp. 290–330.
  * [12] Dinh Dũng and Van Kien Nguyen, Sparse-grid sampling recovery and deep relu neural networks in high-dimensional approximation, arXiv preprint arXiv:2007.08729, (2020).
  * [13] Lawrence C Evans, Partial differential equations, Graduate studies in mathematics, 19 (1998).
  * [14] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio, Generative adversarial nets, in Advances in Neural Information Processing Systems 27, Curran Associates, Inc., 2014, pp. 2672–2680.
  * [15] Ingo Gühring, Gitta Kutyniok, and Philipp Petersen, Error bounds for approximations with deep relu neural networks in w s, p norms, Analysis and Applications, 18 (2020), pp. 803–859.
  * [16] T. J. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Courier Corporation, 2012.
  * [17] Ameya D. Jagtap and George Em Karniadakis, Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations, Communications in Computational Physics, 28 (2020), pp. 2002–2041.
  * [18] Ameya D Jagtap, Ehsan Kharazmi, and George Em Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Computer Methods in Applied Mechanics and Engineering, 365 (2020), p. 113028.
  * [19] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE transactions on neural networks, 9 (1998), pp. 987–1000.
  * [20] Michel Ledoux and Michel Talagrand, Probability in Banach Spaces: isoperimetry and processes, Springer Science & Business Media, 2013.
  * [21] Hyuk Lee and In Seok Kang, Neural algorithm for solving differential equations, Journal of Computational Physics, 91 (1990), pp. 110–131.
  * [22] Randall J LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, 2007.
  * [23] Bo Li, Shanshan Tang, and Haijun Yu, Better approximations of high dimensional smooth functions by deep neural networks with rectified power units, arXiv preprint arXiv:1903.05858, (2019).
  * [24] Jianfeng Lu, Yulong Lu, and Min Wang, A priori generalization analysis of the deep ritz method for solving high dimensional elliptic equations, arXiv preprint arXiv:2101.01708, (2021).
  * [25] Lu Lu, Xuhui Meng, Zhiping Mao, and George E. Karniadakis, Deepxde: A deep learning library for solving differential equations, CoRR, abs/1907.04502 (2019).
  * [26] Tao Luo and Haizhao Yang, Two-layer neural networks for partial differential equations: Optimization and generalization theory, arXiv preprint arXiv:2006.15733, (2020).
  * [27] Siddhartha Mishra and Roberto Molinaro, Estimates on the generalization error of physics informed neural networks (pinns) for approximating pdes, arXiv preprint arXiv:2006.16144, (2020).
  * [28] Hadrien Montanelli and Qiang Du, New error bounds for deep relu networks using sparse grids, SIAM Journal on Mathematics of Data Science, 1 (2019), pp. 78–92.
  * [29] G. Pang, M. D’Elia, M. Parks, and G.E. Karniadakis, npinns: Nonlocal physics-informed neural networks for a parametrized nonlocal universal laplacian operator. algorithms and applications, Journal of Computational Physics, 422 (2020), p. 109760.
  * [30] Guofei Pang, Lu Lu, and George Em Karniadakis, fpinns: Fractional physics-informed neural networks, SIAM Journal on Scientific Computing, 41 (2019), pp. A2603–A2626.
  * [31] Maziar Raissi, Paris Perdikaris, and George E Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics, 378 (2019), pp. 686–707.
  * [32] Deep Ray and Jan S. Hesthaven, An artificial neural network as a troubled-cell indicator, Journal of Computational Physics, 367 (2018), pp. 166–191.
  * [33] Johannes Schmidt-Hieber, Deep relu network approximation of functions on a manifold, arXiv preprint arXiv:1908.00695, (2019).
  * [34] Martin H Schultz, Approximation theory of multivariate spline functions in sobolev spaces, SIAM Journal on Numerical Analysis, 6 (1969), pp. 570–582.
  * [35] Larry Schumaker, Spline functions: basic theory, Cambridge University Press, 2007.
  * [36] Zuowei Shen, Haizhao Yang, and Shijun Zhang, Deep network approximation characterized by number of neurons, arXiv preprint arXiv:1906.05497, (2019).
  * [37] Yeonjong Shin, Jerome Darbon, and George Em Karniadakis, On the convergence and generalization of physics informed neural networks, arXiv preprint arXiv:2004.01806, (2020).
  * [38] Yeonjong Shin, Zhongqiang Zhang, and George Em Karniadakis, Error estimates of residual minimization using neural networks for linear pdes, arXiv preprint arXiv:2010.08019, (2020).
  * [39] Jonathan W Siegel and Jinchao Xu, Approximation rates for neural networks with general activation functions, Neural Networks, 128 (2020), pp. 313–321.
  * [40] Justin A. Sirignano and K. Spiliopoulos, Dgm: A deep learning algorithm for solving partial differential equations, Journal of Computational Physics, 375 (2018), pp. 1339–1364.
  * [41] Charles J Stone, Optimal global rates of convergence for nonparametric regression, The annals of statistics, (1982), pp. 1040–1053.
  * [42] Zheng Sun, Shuyi Wang, Lo-Bin Chang, Yulong Xing, and Dongbin Xiu, Convolution neural network shock detector for numerical solution of conservation laws, Communications in Computational Physics, 28 (2020), pp. 2075–2108.
  * [43] Taiji Suzuki, Adaptivity of deep relu network for learning in besov and mixed smooth besov spaces: optimal rate and curse of dimensionality, arXiv preprint arXiv:1810.08033, (2018).
  * [44] J.W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, vol. 22, Springer Science & Business Media, 2013.
  * [45] Aad W Van Der Vaart and Jon A Wellner, Weak convergence, in Weak convergence and empirical processes, Springer, 1996.
  * [46] Vladimir N Vapnik and A Ya Chervonenkis, On the uniform convergence of relative frequencies of events to their probabilities, in Measures of complexity, Springer, 2015, pp. 11–30.
  * [47] Yufei Wang, Ziju Shen, Zichao Long, and Bin Dong, Learning to discretize: Solving 1d scalar conservation laws via deep reinforcement learning, Communications in Computational Physics, 28 (2020), pp. 2158–2179.
  * [48] E. Weinan and Ting Yu, The deep ritz method: A deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics, 6 (2017), pp. 1–12.
  * [49] Jinchao Xu, The finite neuron method and convergence analysis, arXiv preprint arXiv:2010.01458, (2020).
  * [50] Dmitry Yarotsky, Error bounds for approximations with deep relu networks, Neural Networks, 94 (2017), pp. 103–114.
  * [51] Yaohua Zang, Gang Bao, Xiaojing Ye, and Haomin Zhou, Weak adversarial networks for high-dimensional partial differential equations, Journal of Computational Physics, 411 (2020), p. 109409.