# Convergence and error control of consistent PINNs for elliptic PDEs

Andrea Bonito, Ronald DeVore, Guergana Petrova, and Jonathan W. Siegel

###### Abstract

We provide an a priori analysis of a certain class of numerical methods,
commonly referred to as collocation methods, for solving elliptic boundary
value problems. They begin with information in the form of point values of the
right side $f$ of such equations and point values of the boundary function $g$
and utilize only this information to numerically approximate the solution $u$
of the Partial Differential Equation (PDE). For such a method to provide an
approximation to $u$ with guaranteed error bounds, additional assumptions on
$f$ and $g$, called model class assumptions, are needed. We determine the best
error (in the energy norm) of approximating $u$, in terms of the number of
point samples $m$, under all Besov class model assumptions for the right hand
side $f$ and boundary $g$.

We then turn to the study of numerical procedures and asks whether a proposed
numerical procedure (nearly) achieves the optimal recovery error. We analyze
numerical methods which generate the numerical approximation to $u$ by
minimizing a specified data driven loss function over a set $\Sigma$ which is
either a finite dimensional linear space, or more generally, a finite
dimensional manifold. We show that the success of such a procedure depends
critically on choosing a correct data driven loss function that is consistent
with the PDE and provides sharp error control. Based on this analysis a loss
function ${\cal L}^{*}$ is proposed.

We also address the recent methods of Physics Informed Neural Networks
(PINNs). Minimization of the new loss ${\cal L}^{*}$ over neural network
spaces $\Sigma$ is referred to as consistent PINNs (CPINNs). We prove that
CPINNs provides an optimal recovery of the solution $u$, provided that the
optimization problem can be numerically executed and $\Sigma$ has sufficient
approximation capabilities. Finally, numerical examples illustrating the
benefits of the CPINNs are given.

##  1 Introduction

This paper is concerned with numerical methods for solving elliptic partial
differential equations (PDEs). Our primary goal is to provide a rigorous
analysis of rates of convergence for collocation methods, including PINNs
(Physics Informed Neural Networks), for solving such differential equations.

Let $\Omega$ be a bounded domain in the Euclidean space $\mathbb{R}^{d}$,
$d\geq 2$, and $\overline{\Omega}$ be its closure in $\mathbb{R}^{d}$.
Consider the elliptic Dirichlet problem

$\displaystyle-\nabla\cdot(a({\bf x})\nabla u({\bf x}))=f({\bf x}),\quad{\bf
x}\in\Omega,$ $\displaystyle u({\bf z})=g({\bf z}),\quad{\bf
z}\in\partial\Omega,$ (1.1)

where $\partial\Omega$ is the boundary of $\Omega$. In order to prove the
existence and uniqueness of a solution to (1), one needs to impose conditions
on $f,g,a$, and $\Omega$. A standard set of assumptions that guarantee the
existence and uniqueness of a solution (1) is the following: A1: $\Omega$ is a
Lipschitz (graph) domain, i.e., $\Omega$ is an open connected set in
$\mathbb{R}^{d}$ with a Lipschitz boundary in the sense described in [57]; A2:
$f\in H^{-1}(\Omega)$; A3: $g\in H^{1/2}(\partial\Omega)$; A4: the diffusion
coefficient $a$ satisfies the coercivity condition

$0<r\leq a({\bf x})\leq M,\quad{\bf x}\in\Omega,$ (1.2)

for some constants $r,M>0$. Here the spaces $H^{s}(\Omega)$, $s\in\mathbb{R}$,
are the Sobolev spaces of order $s$ in $L_{2}(\Omega)$.

Under the assumptions A1-A4, the Lax-Milgram theorem (see e.g. [25]) implies
that (1) has a unique solution $u\in H^{1}(\Omega)$ which satisfies

$c\\{\|f\|_{H^{-1}(\Omega)}+\|g\|_{H^{1/2}(\partial\Omega)}\\}\leq\|u\|_{H^{1}(\Omega)}\leq
C\\{\|f\|_{H^{-1}(\Omega)}+\|g\|_{H^{1/2}(\partial\Omega)}\\},$ (1.3)

where the constants $c,C$ depend on $\Omega$, $r$ and $M$ 111Throughout this
paper, we use the notation $c,C,c_{1},C_{1}$, etc, to denote constants. Unless
they are absolute constants, we indicate the quantities on which they depend
when they occur. The value of these constants can change at each occurrence.
We also use the notation $A\asymp B$ to denote that $c_{1}A\leq B\leq c_{2}A$
with $c_{1},c_{2}>0$.. Therefore, if we consider the theoretical loss function

${\cal L}_{T}(v):=\|f+\Delta
v\|_{H^{-1}(\Omega)}+\|g-v\|_{H^{1/2}(\partial\Omega)},$ (1.4)

then minimizing this loss over the whole of $H^{1}(\Omega)$ has $u$ as its
unique solution.

In order to make our presentation as clear as possible, we want to avoid the
technicalities of handling general domains and general diffusion coefficients
$a$. Accordingly, we restrict our presentation to the case $\Omega=(0,1)^{d}$
and $a\equiv 1$. This will allow us to concentrate on the novelties of PINNs
and alternative collocation methods. The treatment of more general settings
for the diffusion $a$ and domains $\Omega$ is left to future works. We also
want to point out that assumptions A2-A3 could be replaced with other
smoothness assumptions on $f$ and $g$, respectively, which will result in a
theory similar to the one described in this paper.

###  1.1 Numerical methods for solving (1)

We are interested in numerical methods for solving (1). The most prominent of
these are the Finite Element Methods (FEMs) and the adaptive variations of
these (AFEMs). Over the last decades, a rigorous theory has been developed for
FEMs and AFEMs which guarantees, a priori, the convergence of these numerical
methods and even provides bounds on their rate of convergence under additional
assumptions on $f$ and $g$. These additional assumptions stipulate extra
smoothness of $f$ and $g$ than those assumed in A2-A3.

A numerical method assumes that $f,g$, and $\Omega$ are known and suggests a
procedure for solving (1) based on that information. The numerical procedure
generates a function $\hat{u}$ which is an approximation to $u$. The
efficiency of the numerical method is then evaluated in the following sense.
One chooses a norm $\|\cdot\|_{X}$ in which to measure performance and then
seeks to quantify how fast the error $\|u-\hat{u}\|_{X}$ between the true
solution $u$ and the output $\hat{u}$ of the numerical procedure tends to zero
in terms of the complexity of the numerical procedure. The classical norm
$\|\cdot\|_{X}$ used to measure error is the energy norm which corresponds to
choosing $X=H^{1}(\Omega)$. Other function norms sometimes used to measure
error correspond to $X=L_{p}(\Omega)$, $1\leq p\leq\infty$. We restrict our
analysis of convergence and rates of convergence to the case $X=H^{1}(\Omega)$
in going further.

An important question is how one can compare the performance of different
numerical methods for solving (1) in order to give a fair competition between
all possible numerical methods. This is typically done by assigning some form
of complexity to the numerical method. In the case of FEMs and AFEMs, this
complexity can be measured in terms of the number $n$ of machine operations
used to compute $\hat{u}_{n}=\hat{u}$ and there are theorems that give a
priori bounds for the error $\|u-\hat{u}_{n}\|_{H^{1}(\Omega)}$ under
additional assumptions on the smoothness of $f$ and $g$ (see e.g. [5]). Some
alternatives to measuring the complexity $n$ in terms of machine operations
are the following. If the numerical approximation $\hat{u}$ to $u$ is an
element from a linear $n$-dimensional space $V_{n}$, then the associated
complexity is typically assigned to be the dimension $n$. This may not be
directly converted to an appropriate number of machine operations because
numerical stability becomes an issue. If the approximation $\hat{u}$ comes
from a nonlinear manifold ${\cal M}_{n}$, then one can use the number of
parameters $n$ determining the manifold as a complexity measure. However, in
the latter case, one has to impose extra conditions on the manifold (or the
numerical procedure) to prevent the use of space filling manifolds and thus
avoid unstable numerical methods, see [9]. In the absence of such
restrictions, the numerical method may require an inordinate amount of
computational resources to achieve a desired accuracy even in the case when
the dimension $n$ is small.

Another issue to consider is in what sense $f,g,\Omega$ are known/given to us.
One setting is to assume that we can query (ask questions about) $f,g,\Omega$,
and receive the answer to such queries. In this setting, one has to decide
which type of questions are allowed and whether these queries are answered
free of charge or should the cost of asking/answering such queries be included
in the numerical cost of the algorithm. For example, in FEMs, the data used
are certain linear functionals (inner products with the basis functions of the
Galerkin space) which are then utilized in the FEM approximation.

In this paper, we are interested in collocation methods as described in the
next section. Accordingly, the linear functionals we consider are point
evaluations of $f$ and $g$ at points from $\overline{\Omega}$. In order for
these to make sense, we assume that $f$ and $g$ are continuous functions,
i.e., $f\in C(\overline{\Omega})$, $g\in C(\partial\Omega)$. The accuracy of
how well $\hat{u}$ approximates $u$ will depend on two components. The first
is the number $m$ of point evaluations and their positioning. The optimal
accuracy that any numerical method can achieve from these $m$ pieces of
information under the model class assumption on the function is called the
optimal recovery rate. A central question in this paper is to determine this
optimal recovery rate under various model class assumptions, provided the $m$
data sites are optimally positioned. This is the subject of §3.

A second item in collocation methods is how one numerically uses the $m$
pieces of information provided to create a numerical approximation to $u$. If
$\Sigma_{n}$ is the underlying space used to create $\hat{u}$ then the
accuracy of the numerical recovery will depend on $n$ and the choice of
$\Sigma_{n}$. In other words, $\hat{u}$ can be viewed as
$\hat{u}=\hat{u}_{n,m}$. If one fixes $m$, then one can study how fast the
accuracy of approximation $\|u-\hat{u}_{n,m}\|$ tends to the optimal recovery
rate as $n\to\infty$. This issue is discussed in §8.

Concerning the sense in which we know $\Omega$, one usually considers
polyhedral domains whose vertices and edges are given to us. As already noted,
for simplicity, we assume that $\Omega=(0,1)^{d}$. However, with some
additional technicality, we could equally well start with polyhedral domain.
The standard way of handling more general domains $\Omega$ is to first find a
polyhedral domain $\widehat{\Omega}$ that approximates $\Omega$, solve the
corresponding PDE on $\widehat{\Omega}$, and then analyze the effect of the
approximation of $\Omega$ by $\widehat{\Omega}$.

###  1.2 Collocation methods and PINNs

Recently, there has been significant interest in using neural networks (NNs)
as a nonlinear manifold to generate the approximation $\hat{u}$ to the
solution $u$ of (1). Let $\Sigma_{n}$ denote the set of outputs of a neural
network with $n$ parameters and fixed architecture. Once $\Sigma_{n}$ is
chosen, the numerical procedure queries $f$ and $g$ and then uses this
information to numerically construct a function $\hat{u}\in\Sigma_{n}$ which
serves as an approximation to $u$ in a specified norm $\|\cdot\|_{X}$. The
queries of $f$ and $g$ are taken as point evaluations at specified points from
$\overline{\Omega}$, thereby providing the values

${\bf f}=(f_{1},\ldots,f_{\widetilde{m}}),\quad f_{i}:=f({\bf x}_{i}),\quad
i=1,\dots,\widetilde{m};\quad{\bf g}=(g_{1},\ldots,g_{{\overline{m}}}),\qquad
g_{i}:=g({\bf z}_{i}),\quad i=1,\dots,\overline{m},$ (1.5)

where ${\bf x}_{i}\in\overline{\Omega}$, $i=1,\dots,\widetilde{m}$, and ${\bf
z}_{i}\in\partial\Omega$, $i=1,\dots,\overline{m}$, are the query sites. We
refer to these points as data sites and denote them by

${\cal X}:=\\{{\bf x}_{1},\dots,{\bf x}_{\widetilde{m}}\\},\quad{\cal
Z}:=\\{{\bf z}_{1},\dots,{\bf z}_{\overline{m}}\\}.$

The performance of such a numerical method will depend on the numbers
$\widetilde{m},\overline{m}$, which we refer to as the query budget, and also
on the positioning of these points.

In PINNs, the numerical procedure to find $\hat{u}$ is to search over
$\Sigma_{n}$ and find a $\hat{u}\in\Sigma_{n}$ which ‘fits the data’. The most
frequently used procedure (there are many variants) is to take $\hat{u}$ as
one of the elements of the set

$\hat{u}\,\in\,\mathop{\rm argmin}_{S\in\Sigma_{n}}{\cal L}_{0}(S),$ (1.6)

where

${\cal L}_{0}(S):=\frac{1}{\widetilde{m}}\sum_{i=1}^{\widetilde{m}}[\Delta
S({\bf x}_{i})+f({\bf
x}_{i})]^{2}+\frac{\lambda}{\overline{m}}\sum_{i=1}^{\overline{m}}[S({\bf
z}_{i})-g({\bf z}_{i})]^{2},$ (1.7)

with $\lambda$ being a tuning parameter. For simplicity, throughout the paper
we will set $\lambda=1$ and consider the loss

${\cal L}(S):=\left[\frac{1}{\widetilde{m}}\sum_{i=1}^{\widetilde{m}}[\Delta
S({\bf x}_{i})+f({\bf
x}_{i})]^{2}\right]^{1/2}+\left[\frac{1}{\overline{m}}\sum_{i=1}^{\overline{m}}[S({\bf
z}_{i})-g({\bf z}_{i})]^{2}\right]^{1/2}$ (1.8)

in our analysis. Note that this is equivalent to the square root of the
typical PINNs loss with $\lambda=1$. The appropriateness of the loss ${\cal
L}_{0}$ and how $\widetilde{m},\overline{m},\lambda$ and the data sites
$({\cal X},{\cal Z})$ should be chosen are major issues and one of the focal
points of this paper.

The PINNs numerical procedures fall into the broad class of collocation
methods for solving PDEs, i.e. numerical methods that use only values of $f$
and $g$ at specified data sites ${\cal X}$ and ${\cal Z}$. Such methods were
studied in the last century (see e.g. [1, 34, 31]) but became less popular
with the advent of FEMs. The novelty of PINNs is to use neural networks
instead of polynomials or splines to build the approximation $\hat{u}$ from
the given data.

Some assumptions are necessary for collocation methods to make sense. Firstly,
to ensure that point values of $f$ and $g$ make sense, we assume that $f$ and
$g$ are continuous function on $\overline{\Omega}$ and $\partial\Omega$,
respectively. Secondly, to have any hope of proving a priori guarantees on
convergence or rates of convergence as
$\widetilde{m},\overline{m},n\to\infty$, we need at a minimum that $f$ comes
from a compact subset ${\cal F}$ of $C(\overline{\Omega})$ and $g$ comes from
a compact subset ${\cal G}$ of $C(\partial\Omega)$. This in turn guarantees
that $u$ is in a compact subset ${\cal U}$ of $H^{1}(\Omega)$.

###  1.3 Overview of this paper

The idea of using neural networks to solve PDEs goes back to the last century
[20, 39], and has recently become extremely popular with the introduction of
physics-informed neural networks (PINNs) [50]. Despite the increasing usage of
PINNs for numerically solving a wide range of PDEs (see for example, [36, 6,
10, 6, 43] and the references therein), a satisfactory analysis of the
convergence and performance of these methods has not been given. Some partial
progress towards this goal has recently been made, see for instance [24, 45,
53, 54, 65].

For example, in [53], convergence is proved in the $C(\overline{\Omega})$-norm
for a modified physics-informed neural network for elliptic PDEs under the
assumption that the right hand side $f$ and boundary data $g$ are Hölder
continuous. Under the additional assumption that the outputs of the neural
network satisfy the boundary conditions exactly, convergence is obtained in
$H^{1}(\Omega)$. In this situation, the boundary values are not enforced via
the loss function and must instead be implemented through the neural network
architecture. We remark that although methods for exactly enforcing boundary
conditions with neural networks have been proposed in [39, 40, 3], this
approach appears unable to rigorously handle arbitrary boundary conditions on
general domains.

Further, in [54] and [65], an abstract framework for analyzing methods based
upon residual minimization, such as PINNs, is presented. This framework has
been used in a variety of follow up works to analyze PINNs (see [2, 58, 22]
for instance). However, the current analysis is, to the best of our knowledge,
not able to obtain convergence rates with respect to the number of collocation
points when a fixed (as opposed to random) set of collocation points are used.

An analysis taking into account the rate of convergence with respect to the
number of collocation points has been obtained in [45]. It applies to the heat
equation (instead of elliptic PDEs) under the assumption that the initial data
lies in the Sobolev space $H^{s}(\Omega)$ for sufficiently large $s$ to ensure
that point values are well-defined.

However, so far, none of the existing literature treats the problem of
obtaining convergence rates in terms of both the number of collocation points
and the network size for elliptic PDEs under general Besov regularity
assumptions on the right hand side $f$ and boundary value $g$. A proper a
priori analysis that addresses this problem would answer the following
questions: Q1: Given the model class assumptions on $f$ and $g$, which query
sites should be used and how large would we have to take $\widetilde{m}$ and
$\overline{m}$ to be able to reach a target accuracy $\varepsilon$ for the
error $\|u-\hat{u}\|_{H^{1}(\Omega)}$? Q2: Given a budget of $m$ queries and
given model class assumption on $f$ and $g$, what is the smallest error that
can be achieved in the recovery of $u$? This is called the optimal recovery
accuracy. Q3: If we use a NN space $\Sigma_{n}$ with $n$ parameters to build
the approximation $\hat{u}$ to $u$, then how large do we need to choose $n$ to
achieve near optimal accuracy? Q4: If we use minimization of a data driven
loss function to generate the numerical approximation $\hat{u}$, then what
form should the loss function take so that small values of this loss function
guaranees small values of the solution approximation error, up to the optimal
recovery rate ?

The goal of the present paper is to answer the above questions for all Besov
space model classes and thereby establish such an a priori analysis of
collocation methods and in particular PINNs. In contrast to previous works, we
also consider the error in the $H^{1}(\Omega)$-norm, and treat the situation
where the boundary values are enforced through the loss function. As already
noted, in order to simplify our presentation, we will only consider the case
of Laplace’s equation (1) with $a\equiv 1$ on $\Omega$ and $\Omega=(0,1)^{d}$.
We leave the problem of handling more general settings to future work. We
remark that an analysis of other methodologies for solving (1) using neural
networks, such as the deep Ritz method [60] and finite neuron method [61], can
be found in [42, 41, 56, 47, 23, 21].

As we have already mentioned, any a priori convergence analysis of collocation
methods requires assumptions on the functions $f$ and $g$ in the form that
$f\in{\cal F}$, $g\in{\cal G}$, where ${\cal F}$ and ${\cal G}$ are compact
subsets of the space of continuous functions on the domains
$\overline{\Omega}$ and $\partial\Omega$, respectively. Such compact sets are
typically described by smoothness conditions. In this paper, we measure
smoothness by membership in Besov classes. For this reason, we begin in §2
with the definition of Besov spaces and the properties we will need in the
sequel. That section will also discuss how well functions in Besov classes can
be approximated (in various norms) by piecewise polynomials. All of the
results in that section are known and are therefore presented in a summary
form. However, for completeness of the paper, and for the specific
presentation of these results, we give proofs of the results we need in the
appendix.

