11institutetext: University of the Basque Country (UPV/EHU), Leioa, Spain
22institutetext: Basque Center for Applied Mathematics (BCAM), Bilbao, Spain
33institutetext: Ikerbasque (Basque Foundation for Sciences), Bilbao, Spain

# On quadrature rules for solving Partial Differential Equations using Neural
Networks

Jon A. Rivera  1122 Jamie M. Taylor  22 Ángel J. Omella  11 David Pardo
112233

###### Abstract

Neural Networks have been widely used to solve Partial Differential Equations.
These methods require to approximate definite integrals using quadrature
rules. Here, we illustrate via 1D numerical examples the quadrature problems
that may arise in these applications and propose different alternatives to
overcome them, namely: Monte Carlo methods, adaptive integration, polynomial
approximations of the Neural Network output, and the inclusion of
regularization terms in the loss. We also discuss the advantages and
limitations of each proposed alternative. We advocate the use of Monte Carlo
methods for high dimensions (above 3 or 4), and adaptive integration or
polynomial approximations for low dimensions (3 or below). The use of
regularization terms is a mathematically elegant alternative that is valid for
any spacial dimension, however, it requires certain regularity assumptions on
the solution and complex mathematical analysis when dealing with sophisticated
Neural Networks.

###### Keywords:

Deep learning Neural Networks Ritz method Least-Squares method Quadrature
rules.

##  1 Introduction

In the last years, the use of Deep Learning (DL) methods has grown
exponentially in multiple areas, including self-driving cars [36, 16], speech
recognition [2, 3], and healthcare [12, 33]. Similarly, the use of DL
algorithms has become popular for solving Partial Differential Equations
(PDEs) – see, e.g., [47, 6, 40, 39, 26, 42, 7, 32].

DL techniques present several advantages and limitations with respect to
traditional PDE solvers based on Finite Elements (FE) [18], Finite Differences
(FD) [24], or Isogeometric Analysis (IGA) [30]. Among the advantages of DL, we
encounter the nonnecessity of generating a grid. In general, DL uses a dataset
in which each datum is independent from others. In contrast, in the linear
system that produces the FEM method there exists a connectivity between the
nodes of the mesh. The independence of the data in DL allows the
parallelization into GPUs for fast computations. Furthermore, DL provides the
possibility of solving certain problems that cannot be solved via traditional
methods, like high-dimensional PDEs [17, 11], some fractional PDEs [4, 31],
and multiple nonlinear PDEs [35, 19].

