# Inf-Sup neural networks for high-dimensional elliptic PDE problems

Xiaokai Huo<EMAIL_ADDRESS>and  Hailiang Liu<EMAIL_ADDRESS>Department
of Mathematics, Iowa State University, Ames, IA 50014

###### Abstract.

Solving high dimensional partial differential equations (PDEs) has
historically posed a considerable challenge when utilizing conventional
numerical methods, such as those involving domain meshes. Recent advancements
in the field have seen the emergence of neural PDE solvers, leveraging deep
networks to effectively tackle high dimensional PDE problems. This study
introduces Inf-SupNet, a model-based unsupervised learning approach designed
to acquire solutions for a specific category of elliptic PDEs. The fundamental
concept behind Inf-SupNet involves incorporating the inf-sup formulation of
the underlying PDE into the loss function. The analysis reveals that the
global solution error can be bounded by the sum of three distinct errors: the
numerical integration error, the duality gap of the loss function (training
error), and the neural network approximation error for functions within
Sobolev spaces. To validate the efficacy of the proposed method, numerical
experiments conducted in high dimensions demonstrate its stability and
accuracy across various boundary conditions, as well as for both semi-linear
and nonlinear PDEs.

###### Key words and phrases:

PDE stability, minimax optimization, machine learning, neural networks, error
bounds

###### 1991 Mathematics Subject Classification:

65K10, 90C26

##  1\. Introduction

In recent years, the burgeoning field of deep neural networks (DNNs) has
propelled the widespread use of machine learning techniques in scientific
computing. Renowned for their status as universal function approximators, DNNs
have gained popularity as ansatz spaces for approximating solutions to diverse
partial differential equations(PDEs) [57, 58, 42, 59, 27, 49, 15, 2]. In
contrast to conventional numerical methods such as finite elements (FEM),
finite differences (FDM), or finite volumes (FVM), which rely on grids or
meshes, neural networks offer a mesh-free, scalable approach that holds
promise for solving high-dimensional PDEs. In scenarios where ample labeled
training data can be obtained through specific numerical methods, supervised
learning becomes a viable option for training DNNs to approximate PDE
solutions, as exemplified in [43, 44, 40, 24, 26]. However, in many contexts,
a learning framework is required to approximate PDE solutions with minimal or
potentially no data. In such situations, the underlying PDE and associated
conditions must be seamlessly integrated into an unsupervised learning
strategy [27, 42]. The method presented here falls into this category.

The idea of incorporating PDEs and associated conditions with neural networks
for inferring PDE solutions can be traced back to a 1994 paper [14]. In this
early work, neural networks were limited to one hidden layer. Subsequently, a
similar idea was explored by [27], where boundary conditions were manually
enforced through a change of variables. These initial concepts evolved into
more sophisticated approaches such as Physics-informed neural networks (PINNs)
[42] and the deep Galerkin method [53], both featuring deeper hidden layers.
Among these, PINNs have gained popularity in the field of scientific machine
learning in recent years. For a differential equation $\mathcal{A}(u)=f$
defined on the domain $\Omega$ with boundary condition $\mathcal{B}u=g$ on
$\partial\Omega$, PINN trains a neural network to minimize the following
functional when evaluated on sampling points,

${\mathcal{L}}(v)=\int_{\Omega}|\mathcal{A}(v)(x)-f|^{2}dx+\lambda\int_{\partial\Omega}|\mathcal{B}(v)(x)-g|^{2}ds,$

where $\lambda>0$ is a weight parameter. However, minimizing such a functional
may lead to a higher order PDE than the original one, posing increased
difficulty when the solution regularity is relatively low. As an improvement,
the deep Least-Squares method introduced in [8] employs a least-squares
functional based on the first order system of scalar second-order elliptic
PDEs, treating boundary conditions in a balanced manner. Various modifications
have been proposed to enhance the accuracy of PINNs or adapt them for specific
problems. Examples include the gradient-enhanced PINNs (gPINNs), which
incorporate higher order solution derivatives [58], and the Variational PINNs
(VPINNs) [22], which incorporate the weak form of the PDE. For a comprehensive
overview of PINNs and their applications across diverse fields, the reader is
directed to [7, 32, 33, 18, 45, 56].

Various advancements have been made in formulating loss functions for solving
PDE problems with deep neural networks (DNNs). The deep Ritz method (DRM) [57]
employs the energy functional of the Poisson equation and is particularly
applicable to problems grounded in a minimization principle. Variational
neural networks (VarNets) [23] utilize the variational (integral) form of
advection-diffusion equations, with test functions sampled from linear spaces.
Using lower-order derivatives with training over space-time regions, the
method using VarNets gets more effective in capturing the PDE solutions. In
the context of hyperbolic conservation laws and the pursuit of entropy
solutions, weak PINNs (wPINNs) were introduced, incorporating the Kruzkov
entropy, as analyzed in [13]. Another notable approach is the Weak Adversarial
Networks (WANs) [59], which also leverage the weak PDE formulation. However,
WANs introduce an additional adversarial network to represent the test
function, enhancing the overall robustness and efficacy of the method.

A notable limitation of weak PDE formulations, such as the Weak Adversarial
Networks (WANs), is their intrinsic lack of numerical stability, as
highlighted in [55]. Addressing this issue, the Deep Double Ritz method,
introduced in the same reference [55], adopts a nested Ritz minimization
strategy to optimize test functions and enhance stability. In this paper, we
propose a simple alternative, namely Inf-SupNet. Our approach originates from
a constraint optimization problem:

$\inf_{u\in X}\int_{\partial\Omega}|\mathcal{B}u(x)-g|^{2}ds\quad\text{subject
to}\quad\mathcal{A}u=f\quad\text{in}\quad\Omega,$

which is demonstrated to be equivalent to the aforementioned PDE problem. In
such formulation, the least square of the boundary data mismatch is the
primary loss, with the underlying PDE serving as a constraint. This constraint
is elevated to the primal-dual space as an inf-sup problem:

$\inf_{u\in X}\sup_{v\in
Y^{*}}\|\mathcal{B}u(x)-g\|_{\partial\Omega}^{2}+(\mathcal{A}u-f,v)_{\Omega},$

where $\mathcal{A}:X\to Y$ represents the underlying differential operator,
and $Y^{*}$ is the dual space of $Y$. InfSupNet employs this loss function,
leveraging the inf-sup form of the PDE problem and employing an adversarial
network for the Lagrangian multiplier $v$. Our method provides flexibility in
handling the dual variable $v$ based on the solution regularity of $u$.

Inspired by a recent work [30] combining primal-dual hybrid gradient methods
with spatial Galerkin methods for approximating conservation laws with
implicit time-discretization, our formulation differs as we introduce neural
networks to represent the trial solution and the Lagrangian multiplier. This
renders the resulting optimization problem non-quadratic in the parameter
space.

Regarding the theoretical accuracy of the InfSupNet algorithm, we derive an
error estimate based on the stability of elliptic problems. The error is shown
to be bounded by the numerical integration error, the neural network
approximation error, and the training error. The min-max structure facilitates
the quantification of the training error through the duality gap between the
loss for minimization and maximization steps.

Contribution. Our key contributions for learning solutions of PDE problems are
as follows.

  * •

Present an InfSupNet, a novel method for learning solutions of PDE problems
using deep neural networks.

  * •

Prove equivalence of the inf-sup problem to the original PDE problem.

  * •

Show the global approximation error is induced by three sources: the training
error, the sampling error, and the approximation error of neural networks; and
each error is shown to be asymptotically small under reasonable assumptions.

  * •

Demonstrate through experiments both stability and accuracy of the method for
high-dimensional elliptic problems posed on both regular and irregular
domains.

### Related work

Here we delve into further discussions on related works, encompassing neural
PDE solvers, sampling methods, and error estimation.

### Mesh Dependent Neural Methods:

Our method, along with other mesh-free approaches like PINNs, DRM, WANs,
relies on auto-differentiation to compute the empirical loss, eliminating the
need for a mesh grid. This enables obtaining solutions at new locations beyond
the sampling points directly after the network is trained. In contrast, mesh-
based methods, as seen in [20, 21, 25, 4], employ conventional numerical
techniques such as finite difference or finite element methods, heavily
relying on mesh grids for computations. While mesh-based methods permit the
use of equation parameters on grids as input for network training and solving
a range of PDEs with varying parameters, changing the mesh resolutions alter
network structures, and this will necessitate new training. In contrast, mesh-
free methods maintain the underlying network when changing the resolution of
numerical integration, making them more adaptable, especially in high
dimensions [10].

### Deep BSDE Solvers:

Establishing a connection between parabolic PDEs and backward stochastic
differential equations (BSDEs), [15] introduces a deep BSDE method to evaluate
solutions at any given space-time locations for high dimensional semi-linear
parabolic PDEs. This method is extended to handle fully nonlinear second-order
PDEs in [2].

### Operator Learning:

A novel approach involves learning infinite dimensional operators with neural
networks [31, 28]. Universal approximation theorems for operator learning have
been established [9], operator learning allows the solution of a family of
PDEs with one neural construction. This approach operates as a form of
supervised learning, requiring only data, devoid of prior knowledge about the
underlying PDE.

### Data Sampling Techniques:

Apart from utilizing deeper hidden layers for solving PDEs, as discussed
earlier, it is also illustrated that the training points can be obtained by a
random sampling of the domain rather than using a mesh, which is beneficial in
higher-dimensional problems [3, 53]. Indeed, different numerical integration
schemes are employed in all mesh-free methods. For an in-depth comparison
study of quadrature rules used in solving PDEs with deep neural networks,
please refer to [47].

### Error Bounds:

Although not as mature as the applications of neural PDE solvers in various
domains, the theoretical development has made significant strides, presenting
rigorous bounds on various error sources [51, 36, 37, 19, 35, 60]. Typically
total error is decomposed into the approximation error of neural networks,
numerical sampling error and training error, each being bounded through
different strategies. For PINNs, as identified in [37], the strategy involves
considering regularity and stability of PDE solutions and quadrature error
bounds to estimate the generalization gap between continuous and discrete
versions of the PDE residual. In deriving error bounds for InfSupNet, we
pursue a related yet distinct approach. Beyond the error on PDE solutions, we
also estimate the error of the Lagrangian multiplier. Moreover, due to the use
of min-max optimization, the training error is no longer solely quantified by
the loss function but rather by the duality gap of the loss function.

### Organization

The paper is organized as follows. In the next section, we describe our
method. First, we reformulate the PDE problem as a constraint optimization in
Section 2.1, and discuss the network structure in Section 2.2. Our algorithm
is given in Section 2.3. Section 3 is devoted to the theoretical analysis of
our algorithm. We prove the equivalence between the optimization problem and
the PDE problem in Section 3.1, Theorem 1. In Section 4, we derive an error
estimate for the proposed algorithm, with a decomposition of the error in
Section 4.1, Theorem 2, with the bounds of each part provided in Section 4.2.
We demonstrate in several numerical experiments the accuracy of our algorithm
in Section 5. Finally, some concluding remarks are given in Section 6.

Notations. We take $L^{2}$ to be the square integrable space with norm defined
by $\|f\|_{L^{2}(\Omega)}=(\int_{\Omega}|f(x)|^{2}dx)^{1/2}$ with
$\Omega\subset\mathbb{R}^{d}$ being an open domain of $\mathbb{R}^{d}$. The
boundary set of the domain $\Omega$ is denoted by $\partial\Omega$. We use
$H^{s}(\Omega)$ to denote the Sobolve space with the up-to $s$-th derivatives
bounded in $L^{2}$. We use $\mathcal{U}_{\theta}$ to denote a U-network, and
$\mathcal{V}_{\tau}$ as V-network, with $\theta,\tau$ as respective network
parameters. For a positive integer $m$, we use $[m]$ to represent
$\\{1,\cdots,m\\}$.

