# A deep first-order system least squares method for solving elliptic PDEs

Francisco M. Bersetche  Departamento de Matemática, Universidad de Buenos
Aires, Buenos Aires, Argentina<EMAIL_ADDRESS>and  Juan Pablo
Borthagaray  Centro de Matemática, Universidad de la República, Montevideo,
Uruguay<EMAIL_ADDRESS>
###### Abstract.

We propose a First-Order System Least Squares (FOSLS) method based on deep-
learning for numerically solving second-order elliptic PDEs. The method we
propose is capable of dealing with either variational and non-variational
problems, and because of its meshless nature, it can also deal with problems
posed in high-dimensional domains. We prove the $\Gamma$-convergence of the
neural network approximation towards the solution of the continuous problem,
and extend the convergence proof to some well-known related methods. Finally,
we present several numerical examples illustrating the performance of our
discretization.

FMB has been supported in part by a PEDECIBA postdoctoral fellowship and the
ANPCyT grant PICT 2018-3017.

##  1\. Introduction

Approximate solution of PDEs using machine learning techniques has been
considered in various forms in the past thirty years. For instance, [15, 16,
17, 21] propose to use neural networks to solve PDEs and ODEs. These articles
compute neural network solutions by using an a priori fixed mesh. In recent
years, there has been an incipient development of mesh-free numerical methods
to solve PDEs by using neural networks. Although the approaches have been
diverse, most of these algorithms aim to train a neural network to approximate
the unknown function, forcing the fulfillment of the PDE and its boundary
conditions through a suitable loss functional. In this regard, among other
works, let us mention [9, 28, 10, 29, 33, 22, 31, 19, 18].

Deep neural networks are not necessarily suitable for solving PDEs in low
dimensions, where they may be outperformed by classical methods specifically
tailored for the problems under consideration. However, neural network methods
have proven to be effective in some circumstances where the application of
classical methods becomes impractical. Such is the case of high-dimensional
PDEs. We refer to [30, 8] for discussion about the suitability of shallow
neural networks for solving high-dimensional PDEs.

The method we propose in this work aims to overcome some disadvantages of the
algorithms available in the literature. Using a first-order formulation we are
able to avoid the computation of second-order derivatives in cost functionals,
thereby saving a significant computational cost in high dimensions. Avoiding
second-order derivatives also allows us to use linear activation functions
and, a priori, gives us the possibility of approximating weak solutions. On
the other hand, counting on explicit representations of the gradients
simplifies the strong imposition of Neumann-type boundary conditions. Namely,
we can impose boundary conditions without adding penalty terms in the loss
function. This results in a reduction in training time. First order
formulations have also been used in [19, 18].

Let $\Omega\subset{\mathbb{R}^{d}}$ be an open domain. In this work, we shall
make use of the spaces

$\displaystyle H^{1}(\Omega)=\\{v\in L^{2}(\Omega)\colon\nabla v\in
L^{2}(\Omega)\\},$ $\displaystyle
H(\operatorname{div};\Omega)=\\{{{\bm{\psi}}}\in[L^{2}(\Omega)]^{d}\colon\operatorname{div}{{\bm{\psi}}}\in
L^{2}(\Omega)\\}.$

We assume there exists a disjoint partition
$\partial\Omega=\Gamma_{\mathcal{D}}\cup\Gamma_{\mathcal{N}}$, with
$|\Gamma_{\mathcal{D}}|>0$, and let ${{\bm{\nu}}}$ denote the outward normal
to $\Omega$. Given sufficiently regular functions
$f,g_{\mathcal{D}},g_{\mathcal{N}}$, we aim to solve the problem

(1.1) $\left\\{\begin{aligned} -\mbox{div}({{\bm{A}}}\nabla u)+Bu&=f&\mbox{ in
}\Omega,\\\ u&=g_{\mathcal{D}}&\mbox{ on }\Gamma_{\mathcal{D}},\\\
{{\bm{A}}}\nabla u\cdot{{\bm{\nu}}}&=g_{\mathcal{N}}&\mbox{ on
}\Gamma_{\mathcal{N}},\\\ \end{aligned}\right.$

where we assume ${{\bm{A}}}\in[L^{\infty}(\Omega)]^{d\times d}$ is a.e.
symmetric and uniformly positive definite: there exist constants
$\lambda,\Lambda$ such that

$0<\lambda\leq\lambda_{\min}({{\bm{A}}}(x))\leq\lambda_{\max}({{\bm{A}}}(x))\leq\Lambda,\quad\mbox{for
a.e. }x\in\Omega,$

where $\lambda_{\min}({{\bm{A}}}(\cdot))$ (resp.
$\lambda_{\max}({{\bm{A}}}(\cdot))$) denotes the minimum (resp. maximum)
eigenvalue of ${{\bm{A}}}(\cdot)$.

We assume the linear operator $B\colon H^{1}(\Omega)\to L^{2}(\Omega)$ in
(1.1) satisfies

$\|Bv\|_{L^{2}(\Omega)}\leq C\|\nabla v\|_{L^{2}(\Omega)}\quad\forall v\in
H^{1}(\Omega)\mbox{ such that }v=0\mbox{ on }\Gamma_{\mathcal{D}}.$

Examples satisfying this condition include
$Bv=\operatorname{div}({{\bm{\beta}}}v)$, with
${{\bm{\beta}}}\in[W^{1,\infty}(\Omega)]^{d}$, and
$Bv={{\bm{\beta}}}\cdot\nabla v+\gamma v$ for some
${{\bm{\beta}}}\in[L^{\infty}(\Omega)]^{d}$, $\gamma\in L^{\infty}(\Omega)$.
We thus remark that (1.1) can accommodate, for example, stationary convection-
reaction-diffusion problems. Following [6], we require problem (1.1) to be
invertible in $H^{1}(\Omega)$, namely, that for every $f\in H^{-1}(\Omega)$
there exists a weak solution $u\in H^{1}(\Omega)$ with $u=0$ on
$\Gamma_{\mathcal{D}}$.

We introduce the flux variable ${{\bm{\phi}}}={{\bm{A}}}\nabla u$ and rewrite
(1.1) as a first-order system:

(1.2) $\left\\{\begin{aligned} {{\bm{\phi}}}-{{\bm{A}}}\nabla u&=0&\mbox{ in
}\Omega,\\\ -\operatorname{div}{{\bm{\phi}}}+Bu-f&=0&\mbox{ in }\Omega,\\\
u&=g_{\mathcal{D}}&\mbox{ on }\Gamma_{\mathcal{D}},\\\
{{\bm{\phi}}}\cdot{{\bm{\nu}}}&=g_{\mathcal{N}}&\mbox{ on
}\Gamma_{\mathcal{N}}.\\\ \end{aligned}\right.$

Our approach is based on seeking minimizers of the loss function

(1.3) $\mathcal{L}(u,{{\bm{\phi}}}):=\|{{\bm{\phi}}}-{{\bm{A}}}\nabla
u\|_{L^{2}(\Omega)}^{2}+\|\operatorname{div}{{\bm{\phi}}}-Bu+f\|_{L^{2}(\Omega)}^{2}$

on a suitable set of admissible functions

$\mathcal{A}:=\\{{{\bm{q}}}=(u,{{\bm{\phi}}})\in H^{1}(\Omega)\times
H(\operatorname{div};\Omega)\colon u=g_{\mathcal{D}}\mbox{ on
}\Gamma_{\mathcal{D}},\ {{\bm{\phi}}}\cdot{{\bm{\nu}}}=g_{\mathcal{N}}\mbox{
on }\Gamma_{\mathcal{N}}\\}.$

Clearly, if (1.1) has a unique solution $u\in H^{1}(\Omega)$, then the unique
minimizer of $\mathcal{L}$ in $\mathcal{A}$ is
${{\bm{q}}}:=~{}(u,{{\bm{A}}}\nabla u)$. Our goal is to compute approximations
to such a minimizer within a suitable space
$\mathcal{A}_{m}\subset\mathcal{A}$. Particularly, in our method we consider a
space $\mathcal{A}_{m}$ composed of neural networks with a fixed architecture
and parameters ${{\bm{\Theta}}}\in\mathbb{R}^{m}$. Some efforts in this
direction include the deep FOSLS method from [5] and the deep mixed residual
method proposed in [20]. Reference [5] proposes the use of a partition of
$\Omega$ and a mid-point quadrature rule for the evaluation of the discrete
loss functional; instead, our algorithm is meshfree and uses random quadrature
points. In more recent work by three of the authors of that work [19], the use
of Monte Carlo integration is discussed albeit not pursued in detail. Such an
approach yields a significant advantage in high-dimensional problems. Our
method can be understood in the setting of the mixed residual methods in [20].
However, a significant difference between our work and [20] is that here we
propose a strong imposition of the boundary conditions instead of the
inclusion of penalization terms in the loss functional. We pre-train neural
networks to accommodate boundary data, which results in a reduction in the
number of iterations required in the solution of the PDE [3, 23].

The error in the approximation of continuous functionals with their discrete
counterparts is usually not taken into account in numerical methods based on
neural networks available in the literature; some recent efforts in this
direction include [24, 25, 34], where convergence rates are proved for a
certain class of elliptic functionals under strong regularity conditions on
the solution of the continuous problem. In other words, the focus is generally
on the convergence of the minimizers of functionals such as $\mathcal{L}$ in
(1.3) over certain neural network spaces towards the minimizer of the same
functional at the continuous level. However, in practice one does not compute
$\mathcal{L}$ exactly but rather approximates it by means of quadrature rules.
Let us call $\mathcal{L}_{N}$ such an approximation to the functional
$\mathcal{L}$, where $N$ is, for example, the number of quadrature points. The
computation of $\mathcal{L}_{N}$ instead of $\mathcal{L}$ can introduce
important changes in the nature of the minimization problem, such as the loss
of convexity of the associated functional [9]. A major contribution of this
work is to present a convergence analysis that considers the discretization of
the functional $\mathcal{L}$. Specifically, we prove the almost-sure
$\Gamma$-convergence of the discrete loss functions towards the continuous
one. As stated in Theorem 3.3, this implies the almost-sure convergence of the
solutions computed numerically to the solution of the continuous problem.

The techniques we develop for this purpose are not only valid for the method
we propose, and we generalize and apply them to the convergence analysis of a
broad class of methods, including the Deep Ritz [9] and the Deep Galerkin [28]
Methods (DRM and DGM, respectively; see Remarks 2 and 3).

### Organization of the paper.

The rest of the paper is organized as follows. Section 2 describes the method
we propose for dealing with (1.2), including the treatment of Dirichlet and
Neumann boundary conditions in strong form, and discusses some aspects
pertaining to its implementation. We perform a convergence analysis for our
method in Section 3. This analysis takes into account the approximation of the
loss functional by means of Monte Carlo integration, and establishes the
convergence of the discrete minimization problem towards the continuous one in
the sense of almost sure $\Gamma$-convergence. Section 4 generalizes the
analysis to include some other well-known methods, thereby establishing their
convergence as well. We illustrate the performance of our method through
computational examples in Section 5, and provide some concluding remarks in
Section 6.

##  2\. Description of the method

The goal of the method we propose is to approximate the unique minimizer
$(u,{{\bm{\phi}}})$ of the functional in (1.3). A natural first approach would
consist in seeking a set of parameters ${{\bm{\Theta}}}_{0}\in\mathbb{R}^{m}$
such that

$\mathcal{L}(u_{{{\bm{\Theta}}}_{0}},{{\bm{\phi}}}_{{{\bm{\Theta}}}_{0}})=\min_{{{\bm{\Theta}}}\in\mathbb{R}^{m}}\mathcal{L}(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}}),$

with the functions $(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})$
belonging to a suitable neural network space. The use of neural networks in
this setting has the advantage that one can easily implement meshfree methods
by randomly sampling collocation points (see [28, 11, 26, 27], for example),
and thereby be able to deal with high-dimensional problems, where most
classical numerical PDE methods become unfeasible.

The enforcement of boundary conditions is a non-trivial aspect to take into
account in this approach. A typical way to tackle this issue is to incorporate
boundary conditions by adding a penalization term [9, 28, 33]. However, in
practice it is observed that enforcing discrete functions to satisfy the
boundary conditions gives rise to a faster training process [3, 23]. We shall
first create suitable auxiliary functions with the purpose of imposing the
boundary conditions in a strong fashion. In this way, we ensure
($u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})\in\mathcal{A}_{m}\subset\mathcal{A}$,
for all ${{\bm{\Theta}}}\in\mathbb{R}^{m}$. Then, the optimization procedure
consists of sampling $N$ points $\\{x_{k}\\}_{k=1}^{N}\subset\Omega$
uniformly, and approximating
$\mathcal{L}(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})\approx\mathcal{L}_{N}(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}}),$
at every step of a gradient descent algorithm, with $\mathcal{L}_{N}$ defined
as

(2.1)
$\mathcal{L}_{N}(u,{{\bm{\phi}}}):=\frac{|\Omega|}{N}\sum_{k=1}^{N}\big{(}{{\bm{\phi}}}(x_{k})-{{\bm{A}}}\nabla
u(x_{k})\big{)}^{2}+\big{(}\operatorname{div}{{\bm{\phi}}}(x_{k})-Bu(x_{k})+f(x_{k})\big{)}^{2}.$

We expose the details below.

###  2.1. Strong imposition of boundary conditions

We follow the ideas from [3] about the imposition of Dirichlet boundary
conditions, and extend the approach to include Neumann boundary conditions.
Instead of trying to compute either $u$ or ${{\bm{\phi}}}$ directly and
incorporate the boundary conditions by a penalization term, we shall enforce
them in the construction of the neural network approximations. For that
purpose, we make use of the following notion.

