# Convergence analysis of unsupervised Legendre-Galerkin neural networks for
linear second-order elliptic PDEs

Seungchan Ko, Seok-Bae Yun, and Youngjoon Hong Department of Mathematics,
Sungkyunkwan University, Suwon, Republic of Korea. Email:
<EMAIL_ADDRESS>of Mathematics, Sungkyunkwan University, Suwon,
Republic of Korea. Email<EMAIL_ADDRESS>of Mathematics,
Sungkyunkwan University, Suwon, Republic of Korea. Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we perform the convergence analysis of unsupervised
Legendre–Galerkin neural networks (ULGNet), a deep-learning-based numerical
method for solving partial differential equations (PDEs). Unlike existing deep
learning-based numerical methods for PDEs, the ULGNet expresses the solution
as a spectral expansion with respect to the Legendre basis and predicts the
coefficients with deep neural networks by solving a variational residual
minimization problem. Since the corresponding loss function is equivalent to
the residual induced by the linear algebraic system depending on the choice of
basis functions, we prove that the minimizer of the discrete loss function
converges to the weak solution of the PDEs. Numerical evidence will also be
provided to support the theoretical result. Key technical tools include the
variant of the universal approximation theorem for bounded neural networks,
the analysis of the stiffness and mass matrices, and the uniform law of large
numbers in terms of the Rademacher complexity.

Keywords: Deep neural network, unsupervised learning, elliptic partial
differential equations, Legendre-Galerkin approximation, spectral element
method, Rademacher complexity

AMS Classification: 68T07, 65N35, 65K10, 65N12, 35J25

##  1 Introduction

In recent years, modern machine learning techniques using deep neural networks
have had tremendous success in various areas from science and engineering:
computer vision [28], natural language processes [31], time series analysis
[46], cognitive science [30], and so on. The notable advance of deep learning
algorithms has offered interesting possibilities for computational
innovations. In particular, based on the rich expressiveness of neural
networks, it has recently begun to gain more attention to solving partial
differential equations (PDEs) using a deep learning algorithm.

PDEs are fundamental mathematical models in studying a wide range of phenomena
naturally arising in science and engineering. Numerical approximation of PDEs
has a long and rich history, and the mathematical theory of the numerical
methods is well-established; see, for example, [7, 10] for the finite element
methods and [50] where the spectral method was studied. On the other hand, in
order to tackle the problems of high-dimensional PDEs, various deep learning-
based techniques solving PDEs were developed such as Deep Ritz Method (DRM)
[15], Deep Galerkin Methods (DGM) [56] and Physics-Informed Neural Network
(PINN) [47, 25]. By utilizing the approximation properties of deep neural
networks and statistical learning theory to numerical PDEs, these new
approaches opened a new area of research which is recently called scientific
machine learning. Since neural networks have a great capability of
approximating nonlinear functions, approximated solutions with parameters
which are updated through machine learning-based optimization techniques have
shown notable performance in solving PDEs numerically. Those methods
introduced certain types of residuals to define a loss function induced from
the given PDEs, which is minimized by an optimization algorithm over the
randomly-selected training points on the domain. In particular, the PINN
exploits randomly-chosen collocation points as training data in the space-time
domain, and hence, the PINN is available for time-dependent high-dimensional
PDEs, on the computational domains of complex geometry [36, 34, 53, 39, 24,
58].

However, these approaches have some drawbacks. For example, only a single
instance from a fixed input data, including initial conditions, boundary
conditions, and external force terms is available at a time. This means that
if the input data is changed, the whole training process needs to be repeated.
Hence, real-time predictions for varying input data are not possible for the
aforementioned methods. Some techniques have been developed to overcome this
problem; data-supervised or data-driven methods (DDM) [9, 2, 61], the Fourier
neural operator (FNO) [32], and the DeepONet (DON) [35]. However, training
these models requires a database, which should be provided in advance by
solving the equations analytically or numerically. To generate a reliable
training dataset, a large number of numerical solutions need to be provided by
massive numerical computations, which is inefficient and often takes a very
long time. Furthermore, using the class of neural networks as the solution
space makes it difficult to impose boundary values. In general, for a function
represented by neural networks, it is not straightforward how to incorporate
certain values on the boundary of the given domain. One common strategy to
deal with this issue is to add penalty terms to loss functions. However, it is
known that the boundary penalty term might cause worse convergence of the
approximate solution and some technical difficulties may arise.

To address the aforementioned issues, various attempts have been made in
recent years. In particular, in [8] the unsupervised Legendre–Galerkin neural
network (ULGNet) was proposed to solve various PDEs without a training
dataset. This approach is based on the classical numerical method, the
Legendre–Galerkin spectral element method for the boundary value problem [19,
50]. By choosing the Legendre polynomials as basis functions, the
corresponding polynomial coefficients are generated via deep neural networks
whose input is data of PDEs such as external forcing terms or initial
conditions. Then the predicted solution of the ULGNet can be obtained by a
linear combination of the coefficients and the Legendre basis functions. For
the training, the loss function is defined as a residual of summation of weak
formulations which take their test function as Legendre basis functions. Then
the random points are sampled from the parametric domain which determines the
source term (input of the neural network) and the neural network is trained in
the direction of minimizing the loss function. Consequently, the ULGNet can
overcome the shortcomings of various models in scientific machine learning.
Since it is based on unsupervised learning process, the expensive training
database is not required while the algorithm is able to predict multiple
instances of PDEs. In addition, the predicted solution can have the exact
boundary values as given in accordance with the choice of spatial basis
functions. Using a weak form as the loss function instead of the strong form
used in other literature (e.g. PINN [47] or DGM [56]) is advantageous in some
aspects, since we can avoid high-order numerical derivatives that cause large
numerical errors, which may be able to enhance the performance of neural
networks. Moreover, using the Legendre basis functions can also improve the
accuracy of the approximation because the numerical integration that appeared
in the ULGNet algorithm can be computed by the Gauss–Legendre quadrature rule
which is highly accurate. From the numerical experiments performed in [8], it
was shown that the ULGNet can compute the approximate solutions with high
accuracy in various situations. Especially, this technique performs very well
for the singularly-perturbed and boundary layer problems, which is mainly due
to the fact that we have freedom in manipulating the basis functions, and
hence we are able to design an enriched scheme in this framework. The more
detailed discussion of ULGNet will be studied in Section 2.

Despite a wide range of successful applications for neural network-based PDE
solvers, the rigorous analysis of such methods is not well developed and only
a few theoretical results are currently known.

Some papers provided evidence of convergence for the neural network-based
method [56, 6]. Note, however, that these papers only proved that there exists
a sequence of neural networks which converges to the true solution, rather
than showing the actual convergence of approximate solutions. Some other
previous research papers on the analysis of neural network-based algorithms
solving PDEs investigated on the representation aspect. In other words, they
focused on whether a solution of the given PDE is well-approximated by neural
networks with quantitative control of numerical errors; see, for example, [20,
23, 42, 37]. On the other hand, if we fix an ansatz space, the generalization
error can be estimated by analyzing complexities using covering numbers [5,
37]. More recently, the analysis for the DRM in a high-dimensional space has
been established; see [37] where the class of two-layer neural networks in the
spectral Barron framework was studied and a more general setting was
considered in [13]. Moreover, some recent papers investigated the
generalization error analysis of the PINN [51, 40, 38] based on the residual
minimization problem [29, 52]. More precisely, [51] focused on the consistency
of the loss function so that the approximate solution converges to the exact
solution as the number of training samples increases with the assumption of
vanishing optimization error. The authors in [40] conducted an a-posteriori
error analysis and showed that the generalization error can be bounded by the
training error and the quadrature error under some assumptions of stability of
the given PDEs. The work in [38] proved a priori generalization error
estimates for a class of second-order linear PDEs by assuming that the true
solutions of PDEs are contained in a Barron-type space, and [52] obtained both
a priori and a posterior error estimates for residual minimization problems in
Sobolev spaces setting. More results on the theoretical analysis of neural
network-based methods of PDEs can be found in [6, 41, 12, 21].

In this perspective, the central aim of the present paper is to perform the
theoretical analysis on the convergence of ULGNet, which supports the
experimental results in [8]. Note that all the aforementioned convergence
results were about the neural network for the single instance prediction. To
the best of our knowledge, this is the first paper which studies convergence
analysis for the unsupervised neural network for the prediction of multiple
instances. More precisely, we consider the following second-order, self-
adjoint linear elliptic PDE: for some $\nu\geq 0$,

$-\Delta u+\nu u=f\quad{\rm{in}}\,\,I^{d}=[-1,1]^{d},$ (1.1)

with either Dirichlet or Neumann boundary conditions. As mentioned earlier, in
the ULGNet setting, the external forcing term $f$ is parametrized by random
samples $\omega$ from a parametric domain which we shall denote by $\Omega$,
and hence the solution $u$ also depends on the parametric variable
$\omega\in\Omega$. Therefore, a proper distance to be measured is the Bochner-
type norm where the distance between the true solution and the approximate
solution is integrated over both $I^{d}$ and $\Omega$. The mathematical
setting and detailed problem formulation will be discussed in a later section.

In general, for the neural network-based methods, the theoretical analysis
heavily depends on the form of given PDEs and the type of boundary conditions.
For example, the authors of [37] established the convergence of the DRM of the
equation (1.1) for both $\nu=0$ and $\nu\neq 0$, and the analysis of the case
$\nu\neq 0$ was much more complicated than the one for the case $\nu=0$. Also,
by referring [37, 13], we note that the analysis with the Dirichlet boundary
condition is quite different from the analysis for the Neumann condition. As
it will be made clear in the later analysis, we will develop a unified
approach that can cover various settings with similar analysis, such as the
analysis for general $\nu\geq 0$ and either Dirichlet or Neumann boundary
conditions.

The rest of the present paper is organized as follows. In Section 2 we will
introduce some preliminaries that will be used throughout the paper, including
notations, a detailed description of the ULGNet, our ansatz class, and the
problem setting under consideration. In Section 3 we shall perform the
convergence analysis of the residual minimization problem under the abstract
setting and this analysis will be applied to concrete PDE problems in Section
4. In Section 5, we provide some numerical experiments which confirm the
theoretical result and a concluding remark will be discussed in 6.

##  2 Preliminaries and mathematical setting

