# When Do Extended Physics-Informed Neural Networks (XPINNs) Improve
Generalization?111Published in SIAM Journal on Scientific Computing (SISC)

Zheyuan Hu Department of Computer Science, National University of Singapore,
Singapore, 119077
(mailto:e0792494@u.nus.edue0792494@u.nus.edu,mailto:kenji<EMAIL_ADDRESS>Ameya D. Jagtap Division of Applied Mathematics, Brown University, Providence,
RI 02912, USA (mailto:ameya_jagtap<EMAIL_ADDRESS>mailto:george_karniadakis<EMAIL_ADDRESS>George Em
Karniadakis33footnotemark: 3 Kenji Kawaguchi22footnotemark: 2

###### Abstract

Physics-informed neural networks (PINNs) have become a popular choice for
solving high-dimensional partial differential equations (PDEs) due to their
excellent approximation power and generalization ability. Recently, Extended
PINNs (XPINNs) based on domain decomposition methods have attracted
considerable attention due to their effectiveness in modeling multiscale and
multiphysics problems and their parallelization. However, theoretical
understanding on their convergence and generalization properties remains
unexplored. In this study, we take an initial step towards understanding how
and when XPINNs outperform PINNs. Specifically, for general multi-layer PINNs
and XPINNs, we first provide a prior generalization bound via the complexity
of the target functions in the PDE problem, and a posterior generalization
bound via the posterior matrix norms of the networks after optimization.
Moreover, based on our bounds, we analyze the conditions under which XPINNs
improve generalization. Concretely, our theory shows that the key building
block of XPINN, namely the domain decomposition, introduces a tradeoff for
generalization. On the one hand, XPINNs decompose the complex PDE solution
into several simple parts, which decreases the complexity needed to learn each
part and boosts generalization. On the other hand, decomposition leads to less
training data being available in each subdomain, and hence such model is
typically prone to overfitting and may become less generalizable. Empirically,
we choose five PDEs to show when XPINNs perform better than, similar to, or
worse than PINNs, hence demonstrating and justifying our new theory.

##  1 Introduction

Deep learning has revolutionized numerous fields in computer science, such as
computer vision and natural language process. Recently, deep neural networks
have also been employed to solve partial differential equations (PDEs) and
integrated into the field of scientific computing, thanks to their unique
optimization [1, 2, 3, 4] and generalization [5] abilities. Physics-informed
neural networks (PINNs) [6] are among the most popular approaches with a wide
variety of successful applications, including heat transfer problems [7],
thrombus material properties [8], nano-optics [9], and fluid mechanics [10,
11]. PINNs are used as surrogates of a target solution for solving PDEs, and a
solution is found by searching for the best parameters of PINNs that satisfy
the physical laws governed by the PDEs. A more recent work [12] proposed
Extended PINNs (XPINNs), which improves on PINNs by employing a domain
decomposition method for partitioning the PDE problem into several sub-
problems on subdomains, where each sub-problem can be solved by individual
networks called as sub-PINNs. In XPINNs, the continuity of the PINN functions
between each subdomain is maintained. XPINNs facilitate parallel computing,
accelerate convergence, and improve generalization empirically. Despite great
progress in applications, currently no theoretical understanding exists on
when and how XPINNs are better than PINNs.

Recently, some works on theoretically understanding of PINNs have emerged [13,
14, 15, 16]. For two-layer networks, Luo and Yang [13] derived prior and
posterior generalization bounds for PINNs based on Barron space and Rademacher
complexity, whereas Lu et al. [14] provided prior error estimates based on
Barron spaces with the softplus activation. For multi-layer networks, Mishra
and Molinaro [15] introduced abstract formalism and stability properties of
the underlying PDEs to derive generalization bounds, and Shin et al. [16] used
the Holder continuity constant to bound the generalization of PINNs. While
these previous bounds significantly advanced our theoretical understanding of
PINNs, we cannot rely on them to study the advantages and disadvantages of
multi-layer XPINNs over PINNs. This is because the previous studies focus on
PINNs, and the previous bounds either only apply to two-layer networks or
depend on variables that are hard to be computed analytically or numerically.
For example, the Holder continuity is often difficult to compute efficiently
and the assumption of Holder continuity regularization is not widely adopted
in practice. Accordingly, it is necessary to employ different approaches to
derive new generalization bounds for the multi-layer XPINN in order to
understand its advantages and limitations.

In this study, we provide an initial step towards understanding how and when
XPINNs improve generalization capabilities of PINNs by proving new
generalization bounds for multi-layer XPINNs and PINNs. Specifically, we first
discuss the Barron space theory for multi-layer networks to define the
function space of neural networks. We then derive a prior generalization bound
for PINNs with the complexity of a target function measured via its Barron
norm without any further assumption. Furthermore, we derive Rademacher
complexity bounds of PINNs via capacity controls based on the spectral norm
and the (2,1) norm, which are then used for our posterior generalization
bounds for PINNs. We then extend these bounds of PINNs to those of XPINNs by
applying the bounds to each of the subdomains in XPINNs and combine them to
form the final results. Overall, our theoretical results predict that in terms
of generalization, the advantages and disadvantages of XPINNs come from the
tradeoff between the reduction in the complexity of decomposed target
functions (within each subdomain) and the increase in the over-fitting due to
less available training data (in each subdomain). That is, the domain
decomposition of XPINNs can make a target function in a subdomain to be less
complex than the whole target function, resulting in a reduction in a
complexity measure, whereas each sub-network tend to utilize less than the
entire available training data. To illustrate when and how XPINNs improve
generalization based on our theory, we first provide analytical examples,
where we mathematically compute and compare the prior bounds of XPINNs and
PINNs. Furthermore, we adopt five PDEs to numerically demonstrate our
posterior bounds via experiments. Both analytical examples and experimental
observations confirm our theoretical prediction and deepen our understanding,
demonstrating that the two factors in our generalization bounds lead to a
tradeoff, leading to different performances of XPINNs over PINNs on various
tasks.

The remainder of this paper is arranged as follows. In Section 2, we provide
properties, background and assumptions on PDEs, PINNs, and the function space
for multi-layer neural networks. In Section 3, our main generalization results
(both prior and posterior bounds) are presented. In Section 4, discussion on
theoretical analysis as well as analytical examples are introduced. In Section
5, extensive experiments are conducted to numerically demonstrate our theory.

##  2 Preliminaries

In this section, we present introductory facts for PDEs, neural networks, as
well as PINNs and XPINNs. We use bold-faced lowercase letters to denote
vectors, and capital letters to denote matrices and network parameters. Given
a vector $\boldsymbol{v}$, we denote its Euclidean norm by
$\|\boldsymbol{v}\|$, while $\|\cdot\|_{p}$ refers to the $p$-norms. For
matrix norms, we denote the spectral norm by $\|\cdot\|_{2}$ and $l_{p,q}$
norms by
$\|\boldsymbol{W}\|_{p,q}=(\sum_{j}(\sum_{k}|W_{j,k}|^{p})^{q/p})^{1/q}$.
Following convention, we define $\inf$ of a set $S$ to be the infimum of the
subset $S$ of $\overline{R}$ (the set of affinely extended real numbers);
e.g., the infimum of the empty set is infinity.

###  2.1 PDE Problem

In this paper, we consider PDEs defined on the bounded domain
$\Omega=(-1,1)^{d}$. More specifically, the PDEs under consideration are in
the form of

$\displaystyle\mathcal{L}u^{*}$ $\displaystyle=f\ \text{in}\ \Omega,\qquad
u^{*}=g\ \text{on}\ \partial\Omega,$ (1)

where $\mathcal{L}$ is the differential operator characterizing the PDE,
$\partial\Omega$ is the boundary of the set $\Omega$,
$f:\boldsymbol{x}=(x_{1},\dots,x_{d})\in{\Omega}\longmapsto
f(\boldsymbol{x})\in\mathbb{R}$ and
$g:\boldsymbol{x}=(x_{1},\dots,x_{d})\in\partial{\Omega}\longmapsto
g(\boldsymbol{x})\in\mathbb{R}$ are given functions, and the function
$u^{*}:\boldsymbol{x}=(x_{1},\dots,x_{d})\in\overline{\Omega}\longmapsto
u^{*}(\boldsymbol{x})\in\mathbb{R}$ is the unknown solution of PDEs with its
domain $\overline{\Omega}=\Omega\cup\partial\Omega$.

###  2.2 PINN and XPINN

In this subsection, we introduce neural network-based PDE solvers PINNs and
XPINNs. Specifically, in PINNs we optimize neural networks via gradient-based
algorithms to enable the network functions to satisfy the data and the
physical laws governed by the PDEs. Given $n_{b}$ boundary training points
$\left\\{\boldsymbol{x}_{b,i}\right\\}_{i=1}^{n_{b}}\subset\partial\Omega$ and
$n_{r}$ residual training points
$\left\\{\boldsymbol{x}_{r,i}\right\\}_{i=1}^{n_{r}}\subset\Omega$, we
approximate the true PDE solution
$u^{*}:\overline{\Omega}\rightarrow\mathbb{R}$ by the PINN function
$u_{\boldsymbol{\theta}}$ parameterized by $\boldsymbol{\theta}$ via
minimizing the empirical loss composed of the boundary loss and the residual
loss, as given below.

$R_{S}(\boldsymbol{\theta})=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}{|u_{\boldsymbol{\theta}}(\boldsymbol{x}_{b,i})-g(\boldsymbol{x}_{b,i})|}^{2}+\frac{1}{n_{r}}\sum_{i=1}^{n_{r}}{|\mathcal{L}u_{\boldsymbol{\theta}}(\boldsymbol{x}_{r,i})-f(\boldsymbol{x}_{r,i})|}^{2},$
(2)

where the first term is included to force the network to satisfy boundary
conditions, while the second term forces the network to satisfy the physical
laws described by the PDEs.

XPINN is an extension of PINN, obtained by decomposing the whole domain
$\overline{\Omega}$ into several subdomains, mapped to several sub-PINNs. The
continuity between each sub-nets is maintained via the interface loss function
and the final solution of XPINN is the combination and ensemble of all sub-
nets, where each of them is responsible for prediction on one subdomain. More
specifically, the original domain $\Omega$ is decomposed into $N_{D}$
subdomains as $\Omega=\cup_{i=1}^{N_{D}}\Omega_{i}$. The loss of XPINN
contains the sum of losses for the sub-nets, which consist of boundary loss
and residual loss, plus the interface loss using points on
$\partial\Omega_{i}\cap\partial\Omega_{j}$, where
$i,j\in\left\\{1,2,...,N_{D}\right\\}$ such that
$\partial\Omega_{i}\cap\partial\Omega_{j}\neq\emptyset$ to maintain the
continuity between the two sub-nets $i$ and $j$. Mathematically, the XPINN
loss for the $i$-th subdomain is

$R_{S}^{i}(\boldsymbol{\theta}^{i})+\lambda_{I}\sum_{i,j:\partial\Omega_{i}\cap\partial\Omega_{j}\neq\emptyset}R_{I}(\boldsymbol{\theta}^{i},\boldsymbol{\theta}^{j}),$
(3)

where $\lambda_{I}\geq 1$ is the weight controlling the strength of interface
loss, $\boldsymbol{\theta}^{i}$ is the parameters for subdomain $i$, and each
$R_{S}^{i}(\boldsymbol{\theta})$ is the PINN loss for subdomain $i$ containing
boundary and residual losses, i.e.

$R_{S}^{i}(\boldsymbol{\theta}^{i})=\frac{1}{n_{b,i}}\sum_{j=1}^{n_{b,i}}{|u_{\boldsymbol{\theta}^{i}}(\boldsymbol{x}^{i}_{b,j})-g(\boldsymbol{x}^{i}_{b,j})|}^{2}+\frac{1}{n_{r,i}}\sum_{j=1}^{n_{r,i}}{|\mathcal{L}u_{\boldsymbol{\theta}^{i}}(\boldsymbol{x}^{i}_{r,j})-f(\boldsymbol{x}^{i}_{r,j})|}^{2},$
(4)

where $n_{b,i}$ and $n_{r,i}$ are the number of boundary points and residual
points in subdomain $i$ respectively, $\boldsymbol{x}^{i}_{b,j}$ and
$\boldsymbol{x}^{i}_{r,j}$ are the $j$-th boundary and residual training
points in subdomain $i$, respectively. Moreover,
$R_{I}(\boldsymbol{\theta}^{i},\boldsymbol{\theta}^{j})$ is the interface loss
between the $i$-th and $j$-th subdomains based on several interface training
points
$\\{\boldsymbol{x}^{ij}_{I,k}\\}_{k=1}^{n_{I,ij}}\subset\partial\Omega_{i}\cap\partial\Omega_{j}$

$\displaystyle R_{I}(\boldsymbol{\theta}^{i},\boldsymbol{\theta}^{j})$
$\displaystyle=\frac{1}{n_{I,ij}}\sum_{k=1}^{n_{I,ij}}[|u_{\boldsymbol{\theta}^{i}}(\boldsymbol{x}^{ij}_{I,k})-\\{\\{u_{\boldsymbol{\theta}^{avg}}\\}\\}|^{2}+$
(5) $\displaystyle\
|(\mathcal{L}u_{\boldsymbol{\theta}^{i}}(\boldsymbol{x}^{ij}_{I,k})-f_{i}(\boldsymbol{x}^{ij}_{I,k}))-(\mathcal{L}u_{\boldsymbol{\theta}^{j}}(\boldsymbol{x}^{ij}_{I,k})-f_{j}(\boldsymbol{x}^{ij}_{I,k}))|^{2}],$

where
$\\{\\{u_{\boldsymbol{\theta}^{avg}}\\}\\}=u_{avg}\coloneqq({u_{\boldsymbol{\theta}^{i}}(\boldsymbol{x}^{ij}_{I,k})+u_{\boldsymbol{\theta}^{j}}(\boldsymbol{x}^{ij}_{I,k})})/{2}$,
$n_{I,ij}$ is the number of interface points between the $i$-th and $j$-th
subdomains, while $\boldsymbol{x}^{ij}_{I,k}$ is the $k$-th interface points
between them. The first term is the average solution continuity between the
$i$-th and the $j$-th sub-nets, while the second term is the residual
continuity condition on the interface given by the $i$-th and the $j$-th sub-
nets.

We also notice a recent paper on improving the training and generalization of
XPINN [17], which includes the following additional interface regularization
term for XPINN:

$R_{A}(\boldsymbol{\theta}^{i},\boldsymbol{\theta}^{j})=\frac{1}{n_{I,ij}}\sum_{k=1}^{n_{I,ij}}\sum_{m=1}^{d}\left|\frac{\partial
u_{\boldsymbol{\theta}^{i}}(\boldsymbol{x}^{ij}_{I,k})}{\partial\boldsymbol{x}_{m}}-\frac{\partial
u_{\boldsymbol{\theta}^{j}}(\boldsymbol{x}^{ij}_{I,k})}{\partial\boldsymbol{x}_{m}}\right|^{2},$
(6)

where $d$ is the problem dimension, i.e., $\boldsymbol{x}\in\mathbb{R}^{d}$.
The additional interface condition forces the continuity of the first order
derivatives between two sub-nets. In our experiment on the Poisson equation in
section 5.4, which includes residual discontinuity, we shall show how this
additional term improves XPINN via decreasing errors near the interface.

For our discussion on generalization, besides the training losses above, the
testing loss evaluating generalization ability is defined as

$R_{{D}}(\boldsymbol{\theta})=\mathbb{E}_{\text{Unif}(\partial\Omega)}{|u_{\boldsymbol{\theta}}(\boldsymbol{x})-g(\boldsymbol{x})|}^{2}+\mathbb{E}_{\text{Unif}(\Omega)}{|\mathcal{L}u_{\boldsymbol{\theta}}(\boldsymbol{x})-f(\boldsymbol{x})|}^{2},$
(7)

where $\text{Unif}(A)$ is the uniform distribution on a set $A$. Note that in
the definition of the population loss, the interface losses of XPINNs are
excluded to compare PINNs and XPINNs with the same quantity – the
generalization bound for the boundary and residual terms. The beneficial
effect of the interface loss is in improving the generalization of the
boundary and residual terms instead of helping the generalization of the
interface term itself. This is because the interface allows the sub-net in the
subdomain $\Omega_{i}$ to implicitly use samples from other subdomain
$\Omega_{j}$ ($i\neq j$) for regularization through the continuity.

Lastly, we denote the boundary empirical loss and the residual empirical loss
by$R_{S\cap\partial\Omega}$ and $R_{S\cap\Omega}$, respectively, and their
population versions by$R_{D\cap\partial\Omega}$ and $R_{D\cap\Omega}$,
inspired by the fact that boundary points are on $\partial\Omega$ and that
residual points are in $\Omega$. Specifically, their mathematical definitions
are:

$\displaystyle R_{S\cap\partial\Omega}$
$\displaystyle=\frac{1}{n_{b}}\sum_{i=1}^{n_{b}}{|u_{\boldsymbol{\theta}}(\boldsymbol{x}_{b,i})-g(\boldsymbol{x}_{b,i})|}^{2}.\quad
R_{S\cap\Omega}=\frac{1}{n_{r}}\sum_{i=1}^{n_{r}}{|\mathcal{L}u_{\boldsymbol{\theta}}(\boldsymbol{x}_{r,i})-f(\boldsymbol{x}_{r,i})|}^{2}.$
(8) $\displaystyle R_{D\cap\partial\Omega}$
$\displaystyle=\mathbb{E}_{\text{Unif}(\partial\Omega)}{|u_{\boldsymbol{\theta}}(\boldsymbol{x})-g(\boldsymbol{x})|}^{2}.\quad
R_{D\cap\Omega}=\mathbb{E}_{\text{Unif}(\Omega)}{|\mathcal{L}u_{\boldsymbol{\theta}}(\boldsymbol{x})-f(\boldsymbol{x})|}^{2}.$

###  2.3 Neural Networks

In this subsection, we define neural networks and their related properties.

######  Definition 2.1.

(Neural Network). A deep neural network (DNN)
$u_{\boldsymbol{\theta}}:\boldsymbol{x}=(x_{1},\dots,x_{d}$
$)\in\overline{\Omega}\longmapsto
u_{\boldsymbol{\theta}}(\boldsymbol{x})\in\mathbb{R}$, parameterized by
$\boldsymbol{\theta}$ of depth $L$ is the composition of $L$ linear functions
with element-wise non-linearity $\sigma$, is expressed as below.

$u_{\boldsymbol{\theta}}(\boldsymbol{x})=\boldsymbol{W}^{L}\sigma(\boldsymbol{W}^{L-1}\sigma(\cdots\sigma(\boldsymbol{W}^{1}\boldsymbol{x})\cdots),$
(9)

where $\boldsymbol{x}\in\mathbb{R}^{d}$ is the input, and
$\boldsymbol{W}^{l}\in\mathbb{R}^{m_{l}\times m_{l-1}}$ is the weight matrix
at $l$-th layer with $d=m_{0}$ and $m_{L}=1$. The parameter vector
$\boldsymbol{\theta}$ is the vectorization of the collection of all
parameters. We denote $h$ as the maximal width of the neural network, i.e.,
$h=\max(m_{L},\cdots,m_{0})$.

We consider DNNs without bias because one can always set
$\boldsymbol{x}\leftarrow[\boldsymbol{x},1]$ to involve the bias term. Note
that the non-linearity $\sigma$ is Lipschitz continuous with Lipschitz
constant 1 in $\Omega$. The widely adopted ReLU activation function
ReLU$(\boldsymbol{x})=\max(0,\boldsymbol{x})$ cannot be used in our setting
due to its non-differentiability.

Because the neural network $u_{\boldsymbol{\theta}}$ is always differentiated
with respect to its input in the residual losses in PINNs, we introduce their
expressions as follows:

$\frac{\partial
u_{\boldsymbol{\theta}}(\boldsymbol{x})}{\partial\boldsymbol{x}}=\boldsymbol{W}^{L}\cdot\boldsymbol{\Phi}^{L-1}\boldsymbol{W}^{L-1}\cdot\dots\cdot\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1}\in\mathbb{R}^{d},$
(10)
$\displaystyle\frac{\partial^{2}u_{\boldsymbol{\theta}}(\boldsymbol{x})}{\partial\boldsymbol{x}^{2}}$
$\displaystyle=\\{\sum_{l=1}^{L-1}(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})\text{diag}(\boldsymbol{\Psi}^{l}\cdots\boldsymbol{\Psi}^{1}\boldsymbol{W}^{1}_{:,j})(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})\\}_{1\leq
j\leq d},$ (11)