######  Definition 2.1 (smooth distance function).

Let $\Gamma_{*}\subset\overline{\Omega}$ be a closed set. We say that a
Lipschitz continuous function $d_{*}\colon\Omega\to\mathbb{R}$ is a smooth
distance function if it satisfies $d_{*}\geq 0$ and $d_{*}(x)=0$ if and only
if $x\in\Gamma_{*}$.

We briefly comment on the use of smooth distance functions in the strong
imposition of Dirichlet and Neumann boundary conditions. In the computation of
$u$ in (1.3), we restrict the class of functions to be

(2.2) $u(x):=G_{\mathcal{D}}(x)+d_{\mathcal{D}}(x)\,v(x),$

where the unknown is the function $v\colon\Omega\to\mathbb{R}$,
$G_{\mathcal{D}}$ is a lifting of the Dirichlet datum, and $d_{\mathcal{D}}$
is a smooth distance function to $\Gamma_{\mathcal{D}}$.

In a similar fashion, we can incorporate normal boundary conditions on the
flux variable ${{\bm{\phi}}}$ in a strong way. We first construct a vector
field ${{\bm{n}}}\colon\Omega\to\mathbb{R}^{d}$ such that
${{\bm{n}}}|_{\Gamma_{\mathcal{N}}}={{\bm{\nu}}}$ and $|{{\bm{n}}}(x)|=1$ for
a.e. $x\in\Omega$, and consider

(2.3)
${{\bm{\phi}}}(x):={{\bm{\psi}}}(x)+\left(G_{\mathcal{N}}(x)-\frac{{{\bm{\psi}}}(x)\cdot{{\bm{n}}}(x)}{1+d_{\mathcal{N}}(x)}\right){{\bm{n}}}(x).$

Above, $G_{\mathcal{N}}$ is a lifting of the Neumann boundary condition,
$d_{\mathcal{N}}$ is a smooth distance function to $\Gamma_{\mathcal{N}}$, and
the unknown is the function ${{\bm{\psi}}}\colon\Omega\to\mathbb{R}^{d}$. By
its definition, the function ${{\bm{\phi}}}$ satisfies the boundary condition
${{\bm{\phi}}}\cdot{{\bm{\nu}}}=g_{\mathcal{N}}$ at $\Gamma_{\mathcal{N}}$. We
remark that we do not require any smoothness on ${{\bm{n}}}$: in particular
this field may be discontinuous at some points in the domain.

Therefore, in the construction of approximate solutions we shall first compute
the vector field ${{\bm{n}}}$ and the scalar functions $d_{\mathcal{D}}$,
$d_{\mathcal{N}}$, $G_{\mathcal{D}}$, $G_{\mathcal{N}}$. Then, we seek
${{\bm{y}}}=(v,{{\bm{\psi}}})$ such that the corresponding pair
$(u,{{\bm{\phi}}})$, given by (2.2) and (2.3), minimizes the loss function
$\mathcal{L}$. The computation of the auxiliary functions ${{\bm{n}}}$,
$d_{\mathcal{D}}$, $d_{\mathcal{N}}$, $G_{\mathcal{D}}$, $G_{\mathcal{N}}$
typically requires fewer degrees of freedom and iterations than the
computation of $(v,{{\bm{\psi}}})$, depending on the complexity of the domain
or the boundary data. Consequently, we shall frequently use a simpler
architecture to represent them. Below, we give details on the computation of
the auxiliary functions.

####  2.1.1. Computation of smooth distance functions

Loosely, for $*\in\\{\mathcal{D},\mathcal{N}\\}$, a smooth distance function
to $\Gamma_{*}$ is a function $d_{*}\colon\Omega\to[0,\infty)$ that
approximates the distance to $\Gamma_{*}$, cf. Definition 2.1. To construct
such functions, we first randomly choose $N_{d}$ points
$\\{x_{i}\\}_{i=1}^{N_{d}}\subset\Omega$ (the same set of points can be used
for either $*=\mathcal{D}$ and $*=\mathcal{N}$) and compute

$d^{(*)}(x_{i})\approx\mbox{dist}(x_{i},\Gamma_{*}).$

This can be done by choosing points on $\Gamma_{*}$ and using efficient
nearest-neighbor search strategies. Once we have computed the quantities
$\\{d^{(*)}(x_{i})\\}_{i=1}^{N_{d}}$, we train a neural network for $d_{*}$ by
using the cost function

$\mathcal{L}_{*}(d)=\frac{1}{N_{d}}\sum_{i=1}^{N}|d(x_{i})-d^{(*)}(x_{i})|^{2}+\frac{1}{N_{d,*}}\sum_{i=1}^{N_{d,*}}|d(x_{*,i})|^{2},$

where $\\{x_{*,i}\\}_{i=1}^{N_{d,*}}$ is a random batch of points on
$\Gamma_{*}$.

In the setting of $d_{\mathcal{D}}$ and $d_{\mathcal{N}}$, we use neural
networks with a single hidden layer and significantly less parameters than the
networks employed in the PDE resolution.

####  2.1.2. Boundary data liftings and normal field

We approximate liftings of the boundary data to $\Omega$ by smooth liftings
[3]: in either (2.2) and (2.3), we require $G_{\mathcal{D}}$ and
$G_{\mathcal{N}}$ to coincide with $g_{\mathcal{D}}$ on $\Gamma_{\mathcal{D}}$
and with $g_{\mathcal{N}}$ on $\Gamma_{\mathcal{N}}$, respectively, and to be
smooth enough so that we can apply the differential operator to them
pointwise. A natural way to enforce the former is to set the $L^{2}$-norms of
the discrepancies on the corresponding boundary subsets as loss functions,
namely

$\mathcal{L}_{\mathcal{D}}(G)=\|G-g_{\mathcal{D}}\|_{L^{2}(\Gamma_{\mathcal{D}})}^{2},\quad\mathcal{L}_{\mathcal{N}}(G)=\|G-g_{\mathcal{N}}\|_{L^{2}(\Gamma_{\mathcal{N}})}^{2}.$

In practice, we consider sets of boundary nodes
$\\{z^{\mathcal{D}}_{i}\\}_{i=1}^{M_{\mathcal{D}}}\subset\Gamma_{\mathcal{D}},$
$\\{z^{\mathcal{N}}_{i}\\}_{i=1}^{M_{\mathcal{N}}}\subset\Gamma_{\mathcal{N}}$
and define the quadratic cost functionals

$\mathcal{L}_{\mathcal{D}}(G)=\frac{1}{M_{\mathcal{D}}}\sum_{i=1}^{M_{\mathcal{D}}}|G(z^{\mathcal{D}}_{i})-g_{\mathcal{D}}(z^{\mathcal{D}}_{i})|^{2},\qquad\mathcal{L}_{\mathcal{N}}(G)=\frac{1}{M_{\mathcal{N}}}\sum_{i=1}^{M_{\mathcal{N}}}|G(z^{\mathcal{N}}_{i})-g_{\mathcal{N}}(z^{\mathcal{N}}_{i})|^{2}.$

In the same fashion as for the smooth distance functions, we consider neural
networks with a single hidden layer to compute the functions $d_{\mathcal{D}}$
and $d_{\mathcal{N}}$.

Analogously, for the computation of the vector field ${{\bm{n}}}$ we start
from the loss function

$\mathcal{L}_{{{\bm{n}}}}({{\bm{m}}})=\|{{\bm{m}}}-{{\bm{\nu}}}\|_{L^{2}(\Gamma_{\mathcal{N}})}^{2}+\||{{\bm{m}}}|^{2}-1\|_{L^{2}(\Omega)}^{2},$

consider a set of randomly selected points
$\\{z^{\mathcal{N},{{\bm{n}}}}_{i}\\}_{i=1}^{M_{\mathcal{N},{{\bm{n}}}}}\subset\Gamma_{\mathcal{N}}$
and $\\{z^{{{\bm{n}}}}_{i}\\}_{i=1}^{M_{{{\bm{n}}}}}\subset\Omega$, and
minimize the cost functional

$\mathcal{L}_{{{\bm{n}}}}({{\bm{m}}})=\frac{1}{M_{\mathcal{N},{{\bm{n}}}}}\sum_{i=1}^{M_{\mathcal{N},{{\bm{n}}}}}|{{\bm{m}}}(z^{\mathcal{N},{{\bm{n}}}}_{i})-{{\bm{\nu}}}|^{2}+\frac{1}{M_{{{\bm{n}}}}}\sum_{i=1}^{M_{{{\bm{n}}}}}||{{\bm{m}}}(z^{{{\bm{n}}}}_{i})|^{2}-1|^{2}.$

We point out that, in practice, the set of auxiliary points
$\\{z^{\mathcal{N},{{\bm{n}}}}_{i}\\}_{i=1}^{M_{\mathcal{N},{{\bm{n}}}}}$ can
be the same as the set $\\{z^{\mathcal{N}}_{i}\\}_{i=1}^{M_{\mathcal{N}}}$
used in the approximation of $\mathcal{L}_{\mathcal{N}}$.

###  2.2. Computational aspects

Once we have built the auxiliary functions, we proceed to compute $u$ and
${{\bm{\phi}}}$. For this purpose, we consider a set of random points
$\\{x_{k}\\}_{k=1}^{N}\subset\Omega$, and seek to minimize the cost functional

(2.4)
$\mathcal{L}_{N}(u,{{\bm{\phi}}}):=\frac{|\Omega|}{N}\sum_{k=1}^{N}\big{(}{{\bm{\phi}}}(x_{k})-{{\bm{A}}}\nabla
u(x_{k})\big{)}^{2}+\big{(}\operatorname{div}{{\bm{\phi}}}(x_{k})-Bu(x_{k})+f(x_{k})\big{)}^{2}.$

From the construction of $u$ and ${{\bm{\phi}}}$ (see (2.2) and (2.3)), the
trainable parameters ${{\bm{\Theta}}}$ arise in the computation of the
auxiliary functions $v$ and ${{\bm{\psi}}}$.

In broad terms, the method we propose can be summarized as follows:

  * •

Stage 1: Train auxiliary functions $d_{\mathcal{D}}$, $d_{\mathcal{N}}$,
$G_{\mathcal{D}}$, $G_{\mathcal{N}}$, and ${{\bm{n}}}$.

  * •

Stage 2: Until some stop criterion is reached, do:

    * –

Select random points $\\{x_{k}\\}_{k=1}^{N}\subset\Omega$.

    * –

For some learning rate $\ell$, do:

${{\bm{\Theta}}}={{\bm{\Theta}}}-\ell\nabla_{{{\bm{\Theta}}}}\mathcal{L}_{N}(u_{{\bm{\Theta}}},{{\bm{\phi}}}_{{\bm{\Theta}}}).$

    * –

Update learning rate.

The computation of $\mathcal{L}_{N}(u,{{\bm{\phi}}})$ requires computing the
derivatives of $u$ and ${{\bm{\phi}}}$ with respect to the input variables,
evaluated at $\\{x_{k}\\}_{k=1}^{N}$. Since we constructed our auxiliary
functions as neural networks, it is possible to compute efficiently these
derivatives by means of the Back-Propagation algorithm. Packages like
TensorFlow allow this kind of computation.

Additionally, our least-squares loss function (2.4) only involves first-order
derivatives in space. We discretize such derivatives by using finite-
difference quotients. Namely, for any function
${{\bm{\varphi}}}\colon{\mathbb{R}^{d}}\to\mathbb{R}^{n}$ we let $h>0$ be a
fixed constant and consider the second-order (with respect to $h$) formula

$\partial_{i}{{\bm{\varphi}}}(x_{k})\simeq\frac{{{\bm{\varphi}}}(x_{k}+h{{\bm{e}}}_{i})-{{\bm{\varphi}}}(x_{k}-h{{\bm{e}}}_{i})}{2h},$

where ${{\bm{e}}}_{i}\in\mathbb{R}^{d}$ is the $i$-th canonical basis vector
in ${\mathbb{R}^{d}}$. We employ this formula for the approximation of $\nabla
u$, $\operatorname{div}{{\bm{\phi}}}$ and the first-order derivatives involved
in $B$.

For the numerical examples we implemented our algorithm by using PyTorch and
discretizing the differential operators by means of finite differences. We
typically use about 10,000 steps of gradient descent, sampling between 1,000
and 5,000 random points in $\Omega$ at each step. A step-type decrease in the
learning rate showed good results in practice. In particular, we start from a
learning rate $\ell=10^{-2}$, which we halve every 1,000-2,500 gradient
descent steps. No particular type of architecture was chosen for the functions
involved. We use three-layer neural networks with linear activation function
(ReLU) for the auxiliary functions, and five-layer networks for the main
variables $v$ and ${{\bm{\psi}}}$. The ADAM [14] optimization algorithm showed
good results in numerical experiments. Further details about the
implementation of the method can be found in Section 5.

Regarding the training of auxiliary functions $d_{\mathcal{D}}$ and
$d_{\mathcal{N}}$, the following procedure showed good results in practice:

  * •

Select $N_{d}$ random points $\\{x_{k}\\}_{k=1}^{N_{d}}\subset\Omega$.

  * •

Initialize a vector ${{\bm{D}}}$ as ${{\bm{D}}}_{i}=\infty$ for
$i=1,...,N_{d}$.

  * •

Until some stop criterion is reached, do:

    * –

Select $M_{*}$ random points $\\{z^{*}_{i}\\}_{i=1}^{M_{*}}\subset\Gamma_{*}$.

    * –