##  2\. Proposed Method

We initiate by presenting an abstract form of the PDE problem, using the
second order linear elliptic equations as a canonical example.

###  2.1. Reformulation as constrained optimization

Consider solving a boundary value problem of the form

$\displaystyle\mathcal{A}u$ $\displaystyle=f,\quad\text{in }\Omega,$ (1)
$\displaystyle\mathcal{B}u$ $\displaystyle=g,\quad\text{on }\partial\Omega,$

where $\mathcal{A}$ is a differential operator with respect to $x$,
$\mathcal{B}$ is a boundary operator which represents Dirichlet, Neumann,
Robin or mixed boundary conditions. Here $\Omega$ is a bounded domain
$\Omega\subset\mathbb{R}^{d}$, for dimension $d\geq 1$, and $f=f(x),g=g(x)$
are given functions defined in $\Omega$ and on the boundary $\partial\Omega$,
respectively. Here $u=u(x):\Omega\mapsto\mathbb{R}$ is the unknown function we
want to seek for. To proceed, we assume that the above boundary value problem
is well-posed. Here as an example, we take

$\displaystyle\mathcal{A}u:=-\sum_{i=1}^{d}\left(\sum_{j=1}^{d}a_{ij}(x)u_{x_{j}}\right)_{x_{i}}+\sum_{i=1}^{d}b_{i}(x)u_{x_{i}}+c(x)u,$
(2)

and

$\displaystyle\mathcal{B}u:=\alpha u+\beta\frac{\partial u}{\partial n},$ (3)

with $\alpha>0$, $\beta\geq 0$, and $\frac{\partial u}{\partial n}$ denoting
the directional derivative of $u$ along the exterior normal direction $n=n(x)$
at the boundary point $x\in\partial\Omega$. Here $a_{ij},b_{i},c_{i}$ are
bounded functions defined in $\Omega$ and $a_{ij}$ are differentiable and
satisfy $\sum_{i,j=1}^{d}a_{ij}(x)\xi_{i}\xi_{j}\geq\lambda\|\xi\|^{2}$ for
some constant $\lambda>0$ for any $\xi\in\mathbb{R}^{d}$ and $x\in\Omega$,
where $a_{ij}=a_{ji}$ for all $i,j\in[d]$ (ellipticity). For the above
elliptic equations, we assume $f\in L^{2}(\Omega),g\in
H^{3/2}(\partial\Omega)$, so the corresponding boundary value problem is well-
posed with solution in $H^{2}(\Omega)$ (see for example [16]).

To solve (1), we reformulate it as a constraint optimization problem

$\inf_{u\in
U}\left\\{\frac{1}{2}\|\mathcal{B}u-g\|_{L^{2}(\partial\Omega)}^{2}\right\\},\quad\text{where
}U=\\{u\in H^{2}(\Omega),\quad\mathcal{A}u=f\;\;\text{ in }\Omega\\}.$ (4)

In this formulation, the objective function is the least square of the
boundary data mismatch, and the PDE itself serves as a constraint. The
equivalence between this optimization problem and equation (1) is obvious, and
existence of problem (4) will be established in the next section.

To remove the PDE constraint, we introduce a Lagrangian multiplier $v\in
L^{2}(\Omega)$, so that the above problem can be converted to an inf-sup
optimization problem

$\inf_{u\in H^{2}(\Omega)}\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v):=\frac{1}{2}\|\mathcal{B}u-g\|_{L^{2}(\partial\Omega)}^{2}+(\mathcal{A}u-f,v),$
(5)

where $(\cdot,\cdot)$ denotes the $L^{2}$ inner product in $\Omega$. Here the
loss function is made up of two terms: the first term takes into consideration
of the boundary condition and the second term accounts for the underlying PDE.
The saddle point of the above Lagrangian can be shown to be the solution to
problem (1), see Theorem 1 in the next section.

###  2.2. Neural networks

Here we use two feedforward neural networks (FNNs) to represent the trial
solution and the Lagrangian multiplier. Given input $x$, an FNN transforms it
to an output, which can be seen as a recursively defined function $u_{\theta}$
on $\mathbb{R}^{d}$ into $\mathbb{R}$:

$u_{\theta}(x)=h_{m}\circ\sigma\cdots h_{2}\circ\sigma\circ h_{1}(x).$

Here $h_{k}$ is an affine map from $\mathbb{R}^{N_{k-1}}$ to
$\mathbb{R}^{N_{k}}$ (here we take $N_{0}=d,N_{m}=1$):

$h_{k}(z)=W_{k}z+b_{k},\quad\theta=\\{W_{k},b_{k}\\}_{k=1}^{m},$

where $W_{k}$ are $N_{k}\times N_{k-1}$ matrices (network weights) and $b_{k}$
are $N_{k}$–vectors (network biases). $\sigma$ is a nonlinear activation
function. As the present method can be combined with other neural network
constructions, FNN serves a base model class to demonstrate our method.
Moreover, by the universal approximation theorem of FNNs [11, 41, 1], with the
use of a proper activation, FFNs of sufficient width and depth can model any
given $H^{2}$ and $L^{2}$ function which covers our solution space and the
function space of the Lagrangian multiplier. Therefore, we use FFNs to model
both trial solutions and the Lagrangian multiplier. Required by assumptions of
the universal approximation theorems ([11, 41, 1]), the activation function
$\sigma$ can be taken as any differentiable non-constant bounded function.
Examples include sigmoid function and $\tanh$.

###  2.3. The Inf-Sup network

The inf-sup reformulation of the PDE problem led us to design the inf-sup
network, which consists of a neural network $\mathcal{U}_{\theta}$
representing the trial solution and another adversarial network
$\mathcal{V}_{\tau}$ representing the Lagrangian multiplier. The subscripts
indicate the parameters of each neural network adopted. We can thus
approximate problem (5) by

$\displaystyle\min_{\theta}\max_{\tau}L(\theta,\tau):=\mathcal{L}(u_{\theta},v_{\tau}),$
(6)

with $u_{\theta}\in\mathcal{U}_{\theta}$ and $v_{\tau}\in\mathcal{V}_{\tau}$.
The differential operator $\mathcal{A}$ and $\mathcal{B}$ in $\mathcal{L}$ are
computed via auto-differentiation of neural networks.

With the above setup, the approximation solution to the inf-sup optimization
problem (6) can be obtained by finding a saddle point for the following min-
max problem:

$\displaystyle\min_{\theta}\max_{\tau}{\hat{L}}(\theta,\tau),$ (7)

where $\hat{L}$ denotes the evaluation of ${\mathcal{L}}(u_{\theta},v_{\tau})$
over samples using the Monte-Carlo integration [48]:

$\displaystyle\hat{L}(\theta,\tau):=\hat{\mathcal{L}}(u_{\theta},v_{\tau})=$
$\displaystyle\frac{|\partial\Omega|}{2N_{b}}\sum_{i=1}^{N_{b}}(\mathcal{B}u_{\theta}(x_{b}^{i})-g(x_{b}^{i}))^{2}+\frac{|\Omega|}{N}\sum_{i=1}^{N}(\mathcal{A}u_{\theta}(x^{i})-f(x^{i}))\cdot
v_{\tau}(x^{i}).$ (8)

Here $N_{b}$ is the number of sampling points $\\{x_{b}^{i}\\}_{i=1}^{N_{b}}$
on the boundary $\partial\Omega$ and $N$ is the number of sampling points
$\\{x^{i}\\}_{i=1}^{N}$ inside the domain $\Omega$. Note that other numerical
integration methods could also be used.

Problem (7) can be solved using techniques of the min-max optimization such as
Stochastic Gradient Descent–Ascent [46] . The algorithm is described in
Algorithm 1.

Algorithm 1 Training of InfSupNet

for number of training iterations do

$\bullet$ Sample $N$ points $\\{x^{1},\dots,x^{N}\\}$ randomly in the domain
$\Omega\subset\mathbb{R}^{d}$.

$\bullet$ Sample $N_{b}$ points $\\{{x}_{b}^{(1)},\dots,{x}_{b}^{(N_{b})}\\}$
randomly from the boundary set $\partial\Omega\subset\mathbb{R}^{d-1}$.

for $k$ steps do

$\bullet$ Update the $\mathcal{V}_{\tau}$-network by ascending its
(stochastic) gradient:
$\tau\leftarrow\tau+\eta\nabla_{\tau}\hat{L}(\theta,\tau)$

end for

for $\ell$ steps do

$\bullet$ Update the $\mathcal{U}_{\theta}$-network by descending its
(stochastic) gradient:
$\theta\leftarrow\theta-\eta\nabla_{\theta}\hat{L}(\theta,\tau)$

end for

end for

return $\theta,\tau$ for the parameters of the neural network
$\mathcal{U}_{\theta},\mathcal{V}_{\tau}$.

##  3\. Theoretical results

In this section, we provide theoretical analysis for our proposed method. For
simplicity of presentation, we present analysis and the result only for the
Poisson equation with $\mathcal{A}=-\Delta$ being the negative Laplacian and
$\mathcal{B}=I$ being the identical map. The results and analysis extend well
to problem (2)-(3) in a straightforward fashion. The Poisson problem reads as

$\displaystyle-\Delta u=f,$ $\displaystyle\quad\text{in }\Omega,$ (9)
$\displaystyle u=g,$ $\displaystyle\quad\text{on }\partial\Omega.$

###  3.1.  Equivalence between the inf-sup problem (5) and the PDE formulation

In this subsection, we prove the following theorem.

######  Theorem 1.

Assume $\Omega$ is a bounded domain with $\partial\Omega$ being smooth, and
$f\in L^{2}(\Omega),{g\in H^{2/3}(\partial\Omega)}$. Then $u$ is the unique
solution to (9) if and only if it solves the optimization problem (5) with
$\mathcal{B}u=u$ and $\mathcal{A}u=-\Delta u$.

In order to prove the above theorem and show the equivalence between the inf-
sup optimization problem (5) and problem (9), we utilize the nice property of
the saddle points.

######  Definition 1.

A pair of functions $(u,v)\in H^{2}(\Omega)\times L^{2}(\Omega)$ is called a
_saddle point_ of the Lagrangian functional $\mathcal{L}(u,v)$ if the
following inequality holds

$\displaystyle\mathcal{L}(u,q)\leq\mathcal{L}(u,v)\leq\mathcal{L}(w,v)\quad\text{for
any }w\in H^{2}(\Omega),q\in L^{2}(\Omega).$ (10)

By Lemma 4 in the appendix, a pair of function $(u^{*},v^{*})$ is a saddle
point of the Lagrangian $\mathcal{L}$ if and only if the inf-sup problem (5)
satisfies the _strong max-min property_

$\displaystyle\inf_{u\in H^{2}(\Omega)}\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v)=\sup_{v\in L^{2}(\Omega)}\inf_{u\in
H^{2}(\Omega)}\mathcal{L}(u,v),$ (11)

and $u^{*}$ and $v^{*}$ are the corresponding primal and dual optimal
solutions. Therefore, to show the equivalence between problem (5) and the PDE
problem (9), it suffices to show the max-min property, and that the solution
to (9) and the saddle point of $\mathcal{L}$ are equivalent.