We use membership in Besov spaces $B_{q}^{s}(L_{p}(\Omega))$ as the smoothness
we impose on $f$ and $g$. Once such assumptions are placed on $f$ and $g$, we
answer questions Q1-Q2 in §3. Questions of this type are commonly referred to
as optimal recovery (OR) of a function $u$ from data. Since our data is given
by point samples, this is also known as optimal sampling. Our general setting
is not usually addressed in the OR literature because our data is not point
evaluations of the target function $u$, but rather of $f$ and $g$.
Nevertheless, we consider the results of §3 to be for the most part known in
the sense that we are simply piecing together various known results and
techniques such as those in [12, 13, 38, 49]. The paper [59], which studies
optimal sampling with respect the negative Sobolev norm $H^{-1}(\Omega)$ on
Lipschitz domains $\Omega$, comes closest to presenting the optimal recovery
results in the form we need. However, note that the cases $0<p,q<1$ for all
$d$, or $p=1$, $0<q\leq\infty$, $d=2$ are not covered in [59]. Notably, we
observe that the optimal sampling rates for Besov spaces in $H^{-1}(\Omega)$
coincide with the optimal sampling rates in a certain space
$L_{\gamma}(\Omega)$, determined by the Sobolev embedding (see Theorem 3.1
(iv)).

The results of §3 do not give a numerical method for finding $\hat{u}$ from
the given data that provides the error bound
$\|u-\hat{u}\|_{X}\leq\varepsilon$. They simply establish the theoretical
benchmark for the optimal accuracy of any numerical procedure for the recovery
of $u$ based on the point samples of $f$ and $g$. This theoretical analysis
does not involve NNs or any other proposed collocation method. NNs and PINNs
only enter the picture when one wants to transform the theoretical analysis
into a numerical procedure that utilizes NNs to find an approximant to $u$.

We next discuss in §4 the use of classical numerical methods such as Finite
Element Methods (FEMs) for optimal recovery. In §5, we turn to the question of
using optimization to solve the PDE. We observe that minimizing the
theoretical loss ${\cal L}_{T}$ over a space $\Sigma_{n}$ will approximate
well $u$, provided $\Sigma_{n}$ has sufficient approximation power. However,
this is not a numerical algorithm per se because it does not provide a
numerical recipe for evaluating ${\cal L}_{T}(S)$ when $S\in\Sigma_{n}$.

The remainder of this paper discusses possible numerical methods to achieve
the optimal recovery rate, i.e., to solve the PDE to the highest accuracy
possible given only the information $({\bf f},{\bf g})$ (see (1.5)). An
important consequence of the analysis in §5 is that the widely used loss
function ${\cal L}_{0}$ (or ${\cal L}$) may not be an appropriate
discretization of the theoretical loss ${\cal L}_{T}$. We show that the
correct loss function ${\cal L}^{*}:H^{1}(\Omega)\to\mathbb{R}$ to use in
collocation methods in our setting is

${\cal
L}^{*}(v)\\!:=\\!\\!\left[\frac{1}{\widetilde{m}}\sum_{i=1}^{\widetilde{m}}|f({\bf
x}_{i})+\Delta v({\bf
x}_{i})|^{\gamma}\right]^{1/\gamma}\\!\\!\\!\\!+\left[\frac{1}{\overline{m}^{2}}\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{\overline{m}}\frac{|[g-v]({\bf z}_{i})-[g-v]({\bf
z}_{j})|^{2}}{|{\bf z}_{i}-{\bf
z}_{j}|^{2}}\right]^{1/2}\\!\\!\\!+\left[\frac{1}{\overline{m}}\sum_{j=1}^{\overline{m}}|g({\bf
z}_{j})-v({\bf z}_{j})|^{2}\right]^{1/2},$ (1.9)

where $\gamma$ is the smallest number (if this is possible) so that
$L_{\gamma}(\Omega)$ embeds in $H^{-1}(\Omega)$. For example,
$\gamma=\frac{2d}{d+2}$ in the case $d\geq 3$. The case $d=2$ is more
complicated and we will discuss it later in the paper.

The loss ${\cal L}^{*}$ and its properties are given in §5, §6 and §7. A large
component of its development centers on how to discretize $L_{p}$
(quasi-)norms and the $H^{1/2}$ norm. Our analysis on this subject may be
useful in the development of other methods for solving PDEs numerically.

The first term in the loss ${\cal L}^{*}$ are shown to bound $\|f+\Delta
v\|_{H^{-1}(\Omega)}$ and the second and third terms are shown to bound
$\|g-v\|_{H^{1/2}(\partial\Omega)}$. These in turn give a bound for
$\|u-v\|_{H^{1}(\Omega)}$ via (1.3). We go on to show how near minimizers of
the suggested loss ${\cal L}^{*}$ give an a priori bound of the approximation
error, provided that the approximation method $\Sigma_{n}$ has suitable
approximation properties. The a priori guarantee provides a near optimal
recovery whenever $\Sigma_{n}$ provides sufficiently good approximation of the
elements in the solution model class ${\cal U}$. Notice that this new loss
function replaces the first discrete $\ell_{2}$ term in (1.7) by a discrete
$\ell_{\gamma}$ loss for a proper choice of $\gamma$ and the second discrete
$\ell_{2}$ loss (used to match boundary values) by a weighted $\ell_{2}$ loss.

We next develop in §8 the properties needed of a set $\Sigma_{n}$ so that
minimizing the loss ${\cal L}^{*}$ over $\Sigma_{n}$ results in a near optimal
recovery of $u$. It turns out that a certain restricted approximation
property, guaranteed to hold when $\Sigma_{n}$ is sufficiently good at
approximating the elements of solution model class ${\cal U}$, is sufficient.
In particular, we show that this is the case for suitable NN spaces
$\Sigma_{n}$, provided $n$ is large enough.

The remainder of this paper focuses on using minimization over neural network
spaces as the numerical method as is the case in PINNs. As noted above, the
loss ${\cal L}^{*}$ is the correct loss to use in such a minimization since it
is consistent with recovering $u$ in the $H^{1}(\Omega)$. When the loss ${\cal
L}^{*}$ is minimized over a neural network space, we call the numerical
procedure consistent PINNs (CPINNs).

The later sections of the paper analyze the theoretical gains in using CPINNs
over PINNS. We give a priori bounds on how well the solution CPINNs
approximates $u$ and, in particular, we provide a sufficient condition on the
size $n$ of the neural network to guarantee optimal recovery. These bounds are
established under the assumption that the optimization problem is solved (or
approximately solved) when optimizing over the NN space $\Sigma_{n}$. As we
know, proving that a particular numerical method of optimization (such as
gradient decent) converges to a minimizer (or near minimizer) of the loss is a
serious issue in optimization, which we do not address here. However, we do
prove that for any function $v\in H^{1}(\Omega)$, the quantity ${\cal
L}^{*}(v)$ always provides an upper bound of the true error
$\|u-v\|_{H^{1}(\Omega)}$ up to the optimal recovery rate for the solution
model class ${\cal U}$. This means that ${\cal L}^{*}(v)$ can be used as an a
posteriori error estimator for any proposed numerical approximation $v$ to
$u$. This error estimate can be used to check whether the output of a PINNs
optimization achieves the desired accuracy $\varepsilon$. In other words,
although we do not guarantee a priori that PINNs will achieve the accuracy
$\varepsilon$ (because of lack of performance analysis of optimization
procedures like gradient descent or stochastic gradient descent), we do give
an a numerically implementable a posteriori bound on performance which may
serve to guarantee optimality.

In summary, a theorem which guarantees that PINNs provide an approximation to
$u$ in the $\|\cdot\|_{H^{1}(\Omega)}$ norm with a prescribed accuracy
$\varepsilon$ requires not only smoothness assumptions on $f$ and $g$ as
described above, assumptions on the spacing and number of elements in the data
sites ${\cal X}$ and ${\cal Z}$, and assumptions on the size of $n$, i.e., the
complexity of $\Sigma_{n}$, but also an efficient numerical method that
properly chooses a good approximation $\hat{u}$ to $u$ in the norm
$\|\cdot\|_{H^{1}(\Omega)}$ through minimizing the loss ${\cal L}^{*}$ over
$\Sigma_{n}$.

### Acknowledgments

We are grateful to Peter Binev, Albert Cohen, Wolfgang Dahmen, and Jinchao Xu
for many insightful conversations about the material in this paper.

This research was supported by the NSF Grants DMS-2409807 (AB), DMS-2424305
(JWS), CCF-2205004 (JWS), and DMS 2134077 (RD and GP), and the MURI ONR Grant
N00014-20-1-278 (RD, GP, and JWS).

##  2 Besov spaces

We start this section by recalling the definition of Besov spaces and their
properties. We confine this paper to the case of functions defined on the
domain $\Omega=(0,1)^{d}$, $d\geq 2$. For the range $s>0$, and
$0<p,q\leq\infty$, the Besov space $B_{q}^{s}(L_{p}(\Omega))$ is a space of
functions with smoothness of order $s>0$ in $L_{p}(\Omega)$. Here $q$ is a
secondary index that gives a fine gradation of these spaces. The material in
this section is taken for the most part from the papers [18, 16, 17] and the
reader will have to refer to those papers for some of the definitions and
proofs. Let us also mention that the univariate case is covered in the book
[15].

If $r$ is a positive integer and $0<p\leq\infty$ and $f\in L_{p}(\Omega)$, we
define the modulus of smoothness $\omega_{r}(f,\cdot)_{p}$ of $f$ by

$\omega_{r}(f,t)_{p}:=\omega_{r}(f,t,\Omega)_{p}:=\sup_{|{\rm\bf h}|\leq
t}\|\Delta_{\rm\bf h}^{r}(f,\cdot)\|_{L_{p}(\Omega_{r{\rm\bf h}})},\quad t>0,$
(2.1)

where

$\Delta_{\rm\bf
h}^{r}(f,\cdot):=(-1)^{r}\sum_{k=0}^{r}(-1)^{k}\binom{r}{k}f(\cdot+k{\rm\bf
h}),$ (2.2)

is the $r$-th difference of $f$ for ${\rm\bf h}\in\mathbb{R}^{d}$ and
$\Omega_{\rm\bf h}:=\\{{\rm\bf x}\in\Omega:[{\rm\bf x},{\rm\bf x}+{\rm\bf
h}]\subset\Omega\\}$. Here $[{\rm\bf x},{\rm\bf x}+{\rm\bf h}]$, ${\rm\bf
x},{\rm\bf h}\in\mathbb{R}^{d}$, denotes the line segment in $\mathbb{R}^{d}$
between ${\rm\bf x}$ and ${\rm\bf x}+{\rm\bf h}$, and $|{\rm\bf h}|$ denotes
the Euclidean norm of ${\rm\bf h}$. If $s>0$ and $0<p\leq\infty$, then
$B_{q}^{s}(L_{p}(\Omega))$ is defined as the set of all functions in
$L_{p}(\Omega)$ for which

$|f|_{B_{q}^{s}(L_{p}(\Omega))}:=\left[\intop\limits_{0}^{1}[t^{-s}\omega_{r}(f,t)_{p}]^{q}\frac{dt}{t}\right]^{1/q}<\infty,\quad
0<q<\infty,$ (2.3)

where $r$ can be taken as any integer strictly bigger than $s$. When
$q=\infty$ we replace the integral by a supremum in the definition. This is a
(quasi-)semi-norm and we obtain the (quasi-)norm for
$B_{q}^{s}(L_{p}(\Omega))$ by adding $\|f\|_{L_{p}(\Omega)}$ to it. An
equivalent (quasi-)semi-norm is given by

$|f|_{B_{q}^{s}(L_{p}(\Omega))}\asymp\left[\sum_{k=0}^{\infty}[2^{ks}\omega_{r}(f,2^{-k})_{p}]^{q}\right]^{1/q},\quad
0<q<\infty,$ (2.4)

with equivalency constants independent of $f$. This is proved by discretizing
the integral in (2.3) and using the monotonicity of $\omega_{r}$ as a function
of $t$. When $q=\infty$, (2.4) uses the $\ell_{\infty}$ norm in place of the
$\ell_{q}$ norm, i.e.,

$|f|_{B_{\infty}^{s}(L_{p}(\Omega))}\asymp\sup_{k\geq
0}2^{ks}\omega_{r}(f,2^{-k})_{p}.$ (2.5)

While different choices of $r$ result in different (quasi-)semi-norms, the
corresponding Besov (quasi-)norms are equivalent, provided $r>s$. To fix
matters, we define the Besov norm with the value of $r=r(s)$ as the smallest
integer strictly larger than $\max\\{s,1\\}$. It follows that $r(s)$ is always
larger than or equal to $2$. This choice is for notational convenience in the
material that follows.

Let us make some remarks on these spaces. Consider the role of $q$ in this
definition. If $q_{2}<q_{1}$, then $B_{q_{2}}^{s}(L_{p}(\Omega))\subset
B_{q_{1}}^{s}(L_{p}(\Omega))$, i.e., these spaces get smaller as $q$ gets
smaller. Thus, all $B_{q}^{s}(L_{p}(\Omega))$, $q>0$, are contained in
$B_{\infty}^{s}(L_{p}(\Omega))$ once $s$ and $p$ are fixed. The effect of $q$
in the definition of the Besov spaces is subtle.

In this paper, the space $B_{\infty}^{s}(L_{p}(\Omega))$ is especially
important, and accordingly, we use the abbreviated notation

$B_{p}^{s}:=B_{p}^{s}(\Omega):=B_{\infty}^{s}(L_{p}(\Omega))$ (2.6)

in going forward. It follows that a function $f\in B_{p}^{s}(\Omega)$ if and
only if

$\omega_{r}(f,t)_{p}\leq Mt^{s},\quad t>0,$ (2.7)

and the smallest $M$ for which (2.7) is valid is the (quasi-)semi-norm
$|f|_{B_{p}^{s}}$. This space is commonly referred to as (generalized) Hölder
smoothness of order $s$ in $L_{p}$. It is important to note that when $s$ is a
positive integer, we take $r>s$ in its definition. Therefore the Besov space
$B_{p}^{s}$ is not a Lipschitz space when $s$ is a positive integer. In view
of (2.5), we have that a function $f$ is in $B_{p}^{s}$ if and only if

$\omega_{r}(f,2^{-k})_{L_{p}(\Omega)}\leq M^{\prime}2^{-ks},\quad
k=0,1,\dots,$ (2.8)

and the smallest $M^{\prime}$ for which this is true is an equivalent
(quasi-)semi-norm.

###  2.1 The Sobolev spaces $H^{s}(\Omega)$

Since the results of this paper heavily use the Sobolev spaces
$H^{s}(\Omega)$, we give here a short review of these spaces and list some of
their properties. The classical Sobolev spaces $H^{s}(\Omega)$, $s>0$, are the
Besov spaces

$H^{s}(\Omega)=B^{s}_{2}(L_{2}(\Omega))\subset
B^{s}_{\infty}(L_{2}(\Omega))=:B^{s}_{2}(\Omega).$

Two of the $H^{s}(\Omega)$ spaces play an important role in what follows,
namely $H^{-1}(\Omega)$ and $H^{1/2}(\partial\Omega)$. Let us recall the
definition of $H^{-1}(\Omega)$. The space $H_{0}^{1}(\Omega)$ is the
collection of functions in $H^{1}(\Omega)$ which vanish on the boundary of
$\Omega$, i.e., it is the closure of smooth functions in $H^{1}(\Omega)$ which
are identically zero on $\partial\Omega$. The space $H^{-1}(\Omega)$ is by
definition the dual space of $H_{0}^{1}(\Omega)$ and is equipped with the norm

$\|f\|_{H^{-1}(\Omega)}:=\sup_{\|v\|_{H_{0}^{1}(\Omega)}=1}\langle
f,v\rangle_{H^{-1}\times H_{0}^{1}}.$ (2.9)

There are two equivalent ways to describe the space $H^{1/2}(\partial\Omega)$
appearing in A3. The first is as the Besov space
$B_{2}^{1/2}(L_{2}(\partial\Omega))$. Here, one needs the concept of Besov
spaces on manifolds. In our case, the boundary of $\Omega$ is quite simple
since it is the union of the faces of $\Omega$. This definition gives the norm
in $H^{1/2}(\partial\Omega)$ through the modulus of smoothness. Using the
averaged modulus of smoothness, see the Appendix, one finds that the semi-norm
of $g\in H^{1/2}(\partial\Omega)$, is

$|g|^{2}_{H^{1/2}(\partial\Omega)}:=\intop\limits_{\partial\Omega}\intop\limits_{\partial\Omega}\frac{|g({\bf
z})-g({\bf z}^{\prime})|^{2}}{|{\bf z}-{\bf z}^{\prime}|^{d}}\,d{\bf z}\,d{\bf
z}^{\prime},$ (2.10)

which is commonly referred to as the intrinsic semi-norm on this space. We
obtain the norm on this space by adding $\|g\|_{L_{2}(\partial\Omega)}$ to
this semi-norm.

The second way to describe $H^{1/2}(\partial\Omega)$ and its norm is through
the trace operator $Tr=Tr_{\partial\Omega}$. Recall that if $g$ is a
continuous function on $\overline{\Omega}$ then the trace $Tr(g)$ is simply
the restriction of $g$ to $\partial\Omega$. While the trace operator has an
extension to certain functions that are not continuous, in our applications
that appear later the function $g$ will always be continuous. Using the trace
operator one defines $H^{1/2}(\partial\Omega)$ as the trace of functions in
$H^{1}(\Omega)$ and its trace norm is

$\|g\|_{H^{1/2}(\partial\Omega)}:=\inf\\{\|v\|_{H^{1}(\Omega)}:\ Tr(v)=g\\}.$
(2.11)

It is well known that the trace norm and the intrinsic norms for
$H^{1/2}(\Omega)$ are equivalent for Lipschitz domains and in particular for
our case of $\Omega=(0,1)^{d}$ (see [27]). We use both of these norms in what
follows while always making clear which norm is being employed.

###  2.2  Piecewise polynomial approximation and interpolation

Recall that one can characterize membership in Besov spaces by piecewise
polynomial approximation. We describe and prove such characterizations in the
appendix. An important additional fact is that the piecewise polynomials in
such characterizations can be described by interpolation. This allows us to
generate near best piecewise polynomial approximations to $f$ and $g$ by using
only the data $(f_{i})$, respectively $(g_{i})$.

Let us begin with the cube $\Omega=(0,1)^{d}$ and discuss polynomial
interpolation on $\overline{\Omega}$ which we will later rescale to any dyadic
cube. We fix $r\in\mathbb{N}$, $r>1$ and let