Update ${{\bm{D}}}$ as: ${{\bm{D}}}_{k}=\min\\{\min_{i=1,\ldots
M_{*}}{|x_{k}-z^{*}_{i}|},{{\bm{D}}}_{k}\\}$

    * –

Define the loss function:

$\mathcal{L}_{*}(d_{*})=\frac{1}{N_{d}}\sum_{k=1}^{N_{d}}|d_{*}(x_{k})-{{\bm{D}}}_{k}|^{2}+\frac{1}{M_{*}}\sum_{i=1}^{M_{*}}|d(z_{i}^{*})|^{2}.$

    * –

For some learning rate $\ell$, do:

${{\bm{\Theta}}}_{d_{*}}={{\bm{\Theta}}}_{d_{*}}-\ell\nabla_{{{\bm{\Theta}}}_{d_{*}}}\mathcal{L}_{*}(d_{*}).$

    * –

Update learning rate.

Here $*\in\\{D,N\\}$, and ${{\bm{\Theta}}}_{d_{*}}$ denotes the trainable
parameters of $d_{*}$.

##  3\. Analysis of the method

In this section, we prove the convergence of our method by using two main
ingredients. First, we put the discretization in a $\Gamma$-convergence
framework. More precisely, the sequence of functionals we consider is related
to the use of meshfree methods in the computation of a regularized version of
the discrete loss functional
$\mathbb{R}^{m}\ni{{\bm{\Theta}}}\mapsto\mathcal{L}(u_{{\bm{\Theta}}},{{\bm{\phi}}}_{{\bm{\Theta}}})$;
see Theorem 3.2 below. Second, we exploit the coercivity of the least-squares
functional and approximation properties of neural networks to conclude that
the sequence of minimizers of the regularized discrete loss functionals
converges to the solution of (1.1) as the number of neural network parameters
$m\to\infty$.

For the sake of simplicity, we consider problem (1.2) with
$g_{\mathcal{D}}=g_{\mathcal{N}}=0$. Otherwise, one could consider
$G_{\mathcal{D}}$ and $G_{\mathcal{N}}$ such that
$G_{\mathcal{D}}=g_{\mathcal{D}}$ on $\Gamma_{\mathcal{D}}$ and
$G_{\mathcal{N}}=g_{\mathcal{N}}$ on $\Gamma_{\mathcal{N}}$, a smooth normal
field ${{\bm{n}}}$ such that ${{\bm{n}}}={{\bm{\nu}}}$ on
$\Gamma_{\mathcal{N}}$, and then the auxiliary functions
$u_{0}=u-G_{\mathcal{D}}$ and
${{\bm{\phi}}}_{0}={{\bm{\phi}}}-G_{\mathcal{N}}{{\bm{n}}}$ would solve the
first-order system

$\left\\{\begin{array}[]{rll}{{\bm{\phi}}}_{0}-{{\bm{A}}}\nabla
u_{0}=&{{\bm{A}}}\nabla G_{\mathcal{D}}-G_{\mathcal{N}}{{\bm{n}}}&\mbox{in
}\Omega,\\\
-\operatorname{div}({{\bm{\phi}}}_{0})+Bu_{0}=&f+\operatorname{div}(G_{\mathcal{N}}{{\bm{n}}})-BG_{\mathcal{D}}&\mbox{in
}\Omega,\\\ u_{0}=&0&\mbox{on }\Gamma_{\mathcal{D}},\\\
{{\bm{\phi}}}_{0}\cdot{{\bm{\nu}}}=&0&\mbox{on }\Gamma_{\mathcal{N}}.\\\
\end{array}\right.$

Naturally, the solution to this system corresponds to the minimum of the
least-squares functional

$(u,{{\bm{\phi}}})\mapsto\|{{\bm{\phi}}}-{{\bm{A}}}\nabla
u+\widetilde{g}\|_{L^{2}(\Omega)}^{2}+\|\operatorname{div}({{\bm{\phi}}})-Bu+\widetilde{f}\|_{L^{2}(\Omega)}^{2},$

with $\widetilde{g}=-{{\bm{A}}}\nabla
G_{\mathcal{D}}+G_{\mathcal{N}}{{\bm{n}}}$ and
$\widetilde{f}=f+\operatorname{div}(G_{\mathcal{N}}{{\bm{n}}})-BG_{\mathcal{D}}$.
This functional can be dealt with by using the same tools as for (1.3), the
only difference being the presence of the zero-order correction term
$\widetilde{g}$ in the first $L^{2}$-norm.

In the following proof of convergence, we restrict ourselves to one hidden
layer neural networks with $n$ neurons. We define the set of discrete
functions

$\mathcal{C}_{m}:=\Big{\\{}(v_{{\bm{\Theta}}},{{\bm{\psi}}}_{{\bm{\Theta}}}):v_{{\bm{\Theta}}}=B_{v}\sigma(A_{v}x+c_{v}),{{\bm{\psi}}}_{{\bm{\Theta}}}=B_{{{\bm{\psi}}}}\sigma(A_{{{\bm{\psi}}}}x+c_{{{\bm{\psi}}}})\Big{\\}},$

with $A_{v},A_{{{\bm{\psi}}}}\in\mathbb{R}^{n\times d}$,
$c_{v},c_{{{\bm{\psi}}}}\in\mathbb{R}^{n\times 1}$,
$B_{v}\in\mathbb{R}^{1\times n}$, $B_{{{\bm{\psi}}}}\in\mathbb{R}^{d\times
n}$, and $\sigma\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$, and $\sigma$ a smooth
and bounded non-constant activation function, applied elementwise. We collect
all the parameters in ${{\bm{\Theta}}}\in\mathbb{R}^{m}$ with $m=3n(d+1)$. We
remark that, whenever we state that $m\to\infty$, we mean that the number of
neurons $n$ is growing to infinity.

Assuming that we are able to construct smooth auxiliary functions
$d_{\mathcal{D}}$, $d_{\mathcal{N}}$ and ${{\bm{n}}}$ as in Section 2.1, we
define the set of discrete admissible functions

(3.1)
$\mathcal{A}_{m}:=\Big{\\{}{{\bm{q}}}_{{{\bm{\Theta}}}}=(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}}):u_{{{\bm{\Theta}}}}=d_{\mathcal{D}}v_{{{\bm{\Theta}}}}\mbox{
and
}{{\bm{\phi}}}_{{{\bm{\Theta}}}}={{\bm{\psi}}}_{{{\bm{\Theta}}}}-\Big{(}\frac{{{\bm{\psi}}}_{{{\bm{\Theta}}}}\cdot{{\bm{n}}}}{1+d_{\mathcal{N}}}\Big{)}{{\bm{n}}},\mbox{
}(v_{{{\bm{\Theta}}}},{{\bm{\psi}}}_{{{\bm{\Theta}}}})\in\mathcal{C}_{m}\Big{\\}}.$

We remark that the fulfillment of the boundary conditions is guaranteed within
the set $\mathcal{A}_{m}$, in the sense that $u_{{\bm{\Theta}}}=0$ if
$d_{\mathcal{D}}=0$ and ${{\bm{\phi}}}_{{\bm{\Theta}}}\cdot{{\bm{n}}}=0$ if
$d_{\mathcal{N}}=0$.

######  Remark 1.

Naturally, when using Montecarlo integration, one is not allowed to take
pointwise evaluations of an arbitrary function $f\in L^{2}(\omega)$. By
density, for every $\epsilon>0$ we can find a continuous function
$f_{\epsilon}$ with $\|f-f_{\epsilon}\|_{L^{2}(\Omega)}<\epsilon$. By the
ellipticity of the functional $\mathcal{L}$ in the $H^{1}(\Omega)\times
H(\operatorname{div};\Omega)$ norm (cf. (3.5) below), if we let
$\mathcal{L}_{\epsilon}$ be the functional (1.3) using $f_{\epsilon}$ instead
of $f$ and ${{\bm{q}}}_{0,\epsilon}$ its minimizer, we then have
$\|{{\bm{q}}}_{0}-{{\bm{q}}}_{0,\epsilon}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}<\epsilon$. We can therefore implement the method
by using $f_{\epsilon}$ instead of $f$ and letting $\epsilon\to 0$ as
$m\to\infty$.

Nevertheless, we emphasize that, for the sake of the theoretical results in
this paper, for any $f\in L^{2}(\Omega)$ we can take any representative of the
equivalence class of $f$ in the definition of the functionals (e.g. in (4.2)).
Our convergence results are not affected because they are stated in an “almost
sure” sense.

###  3.1. Approximation properties of neural networks

Let ${{\bm{q}}}_{0}=(u_{0},{{\bm{\phi}}}_{0})\in\mathcal{A}$ be the unique
minimizer of (1.3). We shall make the assumption that ${{\bm{q}}}_{0}$ can be
approximated by the neural network spaces. Namely, let us assume that

(3.2)
$d({{\bm{q}}}_{0},\mathcal{A}_{m}):=\inf_{{{\bm{q}}}_{{\bm{\Theta}}}\in\mathcal{A}_{m}}\|{{\bm{q}}}_{0}-{{\bm{q}}}_{{\bm{\Theta}}}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}\to 0\quad\mbox{as }m\to\infty.$

We briefly comment on this hypothesis. In first place, there are several by
now classical results [7, 13, 2] regarding the approximation properties of
neural networks, although without the incorporation of boundary conditions. We
additionally point out to [32, 12] for recent results regarding approximation
capabilities of ReLU neural networks, including approximation rates. For deep
ReLU neural networks (with at most $\lceil\log_{2}(d+1)\rceil$ hidden layers),
references [12, 1] establish the capability of networks to represent
simplicial linear finite element functions, which possess good approximation
properties in the $H^{1}$-norm. Therefore, if we use a nonconstant activation
function $\sigma$, then we expect $d({{\bm{q}}},\mathcal{C}_{m})\to 0$ when
$m\to\infty$ for any ${{\bm{q}}}\in H^{1}(\Omega)\times
H(\operatorname{div};\Omega)$.

Condition (3.2) further assumes that the solution $q_{0}$ can be approximated
through the admissible classes $\mathcal{A}_{m}$ that incorporate boundary
conditions. This hypothesis holds, for example, if one assumes certain
regularity of solutions to (1.2). For instance, if $u_{0}\in
C^{1}(\overline{\Omega})$, then it satisfies (recall $g_{\mathcal{D}}=0$)

$\left|\lim_{t\to
0^{+}}\frac{u_{0}(z-t{{\bm{\nu}}})}{t}\right|=\left|\frac{\partial
u_{0}}{\partial{{\bm{\nu}}}}(z)\right|<\infty,\quad z\in\Gamma_{\mathcal{D}}.$

If we write $x=z-t{{\bm{\nu}}}$, then
$t\approx\mbox{dist}(x,\Gamma_{\mathcal{D}})\approx d_{\mathcal{D}}(x)$ and
the finiteness of the limit above essentially means that
$u_{0}/d_{\mathcal{D}}$ is a bounded function. Additionally, if we can
construct auxiliary functions $d_{\mathcal{D}}$, $d_{\mathcal{N}}$, and
${{\bm{n}}}$ in such a way that

(3.3) $\frac{u_{0}}{d_{\mathcal{D}}}\in H^{1}(\Omega),\mbox{ and
}\frac{({{\bm{\phi}}}_{0}\cdot{{\bm{n}}}){{\bm{n}}}}{d_{\mathcal{N}}}\in
H(\operatorname{div};\Omega).$

then there exists a sequence
$\\{(v_{m},{{\bm{\psi}}}_{m})\\}_{m\in\mathbb{N}}$ with
$(v_{m},{{\bm{\psi}}}_{m})\in\mathcal{C}_{m}$ for all $m$, such that

$\left\|v_{m}-\frac{u_{0}}{d_{\mathcal{D}}}\right\|_{H^{1}(\Omega)}\to
0\quad\mbox{and}\quad\left\|{{\bm{\psi}}}_{m}-\sum^{d-1}_{i=1}({{\bm{\phi}}}_{0}\cdot{{\bm{t}}}){{\bm{t}}}-\frac{1+d_{\mathcal{N}}}{d_{\mathcal{N}}}({{\bm{\phi}}}_{0}\cdot{{\bm{n}}}){{\bm{n}}}\right\|_{H(\operatorname{div};\Omega)}\to
0$

as $m\to\infty$. Defining the sequence
$\\{(u_{m},{{\bm{\phi}}}_{m})\\}_{m\in\mathbb{N}}$ as
$u_{m}=d_{\mathcal{D}}v_{m}$ and
${{\bm{\phi}}}_{m}={{\bm{\psi}}}_{m}-\Big{(}\frac{{{\bm{\psi}}}_{m}\cdot{{\bm{n}}}}{1+d_{\mathcal{N}}}\Big{)}{{\bm{n}}}$,
we would have $(u_{m},{{\bm{\phi}}}_{m})\in\mathcal{A}_{m}$ for all $m$, and
$(u_{m},{{\bm{\phi}}}_{m})\to(u_{0},{{\bm{\phi}}}_{0})$ in
$\|\cdot\|_{H^{1}(\Omega)\times H(\operatorname{div};\Omega)}$ and therefore
(3.2) would hold. Clearly, (3.3) is a regularity assumption on the solution of
(1.2), and in turn it translates into its approximability by neural networks.

###  3.2. $\Gamma$-convergence