#### The strong max-min property

Suppose $u$ is the optimal solution to (5), by the linearity of the
constraint, we have

$\displaystyle\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v)=\left\\{\begin{array}[]{cc}+\infty,&\text{if
}-\Delta u-f\neq 0,\\\
\frac{1}{2}\|u-g\|_{L^{2}(\partial\Omega)}^{2},&\text{otherwise}.\end{array}\right.$
(14)

Since for $f,g$ satisfying the assumption of Theorem 1, (9) has a unique
solution, $\inf_{u\in H^{2}(\Omega)}\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v)=0$. On the other hand,

$\displaystyle\inf_{u\in
L^{2}(\Omega)}\mathcal{L}(u,v)=\left\\{\begin{array}[]{cc}-\infty,&\text{if
}v\neq 0,\\\ 0,&\text{otherwise},\end{array}\right.$ (17)

since for $v\neq 0$, one can take $-\Delta u-f=-\lambda v$ and $u=g$ on
$\partial\Omega$ so that $\mathcal{L}(u,v)\to-\lambda\|v\|^{2}$, and let
$\lambda\to\infty$ to get $\mathcal{L}(u,v)\to-\infty$. Therefore,

$\sup_{u\in H^{2}(\Omega)}\inf_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v)=0=\inf_{u\in H^{2}(\Omega)}\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v)$

and the strong max-min property holds.

#### Equivalence between the saddle point and PDE

First we derive the equation of the saddle point. Suppose $(u,v)$ is a saddle
point of $\mathcal{L}$, then by definition (10), $\mathcal{L}(u,v+\varepsilon
q)\leq\mathcal{L}(u,v)\leq\mathcal{L}(u+\varepsilon w,v)$ for any $q\in
L^{2}(\Omega)$, $w\in H^{2}(\Omega)$ and $\varepsilon\in\mathbb{R}$. By the
arbitrariness of $\varepsilon$, one can derive the following equations:

$\displaystyle\int_{\partial\Omega}uwds+\int_{\Omega}(-\Delta
w)vdx=\int_{\partial\Omega}wgds,\quad\text{for any }w\in H^{2}(\Omega),$ (18)
$\displaystyle\int_{\Omega}(-\Delta u-f)qdx=0,\quad\text{for any }q\in
L^{2}(\Omega).$

Note that the converse is also true. If $(u,v)$ satisfies the above equations,
we can compute $\mathcal{L}(u,q)-\mathcal{L}(u,v)=\int_{\Omega}(-\Delta
u-f)(q-v)dx$, which by the second equation of (18), equals zero, and a careful
calculation shows

$\mathcal{L}(w,v)-\mathcal{L}(u,v)=\frac{1}{2}\|u-w\|_{L^{2}(\partial\Omega)}^{2}\geq
0\quad\text{for any}\quad w\in H^{2}(\Omega).$

Therefore, (10) holds.

Next we show that the solution to (9) and the solution to (18) are the same.

######  Lemma 1.

Let the assumptions of Theorem 1 hold. Then $u$ is a solution to (9) if and
only if it is a solution to (18).

###### Proof.

First, by the regularity theory of elliptic equations, there exists a unique
solution $u\in H^{2}(\Omega)$ to (9) (see for example [16]). It is simple to
verify that $(u,v=0)$ is also a solution to (18). Hence it suffices to show
the uniqueness of the solution to (18). By the linearity of the problem (18),
we need to show that for $f=0,g=0$, (18) only has a zero solution. First, from
the second equation of (18) with $f=0$, we get

$\displaystyle-\Delta u=0\quad\text{in }\Omega,$ (19)

almost everywhere. Let $h$ be an arbitrary function in $L^{2}(\partial\Omega)$
and $w\in H^{2}(\Omega)$ be the solution to equation

$\displaystyle-\Delta w=0\quad\text{in }\Omega,\quad w=h\quad\text{on
}\partial\Omega,$

then with (18), we have

$\displaystyle\int_{\partial\Omega}uhds=0\;\text{for any }h\in
L^{2}(\partial\Omega),$

which implies $u=0$ almost everywhere on $\partial\Omega$. This, together with
(19), implies $u(x)=0$ for $x\in\Omega$.

To show $v=0$, let $r\in L^{2}(\Omega)$ be an arbitrary function and we take
the test function $w$ to be the solution to

$\displaystyle-\Delta w=r\quad\text{in }\Omega,\quad w=0\quad\text{on
}\partial\Omega,$

then (18) becomes

$\displaystyle\int_{\Omega}rvdx=0\quad\text{for any }r\in L^{2}(\Omega),$ (20)

which implies $v=0$ almost everywhere. This finishes the proof. ∎

#### Proof of Theorem 1.

If $u$ is a solution to (9), then by Lemma 1, $(u,v=0)$ solves (18) and thus
$(u,v=0)$ is a saddle point of $\mathcal{L}$. Therefore, it solves (11) and
thus problem (5). Conversely, if $(u,v)$ solves (5), then by the strong max-
min property, the primal and dual optimal solutions $u$ and $v$ are the saddle
point of $\mathcal{L}$. Hence $(u,v)$ solves (18). And by Lemma 1, $u$ is the
solution to (9).

###  3.2. The existence of optimization problem (4) and its equivalence to the
PDE formulation

Next, we establish the existence of the optimization problem (4) and
demonstrate that its solution satisfies (9). First, we assert the Brezzi
condition:

$\displaystyle\inf_{v\in L^{2}(\Omega)}\sup_{u\in
H^{2}(\Omega)}\frac{\int_{\Omega}-\Delta
uvdx}{\|u\|_{H^{2}(\Omega)}\|v\|_{L^{2}(\Omega)}}\geq\gamma,$ (21)

where $\gamma>0$ is a positive constant. Specifically, for each $v$, we can
set $-\Delta u=v$ in $\Omega$ and $u=0$ on $\partial\Omega$ and use the
stability estimate $\|u\|_{H^{2}(\Omega)}\leq C\|v\|_{L^{2}(\Omega)}$ for the
Poisson problem. This leads to inequality (21) with $\gamma=1/C$.

######  Lemma 2.

Let the assumptions of Theorem 1 hold. The optimization problem (4) has a
unique solution $u\in H^{2}(\Omega)$ that solves (9).

###### Proof.

Let $V_{f}:=\\{u\in H^{2}(\Omega):-\Delta u=f\\}$ be a subspace of
$H^{2}(\Omega)$. By the theory of elliptic equations, $V_{f}$ is not empty for
$f\in L^{2}(\Omega)$. Taking an element $u_{f}\in V_{f}$ such that $u_{f}\neq
0$. Define $w:=u-u_{f}$ and the problem (4) becomes

$\displaystyle\inf_{w\in
V_{0}}\frac{1}{2}\int_{\partial\Omega}(u_{f}+w-g)^{2}ds.$

Introducing the bilinear form

$a(w_{1},w_{2})=\int_{\partial\Omega}w_{1}\cdot w_{2}ds\quad\text{and}\quad
b(w,g-u_{f})=\int_{\partial\Omega}w(g-u_{f})ds,$

the above optimization problem can be rewritten in the form

$\displaystyle\inf_{w\in V_{0}}\frac{1}{2}a(w,w)-b(w,g-u_{f}).$ (22)

Next we show $a(w_{1},w_{2})$ is a continuous coercive bilinear form on
$V_{0}$. First,

$a(w_{1},w_{2})\leq\|w_{1}\|_{L^{2}(\partial\Omega)}\|w_{2}\|_{L^{2}(\partial\Omega)}\leq
C\|w_{1}\|_{H^{2}(\Omega)}\|w_{2}\|_{H^{2}(\Omega)}$

by the trace theorem. Moreover, by the regularity of elliptic equation

$\displaystyle-\Delta w=0,\text{ in }\Omega,\quad w=w|_{\partial\Omega},\text{
on }\partial\Omega,$

we have $\|w\|_{H^{2}(\Omega)}\leq C\|w\|_{L^{2}(\partial\Omega)}.$ Hence
$a(w_{1},w_{2})$ is coercive on $V_{0}$, satisfying

$\displaystyle a(w,w)\geq\frac{1}{C}\|w\|_{H^{2}(\Omega)}^{2}.$

Due to the above property of $a(\cdot,\cdot)$ and $V_{0}$ being a nonempty
closed and convex set, we apply the Lax-Milgram theorem to conclude that there
exists a unique solution $w^{*}\in V_{0}$ to the problem (22). Hence there
exists a unique solution $u^{*}=w^{*}+u_{f}\in H^{2}(\Omega)$ to the
optimization problem (4). Moreover, the unique solution $u$ satisfies

$\displaystyle a(w^{*},s)=b(s,g-u_{f}),\quad\text{ for any }s\in V_{0}.$

Next we extend the space of $s$ from $V_{0}$ to $H^{2}$. Let $\ell$ be a
linear functional defined on $H^{2}(\Omega)$ by
$\ell(s):=a(w^{*},s)-b(s,g-u_{f})$. By the Brezzi condition (21), there exists
a unique solution $v^{*}$ such that (see [5, Theorem 4.2.1])

$\displaystyle\ell(s)=-\int_{\Omega}(-\Delta s)v^{*}dx\quad\text{ for any
}s\in H^{2}(\Omega).$

Using $u^{*}=w^{*}+u_{f}$, we can get

$\displaystyle a(u^{*}-u_{f},s)-b(s,g-u_{f})=-\int_{\Omega}(-\Delta
s)v^{*}dx,$

which is

$\displaystyle\int_{\partial\Omega}u^{*}sds+\int_{\Omega}(-\Delta
s)v^{*}dx=\int_{\partial\Omega}gsds,$

and hence $(u^{*},v^{*})$ satisfies (18). Moreover, since $\Delta w^{*}=0$ and
$-\Delta u_{f}=f$, $u^{*}$ also satisfies the second equation of (18). By
Lemma 1, $u^{*}$ is a solution to equation (9).

∎

##  4\. Error analysis

In this section, we analyze the generalization error of our method. Let
$(u_{n},v_{n})$ represent the solution obtained from InfSupNet after $n$ steps
training and let $u$ be the precise solution to equation (9). We will derive
an error estimate for the discrepancy between $(u_{n},v_{n})$ and $(u,v)$,
where $v=0$.

Before conducting the error analysis, let’s revisit a well-established result
in the literature (refer to [50]). Such an inequality has previously played a
crucial role in the error analysis of PINNs for elliptics PDEs [51, 60]).

######  Lemma 3.

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain with Lipschitz boundary
$\partial\Omega$. For any $w\in H^{2}(\Omega)$, there exists a constant $C>0$
depending only on $\Omega$ such that

$\displaystyle\|w\|_{H^{1/2}(\Omega)}\leq C(\|\Delta
w\|_{L^{2}(\Omega)}+\|w\|_{L^{2}(\partial\Omega)}).$

###### Proof.

Using the elliptic regularity inequality given in [50, Theorem 2.1]), we have

$\displaystyle\|w\|_{H^{s}(\Omega)}\leq C(\|\Delta
w\|_{H^{s-2}(\Omega)}+\|w\|_{H^{s-1/2}(\partial\Omega)}+\|w\|_{H^{s-2}(\Omega)}),$

for any $s$ real. Taking $s=1/2$ and using the embedding
$\|v\|_{H^{-\alpha}(\Omega)}\leq C\|v\|_{L^{2}(\Omega)}$ for $\alpha\geq 0$,
the above inequality implies