In this section, we introduce some preliminaries related to the numerical
approximation of the given problem (2.1), which will be used throughout the
paper. For $1\leq p\leq\infty$ and $k\in\mathbb{N}$, we denote by
$L^{p}(\Omega)$ the Lebesgue space and by $W^{k,p}(\Omega)$ the standard
Sobolev space, with the abbreviation $H^{k}(\Omega)=W^{k,2}(\Omega)$, and we
write as $H^{1}_{0}(\Omega)$ the subspace of $H^{1}(\Omega)$ consisting of the
functions with zero trace. For a real Banach space $X$, we shall denote by
$C(\Omega;X)$ the subspace of the Bochner space $L^{\infty}(\Omega;X)$, which
consists of the continuous functions from $\Omega$ to $X$. Moreover,
$(\Omega,\mathcal{G},\mathbb{P})$ signifies a probability space, where
$\Omega$ is the sample space consisting of all possible outcomes,
$\mathcal{G}\subset 2^{\Omega}$ is the $\sigma$-algebra, and $\mathbb{P}$ is a
probability measure.

For the sake of simplicity, we write $\|\cdot\|_{p}=\|\cdot\|_{L^{p}(\Omega)}$
and $\|\cdot\|_{k,p}=\|\cdot\|_{W^{k,p}(\Omega)}$, and furthermore, $|\cdot|$
denotes the Euclidean norm. Throughout the paper, for two vectors
$\boldsymbol{a}$ and $\boldsymbol{b}$, $\boldsymbol{a}\cdot\boldsymbol{b}$
means their scalar product. Also, for any Lebesgue measurable set
$Q\subset\mathbb{R}^{d}$, $|Q|$ stands for the standard Lebesgue measure of
$Q$, and $C$ denotes a generic positive constant, which may differ at each
appearance. Finally, for two different quantities $A$ and $B$, the notation
$A\lesssim B$ means that there exists a positive constant $C>0$ such that
$A\leq CB$.

In order to simplify the presentation, we shall restrict ourselves to the
following one-dimensional problem: for some $\nu\geq 0$,

$-u^{\prime\prime}+\nu u=f\quad{\rm{in}}\,\,I=(-1,1),$ (2.1)

with either the homogeneous Dirichlet boundary condition $u(-1)=u(1)=0$ or the
Neumann boundary condition $u^{\prime}(-1)=u^{\prime}(1)=0$. As will be
discussed in Section 4.2, it can be extended to the multi-dimensional problem
with general boundary conditions.

###  2.1 Legendre–Galerkin neural network

In this section, we shall introduce one of the spectral element methods for
solving a two-point boundary value problem, the so-called Legendre–Galerkin
method. Based on this numerical method, we describe the unsupervised
Legendre–Galerkin neural networks (ULGNet) proposed in [8], which we aim to
analyze in the present paper.

To begin, we discuss some basic properties of Legendre polynomials that play a
role in basis functions for our numerical approximation. Legendre polynomials,
$L_{n}(x)$, are mutually orthogonal polynomials with respect to $L^{2}$-inner
product, defined on the interval $[-1,1]$. Here we collect some relevant
results for the Legendre polynomials (see, for example, [50]):

  * •

Three-terms recurrence relation:

$(n+1)L_{n+1}(x)=(2n+1)xL_{n}(x)-nL_{n-1}(x),\quad{\rm{with}}\,\,L_{0}(x)=1\,\,{\rm{and}}\,\,L_{1}(x)=x.$

  * •

The $n$-th Legendre polynomial can be written as the expansion

$L_{n}(x)=\frac{1}{2^{n}}\sum^{[n/2]}_{k=0}(-1)^{k}\frac{(2n-2k)!}{2^{n}k!(n-k)!(n-2k)!}x^{n-2k}.$

  * •

Sturm–Liouville problem:

$\left((1-x^{2})L^{\prime}_{n}(x)\right)^{\prime}+n(n+1)L_{n}(x)=0.$

  * •

Rodriques’ formula:

$L_{n}(x)=\frac{1}{2^{n}n!}\frac{{\rm{d}}^{n}}{{\rm{d}}x^{n}}\left[\left(x^{2}-1\right)^{n}\right].$

To make the stiffness and mass matrices to be sparse in the numerical
approximation, and enable us to impose exact boundary conditions, we adopt the
following compact combination of Legendre Polynomials as basis functions:

$\phi_{k}(x)=L_{k}(x)+a_{k}L_{k+1}(x)+b_{k}L_{k+2}(x),$ (2.2)

where the parameters $\\{a_{k},b_{k}\\}$ are chosen to satisfy the boundary
condition of the given problem. In fact, it is straightforward to show that
for all $k\in\mathbb{N}$, the set $\\{a_{k},b_{k}\\}$ can be uniquely
determined so that the basis function (2.2) satisfies the given boundary
condition either in Dirichlet sense or in Neumann sense.

Let us discuss the Legendre–Galerkin spectral method. The variational form of
the equation (2.1) can be written as

$\int_{I}\left(u^{\prime}(x)v^{\prime}(x)+\nu
u(x)v(x)\right)\,\mathrm{d}x=\int_{I}f(x)v(x)\,\mathrm{d}x\quad\forall v\in
H,$ (2.3)

for a suitable Hilbert space $H$. The existence and uniqueness of the solution
(2.3) follow by the Lax–Milgram theorem [16, 17]. The idea of the
Legendre–Galerkin method is to approximate the solution $u\in H$ by

$u_{N}(x)=\sum^{N-1}_{k=1}\alpha_{k}\phi_{k}(x),$ (2.4)

where $N$ is the number of nodal points on $I$ and $\phi_{k}$ is defined in
(2.2). The coefficients $\\{\alpha_{k}\\}_{k=1}^{N-1}$ can be determined by
solving the discrete approximation of (2.3), namely,

$\int_{I}\left(u^{\prime}_{N}(x)v^{\prime}_{N}(x)+\nu
u_{N}(x)v_{N}(x)\right)\,\mathrm{d}x=\int_{I}f(x)v_{N}(x)\,\mathrm{d}x\quad\forall
v_{N}\in H_{N},$ (2.5)

where $H_{N}={\rm{span}}\,\\{\phi_{1},\phi_{1},\cdots,\phi_{N-1}\\}$. More
precisely, if we write $S=(S_{ij})$,
$M=(M_{ij})\in\mathbb{R}^{(N-1)\times(N-1)}$, $F=(F_{j})\in\mathbb{R}^{N-1}$
with

$S_{ij}=\int_{I}\phi^{\prime}_{i}(x)\phi^{\prime}_{j}(x)\,\mathrm{d}x,\quad
M_{ij}=\int_{I}\phi_{i}(x)\phi_{j}(x)\,\mathrm{d}x,\quad
F_{j}=\int_{I}f(x)\phi_{j}(x)\,\mathrm{d}x,$ (2.6)

(2.5) is equivalent to the linear algebraic system

$(S+\nu M)\alpha=F,$ (2.7)

for $\alpha=(\alpha_{j})\in\mathbb{R}^{N-1}$, and (2.5) can be exactly solved
by

$\alpha=(S+\nu M)^{-1}F.$ (2.8)

Figure 1: Schematic description of the unsupervised Legendre–Galerkin neural
network (ULGNet).

The Unsupervised Legendre–Galerkin neural network, proposed in [8], is based
on the above Legendre–Galerkin spectral element methods. For each given $f$,
instead of computing the coefficients by using (2.7)-(2.8), we approximate the
coefficients $\alpha$ by a deep neural network. The input of the neural
network is external forcing terms, parametrized by the random parameter
$\omega$ contained in the (possibly high-dimensional) probability space
$(\Omega,\mathcal{G},\mathbb{P})$, where $\Omega$ is assumed to be compact.
One typical example is the random external force given by

$f(x;\omega)=\omega_{1}\sin(2\pi\omega_{2}x)+\omega_{3}\cos(2\pi\omega_{4}x),$