where
$\boldsymbol{\Phi}^{l}=\text{diag}[\sigma^{\prime}(\boldsymbol{W}^{l}\sigma(\boldsymbol{W}^{l-1}\sigma(\cdots\sigma(\boldsymbol{W}^{1}\boldsymbol{x}))))]\in\mathbb{R}^{m_{l}\times
m_{l}}$, and
$\boldsymbol{\Psi}^{l}=\text{diag}[\sigma^{\prime\prime}(\boldsymbol{W}^{l}\sigma($
$\boldsymbol{W}^{l-1}\sigma(\cdots\sigma(\boldsymbol{W}^{1}\boldsymbol{x}))))]\in\mathbb{R}^{m_{l}\times
m_{l}}$. In the Appendix, we provide detailed computation of the derivatives.

###  2.4 Generalized Barron Space

In this subsection, we introduce the generalized Barron space [18], which is a
natural building block to construct a function space of multi-layer deep
networks. This will facilitate our study of their approximation and
generalization properties. We begin by presenting some mathematical
background.

Let $X$ be a Banach space such that $X$ embeds continuously into the space
$C^{2,1}(\overline{\Omega})$ of functions on $\overline{\Omega}$. We further
assume that the closed unit ball $B^{X}$ in $X$ is closed in the topology of
$C^{2}(\overline{\Omega})$.

Because $X$ embeds continuously into the space $C^{2,1}(\overline{\Omega})$,
the Lipschitz constants of functions in $B^{X}$ and their derivatives up to
second order are bounded by the same constant and thus the subset $B^{X}$ is
uniformly $2$-equicontinuous (see the supplementary for definition).

The subset $B^{X}$ is pre-compact, i.e., its closure is compact, in the
separable Banach space $C^{2}(\overline{\Omega})$, because $\overline{\Omega}$
is compact, and that $C^{0}(\overline{\Omega})$ is separable. Since $B^{X}$ is
$C^{2}$-closed and pre-compact, it is compact and is a Polish space in
particular.

Let $\mu$ be a finite signed measure on the Borel $\sigma$-algebra of $B^{X}$,
with respect to the $C^{0}$-norm. Then $\mu$ is a signed Radon measure. We
therefore consider the infinite-dimensional vector function version of the
activation $\sigma:B^{X}\rightarrow C^{2}(\overline{\Omega})$ given by
$\sigma:g\longmapsto(\sigma\circ g)\text{ for each }g\in B^{X}$, where
$(a\circ b)$ represents the composition of functions $a$ and $b$. Then, the
infinite-dimensional vector function $\sigma$ is continuous due to the
Lipschtiz continuity of the dimensional wise version of $\sigma$. Thus, the
infinite-dimensional vector function $\sigma$ is strongly measurable (the
preimage of Borel sets in $B^{X}$ are Borel sets in
$C^{0}(\overline{\Omega})$), and $\mu$-integrable in the sense of the Bochner
integral. The above construction of new function class containing all
$\sigma\circ g$ from the class of $g$ in $B^{X}$ can be formalized as follows.

######  Definition 2.2.

(Generalized Barron Space) The generalized Barron space modeled on $X$
associated with the non-linearity $\sigma$ is a normed space
$(\mathcal{B}_{X,\Omega},\|\cdot\|_{X,\Omega})$ with

$\displaystyle\mathcal{B}_{X,\Omega}$ $\displaystyle=\left\\{f\in
C^{2}(\overline{\Omega}):\|f\|_{X,\Omega}<\infty\right\\},\text{ and }$ (12)
$\displaystyle\|f\|_{X,\Omega}$
$\displaystyle=\inf\left\\{\|\mu\|_{\mathcal{M}(B^{X})}:\mu\in\mathcal{M}(B^{X})\
\text{s.t.}\ f=f_{\mu}\ \text{on}\ \overline{\Omega}\right\\},$

where $\mathcal{M}(B^{X})$ denotes the space of Radon measures on $B^{X}$ and
$f_{\mu}=\int_{B^{X}}\sigma(g)d\mu(g)$. Here, the integral represents the
Bochner integral with $g\in B^{X}$.

For example, if $X$ is the space of linear functions from $\mathbb{R}^{d}$ to
$\mathbb{R}$ (which is isomorphic to $\mathbb{R}^{d+1}$), then the generalized
Barron space modeled on $X$ is the (usual) Barron space for two-layer neural
networks [19]. If $X$ is the (usual) Barron space of two-layer neural
networks, then the generalized Barron space modeled on $X$ is the space for
three-layer neural networks. That is, we can construct the generalized Barron
space of $L$-layer networks from that of $(L-1)$-layer networks by recursively
applying its definition. This recursive construction leads to the tree-like
function space (Definition 2.3), which is a function space of multi-layer
neural networks (Theorem 2.1), as given below.

######  Definition 2.3.

(Tree-Like Function Space for Deep Networks) The tree-like function space
$\mathcal{W}^{L}(\Omega)$ of depth $L$ is recursively defined by
$\mathcal{W}^{l}(\Omega)=\mathcal{B}_{\mathcal{W}^{l-1}(\Omega),\Omega}$ for
all $l\in\\{2,3,\dots,L\\}$ where $\mathcal{W}^{1}(\Omega)$ is the space of
linear functions from $\mathbb{R}^{d}$ to $\mathbb{R}$.

The term “tree like function space” is taken from one of the most related
works [18], and a famous paper on neural network approximation [20].
Intuitively, the neural networks resemble the tree structure, where each
neuron in the network is a node in the graph, and the edges of the graph
connect the neurons in the previous and the present layers. Our “tree like
function space” has nothing to do with the tree functions in graph theory. Our
definition contains a vast function space covering the majority of PDE
solutions.

######  Theorem 2.1.

(Embedding of Finite Networks). The tree-like function space
$\mathcal{W}^{L}(\Omega)$ contains all finite multi-layer networks
$u_{\boldsymbol{\theta}}(\boldsymbol{x})$ of depth $L$ satisfying
$\|\boldsymbol{W}^{l}\|_{1,\infty}\leq 1,1\leq l\leq L-1$, and the Barron norm
of networks satisfies
$\|u_{\boldsymbol{\theta}}\|_{\mathcal{W}^{L}(\Omega)}\leq\|\boldsymbol{W}^{L}\|_{1,\infty}$.

Theorem 2.1 shows that the tree-like function space constructed via the
generalized Barron space indeed contains the class of multi-layer neural
networks, the norm of which is controlled by the $1,\infty$ matrix norm of
their parameters. The following is a list of basic properties of the
generalized Barron space, which also holds for the tree-like function space
and justifies our recursive construction of the tree-like function space:

######  Theorem 2.2.

(Property of Generalized Barron Spaces). The following two statements are
true. (1) The generalized Barron space is complete in the metric defined by
the generalized Barron norm $\|\cdot\|_{X,\Omega}$: i.e., the generalized
Barron space is a Banach space. (2) $\mathcal{B}_{X,\Omega}$ embeds
continuously into $C^{2,1}(\overline{\Omega})$ and the closed unit ball of
$\mathcal{B}_{X,\Omega}$ is a closed subset of $C^{2}(\overline{\Omega})$.

The last property indicates that $\mathcal{B}_{X,\Omega}$ satisfies the same
properties, which we imposed on $X$ during the construction, i.e., we can
repeat the construction and consider $B_{\mathcal{B}_{X,\Omega}}$, hence,
ensuring the validity of the recursive construction in the tree-like function
space. As universal approximators, neural networks can also approximate
arbitrary Barron functions accurately.

######  Theorem 2.3.

(Approximation Properties of Tree-Like Functions). Let $\mathbb{P}$ be a
probability measure with compact support in $\Omega$, and $\mathbb{Q}$ be a
probability measure with compact support in $\partial\Omega$. Then for any
$L\geq 1,f\in\mathcal{W}^{L}(\Omega)$ and $m\in\mathbb{N}$, there exists a
neural network $u_{\boldsymbol{\theta}}(\boldsymbol{x})$ of depth $L$, with
width $m_{l}=m^{L-l+1},\ \forall l>1$ such that

$\displaystyle\|u_{\boldsymbol{\theta}}-f\|_{H^{2}(\mathbb{P})}$
$\displaystyle\leq\frac{3L\|f\|_{\mathcal{W}^{L}(\Omega)}}{\sqrt{m}},$ (13)
$\displaystyle\|u_{\boldsymbol{\theta}}-f\|_{L^{2}(\mathbb{Q})}$
$\displaystyle\leq\frac{3C_{\Omega}L\|f\|_{\mathcal{W}^{L}(\Omega)}}{\sqrt{m}},$

where $C_{\Omega}$ is a universal constant only depends on the domain
$\Omega$, and $H^{2}=W^{2,2}$ is the Sobolev space, and
$\|\boldsymbol{\theta}\|_{\mathcal{P}}\leq\|\boldsymbol{W}^{L}\|_{1,\infty}\leq\|f\|_{\mathcal{W}^{L}(\Omega)}$,
where $\|\cdot\|_{\mathcal{P}}$ is the path norm defined as

$\|\boldsymbol{\theta}\|_{\mathcal{P}}=\sum_{i_{L}}\cdots\sum_{i_{0}}|\boldsymbol{W}^{L}_{i_{L}}\cdots\boldsymbol{W}^{1}_{i_{1}i_{0}}|.$
(14)

The path norm is one type of complexity measure of neural networks correlated
to generalization [21]. The above theorem shows that the neural networks can
approximate any target function in the generalized Barron space and its
derivatives well with complexities controlled by the Barron norm of the target
functions, which shows the efficiency of network approximation. Note that
since $\Omega$ is fixed in this paper, the constant $C_{\Omega}$ is actually
universal. This is utilized in the proof of our prior generalization bound in
Theorem 3.1.

Actually, the theoretical result that neural networks can approximate a
function and its derivative is not new [22, 23]. Their proof idea is two-step.
First, they show that polynomials are dense in $C^{k}(\overline{\Omega})$.
Second, they can approximate polynomials in $C^{k}$-norm, using Taylor’s
expansion. In our paper, we adopt a functional analysis approach to adapt to
the Barron space setting, i.e., to show additionally that such networks have
low complexity measured by Barron norm, which is indispensable to our prior
bound.

##  3 Theory

In this section, we introduce our main generalization results, including a
prior bound based on the Barron space and a posterior bound based on the
Rademacher complexity. For both of them, we use the following assumption
adopted from a closely related previous study [13].

######  Assumption 3.1.

(Symmetry and boundedness of $\mathcal{L}$). Throughout the analysis in this
paper, we assume the differential operator $\mathcal{L}$ in the PDE satisfies
the following conditions. The operator $\mathcal{L}$ is a linear second-order
differential operator in a non-divergence form, i.e.,
$(\mathcal{L}u^{*})(\boldsymbol{x})=\sum_{\alpha=1,\beta=1}^{d}\boldsymbol{A}_{\alpha\beta}(\boldsymbol{x})u^{*}_{x_{\alpha}x_{\beta}}(\boldsymbol{x})+\sum_{\alpha=1}^{d}\boldsymbol{b}_{\alpha}(\boldsymbol{x})u^{*}_{x_{\alpha}}(\boldsymbol{x})+c(\boldsymbol{x})u^{*}(\boldsymbol{x})$,
where all
$\boldsymbol{A}_{\alpha\beta},\boldsymbol{b}_{\alpha},c:\Omega\rightarrow\mathbb{R}$
are given coefficient functions and $u^{*}_{x_{\alpha}}$ are the first-order
partial derivatives of the function $u^{*}$ with respect to its $\alpha$-th
argument (the variable $x_{\alpha}$) and $u^{*}_{x_{\alpha}x_{\beta}}$ are the
second-order partial derivatives of the function $u^{*}$ with respect to its
$\alpha$-th and $\beta$-th arguments (the variables $x_{\alpha}$ and
$x_{\beta}$). Furthermore, there exists constant $K>0$ such that for all
$\boldsymbol{x}\in\Omega=[-1,1]^{d}$, and $\alpha,\beta\in[d]$, we have
$A_{\alpha\beta}=A_{\beta\alpha}$ and
$A_{\alpha\beta}(\boldsymbol{x}),b_{\alpha}(\boldsymbol{x}),c(\boldsymbol{x})$
are all $K$-Lipschitz, and their absolute values are not larger than $K$.

Because multiplying the network functions by the coefficients
$\boldsymbol{A},\boldsymbol{b},c$ and differentiation on them influence their
complexities, the universal bound on the coefficients
$\boldsymbol{A},\boldsymbol{b},c$ and the restriction to second order PDEs are
required for our estimation on the Rademacher complexity of the hypothesis
class of PINNs.

###  3.1 A Prior Generalization Bound (Theorem 3.1)

In this subsection, we introduce our prior bound based on the Barron space.

######  Theorem 3.1.

(A prior generalization bound on PINN). Let the Assumption given in 3.1 holds,
then for any $\delta\in(0,1)$ and the depth $L$, suppose that the true
solution $u^{*}(\boldsymbol{x})$ lies in the tree-like function space
$\mathcal{W}^{L}(\Omega)$, and set $\lambda={3(2KC_{\Omega}+1)L^{2}}/{m}$. Let
$\boldsymbol{\theta}^{*}=\arg\min_{\boldsymbol{\theta}}R_{S}(\boldsymbol{\theta})+\lambda\|\boldsymbol{\theta}\|^{2}_{\mathcal{P}}$.
Then, with probability at least $1-\delta$ over the choice of random samples
$S=\left\\{\boldsymbol{x}_{i}\right\\}_{i=1}^{n_{b}+n_{r}}\subset\overline{\Omega}$
with $n_{b}$ boundary points and $n_{r}$ residual points, we obtain the
following generalization bound

$\displaystyle R_{{D}\cap\partial\Omega}(\boldsymbol{\theta}^{*})$
$\displaystyle\leq
R_{S\cap\partial\Omega}(\boldsymbol{\theta}^{*})+8\|u^{*}\|_{\mathcal{W}^{L}(\Omega)}\frac{C(h)\log
n_{b}}{{\sqrt{n_{b}}}}+2\sqrt{\frac{\log(2/\delta)}{n_{b}}},$ (15)
$\displaystyle R_{{D}\cap\Omega}(\boldsymbol{\theta}^{*})\leq$ $\displaystyle
R_{S\cap\Omega}(\boldsymbol{\theta}^{*})+8\left(\|u^{*}\|_{\mathcal{W}^{L}(\Omega)}\right)^{3}\frac{C(h,K)\log
n_{r}}{{\sqrt{n_{r}}}}+2\sqrt{\frac{\log(2/\delta)}{n_{r}}},$

where $C(h)$ and $C(h,K)$ are universal constants depending only on $h$ and
$h,K$, respectively.

In Theorem 3.1, the generalization bounds on the right-hand side of equation
(15) (for both the boundary and the residual points) contain three terms,
where the first term is the empirical training loss, the second term is the
complexity of the model (original network for boundary loss and differentiated
network for residual loss), and the third term is the statistical term.
Moreover, this theorem shows that under certain regularization of the path
norm, the generalization errors of PINNs are controlled by the Barron norm of
the target function $u^{*}$. If the target function $u^{*}$ is more complex
(simpler), i.e., it has larger (smaller) Barron norm, the generalization error
will be larger (smaller). This reflects a data-dependent bound in which neural
networks control their complexity based on those of target functions. The
above advantages are also summarized in the appendix to justify our choice of
Barron space.

###  3.2 A Posterior Generalization Bound (Theorem 3.2)

We now provide a posterior generalization bound based on the optimized network
parameters (which are obtained after optimization). We begin by defining the
Rademacher complexity [24], which is one of key notions in statistical
learning theory.

######  Definition 3.1.

(Rademacher Complexity). Let
$S=\left\\{x_{i}\right\\}_{i=1}^{n}\subset\overline{\Omega}$ be a dataset
containing $n$ samples. The Rademacher complexity of a function class
$\mathcal{F}$ on $S$ is defined as
$\text{Rad}(\mathcal{F};S)=\mathbb{E}_{\epsilon}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}f(x_{i})\right]$,
where $\epsilon_{1},\dots,\epsilon_{n}$ are independent and identically
distributed (i.i.d.) random variables taking values uniformly in $\\{-1,1\\}$.

Intuitively, Rademacher complexity measures the richness of function class
$\mathcal{F}$ by studying its ability to fit random labels of $x_{i}$
generated by $\epsilon_{i}$. Since simpler function classes tend to generalize
better on unseen testing data, we will investigate the Rademacher complexities
of PINNs and XPINNs, which begins with a key lemma on that of neural networks.

######  Lemma 3.1.

[25] For every $L$, and every set of $n$ points $S\subset\overline{\Omega}$,
the hypothesis class $\mathcal{NN}^{L}_{M,N}$ given by the neural networks