We aim to prove the convergence of the neural network approximations computed
by our method towards minimizers of the least-squares functional $\mathcal{L}$
in (1.3). For this purpose, we shall make use of $\Gamma$-convergence theory,
that provides a framework for the convergence of functionals. In particular,
if one has proven the $\Gamma$-convergence of a sequence of functionals and
has a converging sequence of minimizers, then one can guarantee the existence
of solutions to the limit problem, as well as the convergence of either
minimum values and minimizers. We next briefly review the definition and some
basic results pertaining to $\Gamma$-convergence and refer to [4] for further
details.

######  Definition 3.1 (sequential $\Gamma$-convergence).

Let $X$ be a metric space and let $F_{n}$, $F:X\to\overline{\mathbb{R}}$,
where $\overline{\mathbb{R}}:=[-\infty,+\infty]$. We say that $F_{n}$
$\Gamma$-converges to $F$ (and write $F_{n}\xrightarrow[]{\Gamma}F$) if, for
every $x\in X$ we have

  * •

_(lim-inf inequality)_ for every sequence $\\{x_{n}\\}_{n\in\mathbb{N}}\subset
X$ converging to $x$,

$F(x)\leq\liminf_{n\to\infty}F_{n}(x_{n});$

  * •

_(lim-sup inequality)_ there exists a sequence $\\{x_{n}\\}_{n\in\mathbb{N}}$
converging to $x$ such that

$F(x)\geq\limsup_{n\to\infty}F_{n}(x_{n}).$

######  Definition 3.2 (equi-coercivity).

Let $\\{F_{n}\\}_{n\in\mathbb{N}}$ be a sequence of functions
$F_{n}:X\to\overline{\mathbb{R}}$. We say that $\\{F_{n}\\}$ is equi-coercive
if for all $t\in\mathbb{R}$ there exists a compact set $K_{t}\subset X$ such
that $\\{F_{n}\leq t\\}\subset K_{t}$.

######  Theorem 3.1 (fundamental theorem of $\Gamma$-convergence).

Let $(X,d)$ be a metric space, $\\{F_{n}\\}_{n\in\mathbb{N}}$ be an equi-
coercive sequence of functions on $X$, and $F$ be such that
$F_{n}\xrightarrow[]{\Gamma}F$. Then,

$\exists\min_{X}F=\lim_{n\to\infty}\inf_{X}F_{n}.$

Moreover, if $\\{x_{n}\\}_{n\in\mathbb{N}}$ is a precompact sequence in $X$
such that $\lim_{n\to\infty}F_{n}(x_{n})=\lim_{n\to\infty}\inf_{X}F_{n}$, then
every limit of a subsequence of $\\{x_{n}\\}$ is a minimum point for $F$.

We emphasize that the result above guarantees that the equi-coercivity of a
family of functionals combined with their $\Gamma$-convergence yields the
convergence of the minimizers towards the minimizers of the $\Gamma$-limit.

###  3.3. Convergence of the method

We split the proof of convergence of our method into several steps. We start
by proving the following auxiliary lemma, that shows the continuity of the
neural network functions with respect to the parameters.

######  Lemma 3.1 (continuity with respect to neural network parameters).

The map

${{\bm{\Theta}}}\mapsto{{\bm{q}}}_{{\bm{\Theta}}}=(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})\in(\mathcal{A}_{m},\|\cdot\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)})$

is continuous. Moreover, defining the functions
$G_{1},G_{2}:\mathbb{R}^{m}\times\Omega\to\mathbb{R}$,

(3.4)
$G_{1}({{\bm{\Theta}}},x):=|{{\bm{\phi}}}_{{{\bm{\Theta}}}}(x)-{{\bm{A}}}\nabla
u_{{{\bm{\Theta}}}}(x)|^{2},\quad
G_{2}({{\bm{\Theta}}},x):=|\operatorname{div}{{\bm{\phi}}}_{{{\bm{\Theta}}}}(x)-Bu_{{{\bm{\Theta}}}}(x)+f(x)|^{2},$

for any $R>0$ we have $G_{1}\in L^{\infty}(B(0,R)\times\Omega)$ and, assuming
$f\in L^{2}(\Omega)$, there exists a function $s\in
L^{1}(B(0,R)\times\Omega))$, depending on $R$, such that
$|G_{2}({{\bm{\Theta}}},x)|\leq s({{\bm{\Theta}}},x)$ for all
$({{\bm{\Theta}}},x)\in B(0,R)\times\Omega$.

###### Proof.

Let us first focus on a generic neural network
$v_{{{\bm{\Theta}}}}:\mathbb{R}^{d}\to\mathbb{R}$ with one hidden layer,

$v_{{{\bm{\Theta}}}}(x)=B\sigma(Ax+c).$

Above, we assume $\sigma$ is a Lipschitz continuous activation function, and
the parameters $B\in\mathbb{R}^{1\times n}$, $A\in\mathbb{R}^{n\times d}$ and
$c\in\mathbb{R}^{n\times 1}$ are collected in
${{\bm{\Theta}}}\in\mathbb{R}^{m}$, $m=n(d+2)$. Using the fact that
$v_{{{\bm{\Theta}}}}$ and its derivatives depend continuously on the
parameters, one can verify easily that the map $\mathbb{R}^{m}\mapsto
W^{1,\infty}(\Omega)$ such that ${{\bm{\Theta}}}\mapsto v_{{{\bm{\Theta}}}}$
is continuous. Moreover, the function
$G:\mathbb{R}^{m}\times\Omega\to\mathbb{R}$, defined as
$G({{\bm{\Theta}}},x):=v_{{{\bm{\Theta}}}}(x)$ is Lipschitz continuous, and
therefore it is bounded on $B(0,R)\times\Omega$ and its (weak) derivatives are
essentially bounded on the same set as well. Furthermore, if $f\in
L^{2}(\Omega)$ then $|G({{\bm{\Theta}}},x)+f(x)|^{2}\leq
2|G({{\bm{\Theta}}},x)|^{2}+2|f(x)|^{2}\leq
2M+2|f(x)|^{2}=:s({{\bm{\Theta}}},x)$, with $s\in L^{1}(B(0,R)\times\Omega)$.

For arbitrary neural network functions
$(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})$ in the space
$\mathcal{A}_{m}$, defined by (3.1), we exploit the idea above together with
the fact that the auxiliary functions $d_{\mathcal{D}},d_{\mathcal{N}}$ and
${{\bm{n}}}$ are smooth to conclude the desired result. ∎

The following lemma guarantees that, for the loss function $\mathcal{L}$
defined in (1.3), quasi-minimizers over $\mathcal{A}_{m}$ converge towards the
minimizer ${{\bm{q}}}_{0}\in\mathcal{A}$ as $m\to\infty$.

######  Lemma 3.2 (approximation properties of $\mathcal{A}_{m}$).

For every $m\in\mathbb{N}$, let us define the set of neural network quasi-
minimizers

$\mathcal{I}_{m}:=\\{{{\bm{q}}}\in\mathcal{A}_{m}:\mathcal{L}({{\bm{q}}})\leq\mathcal{L}({{\bm{q}}}^{*})+1/m\
\forall{{\bm{q}}}^{*}\in\mathcal{A}_{m}\\}.$

Then, if ${{\bm{q}}}_{0}$ is the unique minimizer of $\mathcal{L}$ in
$\mathcal{A}$, we have

$\sup_{{{\bm{q}}}_{m}\in\mathcal{I}_{m}}\|{{\bm{q}}}_{m}-{{\bm{q}}}_{0}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}\to 0\quad\mbox{as }m\to\infty.$

###### Proof.

From [6], we know that $\mathcal{L}$ is elliptic with respect to the
$H^{1}(\Omega)\times H(\operatorname{div};\Omega)$ norm. Namely, there exist
positive constants $\alpha$ and $\beta$ such that

(3.5) $\alpha\|(u,{{\bm{\phi}}})\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}\leq\|{{\bm{\phi}}}-{{\bm{A}}}\nabla
u\|_{0}^{2}+\|\operatorname{div}({{\bm{\phi}}})-Bu\|_{0}^{2}\leq\beta\|(u,{{\bm{\phi}}})\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)},$

for all $(u,{{\bm{\phi}}})\in H^{1}(\Omega)\times
H(\operatorname{div};\Omega)$.

Let $\varepsilon>0$. By (3.2), we consider $m_{0}>0$ such that
$d({{\bm{q}}}_{0},\mathcal{A}_{m})<\varepsilon$ and $1/m<\varepsilon$ for all
$m>m_{0}$. For every $m>0$, there exists
${{\bm{q}}}^{*}_{m}=(u^{*}_{m},{{\bm{\phi}}}^{*}_{m})\in\mathcal{A}_{m}$ with
$d({{\bm{q}}}_{0},\mathcal{A}_{m})\geq\|{{\bm{q}}}^{*}_{m}-{{\bm{q}}}_{0}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}-\varepsilon$. Then, for all $m>m_{0}$ and every
neural network quasi-minimizer
${{\bm{q}}}_{m}=(u_{m},{{\bm{\phi}}}_{m})\in\mathcal{I}_{m}$, using that the
solution ${{\bm{q}}}_{0}=(u_{0},{{\bm{\phi}}}_{0})$ of (1.2) satisfies the
conditions ${{\bm{\phi}}}_{0}={{\bm{A}}}\nabla u_{0}$ and
$-\operatorname{div}({{\bm{\phi}}}_{0})+Bu_{0}=f$ a.e. in $\Omega$ and
exploiting the upper bound in (3.5), we have

$\displaystyle 0$
$\displaystyle\leq\mathcal{L}({{\bm{q}}}_{m})\leq\mathcal{L}({{\bm{q}}}^{*}_{m})+\varepsilon=\|{{\bm{\phi}}}^{*}_{m}-{{\bm{A}}}\nabla
u^{*}_{m}\|_{0}^{2}+\|\operatorname{div}({{\bm{\phi}}}^{*}_{m})-Bu^{*}_{m}+f\|_{0}^{2}+\varepsilon$
$\displaystyle=\|{{\bm{\phi}}}^{*}_{m}-{{\bm{\phi}}}_{0}-{{\bm{A}}}\nabla(u^{*}_{m}-u_{0})\|_{0}^{2}+\|\operatorname{div}({{\bm{\phi}}}^{*}_{m}-{{\bm{\phi}}}_{0})-B(u^{*}_{m}-u_{0})\|_{0}^{2}+\varepsilon$
$\displaystyle\leq\beta\|{{\bm{q}}}^{*}_{m}-{{\bm{q}}}_{0}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}+\varepsilon\leq\beta(d({{\bm{q}}}_{0},\mathcal{A}_{m})+\varepsilon)+\varepsilon=(2\beta+1)\varepsilon.$

Finally, by combining this estimate with the lower bound in (3.5), and
exploiting the fact that ${{\bm{q}}}_{0}$ satisfies (1.2) a.e. in $\Omega$, we
reach the estimate

$\displaystyle\|{{\bm{q}}}_{m}-{{\bm{q}}}_{0}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}$
$\displaystyle\leq\frac{1}{\alpha}\Big{(}\|{{\bm{\phi}}}_{m}-{{\bm{\phi}}}_{0}-{{\bm{A}}}\nabla(u_{m}-u_{0})\|_{0}^{2}+\|\operatorname{div}({{\bm{\phi}}}_{m}-{{\bm{\phi}}}_{0})+B(u_{m}-u_{0})\|_{0}^{2}\Big{)}$
$\displaystyle\leq\frac{1}{\alpha}\Big{(}\|{{\bm{\phi}}}_{m}-{{\bm{A}}}\nabla
u_{m}\|_{0}^{2}+\|\operatorname{div}({{\bm{\phi}}}_{m})+Bu_{m}+f\|_{0}^{2}\Big{)}=\frac{\mathcal{L}({{\bm{q}}}_{m})}{\alpha}<\frac{(2\beta+1)}{\alpha}\varepsilon,$

for every ${{\bm{q}}}_{m}\in\mathcal{I}_{m}$ and $m>m_{0}$. Since
$\varepsilon$ is arbitrary small, this concludes the proof. ∎

The result above assumes that, given ${{\bm{\Theta}}}\in\mathbb{R}^{m}$, one
can compute $\mathcal{L}(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})$
exactly. This is not the case in general, because we resort to Monte Carlo
integration for the computation of the $L^{2}$ norms in (1.3); cf. the
discrete loss functional (2.4). To deal with this issue, we consider a
regularized version of the loss functions $\mathcal{L}$ and
$\mathcal{L}_{N}:\mathcal{A}_{m}\to\overline{\mathbb{R}}$, using
$\mathbb{R}^{m}$ as domain. Given $R>0$, we define the regularized functional
$L:\mathbb{R}^{m}\to\mathbb{R}$ as

(3.6)
$L({{\bm{\Theta}}}):=\left\\{\begin{array}[]{ll}\mathcal{L}(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})&\mbox{if
}|{{\bm{\Theta}}}|\leq R,\\\ +\infty&\mbox{otherwise. }\\\ \end{array}\right.$

Next, we let $\\{X_{i}\\}_{i\in\mathbb{N}}$ be an i.i.d. sequence of random
variables, defined on a probability space $(\Lambda,\Sigma,P)$ with
$X_{i}:\Lambda\to\Omega\quad\forall i\in\mathbb{N}$, with uniform probability
density on $\Omega$. Given $\lambda\in\Lambda$, $R>0$, and $N\in\mathbb{N}$ we
set $V_{N}(\lambda):=\cup_{i\leq N}\\{X_{i}(\lambda)\\}$, and the regularized
discrete functional $L_{\lambda,N}:\mathbb{R}^{m}\to\overline{\mathbb{R}}$ as