$\displaystyle\|w\|_{H^{1/2}(\Omega)}\leq C(\|\Delta
w\|_{H^{-3/2}(\Omega)}+\|w\|_{L^{2}(\partial\Omega)}+\|w\|_{H^{-3/2}(\Omega})\leq
C(\|\Delta
w\|_{L^{2}(\Omega)}+\|w\|_{L^{2}(\partial\Omega)}+\|w\|_{L^{2}(\Omega)}).$

Furthermore, applying Lemma 2.1 from [34] implies

$\displaystyle\|w\|_{L^{2}(\Omega)}\leq C(\|\Delta
w\|_{H^{-2}(\Omega)}+\|w\|_{H^{-1/2}(\partial\Omega)})\leq C(\|\Delta
w\|_{L^{2}(\Omega)}+\|w\|_{L^{2}(\partial\Omega)}).$

Combining these inequalities yield (23) as asserted. ∎

Applying Lemma 3 to $u_{n}-u$ gives

$\displaystyle\|u_{n}-u\|_{H^{1/2}(\Omega)}$ $\displaystyle\leq
C(\|\Delta(u_{n}-u)\|_{L^{2}(\Omega)}+\|u_{n}-u\|_{L^{2}(\partial\Omega)}).$
(23)

Hence to get the error estimate, it suffices to derive a bound on the terms on
the right hand side of (23).

###  4.1. Error decomposition

We will break down the error into various components, encompassing the Monte
Carlo sampling error, the training loss gap, and the neural network
approximation error. The principal outcome is articulated in the following
theorem.

######  Theorem 2.

Let the assumptions of Theorem 1 hold. Then the following estimate holds:

$\displaystyle\|u_{n}-u\|_{L^{2}(\partial\Omega)}^{2}+\|\Delta(u_{n}-u)\|_{L^{2}(\Omega)}^{2}+\|v_{n}\|_{L^{2}(\Omega)}^{2}\leq
I_{NN}+I_{MC}+I_{GP},$ (24)

where $I_{NN}$ is the neural network approximation error, given by

$\displaystyle
I_{NN}=|\mathcal{L}(u_{n},q_{\tilde{v}})-\mathcal{L}(u_{n},\tilde{v})|+|\mathcal{L}(w_{\tilde{u}},v_{n})-\mathcal{L}(\tilde{u},v_{n})|,$
(25)

where $\tilde{v}=-\Delta u_{n}-f$, and $\tilde{u}$ is the solution to

$\displaystyle-\Delta\tilde{u}$ $\displaystyle=f-v_{n},\quad\text{in }\Omega,$
(26) $\displaystyle\tilde{u}$ $\displaystyle=g,\quad\text{on }\partial\Omega.$

Here $q_{\tilde{v}},w_{\tilde{u}}$ are the best approximation of
$\tilde{v},\tilde{u}$ in $\mathcal{V}_{\tau},\mathcal{U}_{\theta}$ (see (30)
and (33) for the definitions), and $I_{MC}$ is the Mento-Carlo sampling error
defined by

$\displaystyle
I_{MC}=|\mathcal{L}(u_{n},q_{\tilde{v}})-\hat{\mathcal{L}}(u_{n},q_{\tilde{v}})|+|\mathcal{L}(w_{\tilde{u}},v_{n})-\hat{\mathcal{L}}(w_{\tilde{u}},v_{n})|,$
(27)

and $I_{GP}$ is the duality gap of the training loss

$\displaystyle
I_{GP}=\sup_{q\in\mathcal{V}_{\tau}}\hat{\mathcal{L}}(u_{n},q)-\inf_{w\in\mathcal{U}_{\theta}}\hat{\mathcal{L}}(w,v_{n}).$
(28)

Moreover, by (23), inequality (24) implies

$\displaystyle\|u_{n}-u\|_{H^{1/2}(\Omega)}^{2}+\|v_{n}\|_{L^{2}(\Omega)}^{2}\leq
C(I_{NN}+I_{MC}+I_{GP}).$ (29)

###### Proof.

First we estimate $u_{n}-u$. For $u_{n}\in\mathcal{U}_{\theta}\subset
H^{2}(\Omega)$, we define $\tilde{v}=-\Delta u_{n}-f$ so that

$\|u_{n}-u\|_{L^{2}(\partial\Omega)}^{2}+\|\Delta(u_{n}-u)\|_{L^{2}(\Omega)}^{2}\leq
2\mathcal{L}(u_{n},\tilde{v}).$

Let $q_{\tilde{v}}$ be the best approximation of
$\tilde{v}\in\mathcal{V}_{\tau}$ such that

$\displaystyle
q_{\tilde{v}}=\arg\inf_{q\in\mathcal{V}_{\tau}}\|q-\tilde{v}\|_{L^{2}(\Omega)},$
(30)

then

$\displaystyle\mathcal{L}(u_{n},\tilde{v})$
$\displaystyle=\mathcal{L}(u_{n},\tilde{v})-\mathcal{L}(u_{n},q_{\tilde{v}})+\mathcal{L}(u_{n},q_{\tilde{v}})\leq|\mathcal{L}(u_{n},\tilde{v})-\mathcal{L}(u_{n},q_{\tilde{v}})|+\mathcal{L}(u_{n},q_{\tilde{v}}).$
(31)

The first term on the right side gives the first part of the approximation
error (25). Next we consider the second term and rewrite
$\mathcal{L}(u_{n},q_{\tilde{v}})$ using $\hat{\mathcal{L}}$, which is defined
in (8), as

$\displaystyle\mathcal{L}(u_{n},q_{\tilde{v}})$
$\displaystyle=\mathcal{L}(u_{n},q_{\tilde{v}})-\hat{\mathcal{L}}(u_{n},q_{\tilde{v}})+\hat{\mathcal{L}}(u_{n},q_{\tilde{v}})\leq|\mathcal{L}(u_{n},q_{\tilde{v}})-\hat{\mathcal{L}}(u_{n},q_{\tilde{v}})|+\sup_{q\in\mathcal{V}_{\tau}}\hat{\mathcal{L}}(u_{n},q).$

The first term on the right of the last inequality gives the first part of the
Mento-Carlo sampling error (27), and the last term gives the first part of the
duality gap (28). Combing the above inequality with (31), we obtain

$\displaystyle\frac{1}{2}\|$ $\displaystyle
u_{n}-u\|_{L^{2}(\partial\Omega)}^{2}+\|\Delta(u_{n}-u)\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq|\mathcal{L}(u_{n},\tilde{v})-\mathcal{L}(u_{n},q_{\tilde{v}})|+|\mathcal{L}(u_{n},q_{\tilde{v}})-\hat{\mathcal{L}}(u_{n},q_{\tilde{v}})|+\sup_{q\in\mathcal{V}_{\tau}}\hat{\mathcal{L}}(u_{n},q).$
(32)

Next we estimate the difference between $v_{n}$ and $v=0$. For any
$v_{n}\in\mathcal{V}_{\tau}\subset L^{2}(\Omega)$, there exists a solution
$\tilde{u}\in H^{2}(\Omega)$ to (26). Taking $u=\tilde{u}$ in
$\mathcal{L}(u,v_{n})$, we get

$\displaystyle\mathcal{L}(\tilde{u},v_{n})=-\|v_{n}\|_{L^{2}(\Omega)}^{2}.$

Let $w_{\tilde{u}}$ be the best approximation of $\tilde{u}$ in
$\mathcal{U}_{\theta}$, i.e.

$\displaystyle
w_{\tilde{u}}=\arg\inf_{w\in\mathcal{U}_{\theta}}\|w-\tilde{u}\|_{H^{2}(\Omega)},$
(33)

then

$\displaystyle\mathcal{L}(\tilde{u},v_{n})$
$\displaystyle=\mathcal{L}(\tilde{u},v_{n})-\mathcal{L}(w_{\tilde{u}},v_{n})+\mathcal{L}(w_{\tilde{u}},v_{n})\geq-|\mathcal{L}(\tilde{u},v_{n})-\mathcal{L}(w_{\tilde{u}},v_{n})|+\mathcal{L}(w_{\tilde{u}},v_{n}).$
(34)

The first term on the right reflects the approximation error of the neural
network. For the second term, we have

$\displaystyle\mathcal{L}(w_{\tilde{u}},v_{n})$
$\displaystyle=\mathcal{L}(w_{\tilde{u}},v_{n})-\hat{\mathcal{L}}(w_{\tilde{u}},v_{n})+\hat{\mathcal{L}}(w_{\tilde{u}},v_{n})\geq-|\mathcal{L}(w_{\tilde{u}},v_{n})-\hat{\mathcal{L}}(w_{\tilde{u}},v_{n})|+\inf_{w\in\mathcal{U}_{\theta}}\mathcal{L}(w,v_{n}).$

Combing the above inequality with (34), we get the estimate

$\displaystyle\|v_{n}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq|\mathcal{L}(\tilde{u},v_{n})-\mathcal{L}(w_{\tilde{u}},v_{n})|+|\mathcal{L}(w_{\tilde{u}},v_{n})-\hat{\mathcal{L}}(w_{\tilde{u}},v_{n})|-\inf_{w\in\mathcal{U}_{\theta}}\hat{\mathcal{L}}(w,v_{n}).$
(35)

Combing (4.1) and (35) leads to (24) and finishes the proof. ∎

###  4.2. Error bounds

To specify, we choose the hyperbolic tangent (tanh) as the activation function
for both the $\mathcal{V}_{\tau}$ and $\mathcal{U}_{\theta}$ networks.
Initially, our focus is on the neural network approximation error $I_{NN}$ as
defined in (25). This error is connected to the approximation inaccuracies of
neural networks within the function spaces $L^{2}(\Omega)$ and
$H^{2}(\Omega)$.

For the first term in (25), we have

$\displaystyle|\mathcal{L}(u_{n},\tilde{v})-\mathcal{L}(u_{n},q_{\tilde{v}})|$
$\displaystyle=\int_{\Omega}(-\Delta u_{n}-f)\cdot(\tilde{v}-q_{\tilde{v}})dx$
$\displaystyle\leq\int_{\Omega}\tilde{v}\cdot(\tilde{v}-q_{\tilde{v}})dx$
$\displaystyle\leq\|\tilde{v}\|_{L^{2}(\Omega)}\|\tilde{v}-q_{\tilde{v}}\|_{L^{2}(\Omega)}.$

From Proposition 1 (see in Appendix), there exists $\mathcal{V}_{\tau}$ (two
hidden layers with width given by the proposition) such that

$\displaystyle\|\tilde{v}-q_{\tilde{v}}\|_{L^{2}(\Omega)}$
$\displaystyle=\inf_{q\in\mathcal{V}_{\tau}}\|q-\tilde{v}\|_{L^{2}(\Omega)}\leq|\Omega|^{\frac{1}{2}}\inf_{q\in\mathcal{V}_{\tau}}\|q-\tilde{v}\|_{L^{\infty}(\Omega)}\leq|\Omega|^{\frac{1}{2}}C_{1}M^{-s+2}\|\tilde{v}\|_{W^{s-2,\infty}(\Omega)},$
(36)

for any $s>2$, where $C_{1}>0$ is the constant specified in Proposition 1 and
depends only on $d,s,|\Omega|$, and $M$ is a constant depending on the width
of the network. Here we use the inequality
$\|\tilde{v}\|_{L^{2}(\Omega)}\leq|\Omega|^{\frac{1}{2}}\|\tilde{v}\|_{L^{\infty}(\Omega)}$.
Hence

$\displaystyle|\mathcal{L}(u_{n},\tilde{v})-\mathcal{L}(u_{n},q_{\tilde{v}})|\leq
C_{1}|\Omega|M^{-s+2}\|\tilde{v}\|_{W^{s-2,\infty}(\Omega)}^{2}.$ (37)

If ReLU were taken as the activation function for $\mathcal{V}_{\tau}$, one
could use the results from [17, 52] to derive a similar estimate.

The second term in (25) when using equation (26) is bounded by

$\displaystyle|\mathcal{L}(w_{\tilde{u}}$
$\displaystyle,v_{n})-\mathcal{L}(\tilde{u},v_{n})|$
$\displaystyle=\left|\frac{1}{2}\|w_{\tilde{u}}-g\|_{L^{2}(\partial\Omega)}^{2}+\int_{\Omega}(-\Delta
w_{\tilde{u}}+\Delta\tilde{u})v_{n}dx\right|$
$\displaystyle\leq\frac{1}{2}\|w_{\tilde{u}}-\tilde{u}\|_{L^{2}(\partial\Omega)}^{2}+\|v_{n}\|_{L^{2}(\Omega)}\|\Delta(w_{\tilde{u}}-\tilde{u})\|_{L^{2}(\Omega)}$
$\displaystyle\leq
C(\|w_{\tilde{u}}-\tilde{u}\|_{H^{2}(\Omega)}^{2}+\|v_{n}\|_{L^{2}(\Omega)}\|w_{\tilde{u}}-\tilde{u}\|_{H^{2}(\Omega)}).$
(38)

Here $C>0$ is the constant induced in the use of the embedding
$\|h\|_{L^{2}(\Omega)},\|\Delta h\|_{L^{2}(\Omega)}\leq
C\|h\|_{H^{2}(\Omega)}$. By Proposition 1, there exists $\mathcal{U}_{\theta}$
(with two hidden layers whose size given by the proposition) such that

$\displaystyle\|w_{\tilde{u}}-\tilde{u}\|_{H^{2}(\Omega)}$
$\displaystyle=\inf_{w\in\mathcal{U}_{\theta}}\|w-\tilde{u}\|_{H^{2}(\Omega)}\leq|\Omega|^{\frac{1}{2}}\inf_{w\in\mathcal{U}_{\theta}}\|w-\tilde{u}\|_{W^{2,\infty}(\Omega)}$
$\displaystyle\leq C_{2}|\Omega|^{\frac{1}{2}}\max\\{R_{2}^{2},\log^{2}(\beta
M^{s+d+2})\\}M^{-s+2}\|\tilde{u}\|_{W^{s,\infty}(\Omega)},$ (39)

where $C_{2}$ is a positive constant depending on $d,s,|\Omega|$, $R_{2}>1$ is
a constant such that $\tanh^{\prime}x,\tanh^{\prime\prime}x$ are decreasing on
$[R_{2},\infty)$ and $\beta$ is a constant depending on $d,k,\tilde{u}$ and is
given by (52). Since both functions are decreasing on $[1,\infty)$, we can
take $R_{2}=1$. From (26), it follows that

$\displaystyle-\Delta(u-\tilde{u})=v_{n},\quad\text{in }\Omega,\quad
u-\tilde{u}=0,\quad\text{on }\partial\Omega,$

which by the regularity theory of this elliptic problem yields the following
estimate

$\displaystyle\|u-\tilde{u}\|_{W^{s,\infty}(\Omega)}\leq
C\|v_{n}\|_{W^{s-2,\infty}(\Omega)}.$

Hence

$\displaystyle\|\tilde{u}\|_{W^{s,\infty}(\Omega)}\leq\|\tilde{u}-u\|_{W^{s,\infty}(\Omega)}+\|u\|_{W^{s,\infty}(\Omega)}\leq
C\|v_{n}\|_{W^{s-2,\infty}(\Omega)}+\|u\|_{W^{s,\infty}(\Omega)}.$

Therefore, (4.2) implies

$\displaystyle|\mathcal{L}(w_{\tilde{u}},$ $\displaystyle
v_{n})-\mathcal{L}(\tilde{u},v_{n})|$ $\displaystyle\leq~{}$ $\displaystyle
CC_{2}^{2}|\Omega|\max\\{1,\log^{4}(\beta
M^{s+d+2})\\}M^{-2s+4}\|\tilde{u}\|_{W^{s,\infty}(\Omega)}^{2}$
$\displaystyle~{}+CC_{2}|\Omega|^{\frac{1}{2}}\max\\{1,\log^{2}(\beta
M^{s+d+2})\\}M^{-s+2}\|\tilde{u}\|_{W^{s,\infty}}\|v_{n}\|_{L^{2}(\Omega)}$
$\displaystyle\leq~{}$ $\displaystyle
CC_{2}^{2}|\Omega|\max\\{1,\log^{4}(\beta
M^{s+d+2})\\}M^{-2s+4}\|\tilde{u}\|_{W^{s,\infty}}^{2}$
$\displaystyle~{}+CC_{2}|\Omega|^{\frac{1}{2}}\max\\{1,\log^{2}(\beta
M^{s+d+2})\\}M^{-s+2}(\|\tilde{u}\|_{W^{s,\infty}}^{2}+\|v_{n}\|_{L^{2}(\Omega)}^{2})$
$\displaystyle\leq~{}$ $\displaystyle
C(C_{2}|\Omega|^{\frac{1}{2}}\max\\{1,\log^{2}(\beta
M^{s+d+2})\\}M^{-s+2}+C_{2}^{2}|\Omega|\max\\{1,\log^{4}(\beta
M^{s+d+2})\\}M^{-2s+4})$
$\displaystyle\cdot(\|v_{n}\|_{W^{s-2,\infty}(\Omega)}^{2}+\|u\|_{W^{s,\infty}(\Omega)}^{2}).$

Combing the above inequality with (37), we get

$\displaystyle I_{NN}\leq CM^{-s+2}+C\max\\{1,\log^{4}(\beta
M^{s+d+2})\\}M^{-2s+4},$ (40)

where $C>0$ is a constant depending on $d,s,|\Omega|$, and
$\|u\|_{W^{s,\infty}(\Omega)},\|u_{n}\|_{W^{s,\infty}(\Omega)}$ and
$\|v_{n}\|_{W^{s-2,\infty}(\Omega)}$. Considering the fact
$\lim_{M\to\infty}M^{-s+2}\log^{4}(\beta M^{s+d+2})=0$, we can set

$\displaystyle M_{0}=\min\\{M|\log^{4}(\beta M^{s+d+2})M^{-s+2}\leq 1\\},$
(41)

so that (40) implies for $M>M_{0}$,

$\displaystyle I_{NN}\leq 2CM^{-s+2}.$ (42)

Next we derive the estimate for the Monte-Carlo sampling error $I_{MC}$
defined by (27). By [6], we can get that for $N_{b}$ and $N$ large,

$\displaystyle\left|\frac{1}{2}\int_{\partial\Omega}(u_{n}-g)^{2}ds-\frac{|\partial\Omega|}{2N_{b}}\sum_{i=1}^{N_{b}}\left[(u_{n}(x_{b}^{i})-g(x_{b}^{i}))^{2}\right]\right|$
$\displaystyle\leq\frac{\sigma_{1}}{2\sqrt{N_{b}}}|\partial\Omega|\nu_{d-1},$
$\displaystyle\left|\int_{\Omega}((-\Delta
u_{n}-f)q_{\tilde{v}})dx-\frac{|\Omega|}{N}\sum_{i=1}^{N}[(-\Delta
u_{n}(x^{i})-f(x_{i}))q_{\tilde{v}}(x^{i})]\right|$
$\displaystyle\leq\frac{\sigma_{2}}{\sqrt{N}}|\Omega|\nu_{d},$ (43)

where $\nu_{d-1},\nu_{d}$ are the standard normal distributions in $d-1$ and
$d$ dimensions. Here $\sigma_{1}$ is the square root of the variance of
$(u_{n}-g)^{2}$ given by

$\displaystyle\sigma_{1}^{2}=\int_{\partial\Omega}\left|(u_{n}-g)^{2}-\frac{1}{|\partial\Omega|}\int_{\partial\Omega}(u_{n}-g)^{2}ds^{\prime}\right|^{2}ds,$

and $\sigma_{2}$ is the square root of the variance of $-\Delta u_{n}-f$ given
by

$\displaystyle\sigma_{2}^{2}=\int_{\Omega}\left|(-\Delta
u_{n}-f)q_{\tilde{v}}-\frac{1}{|\Omega|}\int_{\Omega}(-\Delta
u_{n}-f)q_{\tilde{v}}dx^{\prime}\right|^{2}dx.$

By $(a-b)^{2}\leq 2(a^{2}+b^{2})$ and Jensen’s inequality
$\left(\frac{1}{|\Omega|}\int_{\Omega}wdx\right)^{2}\leq\frac{1}{|\Omega|}\int_{\Omega}w^{2}dx$
for any $w\in L^{2}(\Omega)$, we get

$\displaystyle\sigma_{1}^{2}$ $\displaystyle\leq
2\|u_{n}-g\|_{L^{4}(\partial\Omega)}^{4}+2|\partial\Omega|\left(\frac{1}{|\partial\Omega|}\int_{\partial\Omega}(u_{n}-g)^{2}ds^{\prime}\right)^{2}\leq
4\|u_{n}-g\|_{L^{4}(\partial\Omega)}^{4}.$ (44)

Similarly,

$\sigma_{2}^{2}\leq 4\|-(\Delta
u_{n}-f)q_{\tilde{v}}\|_{L^{2}(\Omega)}^{2}\leq
4\|q_{\tilde{v}}\|_{L^{\infty}(\Omega)}^{2}\|-\Delta
u_{n}-f)\|_{L^{2}(\Omega)}^{2}.$