$\mathcal{NN}^{L}_{M,N}:=\left\\{\boldsymbol{x}\mapsto
W_{L}\sigma(W_{L-1}\sigma(\cdots\sigma(W_{1}\boldsymbol{x})\cdots)\ |\
\|W_{l}\|_{2}\leq M(l),\frac{\|W_{l}\|_{2,1}}{\|W_{l}\|_{2}}\leq
N(l)\right\\},$ (16)

satisfies the Rademacher complexity bound

$\text{Rad}(\mathcal{NN}^{L}_{M,N};S)\leq\frac{4}{n\sqrt{n}}+\frac{18\sqrt{d\log(2h^{2})}\log
n}{\sqrt{n}}\prod_{l=1}^{L}M(l)\Big{(}\sum_{l=1}^{L}N(l)^{2/3}\Big{)}^{3/2},$
(17)

where $h$ is the maximal width of the neural network, i.e.,
$h=\max(m_{L},\cdots,m_{0})$.

This lemma controls Rademacher complexities of neural networks by the product
of the spectral norms of the network parameter matrices at each layer, i.e.,
$M(l)$. The complexity depends on the network depth via the $N(l)$ term as we
always have $N(l)\geq 1$. We extend this result to the differentiated PINNs in
the following lemma.

######  Lemma 3.2.

(Rademacher Complexity of Differentiated Networks). For every $L$, and every
set of $n$ points $S\subset\overline{\Omega}$, the hypothesis class

$\mathcal{PINN}^{L}_{M,N}=\left\\{\boldsymbol{x}\mapsto\mathcal{L}u(\boldsymbol{x})\
|\ u(\boldsymbol{x})\in\mathcal{NN}^{L}_{M,N}\right\\},$ (18)

satisfies the Rademacher complexity bound

$\displaystyle\text{Rad}(\mathcal{PINN}^{L}_{M,N};S)\leq\frac{8K+4d(L-1)K}{n_{r}\sqrt{n_{r}}}+\frac{18K\sqrt{d\log(2h^{2})}\log
n_{r}}{\sqrt{n_{r}}}\cdot$ (19)
$\displaystyle\quad\prod_{l=1}^{L}M(l)\left(\sum_{l=1}^{L}N(l)^{2/3}\right)^{3/2}\left[1+\sqrt{2}L\prod_{l=1}^{L}M(l)+\sqrt{2}d(L^{2}-1)\left(\prod_{l=1}^{L}M(l)\right)^{2}\right].$

This lemma shows a similar Rademacher complexity bound for PINNs, where the
main differences are due to the first order and second order differentiation
in PINNs, respectively. Using these lemmas, we derive the following posterior
generalization bound.

######  Theorem 3.2.

(A posterior generalization bound on PINN). Let Assumptions 3.1 hold, for any
$\delta\in(0,1)$ and the depth $L$, let the (not regularized) empirical loss
function be
$\boldsymbol{\theta}_{S}=\arg\min_{\boldsymbol{\theta}}R_{S}(\boldsymbol{\theta}).$
Then, with probability at least $1-\delta$ over the choice of random samples
$S=\left\\{\boldsymbol{x}_{i}\right\\}_{i=1}^{n_{b}+n_{r}}\subset\overline{\Omega}$
with $n_{b}$ boundary points and $n_{r}$ residual points, we have the
following generalization bound

$\displaystyle R_{D\cap\partial\Omega}(\boldsymbol{\theta}_{S})$
$\displaystyle\leq
R_{S\cap\partial\Omega}(\boldsymbol{\theta}_{S})+\frac{32}{n_{b}\sqrt{n_{b}}}+\frac{144\sqrt{d\log(2h^{2})}\log
n_{b}}{\sqrt{n_{b}}}\cdot$ (20)
$\displaystyle\quad\prod_{l=1}^{L}M(l)\Big{(}\sum_{l=1}^{L}N(l)^{2/3}\Big{)}^{3/2}+2\sqrt{\frac{\log(2/\delta(M,N))}{2n_{b}}},$
$\displaystyle\quad
R_{D\cap\Omega}(\boldsymbol{\theta}_{S})-R_{S\cap\Omega}(\boldsymbol{\theta}_{S})$
(21)
$\displaystyle\leq\frac{64K+32d(L-1)K}{n_{r}\sqrt{n_{r}}}+2\sqrt{\frac{\log(2/\delta(M,N))}{2n_{r}}}+\frac{144K\sqrt{d\log(2h^{2})}\log
n_{r}}{\sqrt{n_{r}}}\cdot$
$\displaystyle\quad\prod_{l=1}^{L}M(l)\left(\sum_{l=1}^{L}N(l)^{2/3}\right)^{3/2}\left[1+\sqrt{2}L\prod_{l=1}^{L}M(l)+\sqrt{2}d(L^{2}-1)\left(\prod_{l=1}^{L}M(l)\right)^{2}\right].$

where $M(l)=\lceil\|\boldsymbol{W}^{l}\|_{2}\rceil$, and
$N(l)=\lceil\|\boldsymbol{W}^{l}\|_{2,1}/\|\boldsymbol{W}^{l}\|_{2}\rceil$, in
which $\lceil a\rceil$ of $a\in\mathbb{R}$ is the smallest integer that is
greater than or equal to $a$, and

$\delta(M,N)={\delta}/\left[{\prod_{l=1}^{L}M(l)(M(l)+1)N(l)(N(l)+1)}\right].$
(22)

The generalization bounds are mainly controlled by the complexity of networks
measured by the spectral norm $\left\\{M(l)\right\\}_{l=1}^{L}$ and the (2,1)
norm $\left\\{N(l)\right\\}_{l=1}^{L}$, as well as the last statistical term.
When compared to the prior bound in Theorem 3.1, the posterior bound in
Theorem 3.2 is easier to compute numerically, because it only involves terms
related to neural network parameters. Despite the difference, Theorems 3.1 and
3.2 are related. Concretely, Theorem 3.1 shows that if the target function
$u^{*}$ is more complex (simpler), i.e., has larger (smaller) Barron norm, the
generalization error is larger (smaller), which implies that a complex
(simple) neural network has been learnt to fit the target since complex
(simple) network generalizes worse (better). On the other hand, the complexity
of neural networks can also be reflected by the quantities
$\left\\{M(l),N(l)\right\\}_{l=1}^{L}$ in Theorem 3.2 since they are directly
linked to Rademacher complexity in Lemmas 3.1 and 3.2. Hence, we can expect
that PINNs learning more complex target functions have larger
$\left\\{M(l),N(l)\right\\}_{l=1}^{L}$ quantities. And if a trained PINN has
large $\left\\{M(l),N(l)\right\\}_{l=1}^{L}$ quantities, which signifies
higher complexity, the target function fitted should also be complicated
(i.e., it has a larger Barron norm) due to the implicit regularization of
network training, i.e., stochastic gradient training finds out complex
(simple) solution given a complex (simple) target.

###  3.3 Posterior $L_{2}$ Error Generalization Bound

In this subsection, we will bridge the gap between the boundary + residual
generalization error with the $L_{2}$ generalization error. In particular, we
adopt the following assumption widely used in numerical PDE methods [26, 27]
and PINN theory [16, 15, 17, 28]. Intuitively, the assumption states that
minimization of the boundary and residual errors contributes to the
minimization of $L_{2}$ error.

######  Assumption 3.2.

Assume that the PDE satisfies the following norm constraint:

$\displaystyle
C_{1}\|u\|_{L_{2}(\Omega)}\leq\|\mathcal{L}u\|_{L_{2}(\Omega)}+\|u\|_{L_{2}(\partial\Omega)},\qquad\forall
u\in\mathcal{NN}_{L}^{M},\forall L,M,$ (23)

where the positive constant $C_{1}$ does not depend on $u$ but on the domain
and the coefficients of the operators $\mathcal{L}$.

This is Assumption 2.5 in [29] and Assumption 2.1 in [15], which are two
theory papers on PINN. In particular, the function space $(X,\|\cdot\|_{X})$
in [29] becomes the function space of all neural network functions with the
$L_{2}$ norm, $(Y,\|\cdot\|_{Y})$ becomes the $L_{2}$ function space
$(\Omega,\|\cdot\|_{L_{2}(\Omega)})$, and $(Z,\|\cdot\|_{Z})$ becomes the
$L_{2}$ function space $(\partial\Omega,\|\cdot\|_{L_{2}(\partial\Omega)})$.

Intuitively, Assumption 3.2 specifies the well-posedness of the PDE problem,
and justifies the motivation of minimizing the boundary and residual losses.

To show that Assumption 3.2 is realistic, various papers have proved that
Assumption 3.2 holds for various PDEs. [15] proves that our Assumption 3.2
holds for Poisson equation, heat equation, wave equation and Stokes equation.
[28] proves that Assumption 3.2 holds for Kolmogorov equations that include
the heat equation and Black-Scholes equation as special cases. [17] proves
Assumption 3.2 holds for incompressible Navier-Stokes equations. Lastly, [14]
proves that optimization of the Deep Ritz Method objective contributes to
$L_{2}$ error, on Poisson equation and static Schrödinger equation. Therefore,
we refer the readers to [14, 15, 28, 17] for concrete PDE examples satisfying
the assumption.

From the practice side [6], the motivation of the boundary and residual losses
is to embed the physical law governed by the PDE and the data into the neural
networks. Empirically, as long as the PDE problem is well-posed (i.e., it
satisfies Assumption 3.2) and that data are sufficient (which means the
corresponding bound will be small), then PINNs can solve the PDE problem with
small error. Assumption 3.2 is basically the foundation of PINNs.

The following theorem bridges the boundary and residual losses with the L2
error.

######  Theorem 3.3.

Let Assumption 3.2 holds, then for all neural networks parameterized by
$\theta$, we can connect their $L_{2}$ generalization errors with the boundary
+ residual losses as follows:

$\|u_{\theta}-u\|_{L_{2}(\Omega)}\leq\sqrt{2}C_{1}^{-1}\left(R_{D\cap\Omega}(\theta)+R_{D\cap\partial\Omega}(\theta)\right)^{1/2}.$
(24)

Here, we can use both the prior bound (Theorem 3.1) and the posterior bound
(Theorem 3.2) to evaluate the $L_{2}$ distance between solution of PINN or
XPINN and the true solution of PDEs.

###  3.4 Comparing XPINN and PINN by Theorem 3.1

In this subsection, we compare PINNs with XPINNs by using the generalization
bound in Theorem 3.1. We focus on comparing PINN with XPINN on the residual
loss, i.e. $R_{D\cap\Omega}$ and $R_{S\cap\Omega}$, because it is more
representative of the differentiated nets in PINNs. The case for boundary loss
is similar but simpler, which is included in the Appendix. Specifically, the
comparison is performed by computing their respective theoretical bounds. In
particular, the generalization performance of PINN depends on the upper bound
in Theorem 3.1, which is:
$R_{S\cap\Omega}(\boldsymbol{\theta}^{*})+8\left(\|u^{*}\|_{\mathcal{W}^{L}(\Omega)}\right)^{3}{C(h,K)\log
n_{r}}/{{\sqrt{n_{r}}}}+2\sqrt{{2(\log(2/\delta))}/{n_{r}}}$, where $n_{r}$ is
the number of residual training points.

For XPINN’s generalization, we can apply Theorem 3.1 to each of the subdomains
in the XPINN. Specifically, for the $i$-th sub-net in the $i$-th subdomain of
XPINN, i.e., the $\Omega_{i},i\in\left\\{1,2,...,N_{D}\right\\}$, its
generalization performance is upper bounded by
$R_{S\cap\Omega_{i}}(\boldsymbol{\theta}^{*})+8\left(\|u^{*}\|_{\mathcal{W}^{L}(\Omega_{i})}\right)^{3}{C(h,K)\log
n_{r,i}}/{{\sqrt{n_{r,i}}}}+2\sqrt{{2(\log(2/\delta))}/{n_{r,i}}}$, where
$n_{r,i}$ is the number of training boundary points in the $i$-th subdomain.

Hence, since the $i$-th subdomain has $n_{r,i}$ training boundary points and
is in charge of the prediction of $\frac{n_{r,i}}{n_{r}}$ proportion of
testing data, we weight-averaged their generalization errors to get the
generalization error of XPINN:
$\sum_{i=1}^{N_{D}}\frac{n_{r,i}}{n_{r}}(R_{S\cap\Omega_{i}}(\boldsymbol{\theta}^{*})+8\left(\|u^{*}\|_{\mathcal{W}^{L}(\Omega_{i})}\right)^{3}{C(h,K)\log
n_{r,i}}/{{\sqrt{n_{r,i}}}}+2\sqrt{{2(\log(2/\delta))}/{n_{r,i}}})$. If we
omit the last term and assume that their empirical losses are similar, i.e.,
$R_{S\cap\Omega}\approx\sum_{i=1}^{N_{D}}({n_{r,i}}/{n_{r}})R_{S\cap\Omega_{i}}$,
and
$\sqrt{{2\log(2/\delta)}/{n_{r,i}}}\ll\|u^{*}\|_{\mathcal{W}^{L}(\Omega)},\|u^{*}\|_{\mathcal{W}^{L}(\Omega_{i})}$,
then comparing the generalization ability of PINN and XPINN reduces to the
following:

${\|u^{*}\|^{3}_{\mathcal{W}^{L}(\Omega)}}{\text{(PINN)}}\
\qquad\text{versus}\qquad{\sum_{i=1}^{N_{D}}\frac{\log
n_{r,i}\sqrt{n_{r,i}}}{\log
n_{r}\sqrt{n_{r}}}\|u^{*}\|^{3}_{\mathcal{W}^{L}(\Omega_{i})}}{\text{(XPINN)}},$
(25)

where the model having smaller corresponding quantity is more generalizable.

In the next section, we will present three analytic examples and adopt the
above comparison method to illustrate the circumstances under which XPINN is
better or worse than PINN.

###  3.5 Comparing XPINN and PINN by Theorem 3.2

The comparison using Theorem 3.2 is also done via computing their respective
theoretical bounds. The residual loss is considered in the main text while the
case for boundary loss is included in the supplementary material. Concretely,
we denote the upper bound of PINN testing loss as $B_{\text{PINN}}$ and those
of the sub-net $i$ in XPINN as $B_{i,\text{XPINN}}$,
$i\in\left\\{1,2,...,N_{D}\right\\}$, which are provided by the right sides of
Theorem 3.2, i.e., the bounds are

$\displaystyle\quad
B_{\text{PINN}}=\frac{64K+32d(L-1)K}{n_{r}\sqrt{n_{r}}}+2\sqrt{\frac{\log(2/\delta(M,N))}{2n_{r}}}+\frac{144K\sqrt{d\log(2h^{2})}\log
n_{r}}{\sqrt{n_{r}}}\cdot$ (26)
$\displaystyle\quad\prod_{l=1}^{L}M(l)\left(\sum_{l=1}^{L}N(l)^{2/3}\right)^{3/2}\left[1+\sqrt{2}L\prod_{l=1}^{L}M(l)+\sqrt{2}d(L^{2}-1)\left(\prod_{l=1}^{L}M(l)\right)^{2}\right],$

and

$\displaystyle\quad
B_{i,\text{XPINN}}=\frac{64K+32d(L-1)K}{n_{r,i}\sqrt{n_{r,i}}}+2\sqrt{\frac{\log(\frac{2}{\delta(M_{i},N_{i})})}{2n_{r,i}}}+\frac{144K\sqrt{d\log(2h^{2})}\log
n_{r,i}}{\sqrt{n_{r,i}}}\cdot$ (27)
$\displaystyle\quad\prod_{l=1}^{L}M_{i}(l)\left(\sum_{l=1}^{L}N_{i}(l)^{2/3}\right)^{3/2}\left[1+\sqrt{2}L\prod_{l=1}^{L}M_{i}(l)+\sqrt{2}d(L^{2}-1)\left(\prod_{l=1}^{L}M_{i}(l)\right)^{2}\right].$

Specifically, we assume that all sub-PINNs as well as the PINN model use
neural networks with depth $L$ and width $h$. In the bound of PINN, $n_{r}$ is
the total number of residual training samples.
$M(l)=\lceil\|\boldsymbol{W}^{l}\|_{2}\rceil$, and
$N(l)=\lceil\|\boldsymbol{W}^{l}\|_{2,1}/\|\boldsymbol{W}^{l}\|_{2}\rceil$,
where $\boldsymbol{W}^{l}$ is the $l$-th layer parameter matrix in the PINN
model. Moreover, in the bound of XPINN, $n_{r,i}$ is the number of residual
training samples in subdomain $i$.
$M_{i}(l)=\lceil\|\boldsymbol{W}^{l}_{i}\|_{2}\rceil$, and
$N_{i}(l)=\lceil\|\boldsymbol{W}^{l}_{i}\|_{2,1}/\|\boldsymbol{W}^{l}_{i}\|_{2}\rceil$,
where $\boldsymbol{W}^{l}_{i}$ is the $l$-th layer parameter matrix of the
$i$-th subnet in the XPINN model. Because the $i$-th sub-net in XPINN is in
charge of the prediction of $\frac{n_{r,i}}{n_{r}}$ proportion of testing
data, we weight-averaged their bounds to get that of XPINN, i.e.,
$B_{\text{XPINN}}=\sum_{i=1}^{N_{D}}({n_{r,i}}/{n_{r}})B_{i,\text{XPINN}}$
where $B_{\text{XPINN}}$ is the bound for XPINN. Thus, we only need to compare
$B_{\text{PINN}}$ with $B_{\text{XPINN}}$, where the model having smaller
corresponding quantity is more generalizable. These quantities can be directly
measured and calculated from the trained deep nets, which allows for easy
numerical validation. Thus, we validate this comparison method in
computational experiments.

###  3.6 Comparing XPINN and PINN by Theorem 3.3

In this subsection, we compare PINNs with XPINNs on the $L_{2}$ error by using
the generalization bound in Theorem 3.3. Specifically, the comparison is
performed by computing their respective theoretical bounds. In particular, we
have already shown how to compare PINN and XPINN based on the boundary and
residual losses. Denote the boundary and residual bounds of PINN (XPINN) as
$B_{\text{PINN}}^{\text{boundary}}$ ($B_{\text{XPINN}}^{\text{boundary}}$) and
$B_{\text{PINN}}^{\text{residual}}$ ($B_{\text{XPINN}}^{\text{residual}}$).
Then, we only need to compare the following quantites:

$\left(B_{\text{PINN}}^{\text{boundary}}+B_{\text{PINN}}^{\text{residual}}\right)^{1/2}(\text{PINN})\quad\text{versus}\quad\left(B_{\text{XPINN}}^{\text{boundary}}+B_{\text{XPINN}}^{\text{residual}}\right)^{1/2}(\text{XPINN}).$
(28)

##  4 Analytical Examples Based on Theorem 3.1

In this section, we provide analytical examples to further analyze the prior
generalization bound in Theorem 3.1. The examples ensure analytical
expressions of the Barron norm, which results in precise calculation of prior
bounds. Specifically, we show in what cases XPINNs are better than, similar
to, and worse than PINNs in order to demonstrate the tradeoff in XPINN
generalization.

Figure 1: Left: Decomposition causes simpler target in each part (red), but
less training data leads to overfitting. Right: Tradeoff between the two
factors.

###  4.1 Case where XPINN Outperforms PINN

Let us consider the target function $u^{*}(x,y)=2\sin x+\sin y$, on the broken
line
$\Omega=\Omega_{1}\cup\Omega_{2}=[0,1]\times\left\\{0\right\\}\cup\left\\{0\right\\}\times[0,1]$.
Obviously, we have $\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}=3$, recall that
$\mathcal{W}^{2}(\Omega)$ is the natural function space of two-layer sine
networks on $\Omega$. Interestingly, if we restrict $u^{*}(x,y)$ to
$\Omega_{1}$, we have $u^{*}(x,y)=2\sin x$ on $\Omega_{1}$, with a Barron norm
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}=2$. Similarly, if we restrict our
observation to $\Omega_{2}$, we have
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}=1$. Since the lines $\Omega_{1}$ and
$\Omega_{2}$ have the same length, it is natural to assume that the numbers of
train residual data points on $\Omega_{1}$ and $\Omega_{2}$ are the same, i.e.
$n_{r,1}=n_{r,2}=n_{r}/2$ in equation (25). In addition, we also assume that
there is sufficient training data, with $\log n_{r}\gg\log 2$. We compare a
PINN on $\Omega$ and an XPINN with two sub-nets on $\Omega_{1}$ and
$\Omega_{2}$, respectively. Applying Theorem 3.1 and our discussion in section
3.1,
$27=\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}^{3}\gtrsim(\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}^{3}+\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{2})}^{3})/\sqrt{2}={(8+1)}/{\sqrt{2}}={9}/{\sqrt{2}}$.
Thus, XPINN generalizes better than PINN.

The underlying reason is as follows: $u^{*}$ remains simple on $\Omega_{2}$,
while being complex on $\Omega_{1}$. Since XPINN optimizes several sub-nets at
different subdomains, XPINN may learn a simple network on $\Omega_{2}$ where
the solution is simpler, (i.e., very small
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{2})}=1$). Also, it learns a complex
network in other subdomains, where the solution is relatively complicated,
(i.e., relatively large $\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}=2$).
Therefore, compared with PINN, which learns a very complex network on the
whole domain $\Omega$ (extremely large
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}=3$), XPINN tends to have lower overall
complexity, because it is complex on only part of the domain and remains
simple on the rest of the domain, which leads to better generalization.

###  4.2 Case where XPINN is Worse Than PINN

Let us consider the same target function $u^{*}(x,y)=2\sin x+\frac{1}{2}\sin
y$, on a different broken line $\Omega=\Omega_{1}\cup\Omega_{2}$, where
$\Omega_{1}=[0,1]\times\left\\{0\right\\}$ and
$\Omega_{2}=\left\\{(x,y)|y=x,x\in[0,{\sqrt{2}}/{2}]\right\\}$. Obviously, we
have $\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}=2.5$. Moreover, if we restrict
$u^{*}(x,y)$ to $\Omega_{1}$, we have $u^{*}(x,y)=2\sin x$ on $\Omega_{1}$,
with a Barron norm $\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}=2$. However, even
if we restrict our observation to $\Omega_{2}$, we still have
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{2})}=2.5$. Since the lines $\Omega_{1}$
and $\Omega_{2}$ have the same length, it is natural to assume that the
numbers of train residual data points on $\Omega_{1}$ and $\Omega_{2}$ are the
same, i.e. $n_{r,1}=n_{r,2}=n_{r}/2$ in equation (25). We compare a PINN on
$\Omega$ and an XPINN with two sub-nets on $\Omega_{1}$ and $\Omega_{3}$
respectively. Applying Theorem 3.1 and following our discussion in section
3.1, we have
$15.625=\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}^{3}\lesssim(\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}^{3}+\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{2})}^{3})/\sqrt{2}={(2^{3}+2.5^{3})}/{\sqrt{2}}={23.625}/{\sqrt{2}}\approx
16.70$. Thus, in this example XPINN is worse than PINN.