(3.7) $L_{\lambda,N}({{\bm{\Theta}}}):=\left\\{\begin{aligned}
&\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}G_{1}({{\bm{\Theta}}},x)+G_{2}({{\bm{\Theta}}},x)&\mbox{ if
}|{{\bm{\Theta}}}|\leq R,\\\ &+\infty&\mbox{otherwise, }\\\
\end{aligned}\right.$

with $G_{1}$ and $G_{2}$ as in (3.4).

With these definitions, we can prove the pointwise $P$-almost sure convergence
of the sequence $\\{L_{\lambda,N}\\}_{N\in\mathbb{N}}$ towards $L$.

######  Lemma 3.3 (almost sure convergence of regularized discrete loss
functions).

Consider $R>0$, $L$ as in (3.6), $L_{\lambda,N}$ and
$\\{X_{i}\\}_{i\in\mathbb{N}}$ an i.i.d. family of random variables defined in
the probability space $(\Lambda,\Sigma,P)$ as in (3.7). Then
$L_{\lambda,N}({{\bm{\Theta}}})\to L({{\bm{\Theta}}})$ as $N\to\infty$
$P$-almost surely, for all ${{\bm{\Theta}}}\in\mathbb{R}^{m}$.

###### Proof.

Since we are using the same parameter $R$ in the definitions of $L$ and
$L_{\lambda,N}$, if $|{{\bm{\Theta}}}|>R$ we have
$L({{\bm{\Theta}}})=L_{\lambda,N}({{\bm{\Theta}}})=+\infty$ and there is
nothing to be proven. We therefore assume $|{{\bm{\Theta}}}|\leq R$. Recalling
$V_{N}(\lambda)=\cup_{i\leq N}\\{X_{i}(\lambda)\\}$ with $\lambda\in\Lambda$
and the definition of $G_{1}$ and $G_{2}$ (3.4), an application of the strong
law of large numbers yields

$\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}|{{\bm{\phi}}}(x)-{{\bm{A}}}\nabla
u(x)|^{2}\xrightarrow[N\to\infty]{a.s.}\int_{\Omega}|{{\bm{\phi}}}-{{\bm{A}}}\nabla
u|^{2},$

and

$\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}|\operatorname{div}{{\bm{\phi}}}(x)-Bu(x)+f(x)|^{2}\xrightarrow[N\to\infty]{a.s.}\int_{\Omega}|\operatorname{div}{{\bm{\phi}}}-Bu+f|^{2}$

for all $(u,{{\bm{\phi}}})\in\mathcal{A}_{m}$. It follows immdiately that
$L_{\lambda,N}({{\bm{\Theta}}})\to L({{\bm{\Theta}}})$ $P$-almost surely as
$N\to\infty$. ∎

We are now in position to prove the almost sure $\Gamma$-convergence of
$L_{\lambda,N}$ to $L$ as the number of quadrature points $N\to\infty$.

######  Theorem 3.2 (almost sure $\Gamma$-convergence).

Let $R>0$, $L$ be as in (3.6), and $L_{\lambda,N}$ and
$\\{X_{i}\\}_{i\in\mathbb{N}}$ be an i.i.d. family of random variables defined
in the probability space $(\Lambda,\Sigma,P)$ as in (3.7). Then, assuming
$f\in L^{2}(\Omega)$, it holds that $L_{\lambda,N}\xrightarrow[]{\Gamma}L$ as
$N\to\infty$ $P$-almost surely.

###### Proof.

We first observe that the lim-sup inequality is a trivial corollary of Lemma
3.3. Indeed, it suffices to consider the recovery sequence
$\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}\subset\mathbb{R}^{m}$,
${{\bm{\Theta}}}_{N}\equiv{{\bm{\Theta}}}$, and by Lemma 3.3 we have
$L_{\lambda,N}({{\bm{\Theta}}}_{N})\to L({{\bm{\Theta}}})$ with $N\to\infty$
$P$-almost surely.

We next prove the lim-inf inequality. Given
${{\bm{\Theta}}}\in\mathbb{R}^{m}$, let
$\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}\subset\mathbb{R}^{m}$ be a
sequence of parameters such that ${{\bm{\Theta}}}_{N}\to{{\bm{\Theta}}}$. We
aim to prove that

(3.8)
$L({{\bm{\Theta}}})\leq\liminf_{N\to\infty}L_{\lambda,N}({{\bm{\Theta}}}_{N}).$

We observe that, if $|{{\bm{\Theta}}}|>R$ then there exists
$N_{0}=N_{0}(\lambda)$ such that
$L({{\bm{\Theta}}})=L_{\lambda,N}({{\bm{\Theta}}}_{N})=+\infty$ for all
$N>N_{0}$, and (3.8) trivially holds. Therefore, without loss of generality we
assume $\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}\subset\overline{B(0,R)}$.
In that case, we extract a subsequence in such a way that
$L_{\lambda,N}({{\bm{\Theta}}}_{N})\to\liminf_{N\to\infty}L_{\lambda,N}({{\bm{\Theta}}}_{N})$
and, for the sake of simplicity, we omit the relabeling. By Lemma 3.1, the map
${{\bm{\Theta}}}\mapsto(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})\in(\mathcal{A}_{m},\|\cdot\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)})$ is continuous and therefore
$(u_{{{\bm{\Theta}}}_{N}},{{\bm{\phi}}}_{{{\bm{\Theta}}}_{N}})\to(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})$
in the $H^{1}(\Omega)\times H(\operatorname{div};\Omega)$ norm. Because
$\Omega$ is bounded, this implies

$\displaystyle\|u_{{{\bm{\Theta}}}_{N}}-u_{{{\bm{\Theta}}}}\|_{L^{1}(\Omega)}\to
0,$ $\displaystyle\quad\|\nabla u_{{{\bm{\Theta}}}_{N}}-\nabla
u_{{{\bm{\Theta}}}}\|_{L^{1}(\Omega)}\to 0,$
$\displaystyle\|{{\bm{\phi}}}_{{{\bm{\Theta}}}_{N}}-{{\bm{\phi}}}_{{{\bm{\Theta}}}}\|_{L^{1}(\Omega)}\to
0,\quad\mbox{and}$ $\displaystyle\
\|\operatorname{div}{{\bm{\phi}}}_{{{\bm{\Theta}}}_{N}}-\operatorname{div}{{\bm{\phi}}}_{{{\bm{\Theta}}}}\|_{L^{1}(\Omega)}\to
0.$

Then, defining $G_{1}$ and $G_{2}$ as in (3.4), we extract another subsequence
in such a way that
$G_{1}({{\bm{\Theta}}}_{N},x)+G_{2}({{\bm{\Theta}}}_{N},x)\to
G_{1}({{\bm{\Theta}}},x)+G_{2}({{\bm{\Theta}}},x)$ almost everywhere in
$\Omega$, and, as before, we omit the relabeling.

In order to prove (3.8), we are going to show that the latter subsequence
satisfies $L_{\lambda,N}({{\bm{\Theta}}}_{N})\to L({{\bm{\Theta}}})$ with
$N\to\infty$ $P$-almost surely. Let $\varepsilon>0$ be an arbitrary number,
using the triangle inequality, we split

(3.9)
$|L_{\lambda,N}({{\bm{\Theta}}}_{N})-L({{\bm{\Theta}}})|\leq|L_{\lambda,N}({{\bm{\Theta}}}_{N})-L_{\lambda,N}({{\bm{\Theta}}})|+|L_{\lambda,N}({{\bm{\Theta}}})-L({{\bm{\Theta}}})|.$

From Lemma 3.3, it follows that
$|L_{\lambda,N}({{\bm{\Theta}}})-L({{\bm{\Theta}}})|\to 0$ $P$-almost surely.
Thus, there exists $N_{0}=N_{0}(\lambda)$ such that
$|L_{\lambda,N}({{\bm{\Theta}}})-L({{\bm{\Theta}}})|\leq\varepsilon/4$ for all
$N>N_{0}$.

In order to bound the first term in the right hand side in (3.9), we first
observe that Lemma 3.1 shows that $G_{1}$ is uniformly bounded and $G_{2}$ is
bounded above by some integrable function. Thus, there exists $s\in
L^{1}(\Omega)$, depending on $R$, such that

(3.10)
$\big{|}G_{1}({{\bm{\Theta}}}_{N},x)+G_{2}({{\bm{\Theta}}}_{N},x)-G_{1}({{\bm{\Theta}}},x)-G_{2}({{\bm{\Theta}}},x)\big{|}\leq
s(x),$

for all $({{\bm{\Theta}}},x)\in B(0,R)\times\Omega$. Now we apply Egorov’s
Theorem to construct a set $\mathcal{K}\subset\Omega$ such that
$\int_{\mathcal{K}}s(x)dx<\varepsilon/8$ and
$G_{1}({{\bm{\Theta}}}_{N},\cdot)+G_{2}({{\bm{\Theta}}}_{N},\cdot)\to
G_{1}({{\bm{\Theta}}},\cdot)+G_{2}({{\bm{\Theta}}},\cdot)$ uniformly in
$\Omega\setminus\mathcal{K}$. We bound

$|L_{\lambda,N}({{\bm{\Theta}}}_{N})-L_{\lambda,N}({{\bm{\Theta}}})|\leq
A_{1}+A_{2},$

where

$\displaystyle A_{1}$ $\displaystyle=\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)\cap(\Omega\setminus\mathcal{K})}\big{|}G_{1}({{\bm{\Theta}}}_{N},x)+G_{2}({{\bm{\Theta}}}_{N},x)-G_{1}({{\bm{\Theta}}},x)-G_{2}({{\bm{\Theta}}},x)\big{|},$
$\displaystyle A_{2}$ $\displaystyle=\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)\cap\mathcal{K}}\big{|}G_{1}({{\bm{\Theta}}}_{N},x)+G_{2}({{\bm{\Theta}}}_{N},x)-G_{1}({{\bm{\Theta}}},x)-G_{2}({{\bm{\Theta}}},x)\big{|}.$

Using the uniform convergence in $\Omega\setminus\mathcal{K}$, $P$-almost
surely there exists $N_{1}=N_{1}(\lambda)$ such that, if $N>N_{1}$, then
$\big{|}G_{1}({{\bm{\Theta}}}_{N},x)+G_{2}({{\bm{\Theta}}}_{N},x)-G_{1}({{\bm{\Theta}}},x)-G_{2}({{\bm{\Theta}}},x)\big{|}<\frac{\varepsilon}{4|\Omega|}$
for all $x\in\Omega\setminus\mathcal{K}$. Then, it follows that
$A_{1}<\varepsilon/4$ if $N>N_{1}$.

On the other hand, we use (3.10) to derive

$A_{2}\leq\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}\chi_{\mathcal{K}}(x)s(x).$

By the strong law of large numbers, we have

$\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}\chi_{\mathcal{K}}(x)s(x)\xrightarrow[N\to\infty]{a.s.}\int_{\mathcal{K}}s(x)<\frac{\varepsilon}{8}.$

Therefore, $P$-almost surely there exists $N_{2}=N_{2}(\lambda)$ such that, if
$N>N_{2}$ then

$\Big{|}\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}\chi_{\mathcal{K}}(x)-\int_{\mathcal{K}}s(x)\Big{|}<\frac{\varepsilon}{8},$

which implies that $\frac{|\Omega|}{N}\sum_{x\in
V_{N}(\lambda)}\chi_{\mathcal{K}}(x)s(x)<\frac{\varepsilon}{4}$. Consequently,
we have $A_{2}<\frac{\varepsilon}{4}.$

Collecting the estimates above, it follows that $P$-almost surely we can
choose $N^{\prime}=N^{\prime}(\lambda)=\max\\{N_{0},N_{1},N_{2}\\}$ such that

$|L_{\lambda,N}({{\bm{\Theta}}}_{N})-L({{\bm{\Theta}}})|\leq|L_{\lambda,N}({{\bm{\Theta}}}_{N})-L_{\lambda,N}({{\bm{\Theta}}})|+|L_{\lambda,N}({{\bm{\Theta}}})-L({{\bm{\Theta}}})|\leq\varepsilon,$

for all $N>N^{\prime}$. This shows that (3.8) holds, and concludes the proof.
∎

The following theorem is the main result of this section and it roughly states
that, if we have a reasonable procedure for the minimization of
$L_{\lambda,N}$ on $\mathcal{A}_{m}$, then we can expect convergence to the
solution ${{\bm{q}}}_{0}$.

######  Theorem 3.3 (convergence).

Suppose that for any fixed $m\in\mathbb{N}$ and $R>0$ we can construct a
sequence $\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}\subset
B(0,R)\subset\mathbb{R}^{m}$ such that
$\lim_{N\to\infty}L_{\lambda,N}({{\bm{\Theta}}}_{N})=\lim_{N\to\infty}\inf_{{{\bm{\Theta}}}\in\mathbb{R}^{m}}L_{\lambda,N}({{\bm{\Theta}}})$,
with $L_{\lambda,N}$ defined as in (3.7). Let
$(u_{0},{{\bm{\phi}}}_{0})={{\bm{q}}}_{0}=\arg\min_{{{\bm{q}}}\in\mathcal{A}}\mathcal{L}({{\bm{q}}})$.
Given $\varepsilon>0$, there $P$-almost surely exist
$m_{0}=m_{0}(\varepsilon)\in\mathbb{N}$, $R=R(m_{0})>0$ and
$N_{0}=N_{0}(m_{0})\in\mathbb{N}$ such that, if one constructs a sequence
$\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}$ as above, then

$\|(u_{0},{{\bm{\phi}}}_{0})-(u_{{{\bm{\Theta}}}_{N}},{{\bm{\phi}}}_{{{\bm{\Theta}}}_{N}})\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}\leq\varepsilon\quad\mbox{for all }N>N_{0},$