From (36), $q_{\tilde{v}}$ satisfies the inequality

$\displaystyle\|q_{\tilde{v}}-\tilde{v}\|_{L^{\infty}(\Omega)}=\inf_{v\in\mathcal{V}_{\tau}}\|v-\tilde{v}\|_{L^{\infty}(\Omega)}\leq
C_{1}M^{-s+2}\|\Delta u_{n}+f\|_{W^{s-2,\infty}(\Omega)},$

and hence

$\displaystyle\|q_{\tilde{v}}\|_{L^{\infty}(\Omega)}\leq\|q_{\tilde{v}}-\tilde{v}\|_{L^{\infty}(\Omega)}+\|\tilde{v}\|_{L^{\infty}(\Omega)}\leq(1+C_{1}M^{-s+2})\|\Delta
u_{n}+f\|_{W^{s-2,\infty}(\Omega)},$

and thus

$\displaystyle\sigma_{2}^{2}\leq 4(1+C_{1}M^{-s+2})^{2}\|\Delta
u_{n}+f\|_{W^{s-2,\infty}(\Omega)}^{4}.$ (45)

Therefore, combing (44) and (45) leads to

$\displaystyle|\mathcal{L}($ $\displaystyle
u_{n},q_{\tilde{v}})-\hat{\mathcal{L}}(u_{n},q_{\tilde{v}})|\leq\|u_{n}-g\|_{L^{4}(\partial\Omega)}^{2}|\partial\Omega|N_{b}^{-\frac{1}{2}}\nu_{d-1}+2(1+C_{1}M^{-s+2})\|\Delta
u_{n}+f\|_{W^{s-2,\infty}(\Omega)}^{2}|\Omega|N^{-\frac{1}{2}}\nu_{d}.$