where $\omega=(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$ is i.i.d.
uniformly distributed with $\omega_{j}\in[a,b]$ for $a,b\in\mathbb{R}^{+}$ and
$1\leq j\leq 4$. Throughout the whole paper, we will interpret $f(x;\omega)$
as a bivariate function defined on $I\times\Omega$, and assume that

$f(x;\omega)\in C(\Omega;L^{1}(I)).$ (2.9)

This input feature passes through the neural network and generates the
coefficients $\\{\widehat{\alpha}_{k}\\}$ as an output. Then we reconstruct
the solution by

$\widehat{u}(x;\omega)=\sum_{k=1}^{N-1}\widehat{\alpha}_{k}(\omega)\phi_{k}(x).$
(2.10)

In order to train the neural network, we utilize the variational formulation
(2.5) for the residual minimization and define the population loss function by

$\mathcal{L}(\alpha)=\mathbb{E}_{\omega\sim\mathbb{P}_{\Omega}}\left[\sum^{N-1}_{k=1}\bigg{|}\int_{I}\left(\widehat{u}^{\prime}(x;\omega)\phi^{\prime}_{k}(x)+\nu\widehat{u}(x;\omega)\phi_{k}(x)\right)\,\mathrm{d}x-\int_{I}f(x;\omega)\phi_{k}(x)\,\mathrm{d}x\bigg{|}^{2}\right].$
(2.11)

In practice, when we compute the approximate solution, we minimize the
empirical loss function which is the Monte–Carlo integration of (2.11):

$\mathcal{L}^{M}(\alpha)=\frac{|\Omega|}{M}\sum^{M}_{m=1}\sum^{N-1}_{k=1}\bigg{|}\int_{I}\left(\widehat{u}^{\prime}(x;\omega_{m})\phi^{\prime}_{k}(x)+\nu\widehat{u}(x;\omega_{m})\phi_{k}(x)\right)\,\mathrm{d}x-\int_{I}f(x;\omega_{m})\phi_{k}(x)\,\mathrm{d}x\bigg{|}^{2},$
(2.12)

where $M\in\mathbb{N}$ is a number of training data and
$\\{\omega_{m}\\}_{m=1}^{M}$ is an i.i.d. sequence of random variables
distributed according to $\mathbb{P}_{\Omega}$. For each epoch, once the
parameters of the neural network are trained in the direction of minimizing
the empirical loss function $\mathcal{L}^{M}$, the given external force
functions pass through this updated neural network to produce more refined
coefficients, and this process is repeated until the sufficiently small loss
is achieved. Thanks to the intrinsic structure described above, ULGNet can be
trained without training data and incorporate the exact boundary values to the
approximate solutions. A schematic diagram of the ULGNet algorithm is depicted
in Figure 1, and more details with some numerical experiments can be found in
[8].

###  2.2 Feed-forward neural networks

In this section, we define a class of neural networks, which will serve as an
ansatz space throughout the whole paper. For a positive integer
$L\in\mathbb{N}$, we define an $L$-layer feed-forward neural network by a
function $f^{L}(x):\mathbb{R}^{n_{0}}\rightarrow\mathbb{R}^{n_{L}}$, which is
defined recursively by

$f^{\ell}(x)=W^{\ell}\sigma(f^{\ell-1}(x))+b^{\ell}\,\,\,{\rm{for}}\,\,2\leq\ell\leq
L\quad{\rm{and}}\quad f^{1}(x)=W^{1}x+b^{1}.$ (2.13)

Here $W^{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}$ is the weight matrix
and $b^{\ell}\in\mathbb{R}^{n_{\ell}}$ is the bias vector for the $\ell$-th
layer. Also, $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ denotes the activation
function, and for a vector $x=(x_{1},\cdots,x_{k})$, $\sigma(x)$ means the
vector $(\sigma(x_{1}),\cdots,\sigma(x_{k}))$. We shall write the architecture
of a network by the vector ${\vec{\boldsymbol{n}}}=(n_{0},\cdots,n_{L})$, the
family of network parameters by
$\theta:=\theta_{\vec{\boldsymbol{n}}}=\\{(W^{1},b^{1}),\cdots,(W^{L},b^{L})\\}$,
and its realization as a neural network by $\mathcal{R}[\theta](x)$. For a
given network architecture $\vec{\boldsymbol{n}}$, the collection of all
possible parameters is denoted by

$\Theta_{\vec{\boldsymbol{n}}}=\left\\{\\{(W^{\ell},b^{\ell})\\}^{L}_{\ell=1}:W^{\ell}\in\mathbb{R}^{n_{\ell}\times
n_{\ell-1}},\,b^{\ell}\in\mathbb{R}^{n_{\ell}}\right\\}.$ (2.14)

For the index for neural networks, we will follow the notations from [52, 44].
For two different neural networks $\theta_{i}$, $i=1,2$, whose architectures
are denoted by
$\vec{\boldsymbol{n}}_{i}=\left(n^{(i)}_{0},\cdots,n^{(i)}_{L_{i}}\right)$, we
write $\vec{\boldsymbol{n}}_{1}\subset\vec{\boldsymbol{n}}_{2}$ provided that
for all $\theta_{1}\in\Theta_{\vec{\boldsymbol{n}}_{1}}$, there exists
$\theta_{2}\in\Theta_{\vec{\boldsymbol{n}}_{2}}$ satisfying
$\mathcal{R}[\theta_{1}](x)=\mathcal{R}[\theta_{2}](x)$ for any
$x\in\mathbb{R}^{n_{0}}$.

Now let $\\{\vec{\boldsymbol{n}}_{n}\\}_{n\geq 1}$ be a sequence of network
architectures with the inclusion
$\vec{\boldsymbol{n}}_{n}\subset\vec{\boldsymbol{n}}_{n+1}$ for any
$n\in\mathbb{N}$. We shall define the corresponding sequence of neural network
families by

$\mathcal{N}^{*}_{n}=\\{\mathcal{R}[\theta]:\theta\in\Theta_{\vec{\boldsymbol{n}}_{n}}\\}.$
(2.15)

By the definition, it is easy to verify that
$\mathcal{N}^{*}_{n}\subset\mathcal{N}^{*}_{n+1}$ for all $n\in\mathbb{N}$. We
then state the following theorem on the universal approximation property.

######  Theorem 2.1.

Let $K$ be a compact subset of $\mathbb{R}^{d}$ and $g\in
C(K,\mathbb{R}^{D})$. Then we have

$\lim_{n\rightarrow\infty}\inf_{\hat{g}\in\mathcal{N}^{*}_{n}}\|\hat{g}-g\|_{C(K)}=0.$
(2.16)

It can be proved that this statement holds under various scenarios. For
example, the classical form of the universal approximation theorem for
arbitrary width and bounded depth coincides with the above theorem provided
that the activation function is not polynomial [11, 22, 45]. In this case, we
only consider the class of two-layer neural networks ($L=2$), and $n$ can be
interpreted as the index related to the number of neurons in the hidden layer.
On the other hand, in the recent work of [26], the authors discussed the dual
scenario: the networks of bounded width and arbitrary depth. The authors
assumed that the activation function is a non-affine continuous function that
is continuously differentiable at at least one point, with a nonzero
derivative at that point. Under this assumption, they considered the space of
feed-forward neural networks with an arbitrary number of hidden layers each
with $d+D+2$ neurons, and proved (2.16). In this case, $n_{i}=d+D+2$ for all
$i\in\\{1,\cdots,L-1\\}$ and $n$ can be regarded as the index concerning the
number of hidden layers ($n=L$).

In this paper, we will consider the bounded and continuous target functions.
More precisely, in Section 2.3, the function we aim to approximate by neural
networks is $\alpha^{*}(\omega)$ which is the solution of (2.7), where $S$,
$M$ are defined as before, but each component of the load vector
$F_{j}(\omega)=\int_{I}f(x;\omega)\phi_{j}(x)\,\mathrm{d}x$ is the function of
$\omega\in\Omega$. By the assumption (2.9), we know that $\alpha^{*}(\omega)$
is continuous and bounded. Therefore, it suffices to consider the family of
bounded neural networks as our ansatz space.

Let us assume that the target function $\alpha^{*}$ in (2.8) has the bound
$\sup_{\omega\in\Omega}|\alpha^{*}(\omega)|<C_{\alpha}$, and let $\sigma^{*}$
be a bounded activation function with continuous inverse (e.g. sigmoid or
hyperbolic tangent activation functions) and
$\sigma_{-}:=\inf_{x\in\mathbb{R}}\sigma^{*}(x)$ and
$\sigma_{+}:=\sup_{x\in\mathbb{R}}\sigma^{*}(x)$. Next, we define the affine
linear function $h:[\sigma_{-},\sigma_{+}]\rightarrow[-C_{\alpha},C_{\alpha}]$
with $h(\sigma_{-})=-C_{\alpha}$, $h(\sigma_{+})=C_{\alpha}$. Note that the
function $h$ also has a continuous inverse.

We then consider the following class of neural networks:

$\mathcal{N}_{n}=\left\\{h\circ\sigma^{*}(g_{n}):g_{n}\in\mathcal{N}_{n}^{*}\right\\}.$
(2.17)

Then it is obvious that

$\|g\|_{C(\Omega)}\leq C_{\alpha},\quad\forall
g\in\bigcup_{n\in\mathbb{N}}\mathcal{N}_{n}.$ (2.18)

By the definitions of $h$ and $\sigma^{*}$, $(\sigma^{*})^{-1}\circ
h^{-1}(\alpha^{*})$ is well-defined, and is continuous on the compact set
$\Omega$. Then by Theorem 2.1, we have

$\lim_{n\rightarrow\infty}\inf_{\hat{g}\in\mathcal{N}^{*}_{n}}\bigg{\|}\hat{g}-(\sigma^{*})^{-1}\circ
h^{-1}(\alpha^{*})\bigg{\|}_{C(\Omega)}=0.$ (2.19)

Then from (2.19), we obtain by the continuity of the function
$h\circ\sigma^{*}(\cdot)$ that

$\lim_{n\rightarrow\infty}\inf_{\hat{g}\in\mathcal{N}^{*}_{n}}\bigg{\|}h\circ\sigma^{*}(\hat{g})-\alpha^{*}\bigg{\|}_{C(\Omega)}=0.$

Therefore, by the definitions of $\mathcal{N}_{n}$ and $\mathcal{N}^{*}_{n}$,
we have that

$\lim_{n\rightarrow\infty}\inf_{\tilde{g}\in\mathcal{N}_{n}}\|\tilde{g}-\alpha^{*}\|_{C(\Omega)}=\lim_{n\rightarrow\infty}\inf_{\hat{g}\in\mathcal{N}^{*}_{n}}\bigg{\|}h\circ\sigma^{*}(\hat{g})-\alpha^{*}\bigg{\|}_{C(\Omega)}=0.$
(2.20)

The above argument is summarized in the following theorem concerning the
universal approximation properties for bounded neural networks.

######  Theorem 2.2.

Assume that $\Omega$ is compact and
$\alpha^{*}:\Omega\rightarrow\mathbb{R}^{N-1}$ is a continuous function
defined in (2.8). Then we have

$\lim_{n\rightarrow\infty}\inf_{\tilde{g}\in\mathcal{N}_{n}}\|\tilde{g}-\alpha^{*}\|_{C(\Omega)}=0.$
(2.21)

###  2.3 Problem formulation

In this section, we shall state the problem under consideration in detail, and
describe the mathematical setting and the strategy for the analysis.
Henceforth, for given external force $f(x;\omega)$ parametrized by
$\omega\in\Omega$, $u(x;\omega)$ denotes the corresponding true solution (weak
solution) of the equation (2.1) either with Dirichlet boundary condition or
Neumann boundary condition, and $u_{N}$ signifies the Legendre–Galerkin
approximation

$u_{N}(x,\omega)=\sum_{k=1}^{N-1}\alpha^{*}_{k}(\omega)\phi_{k}(x),$ (2.22)

where $\\{\alpha^{*}_{k}(\omega)\\}$ is the set of coefficients computed by
(2.8). Here we shall use the same notations $S$, $M$ and $F(\omega)$ as
defined in (2.6).

Note that both $u$ and $u_{N}$ depend on the external force $f(x;\omega)$ and
hence, depend on the parametric variable $\omega\in\Omega$. In particular,
$\alpha^{*}$ can be regarded as a function contained in
$C(\Omega,\mathbb{R}^{N-1})$ provided that $f(x;\omega)$ satisfies the
assumption (2.9) . In this regard, we can interpret $\alpha^{*}$ as the
solution to the minimization problem which is equivalent to (2.5):

$\alpha^{*}=\operatorname*{arg\,min}_{\alpha\in
C(\Omega,\mathbb{R}^{N-1})}\mathcal{L}(\alpha),$ (2.23)

where $\mathcal{L}$ is the population loss function defined in (2.11).

Next, we define the neural network approximation of $\alpha^{*}$. More
precisely, we look for $\widehat{\alpha}(n):\Omega\rightarrow\mathbb{R}^{N-1}$
solving the continuous residual minimization problem

$\widehat{\alpha}(n)=\operatorname*{arg\,min}_{\alpha\in\mathcal{N}_{n}}\mathcal{L}(\alpha),$
(2.24)

where the minimizer was found over the family of neural networks
$\mathcal{N}_{n}$ for $n\in\mathbb{N}$, and we shall write the associate
solution by

$u_{N,n}(x;\omega)=\sum^{N-1}_{k=1}\widehat{\alpha}(n)_{k}(\omega)\phi_{k}(x).$
(2.25)

Finally, let us define the solution to the discrete residual minimization
problem, namely,

$\widehat{\alpha}(n,M)=\operatorname*{arg\,min}_{\alpha\in\mathcal{N}_{n}}\mathcal{L}^{M}(\alpha),$
(2.26)

where the empirical loss $\mathcal{L}^{M}$ was defined in (2.12). In addition,
we shall write the corresponding solution as

$u_{N,n,M}(x;\omega)=\sum^{N-1}_{k=1}\widehat{\alpha}(n,M)_{k}(\omega)\phi_{k}(x).$
(2.27)

In the present paper, we shall ignore the error caused by the minimization,
and assume that we can always find the exact minimizer for the minimization
problems (2.24) and (2.26). Therefore, we shall employ $\widehat{\alpha}(n,M)$
as the neural network approximation for the exact coefficient $\alpha^{*}$,
and $u_{N,n,M}$ is our approximate solution computed by the numerical scheme
described in Section 2.1. If we were to summarize once more, $N$ stands for
the parameter for the Legendre–Galerkin method (number of basis functions),
$n$ is the parameter for the neural network architecture, and $M$ denotes the
number of sampling points when we compute the discrete loss function.

The main objective of this paper is to show the convergence

$\|u-u_{N,n,M}\|_{L^{2}(\Omega;H^{1}(I))}\rightarrow
0\quad{\rm{as}}\,\,N,n,M\rightarrow\infty.$ (2.28)

For this purpose, we split the error in (2.28) into three parts:

$u-u_{N,n,M}=(u-u_{N})+(u_{N}-u_{N,n})+(u_{N,n}-u_{N,n,M}).$ (2.29)

The first term on the right-hand side of (2.29) is the error for the
Legendre–Galerkin approximation. It is widely known that the Legendre–Galerkin
approximation error is sufficiently small in the numerical analysis, and hence
this error can be ignored in this paper. In fact, it is well known that if the
solution $u$ is in $H^{m}$ for some $m\in\mathbb{N}$, then the corresponding
approximate solution $u_{N}$ is close to $u$ within $\mathcal{O}(N^{-m})$
accuracy. Hence, if the solution $u^{*}$ is sufficiently smooth, the
Legendre–Galerkin method achieves exponential convergence. Furthermore, as $N$
increases, the condition number of the capacitance matrix $S+\nu M$ grows, and
hence, the large $N$ may deteriorate numerical convergence in practice. See,
for example, Table III in [49] where $N$ increased from $64$ to $128$, but the
computed $L^{\infty}$ error increased slightly. Therefore, a suitable choice
of $N$ guarantees that the approximate error is negligible, which justifies
our assumption. In the rest of the paper, we will assume that we have found
sufficiently good $N$ and we will fix that $N$ in the analysis.

The second term is called the approximation error, since it occurs when we
approximate the target function with the functions in our ansatz space. The
third term is often referred to as the generalization error, which is the
error of the neural network-based approximation on predicting unseen data.

We will first perform the convergence analysis of both approximation error and
generalization error for the corresponding coefficients neural networks, i.e.,

$\|\alpha^{*}-\widehat{\alpha}(n)\|_{L^{2}(\Omega)}+\|\widehat{\alpha}(n)-\widehat{\alpha}(n,M)\|_{L^{2}(\Omega)}.$
(2.30)

This error estimate for the neural networks is the crucial step in our
analysis. We will discuss the convergence of the error (2.30) in Section 3,
and based on this result, we will conduct the convergence analysis of the
approximate solution (2.28) in Section 4.

##  3 Abstract convergence analysis

In this section, we shall perform an abstract convergence analysis for the
error (2.30). We begin with an observation on the structure of the loss
functions. Note by the definition of $\mathcal{L}(\alpha)$ and (2.4) that

$\displaystyle\mathcal{L}(\alpha)$
$\displaystyle=\int_{\Omega}\sum^{N-1}_{j=1}\bigg{|}\sum^{N-1}_{k=1}\alpha_{k}(\omega)\left(\int_{I}\phi^{\prime}_{k}(x)\phi^{\prime}_{j}(x)\,\mathrm{d}x\right)-\int_{I}f(x;\omega)\phi_{j}(x)\bigg{|}^{2}\,\mathrm{d}\omega$
(3.1)
$\displaystyle=\int_{\Omega}\sum^{N-1}_{j=1}\bigg{|}\left(A\alpha(\omega)\right)_{j}-\left(F(\omega)\right)_{j}\bigg{|}^{2}\,\mathrm{d}\omega=\left\lVert
A\alpha(\omega)-F(\omega)\right\rVert^{2}_{L^{2}(\Omega)},$

where $A:=S+\nu M$ and $F(\omega)$ are defined in (2.6). Similarly, we also
have

$\mathcal{L}^{M}(\alpha)=\frac{|\Omega|}{M}\sum_{m=1}^{M}|A\alpha(\omega_{m})-F(\omega_{m})|^{2}.$
(3.2)

Here the matter on the loss functions $\mathcal{L}$ and $\mathcal{L}^{M}$ has
been converted into the analysis of the residual induced by the matrix
equation defined by $A$. Note that the matrix $A$ depends on the structure of
the given differential equations, boundary conditions and the choice of basis
functions, and hence the general analysis of the matrix $A$ that can cover
various PDE settings and numerical methods is important. In the next
proposition, we shall present the result on the structure of the matrix $A$
which can cover a wide range of settings of interest.

######  Proposition 3.1.

Assume that $A\in\mathbb{R}^{d\times d}$ is a symmetric, non-singular matrix
and let $\rho_{\min}=\min_{i}\\{|\lambda_{i}|\\}$,
$\rho_{\max}=\max_{i}\\{|\lambda_{i}|\\}$ where $\\{\lambda_{i}\\}$ is the set
of eigenvalues of $A$. Then for any $x\in\mathbb{R}^{d}$ we have

$\rho_{\min}|x|\leq|Ax|\leq\rho_{\max}|x|.$ (3.3)

###### Proof.

By the spectral theorem, we can choose the set $\\{x_{1},\cdots,x_{d}\\}$,
where $x_{j}$ is the corresponding eigenvector for the eigenvalue
$\lambda_{j}$ of $A$, and $\\{x_{1},\cdots,x_{d}\\}$ forms an orthonormal
basis of $\mathbb{R}^{d}$. Hence for arbitrary $x\in\mathbb{R}^{d}$, we can
write $x=\sum^{d}_{j=1}a_{j}x_{j}$. for $a_{j}\in\mathbb{R}$. Then by the
orthonormality and the definition of eigenvalue, we have

$|Ax|^{2}=\bigg{|}A\sum^{d}_{j=1}a_{j}x_{j}\bigg{|}^{2}=\bigg{|}\sum^{d}_{j=1}a_{j}Ax_{j}\bigg{|}^{2}=\bigg{|}\sum^{d}_{j=1}a_{j}\lambda_{j}x_{j}\bigg{|}^{2}=\sum^{d}_{j=1}(a_{j}\lambda_{j})^{2}.$

Since $A$ is non-singular, all the eigenvalues of $A$ are nonzero, and from
the above inequality, we can conclude that

$\rho^{2}_{\min}|x|^{2}\leq|Ax|^{2}\leq\rho^{2}_{\max}|x|^{2},$

which completes the proof. ∎

Throughout the section, we will utilize the matrix $A=S+\nu M$ defined in
(2.6), and assume that $A$ is symmetric and non-singular. As it will be made
clear in the next section, $A$ is indeed symmetric and non-singular in many
cases including the situation under consideration in this paper. Henceforth,
we set $\Omega=[a,b]^{d}$ for $a,b\in\mathbb{R}^{+}$.

###  3.1 Approximation error

This subsection aims to prove that the approximation error for the neural
networks converges to zero, which is encapsulated in the following theorem.
The idea of the proof is to make a connection between the error and the loss
function evaluated at the minimizer, which actually converges to zero as the
network architecture grows.

######  Theorem 3.2.

Assume that (2.9) holds. For $\alpha^{*}\in C(\Omega,\mathbb{R}^{N-1})$
defined in (2.23) and $\widehat{\alpha}(n)\in\mathcal{N}_{n}$ defined from
(2.24), we have that

$\|\alpha^{*}-\widehat{\alpha}(n)\|_{L^{2}(\Omega)}\rightarrow
0\quad{\rm{as}}\,\,n\rightarrow\infty.$

###### Proof.

Since $A$ is symmetric and non-singular, by Proposition 3.1, we have

$\|\alpha^{*}-\widehat{\alpha}(n)\|^{2}_{2}\lesssim\|A\alpha^{*}-A\widehat{\alpha}(n)\|^{2}_{2}\lesssim\|A\alpha^{*}-F\|^{2}_{2}+\|A\widehat{\alpha}(n)-F\|^{2}_{2}=\mathcal{L}(\widehat{\alpha}(n)),$
(3.4)

where we have used (2.7). Since $\widehat{\alpha}(n)$ is the minimizer of
$\mathcal{L}(\cdot)$ over $\mathcal{N}_{n}$, we obtain

$\displaystyle\mathcal{L}(\widehat{\alpha}(n))$
$\displaystyle\leq\inf_{\alpha\in\mathcal{N}_{n}}\mathcal{L}(\alpha)=\inf_{\alpha\in\mathcal{N}_{n}}\|A\alpha-F\|^{2}_{2}$
$\displaystyle\lesssim\inf_{\alpha\in\mathcal{N}_{n}}\left(\|A\alpha-A\alpha^{*}\|^{2}_{2}+\|A\alpha^{*}-F\|^{2}_{2}\right)$
$\displaystyle\lesssim\inf_{\alpha\in\mathcal{N}_{n}}\|\alpha-\alpha^{*}\|^{2}_{2},$

where we have used (2.7) and Proposition 3.1 once more. Then by the universal
approximation property (Theorem 2.2),
$\inf_{\alpha\in\mathcal{N}_{n}}\|\alpha-\alpha^{*}\|^{2}_{2}\rightarrow 0$ as
$n\rightarrow\infty$ and hence $\mathcal{L}(\widehat{\alpha}(n))\rightarrow 0$
as $n\rightarrow\infty$. Together with (3.4), this completes the proof. ∎

######  Remark 3.3.

In many cases, the explicit convergence rate for the approximation error can
be computed. See, for instance, [59, 60] where the approximation of continuous
functions with deep ReLU networks are discussed. However, the approximation
rates for continuous functions deteriorate as the dimension increases. This
can be tackled by assuming that the target function lies in a suitable smaller
function class which has low complexity compared to the functions with
classical regularity. One function class of this type is the so-called Barron
space introduced in the seminal work [3]. If the target function $\alpha^{*}$
belongs to the Barron space, then approximating it by two-layer networks has
the explicit convergence rate $\mathcal{O}(n^{-\frac{1}{2}})$. See also [1,
27, 14, 54, 55, 37] for some variants of Barron spaces and their approximation
properties by two-layer neural networks. In our case, if the external force
$f(x;\omega)$ is, for example, Gaussian, positive definite, radial, or smooth
with respect to $\omega\in\Omega$ (see [3, 14] where the characterization of
Barron functions is discussed), then $F(\omega)$ is the Barron function, and
hence our target function $\alpha^{*}(\omega)=A^{-1}F(\omega)$ becomes the
Barron function.

###  3.2 Generalization error

To handle the generalization error, we need to introduce the quantity so-
called Rademacher complexity [4]. For more information on the Rademacher
complexity and related topics, see [18, 57, 48].

######  Definition 3.4.

For a family $\\{X_{i}\\}_{i=1}^{M}$ of i.i.d. random variables according to
$\mathbb{P}_{\Omega}$, we define the Rademacher complexity of the function
class $\mathcal{F}$ by

$R_{M}(\mathcal{F})=\mathbb{E}_{\\{X_{i},\varepsilon_{i}\\}^{M}_{i=1}}\bigg{[}\sup_{f\in\mathcal{F}}\bigg{|}\frac{1}{M}\sum^{M}_{i=1}\varepsilon_{i}f(X_{i})\bigg{|}\bigg{]},$

where $\varepsilon_{i}$’s are i.i.d. Bernoulli random variables, i.e.,
$\mathbb{P}(\varepsilon_{i}=1)=\mathbb{P}(\varepsilon_{i}=-1)=\frac{1}{2}$.

Note from the above definition that the Rademacher complexity is the average
of the maximum correlation between the vector $(f(X_{1}),\cdots,f(X_{M}))$ and
the random noise $(\varepsilon_{1},\cdots,\varepsilon_{M})$, where the maximum
is chosen over the function class $\mathcal{F}$. Intuitively, the Rademacher
complexity of the function class $\mathcal{F}$ measures the ability of
functions from $\mathcal{F}$ to fit random noise. If we can always find a
function which has a high correlation with a randomly selected noise vector,
the corresponding function class can be interpreted as too large for
statistical purposes. Conversely, if the Rademacher complexity decays as the
sample size $M$ increases, then it is impossible to find a function that
highly correlates in expectation with random noise.

We now make the precise connection between Rademacher complexity and the
generalization error. In particular, we shall use the result that for any
bounded function class $\mathcal{F}$, the condition
$R_{M}({\mathcal{F}})=o(1)$ implies the convergence of the generalization
error, which is called the uniform law of large numbers. In the following
theorem, the function class is assumed to be $b$-uniformly bounded, meaning
that $\|f\|_{\infty}\leq b$ for all $f\in\mathcal{F}$.

######  Theorem 3.5.

[Theorem 4.10 in [57]] Let $\mathcal{F}$ be a $b$-uniformly bounded class of
functions and $M\in\mathbb{N}$. Then for any small number $\delta>0$, we have

$\sup_{f\in\mathcal{F}}\bigg{|}\frac{1}{M}\sum^{M}_{i=1}f(X_{i})-\mathbb{E}[f(X)]\bigg{|}\leq
2R_{M}(\mathcal{F})+\delta,$

with probability at least $1-\exp(-\frac{M\delta^{2}}{2b^{2}})$.

Let us now define the function class

$\mathcal{F}_{n}:=\\{|A\alpha-F|^{2}:\alpha\in\mathcal{N}_{n}\\},$ (3.5)

where $A$ and $F$ are define in (2.6). Then by Proposition 3.1, we have that

$\|A\alpha-F\|_{L^{\infty}(\Omega)}\leq\|A\alpha\|_{L^{\infty}(\Omega)}+\|F\|_{L^{\infty}(\Omega)}\lesssim\|\alpha\|_{L^{\infty}(\Omega)}+\|F\|_{L^{\infty}(\Omega)}.$

Therefore, by (2.18) and (2.9), we observe that for all $n\in\mathbb{N}$,
$\mathcal{F}_{n}$ is $\overline{G}$-uniformly bounded for some positive
constant $\overline{G}>0$. Then the following lemma is the direct consequence
of Theorem 3.5.

######  Lemma 3.6.

Let $\\{\omega_{m}\\}^{M}_{m=1}$ be i.i.d. samples randomly chosen from the
distribution $\mathbb{P}_{\Omega}$ in the definition of the empirical loss
function (2.12). Then for any small $\delta>0$, we have with probability at
least $1-2\exp(-\frac{M\delta^{2}}{32\overline{G}^{2}})$

$\sup_{\alpha\in\mathcal{N}_{n}}\left|\mathcal{L}^{M}(\alpha)-\mathcal{L}(\alpha)\right|\leq
2R_{n}(\mathcal{F}_{n})+\frac{\delta}{2}.$ (3.6)

With Lemma 3.6, we now establish the convergence concerning
$M\rightarrow\infty$ for the generalization error, which is encapsulated in
the following theorem. Here we assume that the Rademacher complexity of
$\mathcal{F}_{n}$ converges to zero for all $n\in\mathbb{N}$, which is the
standard assumption as discussed above.

######  Theorem 3.7.

Assume that (2.9) holds, and suppose that
$\lim_{M\rightarrow\infty}R_{M}(\mathcal{F}_{n})=0$ for all $n\in\mathbb{N}$.
Then we have with probability $1$ over i.i.d. samples that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\widehat{\alpha}(n,M)-\widehat{\alpha}(n)\|^{2}_{2}=0.$

###### Proof.

By Proposition 3.1 and the definition (2.24), we have

$\displaystyle\|\widehat{\alpha}(n)-\widehat{\alpha}(n,M)\|^{2}_{2}$
$\displaystyle\lesssim\|A\widehat{\alpha}(n)-A\widehat{\alpha}(n,M)\|^{2}_{2}\lesssim\|A\widehat{\alpha}(n)-F\|^{2}_{2}+\|A\widehat{\alpha}(n,M)-F\|^{2}_{2}$
(3.7)
$\displaystyle=\mathcal{L}(\widehat{\alpha}(n))+\mathcal{L}(\widehat{\alpha}(n,M))\leq
2\mathcal{L}(\widehat{\alpha}(n,M))$

We shall now apply Lemma 3.6 with $\delta=2M^{-\frac{1}{2}+\varepsilon}$ for
$0<\varepsilon<\frac{1}{2}$. Then with probability
$1-2\exp(-\frac{M^{2\varepsilon}}{8\overline{G}^{2}})$, we have by minimality
of $\widehat{\alpha}(n,M)$,

$\mathcal{L}(\widehat{\alpha}(n,M))\leq\mathcal{L}^{M}(\widehat{\alpha}(n,M))+2R_{M}(\mathcal{F}_{n})+M^{-\frac{1}{2}+\varepsilon}\leq\mathcal{L}^{M}(\widehat{\alpha}(n))+2R_{M}(\mathcal{F}_{n})+M^{-\frac{1}{2}+\varepsilon}$

If we apply Lemma 3.6 once more, we obtain

$\mathcal{L}(\widehat{\alpha}(n,M))\leq\mathcal{L}(\widehat{\alpha}(n))+4R_{M}(\mathcal{F}_{n})+2M^{-\frac{1}{2}+\varepsilon}.$

Therefore, if we take $M\rightarrow\infty$ on (3.7), we have by the assumption
on $R_{M}(\mathcal{F}_{n})$ that

$\lim_{M\rightarrow\infty}\|\widehat{\alpha}(n,M)-\widehat{\alpha}(n)\|^{2}_{2}\lesssim\mathcal{L}(\widehat{\alpha}(n)).$

As before, from Theorem 2.2, we observe by taking $n\rightarrow\infty$ that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\widehat{\alpha}(n,M)-\widehat{\alpha}(n)\|^{2}_{2}\lesssim\lim_{n\rightarrow\infty}\mathcal{L}(\widehat{\alpha}(n))\lesssim\lim_{n\rightarrow\infty}\inf_{\alpha\in\mathcal{N}_{n}}\|\alpha-\alpha^{*}\|^{2}_{2}\rightarrow
0,$

and the proof is completed. ∎

######  Remark 3.8.

The assumption $\lim_{M\rightarrow\infty}R_{M}(\mathcal{F}_{n})=0$ in Theorem
3.7 is common in the statistical learning theory and known for several
function classes [43, 14]. Moreover, in many cases, it was shown that the
Rademacher complexity converges to zero as $M$ goes to infinity with the
convergence rate $\mathcal{O}(M^{-\frac{1}{2}})$, which is same as the rate
for the classical Monte Carlo integration; see, for example, two-layer ReLU
networks [43], spectral Barron setting with the softplus activation function
[37], general Barron class [14] and recent result on the complexity estimate
for general case [21].

##  4 Application to linear second-order elliptic PDEs

In this section, we will derive the convergence of approximate solutions for
the specific PDE problems based on what we have discussed in the previous
sections. In general, one of the most challenging tasks for neural network-
based PDE techniques is to impose the boundary condition. One typical way to
handle this problem is to add penalty terms to loss functions. However, adding
the penalty term to the loss function make it difficult to analyze the
algorithm and may deteriorate the convergence.

As mentioned earlier, one important advantage of using the ULGNet is that we
can impose exact boundary values to the approximate solutions. More precisely,
once the boundary condition is given (either Dirichlet or Neumann),
$\\{a_{k},b_{k}\\}_{k=1}^{N-1}$ in (2.2) are uniquely determined and hence,
the basis function (2.2) is fixed with the exact boundary conditions. Then we
compute the stiffness matrix $S$ and the mass matrix $M$ using (2.6), and
hence the total matrix $A$ is obtained accordingly. Then we can apply the
abstract analysis studied in the previous section. One interesting thing is
that the analysis in Section 3 can cover most of the self-adjoint PDEs. In the
following subsection, we will deal with both the Dirichlet case and the
Neumann case with detailed computation, to highlight the flexibility of the
proposed method. As mentioned earlier, for the sake of simplicity we shall
consider the self-adjoint equation

$-u^{\prime\prime}(x;\omega)+\nu
u(x;\omega)=f(x;\omega)\quad{\rm{for}}\,\,x\in I:=(-1,1),\,\,\omega\in\Omega,$
(4.1)

with the homogeneous boundary conditions

$u(-1;\omega)=u(1;\omega)=0\quad{\rm{or}}\quad
u^{\prime}(-1;\omega)=u^{\prime}(1;\omega)=0\quad\forall\omega\in\Omega.$
(4.2)

Note, however, that our analysis can be extended to the case of any self-
adjoint PDEs and non-homogeneous boundary conditions in a straightforward
manner. After that, we will also discuss the way to extend the analysis to the
higher-dimensional case. Throughout the section, we will use the same
notations, assumptions, and mathematical settings introduced in Section 2.

###  4.1 Main theorem and convergence of approximate solutions

For the input of neural networks, we shall assume that the forcing functions
are generated with respect to the parameter $\omega$ uniformly distributed
over the $d$-dimensional unit cube, i.e. the external force is regarded as a
bivariate function $f(x;\omega)$ where the spatial variable $x$ is defined
over $I$ and the parametric variable $\omega$ is uniformly distributed over
$[a,b]^{d}$ for $a,b\in\mathbb{R}^{+}$. Note that as the dimension $d$
increases, our method can cover a wider range of external force functions. As
mentioned earlier, we assume that (2.9) holds.

We begin with the following lemma, concerning the values of $a_{k}$ and
$b_{k}$, which is a direct consequence of Lemma 4.1 in [50].

######  Lemma 4.1.

Let $L_{k}$ be the $k$-th Legendre polynomial. For each $k\geq 0$, there is a
unique set of $\\{a_{k},b_{k}\\}^{N-1}_{k=1}$ such that
$\phi_{k}(x)=L_{k}(x)+a_{k}L_{k+1}(x)+b_{k}L_{k+2}(x)$ satisfies the given
boundary condition (either Dirichlet or Neumann). In particular,

  * •

homogeneous Dirichlet boundary condition $(\phi_{k}(-1)=\phi_{k}(1)=0$ for all
$k\geq 1)$:

$a_{k}=0\quad{\rm{and}}\quad b_{k}=-1\quad\forall k\geq 1.$ (4.3)

  * •

homogeneous Neumann boundary condition
$(\phi^{\prime}_{k}(-1)=\phi^{\prime}_{k}(1)=0$ for all $k\geq 1)$:

$a_{k}=0\quad{\rm{and}}\quad b_{k}=\frac{-k(k+1)}{(k+2)(k+3)}\quad\forall
k\geq 1.$ (4.4)

We will also use the following lemma, regarding the structure of the stiffness
and the mass matrices, again quoted from [50].

######  Lemma 4.2.

The stiffness matrix $S$ defined in (2.6) is a diagonal matrix with the
diagonal entries

$S_{kk}=-(4k+6)b_{k}\quad\forall k\geq 1.$ (4.5)

Furthermore, the mass matrix $M$ is a symmetric penta-diagonal matrix with the
following non-zero elements:

$M_{jk}=M_{kj}=\begin{cases}\frac{2}{2k+1}+\frac{2a_{k}^{2}}{2k+3}+\frac{2b^{2}_{k}}{2k+5},\quad&j=k,\\\
\frac{2a_{k}}{2k+3}+\frac{2a_{k+1}b_{k}}{2k+5},\quad&j=k+1,\\\
\frac{2b_{k}}{2k+5},\quad&j=k+2.\end{cases}$ (4.6)

Note from Lemma 4.1 and Lemma 4.2 that all entries of the stiffness and mass
matrices are finite, since the basis functions we used are continuous over the
compact interval $[-1,1]$, and hence the integrals in (2.6) are well-defined.

Let us write the stiffness and the mass matrices as $S_{\rm{Dir}}$,
$M_{\rm{Dir}}$ for the Dirichlet case computed by (4.3), (4.5) and (4.6), and
define $A_{\rm{Dir}}=S_{\rm{Dir}}+\nu M_{\rm{Dir}}$. The matrices for the
Neumann case $S_{\rm{Neu}}$, $M_{\rm{Neu}}$ and $A_{\rm{Neu}}$ are similarly
obtained by using (4.4), (4.5) and (4.6). From now on, $A_{\rm{BC}}$ means
either $A_{\rm{Dir}}$ or $A_{\rm{Neu}}$ and the matrix denoted by $A$ in the
previous sections is now considered to be $A_{\rm{BC}}$. By Lemma 4.2,
$A_{\rm{BC}}$ is symmetric. Furthermore for any non-zero
$y=(y_{1},\cdots,y_{N-1})\in\mathbb{R}^{N-1}$,

$y\cdot
A_{\rm{BC}}y=\sum^{N-1}_{i=1}\sum^{N-1}_{j=1}y_{i}y_{j}\int_{I}\left(\phi^{\prime}_{i}(x)\phi^{\prime}_{j}(x)+\nu\phi_{i}(x)\phi_{j}(x)\right)\,\mathrm{d}x\geq\int_{I}\sum^{N-1}_{i=1}y_{i}^{2}|\phi^{\prime}_{i}(x)|^{2}\,\mathrm{d}x=\sum^{N-1}_{i=1}y_{i}^{2}S_{ii}>0,$
(4.7)

which means that $A_{\rm{BC}}$ is positive-definite, and hence is non-
singular. Therefore, we can apply the analysis conducted in Section 3. We now
consider that the loss functions $\mathcal{L}$ and $\mathcal{L}^{M}$ ((3.1)
and (3.2) respectively) where the matrix $A$ is now replaced by $A_{\rm{BC}}$,
and the corresponding coefficients defined through (2.23), (2.24) and (2.26).
Then Theorem 3.2 and Theorem 3.7 together with the triangular inequality give
us that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\alpha^{*}-\widehat{\alpha}(n,M)\|_{L^{2}(\Omega)}=0.$
(4.8)

Now we shall state and prove the main theorem of this paper. Recall that the
definition of approximate solutions (2.22), (2.25) and (2.27).

######  Theorem 4.3.

Assume that (2.9) holds, and suppose that for any $n\in\mathbb{N}$,
$R_{M}(\widetilde{{\mathcal{F}}}_{n})\rightarrow 0$ as $M\rightarrow\infty$,
where
$\widetilde{{\mathcal{F}}}_{n}:=\\{|A_{\rm{BC}}\alpha-F|^{2}:\alpha\in\mathcal{N}_{n}\\}$.
Then for fixed $N\in\mathbb{N}$, we have with probability $1$ over i.i.d.
samples that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|u_{N}-u_{N,n,M}\|_{L^{2}(\Omega;H^{1}(I))}=0.$
(4.9)

###### Proof.

Let us first prove the convergence with respect to the $H^{1}$-seminorm. By
the orthogonality of $\\{\phi^{\prime}_{k}\\}_{k=1}^{N-1}$ (Lemma 4.1, Lemma
4.2), and Theorem 3.2, Theorem 3.7, we note that

$\displaystyle\int_{\Omega}\int_{I}|u^{\prime}_{N}-u^{\prime}_{N,n,M}|^{2}\,\mathrm{d}x\,\mathrm{d}\omega$
$\displaystyle=\int_{\Omega}\int_{I}\bigg{|}\sum^{N-1}_{k=1}\left(\alpha^{*}_{k}-\widehat{\alpha}(n,M)_{k}\right)\phi^{\prime}_{k}\bigg{|}^{2}\,\mathrm{d}x\,\mathrm{d}\omega$
$\displaystyle=\int_{\Omega}\sum^{N-1}_{k=1}|\alpha^{*}_{k}-\widehat{\alpha}(n,M)_{k}|^{2}S_{kk}\,\mathrm{d}\omega$
$\displaystyle\lesssim\|\alpha^{*}-\widehat{\alpha}(n,M)\|^{2}_{L^{2}(\Omega)}\rightarrow
0\quad{\rm{as}}\,\,M,n\rightarrow\infty.$

For the Dirichlet case, the full $H^{1}$-convergence easily follows by the use
of Poincaré’s inequality. Here, however, to cover the Neumann case as well, we
present a more general argument. From the structure of the mass matrix (Lemma
4.1, Lemma 4.2) and Young’s inequality, we obtain that

$\displaystyle\int_{\Omega}\int_{I}|u_{N}-u_{N,n,M}|^{2}\,\mathrm{d}x\,\mathrm{d}\omega$
$\displaystyle=\int_{\Omega}\int_{I}\bigg{|}\sum^{N-1}_{k=1}\left(\alpha^{*}_{k}-\widehat{\alpha}(n,M)_{k}\right)\phi_{k}\bigg{|}^{2}\,\mathrm{d}x\,\mathrm{d}\omega$
$\displaystyle=\int_{\Omega}\sum^{N-1}_{k=1}|\alpha^{*}_{k}-\widehat{\alpha}(n,M)_{k}|^{2}M_{kk}\,\mathrm{d}\omega$
$\displaystyle\,\,\,\,\,\,\,+2\int_{\Omega}\sum^{N-3}_{j=1}|\alpha^{*}_{j}-\widehat{\alpha}(n,M)_{j}||\alpha^{*}_{j+2}-\widehat{\alpha}(n,M)_{j+2}|M_{j\,j+2}\,\mathrm{d}\omega$
$\displaystyle\lesssim\int_{\Omega}\sum^{N-1}_{k=1}|\alpha^{*}_{k}-\widehat{\alpha}(n,M)_{k}|^{2}\,\mathrm{d}\omega+\int_{\Omega}\sum^{N-3}_{j=1}|\alpha^{*}_{j}-\widehat{\alpha}(n,M)_{j}|^{2}\,\mathrm{d}\omega$
$\displaystyle\,\,\,\,\,\,\,+\int_{\Omega}\sum^{N-3}_{j=1}|\alpha^{*}_{j+2}-\widehat{\alpha}(n,M)_{j+2}|^{2}\,\mathrm{d}\omega$
$\displaystyle\lesssim\|\alpha^{*}-\widehat{\alpha}(n,M)\|^{2}_{L^{2}(\Omega)}\rightarrow
0\quad{\rm{as}}\,\,M,n\rightarrow\infty,$

where we have used once more Theorem 3.2 and Theorem (3.7). ∎

###  4.2 Analysis in higher dimensions

In this section, we shall discuss the way to extend the ULGNet scheme and its
analysis to the multi-dimensional case for the spatial variable. Here we
illustrate the two-dimensional case with homogeneous Dirichlet boundary
condition, and the study for general-dimensional cases follows in a
straightforward manner.

Let us consider the following two-dimensional problem:

$\displaystyle-\Delta u+\nu u$
$\displaystyle=f\quad{\rm{in}}\,\,I^{2}=[-1,1]^{2},$ (4.10) $\displaystyle u$
$\displaystyle=0\quad{\rm{on}}\,\,\partial I^{2}.$ (4.11)

Here $f$ is a two-dimensional function, again parametrized by
$\omega\in\Omega$. The weak solution of the equation (4.10)-(4.11) is the
function $u\in H^{1}_{0}$ satisfying the weak formulation

$\int_{I^{2}}\left(\nabla u(x,y)\cdot\nabla v(x,y)+\nu
u(x,y)v(x,y)\right)\,\mathrm{d}x\,\mathrm{d}y=\int_{I^{2}}f(x,y)v(x,y)\,\mathrm{d}x\,\mathrm{d}y,\quad\forall
v\in H^{1}_{0}.$ (4.12)

For the Legendre–Galerkin approximation of (4.10)-(4.11), we define the two-
dimensional basis function by

$\phi_{i,j}(x,y)=\phi_{i}(x)\phi_{j}(y)\quad i,j\in\\{1,\cdots,N-1\\},$ (4.13)

where $\phi_{i}(x)$ and $\phi_{j}(y)$ are one-dimensional basis functions
defined in (2.2). Note that $(a_{k},b_{k})$’s are determined with the same
values computed in (4.3), since $\phi_{i,j}(x,y)$ is reduced to the single-
variable basis function if we fix one of $x$ or $y$, which needs to satisfy
the one-dimensional homogeneous Dirichlet boundary condition as before.

We shall rewrite $\phi_{i,j}(x,y)$, $i,j\in\\{1,\cdots,N-1\\}$ as
$\overline{\phi}_{k}(x,y)$, $k\in\\{1,\cdots,(N-1)^{2}\\}$, by allocating the
multi indices
$\\{(1,1),\cdots,(1,N-1),(2,1),\cdots,(2,N-1),\cdots,(N-1,1),\cdots,(N-1,N-1)\\}$
to $\\{1,\cdots,(N-1)^{2}\\}$. As before, for $i$,
$j\in\\{1,\cdots,(N-1)^{2}\\}$, we define the stiffness matrix, the mass
matrix and the load vector by

$\displaystyle\overline{S}$ $\displaystyle=(\overline{S}_{ij})$
$\displaystyle{\rm{where}}\,\,\overline{S}_{ij}:=\int_{I^{2}}\overline{\phi}^{\prime}_{i}(x,y)\overline{\phi}^{\prime}_{j}(x,y)\,\mathrm{d}x\,\mathrm{d}y,$
(4.14) $\displaystyle\overline{M}$ $\displaystyle=(\overline{M}_{ij})$
$\displaystyle{\rm{where}}\,\,\overline{M}_{ij}:=\int_{I^{2}}\overline{\phi}_{i}(x,y)\overline{\phi}_{j}(x,y)\,\mathrm{d}x\,\mathrm{d}y,$
$\displaystyle\overline{F}$ $\displaystyle=(\overline{F}_{j})$
$\displaystyle{\rm{where}}\,\,\overline{F}_{j}:=\int_{I^{2}}f(x,y)\overline{\phi}^{\prime}_{j}(x,y)\,\mathrm{d}x\,\mathrm{d}y,$

and we let $\overline{A}=\overline{S}+\nu\overline{M}$. Then we define the
loss functions

$\overline{\mathcal{L}}(\alpha)=\|\overline{A}\alpha(\omega)-\overline{F}(\omega)\|^{2}_{L^{2}(\Omega)}\quad{\rm{and}}\quad\overline{\mathcal{L}}^{M}(\alpha)=\frac{|\Omega|}{M}\sum^{M}_{m=1}|\overline{A}\alpha(\omega_{m})-\overline{F}(\omega_{m})|^{2},$
(4.15)

and as before, we further define the Legendre–Galerkin approximation

$\overline{\alpha}^{*}=\operatorname*{arg\,min}_{\alpha\in
C(\Omega,\mathbb{R}^{(N-1)^{2}})}\overline{\mathcal{L}}(\alpha),\quad\overline{u}_{N}=\sum^{(N-1)^{2}}_{k=1}\overline{\alpha}^{*}_{k}\overline{\phi}_{k}(x),$
(4.16)

the continuous residual minimization

$\overline{{\alpha}}(n)=\operatorname*{arg\,min}_{\alpha\in\overline{N}_{n}}\overline{\mathcal{L}}(\alpha),\quad\overline{u}_{N,n}=\sum^{(N-1)^{2}}_{k=1}\overline{{\alpha}}(n)_{k}\overline{\phi}_{k}(x),$
(4.17)

and the discrete residual minimization

$\overline{{\alpha}}(n,M)=\operatorname*{arg\,min}_{\alpha\in\overline{N}_{n}}\overline{\mathcal{L}}^{M}(\alpha),\quad\overline{u}_{N,n,M}=\sum^{(N-1)^{2}}_{k=1}\overline{{\alpha}}(n,M)_{k}\overline{\phi}_{k}(x),\\\
$ (4.18)

where $\overline{N}_{n}$ consists of neural networks with
$(N-1)^{2}$-dimensional output.

Next, by the definition, $\overline{A}$ is symmetric. Moreover, by the same
argument used in (4.7), $\overline{A}$ is positive-definite, and hence non-
singular. Therefore, we can apply Proposition 4.2, from which we can deduce
the convergence

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\overline{\alpha}^{*}-\overline{\alpha}(n,M)\|^{2}_{L^{2}(\Omega)}=0$

with the same arguments proposed in Theorem 3.2 and Theorem 3.7 with
$A=\overline{A}$. Again, by the same argument presented in Theorem 4.3 with
$A_{\rm{BC}}$ replaced by $\overline{A}$, we conclude that

$\lim_{n\rightarrow\infty}\lim_{M\rightarrow\infty}\|\overline{u}_{N}-\overline{u}_{N,n,M}\|_{L^{2}(\Omega;H^{1}(I))}=0.$
(4.19)

##  5 Numerical experiments

In this section, we shall demonstrate some numerical results to verify our
theoretical findings presented in the previous sections. We consider the
equation (2.1) with $\nu=1$, and impose the homogeneous Dirichlet boundary
condition on the interval $[-1,1]$:

$\displaystyle-u^{\prime\prime}+u$ $\displaystyle=f\quad{\rm{in}}\,\,(-1,1),$
(5.1) $\displaystyle u(-1)$ $\displaystyle=u(1)=0.$ (5.2)

For the external forcing term $f$, we randomly generate $M$ samples of the
training data of the ULGNet which is given by

$f(x;\omega_{m})=\omega^{1}_{m}\sin(2\pi\omega^{2}_{m}x)+\omega^{3}_{m}\cos(2\pi\omega^{4}_{m}x),\quad
1\leq m\leq M,$ (5.3)

where the random sample $\omega_{m}$ is selected from the uniform distribution

${\rm{i.i.d.}}\quad\omega_{m}=(\omega^{1}_{m},\omega^{2}_{m},\omega^{3}_{m},\omega^{4}_{m})\sim[0,1]^{4}=:\Omega,\quad
1\leq m\leq M.$

For the spectral element approximation of the equations (5.1)-(5.2), as
presented in Lemma 4.1 and Lemma 4.2, we will utilize the basis function (2.2)
and hence the matrices (4.5), (4.6) with

$a_{k}=0\quad{\rm{and}}\quad b_{k}=1\quad\forall k\geq 1.$

In our experiments, we set N = 32 and use the Legendre-Gauss-Lobatto nodal
points for numerical integrations. Then we solve the algebraic equation (2.7)
to compute $u_{N}$ which will be regarded as true solutions in the present
experiment. For the approximate solutions using the ULGNet, we minimize the
empirical loss function defined by (2.12) to compute the coefficients
$\widehat{\alpha}(n,M)$ and construct the predicted solution $U_{N,n,M}$ as
(2.27). We then measure the distance

$\|u_{N}-u_{N,n,M}\|_{L^{2}(\Omega;L^{2}(I))}$

by using the Monte–Carlo simulation with respect to $\Omega$ with $1000$
random samples. In this setting, we shall conduct the numerical experiments
with two different scenarios. One is to grow the neural network architecture
with fixed $M\in\mathbb{N}$, and the other one is to fix the neural network
architecture and to increase the number of training samples $M\in\mathbb{N}$.

More precisely, for the first experiment, we consider the two-layer neural
networks with $n\in\mathbb{N}$ neurons in the hidden layer, and use the
hyperbolic tangent activation function which satisfies the boundedness
condition assumed in Theorem 2.2. For the training, We fix the number of
training data $M=10000$ and adopt the limited-memory
Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimization algorithm as in [33],
with $20000$ training epochs. Then for each $n\in\mathbb{N}$ varying from $20$
to $300$, we train the model as described above and compute the relative
$L^{2}$-error between $u_{N}$ and $U_{N,n,M}$ evaluated over $1000$ randomly
generated functions of the form (5.3) that were not used in the training.

On the other hand, for the second experiment, we adopt a fully-connected
neural networks with five hidden layers with $128$ neurons for each layer,
again equipped with the hyperbolic tangent activation function. We fix this
neural network architecture and train the models with the same method
described above with the different number of training data $M\in\mathbb{N}$
for each model. In both experiments, all the computations were carried on by
Python software on the GPU cluster in the Department of Mathematics,
Sugnkyunkwan University. The results of the experiments are presented in
Figure 2.

(a) Relative $L^{2}$-error as $n$ increases.

(b) Relative $L^{2}$-error as $M$ increases.

Figure 2: Relative $L^{2}$-error evaluated over $1000$ randomly generated
samples as $n$, $M$ increase.

As shown in Figure 2, we can say that the relative $L^{2}$ error tends to
decrease as $n$ and $M$ increase. Since there are a number of factors in
training machine learning algorithms which may hold us back from obtaining
precise values (e.g. there is randomness for each realization of $\omega_{m}$
and optimization errors may occur), the trend of errors in the figure is not
expressed in a deterministic manner as the classical numerical methods for
PDEs. However, it can be clearly confirmed that from the figure that the error
tends to decrease in both cases, which supports what we theoretically proved
in previous sections.

##  6 Conclusion

In this paper, we have established the convergence of approximate solutions
obtained by the ULGNet algorithm to the weak solution of the linear, second-
order elliptic PDEs. The idea of ULGNet is to represent solutions as a linear
combination of the Legendre–Galerkin basis functions and generate the
coefficients through the deep neural network whose input is an external
forcing term, where the loss function is defined via the variational residual
minimization. The key observation for the proof is to convert the variational
minimization problem to the strong residual minimization solving matrix
equation, which is similar to the framework of classical numerical schemes.
Then we perform the analysis for the derived matrix, enabled by the choice of
our Legendre-Galerkin basis functions. As a consequence, together with the
universal approximation property of neural networks and the convergence of
Rademacher complexity, we showed the strong convergence of the minmizers to
the solutions of the given PDEs.

We expect that the numerical scheme and the analysis studied in the present
paper have considerable potential to be further investigated in various
directions. As the well-known basis functions are used for the scheme, from
the analysis point of view, it is more likely for us to apply various theories
from classical numerical analysis. One interesting future research direction
is to use different basis functions which might lead us to achieve better
computational cost and numerical tractability. On the other hand, it would be
of particular interest to extend the result to nonlinear equations with the
iterative methods, or the time-stepping analysis for the parabolic problem,
which will be addressed in the forthcoming papers.

## References

  * [1] F. Bach.  Breaking the curse of dimensionality with convex neural networks.  Journal of Machine Learning Research, 18(19):1–53, 2017.
  * [2] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner.  Learning data-driven discretizations for partial differential equations.  Proceedings of the National Academy of Sciences, 116(31):15344–15349, 2019.
  * [3] A. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information Theory, 39(3):930–945, 1993.
  * [4] P. L. Bartlett and S. Mendelson.  Rademacher and Gaussian complexities: risk bounds and structural results.  J. Mach. Learn. Res., 3(Spec. Issue Comput. Learn. Theory):463–482, 2002.
  * [5] J. Berner, P. Grohs, and A. Jentzen.  Analysis 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.  SIAM J. Math. Data Sci., 2(3):631–657, 2020.
  * [6] A. Biswas, J. Tian, and S. Ulusoy.  Error estimates for deep learning methods in fluid dynamics.  Numerische Mathematik, 151(3):753–777, Jul 2022.
  * [7] S. C. Brenner and L. R. Scott.  The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics.  Springer, New York, third edition, 2008.
  * [8] J. Choi, N. Kim, and Y. Hong.  Unsupervised legendre-galerkin neural network for stiff partial differential equations.  arXiv:2207.10241 [cs.LG], 2022.
  * [9] B. Chudomelka, Y. Hong, H. Kim, and J. Park.  Deep neural network for solving differential equations motivated by legendre-galerkin approximation.  arXiv:2010.12975 [math.NA], 2020.
  * [10] P. G. Ciarlet.  The finite element method for elliptic problems, volume 40 of Classics in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.  Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58 #25001)].
  * [11] G. Cybenko.  Approximation by superpositions of a sigmoidal function.  Math. Control Signals Systems, 2(4):303–314, 1989.
  * [12] T. De Ryck, A. D. Jagtap, and S. Mishra.  Error estimates for physics informed neural networks approximating the navier-stokes equations, 2022.
  * [13] C. Duan, Y. Jiao, Y. Lai, X. Lu, Q. Quan, and J. Z. Yang.  Analysis of deep ritz methods for laplace equations with dirichlet boundary conditions, 2021.
  * [14] W. E, C. Ma, and L. Wu.  The Barron space and the flow-induced function spaces for neural network models.  Constr. Approx., 55(1):369–406, 2022.
  * [15] W. E and B. Yu.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  Commun. Math. Stat., 6(1):1–12, 2018.
  * [16] L. C. Evans.  Partial differential equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2010.
  * [17] D. Gilbarg and N. S. Trudinger.  Elliptic partial differential equations of second order, volume 224 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, second edition, 1983.
  * [18] G. Gnecco and M. Sanguineti.  Approximation error bounds via Rademacher’s complexity.  Appl. Math. Sci. (Ruse), 2(1-4):153–176, 2008.
  * [19] D. Gottlieb and S. A. Orszag.  Numerical analysis of spectral methods: theory and applications.  CBMS-NSF Regional Conference Series in Applied Mathematics, No. 26. Society for Industrial and Applied Mathematics, Philadelphia, Pa., 1977.
  * [20] P. Grohs, F. Hornung, A. Jentzen, and P. von Wurstemberger.  A proof that artificial neural networks overcome the curse of dimensionality in the numerical approximation of black-scholes partial differential equations, 2018.
  * [21] Q. Hong, J. W. Siegel, and J. Xu.  A priori analysis of stable neural network solutions to numerical pdes, 2021.
  * [22] K. Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural Networks, 4(2):251–257, 1991.
  * [23] M. Hutzenthaler, A. Jentzen, T. Kruse, and T. A. Nguyen.  A proof that rectified deep neural networks overcome the curse of dimensionality in the numerical approximation of semilinear heat equations.  SN Partial Differential Equations and Applications, 1(2):10, Apr 2020.
  * [24] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6), 5 2021.
  * [25] E. Kharazmi, Z. Zhang, and G. E. M. Karniadakis.  $hp$-VPINNs: variational physics-informed neural networks with domain decomposition.  Comput. Methods Appl. Mech. Engrg., 374:Paper No. 113547, 25, 2021\.
  * [26] P. Kidger and T. Lyons.  Universal Approximation with Deep Narrow Networks.  In J. Abernethy and S. Agarwal, editors, Proceedings of Thirty Third Conference on Learning Theory, volume 125 of Proceedings of Machine Learning Research, pages 2306–2327. PMLR, 09–12 Jul 2020.
  * [27] J. M. Klusowski and A. R. Barron.  Approximation by combinations of relu and squared relu ridge functions with $\ell^{1}$ and $\ell^{0}$ controls.  IEEE Transactions on Information Theory, 64(12):7649–7656, 2018\.
  * [28] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  In F. Pereira, C. Burges, L. Bottou, and K. Weinberger, editors, Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
  * [29] I. Lagaris, A. Likas, and D. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [30] B. M. Lake, R. Salakhutdinov, and J. B. Tenenbaum.  Human-level concept learning through probabilistic program induction.  Science, 350(6266):1332–1338, 2015.
  * [31] Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  nature, 521(7553):436, 2015.
  * [32] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Fourier neural operator for parametric partial differential equations, 2020.
  * [33] D. C. Liu and J. Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical Programming, 45(1):503–528, Aug 1989.
  * [34] L. Lu, M. Dao, P. Kumar, U. Ramamurty, G. E. Karniadakis, and S. Suresh.  Extraction of mechanical properties of materials through deep learning from instrumented indentation.  Proceedings of the National Academy of Sciences, 117(13):7052–7062, 2020.
  * [35] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature machine intelligence, 3(3):218–229, 2021.
  * [36] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  DeepXDE: a deep learning library for solving differential equations.  SIAM Rev., 63(1):208–228, 2021.
  * [37] Y. Lu, J. Lu, and M. Wang.  A priori generalization analysis of the deep ritz method for solving high dimensional elliptic partial differential equations.  In M. Belkin and S. Kpotufe, editors, Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 3196–3241. PMLR, 15–19 Aug 2021.
  * [38] T. Luo and H. Yang.  Two-layer neural networks for partial differential equations: Optimization and generalization theory, 2020.
  * [39] X. Meng, Z. Li, D. Zhang, and G. E. Karniadakis.  PPINN: parareal physics-informed neural network for time-dependent PDEs.  Comput. Methods Appl. Mech. Engrg., 370:113250, 16, 2020.
  * [40] S. Mishra and R. Molinaro.  Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs.  IMA J. Numer. Anal., 42(2):981–1022, 2022.
  * [41] J. Müller and M. Zeinhofer.  Deep ritz revisited.  In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2019.
  * [42] J. Müller and M. Zeinhofer.  Error estimates for the deep ritz method with boundary penalty, 2021.
  * [43] B. Neyshabur, Z. Li, S. Bhojanapalli, Y. LeCun, and N. Srebro.  Towards understanding the role of over-parametrization in generalization of neural networks.  CoRR, abs/1805.12076, 2018.
  * [44] P. Petersen and F. Voigtlaender.  Optimal approximation of piecewise smooth functions using deep relu neural networks.  Neural Networks, 108:296–330, 2018.
  * [45] A. Pinkus.  Approximation theory of the MLP model in neural networks.  In Acta numerica, 1999, volume 8 of Acta Numer., pages 143–195. Cambridge Univ. Press, Cambridge, 1999.
  * [46] B. Qian, Y. Xiao, Z. Zheng, M. Zhou, W. Zhuang, S. Li, and Q. Ma.  Dynamic multi-scale convolutional neural network for time series classification.  IEEE Access, 8:109732–109746, 2020.
  * [47] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  J. Comput. Phys., 378:686–707, 2019.
  * [48] S. Shalev-Shwartz and S. Ben-David.  Understanding Machine Learning - From Theory to Algorithms. Cambridge University Press, 2014.
  * [49] J. Shen.  Efficient spectral-Galerkin method. I. Direct solvers of second- and fourth-order equations using Legendre polynomials.  SIAM J. Sci. Comput., 15(6):1489–1505, 1994.
  * [50] J. Shen, T. Tang, and L.-L. Wang.  Spectral methods, volume 41 of Springer Series in Computational Mathematics.  Springer, Heidelberg, 2011.  Algorithms, analysis and applications.
  * [51] Y. Shin, J. Darbon, and G. E. Karniadakis.  On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type PDEs.  Commun. Comput. Phys., 28(5):2042–2074, 2020.
  * [52] Y. Shin, Z. Zhang, and G. E. Karniadakis.  Error estimates of residual minimization using neural networks for linear pdes.  arXiv:2010.08019 [math.NA], 2020.
  * [53] K. Shukla, P. C. D. Leoni, J. L. Blackshire, D. Sparkman, and G. E. Karniadakis.  Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks.  CoRR, abs/2005.03596, 2020.
  * [54] J. W. Siegel and J. Xu.  Approximation rates for neural networks with general activation functions.  Neural Networks, 128:313–321, 2020.
  * [55] J. W. Siegel and J. Xu.  High-order approximation rates for shallow neural networks with cosine and reluk activation functions.  Applied and Computational Harmonic Analysis, 58:1–26, 2022.
  * [56] J. Sirignano and K. Spiliopoulos.  DGM: a deep learning algorithm for solving partial differential equations.  J. Comput. Phys., 375:1339–1364, 2018.
  * [57] M. J. Wainwright.  High-dimensional statistics, volume 48 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge, 2019.  A non-asymptotic viewpoint.
  * [58] S. Wang, S. Sankaran, and P. Perdikaris.  Respecting causality is all you need for training physics-informed neural networks, 2022.
  * [59] D. Yarotsky.  Error bounds for approximations with deep relu networks.  Neural Networks, 94:103–114, 2017.
  * [60] D. Yarotsky.  Optimal approximation of continuous functions by very deep relu networks.  In S. Bubeck, V. Perchet, and P. Rigollet, editors, Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pages 639–649. PMLR, 06–09 Jul 2018.
  * [61] J. Zhuang, D. Kochkov, Y. Bar-Sinai, M. P. Brenner, and S. Hoyer.  Learned discretizations for passive scalar advection in a two-dimensional turbulent flow.  Phys. Rev. Fluids, 6:064605, Jun 2021.