where $(u_{{{\bm{\Theta}}}_{N}},{{\bm{\phi}}}_{{{\bm{\Theta}}}_{N}})$ is the
neural network function defined by the parameters ${{\bm{\Theta}}}_{N}$.

###### Proof.

Let $\varepsilon>0$ and consider the set of neural network quasi-minimizers
introduced in Lemma 3.2,
$\mathcal{I}_{m}=\\{{{\bm{q}}}\in\mathcal{A}_{m}:\mathcal{L}({{\bm{q}}})\leq\mathcal{L}({{\bm{q}}}^{*})+1/m\
\forall{{\bm{q}}}^{*}\in\mathcal{A}_{m}\\}$. By that lemma, there exists
$m_{0}>0$ such that

(3.11) $\|{{\bm{q}}}_{0}-{{\bm{q}}}_{m_{0}}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}<\varepsilon/2,$

for all ${{\bm{q}}}_{m_{0}}\in\mathcal{I}_{m_{0}}$. Next, we fix $R_{0}>0$
large enough so that there exists ${{\bm{\Theta}}}\in B(0,R_{0})$ with
${{\bm{q}}}_{{\bm{\Theta}}}^{*}=(u_{{{\bm{\Theta}}}}^{*},{{\bm{\phi}}}_{{{\bm{\Theta}}}}^{*})\in\mathcal{I}_{m_{0}}$.
For the functional $L$ defined in (3.6), this implies that
$q_{{{\bm{\Theta}}}}\in\mathcal{I}_{m}$ for all
${{\bm{\Theta}}}\in\arg\min_{{{\bm{\Theta}}}\in
B(0,R_{0})}L({{\bm{\Theta}}})$.

For this choice of $m_{0}$ and $R_{0}$, from Theorem 3.2 we have
$L_{\lambda,N}\xrightarrow[]{\Gamma}L$ P-almost surely. From the definition of
$L_{\lambda,N}$ (3.7), it follows immediately that
$\\{L_{\lambda,N}\\}_{N\in\mathbb{N}}$ is an equi-coercive sequence, according
to Definition 3.2. Therefore, we deduce that $P$-almost surely there exists
$N_{0}>0$ such that

(3.12)
$\|(u_{{{\bm{\Theta}}}_{N}},{{\bm{\phi}}}_{{{\bm{\Theta}}}_{N}})-{{\bm{q}}}_{m_{0}}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)}<\varepsilon/2$

for all $N>N_{0}$ for some ${{\bm{q}}}_{m_{0}}\in\mathcal{I}_{m_{0}}$. This
bound follows by Theorem 3.1 because every cluster point of
$\\{{{\bm{\Theta}}}_{N}\\}$ is a minimum point for $L$, and because of the
continuity of the map
${{\bm{\Theta}}}\mapsto(u_{{{\bm{\Theta}}}},{{\bm{\phi}}}_{{{\bm{\Theta}}}})$.

The proof concludes upon combining (3.11) and (3.12). ∎

##  4\. General framework

In this section, we extend the theoretical analysis we performed in Section 3
and put it into an abstract framework. Afterwards, we illustrate how such a
framework applies to some well-established unstructured neural-network methods
for the approximation of PDEs.

Let $\Omega\subset\mathbb{R}^{d}$ and $\gamma\in\mathbb{N}$. We assume our
problem is posed in some admissible vector space

$\mathcal{A}\subset W^{\gamma,1}_{loc}(\Omega;\mathbb{R}^{n}),$

namely, that every function ${{\bm{q}}}\in\mathcal{A}$ has locally integrable
weak derivatives of order up to $\gamma$. The space $\mathcal{A}$ may or may
not include boundary conditions or constraints of any type. In the setting we
described in Section 1, the target dimension is $n=1+d$, the differentiability
index is $\gamma=1$, and we identify
$\mathcal{A}\ni{{\bm{q}}}=(u,{{\bm{\phi}}})$. Additionally, we assume the
space $\mathcal{A}$ is furnished with some norm $\|\cdot\|_{\mathcal{A}}$,
which in our setting corresponds to the $H^{1}(\Omega)\times
H(\operatorname{div};\Omega)$-norm.

We consider $\omega_{1},...,\omega_{K}$ Borel subsets of $\overline{\Omega}$,
each $\omega_{i}$ furnished with a finite Radon measure $\mu_{i}$, and some
given functions $f_{1},...,f_{n_{f}}$ with $f_{i}:\Omega\to\mathbb{R}$. Given
some integrable functions
$F_{i}:\mathbb{R}^{n_{\gamma}n+n_{f}+d}\to\mathbb{R}$, $1\leq i\leq K$, we
define the loss functional

$\mathcal{L}({{\bm{q}}}):=\sum^{K}_{i=1}\int_{\omega_{i}}F_{i}(D^{\alpha_{1}}{{\bm{q}}},\ldots,D^{\alpha_{n_{\gamma}}}{{\bm{q}}},f_{1},\ldots,f_{n_{f}},x)\,d\mu_{i},$

with $F_{i}$ in such a way that all the integrals involved are well defined.
Namely, we assume the loss functional consists of $K$ terms, each of which may
be defined on different subdomains of $\overline{\Omega}$. Each of these terms
involves certain partial derivatives of ${{\bm{q}}}$ of order up to $\gamma$.
The subdomains $\omega_{i}$ need not be open; for example, we could allow for
$\omega_{i}\subset\partial\Omega$ and the corresponding term would be able to
accommodate boundary data. In such a case, the corresponding trace operator
must be bounded on the space $\mathcal{A}$.

Consider now a space $\mathcal{A}_{m}\subset\mathcal{A}$ in such a way that we
have a surjective map ${{\bm{\Theta}}}:\mathbb{R}^{m}\mapsto\mathcal{A}_{m}$.
In the setting from Section 3, this space consists of the functions obtained
through a neural network with a modification to account for boundary
conditions, cf. (3.1). We denote by ${{\bm{q}}}_{{{\bm{\Theta}}}}$ a generic
element of $\mathcal{A}_{m}$. For $1\leq i\leq K$, we define
$G_{i}({{\bm{\Theta}}},x):\mathbb{R}^{n}\times\Omega\to\mathbb{R}$ as

$G_{i}({{\bm{\Theta}}},x)=F_{i}(D^{\alpha_{1}}{{\bm{q}}}_{{{\bm{\Theta}}}},\ldots,D^{\alpha_{n_{\gamma}}}{{\bm{q}}}_{{{\bm{\Theta}}}},f_{1},\ldots,f_{n_{f}},x)$

and, given $R>0$, we define the regularized loss functional
$L:\mathbb{R}^{m}\to\mathbb{R}$

(4.1)
$L({{\bm{\Theta}}}):=\left\\{\begin{array}[]{ll}\mathcal{L}({{\bm{q}}}_{{{\bm{\Theta}}}})&\mbox{if
}|{{\bm{\Theta}}}|\leq R,\\\ +\infty&\mbox{otherwise}.\\\ \end{array}\right.$

Let $\\{X^{1}_{j}\\}_{j\in\mathbb{N}},...,\\{X^{K}_{j}\\}_{j\in\mathbb{N}}$
i.i.d. sequences of random variables, defined in the probability space
$(\Lambda,\Sigma,P)$ with $X^{i}_{j}:\Lambda\to\omega_{i}\quad\forall
j\in\mathbb{N}$, $1\leq i\leq K$, in such a way that the probability density
$\overline{\mu}^{i}$ of $X^{i}_{j}$ is distributed as $\mu_{i}$ on
$\omega_{i}$, that is

$\overline{\mu}^{i}(E)=\frac{\mu_{i}(E)}{\mu_{i}(\omega_{i})}\quad\mbox{for
every Borel set }E\subset\omega_{i}.$

Given $\lambda\in\Lambda$, $R>0$, and $N\in\mathbb{N}$ we define the sampling
nodes $V^{i}_{N}(\lambda):=\cup_{j\leq N}\\{X^{i}_{j}(\lambda)\\}$, and the
regularized discrete loss functional
$L_{\lambda,N}:\mathbb{R}^{m}\to\mathbb{R}$,

(4.2) $L_{\lambda,N}({{\bm{\Theta}}}):=\left\\{\begin{aligned}
&\sum^{K}_{i=1}\frac{\mu_{i}(\omega_{i})}{N}\sum_{x\in
V^{i}_{N}(\lambda)}G_{i}({{\bm{\Theta}}},x)&\mbox{if }|{{\bm{\Theta}}}|\leq
R,\\\ &+\infty&\mbox{otherwise.}\\\ \end{aligned}\right.$

In order to extend our convergence estimates in Section 3 to a general
framework, we consider the following hypotheses:

  * (H1)

The map $\mathbb{R}^{m}\mapsto(\mathcal{A}_{m},\|\cdot\|_{\mathcal{A}})$ with
${{\bm{\Theta}}}\mapsto{{\bm{q}}}_{{\bm{\Theta}}}$ is continuous.

  * (H2)

For all $1\leq i\leq K$ and every convergent sequence
$\\{{{\bm{q}}}_{{{\bm{\Theta}}}_{n}}\\}_{n\in\mathbb{N}}\subset\mathcal{A}_{m}$,
with
${{\bm{q}}}_{{{\bm{\Theta}}}_{n}}\to{{\bm{q}}}_{{{\bm{\Theta}}}}\in\mathcal{A}_{m}$
with respect to the $\mathcal{A}$-norm, there exists a subsequence
$\\{{{\bm{q}}}_{{{\bm{\Theta}}}_{n_{j}}}\\}_{j\in\mathbb{N}}$ such that
$G_{i}({{\bm{\Theta}}}_{n_{j}},x)\to G_{i}({{\bm{\Theta}}},x)$
$\mu_{i}$-almost everywhere.

  * (H3)

For every $R>0$, there exist functions $s_{i}\in L^{1}_{\mu_{i}}(\omega_{i})$
such that $|G_{i}({{\bm{\Theta}}},x)|\leq s_{i}(x)$ for all $1\leq i\leq K$,
for all ${{\bm{\Theta}}}\in B(0,R)$, and $\mu_{i}$-almost every
$x\in\omega_{i}$.

  * (H4)

The loss function $\mathcal{L}$ has a unique minimizer
${{\bm{q}}}_{0}\in\mathcal{A}$.

  * (H5)

Let
$\mathcal{I}_{m}:=\\{{{\bm{q}}}\in\mathcal{A}_{m}:\mathcal{L}({{\bm{q}}})\leq\mathcal{L}({{\bm{q}}}^{*})+1/m\
\forall{{\bm{q}}}^{*}\in\mathcal{A}_{m}\\}$ be the set of neural network
quasi-minimizers. Then,
$\sup_{{{\bm{q}}}_{m}\in\mathcal{I}_{m}}\|{{\bm{q}}}_{m}-{{\bm{q}}}_{0}\|_{\mathcal{A}}\to
0$ as $m\to\infty$.

Let us comment on these assumptions and how they relate to our analysis in the
previous section. Hypothesis (H1) corresponds to the first part in the
conclusion of Lemma 3.1, and guarantees the stability of neural network
functions with respect to the parameters. Hypothesis (H2) roughly states that,
for neural network functions, one can pass from convergence in $\mathcal{A}$
to almost everywhere convergence (up to a subsequence). In our setting, we
showed this condition to hold in the proof of Theorem 3.2. Our assumption (H3)
requires the existence of an $L^{1}$-upper bound for the terms $G_{i}$. This
condition appeared in the second part of Lemma 3.1. The ellipticity of the
loss functional $\mathcal{L}$ guarantees that hypothesis (H4) is satisfied.
Finally, hypothesis (H5) involves the approximability of the solution to the
continuous problem by the neural network quasi-minimizers of $\mathcal{L}$. In
our setting, this appeared in Lemma 3.2, and is a consequence of ellipticity
and assumption (3.2). Clearly, in the definition of the set $\mathcal{I}_{m}$,
one can equivalently use any positive sequence $a_{m}\searrow 0$ instead of
our default choice $a_{m}=1/m$.

The following two results extend Theorem 3.2 and Theorem 3.3, respectively; we
outline the main steps of their proofs. We first address the
$\Gamma$-convergence of the regularized discrete functionals.

######  Theorem 4.1 (almost sure $\Gamma$-convergence, general case).

Let $R>0$, and $L$, $L_{\lambda,N}$ be as in (4.1) and (4.2), respectively.
Then, under assumptions (H1), (H2), and (H3), it holds that
$L_{\lambda,N}\xrightarrow[]{\Gamma}L$ with $N\to\infty$ P-almost surely.

###### Proof.

The arguments used in the proof of Theorem 3.2 can be easily adapted to this
case. Indeed, the lim-sup inequality follows trivially by taking the recovery
sequence $\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}$,
${{\bm{\Theta}}}_{N}\equiv{{\bm{\Theta}}}$ and using a strong law of large
numbers.

To prove the lim-inf inequality, we start from a bounded sequence of
parameters $\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}$ and use (H1) to
extract a converging subsequence
$\\{{{\bm{q}}}_{{{\bm{\Theta}}}_{N}}\\}_{N\in\mathbb{N}}$ in the
$\mathcal{A}$-norm. Then, by (H2) we can extract another subsequence such that
$G_{i}({{\bm{\Theta}}}_{n_{j}},x)\to G_{i}({{\bm{\Theta}}},x)$
$\mu_{i}$-almost everywhere for all $1\leq i\leq K$ and by (H3) we know that
every function $G_{i}$ has an upper bound in $L_{\mu_{i}}^{1}(\Omega)$. The
conclusion then follows by applying Egorov’s Theorem on every subset
$\omega_{1},\ldots\omega_{K}$. ∎