Although XPINN decreases the target function complexity via decomposition, at
least on $\Omega_{1}$ it decreases to 2 from 2.5, it cannot complement the
overfitting of less available training data on generalization, which is
reflected in the $\sqrt{1/n_{r}}$ term in the bound, where $n_{r}$ is the
number of residual training samples, i.e. the complexity grows with less data.
Unfortunately, in this example, the more complexity brought by overfitting due
to less data exceeds the benefit of simpler target function parts after
decomposition. Hence, XPINN performs worse than PINN.

###  4.3 Illustration of a Tradeoff in XPINN generalization

In this section, we summarize the above two examples and derive a tradeoff in
XPINN generalization, which is illustrated in Figure 1. There are two factors
that counter-balance each other to affect XPINN generalization, namely the
simplicity of decomposed target function within each subdomain thanks to
domain decomposition, and the complexity and inclination to overfit due to
less available training data, which counter-balance each other as follows.
When the former is more dominant, XPINN outperforms PINN, as in our example in
section 4.1. Otherwise, PINN outperforms XPINN, as in our example in section
4.2. When the two factors reach a balance, XPINN and PINN perform similarly.

To make the idea clearer, we consider another analytical example. Let us
consider the target function $u^{*}(x,y)=2\sin x+q\sin y$, on the broken line
$\Omega=\Omega_{1}\cup\Omega_{2}$, where
$\Omega_{1}=[0,1]\times\left\\{0\right\\}$ and
$\Omega_{2}=\left\\{(x,y)|y=x,x\in[0,{\sqrt{2}}/{2}]\right\\}$, where
$q\in\mathbb{R}^{+}$ is a fixed constant to be decided. Obviously, we have
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}=2+q$. Further, if we restrict
$u^{*}(x,y)$ to $\Omega_{1}$, we have $u^{*}(x,y)=2\sin x$ on $\Omega_{1}$,
with a Barron norm $\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}=2$. However, even
if we restrict our observation to $\Omega_{2}$, we still have
$\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{2})}=2+q$. Because the lines $\Omega_{1}$
and $\Omega_{2}$ have the same length, it is natural to assume that the
numbers of train residual data points on $\Omega_{1}$ and $\Omega_{2}$ are the
same, i.e. $n_{r,1}=n_{r,2}=n_{r}/2$ in equation (25). We compare a PINN on
$\Omega$ and an XPINN with two sub-nets on $\Omega_{1}$ and $\Omega_{2}$,
respectively. Applying Theorem 3.1 and following our discussion in Section
3.1, we need to compare the following quantities to determine when XPINN
outperforms PINN: $(2+q)^{3}=\|u^{*}\|_{\mathcal{W}^{2}(\Omega)}^{3}\
\text{versus}\
(\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{1})}^{3}+\|u^{*}\|_{\mathcal{W}^{2}(\Omega_{2})}^{3})/\sqrt{2}={(8+(2+q)^{3})}/{\sqrt{2}}$.
When $(2+q)^{3}<8/(\sqrt{2}-1)\approx 19.31$, i.e., $q<0.683$, PINN is better
due to its more obvious effect with less data. When $q>0.683$, XPINN performs
better due to the more obvious effect by decomposing complexity into
simplicity.

In summary, inspired by the two analytical examples, we have shown in this
section that there exists a tradeoff in XPINNs, which results in their
different performance when compared to PINNs. We further demonstrate these
phenomena in computational experiments for various PDEs in the next section.

##  5 Computational Experiments

Figure 2: Data visualization for the KdV experiment. Left: exact solution of
the KdV equation. Right: Training points.

Figure 3: Error visualization for the KdV experiment.

###  5.1 KdV Equation

####  5.1.1 Setup

In this experiment, we consider a one-dimensional KdV equation given by
$u_{t}+uu_{x}=0.0025u_{xxx},x\in[-1,1],t\in[0,1]$, with the boundary condition
of $u(x,0)=\cos(\pi x),x\in[-1,1]$ and the periodic initial condition. The
true solution is visualized in Figure 2 left. The entire dataset for this PDE
is provided by the paper of PINN [6] and CPINN [30].

Following the KdV equation experiment in CPINN [30], the training dataset for
PINN contains 18000 residual points and 914 boundary points. The testing
dataset for PINN contains 102400 points uniformly distributed within the
domain. The backbone model for PINN is a 10-layer neural network with 20
hidden units activated by sine as in our theory. Adam [31] optimizer with 1e-3
learning rate is used for optimization. No regularization is used.

Moreover, from the solution of the KdV equation in Figure 2 left, we observe
that it is complex, and that it fluctuates when $x\geq 0$ corresponding to the
right part. In contrast, the left part is smoother and simpler. Hence, to
reflect our discussion on the prior bound and analytical examples, for XPINN
we partition the whole domain into two subdomains, including (1) the right
domain, corresponding to sub-net named XPINN-R as $x>-0.74$ and (2) the left
domain, corresponding to sub-net named XPINN-L as $x\leq-0.74$. Then, the
target function in subdomain 1 is complicated and it fluctuates significantly,
while that in subdomain 2 is simpler and smoother.

For XPINN, it is given 646 and 268 boundary points, 14000 and 4000 residual
points, in subdomain 1 and 2, respectively. The number of interface points is
10000. The backbone models for XPINN are two 10-layer neural networks with 20
hidden units activated by sine as in our theory. Two Adam [31] optimizers with
1e-3 learning rates are used for optimizations. No regularization is used.

For fair comparison, we keep the same training procedure, i.e., training
epochs, learning rate, model structure, and weight decay, etc. We train each
model for 5000 epochs, and the results reported in the table are those at the
5000-th epochs. For both models, we use unity weight for the residual loss and
zero weight for the residual continuity loss, and also use unity weight for
the boundary loss and the boundary interface loss. For reproducibility, we run
each model for 5 times using fixed random seed 0, 1, 2, 3, 4.

Lastly, we need to pay attention to the choice of $\delta$ in the posterior
bound. Because the bound in Theorem 3.2 holds with probability $1-\delta$, in
the computations we take $\delta=0.1$ so that the bound holds with probability
0.9. At the same time, for XPINN, we need to take $\delta=0.05$ because there
are two sub-nets for a union bound.

####  5.1.2 Results

Table 1: Computational results for KdV equation. | Train Loss |  Relative $L_{2}$ error | Complexity | Bound
---|---|---|---|---
PINN |  3.597e-3$\pm$7.194e-4 |  6.899e-1$\pm$8.015e-3 | 100.00% | 100.00%
XPINN-R | 5.619e-3$\pm$5.056e-3 | 6.955e-1$\pm$9.905e-3 | 101.31% | 121.08%
XPINN-L | 28.50%

We present the experimental results including train loss, and test relative L2
error, as well as the calculated theoretical generalization bound in XPINN and
PINN models in Table 1. Moreover, we also provide the product of norms of the
neural network parameter matrices, which is directly linked to the Rademacher
complexity of neural networks and thus their complexities (see Lemma 3.1). We
compute this quantity to provide an intuitive observation on the complexities
of the optimized neural networks.

In addition, the products of norms of the weight matrices are presented in the
columns “Complexity” where that of PINN is denoted 100% for clear comparison.
The complexity of PINN is similar to that of XPINN-R, while XPINN-R is more
complex than XPINN-L. XPINN-R corresponds to the domain $x>-0.74$, where the
solution is more complicated and XPINN-L corresponds to the simpler domain,
validating the implicit regularization of gradient descent, which learns
simple (complex) function with simple (complex) neural networks. Although
XPINN-R is fitting a less complex function, the available data in its
subdomain 1 is less than those of PINN, so overfitting causes XPINN-R to be as
complex as the PINN fitting the entire function.

Furthermore, the “Bound” columns in Table 1 are the theoretical generalization
bounds for PINN and XPINN, where that of PINN is denoted as 100% for clarity.
The theoretical generalization bound of XPINN is 121.08%, which is slightly
larger than PINN and consistent with their testing performances, i.e., the
error of PINN (6.899e-1) is slightly better than that of XPINN (6.955e-1),
justifying the effectiveness of our generalization bound. In Figure 3, we
visualize the errors of PINN and XPINN. PINN and XPINN have similar error
distributions, which justifies their similar performances and bounds. In
conclusion, for the KdV equation, the positive effect of simplicity of target
functions in every subdomain brought by the domain decomposition is similar to
the negative overfitting effect caused by less available data in each
subdomain, which leads to similar performance of XPINN than PINN overall.

###  5.2 Heat Equation

####  5.2.1 Setup

In this subsection, we consider the heat equation, which is a second order
linear PDE. The one-dimensional heat equation under consideration is
$u_{t}=u_{xx},x\in[-1,1],t\in[0,1]$, where its boundary conditions on
$(x,t)\in\\{-1,1\\}\times[0,1]$, and its initial conditions on
$(x,t)\in[0,1]\times\\{0\\}$, are given by the ground truth solution

$u(x,t)=e^{-\pi^{2}t}\cos(\pi x)+0.6e^{-4\pi^{2}t}\cos(2\pi
x)+0.3e^{4t-4}\cosh(2x)+0.1e^{t-1}\sinh x.$ (29)

Figure 4: Data visualization for the heat experiment. Left: exact solution of
the heat equation. Right: Training points.

Figure 5: Error visualization for the heat experiment. Table 2: Computational results for the heat equation. | Train Loss | Relative $L_{2}$ error | Complexity | Bound
---|---|---|---|---
PINN | 8.589e-5$\pm$2.218e-5 | 1.778e-3$\pm$2.195e-4 | 100.00% | 100.00%
XPINN-T | 2.585e-4$\pm$1.726e-4 | 4.490e-3$\pm$1.517e-3 | 156.24% | 243.22%
XPINN-B | 75.75%

The solution is visualized in Figure 4. The training dataset for PINN contains
2000 residual points and 200 boundary points, whereas the testing dataset for
PINN contains 160801 points within the domain. The backbone model for PINN is
a 9-layer neural network with 20 hidden units activated by tanh. LBFGS with
1e-1 learning rate is used. No regularization is used.

For domain decomposition of XPINN, from Figure 4 left, we observe that the
solution of heat equation is complex near $t=0$ and $t=1$, due to the two
nearby heat sources. The two heat sources are also dissimilar: at $t=0$ the
source is generated by trigonometric functions, while at $t=1$ it is generated
by hyperbolic functions. To design a good XPINN, we should partition the
complexities of the two heterogeneous heat sources into different subdomains.
Thus, we partition the whole domain into a bottom domain $t\leq 0.5$
containing trigonometric heat source, corresponding to XPINN-B, and a top one
$t>0.5$ containing the hyperbolic source, corresponding to XPINN-T.

For fair comparison, we keep the same training procedure, i.e., training
epochs, learning rate, model structure, and weight decay, etc. We train each
model for 20000 epochs, which can train orders of magnitude to $1e-4$ to
$1e-5$, and the results reported in the table are those at the 20000-th
epochs. For both models, we use unity weight for the residual loss and the
residual continuity loss, and use 20 weight for the boundary loss and the
boundary interface loss, which has been adopted in the original code of XPINN
[12]. For reproducibility, we run each model for 5 times using fixed random
seed 0, 1, 2, 3, 4.

####  5.2.2 Results

Table 2 shows the experimental results for the heat equation. In the table,
Train Loss denotes the training boundary plus residual losses, Relative L2
denote the relative L2 error of the model, Complexity denotes the products of
norms of the weight matrices to quantify the network complexity, and Bound
denotes the theoretical bound for the (relative) L2 error.

PINN (1.778e-3) outperforms XPINN (4.490e-3) in relative L2 error, and its
bound (100.00%) is also smaller than that of XPINN (243.22%), i.e., our bound
is consistent with numerical results. For complexity, the order is XPINN-T
(156.24%) $>$ PINN (100%) $>$ XPINN-B (75.75%). XPINN-T is the most complex
because the top heat source is more complicated than the bottom one, and
XPINN-T only has half of the training data. In this case, the negative effect
due to overfitting caused by less data available is more obvious than the
positive influence of less complex target function, so PINN performs better
than XPINN. In Figure 5, we visualize the errors of PINN and XPINN. The errors
of both models are distributed mainly in the bottom part of the domain, where
XPINN-B performs much worse than PINN in the corners of the bottom sub-domain,
which may be due to limit data in this sub-domain, finally deteriorating the
performance in this part.

###  5.3 Advection Equation

####  5.3.1 Setup

Figure 6: Data visualization for the advection experiment. Left: exact
solution of the advection equation. Right: Training points.

Figure 7: Error visualization for the advection equation.

In this subsection, we consider the advection equation to show the difference
between XPINN and PINN, which is given by
$u_{t}+0.5u_{x}=0,x\in[-1,1],t\in[0,1]$, with the initial condition
$u(x,0)=1_{-0.2\leq x\leq 0.2}$. The solution is presented in Figure 6 left.
The training dataset for PINN contains 2000 residual points and 200 boundary
points. The backbone model for PINN is a 6-layer neural network with 20 hidden
units activated by tanh. Adam [31] optimizer with 1e-3 learning rate is used
for optimization. No regularization is used.

For XPINN, from Figure 6 left, we observe that the solution of the advection
equation can be divided into the following two parts,
$\left\\{-0.2<x-0.5t\right\\}\cap\left\\{x-0.5t<0.2\right\\}$ where $u=1$ and
$\left\\{x-0.5t\geq 0.2\right\\}$ and $\left\\{x-0.5t\leq-0.2\right\\}$ where
$u=0$. Hence, in XPINN-LMR (left, middle, and right), we partition the domain
into these three continuous parts mentioned above. XPINN-LMR seems to be a
good partition since in each subdomain the target function is extremely simple
constant function, while the whole function is complex and discontinuous. The
visualization of domain decomposition is provided in Figure 6 right.

For fair comparison, we keep the same training procedure. We train each model
for 5000 epochs, and the results reported in the table are those at the
5000-th epochs. For both models, we use unity weight for the residual loss and
zero weight for the residual continuity loss, and also use unity weight for
the boundary loss and the boundary interface loss. For reproducibility, we run
each model for 5 times using fixed random seed 0, 1, 2, 3, 4.

####  5.3.2 Results

Table 3: Computational results for advection equation. Method | Train Loss | Relative $L_{2}$ error | Complexity | Bound
---|---|---|---|---
PINN | 1.387e-5$\pm$1.298e-5 | 2.052e-1$\pm$1.001e-1 | 100% | 100%
XPINN-L | 4.239e-3$\pm$2.385e-5 | 1.617e-1$\pm$3.582e-2 | 40.53% | 66.59%
XPINN-M | 53.16%
XPINN-R | 79.95%

Table 3 presents all computational results for the advection equation. XPINN
(2.052e-1) performs better than PINN (1.617e-1) in relative L2 error, which is
consistent with their theoretical bounds, i.e., the bound of XPINN (66.59%) is
also smaller than that of PINN (100%). The reason is revealed by the norms.
All sub-nets in XPINN (40.53%, 53.16%,79.95%) are less complicated than PINN
(100%), because in each subdomain of XPINN, the target function is constant,
whose positive influence is more obvious than the overfitting due to less
data. In sum, in the experiment of advection equation, since the target
functions in every sub-domains are extremely simple constant functions, the
positive influence of less complex target functions is much stronger than the
negative effect of overfitting due to less available data in each sub-domains,
thus, XPINN performs better than PINN. Lastly, in Figure 7 we visualize the
errors of PINN and XPINN. The error mainly concentrates in the areas near the
discontinuity part, where XPINN has smaller error than PINN.

###  5.4 Poisson Equation

Figure 8: Data visualization for the Poisson experiment. Left: exact solution
of the Poisson equation. Right: Training points.

####  5.4.1 Setup

In this subsection, we consider a Poisson equation with residual
discontinuity, which is also a second order linear PDE. The equation under
consideration is $u_{xx}+u_{yy}=f$, where $(x,y)\in[0,1]\times[0,1]$, and $f$
is given by $f(x,y)=1$ for $(x,y)\in[0.25,0.75]\times[0.25,0.75]$, and
$f(x,y)=0$ for the rest of the domain. The boundary condition is zero.

The solution is visualized in Figure 8. The training dataset for PINN contains
400 residual points and 80 boundary points, whereas the testing dataset for
PINN contains 1,002,001 points within the domain. The backbone model for PINN
is a 9-layer neural network with 20 hidden units activated by tanh. LBFGS with
1e-1 learning rate is used. No regularization is used.

The weighting strategies for PINN and XPINNs are summarized in Table 4. For
PINN, we use unity weight for residual and 20 weight for boundary. For XPINN1,
i.e., the simplest XPINN model, we use 20 weight for boundary, boundary
interface, and residual interface, while using unity weight for residual. Due
to the existence of residual discontinuity, we need to give stronger
constraint to maintain the residual continuity. So, we use 20 weight for
residual interface rather than unity weight used in the official XPINN code
[12].

In addition, to test the effective of the newly proposed regularization in
[17], we design the second XPINN2, which builds upon XPINN1, i.e., it adopts
the same weights as XPINN1 and use 30 weight for the additional regularization
on the first order derivatives near the interface.

For the last XPINN3, we shall see in the results that XPINN2 does not perform
well on the boundary. So, we try to remedy XPINN2 by XPINN3 through increasing
the weight for the boundary loss, i.e., from 20 to 80 weight on the boundary.

For domain decomposition of XPINN, sub-domain 1 contains the area
$(x,y)\in[0.25,0.75]\times[0.25,0.75]$, while sub-domain 2 contains the rest
of the domain. Specifically, we partition the domain according to the
discontinuity, to force each sub-net in the XPINN to focus on one continuous
part, rather than fitting the entire function containing residual
discontinuity by one network. The sub-net for the sub-domain 1 is called
XPINN-M since it is in the middle of the domain, while the sub-net for sub-
domain 2 is called XPINN-A since it is around the entire domain.

For fair comparison, we keep the same training procedure, i.e., training
epochs, learning rate, model structure, and weight decay, etc. We train each
model for 20000 epochs, and the results reported in the table are those at the
20000-th epochs. For reproducibility, we run each model for 5 times using
fixed random seed 0, 1, 2, 3, 4.

Table 4: Weighting strategies for the Poisson equation. | Residual | Interface R | Additional I | Boundary | Interface B
---|---|---|---|---|---
PINN | 1 | NA | NA | 20 | NA
XPINN1 | 1 | 20 | 0 | 20 | 20
XPINN2 | 1 | 20 | 30 | 20 | 20
XPINN3 | 1 | 20 | 30 | 80 | 20
Table 5: Computational results for the Poisson equation. | Train Loss | Relative $L_{2}$ error | Complexity | Bound
---|---|---|---|---
PINN | 2.688e-4$\pm$3.411e-4 | 5.553e-2$\pm$2.936e-2 | 100.00% | 100.00%
XPINN1-A | 1.181e-2$\pm$4.319e-3 | 4.022e-1$\pm$1.648e-1 | 142.71% | 122.56%
XPINN1-M | 297.91%
XPINN2-A | 1.016e-2$\pm$3.713e-3 | 1.387e-1$\pm$7.030e-3 | 183.44% | 108.57%
XPINN2-M | 292.93%
XPINN3-A | 1.621e-2$\pm$5.222e-3 | 1.108e-1$\pm$1.561e-2 | 195.57% | 106.28%
XPINN3-M | 300.47%

Figure 9: Error visualization in the Poisson equation.

####  5.4.2 Results