$G_{r}:=\left\\{\left(\frac{j_{1}}{r-1},\ldots,\frac{j_{d}}{r-1}\right),\
j_{i}\in\\{0,1,\ldots,r-1\\}\right\\}\subset[0,1]^{d},$

be the tensor product grid of $r^{d}$ equally spaced points in
$\overline{\Omega}$ . Consider the simplicial (Kuhn-Tucker) decomposition of
$\overline{\Omega}$ into simplices $T$,

$\overline{\Omega}=\bigcup\overline{T}.$

Given any one of these simplices $T$, the number $n_{r}$ of grid points in
this simplex $\overline{T}$ is the same as the dimension of ${\cal
P}_{r}:={\cal P}^{d}_{k}$, where ${\cal P}^{d}_{r}$ denotes the linear space
of algebraic polynomials of order $r$ (total degree $r-1$), namely,

${\cal P}^{d}_{r}:=\left\\{\sum_{|{{\bf k}}|_{1}<r}a_{{\bf k}}{\bf x}^{{\bf
k}},\,\,a_{{\bf k}}\in\mathbb{R}\right\\},\quad\hbox{where}\quad{\bf x}^{{\bf
k}}:=x_{1}^{k_{1}}\cdots x_{d}^{k_{d}},\quad{\bf k}:=(k_{1},\dots,k_{d}),\quad
k_{j}\geq 0,\quad|{\bf k}|_{1}:=\sum_{j=1}^{d}k_{j}.$

Let us consider the standard closed simplex

$\overline{T}_{0}:=\\{{\bf
x}=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}:\,x_{1}+\ldots+x_{d}\leq
1,\,\,x_{i}\geq 0,\,i=1,\ldots,d\\}.$

It is well known that polynomial interpolation using elements of ${\cal
P}_{r}$ at the grid points in $G_{r}\cap\overline{T_{0}}$ is well posed (see
e.g. Section 4 in [51] or Chapter 2 in [7]. Let $\phi_{j}:=\phi_{j,T_{0}}$,
$j=1,\dots,n_{r}$, be the Lagrange polynomial basis for ${\cal P}_{r}$
corresponding to the points in $G_{r}\cap\overline{T}_{0}$. Then, the operator

$L_{T_{0}}(f):=\sum_{{\bf x}_{j}\in G_{r}\cap T_{0}}f({\bf x}_{j})\phi_{j},$
(2.12)

is a bounded projector from $C(\overline{\Omega})$ onto ${\cal P}_{r}$, i.e.,

$\|L_{T_{0}}(f)\|_{C(\overline{T_{0}})}\leq\Lambda_{r}\|f\|_{C(\overline{T_{0}})},\quad\hbox{where}\quad\Lambda_{r}:=\|\sum_{j=1}^{n_{r}}|\phi_{j,T_{0}}({{\bf
x}})|\|_{C(\overline{T_{0}})},$ (2.13)

is the Lebesgue constant which depends only on $r$, $d$, and $T_{0}$.

Now, fix $k\geq 0$ and consider the dyadic cubes $I\in{\cal D}_{k}$. If we
rescale to any such dyadic cube $I$ and any simplex $T$ in the simplicial
decomposition of $I$ and denote by $\phi_{j,T}$ the rescaled polynomial
$\phi_{j,T}:=\phi_{j,T_{0}}\circ F_{T}^{-1}$ where
$F_{T}:\overline{T}_{0}\rightarrow\overline{T}$ is an affine map from the
reference simplex $\overline{T}_{0}$ to $\overline{T}$, we obtain the
projector $L_{T}$ onto ${\cal P}_{r}$ which satisfies

$\|L_{T}(f)\|_{C(\overline{T})}\leq\Lambda_{r}(T)\|f\|_{C(\overline{T})},\quad\hbox{where}\quad\Lambda_{r}(T):=\|\sum_{j=1}^{n_{r}}|\phi_{j,T}({{\bf
x}})|\|_{C(\overline{T})},\quad T\subset I\in{\cal D}_{k},\ k\geq 0.$ (2.14)

Notice that $\Lambda_{r}=\Lambda(T)$. Moreover, the interpolating polynomial
$L_{T}(f)$ provides a good approximation to $f$ on $T$. For example, for any
polynomial $P\in{\cal P}_{r}$, we have

$\|f-L_{T}(f)\|_{C(\overline{T})}\leq\|f-P\|_{C(\overline{T})}+\|L_{T}(f-P)\|_{C(\overline{T})}\leq(1+\Lambda_{r})\|f-P\|_{C(\overline{T})}.$
(2.15)

Let us denote by $S^{*}_{k}(f)$ the piecewise polynomial interpolant

$S^{*}_{k}(f):=S^{*}_{k,r}(f):=\sum_{T\in{\cal T}_{k}}L_{T}(f)\chi_{T},$
(2.16)

where ${\cal T}_{k}$ is the collection of all simplicies arising from the
(Kuhn-Tucker) decomposition of the dyadic cubes in ${\cal
D}_{k}(\overline{\Omega})$ and $\chi_{T}$ is the characteristic function of
$T$. The function $S^{*}_{k}(f)$ is continuous on $\overline{\Omega}$ since
the number points from $G_{r}$ that are on a shared face of any two simplicies
from ${\cal T}_{k}$ is the same as the dimension of the polynomial space
${\cal P}^{d-1}_{r}$ .

The following theorem (proved in the appendix) describes the appproximation
accuracy of $S^{*}_{k}$.

######  Theorem 2.1.

Let $s>0$, $0<p\leq\infty$, be fixed with $s>d/p$. Let $r=r(s)$ be the
smallest integer strictly larger than $\max(s,1)$ and let $S_{k}^{*}$ be the
interpolation operator (2.16). Then, for any $f\in
B_{p}^{s}(\Omega)=B_{\infty}^{s}(L_{p}(\Omega))$ and any $\tau\geq p$, we have

$\|f-S_{k}^{*}(f)\|_{L_{\tau}(\Omega)}\leq
C|f|_{B_{p}^{s}(\Omega)}2^{-k(s-d/p+d/\tau)},$ (2.17)

with $C$ independent of $f$ and $k$.

We will also need similar results where the approximation takes place in the
$H^{1}(\Omega)=B_{2}^{1}(L_{2}(\Omega))$ norm.

######  Theorem 2.2.

Let $d\geq 2$. If $f\in B_{q}^{s}(L_{p}(\Omega))$ with $s>d/p$ and $0<p\leq 2$
and $S^{*}_{k}(f)$ is as defined in (2.16), then

$\|f-S_{k}^{*}(f)\|_{H^{1}(\Omega)}\leq
C|f|_{B_{q}^{s}(L_{p}(\Omega))}2^{-k(s-1-d/p+d/2)},\quad k\geq 0.$ (2.18)

There are many results concerning the embeddings of Besov spaces into other
Besov spaces or $L_{p}$ spaces. We are particularly interested in those when a
Besov space is embedded into $C(\Omega)$. For this, we can use the following
theorem proved in the appendix

######  Theorem 2.3.

For every $f\in B_{p}^{s}(\Omega)$, $s>d/p$, $0<p\leq\infty$, there is a
continuous function $\widetilde{f}$ such that $f=\widetilde{f}$ a.e. In fact,
$\widetilde{f}\in$ Lip $\alpha$, with $\alpha:=s-d/p$, i.e.,

$\omega_{r}(\widetilde{f},t)_{C(\overline{\Omega})}\leq
C|f|_{B_{p}^{s}(\Omega)}t^{s-d/p},\quad t>0.$ (2.19)

##  3 Optimal recovery

In this section, we study the question of whether the information $({\bf
f},{\bf g})$ at the fixed data sites ${\cal X}:=\\{{\bf x}_{1},\ldots,{\bf
x}_{\widetilde{m}}\\}$ and ${\cal Z}:=\\{{\bf z}_{1},\ldots,{\bf
z}_{\overline{m}}\\}$ is sufficient to determine $u$ to a prescribed accuracy
$\varepsilon$ in the $H^{1}(\Omega)$ norm. Questions of this type are well
studied and referred to as optimal recovery (OR). The answer to our specific
question depends on the numbers $\widetilde{m},\overline{m}$, the positions of
the points in ${\cal X}$ and ${\cal Z}$, and the assumptions we make on $f$
and $g$. As we have mentioned, we assume that $f$ and $g$ are continuous so
that point evaluation makes sense. It follows that the unknown function $u$ is
continuous as well. If we wish to ensure an error bound on how well we can
recover $u$ in the $H^{1}(\Omega)$ norm, then the conditions we impose on $f$
and $g$ also have to guarantee that $u$ lies in a compact subset of
$H^{1}(\Omega)$. We require that $f,g$ are in certain Besov spaces that
compactly embed into the appropriate space of continuous functions, i.e., $f$
lies in a Besov space that compactly embeds into $C(\Omega)$ and $g$ is in a
Besov space that compactly embeds into $C(\partial\Omega)$.

Before presenting the results on optimal recovery of this section, let us
remark that the results given below are (in essence) all known. For example,
the optimal recovery rates for the Besov model classes of this section have
all been given in the papers [12, 13, 59] save for one small discrepancy on
optimal recovery in $H^{-1}(\Omega)$ (explained below) and for the fact that
they do not directly consider recovery on manifolds like $\partial\Omega$. The
main distinction between the results in [59] and ours given below is in how
the optimal recovery is obtained/proven.The optimal recovery method in [59]
uses linear combinations of certain ‘bump functions’, whereas we use
continuous piecewise polynomial Lagrange interpolation. One can argue that the
fact that optimal recovery can be obtained by interpolation is also known
since it is a prominent method of approximation in Finite Element Methods.
However, the latter community usually considers Sobolev classes, rather than
the more general Besov classes, and also typically restricts $p\geq 1$. In
summary, we could not find a direct reference which proves the results given
below based on continuous piecewise polynomial interpolation. In any case,
those not familiar with optimal recovery may find the results below useful and
the OR community can just do a very light reading of this section. The results
stated in this section are all proved in detail in the appendix.

We assume that $f,g$ belong to the unit ball of a smoothness space that
compactly embeds into the space of continuous functions on the relevant
domain. The function $f$ is defined on the domain
$\Omega\subset\mathbb{R}^{d}$ with continuous extension on
$\overline{\Omega}$. The classical assumptions on function smoothness for a
domain $\Omega$ is membership in a Besov space. Accordingly, we shall assume
that

$f\in{\cal F}:=U({\cal B}),\ \ {\cal B}=B_{q}^{s}(L_{p}(\Omega)),\quad
0<q,p\leq\infty,\quad s>d/p,$ (3.1)

where $U({\cal B})$ is the unit ball of ${\cal B}$ and the restriction on $s$
ensures that ${\cal F}$ compactly embeds into the space of continuous
functions $C(\overline{\Omega})$.

The function $g$ is defined on the manifold $\partial\Omega$. There are two
ways to describe smoothness conditions on such functions $g$. The first is
through the trace operator $Tr:=Tr_{\partial\Omega}$. The second is to place
smoothness conditions directly on $g$ as a function from $C(\partial\Omega)$.
The latter is usually referred to as an intrinsic definition. We choose in
this paper to define the set ${\cal G}$ via the trace operator. Namely, we
define the model class ${\cal G}$ via the trace of functions from the unit
ball of a Besov class
$B_{\overline{q}}^{\overline{s}}(L_{\overline{p}}(\Omega))$ with
$\overline{s}>d/\overline{p}$, to ensure continuity. We also place the
restrictions $0<\bar{p},\bar{q}\leq 2$ in order to simplify the presentation
that follows. In other words, our model class assumptions on the boundary
function $g$ take the form

$g\in{\cal G}:=Tr(U(\overline{\cal B})),\ \ \hbox{where}\ \ \overline{\cal
B}:=B_{\overline{q}}^{\overline{s}}(L_{\overline{p}}(\Omega)),\quad\hbox{with}\
\ \bar{s}>d/\bar{p},\quad 0<\bar{p},\bar{q}\leq 2.$ (3.2)

Then ${\cal G}$ is a compact subset of $H^{1/2}(\partial\Omega)$. Note that in
general the parameters $s,p,q$ used in defining ${\cal F}$ are different from
$\overline{s},\overline{p},\overline{q}$ used in defining ${\cal G}$. These
assumptions imply that the function $u$ we want to numerically approximate is
an element of the set

${\cal U}:=\\{\widetilde{u}\in C(\overline{\Omega}):\ \widetilde{u}\ {\rm
satisfies}\ {\rm(\ref{main})}\ {\rm with}\ f\in{\cal F},\ g\in{\cal G}\\}.$
(3.3)

Given data ${\bf f}:=(f_{1},\ldots,f_{\widetilde{m}})$, the totality of
information we have about $f$ is that it belongs to the set

${\cal F}_{\rm data}:={\cal F}_{{\rm data}(f)}:=\\{\widetilde{f}\in{\cal F}:\
\widetilde{f}_{i}=f_{i},\ i=1,\dots,\widetilde{m}\\}.$

Similarly, if the data ${\bf g}:=(g_{1},\ldots,g_{\overline{m}})$ is coming
from a $g\in{\cal G}$, then, what we know about $g$ is that in lies in the set

${\cal G}_{\rm data}:={\cal G}_{{\rm data}(g)}:=\\{\widetilde{g}\in{\cal G}:\
\widetilde{g}_{i}=g_{i},\ i=1,\dots,\overline{m}\\}.$ (3.4)

Notice that the $g_{i}=u({\bf z}_{i})$, $i=1,\dots,\overline{m}$. Finally, the
totality of information we have about the sought after $u$ is that it is in
the set

${\cal U}_{\rm data}:={{\cal U}_{{\rm data}(u)}:=}\\{\widetilde{u}\in{\cal
U}{:\,-\Delta\widetilde{u}({\bf x}_{i})=f_{i},\
i=1,\dots,\widetilde{m},\quad\widetilde{u}({\bf z}_{i})=g_{i},\ \
i=1,\dots,\overline{m}}\\}.$

or equivalently

${\cal U}_{\rm data}:={\cal U}_{{\rm
data}(u)}:=\\{\widetilde{u}:\widetilde{u}\ \hbox{is a solution to
{\rm(\ref{main})} with }\,f\in{\cal F}_{\rm data},\quad g\in{\cal G}_{\rm
data}\\}.$

We are interested in knowing to what extent the information that $u\in{\cal
U}_{\rm data}$ identifies $u$.

There is a simple theoretical answer to such questions. Given a compact set
$K$ of a Banach space $X$, we let $B(K):=B(K)_{X}$ denote a ball with the
smallest radius in $X$ which contains $K$, called a Chebyshev ball of $K$ in
$X$. If we wish to provide an element from $X$ that simultaneously
approximates all elements in $K$ (in the norm of $X$), then the center of this
ball is the best we can do and its radius

$R(K)_{X}:={\rm rad}(B(K))_{X}$ (3.5)

is the optimal error we can obtain. It is called the error of optimal recovery
of $K$.

If ${\cal U}_{\rm data}$ is a compact subset of $X$, let $B=B({\cal U}_{\rm
data})_{X}$ be a Chebyshev ball of this set. Since all we know about $u$ is
that it is in ${\cal U}_{\rm data}$, any element from $B({\cal U}_{\rm
data})_{X}$ will give a near best approximation to $u$ in $\|\cdot\|_{X}$. The
radius of $B$

$R({\cal U}_{{\rm data}(u)})_{X}:={\rm rad}(B({\cal U}_{{\rm data}(u)}))_{X}$
(3.6)

is then the optimal recovery error. We are interested in bounds for $R({\cal
U}_{\rm data})_{X}$, where $\|\cdot\|_{X}$ is the norm in which we wish to
measure accuracy. We introduce similar notation for the recovery of $f$ and
$g$.

The recovery rate (3.6) will not only depend on the data sites ${\cal X}$, and
${\cal Z}$, but also on the values $({\bf f},{\bf g})$ assigned at these data
sites. In order to obtain uniform estimates, we fix the data sites $({\cal
X},{\cal Z})$ and introduce

$R^{*}({\cal U},{\cal X},{\cal Z})_{X}:=\sup_{u\in{\cal U}}R({\cal U}_{{\rm
data}(u)})_{X},$ (3.7)

where the data values come from any $u\in{\cal U}$. This is measuring the
worst possible performance over the class ${\cal U}$ and is called the uniform
optimal recovery rate at the fixed data sites $({\cal X},{\cal Z})$. If we
prescribe a budget for the number $m=|{\cal X}|+|{\cal Z}|$ of data sites, we
can ask for the optimal error we can achieve under such a budget restriction.
Accordingly, we define

$R_{m}^{*}({\cal U})_{X}:=\inf_{{\cal X}\subset\overline{\Omega},\,{\cal
Z}\subset\partial\Omega,\,|{\cal X}|+|{\cal Z}|=m}R^{*}({\cal U},{\cal
X},{\cal Z})_{X},\quad m\geq 2.$

Similarly,

$R^{*}({\cal F},{\cal X})_{X}:=\sup_{f\in{\cal F}}R({\cal F}_{{\rm
data}(f)})_{X}\quad\hbox{and}\quad R_{\widetilde{m}}^{*}({\cal
F})_{X}:=\inf_{{\cal X}\subset\overline{\Omega},|{\cal
X}|=\widetilde{m}}R^{*}({\cal F},{\cal X})_{X},\quad\widetilde{m}\geq 1,$

where the latter is the minimum error we can achieve over all the possible
choices of ${\cal X}=\\{{\bf x}_{1},\dots,{\bf x}_{\widetilde{m}}\\}$. We can
similarly ask the question of how well can we can recover $g$ from the data
${\bf g}$, and therefore how well we can recover $u$? Each of these questions
has an associated uniform optimal recovery error which is dependent on which
norm $\|\cdot\|_{X}$ we use to measure the error of recovery.

###  3.1 Optimal recovery of $f$

Let ${\cal F}:=U({\cal B})$ where ${\cal B}=B_{q}^{s}(L_{p}(\Omega))$ with
$0\leq q,p\leq\infty$ and $s>\frac{d}{p}$. We define the tensor product grid
$G_{k,r}$ of $\overline{\Omega}=[0,1]^{d}$ for $r\geq 2$,

$G_{k,r}:=\\{{\bf x}_{1},\ldots,{\bf
x}_{\widetilde{m}}\\}\subset\overline{\Omega},\quad\widetilde{m}=(r2^{k})^{d},$
(3.8)

with spacing $h=2^{-k}(r-1)^{-1}$ (the role of $r$ will be made clear below)
and prove that, see Theorem 3.1,

$R^{*}({\cal F},G_{k,r})_{X}\lesssim\widetilde{m}^{-\alpha_{X}}.$

Moreover, we show that for any choice of data cites ${\cal
X}\subset\overline{\Omega}$ with $|{\cal X}|=\widetilde{m}$, we have

$R^{*}({\cal F},{\cal X})_{X}\gtrsim\widetilde{m}^{-\alpha_{X}},$

and therefore $R_{\widetilde{m}}^{*}({\cal
F})_{X}\asymp\widetilde{m}^{-\alpha_{X}}$. The proof of this theorem is
presented in the Appendix.

######  Theorem 3.1.

Let $\Omega=(0,1)^{d}$ and ${\cal F}:=U({\cal B})$, where ${\cal B}$ is the
Besov space $B_{q}^{s}(L_{p}(\Omega))$ with $s>d/p$, and $0<p,q\leq\infty$.
Then the following holds for the optimal recovery rate
$R_{\widetilde{m}}^{*}({\cal F})_{X}$ in the norm of $X$: (i) If
$X=C(\Omega)$, then

$R_{\widetilde{m}}^{*}({\cal F})_{X}\asymp{\widetilde{m}}^{-\alpha_{C}},\
\widetilde{m}\geq 1,\quad{\rm where}\quad\alpha_{C}:=\frac{s}{d}-\frac{1}{p},$
(3.9)

with constants of equivalence independent of $\widetilde{m}$.

(ii) If $X=L_{\tau}(\Omega)$, $\tau>0$, then

$R_{\widetilde{m}}^{*}({\cal F})_{X}\asymp\widetilde{m}^{-\alpha_{\tau}},\
\widetilde{m}\geq 1,\quad{\rm
where}\quad\alpha_{\tau}:=\frac{s}{d}-\left[\frac{1}{p}-\frac{1}{\tau}\right]_{+},$
(3.10)

with constants of equivalence independent of $\widetilde{m}$.

(iii) If $X=H^{1}(\Omega)$, $0<p\leq 2$, then

$R_{\widetilde{m}}^{*}({\cal F})_{X}\asymp\widetilde{m}^{-\alpha_{H}},\
\widetilde{m}\geq 1,\quad{\rm
where}\quad\alpha_{H}:=\frac{s-1}{d}-\left[\frac{1}{p}-\frac{1}{2}\right],$
(3.11)

with constants of equivalence independent of $\widetilde{m}$.

(iv) If $X=H^{-1}(\Omega)$, then

  * •

if $d\geq 3$ and $0<p\leq\infty$, or $d=2$ and $1<p\leq\infty$, we have

$R_{\widetilde{m}}^{*}({\cal F})_{X}\asymp\widetilde{m}^{-\alpha_{-1}},\
\widetilde{m}\geq 1\quad{\rm
where}\quad\alpha_{-1}:=\frac{s}{d}-\left[\frac{1}{p}-\frac{1}{\delta}\right]_{+}\
{\rm and}\ \frac{1}{\delta}:=\frac{1}{2}+\frac{1}{d},$ (3.12)

with constants of equivalence independent of $\widetilde{m}$.

  * •

if $d=2$ and $0<p\leq 1$, we have

$R_{\widetilde{m}}^{*}({\cal
F})_{X}\lesssim\log(\widetilde{m})\widetilde{m}^{-\alpha_{-1}},\qquad
R_{\widetilde{m}}^{*}({\cal F})_{X}\gtrsim\widetilde{m}^{-\alpha_{-1}},\
\qquad\widetilde{m}\geq 1,$ (3.13)

with $\alpha_{-1}$ and $\delta$ as in (3.12).

Moreover, when $\widetilde{m}=(r2^{k})^{d}$, the upper bounds in each of these
estimates are obtained when we take the data sites ${\cal X}=G_{k,r}$,
provided $r>\max(s,1)$. In this case, the approximant to $f$ is given by the
continuous piecewise polynomial

$S_{k}^{*}(f)=\sum_{T\in{\cal T}_{k}}L_{T}(f)\chi_{T},$

described in (2.16), where ${\cal T}_{k}={\cal T}_{k}(\Omega)=\bigcup T$ is
the collection of all simplicies arising from the (Kuhn-Tucker) conforming
decomposition of the dyadic cubes in ${\cal D}_{k}(\Omega)$, the $\chi_{T}$’s
are the characteristic functions of $T$, and the $L_{T}(f)$’s are the
polynomials of order $r$ ($L_{T}(f)\in{\cal P}^{d}_{r}$), gotten by
interpolating $f$ at the data points $T\cap G_{k,r}$.

We remark that in the case $d=2$ and $0<p\leq 1$, the logarithm appears
because of the failure of the Sobolev embedding to provide
$L_{1}(\Omega)\subset H^{-1}(\Omega)$. However, it is actually true that the
Hardy space ${\cal H}^{1}(\Omega)$ embeds into $H^{-1}(\Omega)$. This could
potentially be used to improve our later analysis, but for the sake of
simplicity we do not pursue this here.

We make some further observations that will help explain the above theorem,
especially what is happening for recovery in $H^{-1}(\Omega)$. Notice that
under the model class assumption $f\in{\cal F}=U(B_{q}^{s}(L_{p}(\Omega)))$
with $s>d/p$, the OR for this class is the same whether we measure error in
$H^{-1}(\Omega)$ or in $L_{\delta}(\Omega)$ when $d\geq 3$ or $d=2$,
$1<p\leq\infty$. We also have that

$\|f\|_{H^{-1}(\Omega)}\leq\begin{cases}C(d)\|f\|_{L_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\delta}(\Omega)},\quad
d\geq 3,\\\ \\\
\frac{C}{\tau-1}\|f\|_{L_{\tau}(\Omega)},\quad\,\,d=2,\quad\tau>1,\end{cases}$
(3.14)