Once we have the almost sure $\Gamma$-convergence of the regularized discrete
functionals, the convergence of the neural network minimizers can be proved by
arguing as in Theorem 3.3.

######  Theorem 4.2 (convergence, general case).

Assume hypotheses (H1)–(H5) are satisfied, and suppose that for any fixed
$m\in\mathbb{N}$ and $R>0$ we can construct a sequence
$\\{{{\bm{\Theta}}}_{N}\\}_{N\in\mathbb{N}}\subset
B(0,R)\subset\mathbb{R}^{m}$ such that
$\lim_{N\to\infty}L_{\lambda,N}({{\bm{\Theta}}}_{N})=\lim_{N\to\infty}\inf_{{{\bm{\Theta}}}\in\mathbb{R}^{m}}L_{\lambda,N}({{\bm{\Theta}}})$,
with $L_{\lambda,N}$ defined as in (3.7). Let
${{\bm{q}}}_{0}=\arg\min_{{{\bm{q}}}\in\mathcal{A}}\mathcal{L}({{\bm{q}}})$.
Given $\varepsilon>0$ there exist $m_{0}=m_{0}(\varepsilon)\in\mathbb{N}$,
$R=R(m_{0})>0$ and $N_{0}=N_{0}(m_{0})\in\mathbb{N}$ P-almost surely, such
that

$\|{{\bm{q}}}_{0}-{{\bm{q}}}_{{{\bm{\Theta}}}_{N}}\|_{\mathcal{A}}\leq\varepsilon\quad\mbox{for
all }N>N_{0},$

where ${{\bm{q}}}_{{{\bm{\Theta}}}_{N}}\in\mathcal{A}_{m_{0}}$ is the neural
network function defined by the parameters ${{\bm{\Theta}}}_{N}$.

###### Proof.

We first remark that hypothesis (H4) is needed to guarantee the existence of a
well-defined minimizer ${{\bm{q}}}_{0}\in\mathcal{A}$, and therefore (H5) is
meaningful. Given $\varepsilon>0$, we use hypothesis (H5) to find
$m_{0}\in\mathbb{N}$ such that, if ${{\bm{q}}}_{m_{0}}\in\mathcal{I}_{m_{0}}$
then $\|{{\bm{q}}}_{0}-{{\bm{q}}}_{m_{0}}\|_{\mathcal{A}}<\varepsilon/2$.

Next, we fix $R_{0}>0$ large enough so that there exists ${{\bm{\Theta}}}\in
B(0,R_{0})$ with ${{\bm{q}}}_{{{\bm{\Theta}}}}\in\mathcal{I}_{m_{0}}$, and use
this $R_{0}$ in Theorem 4.1 to deduce that
$L_{\lambda,N}\xrightarrow[]{\Gamma}L$ with $N\to\infty$ P-almost surely. The
result then follows by the equi-coercivity of the sequence
$\\{L_{\lambda,N}\\}_{N\in\mathbb{N}}$ by applying the fundamental theorem of
$\Gamma$-convergence (Theorem 3.1). ∎

We next discuss how two well-known methods fit into the framework in
hypotheses (H1)–(H5), and thus Theorem 4.2 establishes their convergence.

######  Remark 2 (Deep Ritz Method).

The DRM was proposed by E and Yu in [9], and is tailored for numerically
solving variational problems. A prototypical example is the homogeneous
Dirichlet problem, that corresponds to the minimization of the energy
$\mathcal{L}\colon H^{1}_{0}(\Omega)\to\mathbb{R}$,

$\mathcal{L}(u)=\frac{1}{2}\int_{\Omega}|\nabla u|^{2}-\int_{\Omega}fu.$

We assume $\|f\|_{L^{2}(\Omega)}<\infty$, consider
$\mathcal{A}=H^{1}_{0}(\Omega)$, and define the neural network spaces
$\mathcal{A}_{m}$ as in (3.1). Arguing as in Section 3, it is possible to show
that hypotheses (H1)–(H4) hold for this loss function. Indeed, (H1) and (H3)
can be proved in the same fashion as Lemma 3.1, while (H2) follows because for
every bounded sequence in $H^{1}_{0}(\Omega)$ we can extract an almost
everywhere convergent subsequence, and (H4) is a standard PDE result. Finally,
hypothesis (H5) can be obtained from classical approximation results [7, 13,
2, 32, 12].

######  Remark 3 (Deep Galerkin Method).

The DGM was introduced by Sirignano and Spiliopoulos in [28], and uses as loss
functional the $L^{2}$-norm of the PDE residual on the neural network
functions. Within the convergence framework in [28, Section 7], and the
conditions assumed there, we set
$\mathcal{A}:=\mathcal{C}^{0,\delta,\delta/2}(\overline{\Omega}_{T})\cap
L^{2}((0,T];W^{1,2}_{0}(\Omega))\cap W^{(1,2),2}_{0}(\Omega_{T}^{\prime})$,
where $\delta>0$ and $\Omega_{T}^{\prime}$ is any interior subdomain of
$\Omega_{T}$, cf. Theorem 7.3. We furnish this space with the
$\|\cdot\|_{H^{2}(\Omega_{T})}$ norm, and define $\mathcal{A}_{m}$ according
to (3.1).

Then, assumptions (H1) and (H3) can be verified by arguing as in Lemma 3.1 by
requiring suitable regularity assumptions on the initial and boundary data and
parameters of the equation; for example, these hold straightforwardly for
these data and parameters are bounded. Hypothesis (H2) can be proved by using
the boundedness of $\Omega_{T}$ and arguing as in the proof of Theorem 3.2 to
exploit the convergence properties of the $H^{2}(\Omega)$-norm. Finally,
hypotheses (H4) and (H5) are addressed in [28, Theorem 7.3]. We observe that
despite this result is stated for a single minimizing sequence $\\{f^{n}\\}$,
defined in [28, Theorem 7.1], the arguments applies ”uniformly” to any
possible construction of $\\{f^{n}\\}$, and then (H5) is verified. Finally, we
point out that the convergence of discrete minimizers of $\mathcal{L}$ is
proven in the weaker norm $\|\cdot\|_{L^{\rho}(\Omega_{T})}$, with $\rho<2$.
Therefore, our conclusion in Theorem 4.2 is valid if we measure convergence in
such a norm.

##  5\. Numerical experiments

In this section, we present numerical results for the method we proposed in
Section 2. We did not prioritize any particular neural network architecture,
and used between one- and five-layer networks with sigmoidal activation
functions to construct $u_{{{\bm{\Theta}}}}$ and
${{\bm{\phi}}}_{{{\bm{\Theta}}}}$. For the construction of the auxiliary
functions ${{\bm{n}}}$, $d_{\mathcal{D}}$, $d_{\mathcal{N}}$,
$G_{\mathcal{D}}$, $G_{\mathcal{N}}$, we used between one and three-layer
networks with less neurons per layer. In the training process, we used the
ADAM [14] algorithm to update the parameters, with a decaying learning rate
schedule.

We observe an improvement in the method’s performance when explicit
approximations of the auxiliary functions $d_{\mathcal{D}}$ and
$d_{\mathcal{N}}$ are used. These functions, which depend on the geometry of
the domain, are many times explicitly available in practice.

We recall that, as explained in sections 3 and 4, the numerical solution
depends on the number of degrees of freedom $m$ and the number of collocation
points $N$. Both must go to infinity to guarantee convergence. In all the
numerical examples we show below, these quantities remain fixed. Therefore, in
these examples the convergence as a function of the iterations occurs towards
the minimizer of the discrete loss functional $L_{\lambda,N}$ (cf. (3.7))
corresponding to the values of $m$ and $N$ we have set.

######  Example 5.1 (Laplace operator).

We consider the following problem in arbitrary dimension. Let
$\Omega=\\{x\in{\mathbb{R}^{d}}:\ -1<x_{1},...,x_{d}<1\\}$,
$\Gamma_{\mathcal{N}}=[-1,1]^{d-1}\times\\{1\\},$ and $k\in\mathbb{N}$. We
seek $u\colon\Omega\to\mathbb{R}$ such that