Table 5 shows the experimental results for the Poisson equation. PINN
(5.553e-2) outperforms XPINN3 (1.108e-1), and XPINN3 outperforms XPINN2
(1.387e-1), and finally XPINN2 outperforms XPINN1 (4.022e-1), and their bounds
also point to the same result, i.e., PINN (100%) $<$ XPINN3 (106.28%) $<$
XPINN2 (108.57%) $<$ XPINN1 (122.56%). All sub-nets in the XPINNs (142.71%,
297.91%, 183.44%, 292.93%, 195.57%, 300.47%) are more complicated then the
PINN (100%), which is the reason accounting for the failure of XPINNs. In this
case, the lack of available data in each sub-domain of XPINNs impacts the
generalization negatively, whose effect is much more obvious than the
reduction of target function complexity in each sub-domain.

To understand the failure of XPINNs, we visualize the error plots in Figure 9.
For XPINN1, due to the lack of the additional residual interface
regularization on the interface [17], the error is extremely large near the
interface, which causes the largest error of XPINN1. For XPINN2, although the
additional residual regularization does mitigate the error near the interface,
the boundary error becomes much larger, which is because the introduction of
an additional regularization decreases the importance of boundary loss during
the optimization. To solve the problem, in XPINN3 we add more weight to the
boundary loss. While the error on the boundary decreases, the error near the
interface increases, which is due to the tradeoff between different loss
components, namely the tradeoff between the boundary loss and the interface
loss. To conclude, XPINNs perform worse than PINN since they perform bad
either on the boundary or on the interface, which justifies the necessity of
XPINNs with adaptive weight adjustment.

###  5.5 Compressible Euler Equations

####  5.5.1 Setup

Next, we consider the nonlinear inviscid compressible Euler equations, which
govern the physics of high-speed compressible fluid flows. The inviscid
compressible Euler equations admit discontinuous solutions called shock or
contact waves, which are difficult to capture with good accuracy. The two-
dimensional steady-state Euler equations are given as
$F_{x}(U)+F_{y}(U)=0,~{}~{}(x,y)\in[0,1]^{2}$, where fluxes in $x$ and $y$
directions are defined as $F_{x}(U)=(\rho u,p+\rho u^{2},\rho uv,pu+\rho uE)$
and $F_{y}(U)=(\rho v,\rho uv,p+\rho v^{2},pv+\rho vE)$,

Figure 10: Two-dimensional compressible Euler equations: Exact solution (left)
and domain decomposition for XPINN-AM (right).

where $\rho,u,v$ and $p$ are density, velocity components in $x$ and $y$
directions, and pressure, respectively. The total energy $E$ is defined as
$E=\frac{p}{\rho(\gamma-1)}+\frac{1}{2}||\mathbf{u}||_{2}^{2}$, where
$\mathbf{u}=(u,v)$. In this case, we are solving the oblique shock wave
problem on a square domain $[0,1]^{2}$. The bottom boundary is the wall where
slip boundary conditions are applied, whereas left and top boundary are the
inflow boundary where Dirichlet boundary conditions are applied. The right
boundary has extrapolation boundary conditions. A Mach 2 flow is at an angle
of -10 degrees with respect to the bottom wall, which generates an oblique
shock at an angle of 29.3 degrees with the bottom horizontal wall. The exact
solution is given as

$(\rho,u,v,p)=\begin{cases}(1.0,~{}\cos 10^{o},-\sin
10^{o},~{}0.17857)&\text{before shock},\\\
(1.4584,~{}0.8873,~{}0.0,~{}0.3047)&\text{after shock}.\end{cases}$

Among all the primitive variables, we have plotted the fluid density which
accurately shows the position of an oblique shock wave. Figure 10 shows the
exact value of density (left) and domain decomposition (right) for XPINN-AM.
From the domain decomposition figure we observe that the solution of the Euler
equations is divided into the following two parts, $\left\\{y\geq
0.57x+0.1\right\\}\cup\left\\{y\leq 0.5222x-0.0522\right\\}$, where solution
is constant (shown by blue points) and the remaining strip (shown by red
points), where oblique shock wave is present.

####  5.5.2 Results

Figure 11: Two-dimensional compressible Euler equations: Point-wise errors for PINN (left) and XPINN (right). Table 6: Computational results for compressible Euler equations. Method | Train Loss | Relative $L_{2}$ error in $\rho$ | Norms | Bound
---|---|---|---|---
PINN | 1.819e-3$\pm$6.043e-4 | 3.4604e-2$\pm$7.385e-3 | 100.00% | 100.00%
XPINN-A | 9.210e-4$\pm$1.882e-4 | 1.048e-2$\pm$5.3793e-3 | 37.28% | 81.09%
XPINN-M | 64.37 %
XPINN-T | 1.067e-3$\pm$4.829e-4 | 3.5722e-2$\pm$4.290e-3 | 42.37% | 137.63%
XPINN-B | 131.26%

We used a deep net with 10000 residual points, 5 hidden-layers with 20 neurons
in each layer, and 8e-4 learning rate. The activation function is hyperbolic
tangent. Table 6 gives the computational results, and Figure 11 gives the
point-wise error for the density of the fluid. In this case, the XPINN-AM
generalizes better than PINN. Furthermore, the norms of XPINN-AM are much
smaller than that of PINN (100%). We further divide the domain into top
($y\geq 0.5$) and bottom subdomains ($y<0.5$) for XPINN-TB. The XPINN-TB does
not generalize well compared to PINNs, and the complexities of the two sub-
nets are 42.37% and 131.26% for top and bottom subdomains, respectively. These
results prove that the norms and generalization bounds are good indicators for
XPINN based domain decomposition, and can be efficiently used to further
decompose the subdomains.

##  6 Conclusion

In this study, we have investigated the generalization abilities of PINNs and
XPINNs, as well as when and how XPINNs improve generalization. For this
purpose, we have provided both prior and posterior generalization bounds to
explain this from different viewpoints, where for the former we have developed
the Barron space for multi-layer networks, while for the latter we have
derived the complexity of norm-based Rademacher for PINNs.

Through our discussion on theoretical results, analytical examples, and
extensive experiments, we conclude that the domain decomposition method in
XPINNs introduces a tradeoff on generalization. On the one hand, its advantage
is that it decomposes the complex target function into several simple parts,
which lead to the phenomenon that the sum of all parts is smaller than the
whole. However, on the other hand, domain decomposition causes less available
training data in each subdomain, leading to higher empirical Rademacher
complexity and makes models prone to overfitting. When the complexity
reduction brought by XPINN exceeds the increased complexity caused by less
training data, XPINN outperforms PINN, as in our experiment on KdV equation,
heat equation, advection equation, compressible Euler equation and the
analytical example in Section 4.1. When the overfitting caused by insufficient
data is more dominant than the simplicity due to domain decomposition, PINN
outperforms XPINN, as in our experiment on heat equation, wave equation and
our analytical example in Section 4.2. When the two factors reach a balance,
XPINN and PINN perform similarly, as shown in our experiment on advection
equation and our analytical example in Section 4.3.

Our results can also provide a partial explanation for the following
observation. For long-time integration of several PDEs, it has been
empirically observed that only XPINNs are applicable as PINNs tend to be
inaccurate. According to our theory, this is expected as we tend to have a
very high complexity measured by the norms for a whole solution of a long-time
integration, which can be decomposed into less complex sub-solutions in
XPINNs. Our proposed theory can also be useful for adaptive domain
decomposition. Specifically, after initialization of a decomposition, we can
compute the bound during optimization to know how the corresponding XPINN
generalizes. If it does not generalize well, we can diagnose the reason, e.g.,
the negative overfitting effect caused by less training points is more obvious
than the positive effect of less complex target function in each sub-domain,
then we can restructure the decomposition by letting the sub-domain(s) contain
more training points. Overall, the present work provides the first theoretical
understanding on when and how to employ XPINN for better generalization
performances over the vanilla PINNs.

## Acknowledgment

A. D. Jagtap and G. E. Karniadakis would like to acknowledge the funding by
OSD/AFOSR MURI Grant FA9550-20-1-0358, and the US Department of Energy (DOE)
PhILMs project (DE-SC0019453). The authors would like to thank Michael
Penwarden for his insightful comments and reviews on the paper and for his
generous help with the code for experiments.

##  Appendix A Preliminary: Functional Analysis

######  Definition A.1.

(Multi-Index) If $k$ is a positive integer, and $u\in C^{k}(\Omega)$, then we
define the multi-index $\alpha=(\alpha_{1},\cdots,\alpha_{d})$ of order
$|\alpha|=\alpha_{1}+\cdots+\alpha_{d}=k$, and the corresponding derivative
is:

$D^{\alpha}u=\frac{\partial^{\alpha_{1}}}{\partial
x_{1}^{\alpha_{1}}}\cdots\frac{\partial^{\alpha_{n}}}{\partial
x_{d}^{\alpha_{d}}}u.$ (30)

######  Definition A.2.

(Holder space) The Holder space $C^{k,\gamma}(\overline{\Omega})$ consists of
all functions $u\in C^{k}(\overline{\Omega})$ for which the norm

$\|u\|_{C^{k,\gamma}(\overline{\Omega})}:=\sum_{|\alpha|\leq
k}\|D^{\alpha}u\|_{C(\overline{\Omega})}+\sum_{|\alpha|=k}[D^{\alpha}u]_{C^{0,\gamma}(\overline{\Omega})}$
(31)

is finite, where the norm and the semi-norm are defined as

$\displaystyle\|u\|_{C(\overline{\Omega})}$
$\displaystyle:=\sup_{x\in\Omega}|u(x)|.$ (32)
$\displaystyle[u]_{C^{0,\gamma}(\overline{\Omega})}$
$\displaystyle:=\sup_{x,y\in\Omega}\left\\{\frac{|u(x)-u(y)|}{|x-y|^{\gamma}}\right\\}.$

######  Theorem A.1.

The space of functions $C^{k,\gamma}(\overline{\Omega})$ is a Banach space.

######  Definition A.3.

(Completeness) A metric space $X$ is called complete if every Cauchy sequence
in $X$ has a limit that is also in $X$.

######  Definition A.4.

(Relative compactness) A relatively compact subspace $Y$ of a topological
space $X$ is a subset whose closure is compact.

######  Lemma A.1.

Suppose $(f_{n})_{n\in\mathbb{N}}$ is an equicontinuous sequence in
$C(\overline{\Omega})$, then if $(f_{n})_{n\in\mathbb{N}}$ converges to $f$
pointwise, then $f\in C(\overline{\Omega})$ and the convergence is uniform.

###### Proof.

Suppose $x\in\overline{\Omega}$, since $(f_{n})_{n\in\mathbb{N}}$ is an
equicontinuous sequence in $C(\overline{\Omega})$, there exists
$V_{x}\in\mathcal{N}(x)$, where $\mathcal{N}(x)$ denotes the neighborhood of
$x$, such that when $y\in V_{x}$, for all $n\in\mathbb{N}$ we have

$|f_{n}(y)-f_{n}(x)|\leq\epsilon.$ (33)

Therefore,

$\displaystyle|f(y)-f(x)|$
$\displaystyle\leq|f(y)-f_{n}(y)|+|f_{n}(x)-f(x)|+|f_{n}(x)-f_{n}(y)|$ (34)
$\displaystyle\leq|f(y)-f_{n}(y)|+|f_{n}(x)-f(x)|+\epsilon.$

Let $n\rightarrow\infty$, $|f(y)-f(x)|\leq\epsilon$. Thus, $f$ is continuous.

Since the set $\overline{\Omega}$ is compact, there exists
$x_{1},\cdots,x_{m}$ and their neighborhoods $V_{x_{1}},\cdots,V_{x_{m}}$,
such that

$\overline{\Omega}=\cup_{k=1}^{m}V_{x_{k}}.\qquad\forall
x\in\overline{\Omega},\exists x_{k},x\in
V_{x_{k}},|f_{n}(x)-f_{n}(x_{k})|\leq\epsilon,|f(x)-f(x_{k})|\leq\epsilon,$
(35)

where $m$ is large enough for the given $\epsilon$, and we consider the
equicontinuity of $f_{n}$ and the continuity of $f$. For all
$x\in\overline{\Omega}$, there exists $x_{k}$ such that $x\in V_{x_{k}}$.
Given large enough $n$, due to the continuity of $f_{n}$ and $f$, we have

$|f_{n}(x)-f(x)|\leq|f_{n}(x)-f_{n}(x_{k})|+|f_{n}(x_{k})-f(x_{k})|+|f(x_{k})-f(x)|\leq
3\epsilon.$ (36)

Thus the convergence is uniform. ∎

######  Theorem A.2.

(Arzela-Ascoli theorem) Suppose the funciton class $\mathcal{H}\subset
C(\overline{\Omega})$, then $\mathcal{H}$ is relatively compact in
$C(\overline{\Omega})$ if and only if $\mathcal{H}$ is equicontinuous.

###### Proof.

We shall only use the sufficiency part of this theorem. So, only that part
will be proved.

We first prove that $\overline{\mathcal{H}}$ is complete in $C(\Omega)$. We
select arbitrary Cauchy sequence $(f_{n})_{n\geq 1}$ in
$\overline{\mathcal{H}}$. Then for all $\boldsymbol{x}\in\overline{\Omega}$,
$(f_{n}(\boldsymbol{x}))_{n\geq 1}$ is a Cauchy sequence in the Euclidean
space $\mathbb{R}$, which is convergent. By Lemma A.1, $(f_{n})_{n\geq 1}$
uniformly converges to a continuous $f$. Thus, $\overline{\mathcal{H}}$ is
complete.

We only need to show that, for all $\epsilon>0$, $\overline{\mathcal{H}}$ can
be covered by finite balls of radius $\epsilon$. Since $\mathcal{H}$ is
equicontinuous, then for all $\boldsymbol{x}\in\overline{\Omega}$ and
$\epsilon>0$, there exists open set $O_{x}\in\mathcal{N}(x)$, where
$\mathcal{N}(x)$ denotes the neighbourhood of $x$, such that for all $y\in
O_{x}$, we have $|f(y)-f(x)|\leq\epsilon$ for all $f\in\mathcal{H}$. Since
$\overline{\Omega}$ is a compact set, there exists finite
$(O_{x_{i}})_{i=1}^{n}$, such that
$\overline{\Omega}=\cup_{i=1}^{n}O_{x_{i}}$.

Since all $f\in\mathcal{H}$ are continuous, the set
$\cup_{i=1}^{n}\mathcal{H}(x_{i})$ is compact, so there exists finite set
$Z\subset\overline{\Omega}$, such that

$\cup_{i=1}^{n}\mathcal{H}(x_{i})\subset\cup_{z\in Z}B(z,\epsilon),$ (37)

where $\mathcal{H}(x):=\\{f(x):f\in\mathcal{H}\\}$. Now we consider the finite
set $Z^{n}$, for all $\overline{z}=(z_{1},\cdots,z_{n})\in Z^{n}$. Let

$B_{\overline{z}}=\left\\{f\in C(\overline{\Omega}):\sup_{1\leq i\leq
n}\sup_{x\in O_{x_{i}}}|f(x)-z_{i}|<2\epsilon\right\\}.$ (38)

Then we can a finite number of open subsets
$(B_{\overline{z}})_{\overline{z}\in Z^{n}}$.

For all $f,g\in B_{\overline{z}}$, we have

$\|f-g\|_{C(\overline{\Omega})}\leq\sup_{1\leq i\leq n}\sup_{x\in
O_{x_{i}}}\left(|f(x)-z_{i}|+|g(x)-z_{i}|\right)\leq 4\epsilon.$ (39)

Thus $\text{diam}B_{\overline{z}}\leq 4\epsilon$. We want to show
$\mathcal{H}\subset\cup_{\overline{z}\in Z^{n}}B_{\overline{z}}$.

For all $f\in\mathcal{H}$, we know that for all $1\leq i\leq n$, $f(x_{i})$
belongs to one of $B(z_{i},\epsilon)$ where $z_{i}\in Z$, due to the fact that
$\cup_{i=1}^{n}\mathcal{H}(x_{i})\subset\cup_{z\in Z}B(z,\epsilon)$. When
$x\in O_{x_{i}}$, then

$|f(x)-z_{i}|\leq|f(x)-f(x_{i})|+|f(x_{i})-z_{i}|<\epsilon+\epsilon=2\epsilon.$
(40)

Therefore,

$\sup_{1\leq i\leq n}\sup_{x\in O_{x_{i}}}|f(x)-z_{i}|<2\epsilon.$ (41)

∎

######  Lemma A.2.

(Lebesgue’s dominated convergence theorem) Let $(f_{n})_{n\in\mathbb{N}}$ be a
sequence of measurable functions in the space $\overline{\Omega}$ with measure
$\mu$. Suppose that the sequence converges point-wise to a function $f$ and is
dominated by some integrable function $g$ in the sense that $|f_{n}(x)|\leq
g(x)$, for all numbers $n$ in the index set of the sequence and all points
$x\in\overline{\Omega}$. Then $f$ is integrable (in the Lebesgue sense) and

$\lim_{n\rightarrow\infty}\int_{\overline{\Omega}}|f_{n}(x)-f(x)|d\mu(x)=0.$
(42)

######  Definition A.5.

(Sobolev space) The Sobolev space $W^{k,p}(\Omega)$ contains all locally
summable functions $u:\Omega\rightarrow\mathbb{R}$ such that for each multi-
index $\alpha$ with $|\alpha|\leq k$, $D^{\alpha}u$ exists in the weak sense
and belongs to $L^{p}(\Omega)$. Furthermore, if $u\in W^{k,p}(\Omega)$, we
define its norm to be

$\displaystyle\|u\|_{W^{k,p}(\Omega)}$
$\displaystyle:=\left(\sum_{|\alpha|\leq
k}\int_{\Omega}|D^{\alpha}u|^{p}\right)^{1/p}.\qquad p\in[1,\infty).$ (43)
$\displaystyle\|u\|_{W^{k,p}(\Omega)}$ $\displaystyle:=\sum_{|\alpha|\leq
k}\text{ess sup}_{\Omega}|D^{\alpha}u|.\qquad p=\infty.$

######  Theorem A.3.

When $p=2$, the Sobolev space $H^{p}(\Omega)=W^{k,2}(\Omega)$ is a Hilbert
space.

######  Theorem A.4.

(Trace Theorem) Assume $\Omega$ is bounded and $\Omega$ is $C^{1}$. Then there
exists a bounded linear operator:

$T:W^{1,p}(\Omega)\rightarrow L^{p}(\partial\Omega)$ (44)

such that

  * •

$Tu=u|_{\partial\Omega}$ if $u\in W^{1,p}(\Omega)\cap C(\overline{\Omega})$.

  * •

For each $u\in W^{1,p}(\Omega)$ and a constant $C$ depending only on $p$ and
$\Omega$,

$\|Tu\|_{L^{p}(\partial\Omega)}\leq C\|u\|_{W^{1,p}(\Omega)}.$ (45)

##  Appendix B Proofs of the Barron Space

###  B.1 Proof of Theorem 2.2

###### Proof.

(Proof of Theorem 2.2) Because $X$ embeds continuously into $C^{2,1}(\Omega)$,
there exists constants $C_{1},C_{2}>0$ such that

$\|D^{\alpha}g\|_{C^{0}(\Omega)}\leq C_{1}\|g\|_{X},\
[D^{\alpha}g]_{C^{0,1}(\Omega)}\leq C_{2}\|g\|_{X},\ \forall g\in X.$ (46)

Banach Space. By construction, $\mathcal{B}_{X,\Omega}$ is isometric to the
quotient space $\mathcal{M}(B^{X})/N_{K}$ where

$N_{K}=\left\\{\mu\in\mathcal{M}(B^{X})|\int_{B^{X}}\rho(g)d\mu(g))=0,\forall\boldsymbol{x}\in\Omega\right\\}.$
(47)

In particular, $\mathcal{B}_{X,\Omega}$ is a normed vector space with the norm
$\|\cdot\|_{X,\Omega}$. Consider the mapping