However, DL also presents limitations when solving PDEs. For example, in [44]
they show that Fully-Connected Neural Networks suffer from spectral bias [34]
and exhibit different convergence rates for each loss component. In addition,
the convergence of the method is often assumed (see, e.g., [27]) since it
cannot be rigorously guaranteed due to the non-convexity of the loss function.
Another notorious problem is due to quadrature errors. In traditional mesh-
based methods, such as FEM, we first select an approximating solution space,
and then we compute the necessary integrals over each element required to
produce the stiffness matrix. With DL, we first set a quadrature rule and then
construct the approximated function. As we will show throughout this work,
this process may lead to large quadrature errors due to the unknown form of
the integrand. This can also be interpreted as a form of overfitting [48, 49]
that applies to the PDE constraint (e.g., $u"=0$) rather than to the PDE
solution (e.g., $u$). This may have disastrous consequences since it may
result in approximate solutions that are far away from the exact ones at all
points.

There exist different methods to approximate definite integrals by discrete
data. The most common ones used in DL are based on Monte Carlo integration.
These methods compute definite integrals using randomly sampled points from a
given distribution. Monte Carlo integration is suitable for high-dimensional
integrals. However, for low-dimensional integrals (1D, 2D, and 3D),
convergence is slow in terms of the number of integration points in comparison
to other quadrature rules. This produces elevated computational costs.
Examples of existing works that follow a Monte Carlo approach are [45] using a
Deep Ritz Method (DRM) and [43] using a Deep Galerkin Method (DGM). From the
practical point of view, at each iteration of the Stochastic Gradient Descent
(SGD), they consider a mini-batch of randomly selected points to discretize
the domain.

Another existing method to compute integrals in DNNs is the so-called
automatic integration [25]. In this method, the author approximates the
integrand by its high order Taylor series expansion around a given point
within the integration domain. Then, the integrals are computed analytically.
The derivatives needed in the Taylor series expansion are computed via
automatic differentiation (a.k.a autodiff) [15]. Since we are using the
information of the derivatives locally, we are doomed to perform an
overfitting on those conditions.

Another alternative is to use adaptive integration methods. In the Deep Least
Square (DLS) method [9], they use an adaptive mid-point quadrature rule using
local error indicators. These indicators are based on the value of the
residual at randomly selected points, and the resulting quadrature error is
unclear. In Variational Neural Networks (VarNets) [22], they use a Gauss
quadrature rule to evaluate their integrals. As error indicators, they use the
strong form of the residual evaluated over a set of random points that follow
a uniform distribution.

One can also use Gauss-type quadrature rules to evaluate integrals, as they do
in Variational Physics-Informed Neural Networks (VPINNs) [21]. One limitation
of this method is that we cannot select an adequate quadrature order because
we are unaware of the solution properties. In addition, as we are using fixed
quadrature points, the chances of performing overfitting are high.

This works analyzes the problems associated with quadrature rules in DL
methods when solving PDEs. As they remark in [21], there is no proper
quadrature rule in the literature developed for integrals of Deep Neural
Networks (DNNs). In particular, we illustrate the disastrous consequences that
may appear as a result of an inadequate selection of a quadrature rule. For
that, we use a simple one-dimensional (1D) Laplace problem with the Deep Ritz
Method (DRM) [45]. We also show a theoretical example of the Deep Least
Squares method. Then, we propose several alternatives to overcome the
quadrature problems: Monte Carlo integration, exact (Gaussian) integration of
piecewise-polynomial approximations, adaptive integration, and the use of
regularizers. We discuss their advantages and limitations and develop adequate
regularizers for certain problems using similar ideas to those presented in
[27]. Moreover, we illustrate with a numerical example the different
approximated solutions obtained with the different integration strategies.

The remainder of this article is as follows. Section 2 presents our selected
model problem. Section 3 defines the loss functions used in this work for
solving PDEs. Section 4 describes our Neural Network (NN) and explains some
critical implementation aspects. Section 5 illustrates the quadrature problems
via two numerical examples. Section 6 proposes several methods to overcome the
quadrature problems and Section 7 shows numerical results corresponding to the
proposed integration strategies. Finally, Section 8 summarizes the main
findings and possible future lines of research in the area.

##  2 Model problem

Let $\Omega\subset R^{d}$ be a computational domain, and $\Gamma_{D}$ and
$\Gamma_{N}$ two disjoint sections of its boundary, where
$\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$ and the subscripts $D$ and $N$
denote the Dirichlet and Neumann bounds, respectively. We consider the
following boundary value problem:

$\displaystyle[left=\empheqlbrace]-\nabla\cdot(\sigma\nabla u)$
$\displaystyle=f\qquad x\in\Omega,$ (1a) $\displaystyle u$
$\displaystyle=0\qquad x\in\Gamma_{D},$ (1b) $\displaystyle(\sigma\nabla
u)\cdot\bf{n}$ $\displaystyle=g\qquad x\in\Gamma_{N}.$ (1c)

In the above, $\bf{n}$ is the unit normal outward (to the domain) vector, and
we assume $\sigma(x)>0$ for every $x$ and the usual regularity assumptions,
namely, $\sigma\in L^{\infty}(\Omega)$, $f\in L^{2}(\Omega)$, $g\in
H^{1/2}(\Gamma_{N})$, and $u\in V=H^{1}_{0}(\Omega)=\\{v\in
H^{1}(\Omega)\textup{ and }v|_{\Gamma_{D}}=0\\}$, where $H^{1}(\Omega)=\\{v\in
L^{2}(\Omega),\nabla v\in L^{2}(\Omega)\\}$.

In this work we solve two different model problems to illustrate our numerical
results. The solution to model problem 1 is $u(x)=x^{0.7}$ and it satisfies
Eq. (2).

$\left\\{\begin{array}[]{rrcl}-u^{\prime\prime}(x)&=&0.21x^{-1.3}&\qquad
x\in(0,10),\\\ u(0)&=&0,\\\
u^{\prime}(10)&=&\frac{0.7}{10^{0.3}}.\end{array}\right.$ (2)

The solution to model problem 2 is $u(x)=x^{2}$ and it satisfies Eq. (3).

$\left\\{\begin{array}[]{rrcl}u^{\prime\prime}(x)&=&2&\qquad x\in(0,10),\\\
u(0)&=&0,\\\ u^{\prime}(10)&=&20.\end{array}\right.$ (3)

##  3 Loss functions

We introduce the following standard $L^{2}(\Omega)$ inner products:

$(u,v)=\int_{\Omega}u\>v\qquad\textup{and}\qquad(g,v)_{\Gamma_{N}}=\int_{\Gamma_{N}}g\>v.$
(4)

In the following, we consider two methods: the Ritz Method [37], and the Least
Squares Method [29].

###  3.1 Ritz Method

Multiplying the PDE from Eq. (1a) by a test function $v\in V$ (where
$V=H^{1}_{0}(\Omega)$), integrating by parts and incorporating the boundary
conditions, we arrive at the variational formulation:

$\textup{Find }u\in V\textup{ such that }(\sigma\nabla u,\nabla
v)=(f,v)+(g,v)_{\Gamma_{N}}\qquad\forall v\in V.$ (5)

It is easy to prove that the original (strong) and variational (weak)
formulations are equivalent (see, for instance, [20]).

To introduce the Ritz method, we define the energy function
$\mathcal{F}_{R}:V\longrightarrow\mathbb{R}$ given by

$\mathcal{F}_{R}(v)=\frac{1}{2}(\sigma\nabla v,\nabla
v)-(f,v)-(g,v)_{\Gamma_{N}}.$ (6)

As proved in [20], problem (5) is equivalent to the following energy
minimization problem:

$u=\underset{v\in V}{\textup{ arg min }}\mathcal{F}_{R}(v),$ (7)

###  3.2 Least Squares (LS) Method

Reordering the terms of Eq. (1a) and (1c), we define:

$\displaystyle{\cal G}u:=\nabla\cdot(\sigma\nabla u)+f\qquad x\in\Omega,$ (8a)
$\displaystyle{\cal B}u:=(\sigma\nabla u)\cdot\textbf{n}-g\qquad
x\in\Gamma_{N}.$ (8b)

To introduce the Least Squares method, we define the function
$\mathcal{F}_{LS}:V\longrightarrow\mathbb{R}$

$\mathcal{F}_{LS}(v)=|({\cal G}v,{\cal G}v)|+|({\cal B}v,{\cal
B}v)_{\Gamma_{N}}|.$ (9)

We want to minimize the function $\mathcal{F}_{LS}(v)$ subject to the
essential (Dirichlet) BCs. We often find the minimum by taking the derivative
equal to zero and ending up with a linear system of equations. In the context
of DL, we can simply introduce the above loss function $\mathcal{F}_{LS}(v)$
directly in our NN. Therefore, we want to find

$u=\underset{v\in V}{\textup{ arg min }}\mathcal{F}_{LS}(v).$ (10)

##  4 Neural Network Implementation

We train a Neural Network, named $u_{NN}(x;\theta)$, with the following
architecture. In this first part, we define the trainable part of our NN, with
learnable parameters $\theta$. We call it $u_{\theta}$. It is composed by:

  1. 1.

An input layer. This layer receives the data in the form of a $n\times d$
vector, where $n$ is the number of samples and $d$ is the dimension of the
data of the problem.

  2. 2.

$l$ hidden dense layers with $n_{i}$ neurons at each layer
$i\in\\{1,\cdots,l\\}$ and a sigmoid activation function. In our particular
examples, we select $l=1$.

  3. 3.

An output layer that delivers $u_{\theta}$. enumerate

Now, we add more layers to our scheme in order to implement Equations (6) or
(9). For that, we introduce:

enumerate

  4. 4.

A non-trainable layer to impose the Dirichlet boundary conditions. For that,
we select a function $\phi(x)$ that satisfies the Dirichlet conditions of the
problem and its value is nonzero everywhere else [23]. In this work, we select
the following $\phi(x)$ functions for 1D problems in the interval
$\Omega=[a,b]$:

$\phi(x)=\prod_{x_{D}\in\Gamma_{D}}(x-x_{D}).$ (11)

Then, we generate a new output of the NN:
$u_{NN}(x;\theta)=\phi(x)u_{\theta}(x)$ that strongly imposes the homogeneous
Dirichlet boundary conditions.

  5. 5.

A non-trainable layer to compute the loss function $\mathcal{F}_{R}$ or
$\mathcal{F}_{LS}$ following Eqs. (6) or (9). Within this layer, we evaluate
the integrals and the derivatives. We consider different quadrature rules,
being the quadrature points part of the input data of our NN, along with the
physical points of the domain. For computing the derivatives, we use automatic
differentiation, except in some specific cases, where we employ finite
differences. These cases are explicitly indicated throughout the text.

Figure 1 shows a schematic graph of the described NN architecture. Our
software is developed in Python and we use the library Tensorflow 2.0.

To train the NN, we replace in Equations (7) or (10) the search space $V$ by
the learnable parameters $\theta$ included in our NN. The result of the
minimization is a function $u_{NN}(x,\tilde{\theta})$, where $\tilde{\theta}$
are the optimal learnable parameters encountered as a result of the training.
For simplicity, in the following we abuse notation and use the symbol $u_{NN}$
to denote also the solution $u_{NN}(x,\tilde{\theta})$ of our minimization
problem.

Figure 1: Our implementation scheme.

##  5 Quadrature rules

We approximate our integrals from Eqs. (6) and (9) using a quadrature rule of
the form

$\int_{a}^{b}f(x)dx\approx\sum_{i=0}^{n}\omega_{i}f(x_{i}),$ (12)

where $\omega_{i}$ are the weights and $x_{i}$ are the quadrature points.
Examples of quadrature rules that follow the above formula include trapezoidal
rule and Gaussian quadrature rules [28]. We classify these quadrature rules
into two groups: (1) those that only employ points from the interior of the
interval; and (2) those that evaluate the solution at one extreme point or
more (a or b). Integration rules within the later group (e.g., the trapezoidal
rule) are inadequate for our minimization problems because the integrand can
be infinite at the boundary points in the case of singular solutions. Thus, we
focus on quadrature rules that only evaluate the solution at interior points
of the domain, with a special focus on Gaussian quadrature rules.

###  5.1 Illustration of quadrature problems in Neural Networks

####  5.1.1 Ritz method.

We consider the two model problems from Section 2. We approximate $u(x)$ using
the Ritz method. Thus, we search for a NN that minimizes the loss functional
given by Eq. (6). Our NN has one hidden layer with $10$ neurons ($31$
trainable weights). We use automatic differentiation to compute the
derivatives and a three-point Gaussian quadrature rule to approximate the
integrals within each element. We select the Stochastic Gradient Descent (SGD)
optimizer. For model problem 1, we discretize our domain with four equal-size
elements and execute $40,000$ iterations during the optimization process. For
model problem 2, we discretize our domain with ten equal-size elements and
execute $200,000$ iterations during the optimization process.

Figures 2(a) and 2(b) describe the loss evolution of the training process. We
obtain a lower loss than the best possible one (corresponding to the exact
solution), which is due to high quadrature errors.

(a) Model problem 1.

(b) Model problem 2.

Figure 2: Loss evolution of the training process for our two model problems.

Figures 3 and 4 compare the approximate and exact solutions. We observe a
disastrous NN approximation due to quadrature errors. Figures 3(b) and 4(b)
show that the gradient is (almost) zero at the training (quadrature) points of
the first interval. This value minimizes the numerical aproximation of the
integral of $(\nabla u)^{2}$. At the same time, the approximated solution is
large in the first interval, thus maximizing the term $(f,u)$, and minimizing
the total loss:

$\mathcal{F}_{R}(u_{NN})=\frac{1}{2}\underbrace{(\nabla u_{NN},\nabla
u_{NN})}_{\sum_{q_{i}}\omega_{i}(\nabla u_{NN})^{2}\approx
0}-\underbrace{(f,u_{NN})}_{\approx\infty}-(g,u_{NN})_{\Gamma_{N}}\simeq-\infty$

The described quadrature errors can be interpreted as overfitting over a
condition of the problem given by the energy of the solution.

(a) Exact and approximate solutions. (b) Approximate solution in the interval
$[0,2.5]$. The Gauss quadrature points corresponding to the first element are
indicated in blue.

Figure 3: Exact vs approximate Ritz method solutions of model problem 1 using
four elements for evaluation of $\mathcal{F}_{R}(v)$ and a NN with 31 weights.

(a) Exact and approximate solutions. (b) Approximate solution in the interval
$[0,1]$. The Gauss quadrature points corresponding to the first element are
indicated in blue.

Figure 4: Exact vs approximate Ritz method solutions of model problem 2 using
ten elements for evaluation of $\mathcal{F}_{R}(v)$ and a NN with 31 weights.

####  5.1.2 Least Squares Method.

We now consider the following one-dimensional problem:

$\left\\{\begin{array}[]{rrcl}-u^{\prime\prime}(x)&=&0&\qquad x\in(0,1),\\\
u(0)=u^{\prime}(1)&=&0,&\\\ \end{array}\right.$ (13)

where the exact solution is $u(x)=0$. We can easily construct an approximating
function $u_{NN}$ that satisfies Eq. (13) at the three considered Gaussian
points and minimizes Eq. (9), while still being a poor approximation of the
exact solution due to quadrature errors. Figure 5 shows an example.

Figure 5: Exact ($u_{exact}=0$) and approximated solution of problem given by
Eq. (13) and solved with the LS method.

##  6 Integral approximation

We now describe four different methods to improve the integral approximations.

###  6.1 Monte Carlo integration

We consider the following Monte Carlo integral approximation over a set of
points $x_{i}\in(a,b)$,

$\int_{a}^{b}f(x)dx\approx\frac{(b-a)}{n}\sum_{i=1}^{n}f(x_{i}),\hskip
14.22636ptx_{i}\in(a,b)\hskip 5.69046pt\forall i=\\{1,\cdots,n\\}$ (14)

In the above, points $x_{i}$ are randomly selected [1]. While this method is
useful for higher-dimensional integrals, for lower dimensions (1D, 2D, 3D) the
computational cost is high since the value of the integral approximation
converges as $1/\sqrt{n}$ [46].

###  6.2 Piecewise-polynomial approximation

We replace the original NN $u_{NN}$ by a piecewise-polynomial approximation
$u^{*}_{NN}$, which can be exactly differentiated (e.q., via finite
differences) and integrated (via a Gaussian quadrature rule). Figure 6 shows
an example of the aforementioned case when we train a NN using four elements
with linear approximations within each element.

This method controls quadrature errors. However, it is inadequate for high-
dimensional problems as we need a mesh that is difficult to implement and
integration becomes time consuming.

Figure 6: Neural Network approximation $u_{NN}$ and its piecewise-linear
approximation $u^{*}_{NN}$.

###  6.3 Adaptive integration

We first consider a training dataset over the interval $(a,b)$, and a
validation dataset built from a global h-refinement of the training dataset.
Figure 7 shows an example of a training and the corresponding validation
datasets. Then, for each element of the training mesh (e.g., $E_{1}$ in Figure
7), we compare the numerical integral over that element vs the sum of the
integrals over the two corresponding elements on the validation dataset (in
our case, $E_{1}^{1}+E_{1}^{2}$). If the integral values differ by more than a
stipulated tolerance, we h-refine the training element, and we upgrade the
validation dataset so it is built as a global h-refinement of the training
dataset. This process is described in Algorithm 1. Figure 7 also shows a
training and the corresponding validation dataset after refining the first and
third elements.

Generate a training dataset;

Generate the corresponding validation dataset;

Set tolerance $\epsilon$ and maximum iteration number $i_{max}$;

while  _$i <i_{max}$_ do

for  _$j=1,\cdots,n\textup{ (number of elements)}$_ do

Compute integral values $I_{j}$ over the training dataset of elements $E_{j}$;

Compute integral values $I_{j}^{1}$,$I_{j}^{2}$ over the validation dataset of
elements $E_{j}^{1}$,$E_{j}^{2}$ ;

if  _$(I_{j}^{1}+I_{j}^{2})-I_{j} <\epsilon$_ then

$h$-refine the $E_{j}$-th element of the training set;

$h$-refine the $E_{j}^{1}$-th and $E_{j}^{2}$-th elements of the validation
set;

else

continue;

end if

end for

$i=i+1$;

end while

Algorithm 1 Adaptive integration method

(a) Training set (b) Validation set

Figure 7: Black points (dots) correspond to the original (a) training/(b)
validation partitions and blue points (circles) are the points added by the
refinement performed in the first and third elements.

We are able to control the quadrature errors by adding new quadrature points
to the training dataset. However, the simplest way to implement such method is
by using meshes, which posses a limitation on high-dimensional integrals. As
an alternative to generating a mesh, one can randomly add points to the
training set. This entails difficulties when designing an adaptive algorithm.

In the same way that we propose an h-adaptive method, we can also work with
p-adaptivity [5] or a combination of them (e.g., hp-adaptivity [10]).

###  6.4 Regularization methods

We now introduce a problem-specific regularizer designed to control the
quadrature error.

In a one-dimensional setting, we consider the integral functional
$\mathcal{F}_{R}$ as given in (6), and its approximation via a midpoint rule,
$\hat{\mathcal{F}}_{R}$, given by

$\begin{split}\hat{\mathcal{F}}_{R}(u)=&\frac{b-a}{N}\sum\limits_{j=1}^{N}\left(\frac{1}{2}|u^{\prime}(x_{j})|^{2}-f(x_{j})u(x_{j})\right)+g(b)u(b)+g(a)u(a).\end{split}$
(15)

We note that $g=0$, except where the Neumann condition is imposed. While we
focus on the Ritz method, a similar heuristic can be applied to the Least
Squares method.

We introduce a function $\mathcal{R}$ that depends on the learnable parameters
$\theta$ of a given neural network $u_{NN}$, such that for any neural network
with a given architecture,

$|\mathcal{F}_{R}(u_{NN})-\hat{\mathcal{F}}_{R}(u_{NN})|<\mathcal{R}(\theta).$
(16)

If we then consider a loss function $\mathcal{L}$ given by

$\mathcal{L}(\theta)=\hat{\mathcal{F}}_{R}(u_{NN})+\mathcal{R}(\theta),$ (17)

then we may be able to improve the approximation of the quadrature rule, as
the loss contains a term that by design controls the quadrature error.

For simplicity, we consider only the case of a single-layer network with a
one-dimensional input, and the midpoint rule for calculating the integral over
a uniform partition of $(a,b)$. We consider a mid-point rule as in (15), and
define the interval length $\delta=\frac{b-a}{N}$ and intervals
$I_{j}=\left(\frac{\delta}{2}+x_{j},x_{j}+\frac{\delta}{2}\right)$. We
estimate the error of the midpoint rule to integrate $F$ as

$\begin{split}\left|\int_{a}^{b}F(x)\,dx-\sum\limits_{j=1}^{N}F(x_{j})\delta\right|=&\left|\sum\limits_{i=1}^{N}\int_{x_{j}-\frac{\delta}{2}}^{x_{j}+\frac{\delta}{2}}(F(x)-F(x_{j}))\,dx\right|\\\
\leq&\sum\limits_{i=1}^{N}\int_{x_{j}-\frac{\delta}{2}}^{x_{j}+\frac{\delta}{2}}\left|F(x)-F(x_{j})\right|\,dx\\\
\leq&\sum\limits_{i=1}^{N}\max\limits_{t\in
I_{j}}|F^{\prime}(t)|\int_{x_{j}-\frac{\delta}{2}}^{x_{j}+\frac{\delta}{2}}|x-x_{j}|\,dx\\\
=&\frac{\delta^{2}}{4}\sum\limits_{i=1}^{N}\max\limits_{t\in
I_{j}}|F^{\prime}(t)|.\end{split}$ (18)

This estimate scales as $\mathcal{O}(\frac{1}{N})$ for fixed $F$, and thus,
for a large number of integration points, we expect the estimate to be
sufficiently accurate and to avoid “overdamping” of the loss.

With (18) in mind, we estimate the local Lipschitz constants of the integrand
as in (6). The numerical estimation of the Lipschitz constants of NNs has
attracted attention, as they form a way of of estimating the generalizability
of a neural network, and have been used in the training process as a way to
encourage accurate generalization [13, 14, 41]. As we are dealing with loss
functions that involve derivatives of the neural network, we however need
estimates of higher order derivatives of $u_{NN}$. The approach that we employ
is similar in spirit to the work of [27] for obtaining a posteriori error
estimates in PINNs.

Despite the arithmetic complications involved in calculating $\mathcal{R}$,
conceptually the idea reduces to an application of Taylor’s theorem. On a
single interval of integration $I_{j}$, we have that for every $x$, there
exists some $\xi_{x}$ so that

$|F^{\prime}(x)|=|F^{\prime}(x_{i})+(x-x_{i})F^{\prime\prime}(\xi_{x})|\leq|F^{\prime}(x_{i})|+\frac{\delta}{2}||F^{\prime\prime}||_{\infty}.$
(19)

We then find $\mathcal{R}$ using a combination of local and global estimates
for the derivatives of the integrand corresponding to simple pointwise
evaluations at the integration points and global estimates involving the
neural network weights. The exact calculation of the estimates that produce
the regularizer $\mathcal{R}$ are tedious and deferred to Appendix 0.A, with
$\mathcal{R}$ given in (37) via expressions in (26) and (30).

##  7 Numerical Results

In this section, we test some of the alternatives proposed in Section 6 to
overcome the quadrature problems. Specifically, we solve the two model
problems from Section 2 using: (a) a piecewise-linear approximation of the NN,
and (b) an adaptive integration method. We also solve model problem 2 using
regularization methods.

For the cases of piecewise-linear approximation and adaptive integration, we
use a NN architecture composed of one hidden layer with ten neurons and a
sigmoid activation function. Moreover, we select SGD as the optimizer. In the
case of regularization, we use a hyperbolic tangent activation function with
the Adam optimizer.

###  7.1 Piecewise-linear approximation

We select a piecewise-linear approximation of the NN as our approximate
solution. We compute the gradients using Finite Differences and the integrals
using a one-point Gaussian quadrature rule (i.e., the midpoint rule). For
model problem 1, we use two different uniform partitions composed of four and
ten elements, and execute $40,000$ iterations. For model problem 2, we use a
uniform partition composed of ten elements and execute $200,000$ iterations.

Figures 8(a) and 8(b) show that the loss converges to the loss of the exact
solution. We also observe better results as we increase the number of
elements, as physically expected. Figures 9(a) and 9(b) show the corresponding
solutions, which are consistent with the loss evolutions displayed in Figures
8(a) and 8(b).

(a) Model problem 1. (b) Model problem 2.

Figure 8: Loss evolution of the training process for our two model problems
when we use a piecewise-linear approximation of the NN.

(a) Model problem 1. (b) Model problem 2.

Figure 9: Ritz method solution when we use a piecewise-linear approximation of
the NN to solve the problem.

While the use of a piecewise-linear approximation overcomes the quadrature
problems, the convergence is limited to $\mathcal{O}(h)$, where $h$ is the
element size [38]. To increase this speed, it is possible to consider
different piecewise-polynomial approximations, including the use of r-adaptive
algorithms [8].

###  7.2 Adaptive integration

We compute the gradients using automatic differentiation and the integrals
using a three-point Gaussian quadrature rule. As in previous examples, we
start the training with a uniform partition of four elements for model problem
1 and of ten elements for model problem 2. We select a half-size partition of
the training dataset for validation. For model problem 1, we compare the
integral values of the training and validation sets (i.e., we execute
Algorithm 1) every $100$ iterations, with an error tolerance of $0,0001$; for
model problem 2, we compare every $10,000$ iterations, with an error tolerance
of $10$. The tolerance is selected as a small percentage of the value of the
loss.

Figures 10(a) and 10(b) show that the loss converges to the optimum value. As
explained in Section 6, the adaptive integration algorithm refines the
training dataset. For model problem 1, the algorithm refines the first
interval at iteration $25,000$. For model problem 2, two refinements occur in
the first interval. The adaptive integration algorithm automatically selects
the first interval for refinement, where overfitting was taking place in
Figures 3 and 4. Figures 11(a) and 11(b) show the corresponding solutions. We
observe that the approximate solutions properly approximate the exact ones.

(a) Model problem 1. (b) Model problem 2.

Figure 10: Loss evolution of the training process for our two model problems
when using adaptive integration.

(a) Model problem 1. (b) Model problem 2.

Figure 11: Ritz method solution when using adaptive integration.

The extrapolation of this method to higher dimensions (2D or 3D) requires
higher-dimensional discretizations and quadrature rules.

###  7.3 Regularization methods

We do not apply the regularization method to model problem 1 as the method
requires sufficient regularity in order to provide the necessary estimates in
the calculation of $\mathcal{R}$. Since the solution is singular at $x=0$, the
necessary Lipschitz bounds on the integral functional cannot be obtained
within this framework. Instead, we aim to demonstrate that for problems that
are sufficiently regular, our technique can avoid overfitting, and leave open
the question as to how one may adapt the technique to singular problems for
future work. We thus consider model problem 2. We propose the loss defined via

$\mathcal{L}(\theta)=\mathcal{\hat{F}}_{R}(u_{NN})+\mathcal{R}(\theta).$ (20)

Explicitly,

$\mathcal{\hat{F}}_{R}(u_{NN})=\frac{10}{N}\sum\limits_{j=1}^{N}\frac{1}{2}|u^{\prime}_{NN}(x_{j})|^{2}-2u_{NN}(x_{j})-20u_{NN}(20),$
(21)

where $x_{j}=\frac{10}{N}\left(i-\frac{1}{2}\right)$.

####  7.3.1 Experiment 1

We consider $N=50$ points, and a single layer network with $M=10$ neurons. We
use the Adam optimizer with learning rate $10^{-2}$. We solve model problem 2
with two losses: with and without regularization. In both cases, we measure
the metrics $\mathcal{L}$, $\mathcal{R}$, and $\mathcal{\hat{F}}_{R}$. For
validation, we use an equidistant partition of $(0,10)$ with $49$ points, so
that we still use a midpoint rule but with different integration points.

(a) Exact and approximated solutions.

(b) Evolution of $\mathcal{L}$ during training.

(c) Evolution of $\mathcal{\hat{F}}_{R}$ during training.

(d) Evolution of $\mathcal{R}$ during training.

Figure 12: The solution and training information for Experiment 1 without
regularization.

(a) Exact and approximated solutions.

(b) Evolution of $\mathcal{L}$ during training.

(c) Evolution of $\mathcal{\hat{F}}_{R}$ during training.

(d) Evolution of $\mathcal{R}$ during training.

Figure 13: The solution and training information for Experiment 1 with
regularization.

Figure 12 shows the results without regularization. As expected, we see in
Figure 12(a) that the approximation is poor due to overfitting, which is most
notable around $x=0$ and attained within 5000 epochs. Via the provided plots
we can observe the beginning of overfitting in two distinct manners. First, we
observe in Figure 12(c) that the value of $\mathcal{\hat{F}}_{R}$ evaluated
over the validation data begins to diverge from the value on the training
data, becoming apparent at around 1000 epochs. We also see this behaviour
reflected in the evolution of $\mathcal{R}$ in Figure 12(d), with its most
dramatic increase beginning around the same iteration. This rapid increase
also provokes an increase in $\mathcal{L}$, as seen in Figure 12(b). This also
indicates that even if $\mathcal{R}$ is not used as part of the training
process, its increase could be used as a metric to identify overfitting.

Figure 13 describes the results with regularization and we observe a different
behaviour. The approximation is generally good, and we do not see any signs of
overfitting within $10^{5}$ epochs, as shown in Figure 13(a). In particular,
the values of $\mathcal{\hat{F}}_{R}$ at the training and validation data
remain consistent in Figure 13(c). Throughout Figure 13 we see that within
$10^{5}$ epochs all metrics appear to have converged to a limiting value. We
obtain final values $\mathcal{L}\approx-644.22$,
$\mathcal{\hat{F}}_{R}\approx-666.07$, $\mathcal{R}\approx 24.8$. We recall
that the true energy of the exact solution is
$\mathcal{F}_{R}(u_{exact})\approx-666.667$, which suggests the quadrature
rule is accurate. Notice that in the case without regularization, before
overfitting became apparent, $\mathcal{R}$ had already attained values of
around 1000, which is far larger than the value of $\mathcal{R}$ at the
obtained solution when regularization was used.

####  7.3.2 Experiment 2

We now consider a smaller $N$. As we expect $\mathcal{R}$ to scale as
$\frac{1}{N}$, we anticipate a more adverse effect when $N$ is small. To view
this, we consider the same problem of Experiment 1, where we now select $N=20$
integration points. We consider $M=10$ neurons and minimize our problem using
the Adam optimizer with a learning rate of $10^{-2}$. As before, we consider
the cases with and without regularization.

(a) Exact and approximated solutions.

(b) Evolution of $\mathcal{L}$ during training.

(c) Evolution of $\mathcal{\hat{F}}_{R}$ during training.

(d) Evolution of $\mathcal{R}$ during training.

Figure 14: The solution and training information for Experiment 2 without
regularization.

(a) Exact and approximated solutions.

(b) Evolution of $\mathcal{L}$ during training.

(c) Evolution of $\mathcal{\hat{F}}_{R}$ during training.

(d) Evolution of $\mathcal{R}$ during training.

Figure 15: The solution and training information for Experiment 2 with
regularization.

Figure 14 presents the loss evolution without regularization. We observe
overfitting, which is accompanied by divergence of the loss on the validation
dataset, as well as a rapid increase in $\mathcal{R}$, with these features
visible within 5000 epochs.

Figure 15 presents the results with regularization. We observe no signs of
overfitting, with the validation and training loss remaining close in Figure
15(b). All metrics appear to have converged to a limiting value within
$10^{4}$ epochs. However, the large value of $\mathcal{R}$ at the found
solution (approximately 140) has substantially changed the optimization
problem so that the obtained minimizer is far from the desired solution. The
final value of $\mathcal{\hat{F}}_{R}$ is around $-622$, which is far from the
desired value of $-666.67$. This experiment highlights the fact that the
regularizer becomes more effective when a large number of integration points
are used.

##  8 Conclusions & Future work

We first illustrated how quadrature errors can destroy the quality of the
approximated solution when solving PDEs using DL methods. Thus, it is crucial
to select an adequate method to overcome the quadrature problems. Herein, we
proposed four different alternatives: (a) Monte Carlo integration, (b)
piecewise-polynomial approximations of the output of the Neural Network, (c)
adaptive integration, and (d) regularization methods. We discussed the
advantages and limitations of each of these methods, and we illustrated their
performance via simple 1D numerical examples.

In high dimensions, Monte Carlo integration methods are the best choice.
Regularizer methods are another option, but they are problem dependent and
they need to be derived for each different architecture. Moreover, they
require further analysis for highly nonlinear integrands. Furthermore, they
are limited only to sufficiently smooth integral functionals. In addition,
more complex NN architectures (which should be needed in higher dimension)
will hinder the derivation of $\mathcal{R}$.

In low dimensions (three or below), Monte Carlo integration is not competitive
because of its low convergence speed. In these cases, adaptive integration
exhibits faster convergence. In the cases of piecewise-linear approximation
and regularizers, we are also able to overcome the quadrature problems, but
the convergence speed is often slower than with adaptive integration.

Possible future research lines of this work are: (a) to implement adaptive
integration method in 2D and 3D, (b) to apply r-adaptivity methods with
piecewise-polynomial approximations, and (c) to implement regularizers for
high dimensional problems involving more complex NN architectures.

##  9 Acknowledgments

This work has received funding from: the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No 777778 (MATHROCKS); the European Regional Development Fund (ERDF)
through the Interreg V-A Spain-France-Andorra program POCTEFA 2014-2020
Project PIXIL (EFA362/19); the Spanish Ministry of Science and Innovation
projects with references PID2019-108111RB-I00 (FEDER/AEI), PDC
2021-121093-I00, and PID2020-114189RB-I00 and the “BCAM Severo Ochoa”
accrediation of excellence (SEV-2017-0718); and the Basque Government through
the three Elkartek projects 3KIA (KK-2020/00049), EXPERTIA (KK-2021/000 48),
and SIGZE (KK-2021/00095), the Consolidated Research Group MATHMODE
(IT1294-19) given by the Department of Education, and the BERC 2018-2021
program.

##  Appendix 0.A Estimation of the regularizer

Following the heuristic of Section 6.4, we derive an expression for
$\mathcal{R}$ that may be used as a regularizer. The necessary steps are:

  1. 1.

Using the chain rule, we find global upper bounds for the derivatives of a
simple neural network in terms of the weights.

  2. 2.

Via Taylor’s theorem with remainder and using the global estimates for simple
networks, we find local estimates of the derivatives of NNs with a cutoff
function to ensure a homogeneous Dirichlet condition.

  3. 3.

Using local estimates for the derivatives of a NN, we find local Lipschitz
estimates for integrands corresponding to the Ritz method.

We tackle each of these estimations in the following subsections.

###  0.A.1 Global estimates for derivatives of a single layer network

Let $\hat{u}_{NN}$ be a single layer neural network. We write it in the form

$\hat{u}_{NN}(x)=b^{1}+\sum\limits_{i=1}^{M}A_{i}^{1}\sigma(A_{i}^{0}x+b^{0}_{i}).$
(22)

for weights $\theta=(b^{1},A^{1},b^{0},A^{0})$ and an activation function
$\sigma$. We assume that $\sigma$ has globally bounded derivatives, so that
$||\sigma^{(n)}||_{\infty}$ is finite for every $n=0,1,2...$.

A simple application of the triangle quality and that $\sigma$ is bounded
gives that

$\begin{split}|\hat{u}_{NN}(x)|\leq&|b^{1}|+\sum\limits_{i=1}^{M}|A_{i}^{1}\sigma(A_{i}^{0}x+b^{0}_{i})|\\\
\leq&|b^{1}|+\sum\limits_{i=1}^{M}|A_{i}^{1}|||\sigma||_{\infty}\\\
\end{split}$ (23)

The derivatives of $\hat{u}_{NN}$ are

$\hat{u}^{(n)}_{NN}(x)=\sum\limits_{i=1}^{M}A_{i}^{1}(A_{i}^{0})^{n}\sigma^{(n)}(A_{i}^{0}x+b^{0}_{i}),$
(24)

for $n\geq 1$. Similarly, this gives the immediate estimate that

$\left|\hat{u}^{(n)}_{NN}(x)\right|\leq\sum\limits_{i=1}^{M}|A_{i}^{1}||A_{i}^{0}|^{n}||\sigma^{(n)}||_{\infty}.$
(25)

Thus, we define the global upper bounding function $\mathcal{R}^{1}(\theta;n)$
by

$\mathcal{R}^{1}(\theta;n)=\left\\{\begin{array}[]{c
c}|b^{1}|+||\sigma||_{\infty}\sum\limits_{i=1}^{M}|A_{i}^{1}|&n=0,\\\
||\sigma^{(n)}||_{\infty}\sum\limits_{i=1}^{M}|A_{i}^{1}||A_{i}^{0}|^{n}&n\geq
1,\end{array}\right.$ (26)

which gives the global upper bound

$\left|\hat{u}_{NN}^{(n)}(x)\right|\leq\mathcal{R}^{1}(\theta;n)$ (27)

for every $x$.

###  0.A.2 Local derivative estimation of a single layer neural network with
cutoff function

Next, we consider the commonly considered case where our neural network admits
a cutoff function to ensure a Dirichlet boundary condition, and turn to the
case of local estimation. Explicitly, we take
$u_{NN}(x)=\hat{u}_{NN}(x)\phi(x)$, where $\phi$ is zero at the homogeneous
Dirichlet condition, and $\hat{u}_{NN}$ is as in (22). We presume that $\phi$
has bounded derivatives; i.e. $||\phi^{(k)}||_{\infty}$ is finite for
$k=0,1,...$. Let $x$ be in the interval
$I_{j}=\left(x_{j}-\frac{\delta}{2},x_{j}+\frac{\delta}{2}\right)$. Then, we
have

$\begin{split}u^{(n)}_{NN}(x)=&u^{(n)}_{NN}(x_{j})+(x-x_{j})u^{(n+1)}_{NN}(\xi_{x})\\\
=&u^{(n)}_{NN}(x_{j})+(x-x_{j})\sum\limits_{k=0}^{n+1}\binom{n+1}{k}u^{(k)}_{NN}(\xi_{x})\phi^{(n+1-k)}(\xi_{x})\end{split}$
(28)

via Taylor’s theorem with remainder and the product rule for higher order
derivatives. Employing the global upper bounds $\mathcal{R}^{1}(\theta;n)$
from (26), we have that for $x\in I_{j}$,

$\begin{split}\left|u^{(n)}_{NN}(x)\right|\leq&\left|u^{(n)}_{NN}(x_{j})\right|+\frac{\delta}{2}\sum\limits_{k=0}^{n+1}\binom{n+1}{k}|u^{(k)}_{NN}(\xi_{x})||\phi^{(n+1-k)}(\xi_{x})|\\\
\leq&\left|u^{(n)}_{NN}(x_{j})\right|+\frac{\delta}{2}\sum\limits_{k=0}^{n+1}\binom{n+1}{k}\mathcal{R}^{1}(\theta;k)||\phi^{(n+1-k)}||_{\infty}.\end{split}$
(29)

Thus, we define the second intermediate regularizer as

$\mathcal{R}^{2}(\theta;I_{j},n)=\left|u^{(n)}_{NN}(x_{j})\right|+\frac{\delta}{2}\sum\limits_{k=0}^{n+1}\binom{n+1}{k}\mathcal{R}^{1}(\theta;k)||\phi^{(n+1-k)}||_{\infty},$
(30)

giving the estimate that for all $x\in I_{j}$,

$\left|u^{(n)}_{NN}(x)\right|\leq\mathcal{R}^{2}(\theta;I_{j},n).$ (31)

We note that via automatic differentiation, $u^{(n)}$ may be evaluated at the
training data $x_{j}$.

###  0.A.3 Application to integral errors

Our aim is to estimate the error in the integral functional

$\mathcal{F}(u_{NN})=\int_{a}^{b}\frac{1}{2}|u_{NN}^{\prime}(x)|^{2}-f(x)u_{NN}(x)\,dx-g(a)u_{NN}(a)-g(b)u_{NN}(b),$
(32)

when approximated by a simple quadrature rule

$\sum\limits_{j=1}^{N}\left(\frac{1}{2}|u_{NN}^{\prime}(x_{j})|^{2}-f(x_{j})u_{NN}(x_{j})\right)\delta-g(b)u_{NN}(b)-g(a)u_{NN}(a).$
(33)

The boundary terms can be calculated in one dimension without quadrature error
and thus we ignore their contribution. We estimate the error for the
quadrature rule by obtaining Lipschitz bounds of the integrand via

$\begin{split}&\left|\frac{d}{dx}\left(\frac{1}{2}|u_{NN}^{\prime}(x)|^{2}-f(x)u_{NN}(x)\right)\right|\\\
=&\left|u^{\prime}_{NN}(x)u^{\prime\prime}_{NN}(x)-f^{\prime}(x)u_{NN}(x)-f(x)u^{\prime}_{NN}(x)\right|\\\
\leq&|u^{\prime}_{NN}(x)||u^{\prime\prime}_{NN}(x)|+|f^{\prime}(x)||u_{NN}(x)|+|f(x)||u^{\prime}_{NN}(x)|.\end{split}$
(34)

Estimating the (local) Lipschitz constant of the integrand reduces to
estimating (locally) various derivatives of $u_{NN}$. For $x\in I_{j}$, we
estimate the Lipschitz constant of the integrand via

$\left|\frac{d}{dx}\left(\frac{1}{2}|u_{NN}^{\prime}(x)|^{2}-f(x)u_{NN}(x)\right)\right|\leq\mathcal{R}^{3}(\theta;I_{j}),$
(35)

where the regularizer $\mathcal{R}^{3}(\theta,I_{j})$ is given by

$\mathcal{R}^{3}(\theta;I_{j})=\Big{(}\mathcal{R}^{2}(\theta;I_{j},1)\mathcal{R}^{2}(\theta;I_{j},2)+||f||_{\infty}\mathcal{R}^{2}(\theta;I_{j},1)+||f^{\prime}||_{\infty}\mathcal{R}^{2}(\theta;I_{j},0)\Big{)}.$
(36)

We define the final regularizer $\mathcal{R}$ by

$\mathcal{R}(\theta)=\frac{\delta^{2}}{4}\sum\limits_{j=1}^{N}\mathcal{R}^{3}(\theta;I_{j}),$
(37)

which following (18) gives the estimate

$\begin{split}&|\mathcal{F}(u_{NN})-\hat{\mathcal{F}}(u_{NN})|\\\
=&\left|\int_{a}^{b}\frac{1}{2}|u_{NN}^{\prime}(x)|^{2}-f(x)u_{NN}(x)\,dx-\sum\limits_{j=1}^{N}\left(\frac{1}{2}|u_{NN}^{\prime}(x_{j})|^{2}-f(x_{j})u_{NN}(x_{j})\right)\delta\right|\\\
\leq&\mathcal{R}(\theta).\end{split}$ (38)

## References

  * [1] 14 - Monte Carlo integration I: Basic concepts. In: Pharr, M., Humphreys, G. (eds.) Physically Based Rendering, pp. 631–660. Morgan Kaufmann, Burlington (2004). https://doi.org/10.1016/B978-012553180-1/50016-8
  * [2] Afouras, T., Chung, J.S., Senior, A., Vinyals, O., Zisserman, A.: Deep Audio-visual Speech Recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence PP, 1–1 (12 2018). https://doi.org/10.1109/TPAMI.2018.2889052
  * [3] Alam, M., Samad, M., Vidyaratne, L., Glandon, A., Iftekharuddin, K.: Survey on Deep Neural Networks in Speech and Vision Systems. Neurocomputing 417, 302–321 (2020). https://doi.org/10.1016/j.neucom.2020.07.053
  * [4] Antil, H., Khatri, R., Löhner, R., Verma, D.: Fractional deep neural network via constrained optimization. Machine Learning: Science and Technology 2(1) (dec 2020). https://doi.org/10.1088/2632-2153/aba8e7
  * [5] Barros, F.B., Proença, S.P.B., de Barcellos, C.S.: On error estimator and p-adaptivity in the generalized finite element method. International Journal for Numerical Methods in Engineering 60(14), 2373–2398 (2004). https://doi.org/https://doi.org/10.1002/nme.1048
  * [6] Berg, J., Nyström, K.: A unified deep artificial neural network approach to partial differential equations in complex geometries. Neurocomputing 317, 28–41 (2018). https://doi.org/10.1016/j.neucom.2018.06.056
  * [7] Brevis, I., Muga, I., Zee, K.: A machine-learning minimal-residual (ml-mres) framework for goal-oriented finite element discretizations. Computers & Mathematics with Applications 95 (2020). https://doi.org/10.1016/j.camwa.2020.08.012
  * [8] Budd, C.J., Huang, W., Russell, R.D.: Adaptivity with moving grids. Acta Numerica 18, 111–241 (2009). https://doi.org/10.1017/S0962492906400015
  * [9] Cai, Z., Chen, J., Liu, M., Liu, X.: Deep least-squares methods: An unsupervised learning-based numerical method for solving elliptic PDEs. Journal of Computational Physics 420, 109707 (2020). https://doi.org/10.1016/j.jcp.2020.109707
  * [10] Demkowicz, L., Rachowicz, W., Devloo, P.: A fully automatic hp-adaptivity. Journal of Scientific Computing 17 (2002). https://doi.org/10.1023/A:1015192312705
  * [11] Ee, W., Han, J., Jentzen, A.: Deep Learning-Based Numerical Methods for High-Dimensional Parabolic Partial Differential Equations and Backward Stochastic Differential Equations. To appear in Communications in Mathematics and Statistics 5 (06 2017). https://doi.org/10.1007/s40304-017-0117-6
  * [12] Esteva, A., Robicquet, A., Ramsundar, B., Kuleshov, V., DePristo, M., Chou, K., Cui, C., Corrado, G., Thrun, S., Dean, J.: A guide to deep learning in healthcare. Nature Medicine 25, 24–29 (2019). https://doi.org/10.1038/s41591-018-0316-z
  * [13] Fazlyab, M., Robey, A., Hassani, H., Morari, M., Pappas, G.: Efficient and accurate estimation of Lipschitz constants for deep neural networks. Advances in Neural Information Processing Systems 32, 11427–11438 (2019)
  * [14] Gouk, H., Frank, E., Pfahringer, B., Cree, M.J.: Regularisation of neural networks by enforcing lipschitz continuity. Machine Learning 110(2), 393–416 (2021)
  * [15] Güneş Baydin, A., Pearlmutter, B.A., Andreyevich Radul, A., Mark Siskind, J.: Automatic differentiation in machine learning: A survey. Journal of Machine Learning Research 18, 1–43 (2018)
  * [16] Gupta, A., Anpalagan, A., Guan, L., Khwaja, A.S.: Deep learning for object detection and scene perception in self-driving cars: Survey, challenges, and open issues. Array 10, 100057 (2021). https://doi.org/10.1016/j.array.2021.100057
  * [17] Han, J., Jentzen, A., Ee, W.: Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences 115 (07 2017). https://doi.org/10.1073/pnas.1718942115
  * [18] Hauser, J.R. (ed.): Partial Differential Equations: The Finite Element Method, pp. 883–987. Springer Netherlands, Dordrecht (2009). https://doi.org/10.1007/978-1-4020-9920-5_13
  * [19] Huré, C., Pham, H., Warin, X.: Some machine learning schemes for high-dimensional nonlinear PDEs. Math. Comput. 89, 1547–1579 (2020)
  * [20] Johnson, C.: Numerical Solution of Partial Differential Equations by the Finite Element Method. Dover Books on Mathematics Series, Dover Publications, Incorporated (2012)
  * [21] Kharazmi, E., Zhang, Z., Karniadakis, G.E.: VPINNs: Variational Physics-Informed Neural Networks For Solving Partial Differential Equations (2019)
  * [22] Khodayi-Mehr, R., Zavlanos, M.: VarNet: Variational Neural Networks for the Solution of Partial Differential Equations. In: Bayen, A.M., Jadbabaie, A., Pappas, G., Parrilo, P.A., Recht, B., Tomlin, C., Zeilinger, M. (eds.) Proceedings of the 2nd Conference on Learning for Dynamics and Control. Proceedings of Machine Learning Research, vol. 120, pp. 298–307. PMLR, The Cloud (10–11 Jun 2020)
  * [23] Lagaris, I., Likas, A., Fotiadis, D.: Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks 9, 987–1000 (1998). https://doi.org/10.1109/72.712178
  * [24] LeVeque, R.: Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems (Classics in Applied Mathematics Classics in Applied Mathemat). Society for Industrial and Applied Mathematics, USA (2007)
  * [25] Liu, K.: Automatic Integration (2020), https://arxiv.org/abs/2006.15210
  * [26] Lu, L., Meng, X., Mao, Z., Karniadakis, G.E.: DeepXDE: A deep learning library for solving differential equations (2020), https://arxiv.org/abs/1907.04502
  * [27] Mishra, S., Molinaro, R.: Estimates on the generalization error of physics informed neural networks (PINNs) for approximating PDEs. arXiv preprint arXiv:2006.16144 (2020)
  * [28] Moin, P.: Fundamentals of Engineering Numerical Analysis. Cambridge University Press, 2 edn. (2010). https://doi.org/10.1017/CBO9780511781438
  * [29] Mortari, D.: Least-Squares Solution of Linear Differential Equations. Mathematics 5 (02 2017). https://doi.org/10.3390/math5040048
  * [30] Nguyen, V.P., Anitescu, C., Bordas, S.P., Rabczuk, T.: Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation 117, 89–116 (2015). https://doi.org/10.1016/j.matcom.2015.05.008
  * [31] Pang, G., Lu, L., Karniadakis, G.E.: fPINNs: Fractional Physics-Informed Neural Networks. SIAM Journal on Scientific Computing 41(4), A2603–A2626 (2019). https://doi.org/10.1137/18M1229845
  * [32] Paszyński, M., Grzeszczuk, R., Pardo, D., Demkowicz, L.: Deep learning driven self-adaptive hp finite element method. In: Computational Science – ICCS 2021. pp. 114–121. Springer International Publishing (2021). https://doi.org/10.1007/978-3-030-77961-0_11
  * [33] Purushotham, S., Meng, C., Che, Z., Liu, Y.: Benchmarking deep learning models on large healthcare datasets. Journal of Biomedical Informatics 83, 112–134 (2018). https://doi.org/10.1016/j.jbi.2018.04.007
  * [34] Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F.A., Bengio, Y., Courville, A.: On the Spectral Bias of Neural Networks (2019)
  * [35] Raissi, M., Perdikaris, P., Karniadakis, G.E.: Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations (2017)
  * [36] Ranjan, S., Senthamilarasu, S.: Applied Deep Learning and Computer Vision for Self-Driving Cars: Build autonomous vehicles using deep neural networks and behavior-cloning techniques. Packt Publishing (2020), https://books.google.es/books?id=nIX4DwAAQBAJ
  * [37] Ritz, W.: Über eine neue methode zur lösung gewisser variationsprobleme der mathematischen physik. Journal für die reine und angewandte Mathematik 135, 1–61 (1909), http://eudml.org/doc/149295
  * [38] Ruas, V.: An introduction to the mathematical foundations of the finite element method (2002), https://sites.icmc.usp.br/gustavo.buscaglia/cursos/mef/mef.pdf
  * [39] Ruthotto, L., Haber, E.: Deep Neural Networks Motivated by Partial Differential Equations. Journal of Mathematical Imaging and Vision 62 (2020). https://doi.org/10.1007/s10851-019-00903-1
  * [40] Samaniego, E., Anitescu, C., Goswami, S., Nguyen-Thanh, V., Guo, H., Hamdia, K., Zhuang, X., Rabczuk, T.: An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering 362, 112790 (2020). https://doi.org/10.1016/j.cma.2019.112790
  * [41] Scaman, K., Virmaux, A.: Lipschitz regularity of deep neural networks: analysis and efficient estimation. In: Proceedings of the 32nd International Conference on Neural Information Processing Systems. pp. 3839–3848 (2018)
  * [42] Shahriari, M., Pardo, D., Rivera, J.A., Torres-Verdín, C., Picon, A., Del Ser, J., Ossandón, S., Calo, V.: Error Control and Loss Functions for the Deep Learning Inversion of Borehole Resistivity Measurements. International Journal for Numerical Methods in Engineering (2020). https://doi.org/10.1002/nme.6593
  * [43] Sirignano, J., Spiliopoulos, K.: DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics 375, 1339 – 1364 (2018). https://doi.org/https://doi.org/10.1016/j.jcp.2018.08.029
  * [44] Wang, S., Yu, X., Perdikaris, P.: When and why PINNs fail to train: A neural tangent kernel perspective (2020)
  * [45] Weinan, E., Bing, Y.: The Deep Ritz Method: A Deep Learning-Based Numerical Algorithm for Solving Variational Problems. Communications in Mathematics and Statistics 6, 1–12 (2018). https://doi.org/10.1007/s40304-018-0127-z
  * [46] Weinzierl, S.: Introduction to Monte Carlo Methods. ArXiv High Energy Physics \- Phenomenology e-prints (07 2000). https://doi.org/10.1007/978-0-387-87837-9_1
  * [47] Wu, K., Xiu, D.: Data-driven deep learning of partial differential equations in modal space. Journal of Computational Physics 408, 109307 (2020). https://doi.org/10.1016/j.jcp.2020.109307
  * [48] Ying, X.: An Overview of Overfitting and its Solutions. Journal of Physics: Conference Series 1168, 022022 (feb 2019). https://doi.org/10.1088/1742-6596/1168/2/022022
  * [49] Zhang, C., Vinyals, O., Munos, R., Bengio, S.: A Study on Overfitting in Deep Reinforcement Learning abs/1804.06893 (2018)