(5.1) $\left\\{\begin{aligned} -\Delta u&=\prod^{d-1}_{i=1}\sin(k\pi
x_{i})\left((d-1)k^{2}\pi^{2}(1-x_{d}^{2})+2\right)\quad&\mbox{in }\Omega,\\\
u&=0\quad&\mbox{on }\partial\Omega\setminus\Gamma_{\mathcal{N}},\\\ \nabla
u\cdot{{\bm{\nu}}}&=-2\prod^{d-1}_{i=1}\sin(k\pi x_{i})\quad&\mbox{on
}\Gamma_{\mathcal{N}}.\\\ \end{aligned}\right.$

Here, we have ${{\bm{\nu}}}=(0,...,0,1)$ on $\Gamma_{\mathcal{N}}$, and the
solution to (5.1) is

$u=\prod^{d-1}_{i=1}\sin(k\pi x_{i})(1-x^{2}_{d}).$

We point out that the parameter $k$ is a frequency that allows us to choose
how oscillatory the exact solution $u$ is. We first tested the method in a
two-dimensional domain ($d=2$). Figure 5.1 displays the results we obtained
for $k=1$ and by constructing $u_{{{\bm{\Theta}}}}$ and
${{\bm{\phi}}}_{{{\bm{\Theta}}}}$ using neural networks with $15$ sigmoidal
activation functions per layer. At the end of the stochastic gradient descent
algorithm we computed the value $L_{N}({{\bm{\Theta}}})=0.0450$. Taking into
account the ellipticity of the loss function $\mathcal{L}$, arguing as in
Lemma 3.2 we deduce

$\mathcal{L}({{\bm{q}}}_{m})\simeq\|{{\bm{q}}}_{m}-{{\bm{q}}}_{0}\|_{H^{1}(\Omega)\times
H(\operatorname{div};\Omega)},$

and therefore this quantity serves as an error estimator.

Figure 5.2 corresponds to $k=2$, and we used a similar architecture, but with
$18$ sigmoidal activation functions per layer. We observed a fast convergence
in the number of iterations, reaching $L_{N}({{\bm{\Theta}}})=1.89$ by the end
of the minimization algorithm. Finally, Figure 5.3 reports the results we
obtained in case $d=5$, $k=1$. In this case, we used networks with $25$
sigmoidal activation functions per layer and obtained
$L_{N}({{\bm{\Theta}}})=2.32$.

|
---
$u_{{{\bm{\Theta}}}}$
|
---
Loss function vs. iterations.
|
---
$|u_{{{\bm{\Theta}}}}-u|$
|
---
$\|u_{{{\bm{\Theta}}}}-u\|_{L^{2}(\Omega)}$ vs. iterations.
Figure 5.1. Top left: computational solution $u_{{{\bm{\Theta}}}}$ to (5.1) in case $k=1$ and $d=2$. In the computation, we used a learning rate $\ell=0.005$, with 2,000 collocation points, 1,500 optimization steps, and 3603 degrees of freedom (including auxiliary functions). We used one-layer networks both for the main and auxiliary functions. The panel in bottom left exhibits the pointwise discrepancy $|u-u_{{\bm{\Theta}}}|$. We also report the evolution of the loss function (top right) and the $L^{2}$ error (bottom right). |
---
$u_{{{\bm{\Theta}}}}$
|
---
Loss function vs. iterations.
|
---
$|u_{{{\bm{\Theta}}}}-u|$
|
---
$\|u_{{{\bm{\Theta}}}}-u\|_{L^{2}(\Omega)}$ vs. iterations.
Figure 5.2. Computational solution $u_{{{\bm{\Theta}}}}$ (top left), evolution of the loss function (top right), pointwise error (bottom left), and evolution of the $L^{2}$-error (bottom right) for (5.1) with $k=2$ and $d=2$. We employed five-layer networks for the main functions and three-layer networks for the auxiliary functions. We used an initial learning rate $\ell=0.005$, with 5,000 collocation points, 10,000 optimization steps, and 2901 degrees of freedom. We halved the learning rate every 2,500 optimization steps. |
---
${\left.\kern-1.2ptu_{{{\bm{\Theta}}}}\vphantom{\big{|}}\right|_{\\{x_{3},\ldots,x_{5}=0.5\\}}}$
|
---
Loss function vs. iterations.
|
---
$|{\left.\kern-1.2ptu_{{{\bm{\Theta}}}}\vphantom{\big{|}}\right|_{\\{x_{3},\ldots,x_{5}=0.5\\}}}-{\left.\kern-1.2ptu\vphantom{\big{|}}\right|_{\\{x_{3},\ldots,x_{5}=0.5\\}}}|$
|
---
MSE vs. iterations.
Figure 5.3. Slice of the solution $u_{{{\bm{\Theta}}}}$ (top left), evolution
of the loss function (top right), pointwise error (bottom left), and evolution
of the Mean Squared Error (MSE) (bottom right) for (5.1) with $k=1$ and $d=5$.
We employed five-layer networks for the main functions and three-layer
networks for the auxiliary functions. We used an initial learning rate
$\ell=0.005$, with 12,000 collocation points, 20,000 optimization steps, and
5656 degrees of freedom. We halved the learning rate every 4,000 optimization
steps. We estimated the MSE by using 5,000 random points in $\Omega$ (re-
sampled at every step

######  Example 5.2 (singularly perturbed problem).

Let $\varepsilon>0$, $\Omega=(0,1)^{2}$,
${{\bm{b}}}=(-1+2\varepsilon,-1+2\varepsilon)$, $c=2(1-\varepsilon)$, and the
function $f\colon\Omega\to\mathbb{R}$,

$f(x,y)=-\left[x-\left(\frac{1-e^{-x/\varepsilon}}{1-e^{-1/\varepsilon}}\right)+y-\left(\frac{1-e^{-y/\varepsilon}}{1-e^{-1/\varepsilon}}\right)\right]e^{x+y}.$

We consider the singularly perturbed problem: find
$u\colon\Omega\to\mathbb{R}$ such that

(5.2) $\left\\{\begin{aligned} -\varepsilon\Delta u+{{\bm{b}}}\cdot\nabla
u+cu=f\quad&\mbox{in }\Omega,\\\ u=0\quad&\mbox{on
}\partial\Omega.\end{aligned}\right.$

The exact solution to (5.2) is

$u(x,y)=\left(x-\frac{1-e^{-x/\varepsilon}}{1-e^{-1/\varepsilon}}\right)\left(y-\frac{1-e^{-y/\varepsilon}}{1-e^{-1/\varepsilon}}\right)e^{x+y}.$

Figure 5.4 exhibits our computed solutions for this example with
$\varepsilon=0.05$. In that case, we observed a fast convergence towards the
solution, reaching $L_{N}({{\bm{\Theta}}})=0.0112$, as well as a good
adaptation of the discrete solution to the boundary layers.

|
---
$u_{{{\bm{\Theta}}}}(x)$
|
---
Loss function vs. iterations.
|
---
$|u_{{{\bm{\Theta}}}}(x)-u(x)|$
|
---
$\|u_{{{\bm{\Theta}}}}-u\|_{L^{2}(\Omega)}$ vs. iterations.
Figure 5.4.  Computational solution $u_{{{\bm{\Theta}}}}$ (top left),
evolution of the loss function (top right), pointwise error (bottom left), and
evolution of the $L^{2}$-error (bottom right) for (5.2) with
$\varepsilon=0.05$. We used an initial learning rate $\ell=0.005$, with 5,000
collocation points, 10,000 optimization steps, and 3543 degrees of freedom.
Auxiliary functions have been approximated exactly.

##  6\. Concluding remarks

In this work, we have proposed a First-Order System Least Squares (FOSLS)
method based on deep learning for numerically solving second-order elliptic
PDEs. This method is meshless, which is naturally advantageous for high-
dimensional problems, but as a consequence implies that we cannot compute the
loss functions exactly. Taking into account this practical issue, we proved
the almost sure convergence of the neural network minimizers towards the PDE
solutions. We furthermore extended the theoretical framework to incorporate
other methods based on Monte Carlo quadrature.

######  Remark 4 (almost-everywhere solutions).

The convergence proofs in Sections 3 and 4 are based on the use of regularized
versions of the cost functionals and their discretizations. Regularization
consists in restricting the size of the parameters, namely, imposing that
$|{{\bm{\Theta}}}|<R$ for certain $R<\infty$. This ensures that any neural
network function with large derivatives is penalized, thereby preventing
minimizers from approximating non-smooth functions.

Far from being an artificial condition of the proof, regularization mechanisms
of this kind are necessary in the implementation to avoid convergence towards
functions that satisfy the PDE almost everywhere but are not weak solutions of
the target problem. To illustrate this point, consider the following example,
which is just (1.2) in a simplified setting: seek $u,\phi:(0,1)\to\mathbb{R}$
such that

(6.1) $\left\\{\begin{aligned} \phi-u^{\prime}&=0&\mbox{in }(0,1),\\\
\phi^{\prime}&=0&\mbox{in }(0,1),\\\ u(0)&=0,\\\ u(1)&=1.\end{aligned}\right.$

Naturally, the unique minimizer of the least-squares functional (cf. (1.3))

$\mathcal{L}(u,\phi):=\|\phi-u^{\prime}\|_{L^{2}(\Omega)}^{2}+\|\phi^{\prime}\|_{L^{2}(\Omega)}^{2}$

in the corresponding admissible set
$\mathcal{A}=\\{(u,\phi)\in[H^{1}(\Omega)]^{2}\colon u(0)=0,\ u(1)=1\\}$ is
$u(x)=x$ and $\phi(x)=1$. Let $\delta\in(0,1/2)$ be a small number, and
consider the functions

(6.2) $u_{\delta}(x)=\left\\{\begin{aligned} &0&\mbox{in }(0,1/2-\delta)\\\
&\frac{x-1/2+\delta}{2\delta}&\mbox{in }(1/2-\delta,1/2+\delta)\\\ &1&\mbox{in
}(1/2+\delta,1)\end{aligned}\right.,\qquad\phi_{\delta}(x)=\left\\{\begin{aligned}
&0&\mbox{in }(0,1/2-\delta)\\\ &\frac{1}{2\delta}&\mbox{in
}(1/2-\delta,1/2+\delta)\\\ &0&\mbox{in }(1/2+\delta,1)\end{aligned}\right..$

We notice $\phi_{\delta}=u^{\prime}_{\delta}$ a.e. in $(0,1)$ and
$\mathcal{L}(u_{\delta},\phi_{\delta})=0$, although
$(u_{\delta},\phi_{\delta})\notin\mathcal{A}$, because $\phi_{\delta}$ is not
an $H^{1}$ function.

If we utilize the discrete functional (2.1) with collocation points, and none
of these points lies in the interval $(1/2-\delta,1/2+\delta)$, then for these
two functions we would have

$\mathcal{L}_{N}(u_{\delta},\phi_{\delta})=0.$

We remark that, independently of the number of collocation points $N$, one can
always take $\delta>0$ sufficiently small such that the probability of none of
the sampling points lies in $(1/2-\delta,1/2+\delta)$ is significant.
Therefore, if our neural network is capable of producing functions
$(u_{{{\bm{\Theta}}}},\phi_{{{\bm{\Theta}}}})$ approximating
$(u_{\delta},\phi_{\delta})$ in (6.2) (cf. Figure 6.1), then during the
optimization process the descent algorithm may choose to approximate the pair
$(u,\phi)=(\chi_{(1/2,1)},0)$. This function satisfies the differential
equations in (6.1) almost everywhere, but is not a significant solution. The
issue of approximating bad solutions of this kind is mitigated by applying
classic regularization techniques that penalize large parameters, because
$|{{\bm{\Theta}}}|$ must be large in order to $u_{{\bm{\Theta}}}^{\prime}$ be
large at some portion of the domain.

Figure 6.1. In red the function $u_{\delta}$ defined in 6.2. On dashed lines
the solution of problem (6.1).

This difficulty extends to all methods based on the minimization of cost
functionals similar to (2.1), such as DGM [28] or DRM [9]. The issue stems
from the fact that the functional (2.1) is unable to distinguish between
regular solutions (belonging to a suitable Sobolev space) from any other
functions that satisfy the equation almost everywhere. As far as we know, this
problem has not been addressed in the literature, and the question of how to
develop suitable regularization techniques for these approaches remains open.

######  Remark 5 (approximation of non-smooth solutions).

There are, however, problems in which the solution presents large gradients in
regions of the domain. One can typically think of singularly perturbed
problems, such as (5.2), or singularities arising due to poor boundary
regularity, such as for the Poisson problem on an $L$-shaped domain. In those
problems, regularization can limit the approximation capabilities of the
algorithm.

For algebraic boundary singularities, if the boundary conditions are imposed
in a strong fashion, as discussed in Section 2.1.2, one could aim to modify
the rate at which the corresponding auxiliary function $d_{\mathcal{D}}$ or
$d_{\mathcal{N}}$ decreases to zero near the singularity. This could
potentially avoid $v_{{{\bm{\Theta}}}}$ having to approximate a singular
function and lead to a faster convergence. Nevertheless, this requires an a
priori knowledge about the location and behavior of the singularities of the
solution, that is not available in general. We emphasize that the theory we
developed in Section 4 does not make any regularity assumption on the PDE, and
therefore includes the case of non-smooth solutions.

## Acknowledgements

The authors thank Prof. Michael Karkulik and Roberto González for their
insightful comments on a previous version of this manuscript.

## References

  * [1] R. Arora, A. Basu, P. Mianjy, and A. Mukherjee.  Understanding deep neural networks with rectified linear units.  In International Conference on Learning Representations, 2018.
  * [2] A. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information theory, 39(3):930–945, 1993.
  * [3] J. Berg and K. Nyström.  A unified deep artificial neural network approach to partial differential equations in complex geometries.  Neurocomputing, 317:28–41, 2018.
  * [4] A. Braides.  A handbook of $\Gamma$-convergence.  In Handbook of Differential Equations: stationary partial differential equations, volume 3, pages 101–213. Elsevier, 2006.
  * [5] Z. Cai, J. Chen, M. Liu, and X. Liu.  Deep least-squares methods: An unsupervised learning-based numerical method for solving elliptic PDEs.  Journal of Computational Physics, 420:109707, 2020.
  * [6] Z. Cai, R. Lazarov, T. Manteuffel, and S. McCormick.  First-order system least squares for second-order partial differential equations: Part I.  SIAM Journal on Numerical Analysis, 31(6):1785–1799, 1994.
  * [7] G. Cybenko.  Approximation by superpositions of a sigmoidal function.  Mathematics of control, signals and systems, 2(4):303–314, 1989\.
  * [8] W. E and S. Wojtowytsch.  Some observations on high-dimensional partial differential equations with Barron data.  In Mathematical and Scientific Machine Learning, pages 253–269. PMLR, 2022.
  * [9] W. E and B. Yu.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics, 6(1):1–12, 2018.
  * [10] C. He, X. Hu, and L. Mu.  A mesh-free method using piecewise deep neural network for elliptic interface problems.  Journal of Computational and Applied Mathematics, 412:114358, 2022\.
  * [11] J. He, L. Li, and J. Xu.  Relu deep neural networks from the hierarchical basis perspective.  Computers & Mathematics with Applications, 120:105–114, 2022.
  * [12] J. He, L. Li, J. Xu, and C. Zheng.  Relu deep neural networks and linear finite elements.  J. Comput. Math., 38(3):502–527, 2020.
  * [13] K. Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural networks, 4(2):251–257, 1991.
  * [14] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In In Proceedings of the 3rd InternationalConference for Learning Representations—ICLR, pages 7–9, San Diego, CA, 2015.
  * [15] I. E. Lagaris, A. Likas, and D. I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks, 9(5):987–1000, 1998.
  * [16] I. E. Lagaris, A. C. Likas, and D. G. Papageorgiou.  Neural-network methods for boundary value problems with irregular boundaries.  IEEE Transactions on Neural Networks, 11(5):1041–1049, 2000.
  * [17] H. Lee and I. S. Kang.  Neural algorithm for solving differential equations.  Journal of Computational Physics, 91(1):110–131, 1990.
  * [18] M. Liu and Z. Cai.  Adaptive two-layer ReLU neural network: II. Ritz approximation to elliptic PDEs.  Computers & Mathematics with Applications, 113:103–116, 2022.
  * [19] M. Liu, Z. Cai, and J. Chen.  Adaptive two-layer ReLU neural network: I. best least-squares approximation.  Computers & Mathematics with Applications, 113:34–44, 2022.
  * [20] L. Lyu, Z. Zhang, M. Chen, and J. Chen.  Mim: A deep mixed residual method for solving high-order partial differential equations.  Journal of Computational Physics, 452:110930, 2022.
  * [21] A. Malek and R. S. Beidokhti.  Numerical solution for high order differential equations using a hybrid neural network—optimization method.  Applied Mathematics and Computation, 183(1):260–271, 2006.
  * [22] 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.  Journal of Computational physics, 378:686–707, 2019.
  * [23] H. Sheng and C. Yang.  PFNN: a penalty-free neural network method for solving a class of second-order boundary-value problems on complex geometries.  Journal of Computational Physics, 428:110085, 2021.
  * [24] 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.  Communications in Computational Physics, 28(5):2042–2074, 2020\.
  * [25] J. Siegel, Q. Hong, X. Jin, W. Hao, and J. Xu.  A priori analysis of stable neural network solutions to numerical PDEs.  arXiv preprint arXiv:2107.04466, 2022.
  * [26] J. W. Siegel and J. Xu.  High-order approximation rates for neural networks with ReLUk activation functions.  arXiv preprint arXiv:2012.07205, 2020.
  * [27] J. W. Siegel and J. Xu.  Sharp lower bounds on the approximation rate of shallow neural networks.  arXiv preprint arXiv:2106.14997, 2021.
  * [28] J. Sirignano and K. Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [29] Z. Wang and Z. Zhang.  A mesh-free method for interface problems using the deep learning approach.  Journal of Computational Physics, 400:108963, 2020.
  * [30] S. Wojtowytsch and W. E.  Can shallow neural networks beat the curse of dimensionality? A mean field training perspective.  IEEE Transactions on Artificial Intelligence, 1(2):121–129, 2020\.
  * [31] J. Xu.  Finite neuron method and convergence analysis.  Communications in Computational Physics, 28:1707–1745, 2020.
  * [32] D. Yarotsky.  Error bounds for approximations with deep relu networks.  Neural Networks, 94:103–114, 2017.
  * [33] Y. Zang, G. Bao, X. Ye, and H. Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  Journal of Computational Physics, page 109409, 2020.
  * [34] U. Zerbinati.  PINNs and GaLS: a priori error estimates for shallow physics informed neural networks applied to elliptic problems.  IFAC-PapersOnLine, 55(20):61–66, 2022.  10th Vienna International Conference on Mathematical Modelling MATHMOD 2022.