$\mathcal{M}(B^{X})\rightarrow C^{2}(\overline{\Omega}),\ \mu\rightarrow
f_{\mu}=\int_{B^{X}}\rho(g)d(\mu(g)).$ (48)

We prove that $f_{\mu}\in C^{2}(\overline{\Omega})$. For the continuity, we
have

$\displaystyle|f_{\mu}(x)-f_{\mu}(y)|$
$\displaystyle=\int_{B^{X}}|\rho(g(x))-\rho(g(y))|d(\mu(g))$ (49)
$\displaystyle\leq\int_{B^{X}}|g(x)-g(y)|d(\mu(g)).$

For the continuity at point $x$, consider arbitrary sequence
$(x_{n})_{n\in\mathbb{N}}\rightarrow x$. Consider the sequence of function
$h_{n}:B^{X}\rightarrow\mathbb{R}$, $f\mapsto f(x_{n})$ and
$h:B^{X}\rightarrow\mathbb{R}$, $f\mapsto f(x)$. We have that
$|h_{n}(f)|=|f(x_{n})|\leq\|f\|_{C^{0}(\overline{\Omega})}$, which is
integrable. By Lemma A.2,

$\displaystyle\lim_{n\rightarrow\infty}|f_{\mu}(x)-f_{\mu}(x_{n})|$
$\displaystyle\leq\lim_{n\rightarrow\infty}\int_{B^{X}}|g(x)-g(x_{n})|d(\mu(g))$
(50)
$\displaystyle=\lim_{n\rightarrow\infty}\int_{B^{X}}|h(g)-h(g_{n})|d(\mu(g))$
$\displaystyle=\int_{B^{X}}\lim_{n\rightarrow\infty}|h(g)-h(g_{n})|d(\mu(g))$
$\displaystyle=\int_{B^{X}}\lim_{n\rightarrow\infty}|g(x)-g(x_{n})|d(\mu(g))$
$\displaystyle=0,$

due to the continuity of $g$. Since the sequence $(x_{n})_{n\in\mathbb{N}}$ is
arbitrary, we know that $f_{\mu}$ is continuous at all points
$x\in\overline{\Omega}$.

$f_{\mu}$ is first order differentiable due to the following:

$\displaystyle\frac{f_{\mu}(x)-f_{\mu}(x_{n})}{x-x_{n}}-\int_{B^{X}}\frac{\partial\rho(g(x))}{\partial
x_{i}}d(\mu(g))$
$\displaystyle=\int_{B^{X}}\left[\frac{\rho(g(x))-\rho(g(x_{n}))}{x-x_{n}}-\frac{\partial\rho(g(x))}{\partial
x_{i}}\right]d(\mu(g)),$ (51)

where $x_{i}$ is the $i$-th coordinate of $x$, $(x_{n})_{i}\rightarrow x_{i}$
as $n\rightarrow\infty$, and $(x_{n})_{-i}=x_{-i}$, i.e., other coordinates
are the same. Because $X$ embeds continuously into $C^{2,1}(\Omega)$, we have

$\displaystyle\left|\frac{\rho(g(x))-\rho(g(x_{n}))}{x-x_{n}}\right|$ (52)
$\displaystyle\leq\left|\frac{g(x)-g(x_{n})}{x-x_{n}}\right|$
$\displaystyle\leq C_{2}.$

Since the activation function and its derivatives are bounded, we also have

$\displaystyle\left|\frac{\partial\rho(g(x))}{\partial x_{i}}\right|$
$\displaystyle=\left|\rho^{\prime}(g(x))\frac{\partial\rho(g(x))}{\partial
x_{i}}\right|$ (53) $\displaystyle\leq\left|\frac{\partial\rho(g(x))}{\partial
x_{i}}\right|$ $\displaystyle\leq C_{1}.$

By Lemma A.2,

$\displaystyle\lim_{n\rightarrow\infty}\left|\frac{f_{\mu}(x)-f_{\mu}(x_{n})}{x-x_{n}}-\int_{B^{X}}\frac{\partial\rho(g(x))}{\partial
x_{i}}d(\mu(g))\right|$ (54)
$\displaystyle\leq\int_{B^{X}}\lim_{n\rightarrow\infty}\left|\frac{\rho(g(x))-\rho(g(x_{n}))}{x-x_{n}}-\frac{\partial\rho(g(x))}{\partial
x_{i}}\right|d(\mu(g))$ $\displaystyle=0.$

Thus

$\frac{\partial f_{\mu}(x)}{\partial
x_{i}}=\int_{B^{X}}\frac{\partial\rho(g(x))}{\partial x_{i}}d(\mu(g)).$ (55)

Then, by induction we can know that $f_{\mu}$ is $C^{k}$ continuous and

$D^{\alpha}f_{\mu}=\int_{B^{X}}D^{\alpha}gd(\mu(g)).$ (56)

Back to the mapping

$\mathcal{M}(B^{X})\rightarrow C^{2}(\overline{\Omega}),\ \mu\rightarrow
f_{\mu}=\int_{B^{X}}\rho(g)d(\mu(g)).$ (57)

It is continuous because

$\displaystyle\qquad\left\|\int_{B^{X}}\rho(g)d(\mu(g))\right\|_{C^{2}(\Omega)}$
(58) $\displaystyle\leq\int_{B^{X}}\|\rho(g)\|_{C^{2}(\Omega)}d(|\mu|(g))$
$\displaystyle=\int_{B^{X}}\left\\{\|\rho(g)\|_{C^{0}(\Omega)}+\sum_{i=1}^{d}\|\rho^{\prime}(g)\partial_{i}g\|_{C^{0}(\Omega)}\right\\}d(|\mu|(g))$
$\displaystyle\quad+\int_{B^{X}}\left\\{\sum_{i,j=1}^{d}\|\rho^{\prime\prime}(g)\partial_{i}g\partial_{j}g\|_{C^{0}(\Omega)}+\sum_{i,j=1}^{d}\|\rho^{\prime}(g)\partial_{ij}g\|_{C^{0}(\Omega)}\right\\}d(|\mu|(g))$
$\displaystyle\leq\int_{B^{X}}\left\\{1+\sum_{i=1}^{d}\|\partial_{i}g\|_{C^{0}(\Omega)}+\sum_{i,j=1}^{d}\|\partial_{i}g\partial_{j}g\|_{C^{0}(\Omega)}+\sum_{i,j=1}^{d}\|\partial_{ij}g\|_{C^{0}(\Omega)}\right\\}d(|\mu|(g))$
$\displaystyle\leq\left[1+(d+d^{2})C_{1}+d^{2}C_{1}^{2}\right]\|\mu\|_{\mathcal{M}(B^{X})},$

by the definition of Brochner integrals and note that $x\in\mathbb{R}^{d}$,
i.e., $d$ is the input dimension, and that $\partial_{i}$ denotes the first
order derivative with respect to the $i$-th coordinate, and
$\partial^{2}_{ij}$ denotes the second order derivative with respect to the
$i$-th and the $j$-th coordinates. Due to the continuity, $N_{K}$ is the
kernel of a continuous linear map. Therefore, $N_{K}$ is a closed subspace of
$\mathcal{M}_{B^{X}}$. By the theorem in functional analysis, we conclude that
$\mathcal{B}_{X,\Omega}$ is a Banach space.

$\mathcal{B}_{X,\Omega}$ embeds continuously into $C^{2,1}(\Omega)$. In the
proof of statement (1), we already have $\|f_{\mu}\|_{C^{2}(\Omega)}\leq
2C_{1}\|\mu\|_{\mathcal{M}(B^{X})}$. By taking infimum over $\mu$, we have
$\|f\|_{C^{2}(\Omega)}\leq 2C_{1}\|f\|_{\mathcal{B}_{X,\Omega}}$. Furthermore,
for any $x\neq y\in\Omega$, we have

$\displaystyle|f_{\mu}(x)-f_{\mu}(y)|$
$\displaystyle\leq\int_{B^{X}}|\rho(g(x))-\rho(g(y))|d|\mu|(g)$ (59)
$\displaystyle\leq\int_{B^{X}}|g(x)-g(y)|d|\mu|(g)$
$\displaystyle\leq\int_{B^{X}}[g]_{C^{0,1}(\overline{\Omega})}|x-y|d|\mu|(g)$
$\displaystyle\leq C_{2}\|\mu\|_{\mathcal{M}(B^{X})}|x-y|.$

For the derivative, we have the similar conclusion:

$\displaystyle|\partial_{i}f_{\mu}(x)-\partial_{i}f_{\mu}(y)|$
$\displaystyle\leq\int_{B^{X}}|\partial_{i}[\rho(g(x))]-\partial_{i}[\rho(g(y))]|d|\mu|(g)$
(60)
$\displaystyle\leq\int_{B^{X}}|\rho^{\prime}(g(x))\partial_{i}g(x)-\rho^{\prime}(g(y))\partial_{i}g(y)|d|\mu|(g)$
$\displaystyle\leq\int_{B^{X}}\left\\{|\rho^{\prime}(g(x))\partial_{i}g(x)-\rho^{\prime}(g(x))\partial_{i}g(y)|+|\rho^{\prime}(g(x))\partial_{i}g(y)-\rho^{\prime}(g(y))\partial_{i}g(y)|\right\\}d|\mu|(g)$
$\displaystyle\leq\int_{B^{X}}\left\\{[g]_{C^{1,1}(\overline{\Omega})}|x-y|+[g]_{C^{0,1}(\overline{\Omega})}\|\partial_{i}g\|_{C^{0}(\overline{\Omega})}|x-y|\right\\}d|\mu|(g)$
$\displaystyle\leq(C_{2}+C_{1}C_{2})\|\mu\|_{\mathcal{M}(B^{X})}|x-y|,$

where $\partial_{i}$ denotes the first order derivative with respect to the
$i$-th coordinate. Also, for the second order derivatives, we have

$\displaystyle|\partial_{ij}f_{\mu}(x)-\partial_{ij}f_{\mu}(y)|$
$\displaystyle\leq\int_{B^{X}}|\partial_{ij}[\rho(g(x))]-\partial_{ij}[\rho(g(y))]|d|\mu|(g)$
(61) $\displaystyle\leq\int_{B^{X}}(I_{1}+I_{2})d|\mu|(g).$

where $\partial^{2}_{ij}$ denotes the second order derivative with respect to
the $i$-th and the $j$-th coordinates, and we have

$\displaystyle I_{1}$
$\displaystyle=|\rho^{\prime\prime}(g(x))\partial_{i}g(x)\partial_{j}g(x)-\rho^{\prime\prime}(g(y))\partial_{i}g(y)\partial_{j}g(y)|$
(62)
$\displaystyle\leq|\rho^{\prime\prime}(g(x))\partial_{i}g(x)\partial_{j}g(x)-\rho^{\prime\prime}(g(x))\partial_{i}g(y)\partial_{j}g(y)|+|\rho^{\prime\prime}(g(x))\partial_{i}g(y)\partial_{j}g(y)-\rho^{\prime\prime}(g(y))\partial_{i}g(y)\partial_{j}g(y)|$
$\displaystyle\leq|\partial_{i}g(x)\partial_{j}g(x)-\partial_{i}g(y)\partial_{j}g(y)|+C_{1}^{2}|\rho^{\prime\prime}(g(x))-\rho^{\prime\prime}(g(y))|$
$\displaystyle\leq|\partial_{i}g(x)\partial_{j}g(x)-\partial_{i}g(x)\partial_{j}g(y)|+|\partial_{i}g(x)\partial_{j}g(y)-\partial_{i}g(y)\partial_{j}g(y)|+C_{1}^{2}|g(x)-g(y)|$
$\displaystyle\leq 2C_{1}C_{2}|x-y|+C_{1}^{2}C_{2}|x-y|$
$\displaystyle\leq(2C_{1}C_{2}+C_{1}^{2}C_{2})|x-y|.$

and

$\displaystyle I_{2}$
$\displaystyle=|\rho^{\prime}(g(x))\partial_{i}\partial_{j}g(x)-\rho^{\prime}(g(y))\partial_{i}\partial_{j}g(y)|$
(63)
$\displaystyle\leq|\rho^{\prime}(g(x))\partial_{i}\partial_{j}g(x)-\rho^{\prime}(g(x))\partial_{i}\partial_{j}g(y)|+|\rho^{\prime}(g(x))\partial_{i}\partial_{j}g(y)-\rho^{\prime}(g(y))\partial_{i}\partial_{j}g(y)|$
$\displaystyle\leq 2C_{2}|x-y|.$

In sum, we have

$\displaystyle|\partial_{ij}f_{\mu}(x)-\partial_{ij}f_{\mu}(y)|$
$\displaystyle\leq(2C_{1}C_{2}+C_{1}^{2}C_{2}+2C_{2})|x-y|.$ (64)

After taking infimum over $\mu$, we come to the conclusion.

The closed unit ball of $\mathcal{B}_{X,\Omega}$ is a closed subset of
$C^{2}(\Omega)$. We assume that $(f_{n})_{n\in\mathbb{N}}$ is a sequence such
that $\|f_{n}\|_{X,K}\leq 1$. Choose a sequence of measures
$(\mu_{n})_{n\in\mathbb{N}}$ such that $f_{n}=f_{\mu_{n}}$ and
$\|\mu_{n}\|\leq 1+\frac{1}{n}$. These meansures exist because $f_{n}$ are
from the unit ball of $\mathcal{B}(X,\Omega)$. By the compactness theorem of
Radon meansures, there exists a subsequence $\mu_{n_{k}}$ and a Radon measure
$\mu$, such that $\mu_{n_{k}}$ weak converge to $\mu$ with $\|\mu\|\leq 1$.
Since all functions $\rho(g)$ are $C^{2}$ continuous and bounded, we have all
$D^{\alpha}f_{n_{k}}$ converge in the product topology (pointwise convergence)
for all $|\alpha|\leq 2$, since $f_{n}\in C^{2}(\overline{\Omega})$. In
particular, if $f_{\mu_{n}}\rightarrow\hat{f}$ uniformly in the norm of
$C^{2}$, then $\hat{f}=f_{\mu}\in B^{\mathcal{B}_{X,\Omega}}$, i.e. the unit
ball of $\mathcal{B}_{X,\Omega}$ is closed in the $C^{2}(\Omega)$ topology. ∎

###  B.2 Proof of Theorem 2.1

###### Proof.

(Proof of Theorem 2.1). This is immediate by the definition of generalized
Barron space and that of neural networks in Definition 2.1. ∎

###  B.3 Proof of Theorem 2.3

To prove the approximation property, we first prove the following useful
lemma.

######  Lemma B.1.

Let $\mathcal{G}$ be a set in a Hilbert space $H$ such that $\|g\|_{H}\leq R$
for all $g\in\mathcal{G}$. If $f$ is in the closed convex hull of $G$, then
for every $m\in\mathbb{N}$ and $\epsilon>0$, there exist $m$ elements
$g_{1},\cdots,g_{m}\in G$ such that

$\|f-\frac{1}{m}\sum_{i=1}^{m}g_{i}\|_{H}\leq\frac{R+\epsilon}{\sqrt{m}}.$
(65)

###### Proof.

This lemma is proved by using the law of large numbers. See [32] for details.
∎

###### Proof.

(Proof of Theorem 2.3) Due to the choice of $\mathcal{W}^{0}$ as linear
functions, the constants of continuous embedding are $C_{1},C_{2}=1$. By the
fact that $C^{2,1}(\overline{\Omega})$ embeds continuously into
$H^{2}(\Omega)$, we have $\|f\|_{H^{2}}\leq 2\|f\|_{\mathcal{W}^{L}}$.

Recall that the unit ball of $\mathcal{W}^{L}$ is the closed convex hull of
the class $\mathcal{H}=\left\\{\sigma(g):\|g\|_{\mathcal{W}^{L}}=1\right\\}$.
Thus by Lemma, there exists $g_{1},\cdots,g_{m}\in\mathcal{W}^{L-1}$ and
$\epsilon_{1},\cdots,\epsilon_{m}\in\left\\{-1,1\right\\}$ such that

$\|f-\frac{1}{m}\sum_{i=1}^{m}\epsilon_{i}\sigma(g_{i}(\boldsymbol{x}))\|_{H^{2}}\leq\frac{3\|f\|_{\mathcal{W}^{L}}}{\sqrt{m}}.$
(66)

If $L=1$, $g_{i}$ are linear functions and
$u_{m}(\boldsymbol{x})=\sum_{i=1}^{m}\frac{\epsilon_{i}}{m}\sigma(g_{i}(\boldsymbol{x}))$
is a two-layer neural netowrk. Thus the case $L=1$ in the theorem is proved.

We prove the remain by induction. Assume that the theorem has been proved for
the case of $L-1$. Then we note that $\|g_{i}\|_{\mathcal{W}^{L}}=1$, so for
$1\leq i\leq m$ we can find a finite $L-1$-layer network $\hat{g}_{i}$ such
that

$\displaystyle\|f-\frac{1}{m}\epsilon_{i}\sigma(\hat{g}_{i}(\boldsymbol{x}))\|_{H^{2}}$
$\displaystyle\leq\|f-\frac{1}{m}\epsilon_{i}\sigma(g_{i}(\boldsymbol{x}))\|_{H^{2}}+\frac{1}{m}\sum_{i=1}^{m}\|g_{i}-\hat{g}_{i}\|_{H^{2}}$
(67) $\displaystyle\leq\frac{3}{\sqrt{m}}+\frac{m}{m}\frac{3(L-1)}{\sqrt{m}}$
$\displaystyle=\frac{3L}{\sqrt{m}}.$

We merge the $m$ trees associated with $\hat{g}_{i}$ into a single tree,
increasing the width of each layer by a factor of $m$, and add an outer layer
of width $m$ with coefficients $W$.

For the remaining part of the theorem, we only need to apply Theorem A.4. ∎

##  Appendix C Proofs of Rademacher Complexity

In this section, we provide proofs for Rademacher complexity. We first provide
full details for the second order derivatives of neural networks.

$\displaystyle\
\frac{\partial^{2}u_{\boldsymbol{\theta}}(\boldsymbol{x})}{\partial\boldsymbol{x}^{2}}$
(68)
$\displaystyle\overset{(1)}{=}\frac{\partial\text{vec}(\boldsymbol{W}^{L}\cdot\boldsymbol{\Phi}^{L-1}\boldsymbol{W}^{L-1}\cdot\dots\cdot\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})}{\partial\boldsymbol{x}}$
$\displaystyle\overset{(2)}{=}\sum_{l=1}^{L-1}\frac{\partial\text{vec}(\boldsymbol{W}^{L}\cdot\boldsymbol{\Phi}^{L-1}\boldsymbol{W}^{L-1}\cdot\dots\cdot\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})}{\partial\text{vec}(\boldsymbol{\Phi}^{l})}\frac{\partial\text{vec}(\boldsymbol{\Phi}^{l})}{\partial\boldsymbol{x}},$
$\displaystyle\overset{(3)}{=}\sum_{l=1}^{L-1}(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})^{\mathrm{T}}\otimes(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})\frac{\partial\text{vec}(\boldsymbol{\Phi}^{l})}{\partial\boldsymbol{x}}$
$\displaystyle\overset{(4)}{=}\left\\{\sum_{l=1}^{L-1}(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})^{\mathrm{T}}\otimes(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})\frac{\partial\text{vec}(\boldsymbol{\Phi}^{l})}{\partial\boldsymbol{x}_{j}}\right\\}_{1\leq
j\leq d}$
$\displaystyle=\\{\sum_{l=1}^{L-1}(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})^{\mathrm{T}}\otimes(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})$
$\displaystyle\qquad\frac{\partial\text{vec}(\text{diag}[\sigma^{\prime}(\boldsymbol{W}^{l}\sigma(\cdots\boldsymbol{W}^{1}\boldsymbol{x}))])}{\partial\boldsymbol{x}_{j}}\\}_{1\leq
j\leq d}$
$\displaystyle=\\{\sum_{l=1}^{L-1}(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})^{\mathrm{T}}\otimes(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})$
$\displaystyle\qquad\text{vec}(\text{diag}[\frac{\partial\sigma^{\prime}(\boldsymbol{W}^{l}\sigma(\cdots\boldsymbol{W}^{1}\boldsymbol{x}))}{\partial\boldsymbol{x}_{j}}])\\}_{1\leq
j\leq d}$
$\displaystyle=\\{\sum_{l=1}^{L-1}(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})^{\mathrm{T}}\otimes(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})$
$\displaystyle\qquad\text{vec}(\text{diag}[\boldsymbol{\Psi}^{l}\boldsymbol{W}^{l}\cdots\boldsymbol{\Psi}^{1}\boldsymbol{W}^{1}_{:,j})])\\}_{1\leq
j\leq d}$
$\displaystyle\overset{(5)}{=}\left\\{\sum_{l=1}^{L-1}(\boldsymbol{W}^{L}\boldsymbol{\Phi}^{L-1}\cdots\boldsymbol{W}^{l+1})\text{diag}(\boldsymbol{\Psi}^{l}\boldsymbol{W}^{l}\cdots\boldsymbol{\Psi}^{1}\boldsymbol{W}^{1}_{:,j})(\boldsymbol{W}^{l}\cdots\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})\right\\}_{1\leq
j\leq d}.$