with $C$ an absolute constant, as discussed in Lemma 12.4 of the appendix .
These two facts explain the form of the loss function ${\cal L}^{*}$ which we
discuss in detail later.

###  3.2 Optimal recovery of $g$

In this section, we study the optimal recovery rate of functions $g\in{\cal
G}$ in the norm of $H^{1/2}(\partial\Omega)$. Let $Tr=Tr_{\partial\Omega}$ be
the trace operator onto $\partial\Omega$ and let $\overline{\cal
B}:=B_{\bar{q}}^{\bar{s}}(L_{\bar{p}}(\Omega))$, where we fix
$\bar{s},\bar{p},\bar{q}$ and assume that $\overline{s}>d/\overline{p}$ so
that we are sure that the functions in $\overline{\cal B}$ are continuous on
$\overline{\Omega}$. We consider the case $0<\bar{p}\leq 2$, which together
with the restriction $d\geq 2$ implies that $\bar{s}-d(1/\bar{p}-1/2)>1$, and
thus the unit ball of $\overline{\cal B}$ compactly embeds into
$H^{1}(\Omega)$. We define the model class

${\cal G}:=\\{g:\,g=Tr(v):\ \|v\|_{\overline{\cal B}}\leq 1\\},$ (3.15)

and use the trace norm definition of $\|\cdot\|_{H^{1/2}(\partial\Omega)}$
throughout this section.

Recall the tensor product grid $G_{k,r}$ for $\overline{\Omega}=[0,1]^{d}$,
see (3.8). We let

$\overline{G}_{k,r}:=\\{{\bf z}_{i}:\
i=1,\dots,\overline{m}\\}=G_{k,r}\cap\partial\Omega$ (3.16)

denote the set of those grid points of $G_{k,r}$ that are on the boundary
$\partial\Omega$. We will use the points ${\cal Z}=\overline{G}_{k,r}$ to
recover $g$. Note that the number $\bar{m}$ of data sites in
$\overline{G}_{k,r}$ is

$\overline{m}=\\#(\overline{G}_{k,r})=2d[r2^{k}]^{d-1}\asymp 2^{k(d-1)},$
(3.17)

with constants of equivalency depending only on $r$ and $d$.

Given $g\in{\cal G}$, we define the continuous piecewise polynomial

$\overline{S}_{k}(g):=\sum_{i=1}^{\overline{m}}g({\bf
z}_{i})\overline{\phi}_{i},$ (3.18)

where each $\overline{\phi}_{i}$ is the trace of the Lagrange element
$\phi_{i}$ centered at ${\bf z}_{i}\in\overline{G}_{k,r}$. Note that if
$v\in\overline{\cal B}$ is any function whose trace on $\partial\Omega$ is
$g$, and $S_{k}^{*}(v)$ is its continuous piecewise polynomial Lagrange
interpolant at the grid points $G_{k,r}$, see (2.16), then we have

$Tr(S_{k}^{*}(v))=\overline{S}_{k}(g).$

This follows from the fact that each Lagrange interpolant $\phi_{i}$ centered
at ${\bf x}_{i}\in G_{k,r}\setminus\overline{G}_{k,r}$ vanishes at the faces
of the simplex from the corresponding Kuhn-Tucker decomposition that contains
${\bf x}_{i}$. We also want to point out that $\overline{S}_{k}(g)$ does not
depend on $v$.

The proof of the following result is provided in the appendix.

######  Theorem 3.2.

Let $\overline{\cal B}=B_{\overline{q}}^{\overline{s}}(L_{\bar{p}}(\Omega))$
with $\bar{s}>d/\bar{p}$, $0<\bar{p}\leq 2$, $0<q\leq\infty$, and let ${\cal
G}=\\{Tr(v):\|v\|_{\overline{\cal B}}\leq 1\\}$. Then the optimal recovery
rates of the model classe ${\cal G}$ in $H^{1/2}(\partial\Omega)$ is

$R_{\overline{m}}^{*}({\cal
G})_{H^{1/2}(\partial\Omega)}\asymp\overline{m}^{-\beta},\quad\overline{m}\geq
1,\quad\hbox{where}\quad\beta:=\frac{\overline{s}-1}{d-1}-\frac{d}{d-1}\left[\frac{1}{\bar{p}}-\frac{1}{2}\right],$
(3.19)

with constants of equivalency independent of $\overline{m}$.

Moreover, when $\overline{m}=\\#(\overline{G}_{k,r})$, the upper bound in
(3.19) is achieved when we take the data sites ${\cal Z}=\overline{G}_{k,r}$,
provided $r>\max(s,1)$ and the approximant to $g$ is given by the function
$\overline{S}_{k}(g)$.

To understand the exponent $\beta$ in (3.19) better, let $v$ be a function in
$U(B_{q}^{\bar{s}}(L_{\bar{p}}(\Omega)))$ whose trace is $g$. Note that
$\beta$ can be rewritten as

$\beta=\frac{\overline{s}-1-(\frac{d}{\bar{p}}-\frac{d}{2})}{d-1},$

where the numerator is the excess regularity of $v$ in $H^{1}(\Omega)$ and
thereby corresponds to the excess regularity of $g$ in
$H^{1/2}(\partial\Omega)$, and the denominator is $d-1$ since $\partial\Omega$
has dimension $d-1$.

###  3.3 Optimal recovery of $u$

We turn now to determining the optimal $H^{1}(\Omega)$ recovery rate for the
functions $u$ in the model class ${\cal U}$ from data

$-\Delta u({\bf x}_{i})=f_{i},\quad{\bf x}_{i}\in{\cal X},\
i=1,\dots,\widetilde{m};\quad u({\bf z}_{i})=g_{i},\quad{\bf z}_{i}\in{\cal
Z},\ i=1,\dots\overline{m},$ (3.20)

under a given budget $m=\widetilde{m}+\overline{m}$ of data observations.
Recall that ${\cal F}:=U(B_{q}^{s}(L_{p}(\Omega)))$, where $s>d/p$ and
$0<q,p\leq\infty$, ${\cal
G}:=Tr(U(B_{\overline{q}}^{\overline{s}}(L_{\overline{p}}(\Omega)))$, where
$\overline{s}>d/\overline{p}$ and $0<\overline{p}\leq 2$, and

${\cal U}:=\\{u\in C(\overline{\Omega})\ :\ u\ {\rm satisfies}\
{\rm(\ref{main})}\ {\rm with}\ f\in{\cal F},g\in{\cal G}\\}.$

The following theorem gives the optimal recovery rate for ${\cal U}$.

######  Theorem 3.3.

Let $\Omega=(0,1)^{d}$ and let ${\cal F}$, ${\cal G}$, ${\cal U}$ be as above.
Then the following holds for the optimal recovery rate $R^{*}_{m}({\cal
U})_{H^{1}(\Omega)}$:

  * •

if $d\geq 3$, or $d=2$ and $p>1$, we have

$R^{*}_{m}({\cal U})_{H^{1}(\Omega)}\asymp m^{-\min(\alpha,\beta)},\quad m\geq
2,$ (3.21)

  * •

if $d=2$ and $0<p\leq 1$, we have

$m^{-\min(\alpha,\beta)}\lesssim R^{*}_{m}({\cal
U})_{H^{1}(\Omega)}\lesssim\log(m)m^{-\alpha}+m^{-\beta},\quad\ m\geq 2,$
(3.22)

where $\alpha:=\alpha_{-1}$ is given by (3.12), $\beta$ is given by (3.19),
and all constants of equivalence are independent of $m$.

###### Proof.

Let us fix the data cites ${\cal X},{\cal Z}$ with $|{\cal X}|=\widetilde{m}$,
$|{\cal Z}|=\overline{m}$, $m=\overline{m}+\widetilde{m}$, and let
$u_{1},u_{2}\in{\cal U}_{\rm data}$ with corresponding $f_{1},f_{2}\in{\cal
F}_{\rm data}$ and $g_{1},g_{2}\in{\cal G}_{\rm data}$. It follows from (1.3)
that

$\|u_{1}-u_{2}\|_{H^{1}(\Omega)}\asymp\|f_{1}-f_{2}\|_{H^{-1}(\Omega)}+\|g_{1}-g_{2}\|_{H^{1/2}(\Omega)}.$

It follows that

$R^{*}({\cal U},{\cal X},{\cal Z})_{H^{1}(\Omega)}\asymp R^{*}({\cal F},{\cal
X})_{H^{-1}(\Omega)}+R^{*}({\cal G},{\cal Z})_{H^{1/2}(\Omega)}.$ (3.23)

Clearly, if $\overline{m}=\widetilde{m}\asymp m/2$ and using the rates from
Theorem 3.1 and Theorem 3.2, we have that

$R^{*}_{m}({\cal
U})_{H^{1}(\Omega)}\lesssim\begin{cases}m^{-\alpha}+m^{-\beta}\lesssim
m^{-\min(\alpha,\beta)},\quad d\geq 3,\,\hbox{or}\,\,d=2,\,p>1,\\\ \\\
\log(m)m^{-\alpha}+m^{-\beta},\,\,\quad\quad\quad\quad d=2,\,0<p\leq
1.\end{cases}$

On the other hand, using the same theorems and (3.23), we obtain

$\displaystyle
m^{-\min(\alpha,\beta)}\lesssim\inf_{m=\overline{m}+\widetilde{m}}(\widetilde{m}^{-\alpha}+\overline{m}^{-\beta})$
$\displaystyle\lesssim$
$\displaystyle\inf_{m=\overline{m}+\widetilde{m}}\left(R_{\widetilde{m}}^{*}({\cal
F})_{H^{-1}(\Omega)}+R_{\overline{m}}^{*}({\cal G})_{H^{1/2}(\Omega)}\right)$
$\displaystyle\lesssim$ $\displaystyle\inf_{m=|{\cal X}|+|{\cal
Z}|}R^{*}({\cal F},{\cal X})_{H^{-1}(\Omega)}+R^{*}({\cal G},{\cal
Z})_{H^{1/2}(\Omega)}$ $\displaystyle\lesssim$ $\displaystyle\inf_{m=|{\cal
X}|+|{\cal Z}|}R^{*}({\cal U},{\cal X},{\cal
Z})_{H^{1}(\Omega)}=R^{*}_{m}({\cal U})_{H^{1}(\Omega)},$

and the proof is completed. ∎

###  3.4 Final observations on optimal recovery

This section gave the optimal recovery rates for the model classes ${\cal
F},{\cal G},{\cal U}$ from $\widetilde{m},\overline{m},m$, data observations,
respectively. While we began by considering general model class assumptions
${\cal F}$ and ${\cal G}$ for $f$ and $g$ to be the unit balls of Besov spaces
$B_{q}^{s}(L_{p})$, we found that the optimal recovery rate is the same for
many of these model classes. This means that some of the model classes are
superflous in that they are contained in a larger model class with the same
recovery rate. It makes sense only to use these largest model classes in going
forward.

Consider, for example, the model classes ${\cal F}=U({\cal B})$ with ${\cal
B}=B_{q}^{s}(L_{p}(\Omega))$ for $f$. We found that in the case $d\geq 3$ and
$0<p\leq\infty$, or $d=2$ and $1<p\leq\infty$, the optimal recovery rate of
${\cal F}$ in $H^{-1}(\Omega)$ is the same as its optimal recovery rate in
$L_{\delta}(\Omega)$ with ${\delta}=\frac{2d}{d+2}$. Moreover, the model class
${\cal F}$ has an optimal recovery rate $m^{-\alpha/d}$ in $H^{-1}(\Omega)$ if
and only if it is contained in $U(B_{\infty}^{\alpha}(L_{\gamma}(\Omega)))$
and, in addition, ${\cal F}$ embeds into $C(\Omega)$. Largest Model Classes
for $f$:

In the case $d\geq 3$, $0<p\leq\infty$ or $d=2$, $1<p\leq\infty$ the class

${\cal F}:=U({\cal B}),\quad{\cal B}=B_{\infty}^{s}(L_{p}(\Omega)),\ {\rm
with}\ p\geq{\delta}\ {\rm and}\ s>\frac{d}{p}$ (3.24)

has the optimal recovery rate

$R_{m}^{*}({\cal F})_{H^{-1}(\Omega)}\asymp{m}^{-s/d},\quad m\geq 1,$ (3.25)

and any Besov model class that gives the optimal recovery rate $O(m^{-s/d})$
is contained in one of these model classes. In each of these cases, the number
$s$ represents the smoothness of $f\in{\cal F}$ in
$L_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\delta}(\Omega)$.
These largest model classes are pictured in Figure 3.1.

In the case $d=2$, $0<p\leq 1$, the largest class

${\cal F}:=U({\cal B}),\quad{\cal B}=B_{\infty}^{s}(L_{1}(\Omega)),\ {\rm
with}\ s>2,$ (3.26)

has recovery rate

${m}^{-s/2}\lesssim R_{m}^{*}({\cal
F})_{H^{-1}(\Omega)}\lesssim(1+\log(m)){m}^{-s/2},\quad m\geq 1.$ (3.27)