Likewise, we can also bound
$\mathcal{L}(w_{\tilde{u}},v_{n})-\hat{L}(w_{\tilde{u}},v_{n})$,

$\displaystyle|\mathcal{L}($ $\displaystyle
w_{\tilde{u}},v_{n})-\hat{\mathcal{L}}(w_{\tilde{u}},v_{n})|$
$\displaystyle\leq
2\|w_{\tilde{u}}-g\|_{L^{4}(\partial\Omega)}^{2}|\partial\Omega|N_{b}^{-\frac{1}{2}}\nu_{d-1}+4\|v_{n}\|_{L^{2}(\Omega)}\|-\Delta
w_{\tilde{u}}-f\|_{L^{\infty}(\Omega)}|\Omega|N^{-\frac{1}{2}}\nu_{d}$
$\displaystyle\leq
C\|w_{\tilde{u}}-\tilde{u}\|_{L^{\infty}(\Omega)}^{2}|\partial\Omega|N_{b}^{-\frac{1}{2}}\nu_{d-1}+4\|v_{n}\|_{L^{2}(\Omega)}\|-\Delta
w_{\tilde{u}}+\Delta\tilde{u}-v_{n}\|_{L^{\infty}(\Omega)}|\Omega|N^{-\frac{1}{2}}\nu_{d}$
$\displaystyle\leq
C\|w_{\tilde{u}}-\tilde{u}\|_{L^{\infty}(\Omega)}^{2}|\partial\Omega|N_{b}^{-\frac{1}{2}}\nu_{d-1}+C(\|w_{\tilde{u}}-\tilde{u}\|_{W^{2,\infty}(\Omega)}+\|v_{n}\|_{L^{\infty}(\Omega)})\|v_{n}\|_{L^{2}(\Omega)}|\Omega|N^{-\frac{1}{2}}\nu_{d}.$

With (4.2), the above inequality implies

$\displaystyle|$
$\displaystyle\mathcal{L}(w_{\tilde{u}},v_{n})-\hat{L}(w_{\tilde{u}},v_{n})|$
$\displaystyle\leq CC_{2}^{2}\max\\{1,\log^{4}(\beta
M^{s+d+2})\\}M^{-2s+4}(\|v_{n}\|_{W^{s-2,\infty}(\Omega)}^{2}+\|u\|_{W^{s,\infty}}^{2})\cdot(|\partial\Omega|N_{b}^{-\frac{1}{2}}\nu_{d-1}+|\Omega|N^{-\frac{1}{2}}\nu_{d}).$

Combing this with (4.2) and (41), we obtain for $M>M_{0}$,

$\displaystyle I_{MC}\leq
C(1+M^{-s+2})(N_{b}^{-\frac{1}{2}}\nu_{d-1}+N^{-\frac{1}{2}}\nu_{d}),$ (46)

where $C>0$ is a constant depending on $d,s,|\Omega|$,
$\|u\|_{W^{s,\infty}(\Omega)},\|u_{n}\|_{W^{s,\infty}(\Omega)}$,
$\|v_{n}\|_{W^{s-2,\infty}(\Omega)}$.

Note that alternative numerical integration methods can be employed to attain
an enhanced error bound. A judicious selection includes utilizing Gaussian
quadrature in low dimensions and opting for the Quasi-Monte Carlo method or
the Monte-Carlo method when dealing with higher dimensions [48].

The analysis presented above enables us to arrive at a conclusion articulated
in the following theorem.

######  Theorem 3.

Let $\Omega\in\mathbb{R}^{d}$ be a domain with smooth boundary
$\partial\Omega$, and let $u\in H^{2}(\Omega)$ be the solution to the
Dirichlet problem (9), where $f\in W^{s,\infty}(\Omega),g\in
W^{s-2,\infty}(\partial\Omega)$ with $s>2,s\in\mathbb{N}$. For
$M\geq\max\\{5d^{2},M_{0}\\}$ with $M_{0}$ given by (41), there exists two
feedforward neural networks $\mathcal{U}_{\theta},\mathcal{V}_{\tau}$, each
with two hidden layers, activation function $\tanh$ and width given by
$\max\\{3\left[\frac{s}{2}\right]\binom{s+d-1}{s-1}+d(M-1),3\left[\frac{d+2}{2}\right]5^{d}M^{d}\\}$,
such that the result
$(u_{n},v_{n})\in\mathcal{U}_{\theta}\times\mathcal{V}_{\tau}$ after $n$-th
step training of InfSupNet satisfies

$\displaystyle\|$ $\displaystyle
u_{n}-u\|_{H^{1/2}(\Omega)}^{2}+\|v_{n}\|_{L^{2}(\Omega)}^{2}\leq
O(1)(M^{-s+2}+(N_{b}^{-\frac{1}{2}}\nu_{d-1}+N^{-\frac{1}{2}}\nu_{d}))+I_{GP},$
(47)

where $\nu_{d-1},\nu_{d}$ are standard normal random variables in the $d-1$
and $d$ dimensions, $O(1)$ may depend on $d,s,|\Omega|$,
$\|u\|_{W^{s,\infty}(\Omega)},\|u_{n}\|_{W^{s,\infty}(\Omega)}$ and
$\|v_{n}\|_{W^{s-2,\infty}(\Omega)}$, but stay bounded even when $M,N$ and
$N_{b}$ increase.

######  Remark 1.

The error $I_{GP}$ in (28) depends on the algorithm used to train the neural
networks. Suppose in each of the training step, we conduct the Gradient
Descent/Ascent until reaching the local minimum/maximum, i.e.

$u_{n+1}=\inf_{w\in\mathcal{U}_{\theta}}\hat{\mathcal{L}}(w,v_{n})\quad\text{and}\quad
v_{n+1}=\sup_{v\in\mathcal{V}_{\tau}}\hat{\mathcal{L}}(u_{n},v),$

then $I_{GP}$ can be expressed as

$\displaystyle
I_{GP}=\hat{\mathcal{L}}(u_{n+1},v_{n})-\hat{\mathcal{L}}(u_{n},v_{n+1}).$

Hence, if the training algorithm converges,
$(u_{n},v_{n})\to(\hat{u},\hat{v})$ as $n\to\infty$, the above error will
converge to zero. Note that in practice $\hat{\mathcal{L}}(u_{n},v_{n})$ may
not converge to zero, still $I_{GP}$ provides a good quantification on the
training error.

If $\hat{L}(\theta,\tau)$ is assumed to be convex in $\theta$ and concave in
$\tau$, it is proved in [38] that with a single step Gradient Descent-Ascend
(GDA) the training loss gap $I_{GP}$ after $n$ training steps satisfies

$\displaystyle
I_{GP}=\max_{\tau}\hat{L}(\theta_{n},\tau)-\min_{\theta}\hat{L}(\theta,\tau_{n})\leq
C\frac{|\theta^{*}-\theta_{0}|^{2}+|\tau^{*}-\tau_{0}|^{2}}{\eta n},$ (48)

where $C$ is an absolute constant, $\theta_{0},\tau_{0}$ are the initial
parameters, $\eta$ is the step size and $(\theta^{*},\tau^{*})$ is the saddle
point of $\hat{L}(\theta,\tau)$. Therefore, for any small $\epsilon>0$, there
exists $n$ such that $I_{GP}<\epsilon$ after $n$ training steps. For non-
convex $\hat{L}$, it is hard to obtain a sharp convergence rate like (48).
Nonetheless, if $\hat{L}(\theta,\tau)$ is non-convex in $\theta$ but concave
in $\tau$, [29] shows the convergence to a set where
$\nabla_{\theta,\tau}\hat{L}$ becomes small after finite training steps.

For the the multistep GDA approach, convergence rates are established in [39].
For general non-convex non-concave case, convergence analysis is rather
difficult, and we refer to [46] for related discussions.

######  Remark 2.

Based on the estimate (47), it is evident that the neural network
approximation error (comprising the first term on the right) diminishes with
the increasing width of the neural network. Simultaneously, the Monte-Carlo
sampling error reduces with the sampling size $N_{b}$ and $N$. The anticipated
trend is a decrease in the gap of the training loss with the progression of
training iterations.

##  5\. Experiments

We apply our Inf-SupNet algorithm to a set of elliptic problems posed on a
bounded domain. Monte Carlo methods are employed to generate training data
points. To assess the algorithm’s accuracy, we compute the relative $L^{2}$
error $\|u_{\theta}-u^{*}\|_{L^{2}}/\|u^{*}\|_{L^{2}}$, where $u_{\theta}$ is
the computed solution and $u^{*}$ is the true solution $u^{*}$. The norms are
computed through Monte Carlo integration on test data points, separate from
the training set.

The solution $u$ and the Lagrangian multiplier $v$ are both represented using
fully connected neural networks with a $\tanh$ activation function. Training
the Inf-SupNet involved the use of the _RMSprop_ gradient descent algorithm
[54]. To ensure rapid and accurate training, we set the learning rate to
$0.001$ and decay it by $0.5$ every $5000$ training epochs.

###  5.1. Poisson equation with Dirichlet boundary conditions

We first show the capacity of our method in solving high dimensional Poisson
equations. We first consider the domain $\Omega=[-1,1]^{d}$ and examine the
Poisson equation with homogeneous Dirichlet boundary conditions, i.e. $g(x)=0$
for $x\in\partial\Omega$. We consider the source term

$\displaystyle f(x)=\prod_{i=1}^{d}\cos\left(\frac{\pi}{2}x_{i}\right).$

with $d=5$ and use a neural network with 100 neurons in each hidden layer and
4 hidden layers to represent both the $\mathcal{U}$ and $\mathcal{V}$
networks. We use _RMSprop_ with learning rate $0.0005$ and set the learning
rate to scale by 0.5 with a step size 2000. Figure 1 shows that the training
error reaches approximately around 0.6% after 20,000 iterations. Moreover,
according to the bottom left figure, the point-wise error is less than 0.01.
This demonstrates that our algorithm provides an accurate numerical solutions
even in high dimensions. Moreover, the top right figure illustrates that the
Lagrangian multiplier network output $v_{\tau}$ is small, which agrees with
the theoretical analysis.

Examining the loss curve in Figure 1, it is evident that the loss does not
exhibit monotonic decay throughout the training iterations. This behavior is
inherent in the min-max optimizations, where the gap between the training loss
during gradient descent and ascent that converges to zero, rather than the
loss itself. However, although the loss may increase, the error consistently
decreases until the loss stablizes.

Figure 1. Numerical results for the Poisson equation with Dirichlet BC. Top
left: the true solution $u^{*}$, top middle: the computed solution
$u_{\theta}$; top right: the Lagrangian multiplier $v_{\tau}$; bottom left:
the difference $u_{\theta}-u^{*}$; bottom middle the training loss over
iterations; bottom right: the test error (relative $L^{2}$ error with respect
to the true solution) over iterations. The solutions are plotted over
$x_{2},x_{4}$ with $x_{1}=x_{3}=x_{5}=0$. The EMAs of the loss and error are
plotted at bottom middle and right.

To show the convergence of the algorithm, we compute the relative error with
true solutions for different numbers of sampling points with $d=3$, and the
relative error, presenting the results in Table 1. The table shows that the
error decreases as the number of sampling points increases. The number of
sampling points not only influences the Monte-Carlo integration error, but
also impacts the training error. Hence, to get a more accurate numerical
solution, we can increase the number of sampling points. However, this comes
at the expense of additional computational costs.