In (1), note that
$\text{vec}(\boldsymbol{W}^{L}\cdot\boldsymbol{\Phi}^{L-1}\boldsymbol{W}^{L-1}\cdot\dots\cdot\boldsymbol{\Phi}^{1}\boldsymbol{W}^{1})\in\mathbb{R}^{d},\boldsymbol{x}\in\mathbb{R}^{d}$,
thus the result of (1) is in $\mathbb{R}^{d\times d}$. In (2), we apply the
chain rule, with the first term $\in\mathbb{R}^{d\times m_{l}^{2}}$ and the
second term $\in\mathbb{R}^{m_{l}^{2}\times d}$. In (3) we use the formula
$\frac{\partial\text{vec}(AXB)}{\partial\text{vec}(X)}=B^{\mathrm{T}}\otimes
A$, and the first term $\in\mathbb{R}^{d\times m_{l}}$, second
$\in\mathbb{R}^{1\times m_{l}}$, third $\in\mathbb{R}^{m_{l}^{2}\times d}$. In
(4), we decompose the calculation into dimensional-wise with
$\frac{\partial\text{vec}(\boldsymbol{\Phi}^{l})}{{\partial\boldsymbol{x}_{j}}}\in\mathbb{R}^{m_{l}^{2}}$.
In (5), we use the fact that $B^{\mathrm{T}}\otimes
A\text{vec}(X)=\text{vec}(AXB)$.

###  C.1 Spectral Norm for Complexity

In this subsection, we shall use a covering number approach to Rademacher
complexity. The following lemma is the key to connect them.

######  Lemma C.1.

[25] Let $\mathcal{F}$ be a real-valued function class taking values in
$[0,1]$, and assume that $\mathbf{0}\in\mathcal{F}$. Then

$\text{Rad}\left(\mathcal{F};{S}\right)\leq\inf_{\alpha>0}\left(\frac{4\alpha}{\sqrt{n}}+\frac{12}{n}\int_{\alpha}^{\sqrt{n}}\sqrt{\log\mathcal{N}\left(\mathcal{F}_{S},\varepsilon,\|\cdot\|_{2,2}\right)}d\varepsilon\right),$
(69)

where $S$ is the dataset, and $\mathcal{F}_{S}$ is the set containing the
image of the dataset $S$ under all mappings in $\mathcal{F}$.

We note that this lemma requires that the hypothesis is in the interval
$[0,1]$. In practice, we can consider the class of truncated neural networks.
More specifically, if the class of neural network is denoted $\mathcal{F}$,
then we consider the following class:

$\widehat{\mathcal{F}}=\mathcal{F}_{+}+\mathcal{F}_{-},$ (70)

where $\mathcal{F}_{+}=\left\\{f:f\cap[0,1],f\in\mathcal{F}\right\\}$ and
$\mathcal{F}_{-}=\left\\{f:f\cap[-1,0],f\in\mathcal{F}\right\\}$. Then, the
function class $\widehat{F}$ is bounded in the interval $[-1,1]$, which is
suitable for prediction of the target function $u^{*}({\boldsymbol{x}})$ which
is bounded by 1. In addition, their Rademacher complexity have the
relationship:
$\text{Rad}(\widehat{\mathcal{F}})\leq\text{Rad}(\mathcal{F}_{+})+\text{Rad}(\mathcal{F}_{-})$.
Throughout this paper, we will adopt the truncated neural network function
class unless specified.

######  Definition C.1.

(Matrix Covering) We use $\mathcal{N}(U,\epsilon,\|\cdot\|)$ to denote the
least cardinality of any subset $V\subset U$ that covers $U$ at scale
$\epsilon$ with norm $\|\cdot\|$, i.e.,

$\sup_{A\in U}\min_{B\in V}\|A-B\|\leq\epsilon.$ (71)

######  Lemma C.2.

Consider the Hilbert space $\mathcal{H}$ with the norm
$\|\cdot\|_{\mathcal{H}}$. Let $U\in\mathcal{H}$ has the reresentation
$U=\sum_{i=1}^{d}\alpha_{i}V_{i}$, where all $V_{i}\in\mathcal{H}$, and
$\alpha_{i}>0$. Then for any positive integer $k$, there exists a choice of
non-negative integers $(k_{1},k_{2},\cdots,k_{d})$, such that
$\sum_{i=1}^{d}k_{i}=k$, and

$\|U-\frac{\|\alpha\|_{1}}{k}\sum_{i=1}^{d}k_{i}V_{i}\|^{2}_{\mathcal{H}}\leq\frac{\|\alpha\|_{1}}{k}\sum_{i=1}^{d}\alpha_{i}\|V_{i}\|^{2}_{\mathcal{H}}\leq\frac{\|\alpha\|_{1}^{2}}{k}\max_{1\leq
i\leq d}\|V_{i}\|^{2}_{\mathcal{H}}.$ (72)

###### Proof.

See [33, 25] for details. ∎

######  Lemma C.3.

[25] Let conjugate exponents $(p,q)$ and $(r,s)$ be given with $p\leq 2$, as
well as positive reals $(a,b,\epsilon)$, and positive integer $m$. Let matrix
$X\in\mathbb{R}^{n\times d}$ be given with $\|X\|_{p,p}\leq b$, where $n$ is
the number of training data, and $d$ is the input dimension. Then,

$\log\mathcal{N}\Big{(}\left\\{XA:A\in\mathbb{R}^{d\times
m},\|A^{\mathrm{T}}\|_{q,s}\leq
a\right\\},\epsilon,\|\cdot\|_{2,2}\Big{)}\leq\lceil\frac{a^{2}b^{2}m^{2/r}}{\epsilon^{2}}\rceil\log(2dm),$
(73)

where $m$ can be interpreted as the hidden dimension of the model.

###### Proof.

Fix the dataset matrix $X$, and construct $Y\in\mathbb{R}^{n\times d}$ by
$Y_{:,j}:=X_{:,j}/\|X_{:,j}\|_{p}$. We set the following quantites:

$\displaystyle N$ $\displaystyle:=2dm\in\mathbb{N}^{*},$ (74) $\displaystyle
k$
$\displaystyle:=\lceil\frac{a^{2}b^{2}m^{2/r}}{\epsilon^{2}}\rceil\in\mathbb{N}^{*},$
$\displaystyle\overline{a}$ $\displaystyle:=am^{1/r}\|X\|_{p}.$

Then, we define the following matrix sets:

$\displaystyle\\{V_{1},V_{2},\cdots,V_{N}\\}$
$\displaystyle:=\left\\{gY\boldsymbol{e}_{i}\boldsymbol{e}_{j}^{\mathrm{T}:g\in\\{-1,1\\}},i\in\\{1,2,\cdots,d\\},j\in\\{1,2,\cdots,m\\}\right\\},$
(75) $\displaystyle\mathcal{C}$
$\displaystyle:=\left\\{\frac{\overline{c}}{k}\sum_{i=1}^{N}k_{i}V_{i}:k_{i}\geq
0,\sum_{i=1}^{N}k_{i}=k\right\\}=\left\\{\frac{\overline{c}}{k}\sum_{i=1}^{N}k_{i}V_{i}:k_{i}\geq
0,\sum_{i=1}^{N}k_{i}=k\right\\}.$

We first note that since $p\leq 2$,

$\max_{i}\|V_{i}\|_{2}\leq\max_{i}\|Y\boldsymbol{e}_{i}\|_{2}=\max_{i}\frac{\|X\boldsymbol{e}_{i}\|_{2}}{\|Y\boldsymbol{e}_{i}\|_{p}}\leq
1,$ (76)

where we use the fact that the matrix
$\boldsymbol{e}_{i}\boldsymbol{e}_{j}^{\mathrm{T}}$ has only one non-zero
entry at $(i,j)$.

We then show that $\mathcal{C}$ is the desired cover. Due to the construction
of $\mathcal{C}$, its cardinal number $|\mathcal{C}|\leq N^{k}$, as we can
interpret its construction as choose one from $V_{i}$ at each time and there
are $k$ such choices to form an element in $\mathcal{C}$.

We consider a matrix $\|A\|_{q,s}\leq a$, and construct a covering element in
$\mathcal{C}$ as follows. Let the matrix $\alpha\in\mathbb{R}^{d\times m}$,
whose elements in the $j$-th row are all equal to $\|X_{:,j}\|_{p}$, then
$XA=Y(\alpha\odot A)$.

$\displaystyle\|\alpha\|_{p,r}$
$\displaystyle=\left\|\left(\left\|\alpha_{:,1}\right\|_{p},\ldots,\left\|\alpha_{:,m}\right\|_{p}\right)\right\|_{r}$
(77)
$\displaystyle=\left\|\left(\left\|\left(\left\|X_{:,1}\right\|_{p},\ldots,\left\|X_{:,d}\right\|_{p}\right)\right\|_{p},\ldots,\left\|\left(\left\|X_{:,1}\right\|_{p},\ldots,\left\|X_{:,d}\right\|_{p}\right)\right\|_{p}\right)\right\|_{r}$
$\displaystyle=m^{1/r}\left\|\left(\left\|X_{:,1}\right\|_{p},\ldots,\left\|X_{:,d}\right\|_{p}\right)\right\|_{p}=m^{1/r}\left(\sum_{j=1}^{d}\left\|X_{:,j}\right\|_{p}^{p}\right)^{1/p}$
$\displaystyle=m^{1/r}\left(\sum_{j=1}^{d}\sum_{i=1}^{n}X_{i,j}^{p}\right)^{1/p}=m^{1/r}\|X\|_{p}.$

Define $B:=\alpha\odot A$, whereby using conjugacy of $\|\cdot\|_{p,r}$ and
$\|\cdot\|_{q,s}$ gives

$\|B\|_{1}\leq\langle\alpha,|A|\rangle\leq\|\alpha\|_{p,r}\|A\|_{q,s}\leq
m^{1/r}\|X\|_{p}a=\bar{a}.$ (78)

Consequently, $XA$ is equal to

$YB=Y\sum_{i=1}^{d}\sum_{j=1}^{m}B_{ij}\mathbf{e}_{i}\mathbf{e}_{j}^{\top}=\|B\|_{1}\sum_{i=1}^{d}\sum_{j=1}^{m}\frac{B_{ij}}{\|B\|_{1}}\left(Y\mathbf{e}_{i}\mathbf{e}_{j}^{\top}\right)\in\bar{a}\cdot\operatorname{conv}\left(\left\\{V_{1},\ldots,V_{N}\right\\}\right)$

where $\operatorname{conv}\left(\left\\{V_{1},\ldots,V_{N}\right\\}\right)$ is
the convex hull of $\left\\{V_{1},\ldots,V_{N}\right\\}$.

Combining the preceding constructions with Lemma A.6 there exist nonnegative
integers $\left(k_{1},\ldots,k_{N}\right)$ with $\sum_{i}k_{i}=k$ with

$\left\|XA-\frac{\bar{a}}{k}\sum_{i=1}^{N}k_{i}V_{i}\right\|_{2}^{2}=\left\|YB-\frac{\bar{a}}{k}\sum_{i=1}^{N}k_{i}V_{i}\right\|_{2}^{2}\leq\frac{\bar{a}^{2}}{k}\max_{i}\left\|V_{i}\right\|_{2}^{2}\leq\frac{a^{2}m^{2/r}\|X\|_{p}^{2}}{k}\leq\epsilon^{2}$

The desired cover element is thus
$\frac{a}{k}\sum_{i}k_{i}V_{i}\in\mathcal{C}$. ∎

We first revisite the definition of neural networks. Given weight matrices
$\mathcal{W}=\left(W_{1},\ldots,W_{L}\right)$, where $L$ is the network depth,
we define the mapping $F_{\mathcal{W}}$, and more generally for $i\leq L$
define $\mathcal{W}_{1}^{i}:=\left(W_{1},\ldots,W_{i}\right)$ and

$F_{\mathcal{W}}(Z)=W_{L}\sigma\left(W_{L-1}\cdots\sigma\left(W_{1}Z\right)\cdots\right),$
(79)
$F_{\mathcal{W}_{1}^{i}}(Z):=\sigma\left(W_{i}\sigma\left(W_{i-1}\cdots\sigma\left(W_{1}Z\right)\cdots\right)\right),$
(80)

with the convention $F_{\emptyset}(Z)=Z$, where $Z$ is the input.

Define two sequences of matrix spaces $\mathcal{V}_{1},\ldots,\mathcal{V}_{L}$
and $\mathcal{W}_{2},\ldots,\mathcal{W}_{L+1}$, where $\mathcal{V}_{i}$ has a
norm $|\cdot|_{i}$ and $\mathcal{W}_{i}$ has norm $|||\cdot|||_{i}$.
Specifically, we choose all vector spaces as Euclidean spaces, and choose
$|\cdot|_{i}=\|\cdot\|_{2,2}$ and $|||\cdot|||_{i}=\|\cdot\|_{2,2}$. The
inputs $Z\in\mathcal{V}_{1}=\mathbb{R}^{d\times n}$ satisfy a norm constraint
$\|Z\|_{\infty.\infty}\leq 1$, which means the absolute values of all entries
in all data is not larger than $B$. Specifically, we are using
$Z=X^{\mathrm{T}}$.

The linear operators $A_{i}:\mathcal{V}_{i}\rightarrow\mathcal{W}_{i+1}$ are
associated with some operator norm $|A_{i}|_{i\rightarrow i+1}\leq c_{i}$,
where we use $|A_{i}|_{i\rightarrow i+1}=\|A_{i}\|_{2}$, which is the spectral
norm and satisfies:

$\|A_{i}\|_{2}=\sup_{\|Z\|_{2,2}\leq 1}\|A_{i}Z\|_{2,2}:=c_{i}.$ (81)

For the activation function, we consider the sine activation function, which
is an 1-Lipschitz mappings
$\sigma:\mathcal{W}_{i+1}\rightarrow\mathcal{V}_{i+1}$, having the Lipschitz
constant $\rho_{i}=1$, measured with respect to norms
$|\cdot|_{i}=\|\cdot\|_{\infty}$ and $|||\cdot|||_{i}=\|\cdot\|_{\infty}$. In
other words, for any $z,z^{\prime}\in\mathcal{W}_{i+1}$, we have

$|\sigma(z)-\sigma(z^{\prime})|_{i+1}=\|\sigma(z)-\sigma(z^{\prime})\|_{2,2}\leq\|z-z^{\prime}\|_{2,2}=|||z-z^{\prime}|||_{i+1}.$
(82)

We will prove the following lemma on covering number.

######  Lemma C.4.

Let the resolutions for covering $(\epsilon_{1},\ldots,\epsilon_{L})$ be
given, and given the operator norm bounds $(c_{1},\ldots,c_{L})$. Suppose the
matrices $\mathcal{W}=\left(W_{1},\ldots,W_{L}\right)$ lie within the set
$\mathcal{B}_{1}\times\cdots\times\mathcal{B}_{L}$, where $\mathcal{B}_{i}$
are arbitrary classes with the property that each $W_{i}\in\mathcal{B}_{i}$
has $|W_{i}|_{i\rightarrow i+1}=\|W_{i}\|_{2}\leq c_{i}$. Lastly, let the
dataset $Z=X^{\mathrm{T}}$ be given with $|Z|_{1}=\|Z\|_{\infty,\infty}\leq
1$. Then, letting

$\tau:=\sum_{j=1}^{L}\epsilon_{j}\rho_{j}\prod_{l=j+1}^{L}\rho_{l}c_{l}=\sum_{j=1}^{L}\epsilon_{j}\prod_{l=j+1}^{L}c_{l},$
(83)

then the neural net images
$\mathcal{H}_{Z}:=\left\\{F_{\mathcal{W}}(Z):\mathcal{W}\in\mathcal{B}_{1}\times\cdots\times\mathcal{B}_{L}\right\\}$
have covering number bound:

$\displaystyle\quad\mathcal{N}(\mathcal{H}_{Z},\tau,|\cdot|_{L+1})$ (84)
$\displaystyle=\mathcal{N}(\mathcal{H}_{Z},\tau,\|\cdot\|_{2,2})$
$\displaystyle\leq\prod_{i=1}^{L}\sup_{W_{j}\in\mathcal{B}_{j},\forall
j<i}\mathcal{N}\left(\left\\{W_{i}F_{\mathcal{W}_{1}^{i-1}}(Z):W_{i}\in\mathcal{B}_{i}\right\\},\epsilon_{i},|||\cdot|||_{i+1}\right)$
$\displaystyle=\prod_{i=1}^{L}\sup_{W_{j}\in\mathcal{B}_{j},\forall
j<i}\mathcal{N}\left(\left\\{W_{i}F_{\mathcal{W}_{1}^{i-1}}(Z):W_{i}\in\mathcal{B}_{i}\right\\},\epsilon_{i},\|\cdot\|_{2,2}\right).$

###### Proof.

In this paper, we consider the case when all $c_{1}\geq 1$. We first denote
$\mathcal{F}_{i}$ as the covering set of the image set of all $i$-layer neural
networks, which are inductively constructed as follows. Choose an
$\epsilon_{1}$-cover $\mathcal{F}_{1}$ of the set of one-layer neural networks
$\left\\{W_{1}Z:W_{1}\in\mathcal{B}_{1}\right\\}$, then its cardinality
satisfies

$\left|\mathcal{F}_{1}\right|\leq\mathcal{N}\left(\left\\{W_{1}Z:W_{1}\in\mathcal{B}_{1}\right\\},\epsilon_{1},|||\cdot|||_{2}\right)=:N_{1},$
(85)

by the definition of covering number. For every element $F\in\mathcal{F}_{i}$,
which is a covering vector, we construct an $\epsilon_{i+1}$-cover
$\mathcal{G}_{i+1}(F)$ of the set

$\left\\{W_{i+1}\sigma_{i}(F):W_{i+1}\in\mathcal{B}_{i+1}\right\\},$ (86)

where $F$ is chosen and fixed. Since the covers are proper, i.e.,

$F\in\mathcal{F}_{i}\subset\left\\{W_{i}F_{\mathcal{W}_{1}^{i-1}}(Z):W_{j}\in\mathcal{B}_{j},\forall
j\leq i\right\\},$ (87)

meaning that $F=W_{i}F_{W_{1}^{i-1}}(Z)$ for some matrices
$\left(W_{1},\ldots,W_{i}\right)\in$
$\mathcal{B}_{1}\times\cdots\times\mathcal{B}_{i}$. Then, we obtain

$\left|\mathcal{G}_{i+1}(F)\right|\leq\sup_{\forall j\leq
i.W_{j}\in\mathcal{B}_{j}}\mathcal{N}\left(\left\\{W_{i+1}F_{W_{1},\ldots,W_{i}}(Z):W_{i+1}\in\mathcal{B}_{i+1}\right\\},\epsilon_{i+1},|||\cdot|||_{i+2}\right)=:N_{i+1}.$
(88)