Largest Model Classes for $g$: We have considered the model classes ${\cal
G}=Tr(B^{\bar{s}}_{\infty}(L_{\bar{p}}(\Omega))$ with $0<\bar{p}\leq 2$ and
$\bar{s}>\max\\{\frac{d}{\bar{p}},1\\}$. Thus, all of the model classes which
give a given rate $O(\overline{m}^{-\alpha})$ will all be contained in one
model class

${\cal G}:=Tr(U(\overline{\cal B})),\quad\overline{\cal
B}=B_{\infty}^{\bar{s}}(L_{2}(\Omega)),\quad\ {\rm with}\ \bar{s}>d/2.$ (3.28)

This class has optimal recovery rate in $H^{1/2}(\partial\Omega)$

$R_{m}^{*}({\cal G})_{H^{1/2}(\partial\Omega)}\asymp
m^{-(\bar{s}-1)/(d-1)},\quad m\geq 1.$ (3.29)

Moreover, any Besov model class that gives this recovery rate is contained in
one of these largest model classes. For each of these largest model classes
$\bar{s}-1$ is the excess smoothness of $v$ in $H^{1}(\Omega)$ and also the
excess smoothness of $g=Tr(v)$ in $H^{1/2}(\Omega)$.

######  Remark 3.4.

Because of the above remarks, in going further, we always take ${\cal F}$ as
in (3.24) and ${\cal G}$ as in (3.28).

s$1/\delta$$1/p$$0$$s{=}d/p$s$1/2$$1/p$$0$$\frac{d}{2}$$s{=}d/p$ Figure 3.1:
The spaces $B^{s}_{\infty}(L_{p}(\Omega))$ are represented by the point
$(1/p,s)$ in the first quadrant. The arrows indicate embedding with a slope of
$d$ for the oblique lines. In particular, the spaces above the line $s=d/p$
embed into $C(\overline{\Omega})$. (Left) The largest model classes for $f$
are the unit balls of spaces just above the thick and dashed segment or on the
on the thick vertical half-line (excluding the point
$(1/{\delta},d/{\delta})$). (Right) The largest model classes for $g$ are the
traces of the unit balls of spaces on the thick vertical half-line (excluding
the point $(1/2,d/2)$).

Model Classes for $u$: If we use one of the model classes from (3.24) for
${\cal F}$ and one of the model classes from (3.28) for ${\cal G}$, we obtain
a model class ${\cal U}$ for $u$. The optimal recovery rate for this model
class is

$R_{m}^{*}({\cal U})_{H^{1}(\Omega)}\asymp
m^{-\min\\{s/d,(\bar{s}-1)/(d-1)\\}},\quad m\geq 2,$

which could be obtained by assigning $\widetilde{m}\asymp m/2$ points in
$\overline{\Omega}$ and $\overline{m}\asymp m/2$ points on the boundary
$\partial\Omega$.

##  4 Numerical algorithms based on linear approximation

The above analysis of optimal recovery tells us to what extent the information
$({\bf f},{\bf g})$ determines $u$. That is, the optimal recovery analysis
tells us that if two functions $u,\widetilde{u}$ from our model class ${\cal
U}$ satisfy the same data, then they are close in the $H^{1}(\Omega)$ norm,
that is

$\|u-\widetilde{u}\|_{H^{1}(\Omega)}\lesssim R_{m}^{*}({\cal
U})_{H^{1}(\Omega)}.$

However, it does not give a numerical algorithm that takes the data and
creates an approximation $\widehat{u}$ to $u$ with accuracy of the optimal
recovery rate. In this section, we discuss two numerical methods which would
accomplish the latter task. The methods discussed in this section are not in
the form of PINNs. The latter will be discussed in §9.

We assume in this section that the data sites form a tensor product grid $G$
of $\overline{\Omega}$. We know such data sites $G$ provide a near optimal
recovery. We have shown that if we use the data ${\bf f}$ at $G$, then we can
numerically construct a piecewise polynomial interpolant $\widehat{f}$ to $f$
at the data sites which is a near optimal (in terms of the budget
$\widetilde{m}$) approximation to $f$ in the $H^{-1}(\Omega)$ norm. Similarly,
we can numerically construct a $\widehat{g}$ from the data ${\bf g}$ which is
a near optimal recovery of $g$ in the $H^{1/2}(\partial\Omega)$ norm.

We define $\widetilde{u}$ as the solution to (1) with right hand side
$\widehat{f}$ and boundary value $\widehat{g}$. We know that $\widetilde{u}$
is a near optimal recovery of $u$ in the $H^{1}(\Omega)$ norm. Then, in view
of (1.3), we have

$\|u-\widetilde{u}\|_{H^{1}(\Omega)}\lesssim R_{m}^{*}({\cal
U})_{H^{1}(\Omega)}.$ (4.1)

Note that the constructions of $\widehat{f}$ and $\widehat{g}$ are numerical
and we can control their complexity. In order to do that for $\widetilde{u}$,
we use an existing numerical method for solving (1). This will, of course,
incur a numerical error depending on the algorithm we choose. We mention two
natural possibilities.

###  4.1 Using FEMs or AFEMs

Let $\widehat{u}_{n}$ be the approximation to $\widetilde{u}$ obtained from
the Galerkin projection onto a standard finite element space $V_{n}$ of
dimension $n$. That is, $\widehat{u}_{n}$ is the Galerkin solution to (1) with
right side $\widehat{f}$ and boundary value $\widehat{g}$. Then, we have

$\|u-\widehat{u}_{n}\|_{H^{1}(\Omega)}\lesssim R_{m}^{*}({\cal
U})_{H^{1}(\Omega)}+\varepsilon_{n}({\cal U})_{H^{1}(\Omega)},$ (4.2)

where

$\varepsilon_{n}({\cal U})_{H^{1}(\Omega)}:=\mathop{\rm dist}({\cal
U},V_{n})_{H^{1}(\Omega)},\quad n\geq 1.$ (4.3)

The term $\varepsilon_{n}$ will tend to zero as $n\to\infty$ and we would
obtain the optimal recovery error on the right side of (4.2) if $n$ is chosen
suitably large. The actual decay rate of $\varepsilon_{n}$ depends on the
regularity of the model class ${\cal U}$ in $H^{1}(\Omega)$. Therefore, the
issue becomes what do our model class assumptions $f\in{\cal F}$, $g\in{\cal
G}$, say about the regularity of $u$ in the $H^{1}(\Omega)$ norm. If we use
linear spaces such as standard finite element spaces, then we would want to
know the regularity of $u$ in the scale of $H^{t}(\Omega)$ spaces, $t>1$.
There are several theorems that obtain such regularity results (see e.g. [11,
35, 46] and the references therein) and thereby obtain concrete bounds on the
decay rate $\varepsilon_{n}$. We do not elaborate on this further but refer
the interested reader to the papers cited above.

If we use nonlinear numerical methods such as AFEMs from for solving (1) with
right side $\widehat{f}$, and boundary value $\widehat{g}$, then the decay
rate for $\varepsilon_{n}({\cal U})$ can be improved (see [5] for a summary of
such results).

###  4.2 Reduced models

If one is faced with solving several PDEs with the same data sites $({\cal
X},{\cal Z})$ but with different data $({\bf f},{\bf g})$, one can build a
numerical algorithm as follows. Suppose that ${\cal X}=G_{k,r}$,
$|G_{k,r}|=\widetilde{m}=[r2^{k}]^{d}$, see (3.8) and ${\cal
Z}=\overline{G}_{k,r}$, $|\overline{G}_{k,r}|=\overline{m}=2d[r2^{k}]^{d-1}$,
see (3.16) $\overline{\Omega}$ with spacing $h=r^{-1}2^{-k}$.

We know that the simplicial interpolation operator $S_{k}^{*}$ defined in
(2.16) provides a near optimal recovery for ${\cal F}$. Namely,

$\|f-S_{k}^{*}(f)\|_{H^{-1}(\Omega)}\lesssim R_{{\widetilde{m}}}^{*}({\cal
F})_{H^{-1}(\Omega)},\quad\widetilde{m}=[r2^{k}]^{d}.$ (4.4)

We write

$S_{k}^{*}(f)=\sum_{i=1}^{\widetilde{m}}f({\bf x}_{i})\phi_{i},$ (4.5)

where the $\phi_{i}$ are the local Lagrange functions with $\phi_{i}({\bf
x}_{i})=1$ and $\phi_{i}({\bf x}_{j})=0$ when $j\neq i$, supported on the
simplex containing the point ${\bf x}_{i}$. Let $\phi_{i}^{*}$ be the solution
to (1) with right side $\phi_{i}$ and zero boundary conditions.

Similarly, the interpolant $\overline{S}_{k}(g)$ to $g\in{\cal G}$ used for
optimal recovery, see (3.18), can be written as

$\overline{S}_{k}(g)=\sum_{i=1}^{\overline{m}}g({\bf
z}_{i})\psi_{i},\quad\overline{m}=2d[r2^{k}]^{d-1}.$ (4.6)

Let $\psi_{i}^{*}$ be the solution to (1) with $g=\psi_{i}$ and zero right
side.

Given any data $({\bf f},{\bf g})$, the function

$\widetilde{u}:=\sum_{i=1}^{\widetilde{m}}f({\bf
x}_{i})\phi_{i}^{*}+\sum_{i=1}^{\overline{m}}g({\bf z}_{i})\psi_{i}^{*}$ (4.7)

is a near optimal recovery of the solution $u$ to (1) with right side $f$ and
boundary values $g$. Note that $\widetilde{u}$ is the solution to (1) with
right side $S_{k}^{*}(f)$ and boundary value $\overline{S}_{k}(g)$.

Given a budget $m$, one computes offline to a sufficiently high accuracy
suitable approximations $\widehat{\phi}_{i}$ and $\widehat{\psi}_{i}$ in
$H^{1}(\Omega)$ to the $\phi_{i}^{*}$ and $\psi_{i}^{*}$, respectively. Then
the function

$\widehat{u}:=\sum_{i=1}^{\widetilde{m}}f({\bf
x}_{i})\widehat{\phi}_{i}+\sum_{i=1}^{\overline{m}}g({\bf
z}_{i})\widehat{\psi}_{i}$ (4.8)

is a near optimal recovery for the solution $u$ to (1) with $f\in{\cal F}$,
$g\in{\cal G}$, where $\left.f\right|_{{\cal X}}={\bf f}$ and
$\left.g\right|_{{\cal Z}}={\bf g}$.

Note also that the near optimality can be guaranteed for all model class
assumptions provided $s$ and $\overline{s}$ are all less than a fixed number
$s_{0}$.

##  5 Numerical methods based on optimization

Another approach to finding a near optimal approximation to $u$ from the given
data is through optimization, such as the one used in PINNS. To analyze this
approach, we begin this section by first considering a theoretical
optimization algorithm which we shall later turn into a numerical procedure.
We begin by introducing the functional ${\cal
L}_{T}:H^{1}(\Omega)\to\mathbb{R}$, defined as

${\cal L}_{T}(v):=\|f+\Delta
v\|_{H^{-1}(\Omega)}+\|v-g\|_{H^{1/2}(\partial\Omega)},$ (5.1)

for fixed $f\in{\cal F}\subset H^{-1}(\Omega)$ and $g\in{\cal G}\subset
H^{1/2}(\partial\Omega)$. Clearly, the unique solution $u\in H^{1}(\Omega)$ to
(1) with right side $f$ and boundary condition $g$ is the unique minimizer of
this functional since ${\cal L}_{T}(u)=0$. In other words, we have

$u=\mathop{\rm argmin}_{v\in H^{1}(\Omega)}{\cal L}_{T}(v).$ (5.2)

Note also that because of (1.3), we have for every $v\in H^{1}(\Omega)$

$\|u-v\|_{H^{1}(\Omega)}\asymp\|f+\Delta
v\|_{H^{-1}(\Omega)}+\|g-v\|_{H^{1/2}(\partial\Omega)}={\cal L}_{T}(v),$ (5.3)

with absolute constants in the equivalence. In other words, the loss ${\cal
L}_{T}(v)$ for any $v$ tells us how close $v$ is to the true solution $u$.

The minimization in (5.2) is over all $v\in H^{1}(\Omega)$. To make this
minimization more amenable to numerical implementation, we take this
minimization over a set $\Sigma=\Sigma_{n}$, where $\Sigma$ is either a linear
space of dimension $n$ or a nonlinear manifold of order $n$ (i.e., depending
on $n$ parameters), and study any function

$S_{\Sigma}\in\\{\mathop{\rm argmin}_{S\in\Sigma}{\cal L}_{T}(S)\\}.$ (5.4)

We want to understand how close such an $S_{\Sigma}$ is to $u$ with error
measured in the norm of $H^{1}(\Omega)$. Therefore, we need to require that
$\Sigma\subset H^{1}(\Omega)$.

Let us fix our model class ${\cal U}$ as given in (3.3). If we want
$S_{\Sigma}$ to be a good approximation to $u$, then we need $\Sigma$ to be
good at approximating the elements of ${\cal U}$. We define the error

$E({\cal U},\Sigma):=\sup_{v\in{\cal
U}}\inf_{S\in\Sigma}\|v-S\|_{H^{1}(\Omega)}.$

If we return to our specific $u\in{\cal U}$ determined by $f$ and $g$, then
for any $S_{\Sigma}$ from (5.4) and any $S\in\Sigma$, we have from (5.3) that

$\|u-S_{\Sigma}\|_{H^{1}(\Omega)}\lesssim{\cal L}_{T}(S_{\Sigma})\lesssim{\cal
L}_{T}(S)\lesssim\|u-S\|_{H^{1}(\Omega)},$ (5.5)

where the last inequality is a consequence of the left inequality in (1.3).
Since (5.5) holds for all $S\in\Sigma$, we have that

$\|u-S_{\Sigma}\|_{H^{1}(\Omega)}\lesssim\inf_{S\in\Sigma}\|u-S\|_{H^{1}(\Omega)}\lesssim
E({\cal U},\Sigma).$ (5.6)

The discussion that we have just given shows that if we solve the continuous
minimization problem (5.2) over $\Sigma$ instead of all of $H^{1}(\Omega)$,
then the solution will get as close to $u$ as the efficiency of $\Sigma$ in
approximating ${\cal U}$, i.e., $E({\cal U},\Sigma)$. This error can be made
as small as we wish by taking finer spaces for $\Sigma$, i.e., letting $n$
increase. The error we incur is then proportional to the error in $E({\cal
U},\Sigma)$ which tell us that the best candidates for $\Sigma$ are those that
approximate our model class ${\cal U}$ well.

Several issues arise that prevent the direct implementation of the above loss
in our setting. The first of these is that we do not know neither $f$ nor $g$.
We only know these functions through the given data. The remedy for this is to
introduce a surrogate for the $H^{-1}(\Omega)$ and $H^{1/2}(\partial\Omega)$
norms that use only the data information $({\bf f},{\bf g})$ we have about $f$
and $g$. The next section addresses this issue.

##  6 Discretizing norms

The optimization procedure of the previous section is not a numerical
algorithm since it does not incorporate numerical methods for estimating the
norms appearing in the loss ${\cal L}_{T}$, namely
$\|\cdot\|_{H^{-1}(\Omega)}$ and $\|\cdot\|_{H^{1/2}(\partial\Omega)}$. In
this section, we address this issue by introducing discrete norms which
involve only the values of functions at the data sites. Replacing the norms in
${\cal L}$ by these discrete norms leads to a loss which can be numerically
computed. The equivalence between the norms in ${\cal L}_{T}$ and their
discrete counterparts holds (modulo corresponding optimal recovery rates) for
functions whose Laplacians are uniformly bounded in $\cal B$ and whose traces
are uniformly bounded in $\overline{\cal B}$.

###  6.1 A discrete $L_{\tau}$ norm

In this section, we study the discretization of the $L_{\tau}$ norm for any
$1\leq\tau\leq\infty$. Let ${\cal F}=U(B_{q}^{s}(L_{p}(\Omega)))$, $s>d/p,$ be
a model class assumption on $f$. We consider the uniform grid
$G_{k,r}\subset\overline{\Omega}$, see (3.8), consisting of $m$ points,
$m=|G_{k,r}|=[r2^{k}]^{d}$, $k\geq 0$ and $r>\max\\{s,1\\}$. For any
continuous function $f$ and any $1\leq\tau<\infty$, we define

$\|f\|^{*}_{L_{\tau}}:=\left[\frac{1}{m}\sum_{j=1}^{m}|f({\rm\bf
x}_{j})|^{\tau}\right]^{1/\tau},\quad\hbox{where}\quad{\bf x}_{j}\in G_{k,r}.$
(6.1)

When $\tau=\infty$ we make the usual modification to define the
$\|f\|^{*}_{L_{\infty}}$ norm.

The following lemma shows that the discrete $L_{\tau}$ norm is close to the
true $L_{\tau}$ norm for functions in the model class ${\cal F}$.

######  Lemma 6.1.

Let ${\cal B}=B_{q}^{s}(L_{p}(\Omega))$, $\Omega=(0,1)^{d}$, be a Besov space
with $s>d/p$. If $f\in{\cal B}$, then for any $1\leq\tau\leq\infty$, we have

$\|f\|_{L_{\tau}(\Omega)}\lesssim\|f\|^{*}_{L_{\tau}}+\|f\|_{{\cal
B}}m^{-\alpha_{\tau}},\quad\hbox{and}\quad\|f\|^{*}_{L_{\tau}}\lesssim\|f\|_{L_{\tau}(\Omega)}+\|f\|_{{\cal
B}}m^{-\alpha_{\tau}},$

with $\alpha_{\tau}:=\frac{s}{d}-\left(\frac{1}{p}-\frac{1}{\tau}\right)_{+}$.
Recall that $m^{-\alpha_{\tau}}$ is the uniform rate of optimal recovery of
the model class ${\cal F}$ in $L_{\tau}(\Omega)$. The constant in $\lesssim$
depends only on $r$ and $d$.

###### Proof.

We prove the lemma for $1\leq\tau<\infty$. A simple modification of the proof
handles the case $\tau=\infty$. We use the simplicial interpolation operator
$S^{*}_{k}$ given in (2.16). Let ${\cal T}_{k}:={\cal T}_{k}(\Omega)$ be the
set of simplices that make up the simplicial decomposition of the cubes in
${\cal D}_{k}(\overline{\Omega})$. Then, from (2.16) we have

$\|S_{k}^{*}(f)\|^{\tau}_{L_{\tau}(\Omega)}=\sum_{T\in{\cal
T}_{k}}\|L_{T}(f)\|^{\tau}_{L_{\tau}(T)},$ (6.2)

where $L_{T}(f)$ is the polynomial in ${\cal P}^{d}_{r}$ that interpolates the
data $\\{f({\bf x}_{i}):\,{\bf x}_{i}\in\overline{T}\\}$. By equivalence of
norms on ${\cal P}^{d}_{r}$, we know that

$\|L_{T}(f)\|_{L_{\tau}(T)}\asymp\Bigl{[}\frac{1}{m}\sum_{{\bf
x}_{i}\in{\overline{T}}}|f({\bf x}_{i})|^{\tau}\Bigr{]}^{1/\tau},$ (6.3)

with constants of equivalency depending only on $d,r$. Summing over $T\in{\cal
T}_{k}$ the quantities in (6.3) gives

$\|S_{k}^{*}(f)\|_{L_{\tau}(\Omega)}\asymp\Bigl{[}\frac{1}{m}\sum_{i=1}^{m}|f({\bf
x}_{i})|^{\tau}\Bigr{]}^{1/\tau}=\|f\|^{*}_{L_{\tau}(\Omega)}.$ (6.4)

For $f\in{\cal B}=B_{q}^{s}(L_{p}(\Omega)$, the piecewise polynomial
interpolant $S_{k}^{*}(f)$ provides a near optimal recovery of $f$ from the
data ${\bf f}$, that is

$\|f-S_{k}^{*}(f)\|_{L_{\tau}(\Omega)}\leq C\|f\|_{{\cal
B}}m^{-\alpha_{\tau}}.$ (6.5)

It follows then from (6.4) that

$\left|\|f\|_{L_{\tau}(\Omega)}-\|S_{k}^{*}(f)\|_{L_{\tau}(\Omega)}\right|\leq
C\|f\|_{{\cal B}}m^{-\alpha_{\tau}},$ (6.6)

and the proof is completed. ∎

###  6.2 Discrete $H^{1/2}(\partial\Omega)$ norms

In this section, we introduce discrete $H^{1/2}(\partial\Omega)$ norms and
discuss their accuracy in computing the true $H^{1/2}(\partial\Omega)$ norm
for functions $g$ in our model class ${\cal G}$. The discrete norm we
introduce will only uses the values of $g$ at the data sites and can therefore
be computed from the given data. We fix any $k\geq 1$ and let ${\cal
Z}:=\overline{G}_{k,r}=\\{{\bf
z}_{j}\\}_{j=1}^{\bar{m}}=G_{k,r}\cap\partial\Omega$, see (3.16), be the grid
points on the boundary $\partial\Omega$. The number of grid points in ${\cal
Z}$ is $\overline{m}=2d[r2^{k}]^{d-1}$.

We fix a model class ${\cal G}=Tr(U(\overline{\cal B}))$ of the form (3.28).
The optimal recovery rate for this class is

$R_{m}^{*}({\cal G})_{H^{1/2}(\partial\Omega)}\asymp m^{-\beta},\quad m\geq
1,\quad\beta=\frac{\bar{s}-1}{d-1}.$ (6.7)

This also serves to fix $r$ which we recall is the smallest integer satisfying
$r>\max(\bar{s},1)$. For any continuous function $g$ which is the trace of a
$v\in\overline{B}$, we use the trace norm

$\|g\|_{Tr(\overline{\cal B})}:=\inf_{Tr(v)=g}\|v\|_{\overline{B}}.$ (6.8)

Recall that the intrinsic semi-norm for $H^{1/2}(\partial\Omega)$ is given by

$|g|_{H^{1/2}(\partial\Omega)}^{2}:=\intop\limits_{\partial\Omega\times\partial\Omega}\frac{|g({\bf
x})-g({\bf y})|^{2}}{|{\bf x}-{\bf y}|^{d}}\,d{\bf x}d{\bf y}.$ (6.9)

We obtain the intrinsic $H^{1/2}(\partial\Omega)$ norm by adding
$\|g\|_{L_{2}(\partial\Omega)}$ to this semi-norm. It is known that this
definition is equivalent to the trace norm definition (6.8) for this domain
$\Omega$ (see [27]). We shall use (6.9) through out this section.

Let $\overline{S}_{k}$ be as defined in (3.18). We have proven in Theorem 3.2
that

$\|g-\overline{S}_{k}(g)\|_{H^{1/2}(\partial\Omega)}\leq
C\|g\|_{Tr(\overline{\cal B})}\bar{m}^{-\beta},$ (6.10)

with $C$ independent of $\overline{m}$ and $g$. This gives the comparison

$\left|\|g\|_{H^{1/2}(\partial\Omega)}-\|\overline{S}_{k}(g)\|_{H^{1/2}(\partial\Omega)}\right|\leq\|g-\overline{S}_{k}(g)\|_{H^{1/2}(\partial\Omega)}\leq
C\|g\|_{Tr(\overline{\cal B})}\bar{m}^{-\beta}.$

We concentrate therefore on finding a discrete $H^{1/2}(\partial\Omega)$ norm
for the functions $\overline{S}_{k}(g)$, $g\in{\cal G}$. Note that
$\overline{S}_{k}(g)$ is completely determined by the values of $g$ at the
data sites $\overline{G}_{k,r}$.

Given $k\geq 0$, let ${\cal T}_{k}:={\cal T}_{k}(\partial\Omega)$ denote the
collection of all $d-1$ dimensional simplicies given by the Kuhn-Tucker
partition of the dyadic cubes $I\in{\cal D}_{k}(\partial\Omega)$ of the
boundary $\partial\Omega$. For each dyadic cube $I\in{\cal
D}_{k}(\partial\Omega)$, there are $(d-1)!$ simplices contained in $I$

Let us further introduce the finite dimensional linear space $V^{r}({\cal
T}_{k})$ of functions that are continuous on $\partial\Omega$ and piecewise
polynomial (subordinate to ${\cal T}_{k}$) of order $r$ on $\partial\Omega$.
The linear operator $\overline{S}_{k}$ is a projector onto $V^{r}({\cal
T}_{k})$.

If $g$ is in $C(\partial\Omega)$, we define

$|g|_{H^{1/2}(\partial\Omega)}^{*}:=\left[\frac{1}{\overline{m}^{2}}\sum_{\begin{subarray}{c}i\not=j\\\
\end{subarray}}\frac{|g({\bf z}_{i})-g({\bf z}_{j})|^{2}}{|{\bf z}_{i}-{\bf
z}_{j}|^{d}}\right]^{1/2}$ (6.11)

and prove that the the two semi-norms $|\cdot|_{H^{1/2}(\partial\Omega)}^{*}$
and $|\cdot|_{H^{1/2}(\partial\Omega)}$ are equivalent on $V^{r}({\cal
T}_{k})$.

We begin with the following lemma.

######  Lemma 6.2.

For each pair $T,T^{\prime}\in{\cal T}_{k}$ and each $S\in V^{r}({\cal
T}_{k})$, we have

$\left[\frac{1}{\bar{m}^{2}}\sum_{{\begin{subarray}{c}i\not=j\\\ {\bf
z}_{i}\in\overline{T},\ {\bf
z}_{j}\in\overline{T}^{\prime}\end{subarray}}}\frac{|S({{\bf z}_{i}})-S({{\bf
z}_{j}})|^{2}}{|{{\bf z}_{i}-{\bf
z}_{j}}|^{d}}\right]^{1/2}{\asymp}\left[\intop\limits_{T\times
T^{\prime}}\frac{|S({\bf z})-S({\bf z}^{\prime})|^{2}}{|{\bf z}-{\bf
z}^{\prime}|^{d}}d{\bf z}\ d{\bf z}^{\prime}\right]^{1/2},$ (6.12)

where the constants in $\asymp$ depend only on $r$ and $d$.

###### Proof.

We fix $d,\,r$, and a pair $T,T^{\prime}\in{\cal T}_{k}:={\cal
T}_{k}(\partial\Omega)$. Let us first consider the case when $k\leq 4$. Let
${\cal N}(S):={\cal N}_{T,T^{\prime}}(S)$ denote the expression in (6.12)
involving the sum, and ${\cal N}^{\prime}(S)$ the expression involving the
integral. Both are semi-norms on the finite dimensional linear space
$X=X(T,T^{\prime})$ of functions from $V^{r}({\cal T}_{k})$ restricted to
$T\cup T^{\prime}$. A function $S\in X$ satisfies ${\cal N}(S)=0$ if and only
if $S$ is constant on $T$, $T^{\prime}$ and, additionally, $S$ is constant on
$T\cup T^{\prime}$ when $T$ and $T^{\prime}$ touch. A similar statement holds
for the semi-norm ${\cal N}^{\prime}$. It follows that for any pair
$T,T^{\prime}$, we have

$c{\cal N}(S)\leq{\cal N}^{\prime}({S})\leq C{\cal N}(S),\quad S\in
X(T,T^{\prime}),$ (6.13)

where the constants $c,C$ depend on $T$, $T^{\prime}$, $r$ and $d$, but are
independent of $S$. Since there are only a finite number (depending on $d$) of
pairs of simplicies in ${\cal T}_{k}(\partial\Omega)$, $k\leq 4$, there are
constants $c_{1},C_{1}$ which depend only on $r,d$ for which (6.13) is valid
for all pairs of simplices from ${\cal T}_{k}(\partial\Omega)$, $k\leq 4$.
Then (6.13) holds for all such pairs with $c=c_{1}$, $C=C_{1}$, and we have
completed the proof in this case.

We now consider the case $k>4$. Given a pair $T,T^{\prime}$, we let
$c(T,T^{\prime})$ be the largest number such that the lower inequality in
(6.13) holds uniformly for all $S\in X(T,T^{\prime})$. Similarly
$C(T,T^{\prime})$ is the smallest constant so that the upper inequality of
(6.13) holds uniformly on $X(T,T^{\prime})$. The above argument using
equivalence of semi-norms shows that there always exist positive constants
$c(T,T^{\prime}),C(T,T^{\prime})$. We are left to show that there are
constants $0<c^{*}<C^{*}<\infty$, depending only on $d$ and $r$, for which
$c(T,T^{\prime})\geq c^{*}$ and $C(T,T^{\prime})\leq C^{*}$ holds for all
pairs $T,T^{\prime}$ such that $T,T^{\prime}\in{\cal T}_{k}$ and $k>4$.

We consider two cases. Case 1: $T\cap T^{\prime}\neq\emptyset$. In this case
there is a pair $T_{0},T_{0}^{\prime}\in{\cal T}_{4}$, and a linear mapping
consisting of a translation and a dilation with factor $2^{-k+4}$ that rigidly
transforms $T_{0},T_{0}^{\prime}$ to the pair $T,T^{\prime}$. If we use this
linear transformation to change variables in (6.2) we obtain the validity of
(6.2) for this case with the same constants $c_{1},C_{1}$.

Case 2: $T\cap T^{\prime}=\emptyset$. First note that from the $L_{2}$
discretization of the previous section, we have that there are constants
$c_{2},C_{2}$ depending only on $r,d$ such that for each $S\in V^{r}({\cal
T}_{k}(\partial\Omega))$

$c_{2}\Bigl{[}\intop\limits_{T\times T^{\prime}}|S({\bf z})-S({\bf
z}^{\prime})|^{2}d{\bf z}\ d{\bf
z}^{\prime}\Bigr{]}^{1/2}\leq\Bigl{[}\frac{1}{\overline{m}^{2}}\sum_{\begin{subarray}{c}{\bf
z}\not={\bf z}^{\prime}\in\mathcal{Z}\\\ {\bf z}\in T,\ {\bf z}^{\prime}\in
T^{\prime}\end{subarray}}|S({\bf z})-S({\bf
z}^{\prime})|^{2}\Bigr{]}^{1/2}\leq C_{2}\Bigl{[}\intop\limits_{T\times
T^{\prime}}|S({\bf z})-S({\bf z}^{\prime})|^{2}d{\bf z}\ d{\bf
z}^{\prime}\Bigr{]}^{1/2}.$ (6.14)

Now, on $T\times T^{\prime}$, the expressions $|{\bf z}-{\bf z}^{\prime}|$
that appear in (6.12) are all comparable with absolute constants. Therefore,
we obtain (6.12) holds in this case as well with constants depending only on
$r,d$. ∎

We can now prove that the discrete semi-norm and the true semi-norm of
functions in ${\cal G}$ are comparable.

######  Theorem 6.3.

For any $g\in{\cal G}=Tr(U(\overline{B}))$, where
$\overline{B}=B^{\bar{s}}_{\infty}(L_{2}(\Omega))$ with $\bar{s}>d/2$, and any
$\overline{m}\geq 1$, we have

$|g|_{H^{1/2}(\partial\Omega)}\lesssim|g|^{*}_{H^{1/2}(\partial\Omega)}+{\|g\|_{Tr(\overline{\cal
B})}}\overline{m}^{-\beta},\quad\hbox{and}\quad|g|^{*}_{H^{1/2}(\partial\Omega)}\lesssim|g|_{H^{1/2}(\partial\Omega)}+{\|g\|_{Tr(\overline{\cal
B})}}\overline{m}^{-\beta},$ (6.15)

with

$\beta=\frac{\bar{s}-1}{d-1},$ (6.16)

and constants in $\lesssim$ depending only on $r,d$. Note that the optimal
recovery rate of ${\cal G}$ in the $H^{1/2}(\partial\Omega)$ is
$\bar{m}^{-\beta}$.

###### Proof.

For any $S\in V^{r}({\cal T}_{k})$, we have

$|S|_{H^{1/2}(\partial\Omega)}=\left[\sum_{T,T^{\prime}\in{\cal
T}_{k}}\intop\limits_{T}\intop\limits_{T^{\prime}}\frac{|S({\bf z})-S({\bf
z}^{\prime})|^{2}}{|{\bf z}-{\bf z}^{\prime}|^{d}}\ d{\bf z}\ d{\bf
z}^{\prime}\right]^{1/2}$

and

$|S|_{H^{1/2}(\partial\Omega)}^{*}=\left[\sum_{{\bf z}\not={\bf
z}^{\prime}\in\overline{G}_{k,r}}\frac{1}{\overline{m}^{2}}\frac{|S({\bf
z})-S({\bf z}^{\prime})|^{2}}{|{\bf z}-{\bf
z}^{\prime}|^{d}}\right[^{1/2}\asymp\left[\sum_{T,T^{\prime}\in{\cal
T}_{k}}\sum_{\begin{subarray}{c}{\bf z}\not={\bf
z}^{\prime}\in\overline{G}_{k,r}\\\ {\bf z}\in T,\ {\bf z}^{\prime}\in
T^{\prime}\end{subarray}}\frac{1}{\overline{m}^{2}}\frac{|S({\bf z})-S({\bf
z}^{\prime})|^{2}}{|{\bf z}-{\bf z}^{\prime}|^{d}}\right]^{1/2}.$

Here, in the second expression, we only have equivalence since a given data
site may be used for more than one simplex. Note, however that the constant in
the second equivalence depends only on $d$ and $r$.

In view of Lemma 6.2, we know that there exist constants $0<c<C<\infty$ that
depend only on $d$ and $r$ such that for each $S\in V^{r}({\cal T}_{k})$, we
have

$c|S|_{H^{1/2}(\partial\Omega)}\leq|S|_{H^{1/2}(\partial\Omega)}^{*}\leq
C|S|_{H^{1/2}(\partial\Omega)}.$ (6.17)

In particular this holds for $S=\overline{S}_{k}(g)$ whenever $g\in{\cal G}$.
If we combine (6.17) with (6.10) and observe that
$|g|^{*}_{H^{1/2}(\partial\Omega)}=|\overline{S}_{k}(g)|^{*}_{H^{1/2}(\partial\Omega)}$
we obtain the theorem. ∎

In order to define a discrete norm for $H^{1/2}(\partial\Omega)$, we will also
need a discrete $L_{2}(\partial\Omega)$ norm for the functions in our model
class ${\cal G}$. Using the same $\overline{m}$ data sites ${\cal
Z}:=\overline{G}_{k,r}$, we define

$\|g\|^{*}_{L_{2}(\partial\Omega)}:=\left[\frac{1}{\overline{m}}\sum_{j=1}^{\overline{m}}|g({\bf
z}_{j})|^{2}\right]^{1/2},\quad g\in{\cal G},\quad{\bf
z}_{j}\in\overline{G}_{k,r}.$ (6.18)

Arguing as we have done in §6.1, we have for
$\alpha=\frac{\bar{s}-1/2}{d-1}>\beta$ and any $g\in Tr(\overline{\cal B})$

$\|g\|_{L_{2}(\partial\Omega)}\lesssim\|g\|^{*}_{L_{2}(\partial\Omega)}+\|g\|_{Tr(\overline{\cal
B})}\overline{m}^{-\alpha},\quad\hbox{and}\quad\|g\|^{*}_{L_{2}(\partial\Omega)}\lesssim\|g\|_{L_{2}(\partial\Omega)}+\|g\|_{Tr(\overline{\cal
B})}\overline{m}^{-\alpha}.$ (6.19)

We now define for any continuous function $g$ the discrete
$H^{1/2}(\partial\Omega)$ norm

$\|g\|^{*}_{H^{1/2}(\partial\Omega)}:=\|g\|^{*}_{L_{2}(\partial\Omega)}+|g|^{*}_{H^{1/2}(\partial\Omega)}.$
(6.20)

######  Theorem 6.4.

For any $g\in{\cal G}=Tr(U(\overline{B}))$, where
$\overline{B}=B^{\bar{s}}_{\infty}(L_{2}(\Omega)$ with $\bar{s}>d/2$, and any
$\overline{m}\geq 1$, we have

$\|g\|_{H^{1/2}(\partial\Omega)}\lesssim\|g\|^{*}_{H^{1/2}(\partial\Omega)}+\|g\|_{Tr(\overline{\cal
B})}\overline{m}^{-\beta},\quad\hbox{and}\quad\|g\|^{*}_{H^{1/2}(\partial\Omega)}\lesssim\|g\|_{H^{1/2}(\partial\Omega)}+\|g\|_{Tr(\overline{\cal
B})}\overline{m}^{-\beta},$ (6.21)

where $\beta$ is given in (6.16), and the constants in $\lesssim$ depend only
on $r$ and $d$.

###### Proof.

The theorem follows from the comparisons (6.15), (6.19), and the fact that
$\alpha>\beta$. ∎

##  7 A discrete loss function with error control

In this section, we introduce a discrete loss function ${\cal L}^{*}$ which is
a surrogate for the theoretical loss function ${\cal L}_{T}$ of (5.1). The
advantage of ${\cal L}^{*}$ is that it can be computed directly from the data
$({\bf f},{\bf g})$. Let $k,r$ be fixed and let ${\cal
X}=G_{k,r}\subset\overline{\Omega}$ and ${\cal Z}=\overline{G}_{\bar{k},r}$ be
the data sites of $\partial\Omega$. Note that $\bar{k}$ could be different
from $k$. This fixes $\widetilde{m}=\\#(G_{k,r})$ and
$\overline{m}:=\\#(\overline{G}_{\bar{k},r})$ for the remainder of this
section.

######  Remark 7.1.

In this and the following two sections, the results and the proofs of the
results are simplest in the case $d\geq 3$. When $d=2$, the results and proofs
are clouded by the appearance of log factors. For this reason, we state the
results in correct form (including the log factors) for all cases $d\geq 2$
but give the exposition and proofs only for $d\geq 3$. The proofs for the case
$d=2$ are given in the Appendix, §12.7.

Let $\|\cdot\|^{*}_{L_{\gamma}}$ be the discrete norm, see (6.3), introduced
in §6.1 for the following choices of

$\gamma=\begin{cases}\frac{2d}{d+2},\quad\quad&d\geq 3,\\\ \\\
1+[\log(\widetilde{m})]^{-1},\,\,\,\quad&d=2,\end{cases}$

and let $\|\cdot\|^{*}_{H^{1/2}}$ be the discrete norm, see (6.11), for
$H^{1/2}(\partial\Omega)$ introduced in §6.2. Given the data vectors for
$({\bf f},{\bf g})$, we introduce the discrete loss function

${\cal L}^{*}(v):=\begin{cases}\|f+\Delta
v\|^{*}_{L_{\gamma}(\Omega)}+\|g-v\|^{*}_{H^{1/2}(\partial\Omega)},\quad\quad\quad\quad&d\geq
3,\\\ \\\ [1+\log(\widetilde{m})]\|f+\Delta
v\|^{*}_{L_{\gamma}(\Omega)}+\|g-v\|^{*}_{H^{1/2}(\partial\Omega)},\quad&d=2,\end{cases}$
(7.1)

which is defined whenever $v$ and $\Delta v$ are continuous on
$\overline{\Omega}$. Notice that to compute ${\cal L}^{*}(v)$ we use only the
values of $v$ at the data sites.

We want to show that making ${\cal L}^{*}(v)$ small guarantees that $v$ is a
good approximation to $u$ in the $H^{1}(\Omega)$ norm. As in [4], such a
result requires model class assumptions on $f$ and $g$. Keeping in mind the
remarks of §3.4, we make the model class assumptions ${\cal F}=U({\cal B})$ of
(3.24) and ${\cal G}=Tr(U(\overline{\cal B}))$ of (3.28).

Recall that, in the case $d\geq 3$ and all $0<p\leq\infty$, or $d=2$ and
$1<p\leq\infty$, for a model class ${\cal F}$ of the form (3.24), the uniform
optimal recovery rate for ${\cal F}$ in $H^{-1}(\Omega)$ is
$\asymp\widetilde{m}^{-s/d}$, $\widetilde{m}\geq 1$ where $s$ is the excess
regularity of ${\cal F}$ in $L_{\delta}(\Omega)$. For a model class ${\cal G}$
of the form (3.28), its uniform optimal recovery rate in
$H^{1/2}(\partial\Omega)$ is $\asymp\overline{m}^{-(\bar{s}-1)/(d-1)}$ ,where
$\bar{s}-1$ is the excess regularity in $H^{1/2}(\partial\Omega)$. These model
class assumptions on $f,g$ imply that the solution $u$ to (1) is in a model
class ${\cal U}$ which has a uniform optimal recovery rate
$\max\\{\widetilde{m}^{-\frac{s}{d}},\overline{m}^{-(\bar{s}-1)/(d-1)}\\}$. We
introduce the following notation for a function $v$

$\|v\|_{{\cal U}}:=\max\\{\|\Delta v\|_{{\cal B}},\|Tr(v)\|_{Tr(\overline{\cal
B})}\\}.$ (7.2)

The following theorem bounds the error $\|u-v\|_{H^{1}(\Omega)}$ in terms of
the discrete loss ${\cal L}^{*}(v)$, provided $\Delta v\in{\cal B}$ and
$v\in\overline{\cal B}$.

######  Theorem 7.2.

Let $u$ be the solution to (1) with $f\in{\cal F}=U({\cal B})$ of (3.24) and
$g\in{\cal G}=Tr(U(\overline{\cal B}))$ of (3.28). Given the data $({\bf
f},{\bf g})$ of $f$ and $g$ at grid data sites
$(G_{k,r},\overline{G}_{\bar{k},r})$ with $|G_{k,r}|=\widetilde{m}$ and
$|\overline{G}_{\bar{k},r}|=\overline{m}$, consider the functional ${\cal
L}^{*}$, defined in (7.1). If $v$ is any continuous function in
$H^{1}(\Omega)$, then

$\|u-v\|_{H^{1}(\Omega)}\lesssim{\cal L}^{*}(v)+\left[1+\|v\|_{{\cal
U}}\right]{\cal R}_{{\cal U}}(\widetilde{m},\overline{m}),$ (7.3)

with the constants in $\lesssim$ independent of $u,v$, $\widetilde{m}$ and
$\overline{m}$, where

${\cal R}_{{\cal
U}}(\widetilde{m},\overline{m}):=\begin{cases}\max\\{{\widetilde{m}}^{-\frac{s}{d}},\overline{m}^{-\frac{\bar{s}-1}{d-1}}\\},&d\geq
3,\\\
\max\\{\log({\widetilde{m}}){\widetilde{m}}^{-\frac{s}{2}},\overline{m}^{-(\bar{s}-1)}\\},&d=2.\end{cases}$
(7.4)

###### Proof.

We only prove this for $d\geq 3$. From (1.3), we have

$\displaystyle\|u-v\|_{H^{1}(\Omega)}$ $\displaystyle\lesssim$
$\displaystyle\|f+\Delta
v\|_{H^{-1}(\Omega)}+\|g-Tr(v)\|_{H^{1/2}(\partial\Omega)}\lesssim\|f+\Delta
v\|_{L_{\gamma}(\Omega)}+\|g-Tr(v)\|_{H^{1/2}(\partial\Omega)}$
$\displaystyle\lesssim$ $\displaystyle\left[\|f+\Delta
v\|^{*}_{L_{\gamma}(\Omega)}+\|g-Tr(v)\|^{*}_{H^{1/2}(\Omega)}\right]+\left[\|f+\Delta
v\|_{{\cal B}}\widetilde{m}^{-\frac{s}{d}}+\|g-Tr(v)\|_{Tr(\overline{\cal
B})}\overline{m}^{-\frac{(\bar{s}-1)}{d-1}}\right],$

where in the second inequality we used the continuous embedding of
$L_{\gamma}(\Omega)$ into $H^{-1}(\Omega)$ and in the third inequality we used
the comparisons between continuous and discrete norms (see Lemma 6.1 and
Theorem 6.3). Finally, if we use the facts that $\|f\|_{{\cal B}}\leq 1$ and
$\|g\|_{Tr(\overline{\cal B})}\leq 1$, we complete the proof of the theorem. ∎

Let us make some remarks on this theorem.

######  Remark 7.3.

If we knew in advance the smoothness class of $f$, or at least the value of
$p$, and if this value were $1<p\leq\infty$, then in the case $d=2$ we could
modify appropriately the the loss ${\cal L}^{*}$ so that it is independent of
$\widetilde{m}$. This choice would lead to the actual optimal recovery rates
$\max\\{{\widetilde{m}}^{-\frac{s}{2}},\overline{m}^{-(\bar{s}-1)}\\}$ for
this class (i.e the logarithm will not be present), see the Appendix, §12.7.

######  Remark 7.4.

It is natural that $\|v\|_{{\cal U}}$ enters into the error bounds since we
need some control on $v$ away from the data. However, this means that when we
try to obtain a good approximation to $u$ (in $H^{1}(\Omega)$), through
optimization, then the norm $\|\cdot\|_{{\cal U}}$ will have to enter the
picture [4].

######  Remark 7.5.

This error estimator can be used to monitor the error in any optimization
scheme that attempts to minimize the loss ${\cal L}^{*}(v)$. For example,
suppose that $\Sigma=\Sigma_{n}$ is a linear space of dimension $n$ or a
nonlinear manifold determined by $n$ parameters such that each $S\in\Sigma$ is
continuous and has a continuous Laplacian $\Delta S$. If at any stage of the
optimization procedure we have an $S\in\Sigma_{n}$ then we can use (7.3) to
bound the error of $\|u-S\|_{H^{1}(\Omega)}$. Of course, a good estimate
requires a bound for $\|S\|_{\cal U}$.

######  Remark 7.6.

Note that for fixed values of $\widetilde{m}$ and $\overline{m}$, the losses
${\cal L}$ and ${\cal L}^{*}$ are equivalent up to a constant depending upon
the parameters $\widetilde{m}$ and $\overline{m}$. This means that if the
original loss ${\cal L}$ is driven to $0$, i.e. we are perfectly interpolating
the data, then both losses give the same error control. However, in the
situation where the data is not interpolated, the new loss ${\cal L}^{*}$
gives control on the error according to Theorem 7.2, while the original PINNs
loss ${\cal L}$ does not.

In view of the above remarks, we see that the effectiveness of using
$\Sigma_{n}$ together with the discrete loss ${\cal L}^{*}$ to provide an
approximation to $u$ is not simply governed by the error $E({\cal
U},\Sigma_{n})_{H^{1}(\Omega)}$ of approximating the class ${\cal U}$ by
$\Sigma_{n}$ in the $H^{1}(\Omega)$ but rather by a form of restricted
approximation which involves $\|\cdot\|_{\cal U}$ [4].

##  8 A numerical optimal recovery algorithm

Let ${\cal U}$ be the model class (3.3) determined by the model classes ${\cal
F}$ and ${\cal G}$, where ${\cal F}$ is a maximal model class (3.24), (3.26)
and ${\cal G}$ is a maximal model class (3.28). To build a numerical algorithm
for the recovery of $u\in{\cal U}_{\rm data}$, we want to minimize the
discrete loss ${\cal L}^{*}$ over a suitable set $\Sigma_{n}$. Let $M\geq 1$
and let $\widetilde{u}$ be such that

$\widetilde{u}\in\\{\mathop{\rm argmin}_{\|v\|_{{\cal U}}\leq M}{\cal
L}^{*}(v)\\}.$ (8.1)

Clearly, $u$ is a solution to this problem since $\|u\|_{{\cal U}}\leq 1$ and
${\cal L}^{*}(u)=0$. Thus, ${\cal L}^{*}(\widetilde{u})=0$ and it follows from
(7.3) that

$\|u-\widetilde{u}\|_{H^{1}(\Omega)}\lesssim M{\cal R}_{{\cal
U}}(\tilde{m},\bar{m}),\quad m\geq 1.$ (8.2)

with the implied constant depending only on $r$ and $d$. In other words,
$\widetilde{u}$ provides a near optimal recovery when $d\geq 3$ and a near
optimal recovery up to a $\log m$ factor when $d=2$. Of course, the
minimization problem (8.1) is taken over too large of a set to be numerically
viable and so we would like to utilize minimization over a smaller set.

Let $\Sigma_{n}\subset H^{1}(\Omega)$ be an approximating set. Here, the
primary examples for $\Sigma_{n}$ are linear spaces of finite dimension $n$ or
a nonlinear manifolds depending on $n$ parameters such as NNs. For $M\geq 1$
(typically taken as some fixed number not dependent on $n$) we define the set

$\Sigma_{n}(M):=\\{S\in\Sigma_{n}:\|S\|_{{\cal U}}\leq M\\}\subset\Sigma_{n},$
(8.3)

and consider the solution $\widehat{u}$ of the minimization problem

$\widehat{S}\in\\{\mathop{\rm argmin}_{S\in\Sigma_{n}(M)}{\cal L}^{*}(S)\\}.$
(8.4)

We assume for the moment that $d\geq 3$ and show that $\widehat{S}$ is a near
optimal recovery of $u$ from the given data provided the set $\Sigma_{n}(M)$
is sufficiently good at approximating the elements of our model class ${\cal
U}$. Later, we explain the changes needed for the case $d=2$.

In order to present our results, we denote by

$E(v,\Sigma_{n}(M)):=\begin{cases}\inf_{S\in\Sigma_{n}(M)}\left(\|\Delta
v-\Delta
S\|_{L_{\gamma}(\Omega)}+\|Tr(v-S)\|_{H^{1/2}(\partial\Omega)}\right)&d\geq
3,\\\ \\\ \inf_{S\in\Sigma_{n}(M)}\left((1+\log(\widetilde{m}))\|\Delta
v-\Delta
S\|_{L_{\gamma}(\Omega)}+\|Tr(v-S)\|_{H^{1/2}(\partial\Omega)}\right)&d=2.\end{cases}$
(8.5)

the error of simultaneously approximating $v$ and $\Delta v$. We remark that
in the case $d=2$ and $1<p\leq\infty$, then this analysis can be modified to
remove the logarithm. Then, for the model class ${\cal U}$ of $H^{1}(\Omega)$,
we denote by

$E({\cal U},\Sigma_{n}(M)):=\sup_{v\in{\cal U}}E(v,\Sigma_{n}(M)),$

the error of restricted approximation of the class ${\cal U}$ by elements from
$\Sigma_{n}(M)$.

We have the following theorem.

######  Theorem 8.1.

Let $\Sigma_{n}$ be a set in $H^{1}(\Omega)$ and let $\Sigma_{n}(M)$ be
defined by (8.3)for a fixed value $M\geq 1$. If $\widehat{S}$ is any function
from (8.4) then

$\|u-\widehat{S}\|_{H^{1}(\Omega)}\lesssim M{\cal R}_{{\cal
U}}(\widetilde{m},\overline{m})+E(u,\Sigma_{n}(M)),$ (8.6)

with the constant in $\lesssim$ depending only on $r,d$. Moreover, if the
approximation set $\Sigma_{n}(M)$ is such that $E({\cal U},\Sigma_{n}(M))\leq
CR_{{\cal U}}(\widetilde{m},\overline{m})$ for some $C$ depending only on
$r,d$, then $\widehat{S}$ provides a uniform near optimal recovery (up to a
logarithmic factor in the case $d=2$) for the model class ${\cal U}$.

###### Proof.

We only prove this in the case $d\geq 3$. In the case $d=2$ the same
modifications can be made as before. We use (7.3) for $\widehat{S}$ to obtain

$\|u-\widehat{S}\|_{H^{1}(\Omega)}\lesssim{\cal L}^{*}(\widehat{S})+(1+M){\cal
R}_{{\cal U}}(\widetilde{m},\overline{m}).$ (8.7)

For all $S\in\Sigma_{n}(M)$, Lemma 6.1 and Theorem 6.4 give the bound

${\cal L}^{*}(\widehat{S})\leq{\cal L}^{*}(S)\lesssim\|\Delta u-\Delta
S\|_{L_{\gamma}(\Omega)}+\|Tr(u-S)\|_{H^{1/2}(\partial\Omega)}+CM{\cal
R}_{{\cal U}}(\widetilde{m},\overline{m}),$

with the constant in $\lesssim$ depending only on $r,d$. Taking an infimum
over all $S\in\Sigma_{M}$ gives

${\cal L}^{*}(\widehat{S})\lesssim E(u,\Sigma_{n}(M))+M{\cal R}_{{\cal
U}}(\widetilde{m},\overline{m}).$ (8.8)

When this is inserted in (8.7) we obtain the theorem. ∎

##  9 CPINNs

The previous section has shown that any minimizer of the loss ${\cal L}^{*}$
over $\Sigma_{n}(M)$ provides a near optimal recovery of $u$ provided
$\Sigma_{n}(M)$ is sufficiently good at approximating the elements of the
model class ${\cal U}$. In this section, we specialize these results to the
case where $\Sigma:=\Sigma_{n}$ is a space generated by neural networks with
$n$ parameters and a prescribed architecture. As a primary example, we can
take the activation function $\sigma$ to be ${\rm ReLU}^{k}$ with $k\geq 3$.
This guarantees that the elements in $\Sigma_{n}$ have continuous Laplacians,
i.e., $\Delta S$ is continuous whenever $S\in\Sigma_{n}$. This also guarantees
that ${\cal L}^{*}(S)$ is well defined for each $S\in\Sigma_{n}$. Regarding
the architecture of the network, we consider the case of deep networks which
is the space usually utilized in PINNs. A similar discussion applies to
shallow networks.

Consistent PINNs (CPINNs): We call any algorithm which minimizes the loss
${\cal L}^{*}$ over a neural network space $\Sigma_{n}$ consistent PINNs and
denote it by the acronym CPINNs. The discussion in Remark 7.6 has given the
advantages of using CPINNs over standard PINNs with the loss ${\cal L}$. Here,
we do not necessarily require the minimization take place over the constrained
space $\Sigma_{n}(M)$, unless explicitly stated since the latter invokes
numerical difficulties in implementation.

We do not address the question of designing an optimization algorithm that is
guaranteed to come close to minimizing the loss ${\cal L}^{*}$ over
$\Sigma_{n}$ or $\Sigma_{n}(M)$. We know that the lack of such theoretical
guarantees is a bottleneck in rigorously proving convergence results. We
analyze the performance of any minimizer or near minimizer of the loss without
addressing how such a minimizer is found.

Let $\widehat{u}_{n}$ denote any minimizer of (8.4). Theorem 8.1 guarantees
that $\widehat{u}_{n}$ is a near optimal recovery of $u$ provided that for the
model class ${\cal U}$, we can make $E({\cal U},\Sigma_{n}(M))$ small by
taking $n$ large. Thus, the question is whether $\Sigma_{n}$ has favorable
restricted approximation properties. There is a large literature of
approximation theoretic results concerning how efficiently deep neural
networks can approximate various classes of functions (see [19] and the
references therein). Let us briefly discuss the specific case of approximating
functions from Besov spaces using deep ReLU neural networks. When error is
measured in the $L_{p}$-norm, optimal rates have been obtained [55, 64, 52].
There has also been recent work on the problem when error is measured in
another Besov space, although in this case a smoother activation function must
be taken [63, 62].

On the other hand, the question of restricted approximation that is used in
the present paper has seemingly not been studied. We conjecture that the
techniques developed in [55, 64, 52] can be appropriately modified to give the
same rates for restricted approximation with a fixed suitably large value for
$M$ (independent of $n$). Since the proofs are already quite difficult, we
leave the proof of this conjecture to future work.

##  10 Numerical Illustrations

In this section, we present some numerical experiments in support of the
theoretical results of the previous sections. The main take away from these
theoretical results is that when using optimization of a prescribed loss
function over a set $\Sigma_{n}$ to numerically compute an approximation
$\widehat{u}$ to the solution $u$ of (1), then the loss function ${\cal
L}^{*}$ given by (1.9) as in CPINNs is the proper choice. The purpose of the
numerical experiments of this section is to analyze the effect of using ${\cal
L}^{*}$ rather than the standard loss ${\cal L}$ when taking $\Sigma_{n}$ as a
neural network space.

As a prelude to presenting our numerical experiments, we make some remarks on
minimization of loss functions over $\Sigma_{n}$ that are well known to
practitioners. Since both losses are non convex when viewed as a function of
the parameters of the neural network space $\Sigma_{n}$, there is not a known
mathematical procedure which is proven to converge to a minimizer of the loss.
Typical approaches used in practice, try various initial guesses of the
parameters and use variations of gradient descent with variable step sizes,
i.e., learning rates, in an attempt to find a minimizer of the loss. There is
no a priori guarantee that any specific recipe will work. In other words, one
cannot blindly apply any prescribed procedure with an a priori guarantee of
success.

Another drawback of this approach to solving the PDE is that one cannot test
its efficacy unless the true solution is in some sense known so that one can
measure performance a posteriori. One advantage of using the loss ${\cal
L}^{*}$ in place of ${\cal L}$ is that our theory provides an a posteriori
bound of the error (see (7.3)) that does not require knowledge of the solution
$u$. In other words, when using ${\cal L}^{*}$, one can guarantee an a
posterior bound of the error by computing ${\cal L}^{*}(S)$, where $S$ is the
output of the numerical optimization procedure. Indeed, one would compute
$\|S\|_{\mathcal{U}}$ and from (7.3) we know

$\|u-S\|_{H^{1}(\Omega)}\lesssim\mathcal{L}^{*}(S)+(1+\|S\|_{{\cal
U}})R_{{\cal U}}(\widetilde{m},\overline{m}).$ (10.1)

Here the second term is the recovery error. All quantities on the right side
of (10.1) only involve $S$ and so one has a bound for the error without
knowing $u$.

Let us also note that in practice, for the PDE setting and in other learning
scenarios, one does not try to impose the model class restriction in the
optimization routine since it is computationally too expensive. We follow this
paradigm in our numerical procedures.

Finally, we want to point out that the minimization of the new loss
$\mathcal{L}^{*}$ differs from the standard loss $\mathcal{L}$ in advocating
the use of:

  * •

a discrete $L_{\gamma}(\Omega)$ norm (we take $\gamma=1.1$ when $d=2$) rather
than a discrete $L_{2}(\Omega)$ for the PDE residual;

  * •

a discrete $H^{1/2}(\partial\Omega)$ norm rather than a discrete
$L_{2}(\partial\Omega)$ norm for the boundary condition residual.

Although it is not possible to completely separate the effect of these two
modifications, we consider two different test cases corresponding to a
harmonic solution $u$ with non-vanishing boundary conditions and a singular
solution with vanishing boundary conditions; see Sections 10.2. In each case,
we compare the results of the minimization of the new loss function (referred
as CPINNs) and the standard loss (referred as PINNs).

###  10.1 Practical Algorithm

We now describe the practical algorithm implemented for the minimization of
the losses ${\cal L}^{*}$ and ${\cal L}$. We remark that PINNs algorithms rely
on the inherent assumption that the minimization processes are successful. As
noted above, there is no prescribed recipe for the initial choice of
parameters, the activation functions, the optimization algorithm and the
hyperparameters that will guarantee convergence. Rather, one is left to
experimenting with different choices and employing different heuristics until
satisfactory convergence of this method is obtained.

There is no way to explore all of the possible minimization algorithms.
Moreover, once a particular numerical optimization procedure is chosen to
implement, one is always subjected to the possibility that some other choice
would have been better and demonstrated different results. So,our numerical
results must be viewed from this perspective.

Our approach in this section is to find one procedure that works reasonably
well in minimizing the loss (1.9). We remark that in order to conform the
common practice of minimizing the squared residual in PINNs, we actually used
a loss function obtained by squaring each of the terms in (1.9) and (1.8) in
our experiments, i.e. we minimize the losses ($sq$ is shorthand for squared)

${\cal
L}_{sq}^{*}:=\\!\\!\left[\frac{1}{\widetilde{m}}\sum_{i=1}^{\widetilde{m}}|f({\bf
x}_{i})+\Delta v({\bf
x}_{i})|^{\gamma}\right]^{2/\gamma}\\!\\!\\!\\!+\left[\frac{1}{\overline{m}^{2}}\sum_{\begin{subarray}{c}i,j=1\\\
i\neq j\end{subarray}}^{\overline{m}}\frac{|[g-v]({\bf z}_{i})-[g-v]({\bf
z}_{j})|^{2}}{|{\bf z}_{i}-{\bf
z}_{j}|^{2}}\right]\\!\\!\\!+\left[\frac{1}{\overline{m}}\sum_{j=1}^{\overline{m}}|g({\bf
z}_{j})-v({\bf z}_{j})|^{2}\right],$

and

${\cal L}_{sq}:=\left[\frac{1}{\widetilde{m}}\sum_{i=1}^{\widetilde{m}}[\Delta
v({\bf x}_{i})+f({\bf
x}_{i})]^{2}\right]+\left[\frac{1}{\overline{m}}\sum_{i=1}^{\overline{m}}[v({\bf
z}_{i})-g({\bf z}_{i})]^{2}\right].$

Here we take a uniform grid of collocation points and $\widetilde{m}$ is the
total number of collocation points (where we penalize the residual of the
equation), while $\overline{m}$ is the number of collocation points which lie
on the boundary (where we penalize the residual of the boundary values).

We then apply the exact same procedure to minimize both the original PINNs
loss ${\cal L}_{sq}$ and the CPINNs loss ${\cal L}_{sq}^{*}$. We observe from
the final loss values that our procedure obtains sufficient minimization of
both loss functions, so that our results enable us to compare the effect of
using the two different loss functions ${\cal L}_{sq}$ and ${\cal
L}_{sq}^{*}$. These results demonstrate a significant improvement in accuracy
when minimizing the loss ${\cal L}_{sq}^{*}$ instead of ${\cal L}_{sq}$.

Finally, we remark that the running time of the training the we have
implemented is very high and is not remotely competitive with traditional
numerical methods for this problem. Improving the efficiency of the training
of PINNs on both ${\cal L}_{sq}$ and ${\cal L}_{sq}^{*}$ is an important and
active area of research (see for instance [44, 14, 37, 48, 66]), and
incorporating these recently developed techniques into our algorithm is an
area of future work.

In our experiments, we use a deep neural network with the following structure.
We use a deep network with $L=8$ hidden layers of width $W=100$. The first
layer uses a $\tanh$ activation function, while the remaining layers are
residual layers [33] with the ReLU3 activation function. The network is
initialized as follows:

  * •

The parameters in the first layer are randomly initialized from a normal
distribution with variance $1$.

  * •

For the residual layers with ReLU3 activation function, the parameters are
randomly initialized from a normal distribution with variance

$\sqrt{2}\left(\frac{2}{15}\right)^{1/6}\frac{1}{\sqrt{LW}}.$

In analogy with Xavier and He initialization [26, 32], this ensures that the
variance of the output of each layer gets multiplied by a factor of $(1+2/L)$
in each residual layer, which guarantees that the expected magnitude of the
neural network function at initialization is bounded.

  * •

In the final layer the weights are initialized randomly from a normal
distribution with variance $1/\sqrt{W}$, which ensures that the network output
has the same variance as the output of the final hidden layer.

We train the networks with gradient descent with an initial step size of
$10^{-3}$ decreasing by a factor of $2$ every $4000$ steps and momentum $0.9$
for a total of $40,000$ steps. We also rescale the velocities (i.e. weighted
averages of the previous gradients according to the momentum) in each layer to
have the same maximum as the parameters in that layer, we have found that this
makes the training somewhat more stable and improves the final loss values
obtained. The full code and all experimental details can be found at
https://github.com/jwsiegel2510/consistent-PINNs.

We also note that although the $H^{1/2}(\partial\Omega)$ norm is nonlocal, the
implementation of the discrete $H^{1/2}(\partial\Omega)$-norm present in
${\cal L}_{sq}^{*}$ can easily be vectorized and is fairly efficient in our
practical implementation. Of course, the computation and optimization of the
discrete $H^{1/2}(\partial\Omega)$-norm scales quadratically with the number
of collocation points on the boundary, while the norm used in the standard
loss ${\cal L}_{sq}$ scales only linearly. In our examples, the number of
collocation points in the interior scales quadratically with the number of
points on the boundary, so that using the discrete
$H^{1/2}(\partial\Omega)$-norm on the boundary is comparably expensive to the
penalty on the interior. In addition, the evaluation of the neural network at
the collocation points typically dominates the calculation, rather than the
calculation of the norms or their gradients from these values. However, when
using a large number of collocation points on the boundary optimizing the new
loss function ${\cal L}_{sq}^{*}$ may be less efficient than optimizing the
original ${\cal L}_{sq}$ loss. This can potentially be mitigated by using a
truncated version of the discrete $H^{1/2}(\partial\Omega)$-norm, but we leave
this to future work.

###  10.2 Results

In our first experiment, we solve the Poisson equation

$-\Delta u=f~{}\text{in $\Omega$},~{}u=g~{}\text{on $\partial\Omega$},$ (10.2)

where $f$ and $g$ are chosen so that the exact solution is given by

$u(x,y)=e^{x}\cos(\pi y),\qquad(x,y)\in\Omega=(0,1)^{2}.$

For this solution $u$, we have $f=0$ and the boundary values $g$ are non-zero.
The results of this experiment with different numbers of collocation points
$\widetilde{m}$ can be found in Table 1.

In our second experiment, we solve the Poisson equation (10.2) where $f$ and
$g$ are chosen so that the exact solution is given by

$u(x,y)=1000x(1-x)y(1-y)\left((x-0.5)^{2}+(y-0.5)^{2}\right)^{9/4},\qquad(x,y)\in(0,1)^{2}.$

This solution is chosen such that it is not smooth (although its second
derivatives are still continuous) to test the ability of the algorithm to
handle lower regularity. The results of this experiment can be found in Table
2.

We remark that the $H^{1}$ relative errors in both experiments are reported
with respect to the discrete $H^{1}(\Omega)$-norm

$\|v\|_{H^{1}(\Omega)}^{*}:=\sqrt{\left(\|v\|_{L_{2}(\Omega)}^{*}\right)^{2}+\left(\||\nabla
v|\|_{L_{2}(\Omega)}^{*}\right)^{2}}$

using $500$ points in each direction.

$\widetilde{m}$ | $\overline{m}$ | CPINNs | CPINNs | PINNs | PINNs
---|---|---|---|---|---
$H^{1}$ Rel error | Loss ${\cal L}_{sq}^{*}$ |  $H^{1}$ Rel error | Loss ${\cal L}_{sq}$
25 | 16 | 2.02% | $2.4\cdot 10^{-5}$ | 1.77% | $9\cdot 10^{-6}$
100 | 36 | 0.15% | $3\cdot 10^{-6}$ |  $0.39$% | $5\cdot 10^{-6}$
225 | 56 | 0.06% | $2\cdot 10^{-6}$ |  $0.55$% | $1.9\cdot 10^{-5}$
400 | 76 | 0.07% | $6\cdot 10^{-6}$ |  $0.58$% | $2.5\cdot 10^{-5}$
625 | 96 | 0.11% | $1.1\cdot 10^{-5}$ |  $0.48$% | $2\cdot 10^{-5}$
900 | 116 | 0.07% | $5\cdot 10^{-6}$ | $0.34\%$ | $1.3\cdot 10^{-5}$
Table 1: Results of the numerical experiments for the harmonic solution $u$.
We see that the CPINNs loss function ${\cal L}_{sq}^{*}$ results in about
$3-5$ times lower error compared with the original PINNs loss ${\cal L}_{sq}$.

From these results, we see that using the new loss function ${\cal
L}^{*}_{sq}$ consistently gives errors that are about $3-10$ times smaller
compared with the original PINNs loss. From the final loss values achieved, we
conclude that this is not due to the training algorithm performing better on
the new loss, but is rather due to the fact that ${\cal L}_{sq}^{*}$ gives a
better control on the solution error than ${\cal L}_{sq}$. We remark that by
weighting the terms in ${\cal L}_{sq}$ differently (i.e. by choosing different
values of $\lambda$ in (1.7)), one may be able to improve the results for the
original PINNs loss, as we have not rigorously tested this. However, our
method does not require such tuning as all terms in ${\cal L}_{sq}^{*}$ should
naturally carry the same weight.

We also remark that in our first experiments, the solution is very smooth.
This has the consequence that the optimization error is the dominant source of
error once $\widetilde{m}\geq 100$, i.e. in this regime the sampling error is
already significantly smaller than the optimization error. This is why the
error does not continue to decrease as the number of collocation points is
increased beyond this point. By using a more effective optimization method
which drives the loss ${\cal L}_{sq}^{*}$ this balance can be shifted and
lower error can be obtained. On the other hand, in our second experiment, the
solution $u$ is not particularly smooth. In this case, the sampling error is
dominant all the way up to $\widetilde{m}=900$ collocation points and we see
the error of our method decrease all the way down to this point. This is due
to the fact that the optimal sampling rate decreases much more slowly because
of the non-smoothness of $u$.

$\widetilde{m}$ | $\overline{m}$ | CPINNs | CPINNs | CPINNs | CPINNs
---|---|---|---|---|---
$H^{1}$ Rel error | Loss ${\cal L}_{sq}^{*}$ |  $H^{1}$ Rel error | Loss ${\cal L}_{sq}$
25 | 16 | 52.4% | $1.4\cdot 10^{-4}$ | 46.7% | $2.6\cdot 10^{-3}$
100 | 36 | 22.3% | $1.4\cdot 10^{-3}$ | 28.0% | $6.6\cdot 10^{-3}$
225 | 56 | 9.6% | $4.2\cdot 10^{-3}$ | 26.7% | $1.6\cdot 10^{-2}$
400 | 76 | 7.2% | $5.9\cdot 10^{-3}$ | 19.2% | $7.7\cdot 10^{-3}$
625 | 96 | 5.5% | $6.8\cdot 10^{-3}$ | 25.4% | $2.3\cdot 10^{-2}$
900 | 116 | 3.9% | $3.5\cdot 10^{-3}$ | 23.6% | $2.8\cdot 10^{-2}$
Table 2: Results of the numerical experiments for the non-smooth solution $u$.
We see that the CPINNs loss function ${\cal L}_{sq}^{*}$ results in about
$5-10$ times lower error compared with the original PINNs loss ${\cal
L}_{sq}$.

##  11 Concluding Remarks

We have investigated collocation methods for solving elliptic PDEs of the form
(1). A collocation method seeks to generate an approximation to the true
solution $u$ of (1) from point value information of the right side $f$ and the
boundary value $g$. Under model class assumptions of the form $f\in{\cal F}$
and $g\in{\cal G}$ which are both unit balls of Besov spaces, we showed that
there is an optimal error (called the error of optimal recovery) in recovering
$u$ in the $H^{1}(\Omega)$ norm. We determined this optimal recovery error for
the various model classes.

We then turned to the study of theoretical and numerical algorithms which
would yield an optimal recovery. The most prominent of these methods is to
minimize a certain loss function ${\cal L}$ that measures how well we have fit
the data observations. We have shown that the typical loss function used in
practice (based on least squares minimization) is not the correct loss and
have proposed the correct loss function ${\cal L}^{*}$ defined in (7.1). We
then showed that any minimizer of the correct loss ${\cal L}^{*}$ over a
suitable set $\Sigma_{n}(M)$ yields a near optimal recovery of $u$.

We have then gone on to show that a modfication of PINNs which we call
Consistent PINNs (CPINNs) leads to a provable near optimal recovery algorithm.
The modifications needed are to use the loss ${\cal L}^{*}$ in place of the
loss ${\cal L}$ usually used in PINNs. Secondly, the minimization needs to be
taken over the constrained set $\Sigma_{n}(M)$ rather than all of
$\Sigma_{n}$. Finally, we provide numerical experiments which demonstrate that
the CPINNs loss function significantly outperforms the original PINNs loss
function on elliptic PDEs.

We now wish to make some further comments that one can view as serious caveats
to the above theory. The first and most serious objection to the above
development is that we do not provide a numerical algorithm for finding a
minimizer or near minimizer of the loss ${\cal L}^{*}$ over $\Sigma_{n}(M)$.
This is the same objection that can be made for almost all learning problems
and certainly for PINNs. The usual numerical method used to try and solve the
minimization problem in PINNs, as well as other learning problems, is to use
gradient descent or some modification of gradient descent, to solve the
minimization of the loss. Unfortunately, there is no proof of convergence or
accuracy when applying gradient descent to non-convex losses such as ${\cal
L}$ or ${\cal L}^{*}$. This has not impeded the use of gradient descent in
practice with much empirical success. We have utilized such techniques as well
in our numerical examples.

We do not want to enter into a discussion of the art of using gradient descent
in optimization. However, we want to address two issues that arise relevant to
this paper. The first of these is the fact that we must optimize the loss over
$\Sigma_{n}(M)$ rather than $\Sigma_{n}$ itself to gain the proof of
optimality. This is a serious numerical burden. We have ignored this issue in
our numerical experiments much like one ignores the lack of a provable
gradient descent algorithm. One could implement the restricted approximation
term by adding a penalty to the loss ${\cal L}^{*}$ to guarantee that at each
step of the optimization one remains in the constrained set $\Sigma_{n}(M)$
(see [4]) but this does not ease the computational burden and one is still
faced with a lack of convergence of the numerical optimization.

Finally, let us point out another issue which is understood by practitioners
but not usually pointed out in the theory. Neural network spaces are a very
unstable manifold. By this we mean that changing parameters very slightly can
have a huge effect on the output of the neural network. That is, the mapping
from parameters to output in neural networks is very unstable. This
instability appears more in deep networks than in shallow networks. Avoiding
this lack of stability is an art in practice that manifests itself in
judiciously choosing the starting parameters and the learning rates (step size
in gradient descent). The theoretical aspects of this lack of stability, and
how to avoid the instability theoretically, is addressed in a series of recent
papers [8, 29, 28, 30]. There remains a serious gap between theoretical
algorithms for learning and their numerical implementation.

##  12 Appendix

In this appendix, we provide proof of the results stated in sections
$\S\ref{S:Besov}$ and §3.

###  12.1 Local approximation by polynomials

There are many important results on the approximation of functions in Besov
classes. We will use approximation by piecewise polynomials. We begin by
describing local polynomial approximation.

For any integer $r\geq 1$, we let ${\cal P}_{r}^{d}:={\cal P}_{r}$ denote the
linear space of algebraic polynomials of order $r$ (total degree $r-1$),
namely,

${\cal P}_{r}:=\left\\{\sum_{|{{\bf k}}|_{1}<r}a_{{\bf k}}{\bf x}^{{\bf
k}},\,\,a_{{\bf k}}\in\mathbb{R}\right\\},\quad\hbox{where}\quad{\bf x}^{{\bf
k}}:=x_{1}^{k_{1}}\cdots x_{d}^{k_{d}},\quad{\bf k}:=(k_{1},\dots,k_{d}),\quad
k_{j}\geq 0,\quad|{\bf k}|_{1}:=\sum_{j=1}^{d}k_{j}.$

Note that $\omega_{r}(P,t)_{p}=0$, $t\geq 0$, for all $P\in{\cal P}_{r}$. If
$I$ is any cube in $\mathbb{R}^{d}$ and $f\in L_{p}(I)$, $0<p\leq\infty$, we
denote by

$E_{r}(f,I)_{p}:=\inf_{P\in{\cal P}_{r}}\|f-P\|_{L_{p}(I)},$ (12.1)

the error of approximation of $f$ on $I$ in the $\|\cdot\|_{L_{p}(I)}$ norm by
the elements of ${\cal P}_{r}$. A well known result in approximation theory,
commonly referred to as Whitney’s theorem, says that for any $f\in L_{p}(I)$
with $I$ a cube with sidelength $\ell_{I}$, we have

$cE_{r}(f,I)_{p}\leq\omega_{r}(f,\ell_{I})_{L_{p}(I)}\leq CE_{r}(f,I)_{p},$
(12.2)

with the constants $c,C$ depending only on $r,d$ and also $p$ if $p$ is close
to $0$. Whitney’s theorem usually only refers to the lower inequality in
(12.2). However, the upper inequality follows trivially since

$\omega_{r}(f,\ell_{I})_{L_{p}(I)}=\omega_{r}(f-P,\ell_{I})_{L_{p}(I)}\leq
C\|f-P\|_{L_{p}(I)},$

holds for any polynomial $P\in{\cal P}_{r}$.