Table 1. Error table for the Poisson equation for different number of sampling points(d=3) $N_{in}$ | $N_{bd}$ | relative $L^{2}$ error
---|---|---
2500 | 300 | $9.1638\times 10^{-3}$
5000 | 432 | $4.3924\times 10^{-3}$
7500 | 516 | $3.6575\times 10^{-3}$
10000 | 600 | $2.6667\times 10^{-3}$

The selection of the neural network structure is crucial for algorithm
implementation. Table 2 presents the error corresponding to various widths and
depths (where depth represents the number of hidden layers). Elevating both
the width and depth enhances the network’s capability to approximate the true
solution, subsequently reducing the approximation error and yielding a closer
alignment with the true solution.

Table 2. Error table for the Poisson equation for different neural network structures (d=3) width | error (depth=1) | error (depth=2) | error (depth=3)
---|---|---|---
10 | $5.7629\times 10^{-2}$ | $5.7293\times 10^{-2}$ | $2.3362\times 10^{-2}$
20 | $1.6861\times 10^{-2}$ | $8.4221\times 10^{-3}$ | $5.3472\times 10^{-3}$
40 | $5.9664\times 10^{-3}$ | $2.6667\times 10^{-3}$ | $2.3664\times 10^{-3}$

###  5.2. Irregular domains

We consider the domain $\Omega=[-1,1]^{2}\setminus[0,1]^{2}$ and examine the
Poisson equation with homogeneous Dirichlet boundary conditions, where

$\displaystyle f(x)=e^{-(x+\frac{1}{2})^{2}-(y+\frac{1}{2})^{2}}.$

Taking $4$ hidden layers, each with $40$ neurons in both the $\mathcal{U}$ and
$\mathcal{V}$ networks, we approximate the solution since no analytical
solution is available. Also, we use the numerical solution through finite
element methods as as a reference for the true solution $u^{*}$. Figure 2
illustrates that the relative $L^{2}$ error reaches approximately $0.01$ after
$20,000$ training iterations, with pointwise errors being less than $0.01$.
This observation shows that our algorithms work well in handling irregular
domains.

Figure 2. Numerical results on L-shaped domain. Left: solution computed by
InfSupNet; middle: pointwise error with the FEM solution; Right: relative
$L^{2}$ error during training

###  5.3. Mixed boundary conditions

We consider the domain $\Omega=[0,1]\times[0,1]$ and the boundary set
$\Gamma_{D}=\\{(0,y)\cup(1,y)\in\partial\Omega\\}$,
$\Gamma_{N}=\\{(x,0)\cup(x,1)\in\partial\Omega\\}$, with boundary conditions

$\displaystyle u=0\quad\text{ on }\Gamma_{D},$ $\displaystyle\frac{\partial
u}{\partial n}=g\quad\text{ on }\Gamma_{N},$

where $g=\sin(5x)$. The source term in the equation
$f=10\exp(-((x-0.5)^{2}+(y-0.5)^{2})/0.02)$. The results shown in Figure 3
affirm the effectiveness of our approach in adeptly managing mixed boundary
conditions.

Figure 3. Numerical results for mixed boundary conditions. Left: solution
computed by InfSupNet; middle: pointwise error with the FEM solution; Right:
relative $L^{2}$ error during training

###  5.4. Semilinear elliptic PDE

We address the semilinear elliptic PDE giving by

$\displaystyle-\nabla\cdot(a(x)\nabla u)=f,\quad\text{ in }\Omega,$
$\displaystyle u=0,\quad\text{ on }\partial\Omega,$

where $a(x)=\sum_{i=1}^{d}x_{i}$, $f=f(x)=-(2d+2)\sum_{i=1}^{d}x_{i}$, and the
true solution is $u^{*}(x)=1+|x|^{2},$ with $d=5$. We randomly select 10,000
sampling points inside $\Omega$ and 400 points on the boundary
$\partial\Omega$. During training, we use a fixed learning rate of 0.0002. The
relative $L^{2}$ error between the computed solution and the true solution
$u^{*}$ is plotted in Figure 4, revealing that the relative error reaches
approximately $0.01$ after $10000$ training iterations.

###  5.5. Nonlinear elliptic PDE

We solve the nonlinear PDE

$\displaystyle-\nabla\cdot(q(u)\nabla u(x))=0,\quad\text{in }\Omega,$
$\displaystyle u=0,\quad\text{on }\\{x_{1}=0\\}\cap\partial\Omega,$
$\displaystyle u=1,\quad\text{on }\\{x_{1}=1\\}\cap\partial\Omega,$
$\displaystyle\frac{\partial u}{\partial n}=0,\text{ on
}\partial\Omega\backslash\\{x_{1}=0\text{ or }1\\},$

where, $q(u)=(1+u)^{m}$. The true solution
$u(x)=((2^{m+1}-1)x_{1}+1)^{\frac{1}{m+1}}-1$. We consider $m=2$ and $d=5$.
The relative $L^{2}$ error between the computed solution and the true solution
$u^{*}$ is plotted in Figure 4, showing that the relative error reaches
approximately $0.01$ after $10000$ training iterations.

Figure 4. Numerical results for semilinear and nonlinear elliptic PDEs with
$d=5$. Left: Relative error with iterations for the semilinear PDE; right:
Relative error with iterations for the nonlinear PDE

##  6\. Concluding remarks

We introduce Inf-SupNet, a novel method for learning solutions to elliptic PDE
problems employing two deep neural networks. Our approach is grounded in two
fundamental concepts. Firstly, we establish an inf-sup reformulation of the
underlying PDE problem, incorporating both the unknown variable and a
Lagrangian multiplier. Secondly, we address the saddle point problem by
leveraging two deep neural networks. Theoretically, we demonstrate that the
global approximation error is governed by three factors: the training error,
the Monte-Carlo sampling error, and the universal approximation error inherent
in neural networks.

In our experimental evaluations, we showcase the stability and accuracy of the
method when applied to high-dimensional elliptic problems on both regular and
irregular domains. However, certain limitations of Inf-SupNet necessitate
careful consideration. Firstly, while the method is robust with guaranteed
error bounds, it may not necessarily outperform some state-of-the-art methods.
Secondly, akin to most min-max problems, the selection of the optimizer is
problem-dependent. Further investigation of training error would be highly
beneficial. Thirdly, the error bounds we obtained are based on a posterior
estimate, deriving a prior estimate would enhance the quantification of error.
Future endeavours may focus on improving the efficiency of Inf-SupNet and
extending its applicability to a broader class of time-dependent problems.

## Acknowledgements

This research was partially supported by the National Science Foundation under
Grant DMS1812666. We would like to thank Marius Zeinhofer in pointing out a
mistake in the initial draft.

## Appendix

We present a lemma along with its proof.

######  Lemma 4.

Let $(u^{*},v^{*})\in H^{2}(\Omega)\times L^{2}(\Omega)$. Then $(u^{*},v^{*})$
is a saddle point of $\mathcal{L}(u,v)$ (defined in (5)) if and only if the
strong min-max property (11) holds.

###### Proof.

First, assume $(u^{*},v^{*})$ is a saddle point of $\mathcal{L}(u,v)$, we
prove (11). For any $u,q$ we have

$\displaystyle\mathcal{L}(u,q)\leq\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v),\text{ for any }q\in L^{2}(\Omega).$

Taking the infimum gives

$\displaystyle\inf_{u\in H^{2}(\Omega)}\mathcal{L}(u,q)\leq\inf_{u\in
H^{2}(\Omega)}\sup_{v\in L^{2}(\Omega)}\mathcal{L}(u,v)$

holds for any $q\in L^{2}(\Omega)$. Hence

$\displaystyle\sup_{q\in L^{2}(\Omega)}\inf_{u\in
H^{2}(\Omega)}\mathcal{L}(u,q)\leq\inf_{u\in H^{2}(\Omega)}\sup_{v\in
L^{2}(\Omega)}\mathcal{L}(u,v).$ (49)

From (10), we get

$\displaystyle\sup_{q\in
L^{2}(\Omega)}\mathcal{L}(u^{*},q)\leq\mathcal{L}(u^{*},v^{*})\leq\inf_{w\in
H^{2}(\Omega)}\mathcal{L}(w,v^{*}).$

Therefore,