Lastly we construct the cover

$\mathcal{F}_{i+1}:=\bigcup_{F\in\mathcal{F}_{i}}\mathcal{G}_{i+1}(F),$ (89)

whose cardinality satisfies

$\left|\mathcal{F}_{i+1}\right|\leq\left|\mathcal{F}_{i}\right|\cdot
N_{i+1}\leq\prod_{l=1}^{i+1}N_{l}.$ (90)

Define $\mathcal{F}:=\left\\{\sigma(F):F\in\mathcal{F}_{L}\right\\}$. By
construction, $\mathcal{F}$ satisfies the desired cardinality constraint. To
show that it is indeed a cover, fix any $\left(W_{1},\ldots,W_{L}\right)$
satisfying the above constraints, and for convenience define recursively the
mapped elements

$F_{1}=W_{1}X\in\mathcal{W}_{2},\quad
G_{i}=\sigma\left(F_{i}\right)\in\mathcal{V}_{i+1}\quad
F_{i+1}=W_{i+1}G_{i}\in\mathcal{W}_{i+2}.$ (91)

The goal is to show the existence of $\hat{G}_{L}\in\mathcal{F}$ satisfying:

$|G_{L}-\hat{G}_{L}|_{L+1}=\|G_{L}-\hat{G}_{L}\|_{2,2}\leq\tau$ (92)

To this end, inductively construct approximating elements
$\left(\hat{F}_{i},\hat{G}_{i}\right)$ as follows. The Base case: set
$\hat{G}_{0}=X$ since $G_{0}=X$. For other cases, choose
$\hat{F}_{i}\in\mathcal{F}_{i}$ with
$|||W_{i}\hat{G}_{i-1}-\hat{F}_{i}|||_{i+1}\leq\epsilon_{i}$, since
$\mathcal{F}_{i}$ is an $\epsilon_{i}$ cover of the following set:

$\left\\{W_{i}F_{\mathcal{W}_{1}^{i-1}}(Z):W_{j}\in\mathcal{B}_{j},\forall
j\leq i\right\\}.$ (93)

And set $\hat{G}_{i}:=\sigma(\hat{F}_{i})$. To complete the proof, it will be
shown inductively that

$\left|G_{i}-\widehat{G}_{i}\right|_{i+1}\leq\sum_{1\leq j\leq
i}\epsilon_{j}\rho_{j}\prod_{l=j+1}^{i}\rho_{l}c_{l}$ (94)

For the base case,

$\left|G_{0}-\widehat{G}_{0}\right|_{1}=0$ (95)

For the inductive step, we obtain

$\displaystyle\left|G_{i+1}-\widehat{G}_{i+1}\right|_{i+2}$
$\displaystyle\leq\rho_{i+1}|||F_{i+1}-\widehat{F}_{i+1}|||_{i+2}$ (96)
$\displaystyle\leq\rho_{i+1}|||F_{i+1}-W_{i+1}\widehat{G}_{i}|||_{i+2}+\rho_{i+1}|||W_{i+1}\widehat{G}_{i}-\widehat{F}_{i+1}|||_{i+2}$
$\displaystyle\leq\rho_{i+1}\left|W_{i+1}\right|_{i+1\rightarrow
i+2}\left|G_{i}-\widehat{G}_{i}\right|_{i+1}+\rho_{i+1}\epsilon_{i+1}$
$\displaystyle\leq\rho_{i+1}c_{i+1}\left(\sum_{j\leq
i}\epsilon_{j}\rho_{j}\prod_{l=j+1}^{i}\rho_{l}c_{l}\right)+\rho_{i+1}\epsilon_{i+1}$
$\displaystyle=\sum_{j\leq
i+1}\epsilon_{j}\rho_{j}\prod_{l=j+1}^{i+1}\rho_{l}c_{l},$

where we note that $\prod_{l=i+1}^{i}\rho_{l}c_{l}=1$. ∎

The whole-network covering bound in terms of spectral and $(2,1)$ norms now
follows by the general norm covering number and the matrix covering lemma.

######  Theorem C.1.

Let spectral norm bounds $\left(s_{1},\ldots,s_{L}\right)$, and matrix $(2,1)$
norm bounds $\left(b_{1},\ldots,b_{L}\right)$ be given. Let data matrix
$X\in\mathbb{R}^{n\times d}$ be given, where the $n$ rows correspond to data
points. Let $\mathcal{H}_{X}$ denote the family of matrices obtained by
evaluating $X$ with all choices of network $F_{\mathcal{W}}$, i.e.,

$\mathcal{H}_{X}:=\left\\{F_{\mathcal{W}}\left(X^{\mathrm{T}}\right):\mathcal{W}=\left(W_{1},\ldots,W_{L}\right),\left\|W_{i}\right\|_{2}\leq
s_{i},\left\|W_{i}\right\|_{2,1}\leq b_{i}\right\\}.$ (97)

Then for any $\epsilon>0$, we have the covering number bound:

$\log\mathcal{N}\left(\mathcal{H}_{X},\epsilon,\|\cdot\|_{2,2}\right)\leq\frac{nd\log\left(2h^{2}\right)}{\epsilon^{2}}\left(\prod_{j=1}^{L}s_{j}^{2}\right)\left(\sum_{i=1}^{L}\left(\frac{b_{i}}{s_{i}}\right)^{2/3}\right)^{3},$
(98)

where the network width is denoted $h$.

###### Proof.

We set the matrix constraint sets as
$\mathcal{B}_{i}=\left\\{W_{i}:\left\|W_{i}\right\|_{2}\leq
s_{i},\left\|W_{i}\right\|_{2,1}\leq b_{i}\right\\}$, and lastly the per-layer
cover resolutions $\left(\epsilon_{1},\ldots,\epsilon_{L}\right)$ set
according to

$\epsilon_{i}:=\frac{\alpha_{i}\epsilon}{\rho_{i}\prod_{j>i}\rho_{j}s_{j}}=\frac{\alpha_{i}\epsilon}{\prod_{j=i+1}^{L}s_{j}},\quad\text{where}\quad\alpha_{i}:=\frac{1}{\bar{\alpha}}\left(\frac{b_{i}}{s_{i}}\right)^{2/3}\quad,\quad\bar{\alpha}:=\sum_{j=1}^{L}\left(\frac{b_{j}}{s_{j}}\right)^{2/3}.$
(99)

By this choice, it follows that the final cover resolution $\tau$ provided by
Lemma satisfies

$\tau\leq\sum_{j=1}^{L}\epsilon_{j}\prod_{l=j+1}^{L}s_{l}=\sum_{j=1}^{L}\alpha_{j}\epsilon=\epsilon$
(100)

To start, the covering number estimate from Lemma C.4 can be combined with
Lemma C.3 with $p=2,s=1$ to get

$\displaystyle\quad\log\mathcal{N}\left(\mathcal{H}_{X},\epsilon,\|\cdot\|_{2,2}\right)$
(101) $\displaystyle=\sum_{i=1}^{L}\sup_{W_{j}\in\mathcal{B}_{j},\forall
j<i}\log\left(\mathcal{N}\left(\left\\{W_{i}F_{\mathcal{W}_{1}^{i-1}}(Z):W_{i}\in\mathcal{B}_{i}\right\\},\epsilon_{i},\|\cdot\|_{2,2}\right)\right)$
$\displaystyle=\sum_{i=1}^{L}\sup_{W_{j}\in\mathcal{B}_{j},\forall
j<i}\log\left(\mathcal{N}\left(\left\\{F_{\mathcal{W}_{1}^{i-1}}(Z)^{\mathrm{T}}W_{i}^{\mathrm{T}}:W_{i}\in\mathcal{B}_{i}\right\\},\epsilon_{i},\|\cdot\|_{2,2}\right)\right)$
$\displaystyle\leq\sum_{i=1}^{L}\sup_{W_{j}\in\mathcal{B}_{j},\forall
j<i}\frac{b_{i}^{2}\|F_{\mathcal{W}_{1}^{i-1}}(Z)^{\mathrm{T}}\|_{2,2}^{2}\log(2h^{2})}{\epsilon_{i}^{2}}$

We can bound the intermediate outputs of the neural network as follows:

$\displaystyle\left\|F_{\left(W_{1},\ldots,W_{i-1}\right)}\left(Z\right)^{\top}\right\|_{2,2}$
$\displaystyle=\left\|F_{\left(W_{1},\ldots,W_{i-1}\right)}\left(Z\right)\right\|_{2,2}$
(102)
$\displaystyle=\|\sigma\left(W_{i-1}F_{\left(W_{1},\ldots,W_{i-2}\right)}\left(Z\right)\|_{2,2}\right.$
$\displaystyle\leq\left\|W_{i-1}F_{\left(W_{1},\ldots,W_{i-2}\right)}\left(Z\right)\right\|_{2,2}$
$\displaystyle\leq\left\|W_{i-1}\right\|_{2}\left\|F_{\left(W_{1},\ldots,W_{i-2}\right)}\left(Z\right)\right\|_{2,2},$

which by induction gives

$\max_{j}\left\|F_{\left(W_{1},\ldots,W_{i-1}\right)}\left(Z\right)^{\mathrm{T}}\right\|_{2,2}\leq\|Z\|_{2,2}\prod_{j=1}^{i-1}\left\|W_{j}\right\|_{2}.$
(103)

Combining all equations, and then expanding the choice of $\epsilon_{i}$ and
collecting terms, we attain

$\displaystyle\log\mathcal{N}\left(\mathcal{H}_{X},\epsilon,\|\cdot\|_{2}\right)$
$\displaystyle\leq\sum_{i=1}^{L}\sup_{\left(W_{1},\ldots,W_{i-1}\right),j<i,W_{j}\in\mathcal{B}_{j}}\frac{b_{i}^{2}\|X\|_{2,2}^{2}\prod_{j<i}\left\|W_{j}\right\|_{2}^{2}}{\epsilon_{i}^{2}}\log\left(2h^{2}\right)$
(104)
$\displaystyle\leq\|X\|_{2,2}^{2}\sum_{i=1}^{L}\frac{b_{i}^{2}\prod_{j<i}s_{j}^{2}}{\epsilon_{i}^{2}}\log\left(2h^{2}\right)$
$\displaystyle=\|X\|_{2,2}^{2}\frac{\log\left(2h^{2}\right)\prod_{j=1}^{L}s_{j}^{2}}{\epsilon^{2}}\sum_{i=1}^{L}\frac{b_{i}^{2}}{\alpha_{i}^{2}s_{i}^{2}}$
$\displaystyle=\|X\|_{2,2}^{2}\frac{\log\left(2h^{2}\right)\prod_{j=1}^{L}s_{j}^{2}}{\epsilon^{2}}\sum_{i=1}^{L}\frac{\bar{\alpha}^{2}b_{i}^{2}}{s_{i}^{2}}\cdot\left(\frac{s_{i}^{2}}{b_{i}^{2}}\right)^{2/3}$
$\displaystyle=\|X\|_{2,2}^{2}\frac{\log\left(2h^{2}\right)\prod_{j=1}^{L}s_{j}^{2}}{\epsilon^{2}}\bar{\alpha}^{2}\sum_{i=1}^{L}\left(\frac{b_{i}^{2}}{s_{i}^{2}}\right)^{1/3}$
$\displaystyle=\|X\|_{2,2}^{2}\frac{\log\left(2h^{2}\right)\prod_{j=1}^{L}s_{j}^{2}}{\epsilon^{2}}\left(\bar{\alpha}^{3}\right)$
$\displaystyle\leq
nd\|X\|_{\infty,\infty}^{2}\frac{\log\left(2h^{2}\right)\prod_{j=1}^{L}s_{j}^{2}}{\epsilon^{2}}\left(\bar{\alpha}^{3}\right)$

where we have used:

$\epsilon_{i}=\frac{\alpha_{i}\epsilon}{\prod_{j=i+1}^{L}s_{j}},\quad\alpha_{i}=\frac{1}{\bar{\alpha}}\left(\frac{b_{i}}{s_{i}}\right)^{2/3},\quad\bar{\alpha}:=\sum_{j=1}^{L}\left(\frac{b_{j}}{s_{j}}\right)^{2/3}.$
(105)

∎

######  Lemma C.5.

For every $L$, and every set of $n$ points $S\subset\overline{\Omega}$, the
hypothesis class $\mathcal{NN}^{L}_{M,N}$ given by the neural networks

$\mathcal{NN}^{L}_{M,N}:=\left\\{\boldsymbol{x}\mapsto
W_{L}\sigma(W_{L-1}\sigma(\cdots\sigma(W_{1}\boldsymbol{x})))\ |\
\|W_{l}\|_{2}\leq M(l),\frac{\|W_{l}\|_{2,1}}{\|W_{l}\|_{2}}\leq N(l),\forall
l\right\\},$ (106)

satisfies the Rademacher complexity bound

$\text{Rad}(\mathcal{NN}^{L}_{M,N};S)\leq\frac{4}{n\sqrt{n}}+\frac{18\sqrt{d\log(2h^{2})}\log
n}{\sqrt{n}}\prod_{l=1}^{L}M(l)\Big{(}\sum_{l=1}^{L}N(l)^{2/3}\Big{)}^{3/2},$
(107)

where $h$ is the maximal width of the neural network, i.e.,

$h=\max(m_{L},\cdots,m_{0}).$ (108)

###### Proof.

Consider the covering number bound:

$\log\mathcal{N}\left((\mathcal{NN}_{M,N}^{L})_{S},\epsilon,\|\cdot\|_{2}\right)\leq\frac{nd\log\left(2h^{2}\right)}{\epsilon^{2}}\left(\prod_{j=1}^{L}s_{j}^{2}\right)\left(\sum_{i=1}^{L}\left(\frac{b_{i}}{s_{i}}\right)^{2/3}\right)^{3}=:\frac{R}{\epsilon^{2}}.$
(109)

What remains is to relate covering numbers and Rademacher complexity via a
Dudley entropy integral:

$\text{Rad}(\mathcal{NN}^{L}_{M,N};S)\leq\inf_{\alpha>0}\left(\frac{4\alpha}{\sqrt{n}}+\frac{12}{n}\int_{\alpha}^{\sqrt{n}}\sqrt{\frac{R}{\epsilon^{2}}}\mathrm{~{}d}\epsilon\right)=\inf_{\alpha>0}\left(\frac{4\alpha}{\sqrt{n}}+\log(\sqrt{n}/\alpha)\frac{12\sqrt{R}}{n}\right).$
(110)

The inf is uniquely minimized at $\alpha:=3\sqrt{R/n}$, but the desired bound
may be obtained by the simple choice $\alpha:=1/n$, and plugging the resulting
Rademacher complexity estimate:

$\displaystyle\text{Rad}(\mathcal{NN}^{L}_{M,N};S)$
$\displaystyle\leq\frac{4}{n\sqrt{n}}+\log(n^{3/2})\frac{12\sqrt{R}}{n}$ (111)
$\displaystyle=\frac{4}{n\sqrt{n}}+\frac{18\log n}{n}\sqrt{R}$
$\displaystyle\leq\frac{4}{n\sqrt{n}}+\frac{18\sqrt{d\log(2h^{2})}\log
n}{\sqrt{n}}\prod_{l=1}^{L}M(l)\Big{(}\sum_{l=1}^{L}N(l)^{2/3}\Big{)}^{3/2}.$

∎

Recall the expressions of the differentiated PINNs:

$\frac{\partial
u_{\boldsymbol{\theta}}(\boldsymbol{x})}{\partial\boldsymbol{x}}=W_{L}\Phi_{L-1}W_{L-1}\dots\Phi_{1}W_{1}\in\mathbb{R}^{d},$
(112)
$\frac{\partial^{2}u_{\boldsymbol{\theta}}(\boldsymbol{x})}{\partial\boldsymbol{x}^{2}}=\left\\{\sum_{l=1}^{L-1}({W}_{L}{\Phi}_{L-1}\cdots{W}_{l+1})\text{diag}({\Psi}_{l}{W}_{l}\cdots{\Psi}_{1}({W}_{1})_{:,j})({W}_{l}\cdots{\Phi}_{1}{W}_{1})\right\\}_{1\leq
j\leq d}.$ (113)

where

$\Phi_{l}=\text{diag}[\sigma^{\prime}(W_{l}\sigma(W^{l-1}\sigma(\cdots\sigma(W^{1}\boldsymbol{x})\cdots)]\in\mathbb{R}^{m_{l}\times
m_{l}},$ (114)
$\Psi_{l}=\text{diag}[\sigma^{\prime\prime}(W^{l}\sigma(W^{l-1}\sigma(\cdots\sigma(W^{1}\boldsymbol{x})\cdots)]\in\mathbb{R}^{m_{l}\times
m_{l}}.$ (115)

The forwards passes for one input $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$ of the
PINN model are:

$W_{L}\Phi_{L-1}W_{L-1}\dots\Phi_{1}W_{1}b(\boldsymbol{x}_{i})\in\mathbb{R},$
(116)
$\sum_{j=1}^{d}\sum_{l=1}^{L-1}{W}_{L}{\Phi}_{L-1}\cdots{W}_{l+1}\text{diag}({\Psi}_{l}{W}_{l}\cdots{\Psi}_{1}({W}_{1})_{:,j}){W}_{l}\cdots{\Phi}_{1}{W}_{1}A_{:,j}\in\mathbb{R}.$
(117)

Now we begin to prove their Rademacher complexities.

######  Lemma C.6.

The hypothesis class $\mathcal{F}_{2}$ given by

$\mathcal{F}_{2}=\left\\{\boldsymbol{x}\mapsto
W_{L}\Phi_{L-1}W_{L-1}\dots\Phi_{1}W_{1}b(\boldsymbol{x})\ |\
\|W_{l}\|_{2}\leq M(l),\frac{\|W_{l}\|_{2,1}}{\|W_{l}\|_{2}}\leq N(l),\forall
l\right\\},$ (118)

satisfies the Rademacher complexity bound

$\text{Rad}(\mathcal{F}_{2};S)\leq\frac{4}{n\sqrt{n}}+\frac{18L\sqrt{2d\log(2h^{2})}\log
n}{\sqrt{n}}\left(\prod_{l=1}^{L}M(l)\right)^{2}\left(\sum_{l=1}^{L}N(l)^{2/3}\right)^{3/2}.$
(119)

###### Proof.

We only need to consider how to cover the set of PINNs. We first consider the
first-order derivatives, i.e., how to cover:

$\boldsymbol{x}\in\mathbb{R}^{d}\mapsto
W_{L}\Phi_{L-1}W_{L-1}\dots\Phi_{1}W_{1}b(\boldsymbol{x})\in\mathbb{R}.$ (120)

To do so, we consider the covering in Lemma C.4. Specifically, denote
$\mathcal{F}_{i}$ as the $\epsilon_{i}$ covering of the following set:

$\left\\{W_{i}F_{\mathcal{W}_{1}^{i-1}}(Z):W_{j}\in\mathcal{B}_{j},\forall
j\leq i\right\\},$ (121)

which is the set of all $i$-layer neural networks, with weight matrices in the
set $\mathcal{B}_{j}$. In Lemma C.4, we have inductively constructed covers
$\mathcal{F}_{1},\ldots,\mathcal{F}_{L}$. Specifically, $\mathcal{F}_{1}$ is
an $\epsilon_{1}$-cover of the set of one-layer neural networks
$\left\\{W_{1}Z:W_{1}\in\mathcal{B}_{1}\right\\}$, and its cardinality
satisfies

$\left|\mathcal{F}_{1}\right|\leq\mathcal{N}\left(\left\\{W_{1}Z:W_{1}\in\mathcal{B}_{1}\right\\},\epsilon_{1},\|\cdot\|_{2,2}\right).$
(122)

For every element $F\in\mathcal{F}_{i}$, which is a covering vector of the set
of image of all $i$-layer neural networks, we construct an
$\epsilon_{i+1}$-cover $\mathcal{G}_{i+1}(F)$ of the set