$\displaystyle\inf_{u\in H^{2}(\Omega)}\sup_{v\in
L^{2}(\Omega}\mathcal{L}(u,v)$ $\displaystyle\leq\sup_{q\in
L^{2}(\Omega)}\mathcal{L}(u^{*},q)\leq\mathcal{L}(u^{*},v^{*})\leq\inf_{w\in
H^{2}(\Omega)}\mathcal{L}(w,v^{*})$ $\displaystyle\leq\sup_{q\in
L^{2}(\Omega)}\inf_{w\in H^{2}(\Omega)}\mathcal{L}(w,q),$

which together with (49) implies (11).

Next, assume the strong min-max property (11) holds, $u^{*}$ is the optimal
solution to the primal problem and $v^{*}$ is the optimal solution to the dual
problem. From (14) and (17) we can see that $u^{*}$ satisfies the equation
$-\Delta u^{*}-f=0$, and $v^{*}=0$. Therefore, for any $u\in H^{2}(\Omega)$,

$\displaystyle\mathcal{L}(u^{*},v^{*})=\frac{1}{2}\|u^{*}-g\|_{L^{2}(\partial\Omega)}^{2}=\inf_{u\in
H^{2}(\Omega)}\frac{1}{2}\|u-g\|_{L^{2}(\partial\Omega)}^{2}\leq\frac{1}{2}\|u-g\|_{L^{2}(\Omega)}=L(u,v^{*}).$

And for any $v\in L^{2}(\Omega)$:

$\mathcal{L}(u^{*},v^{*})=\frac{1}{2}\|u^{*}-g\|_{L^{2}(\partial\Omega)}^{2}+(-\Delta
u^{*}-f,v)\geq\mathcal{L}(u^{*},v).$

These two inequalities imply that $(u^{*},v^{*})$ is a saddle point of
$\mathcal{L}$. ∎

Here, we present a quantitative universal approximation result; refer to [12]
for a comprehensive proof.

######  Proposition 1 ([12]).

For $h\in W^{s,\infty}(\Omega)$, and any $M\geq 5d^{2}$, there exists a
feedforward neural network $\mathcal{U}$ with two hidden layers, having a
width given by

$\max\\{3\left[\frac{s}{2}\right]\binom{s+d-1}{s-1}+d(M-1),3\left[\frac{d+2}{2}\right]5^{d}M^{d}\\},$

such that

$\displaystyle\inf_{\hat{h}\in\mathcal{U}}\|h-\hat{h}\|_{L^{\infty}(\Omega)}\leq\frac{C_{1}(d,s,|\Omega|)}{M^{s}}\|h\|_{W^{s,\infty}(\Omega)}$
(50)

for some constant $C_{1}=C_{1}(d,s,|\Omega|)>0$, and for $k=1,\ldots,s-1$, and

$\displaystyle\inf_{\hat{h}\in\mathcal{U}}\|h-\hat{h}\|_{W^{k,\infty}(\Omega)}\leq
C_{2}(d,k,s,|\Omega|)\frac{\max\\{R_{k}^{k},\log^{k}(\beta
M^{s+d+2})\\}}{M^{s-k}}\|h\|_{W^{s,\infty}(\Omega)},$ (51)

where $C_{2}=C_{2}(d,k,s,|\Omega|)>0$ is a positive constant only depending on
$d,k,s,|\Omega|$,

$\displaystyle\beta=C(d,k,s)\frac{\max\\{1,\|h\|_{W^{k,\infty}(\Omega)}^{\frac{1}{2}}\\}}{\min\\{1,{\|h\|_{W^{s,\infty}(\Omega)}}^{\frac{1}{2}}\\}},$
(52)

with $C=C(d,k,s)$ being a positive constant depending on $d,k,s$, and $R_{k}$
is the lowest value such that $\tanh^{(m)}(x)$ is decreasing on
$[R_{k},\infty)$ for all $1\leq m\leq k$.

## References

  * [1] Andrew R Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information theory, 39(3):930–945, 1993.
  * [2] Christian Beck, Weinan E, and Arnulf Jentzen.  Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations.  Journal of Nonlinear Science, 29:1563–1619, 2019.
  * [3] Jens Berg and Kaj Nyström.  A unified deep artificial neural network approach to partial differential equations in complex geometries.  Neurocomputing, 317:28–41, 2018.
  * [4] Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik.  Prediction of aerodynamic flow fields using convolutional neural networks.  Computational Mechanics, 64:525–545, 2019.
  * [5] Daniele Boffi, Franco Brezzi, and Michel Fortin.  Mixed finite element methods and applications, volume 44.  Springer, 2013.
  * [6] Russel E Caflisch.  Monte-Carlo and quasi-Monte-Carlo methods.  Acta Numerica, 7:1–49, 1998.
  * [7] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (PINNs) for fluid mechanics: A review.  Acta Mechanica Sinica, 37(12):1727–1738, 2021.
  * [8] Zhiqiang Cai, Jingshuang Chen, Min Liu, and Xinyu Liu.  Deep least-squares methods: An unsupervised learning-based numerical method for solving elliptic PDEs.  Journal of Computational Physics, 420:109707, 2020.
  * [9] Tianping Chen and Hong Chen.  Approximations of continuous functionals by neural networks with application to dynamic systems.  IEEE Transactions on Neural networks, 4(6):910–918, 1993.
  * [10] Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli.  Scientific machine learning through physics–informed neural networks: Where we are and what’s next.  Journal of Scientific Computing, 92(3):88, 2022.
  * [11] George Cybenko.  Approximation by superpositions of a sigmoidal function.  Mathematics of control, signals and systems, 2(4):303–314, 1989\.
  * [12] Tim De Ryck, Samuel Lanthaler, and Siddhartha Mishra.  On the approximation of functions by tanh neural networks.  Neural Networks, 143:732–750, 2021.
  * [13] Tim De Ryck, Siddhartha Mishra, and Roberto Molinaro.  Weak physics informed neural networks for approximating entropy solutions of hyperbolic conservation laws.  In Seminar für Angewandte Mathematik, Eidgenössische Technische Hochschule, Zürich, Switzerland, Rep, volume 35, page 2022, 2022\.
  * [14] MWMG Dissanayake and Nhan Phan-Thien.  Neural-network-based approximations for solving partial differential equations.  Communications in Numerical Methods in Engineering, 10(3):195–201, 1994.
  * [15] Weinan E, Jiequn Han, and Arnulf Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations.  Communications in Mathematics and Statistics, 5(4):349–380, 2017\.
  * [16] Pierre Grisvard.  Elliptic problems in nonsmooth domains.  SIAM, 2011.
  * [17] Ingo Gühring, Gitta Kutyniok, and Philipp Petersen.  Error bounds for approximations with deep ReLU neural networks in $W^{s,p}$ norms.  Analysis and Applications, 18(05):803–859, 2020.
  * [18] Ehsan Haghighat, Maziar Raissi, Adrian Moure, Hector Gomez, and Ruben Juanes.  A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics.  Computer Methods in Applied Mechanics and Engineering, 379:113741, 2021.
  * [19] Yuling Jiao, Yanming Lai, Yisu Lo, Yang Wang, and Yunfei Yang.  Error analysis of deep Ritz methods for elliptic equations.  Analysis and Applications, 22(1):57–87, 2024.
  * [20] Biswajit Khara, Aditya Balu, Ameya Joshi, Soumik Sarkar, Chinmay Hegde, Adarsh Krishnamurthy, and Baskar Ganapathysubramanian.  Neufenet: Neural finite element solutions with theoretical bounds for parametric PDEs.  arXiv preprint arXiv:2110.01601, 2021.
  * [21] Biswajit Khara, Ethan Herron, Zhanhong Jiang, Aditya Balu, Chih-Hsuan Yang, Kumar Saurabh, Anushrut Jignasu, Soumik Sarkar, Chinmay Hegde, Adarsh Krishnamurthy, and Baskar Ganapathysubramanian.  Neural PDE solvers for irregular domains.  arXiv preprint arXiv:2211.03241, 2022.
  * [22] Ehsan Kharazmi, Zhongqiang Zhang, and George Em Karniadakis.  Variational physics-informed neural networks for solving partial differential equations.  arXiv preprint arXiv:1912.00873, 2019.
  * [23] Reza Khodayi-Mehr and Michael Zavlanos.  Varnet: Variational neural networks for the solution of partial differential equations.  In Learning for Dynamics and Control, pages 298–307. PMLR, 2020\.
  * [24] Yuehaw Khoo, Jianfeng Lu, and Lexing Ying.  Solving parametric PDE problems with artificial neural networks.  European Journal of Applied Mathematics, 32(3):421–435, 2021.
  * [25] Yuehaw Khoo, Jianfeng Lu, and Lexing Ying.  Solving parametric PDE problems with artificial neural networks.  European Journal of Applied Mathematics, 32(3):421–435, 2021.
  * [26] Gitta Kutyniok, Philipp Petersen, Mones Raslan, and Reinhold Schneider.  A theoretical analysis of deep neural networks and parametric PDEs.  Constructive Approximation, 55(1):73–125, 2022.
  * [27] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [28] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  International Conference on Learning Representations, 2021.
  * [29] Tianyi Lin, Chi Jin, and Michael Jordan.  On gradient descent ascent for nonconvex-concave minimax problems.  In International Conference on Machine Learning, pages 6083–6093. PMLR, 2020.
  * [30] Siting Liu, Stanley Osher, Wuchen Li, and Chi-Wang Shu.  A primal-dual approach for solving conservation laws with implicit in time approximations.  Journal of Computational Physics, 472:111654, 2023.
  * [31] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature Machine Intelligence, 3(3):218–229, 2021.
  * [32] Mohammadamin Mahmoudabadbozchelou, George Em Karniadakis, and Safa Jamali.  nn-PINNs: Non-newtonian physics-informed neural networks for complex fluid modeling.  Soft Matter, 18(1):172–185, 2022.
  * [33] Zhiping Mao, Ameya D Jagtap, and George Em Karniadakis.  Physics-informed neural networks for high-speed flows.  Computer Methods in Applied Mechanics and Engineering, 360:112789, 2020.
  * [34] Sandra May, Rolf Rannacher, and Boris Vexler.  Error analysis for a finite element approximation of elliptic dirichlet boundary control problems.  SIAM Journal on Control and Optimization, 51(3):2585–2611, 2013\.
  * [35] Piotr Minakowski and Thomas Richter.  A priori and a posteriori error estimates for the Deep Ritz method applied to the Laplace and Stokes problem.  Journal of Computational and Applied Mathematics, 421:114845, 2023\.
  * [36] Siddhartha Mishra and Roberto Molinaro.  Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs.  IMA Journal of Numerical Analysis, 42(2):981–1022, 2022.
  * [37] Siddhartha Mishra and Roberto Molinaro.  Estimates on the generalization error of physics-informed neural networks for approximating PDEs.  IMA Journal of Numerical Analysis, 43(1):1–43, 2023.
  * [38] Arkadi Nemirovski.  Prox-method with rate of convergence $o(1/t)$ for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems.  SIAM Journal on Optimization, 15(1):229–251, 2004.
  * [39] Maher Nouiehed, Maziar Sanjabi, Tianjian Huang, Jason D Lee, and Meisam Razaviyayn.  Solving a class of non-convex min-max games using iterative first order methods.  Advances in Neural Information Processing Systems, 32, 2019.
  * [40] Houman Owhadi.  Bayesian numerical homogenization.  Multiscale Modeling & Simulation, 13(3):812–828, 2015.
  * [41] Allan Pinkus.  Approximation theory of the MLP model in neural networks.  Acta Numerica, 8:143–195, 1999.
  * [42] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707, 2019.
  * [43] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Inferring solutions of differential equations using noisy multi-fidelity data.  Journal of Computational Physics, 335:736–746, 2017.
  * [44] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Machine learning of linear differential equations using Gaussian processes.  Journal of Computational Physics, 348:683–693, 2017.
  * [45] Majid Rasht-Behesht, Christian Huber, Khemraj Shukla, and George Em Karniadakis.  Physics-informed neural networks (PINNs) for wave propagation and full waveform inversions.  Journal of Geophysical Research: Solid Earth, 127(5):e2021JB023120, 2022.
  * [46] Meisam Razaviyayn, Tianjian Huang, Songtao Lu, Maher Nouiehed, Maziar Sanjabi, and Mingyi Hong.  Nonconvex min-max optimization: Applications, challenges, and recent theoretical advances.  IEEE Signal Processing Magazine, 37(5):55–66, 2020.
  * [47] Jon A Rivera, Jamie M Taylor, Ángel J Omella, and David Pardo.  On quadrature rules for solving partial differential equations using neural networks.  Computer Methods in Applied Mechanics and Engineering, 393:114710, 2022.
  * [48] Christian P Robert, George Casella, and George Casella.  Monte Carlo statistical methods, volume 2.  Springer, 1999.
  * [49] Keith Rudd and Silvia Ferrari.  A constrained integration (Ritz) approach to solving partial differential equations using artificial neural networks.  Neurocomputing, 155:277–285, 2015.
  * [50] Martin Schechter.  On $L^{p}$ estimates and regularity II.  Mathematica Scandinavica, 13.1: 47–69, 1963.
  * [51] Yeonjong Shin, Zhongqiang Zhang, and George Em Karniadakis.  Error estimates of residual minimization using neural networks for linear PDEs.  Journal of Machine Learning for Modeling and Computing, 4(4), 2023\.
  * [52] Jonathan W Siegel.  Optimal approximation rates for deep ReLU neural networks on Sobolev and Besov spaces.  Journal of Machine Learning Research, 24(357):1–52, 2023.
  * [53] Justin Sirignano and Konstantinos Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  Journal of computational physics, 375:1339–1364, 2018.
  * [54] Tijmen Tieleman and Geoffrey Hinton.  RMSPROP: divide the gradient by a running average of its recent magnitude. Coursera: Neural networks for machine learning.  COURSERA Neural Networks Mach. Learn, 17, 2012.
  * [55] Carlos Uriarte, David Pardo, Ignacio Muga, and Judit Muñoz-Matute.  A deep double deeponet method (D2RM) for solving partial differential equations using neural networks.  Computer Methods in Applied Mechanics and Engineering, 405:115892, 2023.
  * [56] Liu Yang, Xuhui Meng, and George Em Karniadakis.  B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data.  Journal of Computational Physics, 425:109913, 2021.
  * [57] Bing Yu et al.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics, 6(1):1–12, 2018.
  * [58] Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis.  Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems.  Computer Methods in Applied Mechanics and Engineering, 393:114823, 2022.
  * [59] Yaohua Zang, Gang Bao, Xiaojing Ye, and Haomin Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  Journal of Computational Physics, 411:109409, 2020.
  * [60] Marius Zeinhofer, Rami Masri, and Kent-André Mardal.  A unified framework for the error analysis of physics-informed neural networks. arXiv preprint, arXiv:2311.00529, 2023.