standard global,type,entry[boolean]subentry[true]

# Generalization in Graph Neural Networks: Improved PAC-Bayesian Bounds on
Graph Diffusion

Haotian Ju  Dongyue Li  Aneesh Sharma  Hongyang R. Zhang***The paper is
presented at Artificial Intelligence and Statistics (AISTATS), 2023. Email
correspondence: {ju.h, li.dongyu<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
###### Abstract

Graph neural networks are widely used tools for graph prediction tasks.
Motivated by their empirical performance, prior works have developed
generalization bounds for graph neural networks, which scale with graph
structures in terms of the maximum degree. In this paper, we present
generalization bounds that instead scale with the largest singular value of
the graph neural network’s feature diffusion matrix. These bounds are
numerically much smaller than prior bounds for real-world graphs. We also
construct a lower bound of the generalization gap that matches our upper bound
asymptotically. To achieve these results, we analyze a unified model that
includes prior works’ settings (i.e., convolutional and message-passing
networks) and new settings (i.e., graph isomorphism networks). Our key idea is
to measure the stability of graph neural networks against noise perturbations
using Hessians. Empirically, we find that Hessian-based measurements correlate
with the observed generalization gaps of graph neural networks accurately.
Optimizing noise stability properties for fine-tuning pretrained graph neural
networks also improves test performance on several graph-level classification
tasks.

##  1 Introduction

A central measure of success for a machine learning model is the ability to
generalize well from training data to test data. For linear and shallow
models, the generalization gap between their training performance and test
performance can be quantified via complexity notions such as the
Vapnik–Chervonenkis dimension and Rademacher complexity. However, formally
explaining the empirical generalization performance of deep models remains a
challenging problem and an active research area (see, e.g., recent textbook by
[116]). There are by now many studies for fully-connected and convolutional
neural networks that provide an explanation for their superior empirical
performance [72, 95]. Our work seeks to formally understand generalization in
graph neural networks (GNN) [63], which are commonly used for learning on
graphs [75].

As a concrete example for motivating the study of generalization performance,
we consider the fine-tuning of pretrained graph neural networks [105]. Given a
pretrained GNN learned on a diverse range of graphs, fine-tuning the
pretrained GNN on a specific prediction task is a common approach for transfer
learning on graphs. An empirical problem with fine-tuning is that, on one
hand, pretrained GNNs use lots of parameters to ensure representational power.
On the other hand, fine-tuning a large GNN would overfit the training data and
suffer poor test performance without proper algorithmic intervention. Thus, a
better understanding of generalization in graph neural networks can help us
identify the cause of overfitting and, consequently, inspire designing robust
fine-tuning methods for graph neural networks.

A naive application of the generalization bounds from fully-connected
feedforward networks [72, 84] to GNNs would imply an extra term in the
generalization bound that scales with $n^{l-1}$, where $n$ is the number of
nodes in the graph, hence rendering the error bounds vacuous. Besides, [85]
shows that the VC dimension of GNN scales with $n$. Thus, although the VC
dimension is a classical notion for deriving learning bounds, it is oblivious
to the graph structure. Recent works have taken a step towards addressing this
issue with better error analysis. [98] (verma2019stability) find that one-
layer graph neural networks satisfy uniform stability properties [98],
following the work of [71]. The generalization bound of [98] scales with the
largest singular value of the graph diffusion matrix of the model. However,
their analysis only applies to a single layer and node prediction. [104]
(garg2020generalization) analyze an $l$ layer message-passing neural network —
with $l-1$ graph diffusion layers and $1$ pooling layer — for graph prediction
tasks [104]. Their result scales with $d^{l-1}$, where $d$ is the maximum
degree of the graph. Subsequently, [119] (liao2020pac) develop a tighter bound
but still scales with $d^{l-1}$ [119]. For both results, the graph’s maximum
degree is used to quantify the complexity of node aggregation in each
diffusion step.

(a) Spectral norm vs. max degree for five graphs

(b) Comparing generalization bounds

Figure 1: The spectral norm bounds for graph diffusion matrices are orders of
magnitude smaller than maximum degree bounds on real-world graphs; In Figure
1(a), we measure the spectral norm and the max degree for five graph datasets.
In Figure 1(b), the Hessian-based generalization measure (plotted in green,
scaled according to the _left axis_) matches the empirically observed
generalization gaps of graph neural networks (plotted in yellow, scaled
according to the _left axis_). The blue line shows a uniform-convergence bound
(scaled according to the _right axis_) that is orders of magnitude larger than
the observed gaps.

Our main contribution is to show generalization bounds for graph neural
networks by reducing the max degree to the spectral norm of the graph
diffusion matrix. We achieve this by analyzing the stability of a graph neural
network against noise injections. To illustrate, denote an $l$-layer GNN by
$f$. By PAC-Bayesian analysis [65], the generalization gap of $f$ will be
small if $f$ remains stable against noise injections; otherwise, the
generalization gap of $f$ will be large. Empirically, quantifying the noise
stability of $f$ via Lipschitz-continuity properties of its activation
functions leads to nonvacuous bounds for feedforward networks that correlate
with their observed generalization gaps [79, 123, 73, 106]. Our theoretical
analysis rigorously formalizes this with refined stability analysis of graph
neural networks through Hessians, leading to tight generalization bounds on
the graph diffusion matrix.

Our Contribution. The goal of this work is to improve the theoretical
understanding of generalization in graph neural networks, and in that vein, we
highlight two results below:

  * •

First, we prove sharp generalization bounds for message-passing neural
networks [70, 74, 82], graph convolutional networks [77], and graph
isomorphism networks [99]. Our bounds scale with the spectral norm of
$P_{{}_{G}}^{l-1}$ for an $l$-layer network, where $P_{{}_{G}}$ denotes a
diffusion matrix on a graph $G$ and varies between different models (see
Theorem 3.1 for the full statement). We then show a matching lower bound
instance where the generalization gap scales with the spectral norm of
$P_{{}_{G}}^{l-1}$ (see Theorem 3.2).

  * •

Second, our stability analysis of graph neural networks provides a practical
tool for measuring generalization. Namely, we show that the trace of the loss
Hessian matrix can measure the noise stability of GNN. The formal statement is
given in Lemma 4.3, and our techniques, which include a uniform convergence of
the Hessian matrix, may be of independent interest. We note that the proof
applies to twice-differentiable and Lipschitz-continuous activations (e.g.,
tanh and sigmoid).

Taken together, these two results provide a sharp understanding of
generalization in terms of the graph diffusion matrix for graph neural
networks. We note that the numerical value of our bounds in their dependence
on the graph is much smaller than prior results [104, 119], as is clear from
Figure 1(a). Moreover, the same trend holds even after taking weight norms
into account (see Figure 2, Section 3.3). Further, the Hessian-based bounds
(see Lemma 4.3, Section 4) are non-vacuous, matching the scale of empirically
observed generalization gaps in Figure 1(b).

Finally, motivated by the above analysis, we also present an algorithm that
performs gradient updates on perturbed weight matrices of a graph neural
network. The key insight is that minimizing the average loss of multiple
perturbed models with independent noise injections is equivalent to
regularizing $f$’s Hessian in expectation. We conduct experiments on several
graph classification tasks with Molecular graphs that show the benefit of this
algorithm in the fine-tuning setting.

##  2 Related Work

Generalization Bounds: An article by [78] (zhang2016understanding) finds that
deep nets have enough parameters to memorize real images with random labels,
yet they still generalize well if trained with true labels. This article
highlights the overparametrized nature of modern deep nets (see also a recent
article by [110]), motivating the need for complexity measures beyond
classical notions. In the case of two-layer ReLU networks, [95]
(neyshabur2018towards) show that (path) norm bounds better capture the
“effective number of parameters” than VC dimension—which is the number of
parameters for piecewise linear activations [90].

For multilayer networks, subsequent works have developed norm, and margin
bounds, either via Rademacher complexities [72, 80, 107], or PAC-Bayesian
bounds [84, 79, 118, 123]. All of these bounds apply to the fine-tuning
setting following the distance from the initialization perspective. Our
analysis approach builds on the work of [79] (arora2018stronger) and [123]
(ju2022robust). The latter work connects perturbed losses and Hessians for
feedforward neural networks, with one limitation Hessians do not show any
explicit dependence on the data. This is a critical issue for GNN as we need
to incorporate the graph structure in the generalization bound. Our result
instead shows an explicit dependence on the graph and applies to message-
passing layers that involve additional nonlinear mappings. We will compare our
analysis approach and prior analysis in more detail when we present the proofs
in Section 4 (see Remark 4.4).

Graph Representation Learning: Most contemporary studies of learning on graphs
consider either node-level or graph-level prediction tasks. Our result applies
to graph prediction while permitting an extension to node prediction: see
Remark 4.2 in Section 4. Most graph neural networks follow an information
diffusion mechanism on graphs [63]. Early work takes inspiration from ConvNets
and designs local convolution on graphs, e.g., spectral networks [66], GCN
[77], and GraphSAGE [76] (among others). Subsequent works have designed new
architectures with graph attention [86] and isomorphism testing [99]. [74]
(gilmer2017neural) synthesize several models into a framework called message-
passing neural networks. Besides, one could also leverage graph structure in
the pooling layer (e.g., differentiable pooling and hierarchical pooling [89,
88]). It is conceivable that one can incorporate the model complexity of these
approaches into our analysis. Recent work applies pretraining to large-scale
graph datasets for learning graph representations [105]. Despite being an
effective transfer learning approach, few works have examined the
generalization of graph neural nets in the fine-tuning step.

Besides learning on graphs, GNNs are also used for combinatorial optimization
[97] and causal reasoning [109]. There is another approach for graph
prediction using graph kernels [64] such as the graph neural tangent kernel
[91]. Lastly, we remark that graph diffusion processes have been studied in
earlier literature on social and information networks [67, 68, 100]. For
further references about different applications of GNN, see review articles
[75, 103, 108, 101].

Generalization in Graph Neural Networks: Recent work explores generalization
by formalizing the role of the algorithm and the alignment between networks
and tasks [120]. [113] finds that transductive Rademacher complexity-based
bound provides insights into the behavior of GNNs in the stochastic block
model. Besides, there are works about size generalization, which refer to
performance degradation when models extrapolate to graphs of different sizes
from the input [97, 121]. It is conceivable that the new tools we have
developed may be useful for studying extrapolation.

Expressivity of Graph Neural Networks: The expressivity of GNN for graph
classification can be related to graph isomorphism tests and has connections
to one-dimensional Weisfeiler-Lehman testing of graph isomorphism [94, 99].
This implies limitations of GNN for expressing tasks such as counting cycles
[96, 102, 111]. The expressiveness view seems orthogonal to generalization,
which instead concerns the sample efficiency of learning. For further
discussions and references, see a recent survey by [122].

##  3 Sharp Generalization Bounds for Graph Neural Networks

We first introduce the problem setup for analyzing graph neural networks.
Then, we state our generalization bounds for graph neural networks and compare
them with the prior art. Lastly, we construct an example to argue that our
bounds are tight.

###  3.1 Problem setup

Consider a graph-level prediction task. Suppose we have $N$ examples in the
training set; each example is an independent sample from a distribution
denoted as $\mathcal{D}$, which is jointly supported on the feature space
$\mathcal{X}$ times the label space $\mathcal{Y}$. In each example, we have an
undirected graph denoted as $G=(V,E)$, which describes the connection between
$n$ entities, represented by nodes in $V$. For example, a node could represent
a molecule, and an edge between two nodes is a bond between two molecules.
Each node also has a list of $d$ features. Denote all node features as an $n$
by $d$ matrix $X$. For graph-level prediction tasks, the goal is to predict a
graph label $y$ for every example. We will describe a few examples of such
tasks later in Section 5.2.

Message-passing neural networks (MPNN). We study a model based on several
prior works for graph-level prediction tasks [70, 74, 104, 119]. Let $l$ be
the number of layers: the first $l-1$ layers are diffusion layers, and the
last layer is a pooling layer. Let $d_{t}$ denote the width of each layer for
$t$ from $1$ up to $l$. There are several nonlinear mappings in layer $t$,
denoted as $\phi_{t},\rho_{t}$, and $\psi_{t}$; further, they are all centered
at zero. There is a weight matrix $W^{(t)}$ of dimension $d_{t-1}$ by $d_{t}$
for transforming neighboring features, and another weight matrix $U^{(t)}$ of
dimension $d$ by $d_{t}$ for transforming the anchor node feature.

For the first $l-1$ layers, we recursively compute the node embedding from the
input features $H^{(0)}=X$:

$\displaystyle
H^{(t)}=\phi_{t}\Big{(}{XU^{(t)}+\rho_{t}\big{(}{P_{{}_{G}}\psi_{t}(H^{(t-1)})}\big{)}W^{(t)}}\Big{)},\text{
for }t=1,2,\dots,l.$ (1)

For the last layer $l$, we aggregate the embedding of all nodes: let
$\bm{1}_{n}$ be a vector with $n$ values of one:

$\displaystyle H^{(l)}$
$\displaystyle=\frac{1}{n}\bm{1}_{n}^{\top}H^{(l-1)}W^{(l)}.$ (2)

Note that this setting subsumes many existing GNNs. Several common designs for
the graph diffusion matrix $P_{{}_{G}}$ would be the adjacency matrix of the
graph (denoted as $A$). $P_{{}_{G}}$ can also be the normalized adjacency
matrix, $D^{-1}A$, with $D$ being the degree-diagonal matrix. Adding an
identity matrix in $A$ is equivalent to adding self-loops in $G$. For GCN, we
set $U^{(t)}$ as zero, $\rho_{t}$ and $\psi_{t}$ as identity mappings.

Notations. For any matrix $X$, let $\left\|X\right\|$ denote the largest
singular value (or spectral norm) of $X$. Let $\left\|X\right\|_{{F}}$ denote
the Frobenius norm of $X$. We use the notation $f(N)\lesssim g(N)$ to indicate
that there exists a fixed constant $c$ that does not grow with $N$ such that
$f(N)\leq c\cdot g(N)$ for large enough values of $N$. Let $\bm{W}$ and
$\bm{U}$ denote the union of the $W$ and $U$ matrices in a model $f$,
respectively.

###  3.2 Main results

Given a message-passing neural network denoted as $f$, what can we say about
its generalization performance, i.e., the gap between population and empirical
risks? Let $f(X,G)$ denote the output of $f$, given input with graph $G$, node
feature matrix $X$, and label $y$. The loss of $f$ for this input example is
denoted as $\ell(f(X,G),y)$. Let $\hat{{L}}(f)$ denote the empirical loss of
$f$ over the training set. Let ${L}(f)$ denote the expected loss of $f$ over a
random drawn of $\mathcal{D}$. We are interested in the generalization gap of
$f$, i.e., ${L}(f)-\hat{{L}}(f)$. How would the graph diffusion matrix
$P_{{}_{G}}$ affect the generalization gap of graph neural networks?

To motivate our result, we examine the effect of incorporating graph diffusion
in a one-layer linear neural network. That is, we consider $f(X,G)$ to be
$\frac{1}{n}\bm{1}_{n}^{\top}P_{{}_{G}}XW^{(1)}$, which does not involve any
nonlinear mapping for simplicity of our discussion. In this case, by standard
spectral norm inequalities for matrices, the Euclidean norm of $f$ (which is a
vector) satisfies:

$\displaystyle\left\|f(X,G)\right\|=$
$\displaystyle\left\|\frac{1}{n}\bm{1}_{n}^{\top}P_{{}_{G}}XW^{(1)}\right\|$
$\displaystyle\leq$
$\displaystyle\left\|\frac{1}{n}\bm{1}_{n}^{\top}\right\|\cdot\left\|P_{{}_{G}}\right\|\cdot\left\|X\right\|\cdot\left\|W^{(1)}\right\|$
(3)

Thus, provided that the loss function $\ell(\cdot,y)$ is Lipschitz-continuous,
standard arguments imply that the generalization gap of $f$ scales with the
spectral norm of ${P_{{}_{G}}}$ (divided by ${\sqrt{N}}$) [83]. Let us compare
this statement with a fully-connected neural net that averages the node
features, i.e., the graph diffusion matrix $P_{{}_{G}}$ is the identity
matrix. The spectral norm of $P_{{}_{G}}$ becomes one. Together, we conclude
that the graph structure affects the generalization bound of a single layer
GNN by adding the spectral norm of ${P_{{}_{G}}}$.

Our main result is that incorporating the spectral norm of the _graph
diffusion matrix_ $P_{{}_{G}}^{l-1}$ is sufficient for any $l$ layer MPNN. We
note that the dependence is a power of $l-1$ because there are $l-1$ graph
diffusion layers: see equation (1). Let $f$ be an $l$-layer network whose
weights $\bm{W},\bm{U}$ are defined within a hypothesis set $\mathcal{H}$: For
every layer $i$ from $1$ up to $l$, we have that

$\displaystyle\left\|W^{(i)}\right\|\leq$ $\displaystyle
s_{i},~{}~{}\left\|W^{(i)}\right\|_{{F}}\leq s_{i}r_{i},$
$\displaystyle\left\|U^{(i)}\right\|\leq$ $\displaystyle
s_{i},~{}~{}\left\|U^{(i)}\right\|_{{F}}\,\,\leq s_{i}r_{i},$ (4)

where $s_{1},s_{2},\dots,s_{l}$ and $r_{1},r_{2},\dots,r_{l}$ are bounds on
the spectral norm and stable rank and are all greater than or equal to one,
without loss of generality. We now present the full statement.

######  Theorem 3.1.

Suppose all of the nonlinear activations in
$\\{\phi_{t},\rho_{t},\psi_{t}:\forall\,t\\}$ and the loss function
$\ell(\cdot,y)$ (for any fixed label $y\in\mathcal{Y}$) are twice-
differentiable, Lipschitz-continuous and their first-order and second-order
derivatives are both Lipschitz-continuous.

With probability at least $1-\delta$ over the randomness of $N$ independent
samples from $\mathcal{D}$, for any $\delta>0$, and any $\epsilon>0$ close to
zero, any model $f$ with weight matrices in the set $\mathcal{H}$ satisfies:

$\displaystyle{L}(f)\leq(1+\epsilon)\hat{{L}}(f)+{\sum_{i=1}^{l}\sqrt{\frac{CBd_{i}\Big{(}{\max\limits_{(X,G,y)\sim\mathcal{D}}\left\|X\right\|^{2}\left\|P_{G}\right\|^{2(l-1)}}\Big{)}\Big{(}{r_{i}^{2}\prod\limits_{j=1}^{l}s_{j}^{2}}\Big{)}}{N}}}+\mathcal{\textup{O}}\left(\frac{\log(\delta^{-1})}{N^{3/4}}\right),$
(5)

where $B$ is an upper bound on the value of the loss function $\ell(x,y)$ for
any $(x,y)\sim\mathcal{D}$, $C$ is a fixed constant depending on the
activation, and the loss function (see Eq. (44), Appendix A.2.4).

As a remark, prior works by [104] and [119] consider an MPNN with $W^{(t)}$
and $U^{(t)}$ being the same for $t$ from $1$ up to $l$, motivated by
practical designs [74, 82]. Thus, their analysis is conducted separately for
GCN and MPNN with weight tying. By contrast, our result allows $W^{(t)}$ and
$U^{(t)}$ to be arbitrarily different across different layers. This unifies
GCN and MPNN without weight tying in the same framework so that we can unify
their analysis. We defer the proof sketch of our result and a discussion to
Section 4.

###  3.3 Comparison with prior art

Table 1: How does the generalization gap of graph neural networks scale with graph properties? In this work, we show spectrally-normalized bounds on $P_{{}_{G}}$ and compare our results with prior results in the following table. We let $A$ denote the adjacency matrix, $D$ be the degree-diagonal matrix of $A$, and $l$ be the depth of the GNN. Previous generalization bounds scale with the graph’s maximum degree denoted as $d$. Our result instead scales with the spectral norm of $P_{{}_{G}}$ and applies to graph isomorphism networks (GIN) [99] and GraphSAGE with mean aggregation [76].  Graph Dependence | GCN | MPNN | GIN | GraphSAGE-Mean
---|---|---|---|---
[104] (garg2020generalization) | $d^{l-1}$ | $d^{l-1}$ | - | -
[119] (liao2020pac) | $d^{\frac{l-1}{2}}$ | $d^{l-1}$ | - | -
Ours (Theorems 3.1 and 4.5) | $1$ | $\left\|A\right\|^{l-1}$ | $\sum_{i=1}^{l-1}\frac{\left\|A\right\|^{i}}{l-1}$ | $\left\|D^{-1}A\right\|^{l-1}$

In Table 1, we compare our result with prior results. We first illustrate the
effects of graph properties on the generalization bounds. Then we will also
show a numerical comparison to incorporate the other components of the bounds.

  * •

Suppose $P_{{}_{G}}$ is the adjacency matrix of $G$. Then, one can show that
for any undirected graph $G$, the spectral norm of ${P_{{}_{G}}}$ is less than
the maximum degree $d$ (cf. Fact A.1, Appendix A for a proof). This explains
why our result is strictly less than prior results for MPNN in Table 1.

  * •

Suppose $P_{{}_{G}}$ is the normalized and symmetric adjacency matrix of $G$:
$P_{{}_{G}}=\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}$, where $\tilde{A}$ is
$A+\operatorname{Id}$ and $\tilde{D}$ is the degree-diagonal matrix of
$\tilde{A}$. Then, the spectral norm of $P_{{}_{G}}$ is at most one (cf. Fact
A.1, Appendix A for a proof). This fact explains why the graph dependence of
our result for GCN is $1$ in Table 1. Thus, we can see that this provides an
exponential improvement compared to the prior results.

Thus, for the above diffusion matrices, we conclude that the spectral norm of
${P_{{}_{G}}}$ is strictly smaller than the maximum degree of graph $G$
(across all graphs in the distribution $\mathcal{D}$).

Next, we conduct an empirical analysis to compare our results and prior
results numerically. Following the setting of prior works, we use two types of
models that share their weight matrices across different layers, including GCN
[77] and the MPNN specified in [119]. For both models, we evaluate the
generalization bounds by varying the network depth $l$ between $2,4$, and $6$.

We consider graph prediction tasks on three collaboration networks, including
IMDB-B, IMDB-M, and COLLAB [69]. IMDB-B includes a collection of movie
collaboration graphs. In each graph, a node represents an actor or an actress,
and an edge denotes a collaboration in the same movie. The task is to classify
each graph into the movie genre as Action or Romance. The IMDB-M is a multi-
class extension with the movie graph label Comedy, Romance, or Sci-Fi. COLLAB
includes a list of ego-networks of scientific researchers. Each graph includes
a researcher and her collaborators as nodes. An edge in the graph indicates a
collaboration between two researchers. The task is to classify each ego-
network into the field of the researcher, including High Energy, Condensed
Matter, and Astro Physics.

We report the numerical comparison in Figure 2, averaged over three random
seeds. Our results are consistently smaller than previous results. As
explained in Table 1, the improvement comes from the spectral norm bounds on
graphs compared with the max degree bounds.

(a) Two-layer GCN

(b) Four-layer GCN

(c) Six-layer GCN

(d) Two-layer MPNN

(e) Four-layer MPNN

(f) Six-layer MPNN

Figure 2: Comparing our result and prior results [104, 119] on three graph
classification tasks conducted on GCNs and MPNNs, respectively.

###  3.4 A matching lower bound

Next, we show an instance with the same dependence on the graph diffusion
matrix as our upper bound. In our example:

  * •

The graph is the complete graph with self-loops inserted in each node. Thus,
the adjacency matrix of $G$ is a square matrix with all ones. We will set
$P_{{}_{G}}$ as the adjacency matrix of $G$.

  * •

In the first $l-1$ graph diffusion layers, the activation functions
$\phi,\rho,\psi$ are all linear functions. Further, we fix all the parameters
of $\bm{U}$ as zero.

  * •

The loss function $\ell$ is the logistic loss.

Then, we demonstrate a data distribution such that there always exists some
weight matrices within $\mathcal{H}$ whose generalization gap must increase in
proportion to the spectral norm of $P_{{}_{G}}^{l-1}$ and the product of the
spectral norm of every layer $s_{1},s_{2},\dots,s_{l}$.

######  Theorem 3.2.

Let $N_{0}$ be a sufficiently large value. For any norms
$s_{1},s_{2},\dots,s_{n}$, there exists a data distribution $\mathcal{D}$ on
which with probability at least $0.1$ over the randomness of $N$ independent
samples from $\mathcal{D}$, for any $N\geq N_{0}$, the generalization gap of
$f$ is greater than the following:

$\displaystyle{\left|{{L}(f)-\hat{L}(f)}\right|}\gtrsim\sqrt{\frac{\Big{(}{\max\limits_{(X,G,y)\sim\mathcal{D}}\left\|P_{{}_{G}}\right\|^{2(l-1)}}\Big{)}\Big{(}{\prod\limits_{i=1}^{l}s_{i}^{2}}\Big{)}}{N}}.$
(6)

Notice that the lower bound in (6) exhibits the same scaling in terms of
$G$—$\left\|P_{{}_{G}}\right\|^{l-1}$—as our upper bound from equation (5).
Therefore, we conclude that our spectral norm bound is tight for multilayer
MPNN. The proof of the lower bound can be found in Appendix A.3.

######  Remark 3.3.

Our results from Theorem 3.1 and 3.2 together suggest the generalization error
bound scales linearly in $l$. To verify whether this is the case, we conducted
an empirical study on three architectures (GCN, GIN-Mean, and GIN-Sum) that
measured the growth of generalization errors as the network depth $l$ varies.
We find that the generalization error grows sublinearly with $l$ to
$\left\|P_{{}_{G}}\right\|$. We also note that this sublinear growth trend has
been captured by our Hessian-based generalization bound (cf. Figure 1(a)). It
would be interesting to understand better why the sublinear trend happens and
further provide insight into the behavior of GNN.

######  Remark 3.4.

Theorem 3.2 suggests that in the worst case, the generalization bound would
have to scale with the spectral norms of the graph and the weight matrices.
Although this is vacuous for large $l$, later in Lemma 4.1, we show a data-
dependent bound using the trace of the Hessians, which is non-vacuous. As
shown in Figure 1(a), Hessian-based measurements match the scale of actual
generalization errors: the green line, calculated based on the trace of the
loss Hessian matrix (cf. equation (7)), matches the scale of actual
generalization gaps plotted in the yellow line.

##  4 Proof Techniques and Extensions

Our analysis for dealing with the graph structure seems fundamentally
different from the existing analysis. In the margin analysis of [119], the
authors also incorporate the graph structure in the perturbation error. For
bounding the perturbation error, the authors use a triangle inequality that
results in a $(1,\infty)$ norm of the matrix ${P_{{}_{G}}}$ (see Lemma 3.1 of
[119] for GCN). We note that this norm can be larger than the spectral norm by
a factor of $\sqrt{n}$, where $n$ is the number of nodes in $G$: in the case
of a star graph, this norm for the graph diffusion matrix of GCN is
$\sqrt{n}$. By comparison, the spectral norm of the same matrix is less than
one (see Fact A.1, Appendix A).

How can we tighten the perturbation error analysis and the dependence on
$P_{{}_{G}}$ in the generalization bounds, then? Our proof involves two parts:

  * •

Part I: By expanding the perturbed loss of a GNN, we prove a bound on the
generalization gap using the trace of the Hessian matrix associated with the
loss.

  * •

Part II: Then, we explicitly bound the trace of the Hessian matrix with the
spectral norm of the graph using the Lipschitzness of the activation
functions.

Part I: Measuring noise stability using the Hessian. We first state an
implicit generalization bound that measures the trace of the Hessian matrix.
Let $\mathbf{H}^{(i)}$ denote the Hessian matrix of the loss $\ell(f(X,G),y)$
with respect to layer $i$’s parameters, for each $i$ from $1$ up to $l$.
Particularly, $\mathbf{H}^{(i)}$ is a square matrix whose dimension depends on
the number of variables within layer $i$. Let $\mathbf{H}$ denote the Hessian
matrix of the loss $\ell(f(X,G),y)$ over all parameters of $f$.

######  Lemma 4.1.

In the setting of Theorem 3.1, with probability at least $1-\delta$ over the
randomness of the $N$ training examples, for any $\delta>0$ and $\epsilon$
close to $0$, we get:

$\displaystyle{L}(f)\leq(1+\epsilon)\hat{{L}}(f)+(1+\epsilon)\sum_{i=1}^{l}\sqrt{\frac{B\cdot\left(\max\limits_{(X,G,y)\sim\mathcal{D}}\operatorname{Tr}\big{[}\mathbf{H}^{(i)}[\ell(f(X,G),y)]\big{]}\right)s_{i}^{2}r_{i}^{2}}{N}}+\mathcal{\textup{O}}\left(\frac{\log(\delta^{-1})}{N^{3/4}}\right).$
(7)

##### Proof Sketch.

At a high level, the above result follows from Taylor’s expansion of the
perturbed loss. Suppose each parameter of $f$ is perturbed by an independent
noise drawn from a Gaussian distribution with mean zero and variance
$\sigma^{2}$. Let $\tilde{\ell}(f(X,G),y)$ be the perturbed loss value of an
input example $X,G$ with label $y$. Let $\mathcal{E}$ denote the noise
injections organized in a vector. Using Taylor’s expansion of the perturbed
loss $\tilde{\ell}$, we get:

$\displaystyle\tilde{\ell}(f(X,G),y)-\ell(f(X,G),y)=\mathcal{E}^{\top}\nabla\ell(f(X,G),y)+\frac{1}{2}{\mathcal{E}}^{\top}\mathbf{H}\big{[}\ell(f(X,G),y)\big{]}{\mathcal{E}}+\mathcal{\textup{O}}(\sigma^{3}).$
(8)

Notice that the expectation of the first-order expansion term above equals
zero. The expectation of the second-order expansion term becomes $\sigma^{2}$
times the trace of the loss Hessian. To derive equation (7), we use a PAC-
Bayes bound of [65, Theorem 2 ]. There are two parts to this PAC-Bayes bound:

  * •

The expectation of the noise perturbation in equation (8), taken over the
injected noise $\mathcal{E}$;

  * •

The KL divergence between the prior and the posterior, which is at most
$s_{i}^{2}r_{i}^{2}$ for layer $i$, for $i$ from $1$ up to $l$, within the
hypothesis set $\mathcal{H}$.

Thus, one can balance the two parts by adjusting the noise variance at each
layer—this leads to the layerwise Hessian decomposition in equation (7).

A critical step is showing the uniform convergence of the Hessian matrix. We
achieve this based on the Lipschitz-continuity of the first and twice
derivatives of the nonlinear activation mappings. With these conditions, we
prove the uniform convergence with a standard $\epsilon$-cover argument. The
complete proof can be found in Appendix A.2.1.

######  Remark 4.2.

Our argument in Lemma 4.1 applies to graph-level prediction tasks, which
assume an unknown distribution of graphs. A natural question is whether the
analysis applies to node-level prediction tasks, which are often treated as
semi-supervised learning problems. The issue with directly applying our
analysis to semi-supervised learning is that the size of a graph is only
finite. Instead, a natural extension would be to think about our graph as a
random sample from some population and then argue about generalization in
expectation of the random sample. It is conceivable that one can prove a
similar spectral norm bound for node prediction in this extension. This would
be an interesting question for future work.

Part II: Spectral norm bounds of the trace of the Hessian. Next, we explicitly
analyze the trace of the Hessian at each layer. We bound the trace of the
Hessian using the spectral norm of the weight matrices and the graph based on
the Lipschitz-continuity conditions from Theorem 3.1. Notice that the last
layer is a linear pooling layer, which can be deduced from layer $l-1$. Hence,
we consider the first $l-1$ layers below.

######  Lemma 4.3.

In the setting of Theorem 3.1, the trace of the loss Hessian matrix
$\mathbf{H}^{(i)}$ taken over $W^{(i)}$ and $U^{(i)}$ satisfies the following,
for any $i=1,2,\cdots,l-1$,

$\displaystyle{\left|{\operatorname{Tr}\left[\mathbf{H}^{(i)}\big{[}{\ell\big{(}{f(X,G),y}\big{)}}\big{]}\right]}\right|}$
$\displaystyle\lesssim~{}$ $\displaystyle
s_{l}^{2}\Bigg{(}\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}\right\|_{{F}}+\sum_{p=1}^{d_{0}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}+\left\|\frac{\partial
H^{(l-1)}}{\partial W^{(i)}}\right\|_{{F}}^{2}+\left\|\frac{\partial
H^{(l-1)}}{\partial U^{(i)}}\right\|_{{F}}^{2}\Bigg{)}$ (9)
$\displaystyle\lesssim~{}$
$\displaystyle{\left\|X\right\|^{2}\left\|P_{{}_{G}}\right\|^{2(l-1)}}{\prod\nolimits_{j=1:\,j\neq
i}^{l}s_{j}^{2}}.$ (10)

##### Proof Sketch.

Equation (9) uses the chain rule to expand out the trace of the Hessian and
then applies the Lipschitzness of the loss function. Based on this result,
equation (10) then bounds the first and second derivatives of $H^{(l-1)}$.
This step is achieved via an induction of $\partial H^{(j)}$ and
$\partial^{2}H^{(j)}$ over $W^{(i)}$ and $U^{(i)}$, for $j=1,\dots,l-1$ and
$i=1,\dots,j$. The induction relies on the feedforward architecture and the
Lipschitzness of the first and second derivatives. We leave out a few details,
such as the constants in equations (9) and (10) that can be found in Appendix
A.2.2 and A.2.3. Combining both parts together, we get equation (1).

######  Remark 4.4.

We compare our analysis with the approach of [119]. Both our analysis and
[119] follow the PAC-Bayesian framework. But additionally, we explore
Lipschitz-continuity properties of the first and second derivatives of the
activation functions (e.g., examples of such activations include tanh and
sigmoid). This allows us to measure the perturbation loss with Hessians, which
captures data-dependent properties much more accurately than the margin
analysis of [119]. It would be interesting to understand if one could still
achieve spectral norm bounds on graphs under weaker smoothness conditions.
This is left for future work.

###  4.1 Extensions

Graph isomorphism networks. This architecture concatenates every layer’s
embedding together for more expressiveness [99]. A classification layer is
used after the layers. Let $V^{(i)}$ denote a $d_{i}$ by $k$ matrix (recall
$k$ is the output dimension). Denote the set of these matrices by
$\mathcal{V}$. We average the loss of all of the classification layers. Let
$\hat{L}_{{}_{GIN}}(f)$ denote the average loss of $f$ over $N$ independent
samples of $\mathcal{D}$. Let ${{L}}_{{}_{GIN}}(f)$ denote the expected loss
of $f$ over a random sample of $\mathcal{D}$. See also equation (45) in
Appendix A.4 for their precise definitions.

Next, we state a generalization bound for graph isomorphism networks. Let $f$
be any $l$-layer MPNN with weights defined in a hypothesis space
$\mathcal{H}$: the parameters of $f$ reside within the constraints from
equation (4); further, for every $i$ from $1$ up to $l$, the spectral norm of
$V^{(i)}$ is less than $s_{l}$. Building on Lemma 4.3, we show a bound that
scales with the spectral norm generalization of the averaged graph diffusion
matrices. Let $P_{{}_{GIN}}$ denote the average of $l-1$ matrices:
$P_{{}_{G}},P_{{}_{G}}^{2},\dots,P_{{}_{G}}^{l-1}$. We state the result below.

######  Corollary 4.5.

Suppose the nonlinear activation mappings and the loss function satisfy the
conditions stated in Theorem 3.1. With probability at least $1-\delta$ for any
$\delta\geq 0$, and any $\epsilon$ close to zero, any $f$ in $\mathcal{H}$
satisfies:

$\displaystyle{L}_{{}_{GIN}}(f)\leq(1+\epsilon)\hat{L}_{{}_{GIN}}(f)+{\sum_{i=1}^{l}\sqrt{\frac{{BCd_{i}\cdot\left({\max\limits_{(X,G,y)\sim\mathcal{D}}\left\|X\right\|^{2}\left\|P_{{}_{GIN}}\right\|^{2}}\right)\Big{(}{r_{i}^{2}\prod\limits_{j=1}^{l}s_{j}^{2}}\Big{)}}}{N}}}+\mathcal{\textup{O}}\left(\frac{\log(\delta^{-1})}{N^{3/4}}\right),$
(11)

where $B$ is an upper bound on the value of the loss function $\ell$ across
the data distribution $\mathcal{D}$, $C$ is a fixed constant that only depends
on the Lipschitz-continuity of the activation mappings and the loss.

The proof can be found in Appendix A.4. In particular, we apply the trace norm
bound over the model output of every layer. The classification layer, which
only uses a linear transformation, can also be incorporated.

Fine-tuning Graph Neural Networks. We note that all of our bounds can be
applied to the fine-tuning setting, where a graph neural network is
initialized with pretrained weights and then fine-tuned on the target task.
The results can be extended to this setting by setting the norm bounds within
equation (4) as the distance between the pretrained and fine-tuned model.

##  5 An Algorithm for Fine-tuning Graph Neural Networks

A common practice to apply deep networks on graphs is to fit pretrained GNNs,
which can be fine-tuned on a target task. Typically, only a small amount of
data is available for fine-tuning. Thus, the fine-tuned model may overfit the
training data, incurring a large generalization gap. A central insight from
our analysis is that maintaining a small perturbed loss ensures lower
generalization gaps. Motivated by this observation, we design an algorithm to
minimize the perturbed risk of GNN.

Let $f$ denote a GNN and $\tilde{\ell}(f)$ be the perturbed loss of $f$, with
noise injected inside $f$’s weight matrices. Recall from step (8) that
$\tilde{\ell}(f)$ is equal to $\ell(f)$ plus several expansion terms. In
particular, minimizing the expectation of $\tilde{\ell}(f)$ is equivalent to
minimizing $\hat{L}(f)$ plus the trace of the Hessian matrix. To estimate this
expectation, we sample several noise perturbations independently. Because
Taylor’s expansion of $\tilde{\ell}(f)$ also involves the gradient, we cancel
this out by computing the perturbed loss with the negated perturbation.
Algorithm 1 describes the complete procedure.

We evaluate the above algorithm for fine-tuning pretrained GNNs. Empirical
results reveal that this algorithm achieves better test performance compared
with existing regularization methods for five graph classification tasks. The
code repository for reproducing the experiments is at
https://github.com/VirtuosoResearch/Generalization-in-graph-neural-networks.

###  5.1 Experimental setup

We focus on graph classification tasks, including five datasets from the
MoleculeNet benchmark [87]. Each dataset aims to predict whether a molecule
has a certain chemical property given its graph representation. We use
pretrained GINs from [105] and fine-tune the model on each downstream task.
Following their experimental setup, we use the scaffold split for the dataset,
and the model architecture is fixed for all five datasets. Each model has $5$
layers; each layer has $300$ hidden units and uses average pooling in the
readout layer. We set the parameters, such as the learning rate and the number
of epochs following their setup.

Algorithm 1 Algorithms for fine-tuning graph neural networks

Input: A training dataset $\\{(X_{i},G_{i},y_{i})\\}_{i=1}^{N}$ with node
feature $X_{i}$, graph $G_{i}$, and graph-level label $y_{i}$, for
$i=1,\dots,N$.

Require: Number of perturbations $m$, noise variance $\sigma^{2}$, learning
rate $\eta$, and number of epochs $T$.

Output: A trained model $f^{(T)}$.

1: At $t=0$, initialize the parameters of $f^{(0)}$ with pretrained GNN weight
matrices.

2: for $1\leq t\leq T$ do

3: for $1\leq i\leq m$ do

4: Add perturbation $\mathcal{E}_{i}$ drawn from a normal distribution with
mean zero and variance $\sigma^{2}$.

5: Let $\tilde{L}_{i}(f^{(t-1)})$ be the training loss of the model
$f^{(t-1)}$ with weight matrix perturbed by $\mathcal{E}_{i}$.

6: Let $\tilde{L}_{i}^{{}^{\prime}}(f^{(t-1)})$ be the training loss of the
model $f^{(t-1)}$ with weight matrix perturbed by $-\mathcal{E}_{i}$.

7: end for

8: Use stochastic gradient descent to update $f^{(t)}$ as
$f^{(t-1)}-\frac{\eta}{2m}\sum_{i=1}^{m}\big{(}\nabla\tilde{{L}}_{i}\big{(}{f^{(t-1)}}\big{)}+\nabla\tilde{{L}}_{i}^{{}^{\prime}}\big{(}{f^{(t-1)}}\big{)}\big{)}$.

9: end for

We compare our algorithm with previous regularization methods that serve as
benchmark approaches for improving generalization. This includes early
stopping, weight decay, dropout, weight averaging [81], and distance-based
regularization [115]. For implementing our algorithm, we set the number of
perturbations as $10$ and choose the noise standard deviation $\sigma$ with a
grid search in $\\{0.01,0.02,0.05,0.1,0.2,0.5\\}$.

###  5.2 Experimental results

Table 2 reports the test ROC-AUC performance averaged over multiple binary
prediction tasks in each dataset. Comparing the average ranks of methods
across datasets, our algorithm outperforms baselines on all five molecular
property prediction datasets. The results support our theoretical analysis
that the noise stability property of GNN is a strong measure of empirical
generalization performance. Next, we provide details insights from applying
our algorithm.

First, we hypothesize that our algorithm is particularly effective when the
empirical generalization gap is large. To test the hypothesis, we vary the
size of the training set in the BACE dataset; we compare the performance of
our algorithm with early stopping until epoch $100$. We plot the
generalization gap between the training and test losses during training, shown
in Figure 3(a)-3(b). As the trend shows, our algorithm consistently reduces
the generalization gap, particularly when the training set size $N$ is $600$.

Second, we hypothesize that our algorithm helps reduce the trace of the
Hessian matrix (associated with the loss). We validate this by plotting the
trace of the Hessian as the number of epochs progresses during training, again
using the BACE dataset as an example. Specifically, we average the trace over
the training dataset. Figure 3(c) shows the averaged trace values during the
fine-tuning process. The results confirm that noise stability optimization
reduces the trace of the Hessian matrix (more significantly than early
stopping). We note that noise stability optimization also reduces the largest
eigenvalue of the Hessian matrix, along with reducing the trace. This can be
seen in Figure 3(d).

Lastly, we study the number of perturbations used in our algorithm. While more
perturbations would lead to a better estimation of the noisy stability, we
observe that using $10$ perturbations is sufficient for getting the most gain.
We also validate that using negated perturbations consistently performs better
than not using them across five datasets. This is because the negated
perturbation cancels out the first-order term in Taylor’s expansion. In our
ablation study, we find that adding the negated perturbation performs better
than not using it by 1% on average over the five datasets.

######  Remark 5.1.

We note that noise stability optimization is closely related to sharpness-
aware minimization (SAM) [114]. Noise stability optimization differs in two
aspects compared with SAM. First, SAM requires solving constrained minimax
optimization, which may not even be differentiable [112]. Our objective
remains the same after perturbation. Second, SAM reduces the largest
eigenvalue of the Hessian matrix, which can be seen from Taylor’s expansion of
$\tilde{\ell}(f)$. We reduce the trace of the Hessian matrix, thus reducing
the largest eigenvalue. There is another related work that regularizes noise
stability in NLP [117]. Their approach adds noise perturbation to the input
and regularizes the loss change in the output. Our approach adds perturbation
to weight matrices.

Table 2: Test ROC-AUC (%) score for five molecular property prediction datasets with different regularization methods. The reported results are averaged over five random seeds. Dataset | SIDER | ClinTox | BACE | BBBP | Tox21
---|---|---|---|---|---
# Molecule Graphs | 1,427 | 1,478 | 1,513 | 2,039 | 7,831
# Binary Prediction Tasks | 27 | 2 | 1 | 1 | 12
Early Stopping |  61.06$\pm$1.48 |  68.25$\pm$2.63 |  82.86$\pm$0.95 |  67.80$\pm$1.05 |  77.52$\pm$0.23
Weight Decay |  61.30$\pm$0.21 |  67.43$\pm$2.88 |  83.72$\pm$0.99 |  67.98$\pm$2.41 |  78.23$\pm$0.35
Dropout |  63.90$\pm$0.90 |  73.70$\pm$2.80 |  84.50$\pm$0.70 |  68.07$\pm$1.30 |  78.30$\pm$0.30
Weight Averaging |  63.67$\pm$0.34 |  78.78$\pm$1.49 |  83.93$\pm$0.36 |  70.26$\pm$0.24 |  77.59$\pm$0.11
Distance-based Reg. |  64.36$\pm$0.48 |  76.68$\pm$1.19 |  84.65$\pm$0.48 |  70.37$\pm$0.44 |  78.62$\pm$0.24
Ours (Alg. 1) | 65.13$\pm$0.18 | 80.18$\pm$0.82 | 85.07$\pm$0.43 | 71.22$\pm$0.36 | 79.31$\pm$0.24
Figure 3: In Figures 3(a) and 3(b), we show that our algorithm is particularly
effective at reducing the generalization gap for small training dataset sizes
$N$. In Figures 3(c) and 3(d), we find that both the trace and the largest
eigenvalue of the loss Hessian matrix decreased during training.

(a) $N=1200$

(b) $N=600$

(c) Trace

(d) Largest Eigenvalue

##  6 Conclusion

This work develops generalization bounds for graph neural networks with a
sharp dependence on the graph diffusion matrix. The results are achieved
within a unified setting that significantly extends prior works. In
particular, we answer an open question mentioned in [119]: a refined PAC-
Bayesian analysis can improve the generalization bounds for message-passing
neural networks. These bounds are obtained by analyzing the trace of the
Hessian matrix with the Lipschitz-continuity of the activation functions.
Empirical findings suggest that the Hessian-based bound matches observed gaps
on real-world graphs. Thus, our work also develops a practical tool to measure
the generalization performance of graph neural networks. The algorithmic
results with noise stability optimization further demonstrate the practical
implication of our findings.

Our work opens up many interesting questions for future work. Could the new
tools we have developed be used to study generalization in graph attention
networks [86]? Could Hessians be used for measuring out-of-distribution
generalization gaps of graph neural networks? We remark that in a subsequent
paper [124], we provide a rigorous analysis of the convergence rates of our
algorithm presented in the experiments, showing matching upper and lower
bounds given Lipschitz-continuous conditions of the gradients.

## Acknowledgement

Thanks to Renjie Liao, Haoyu He, and the anonymous referees for providing
constructive feedback on our work. Thanks to Yang Yuan for the helpful
discussions. H. Z. is grateful to Stefanie Jegelka for an invitation to
present this paper at her lab meeting. H. J. and D. L. acknowledge financial
support from the startup fund of Khoury College of Computer Sciences,
Northeastern University. [sorting=nyt]

## References

  * [1] Sanjeev Arora “Technical perspective: Why don’t today’s deep nets overfit to their training data?” In _Communications of the ACM_ , 2021
  * [2] Sanjeev Arora, Rong Ge, Behnam Neyshabur and Yi Zhang “Stronger generalization bounds for deep nets via a compression approach” In _ICML_ , 2018
  * [3] Waiss Azizian and Marc Lelarge “Expressive power of invariant and equivariant graph neural networks” In _ICLR_ , 2021
  * [4] Peter Bartlett, Dylan J Foster and Matus Telgarsky “Spectrally-normalized margin bounds for neural networks” In _NeurIPS_ , 2017
  * [5] Peter Bartlett, Nick Harvey, Christopher Liaw and Abbas Mehrabian “Nearly-tight VC-dimension and pseudodimension bounds for piecewise linear neural networks” In _JMLR_ , 2019
  * [6] Joan Bruna, Wojciech Zaremba, Arthur Szlam and Yann LeCun “Spectral networks and locally connected networks on graphs” In _ICLR_ , 2014
  * [7] Ines Chami, Sami Abu-El-Haija, Bryan Perozzi, Christopher Ré and Kevin Murphy “Machine learning on graphs: A model and comprehensive taxonomy” In _arXiv preprint arXiv:2005.03675_ , 2020, pp. 1
  * [8] Zhengdao Chen, Lei Chen, Soledad Villar and Joan Bruna “Can graph neural networks count substructures?” In _NeurIPS_ , 2020
  * [9] Hanjun Dai, Bo Dai and Le Song “Discriminative embeddings of latent variable models for structured data” In _ICML_ , 2016
  * [10] Constantinos Daskalakis, Stratis Skoulakis and Manolis Zampetakis “The complexity of constrained min-max optimization” In _STOC_ , 2021
  * [11] Simon S Du, Kangcheng Hou, Russ R Salakhutdinov, Barnabas Poczos, Ruosong Wang and Keyulu Xu “Graph neural tangent kernel: Fusing graph neural networks with graph kernels” In _NeurIPS_ , 2019
  * [12] Gintare Karolina Dziugaite and Daniel M Roy “Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data” In _UAI_ , 2017
  * [13] Federico Errica, Marco Podda, Davide Bacciu and Alessio Micheli “A fair comparison of graph neural networks for graph classification” In _ICLR_ , 2020
  * [14] Pascal Esser, Leena Chennuru Vankadara and Debarghya Ghoshdastidar “Learning theory can (sometimes) explain generalisation in graph neural networks” In _NeurIPS_ , 2021
  * [15] Pierre Foret, Ariel Kleiner, Hossein Mobahi and Behnam Neyshabur “Sharpness-aware minimization for efficiently improving generalization” In _ICLR_ , 2021
  * [16] Vikas Garg, Stefanie Jegelka and Tommi Jaakkola “Generalization and representational limits of graph neural networks” In _ICML_ , 2020
  * [17] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals and George E Dahl “Neural message passing for quantum chemistry” In _ICML_ , 2017
  * [18] Ashish Goel, Sanjeev Khanna, Sharath Raghvendra and Hongyang Zhang “Connectivity in random forests and credit networks” In _SODA_ , 2014
  * [19] Ashish Goel, Kamesh Munagala, Aneesh Sharma and Hongyang Zhang “A note on modeling retweet cascades on Twitter” In _Workshop on Algorithms and Models for the Web-Graph, 2015_ , 2015
  * [20] Noah Golowich, Alexander Rakhlin and Ohad Shamir “Size-independent sample complexity of neural networks” In _COLT_ , 2018
  * [21] Henry Gouk, Timothy M Hospedales and Massimiliano Pontil “Distance-Based Regularisation of Deep Networks for Fine-Tuning” In _ICLR_ , 2021
  * [22] Benjamin Guedj “A Primer on PAC-Bayesian Learning” In _Proceedings of the French Mathematical Society_ , 2019
  * [23] Will Hamilton, Rex Ying and Jure Leskovec “Representation learning on graphs: Methods and applications” In _arXiv preprint arXiv:1709.05584_ , 2017
  * [24] Will Hamilton, Zhitao Ying and Jure Leskovec “Inductive representation learning on large graphs” In _NeurIPS_ , 2017
  * [25] Moritz Hardt and Benjamin Recht “Patterns, predictions, and actions: A story about machine learning” In _arXiv preprint arXiv:2102.05242_ , 2021
  * [26] Moritz Hardt, Benjamin Recht and Yoram Singer “Train faster, generalize better: Stability of stochastic gradient descent” In _ICML_ , 2016
  * [27] Weihua Hu, Bowen Liu, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay Pande and Jure Leskovec “Strategies for pre-training graph neural networks” In _ICLR_ , 2020
  * [28] Hang Hua, Xingjian Li, Dejing Dou, Cheng-Zhong Xu and Jiebo Luo “Noise stability regularization for improving BERT fine-tuning” In _ACL_ , 2021
  * [29] Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov and Andrew Gordon Wilson “Averaging weights leads to wider optima and better generalization” In _UAI_ , 2018
  * [30] Stefanie Jegelka “Theory of Graph Neural Networks: Representation and Learning” In _arXiv preprint arXiv:2204.07697_ , 2022
  * [31] Yiding Jiang, Behnam Neyshabur, Hossein Mobahi, Dilip Krishnan and Samy Bengio “Fantastic generalization measures and where to find them” In _ICLR_ , 2020
  * [32] Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M Kakade and Michael I Jordan “A short note on concentration inequalities for random vectors with subgaussian norm” In _arXiv preprint arXiv:1902.03736_ , 2019
  * [33] Wengong Jin, Regina Barzilay and Tommi Jaakkola “Junction tree variational autoencoder for molecular graph generation” In _ICML_ , 2018
  * [34] Haotian Ju, Dongyue Li and Hongyang R Zhang “Robust Fine-Tuning of Deep Neural Networks with Hessian-based Generalization Guarantees” In _ICML_ , 2022
  * [35] Haotian Ju, Dongyue Li and Hongyang R Zhang “Noise Stability Optimization for Flat Minima with Tight Rates” In _arXiv preprint arXiv:2306.08553_ , 2023
  * [36] Thomas N Kipf and Max Welling “Semi-supervised classification with graph convolutional networks” In _ICLR_ , 2017
  * [37] Dongyue Li and Hongyang R Zhang “Improved regularization and robustness for fine-tuning in neural networks” In _NeurIPS_ , 2021
  * [38] Renjie Liao, Raquel Urtasun and Richard Zemel “A PAC-Bayesian Approach to Generalization Bounds for Graph Neural Networks” In _ICLR_ , 2021
  * [39] Philip M Long and Hanie Sedghi “Generalization bounds for deep convolutional neural networks” In _ICLR_ , 2020
  * [40] David McAllester “A PAC-Bayesian tutorial with a dropout bound” In _arXiv preprint arXiv:1307.2118_ , 2013
  * [41] Mehryar Mohri, Afshin Rostamizadeh and Ameet Talwalkar “Foundations of machine learning” MIT press, 2018
  * [42] Christopher Morris, Martin Ritzert, Matthias Fey, William L Hamilton, Jan Eric Lenssen, Gaurav Rattan and Martin Grohe “Weisfeiler and leman go neural: Higher-order graph neural networks” In _AAAI_ , 2019
  * [43] Behnam Neyshabur, Srinadh Bhojanapalli and Nathan Srebro “A pac-bayesian approach to spectrally-normalized margin bounds for neural networks” In _ICLR_ , 2018
  * [44] Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun and Nathan Srebro “Towards understanding the role of over-parametrization in generalization of neural networks” In _ICLR_ , 2019
  * [45] Ryoma Sato, Makoto Yamada and Hisashi Kashima “Approximation ratios of graph neural networks for combinatorial problems” In _NeurIPS_ , 2019
  * [46] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner and Gabriele Monfardini “The graph neural network model” In _IEEE transactions on neural networks_ , 2008
  * [47] Franco Scarselli, Ah Chung Tsoi and Markus Hagenbuchner “The vapnik–chervonenkis dimension of graph and recursive neural networks” In _Neural Networks_ , 2018
  * [48] Daniel Selsam, Matthew Lamm, Benedikt Bünz, Percy Liang, Leonardo Moura and David L Dill “Learning a SAT solver from single-bit supervision” In _ICLR_ , 2019
  * [49] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio and Yoshua Bengio “Graph attention networks” In _ICLR_ , 2018
  * [50] Saurabh Verma and Zhi-Li Zhang “Stability and generalization of graph convolutional neural networks” In _KDD_ , 2019
  * [51] S. Vishwanathan, Nicol Schraudolph, Risi Kondor and Karsten Borgwardt “Graph kernels” In _JMLR_ , 2010
  * [52] Zhenqin Wu, Bharath Ramsundar, Evan N Feinberg, Joseph Gomes, Caleb Geniesse, Aneesh S Pappu, Karl Leswing and Vijay Pande “MoleculeNet: a benchmark for molecular machine learning” In _Chemical science_ , 2018
  * [53] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang and Yu Philip “A comprehensive survey on graph neural networks” In _IEEE transactions on neural networks and learning systems_ , 2020
  * [54] Keyulu Xu, Weihua Hu, Jure Leskovec and Stefanie Jegelka “How powerful are graph neural networks?” In _ICLR_ , 2019
  * [55] Keyulu Xu, Jingling Li, Mozhi Zhang, Simon S Du, Ken-ichi Kawarabayashi and Stefanie Jegelka “What can neural networks reason about?” In _ICLR_ , 2020
  * [56] Keyulu Xu, Mozhi Zhang, Jingling Li, Simon S Du, Ken-ichi Kawarabayashi and Stefanie Jegelka “How neural networks extrapolate: From feedforward to graph neural networks” In _ICLR_ , 2021
  * [57] Pinar Yanardag and S. Vishwanathan “Deep graph kernels” In _KDD_ , 2015
  * [58] Gilad Yehudai, Ethan Fetaya, Eli Meirom, Gal Chechik and Haggai Maron “From local structures to size generalization in graph neural networks” In _ICML_ , 2021
  * [59] Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton and Jure Leskovec “Hierarchical graph representation learning with differentiable pooling” In _NeurIPS_ , 2018
  * [60] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht and Oriol Vinyals “Understanding deep learning requires rethinking generalization” In _ICLR_ , 2017
  * [61] Hongyang Zhang, Huacheng Yu and Ashish Goel “Pruning based distance sketches with provable guarantees on random graphs” In _WWW_ , 2019
  * [62] Muhan Zhang, Zhicheng Cui, Marion Neumann and Yixin Chen “An end-to-end deep learning architecture for graph classification” In _AAAI_ , 2018

## References

  * [63] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner and Gabriele Monfardini “The graph neural network model” In _IEEE transactions on neural networks_ , 2008
  * [64] S. Vishwanathan, Nicol Schraudolph, Risi Kondor and Karsten Borgwardt “Graph kernels” In _JMLR_ , 2010
  * [65] David McAllester “A PAC-Bayesian tutorial with a dropout bound” In _arXiv preprint arXiv:1307.2118_ , 2013
  * [66] Joan Bruna, Wojciech Zaremba, Arthur Szlam and Yann LeCun “Spectral networks and locally connected networks on graphs” In _ICLR_ , 2014
  * [67] Ashish Goel, Sanjeev Khanna, Sharath Raghvendra and Hongyang Zhang “Connectivity in random forests and credit networks” In _SODA_ , 2014
  * [68] Ashish Goel, Kamesh Munagala, Aneesh Sharma and Hongyang Zhang “A note on modeling retweet cascades on Twitter” In _Workshop on Algorithms and Models for the Web-Graph, 2015_ , 2015
  * [69] Pinar Yanardag and S. Vishwanathan “Deep graph kernels” In _KDD_ , 2015
  * [70] Hanjun Dai, Bo Dai and Le Song “Discriminative embeddings of latent variable models for structured data” In _ICML_ , 2016
  * [71] Moritz Hardt, Benjamin Recht and Yoram Singer “Train faster, generalize better: Stability of stochastic gradient descent” In _ICML_ , 2016
  * [72] Peter Bartlett, Dylan J Foster and Matus Telgarsky “Spectrally-normalized margin bounds for neural networks” In _NeurIPS_ , 2017
  * [73] Gintare Karolina Dziugaite and Daniel M Roy “Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data” In _UAI_ , 2017
  * [74] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals and George E Dahl “Neural message passing for quantum chemistry” In _ICML_ , 2017
  * [75] Will Hamilton, Rex Ying and Jure Leskovec “Representation learning on graphs: Methods and applications” In _arXiv preprint arXiv:1709.05584_ , 2017
  * [76] Will Hamilton, Zhitao Ying and Jure Leskovec “Inductive representation learning on large graphs” In _NeurIPS_ , 2017
  * [77] Thomas N Kipf and Max Welling “Semi-supervised classification with graph convolutional networks” In _ICLR_ , 2017
  * [78] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht and Oriol Vinyals “Understanding deep learning requires rethinking generalization” In _ICLR_ , 2017
  * [79] Sanjeev Arora, Rong Ge, Behnam Neyshabur and Yi Zhang “Stronger generalization bounds for deep nets via a compression approach” In _ICML_ , 2018
  * [80] Noah Golowich, Alexander Rakhlin and Ohad Shamir “Size-independent sample complexity of neural networks” In _COLT_ , 2018
  * [81] Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov and Andrew Gordon Wilson “Averaging weights leads to wider optima and better generalization” In _UAI_ , 2018
  * [82] Wengong Jin, Regina Barzilay and Tommi Jaakkola “Junction tree variational autoencoder for molecular graph generation” In _ICML_ , 2018
  * [83] Mehryar Mohri, Afshin Rostamizadeh and Ameet Talwalkar “Foundations of machine learning” MIT press, 2018
  * [84] Behnam Neyshabur, Srinadh Bhojanapalli and Nathan Srebro “A pac-bayesian approach to spectrally-normalized margin bounds for neural networks” In _ICLR_ , 2018
  * [85] Franco Scarselli, Ah Chung Tsoi and Markus Hagenbuchner “The vapnik–chervonenkis dimension of graph and recursive neural networks” In _Neural Networks_ , 2018
  * [86] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio and Yoshua Bengio “Graph attention networks” In _ICLR_ , 2018
  * [87] Zhenqin Wu, Bharath Ramsundar, Evan N Feinberg, Joseph Gomes, Caleb Geniesse, Aneesh S Pappu, Karl Leswing and Vijay Pande “MoleculeNet: a benchmark for molecular machine learning” In _Chemical science_ , 2018
  * [88] Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton and Jure Leskovec “Hierarchical graph representation learning with differentiable pooling” In _NeurIPS_ , 2018
  * [89] Muhan Zhang, Zhicheng Cui, Marion Neumann and Yixin Chen “An end-to-end deep learning architecture for graph classification” In _AAAI_ , 2018
  * [90] Peter Bartlett, Nick Harvey, Christopher Liaw and Abbas Mehrabian “Nearly-tight VC-dimension and pseudodimension bounds for piecewise linear neural networks” In _JMLR_ , 2019
  * [91] Simon S Du, Kangcheng Hou, Russ R Salakhutdinov, Barnabas Poczos, Ruosong Wang and Keyulu Xu “Graph neural tangent kernel: Fusing graph neural networks with graph kernels” In _NeurIPS_ , 2019
  * [92] Benjamin Guedj “A Primer on PAC-Bayesian Learning” In _Proceedings of the French Mathematical Society_ , 2019
  * [93] Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M Kakade and Michael I Jordan “A short note on concentration inequalities for random vectors with subgaussian norm” In _arXiv preprint arXiv:1902.03736_ , 2019
  * [94] Christopher Morris, Martin Ritzert, Matthias Fey, William L Hamilton, Jan Eric Lenssen, Gaurav Rattan and Martin Grohe “Weisfeiler and leman go neural: Higher-order graph neural networks” In _AAAI_ , 2019
  * [95] Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun and Nathan Srebro “Towards understanding the role of over-parametrization in generalization of neural networks” In _ICLR_ , 2019
  * [96] Ryoma Sato, Makoto Yamada and Hisashi Kashima “Approximation ratios of graph neural networks for combinatorial problems” In _NeurIPS_ , 2019
  * [97] Daniel Selsam, Matthew Lamm, Benedikt Bünz, Percy Liang, Leonardo Moura and David L Dill “Learning a SAT solver from single-bit supervision” In _ICLR_ , 2019
  * [98] Saurabh Verma and Zhi-Li Zhang “Stability and generalization of graph convolutional neural networks” In _KDD_ , 2019
  * [99] Keyulu Xu, Weihua Hu, Jure Leskovec and Stefanie Jegelka “How powerful are graph neural networks?” In _ICLR_ , 2019
  * [100] Hongyang Zhang, Huacheng Yu and Ashish Goel “Pruning based distance sketches with provable guarantees on random graphs” In _WWW_ , 2019
  * [101] Ines Chami, Sami Abu-El-Haija, Bryan Perozzi, Christopher Ré and Kevin Murphy “Machine learning on graphs: A model and comprehensive taxonomy” In _arXiv preprint arXiv:2005.03675_ , 2020, pp. 1
  * [102] Zhengdao Chen, Lei Chen, Soledad Villar and Joan Bruna “Can graph neural networks count substructures?” In _NeurIPS_ , 2020
  * [103] Federico Errica, Marco Podda, Davide Bacciu and Alessio Micheli “A fair comparison of graph neural networks for graph classification” In _ICLR_ , 2020
  * [104] Vikas Garg, Stefanie Jegelka and Tommi Jaakkola “Generalization and representational limits of graph neural networks” In _ICML_ , 2020
  * [105] Weihua Hu, Bowen Liu, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay Pande and Jure Leskovec “Strategies for pre-training graph neural networks” In _ICLR_ , 2020
  * [106] Yiding Jiang, Behnam Neyshabur, Hossein Mobahi, Dilip Krishnan and Samy Bengio “Fantastic generalization measures and where to find them” In _ICLR_ , 2020
  * [107] Philip M Long and Hanie Sedghi “Generalization bounds for deep convolutional neural networks” In _ICLR_ , 2020
  * [108] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang and Yu Philip “A comprehensive survey on graph neural networks” In _IEEE transactions on neural networks and learning systems_ , 2020
  * [109] Keyulu Xu, Jingling Li, Mozhi Zhang, Simon S Du, Ken-ichi Kawarabayashi and Stefanie Jegelka “What can neural networks reason about?” In _ICLR_ , 2020
  * [110] Sanjeev Arora “Technical perspective: Why don’t today’s deep nets overfit to their training data?” In _Communications of the ACM_ , 2021
  * [111] Waiss Azizian and Marc Lelarge “Expressive power of invariant and equivariant graph neural networks” In _ICLR_ , 2021
  * [112] Constantinos Daskalakis, Stratis Skoulakis and Manolis Zampetakis “The complexity of constrained min-max optimization” In _STOC_ , 2021
  * [113] Pascal Esser, Leena Chennuru Vankadara and Debarghya Ghoshdastidar “Learning theory can (sometimes) explain generalisation in graph neural networks” In _NeurIPS_ , 2021
  * [114] Pierre Foret, Ariel Kleiner, Hossein Mobahi and Behnam Neyshabur “Sharpness-aware minimization for efficiently improving generalization” In _ICLR_ , 2021
  * [115] Henry Gouk, Timothy M Hospedales and Massimiliano Pontil “Distance-Based Regularisation of Deep Networks for Fine-Tuning” In _ICLR_ , 2021
  * [116] Moritz Hardt and Benjamin Recht “Patterns, predictions, and actions: A story about machine learning” In _arXiv preprint arXiv:2102.05242_ , 2021
  * [117] Hang Hua, Xingjian Li, Dejing Dou, Cheng-Zhong Xu and Jiebo Luo “Noise stability regularization for improving BERT fine-tuning” In _ACL_ , 2021
  * [118] Dongyue Li and Hongyang R Zhang “Improved regularization and robustness for fine-tuning in neural networks” In _NeurIPS_ , 2021
  * [119] Renjie Liao, Raquel Urtasun and Richard Zemel “A PAC-Bayesian Approach to Generalization Bounds for Graph Neural Networks” In _ICLR_ , 2021
  * [120] Keyulu Xu, Mozhi Zhang, Jingling Li, Simon S Du, Ken-ichi Kawarabayashi and Stefanie Jegelka “How neural networks extrapolate: From feedforward to graph neural networks” In _ICLR_ , 2021
  * [121] Gilad Yehudai, Ethan Fetaya, Eli Meirom, Gal Chechik and Haggai Maron “From local structures to size generalization in graph neural networks” In _ICML_ , 2021
  * [122] Stefanie Jegelka “Theory of Graph Neural Networks: Representation and Learning” In _arXiv preprint arXiv:2204.07697_ , 2022
  * [123] Haotian Ju, Dongyue Li and Hongyang R Zhang “Robust Fine-Tuning of Deep Neural Networks with Hessian-based Generalization Guarantees” In _ICML_ , 2022
  * [124] Haotian Ju, Dongyue Li and Hongyang R Zhang “Noise Stability Optimization for Flat Minima with Tight Rates” In _arXiv preprint arXiv:2306.08553_ , 2023

Organization. In Appendix A, we state the complete proofs for our results. In
Appendix B, we describe extra experiment details to complement our results.

##  Appendix A Proofs

This section provides the complete proofs for our results in Section 3. First,
we state several notations and facts needed in the proofs. Then we provide the
proof of the Hessian-based generalization bound for MPNN, as stated in Lemma
4.1. After that, in Appendix A.2, we provide the proof of Theorem 3.1, a key
step of which is the proof of Lemma 4.3. Next, in Appendix A.3, we state the
proof of the lower bound. Lastly, in Appendix A.4, we will provide proof for
the case of graph isomorphism networks.

First, we state several facts about graphs and provide a short proof of them.

######  Fact A.1.

Let $G=(V,E)$ be an undirected graph. Let $d_{G}$ be the maximum degree of
$G$.

  1. a)

Let $A$ be the adjacency matrix of $G$. Then, the adjacency matrix satisfies:
$\sqrt{d_{G}}\leq\left\|A\right\|\leq d_{G}$.

  2. b)

The symmetric and degree-normalized adjacency matrix satisfies
$\left\|D^{-1/2}AD^{-1/2}\right\|\leq 1$.

###### Proof.

Based on the definition of the spectral norm, we get

$\displaystyle\left\|A\right\|=\max_{\left\|x\right\|=1}x^{\top}Ax=\max_{\left\|x\right\|=1}\sum_{(i,j)\in
E}x_{i}x_{j}\leq\max_{\left\|x\right\|=1}\sum_{(i,j)\in
E}\frac{1}{2}(x_{i}^{2}+x_{j}^{2})\leq d_{G}\sum_{i\in V}x_{i}^{2}=d_{G}.$

Assume that node $i$ has the maximum degree $d_{G}$. Denote edges set
$E_{i}=\\{(i,i_{k})\\}_{k=1}^{d_{G}}\subseteq E$. Let
$x_{i}=\frac{1}{\sqrt{2}}$, $x_{i_{k}}=\frac{1}{\sqrt{2d_{G}}}$ for all
$k=1,\dots,d_{G}$. The rest entries of $x$ are equal to zero. Thus, $x$ is a
normalized vector. Next, we have

$\displaystyle\left\|A\right\|=\max_{\left\|x\right\|=1}\sum_{(i,j)\in
E}x_{i}x_{j}\geq\max_{\left\|x\right\|=1}2\sum_{(i,j)\in
E_{i}}x_{i}x_{j}=2d_{G}\cdot\frac{1}{\sqrt{2}}\frac{1}{\sqrt{2d_{G}}}=\sqrt{d_{G}}.$

An example in which $\left\|P_{G}\right\|$ gets close to $\sqrt{d_{G}}$ is the
star graph. An example in which $\left\|P_{G}\right\|$ gets close to $d_{G}$
is the complete graph.

Next, we focus on case b). From the definition of the spectral norm, we know

$\displaystyle\left\|D^{-1/2}AD^{-1/2}\right\|$
$\displaystyle=\max_{\left\|x\right\|=1}x^{\top}\big{(}{D^{-1/2}AD^{-1/2}}\big{)}x=\max_{\left\|x\right\|=1}\sum_{(i,j)\in
E}\frac{x_{i}x_{j}}{\sqrt{d_{i}d_{j}}}$
$\displaystyle\leq\max_{\left\|x\right\|=1}\sum_{(i,j)\in
E}\frac{x_{i}^{2}}{2d_{i}}+\frac{x_{j}^{2}}{2d_{j}}=\sum_{i\in V}x_{i}^{2}=1.$

During the middle of the above step, we used the Cauchy-Schwartz inequality.
The proof of this result is now completed. ∎

Notations: For two matrices $X$ and $Y$ that are both of dimension $d_{1}$ by
$d_{2}$, the Hadamard product of $X$ and $Y$, denoted as $X\odot Y$, is equal
to the entrywise product of $X$ and $Y$.

###  A.1 Proof of our PAC-Bayesian bound (Lemma 4.1)

To be precise, we will restate the conditions required in Theorem 3.1
separately below. The conditions are exactly the same as stated in Section 3.

######  Assumption A.2.

Assume that all the activation functions
$\phi_{i}(\cdot),\rho_{i}(\cdot),\psi_{i}(\cdot)$ for any $1\leq i\leq l-1$
and the loss function $\ell(x,y)$ over $x$ are twice-differentiable and
$\kappa_{0}$-Lipschitz. Their first-order derivatives are
$\kappa_{1}$-Lipschitz and their second-order derivatives are
$\kappa_{2}$-Lipschitz.

Based on the above assumption, we provide the precise statement for Taylor’s
expansion, used in equation (8).

######  Proposition A.3.

In the setting of Theorem 3.1, suppose each parameter in layer $i$ is
perturbed by an independent noise drawn from $\mathcal{N}(0,\sigma_{i}^{2})$.
Let $\tilde{\ell}(f(X,G),y)$ be the perturbed loss function with noise
perturbation injection vector $\mathcal{E}$ on all parameters $\bm{W}$ and
$\bm{U}$. There exist some fixed value $C_{1}$ that do not grow with $N$ and
$1/\delta$ such that

$\displaystyle{\left|{\tilde{\ell}(f(X,G),y)-\ell(f(X,G),y)-\frac{1}{2}\sum_{i=1}^{l}\sigma_{i}^{2}\operatorname{Tr}\left[\mathbf{H}^{{(i)}}[\ell(f(X,G),y)]\right]}\right|}\leq
C_{1}\sum_{i=1}^{l}\sigma_{i}^{3}.$

###### Proof.

By Taylor’s expansion, the following identity holds

$\displaystyle\tilde{\ell}(f(X,G),y)-\ell(f(X,G),y)=\mathop{\mathbb{E}}_{\mathcal{E}}\left[\mathcal{E}^{\top}\nabla\ell(f)+\frac{1}{2}{\mathcal{E}}^{\top}\mathbf{H}[\ell(f)]{\mathcal{E}}+R(\ell(f),\mathcal{E})\right].$

where $R(\ell(f),\mathcal{E})$ is the rest of the first-order and the second-
order terms. Since each entry in $\mathcal{E}$ follows the normal
distribution, we have
$\mathop{\mathbb{E}}_{\mathcal{E}}\left[\mathcal{E}^{\top}\nabla\ell(f)\right]=0$.
The Hessian term turns to

$\displaystyle{\mathcal{E}}^{\top}\mathbf{H}[\ell(f)]{\mathcal{E}}=\sum_{i=1}^{l}\sigma_{i}^{2}\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X,G),y)]\right].$

Since the readout layer is linear, by Proposition A.3, there exists a fixed
constant $\bar{C}$ that does not grow with $N$ and $\delta^{-1}$ such that
$|R(\ell(f),\mathcal{E})|\leq\bar{C}\left\|\mathcal{E}\right\|^{3}$. Based on
[93, Lemma 2 ], for any $x$ drawn from a normal distribution
$\mathcal{N}(0,\sigma^{2})$, we have
${\mathop{\mathbb{E}}\left[x^{3}\right]}\leq 6\sigma^{3}$. Hence, we get
$\mathop{\mathbb{E}}\left[R(\ell(f),\mathcal{E})\right]\leq
C_{1}\sum_{i=1}^{l}{\sigma_{i}^{3}}$, where
$C_{1}=\mathcal{\textup{O}}(h^{2}\bar{C})$ is a fixed constant. Thus, we have
finished the proof. ∎

Next, we state a Lipschitz-continuity upper bound of the network output at
each layer. This will be needed in the $\epsilon$-covering argument later in
the proof of Theorem 3.1. To simplify the notation, we will abbreviate
explicit constants that do not grow with $N$ and $1/\delta$ in the notation
$\lesssim$; more specifically, we use $A(n)\lesssim B(n)$ to indicate that
there exists a function $c$ that does not depend on $N$ and $1/\delta$ such
that $A(n)\leq c\cdot B(n)$ for large enough values of $n$.

######  Proposition A.4.

In the setting of Theorem 3.1, for any $j=1,\dots,l-1$, the change in the
Hessian of output of the $j$ layer network $H^{(j)}$ with respect to $W_{i}$
and $U_{i}$ under perturbation on $W$ and $U$ can be bounded as follows:

$\displaystyle\left\|\mathbf{H}_{\bm{W}}^{(i)}[\tilde{H}^{(j)}]-\mathbf{H}_{\bm{W}}^{(i)}[H^{(j)}]\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (12)
$\displaystyle\left\|\mathbf{H}_{\bm{U}}^{(i)}[\tilde{H}^{(j)}]-\mathbf{H}_{\bm{U}}^{(i)}[H^{(j)}]\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (13)

Above, the notation $\mathbf{H}_{\bm{W}}^{(i)}[\tilde{H}^{(j)}]$ is the
perturbation of the Hessian matrix of $H^{(j)}$ by $\Delta\bm{W}$ and
$\Delta\bm{U}$, specific to the variables of $\bm{W}$; likewise,
$\mathbf{H}_{\bm{U}}^{(i)}[\tilde{H}^{(j)}]$ is the perturbation of the
Hessian matrix specific to the variables of $\bm{U}$.

The proof of Proposition A.4 will be deferred until Appendix A.1.1. Based on
Propositions A.3 and A.4, now we are ready to present the proof of Lemma 4.1.

###### Proof of Lemma 4.1.

First, we separate the gap of ${L}(f)$ and $\frac{1}{\beta}\hat{{L}}(f)$ into
three parts:

$\displaystyle{L}(f)-\frac{1}{\beta}\hat{{L}}(f)=\underbrace{\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\ell(f(X,G),y)\right]-\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\tilde{\ell}(f(X,G),y)\right]}_{E_{1}}+\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\tilde{\ell}(f(X,G),y)\right]$
$\displaystyle-\frac{1}{\beta}\Big{(}{\frac{1}{N}\sum_{i=1}^{N}\tilde{\ell}(f(X_{i},G_{i}),y_{i})}\Big{)}+\underbrace{\frac{1}{\beta}\Big{(}{\frac{1}{N}\sum_{i=1}^{N}\tilde{\ell}(f(X_{i},G_{i}),y_{i})}\Big{)}-\frac{1}{\beta}\Big{(}{\frac{1}{N}\sum_{i=1}^{N}\ell(f(X_{i},G_{i}),y_{i})}\Big{)}}_{E_{2}}.$

for any $\beta\in(0,1)$. Above, $\tilde{\ell}(f(X,G),y)$ is the perturbed loss
from $\ell(f(X,G),y)$ with noise injections $\mathcal{E}$ added to all the
parameters in $\bm{W}$ and $\bm{U}$. By Taylor’s expansion from Proposition
A.3, we can bound the difference between $\tilde{\ell}(f(X,G),y)$ and
$\ell(f(X,G)$ with the trace of the Hessian. Therefore

$\displaystyle{L}(f)-\frac{1}{\beta}\hat{{L}}(f)\leq\,$
$\displaystyle-\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\frac{1}{2}\sum_{i=1}^{l}\sigma_{i}^{2}\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X,G),y)]\right]\right]+\sum_{i=1}^{l}C_{1}\sigma_{i}^{3}$
(by Prop. A.3 for $E_{1}$)
$\displaystyle+\Bigg{(}\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\tilde{\ell}(f(X,G),y)\right]-\frac{1}{\beta}\Big{(}{\frac{1}{N}\sum_{i=1}^{N}\tilde{\ell}(f(X_{i},G_{i}),y_{i})}\Big{)}\Bigg{)}$
$\displaystyle+\frac{1}{2\beta}\sum_{i=1}^{l}\sigma_{i}^{2}\Big{(}{\frac{1}{N}\sum_{j=1}^{N}\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X_{j},G_{j}),y_{j})]\right]}\Big{)}+\frac{1}{\beta}\sum_{i=1}^{l}C_{1}\sigma_{i}^{3}.$
(by Prop. A.3 for $E_{2}$)

By rearranging the above equation, we get the following:

$\displaystyle{L}(f)-\frac{1}{\beta}\hat{{L}}(f)$
$\displaystyle\leq\frac{1}{2}\sum_{i=1}^{l}\sigma_{i}^{2}\Bigg{(}\underbrace{\frac{1}{N}\sum_{j=1}^{N}\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X_{j},G_{j}),y_{j})]\right]-\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X,G),y)]\right]\right]}_{E_{3}}\Bigg{)}$
$\displaystyle+\frac{1}{2}\Big{(}{\frac{1}{\beta}-1}\Big{)}\underbrace{\sum_{i=1}^{l}\frac{\sigma_{i}^{2}}{N}{\sum_{j=1}^{N}\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X_{j},G_{j}),y_{j})]\right]}}_{E_{4}}$
$\displaystyle+\Big{(}{1+\frac{1}{\beta}}\Big{)}C_{1}\sum_{i=1}^{l}\sigma_{i}^{3}+\underbrace{\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\tilde{\ell}(f(X,G),y)\right]-{\frac{1}{\beta
N}\sum_{i=1}^{N}\tilde{\ell}(f(X_{i},G_{i}),y_{i})}}_{E_{5}}.$

Based on Proposition A.4, the Hessian operator $\mathbf{H}^{(i)}$ is
Lipschitz-continuous for some parameter that does not depend on $N$ and
$1/\delta$, for any $i=1,2\dots,l$. Therefore, from [123, Lemma 2.4 ], there
exist some fixed values $C_{2}$, $C_{3}$ that do not grow with $N$ and
$1/\delta$, such that with probability at least $1-\delta$ over the randomness
of the training set. Therefore, the matrix inside the trace of $E_{3}$
satisfies

$\displaystyle\left\|\frac{1}{N}\sum_{j=1}^{N}\mathbf{H}^{(i)}[\ell(f(X_{j},G_{j}),y_{j})]-\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\mathbf{H}^{(i)}[\ell(f(X,G),y)]\right]\right\|_{{F}}\leq\frac{C_{2}\sqrt{\log(C_{3}N/\delta)}}{\sqrt{N}},$
(14)

for any $i=1,\dots,l.$ Thus, by the Cauchy-Schwartz inequality, $E_{3}$ is
less than $\sqrt{2h^{2}}$ times the RHS of equation (14). Suppose the loss
function $\ell(f(X,G),y)$ lies in a bounded range $[0,B]$ given any
$(X,G,y)\sim\mathcal{D}$. By the PAC-Bayes bound of [65, Theorem 2 ] (see also
[92]), we choose $\bm{U}$ as a prior distribution and $\bm{W}+\bm{U}$ as a
posterior distribution. For any $\beta\in(0,1)$ and $\delta\in[0,1)$, with
probability at least $1-\delta$, $E_{5}$ satisfies:

$\displaystyle\mathop{\mathbb{E}}_{(X,G,y)\sim\mathcal{D}}\left[\tilde{\ell}(f(X,G),y)\right]-{\frac{1}{\beta
N}\sum_{i=1}^{N}\tilde{\ell}(f(X_{i},G_{i}),y_{i})}\leq$
$\displaystyle\frac{B}{2\beta(1-\beta)N}\Big{(}{\sum_{i=1}^{l}{\frac{\left\|W^{(i)}\right\|_{{F}}^{2}+\left\|U^{(i)}\right\|_{{F}}^{2}}{2\sigma_{i}^{2}}}+\log\frac{1}{\delta}}\Big{)}$
$\displaystyle\leq$
$\displaystyle\frac{B}{2\beta(1-\beta)N}\Big{(}{\sum_{i=1}^{l}\frac{s_{i}^{2}r_{i}^{2}}{\sigma_{i}^{2}}+\log\frac{1}{\delta}}\Big{)}.$
(15)

The above is because ${\bm{W}}$ and ${\bm{U}}$ are inside the hypothesis set
$\mathcal{H}$. For any $i=1,\dots,l$, let

$\alpha_{i}=\max\limits_{(X,G,y)\sim\mathcal{D}}\operatorname{Tr}\left[\mathbf{H}^{(i)}[\ell(f(X,G),y)]\right].$

Lastly, we use $\sigma_{i}^{2}\alpha_{i}$ above to upper bound $E_{4}$.
Combined with equations (14) and (15), with probability at least $1-2\delta$,
we get

$\displaystyle{L}(f)-\frac{1}{\beta}\hat{{L}}(f)\leq\,$
$\displaystyle\frac{C_{2}\sqrt{2h^{2}\log(C_{3}N/\delta)}}{\sqrt{N}}\sum_{i=1}^{l}\sigma_{i}^{2}+\Big{(}{1+\frac{1}{\beta}}\Big{)}C_{1}\sum_{i=1}^{l}\sigma_{i}^{3}$
$\displaystyle+{\frac{1}{2}\Big{(}{\frac{1}{\beta}-1}\Big{)}\sum_{i=1}^{l}\alpha_{i}\sigma_{i}^{2}+\frac{B}{2\beta(1-\beta)N}\Bigg{(}\sum_{i=1}^{l}{\frac{s_{i}^{2}r_{i}^{2}}{\sigma_{i}^{2}}}+\log\frac{1}{\delta}\Bigg{)}}.$

Next, we will select $\sigma_{i}$ to minimize the last line above. One can
verify that this is achieved when

$\displaystyle\sigma_{i}^{2}=\frac{s_{i}r_{i}}{1-\beta}\sqrt{\frac{B}{\alpha_{i}N}},\text{
for every }i=1,2,\dots,l.$

With this setting of the noise variance, the gap between ${L}(f)$ and
$\hat{{L}}(f)/{\beta}$ becomes:

$\displaystyle{L}(f)-\frac{1}{\beta}\hat{{L}}(f)$ $\displaystyle\leq$
$\displaystyle\frac{1}{\beta}\sum_{i=1}^{l}\sqrt{\frac{B\alpha_{i}s_{i}^{2}r_{i}^{2}}{N}}+\frac{C_{2}\sqrt{2h^{2}\log(C_{3}N/\delta)}}{\sqrt{N}}\sum_{i=1}^{L}\sigma_{i}^{2}+\Big{(}{1+\frac{1}{\beta}}\Big{)}C_{1}\sum_{i=1}^{l}\sigma_{i}^{3}+\frac{C}{2\beta(1-\beta)N}\log\frac{1}{\delta}.$

Let $\beta$ be a fixed value close to $1$ and independent of $N$ and
$\delta^{-1}$; let $\epsilon=(1-\beta)/\beta$. We get

$\displaystyle{L}(f)\leq\,(1+\epsilon)\hat{{L}}(f)+(1+\epsilon)\sum_{i=1}^{l}\sqrt{\frac{B\alpha_{i}r_{i}^{2}s_{i}^{2}}{N}}+\xi,\text{
where }$
$\displaystyle\xi=\frac{C_{2}\sqrt{2h^{2}\log(C_{3}N/\delta)}}{\sqrt{N}}\sum_{i=1}^{L}\sigma_{i}^{2}+\Big{(}{1+\frac{1}{\beta}}\Big{)}C_{1}\sum_{i=1}^{l}{\sigma_{i}^{3}}+\frac{C}{2\beta(1-\beta)N}\log\frac{1}{\delta}.$

Notice that $\xi$ is of order $\mathcal{\textup{O}}(N^{-3/4}\
+\log(\delta^{-1})N^{-1})\leq\mathcal{\textup{O}}(\log(\delta^{-1})/N^{3/4})$.
Therefore, we have finished the proof of equation (7). ∎

####  A.1.1 Proof of Proposition A.4

For any $j=1,2,\dots,l$, let $\tilde{H}^{(j)}$ be the perturbed network output
after layer $j$, with perturbations given by $\Delta\bm{W}$ and
$\Delta\bm{U}$. We show the following Lipschitz-continuity property for
$H^{(j)}$.

######  Claim A.5.

Suppose that Assumption A.2 holds. For any $j=1,\dots,l-1$, the change in the
output of the $j$ layer network $H^{(j)}$ with perturbation added to $\bm{W}$
and $\bm{U}$ can be bounded as follows:

$\displaystyle\left\|\tilde{H}^{(j)}-H^{(j)}\right\|_{{F}}\lesssim\sum_{t=1}^{j}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (16)

###### Proof.

We will prove using induction with respect to $j$. If $j=1$, we have

$\displaystyle\left\|\phi_{1}\Big{(}{X\big{(}{U^{(1)}+\Delta
U^{(1)}}\big{)}+\rho_{1}\big{(}{P_{{}_{G}}\psi_{1}\big{(}{X}\big{)}}\big{)}\big{(}{W^{(1)}+\Delta
W^{(1)}}\big{)}}\Big{)}-\phi_{1}\Big{(}{XU^{(1)}+\rho_{1}\big{(}{P_{{}_{G}}\psi_{1}\big{(}{X}\big{)}}\big{)}W^{(1)}}\Big{)}\right\|_{{F}}$
$\displaystyle\leq$ $\displaystyle\kappa_{0}\left\|X\Delta
U^{(1)}+\rho_{1}\big{(}{P_{{}_{G}}\psi_{1}\big{(}{X}\big{)}}\big{)}\Delta
W^{(1)}\right\|_{{F}}\lesssim\left\|\Delta U^{(1)}\right\|+\left\|\Delta
W^{(1)}\right\|.$

Hence, we know that equation (16) will be correct when $j=1$. Assuming that
equation (16) is correct for any $j\geq 1$, the perturbation of layer $j+1$’s
network output $H^{(j+1)}$ is less than

$\displaystyle\left\|\tilde{H}^{(j+1)}-H^{(j+1)}\right\|_{{F}}$
$\displaystyle\leq$ $\displaystyle\,\kappa_{0}\left\|X\Delta
U^{(j+1)}+\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}\big{(}{\tilde{H}^{(j)}}\big{)}}\big{)}\big{(}{W^{(j+1)}+\Delta
W^{(j+1)}}\big{)}-\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}\big{(}{H^{(j)}}\big{)}}\big{)}W^{(j+1)}\right\|_{{F}}$
$\displaystyle\lesssim$ $\displaystyle\left\|\Delta
U^{(j+1)}\right\|+\left\|\Delta
W^{(j+1)}\right\|+\left\|\tilde{H}^{(j)}-H^{(j)}\right\|_{{F}}.$

Thus, we have finished the proof of the induction step. ∎

Next, for any $i$ and $j$, let $\frac{\partial\tilde{H}^{(j)}}{\partial
W^{(i)}}$ be the perturbation of the partial derivative of $H^{(j)}$ with
perturbations given by $\Delta\bm{W}$ and $\Delta\bm{U}$.

######  Claim A.6.

Suppose that Assumption A.2 holds. For any $j=1,\dots,l-1$, the change in the
Jacobian of the $j$-th layer’s output $H^{(j)}$ with respect to $W^{(i)}$ and
$U^{(i)}$ satisfies:

$\displaystyle\left\|\frac{{{\partial\tilde{H}}}^{(j)}}{\partial
W^{(i)}}-\frac{\partial H^{(j)}}{\partial W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (17)
$\displaystyle\left\|\frac{{\partial\tilde{H}}^{(j)}}{\partial
U^{(i)}}-\frac{\partial H^{(j)}}{\partial U^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (18)

###### Proof.

We will consider a fixed $i=1,\dots,l-1$ and take induction over
$j=i,\dots,l-1$. We focus on the proof of equation (17), while the proof of
equation (18) will be similar. To simplify the derivation, we use two
notations for brevity. Let

$\displaystyle F_{j}=P_{{}_{G}}\psi_{j}\big{(}{H^{(j-1)}}\big{)}W^{(j)}\text{
and }E_{j}=XU^{(j)}+\rho_{j}\big{(}{F_{j}}\big{)}.$

First, we consider the base case when $j=i$. By the chain rule, we have:

$\displaystyle\left\|\frac{\partial\tilde{H}^{(i)}}{\partial
W^{(i)}}-\frac{\partial H^{(i)}}{\partial W^{(i)}}\right\|_{{F}}=$
$\displaystyle\left\|\phi^{\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}\odot\frac{\partial\tilde{E}_{i}}{\partial
W^{(i)}}-\phi^{\prime}_{i}\big{(}{E_{i}}\big{)}\odot\frac{\partial
E_{i}}{\partial W^{(i)}}\right\|_{{F}}$ $\displaystyle\lesssim$
$\displaystyle\left\|\phi^{\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{E}_{i}}{\partial
W^{(i)}}-\frac{\partial E_{i}}{\partial W^{(i)}}\right\|_{{F}}.$

From Claim A.5, we know

$\displaystyle\left\|\phi^{\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}\leq\kappa_{1}\left\|\tilde{E}_{i}-E_{i}\right\|_{{F}}\lesssim\left\|\Delta
W^{(i)}\right\|+\left\|\Delta U^{(i)}\right\|.$

By the chain rule again, we get:

$\displaystyle\left\|\frac{\partial\tilde{E}_{i}}{\partial
W^{(i)}}-\frac{\partial E_{i}}{\partial
W^{(i)}}\right\|_{{F}}\lesssim\left\|\rho^{\prime}_{i}\big{(}{\tilde{F}_{i}}\big{)}-\rho^{\prime}_{i}\big{(}{F_{i}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{F}_{i}}{\partial
W^{(i)}}-\frac{\partial F_{i}}{\partial
W^{(i)}}\right\|_{{F}}\lesssim{\left\|\Delta W^{(i)}\right\|+\left\|\Delta
U^{(i)}\right\|}.$ (by Claim A.5 again)

Hence, we know that equation (17) will be correct when $j=i$. Assuming that
equation (17) will be correct for any $j$ up to $j\geq i$, we have

$\displaystyle\left\|\frac{\partial\tilde{H}^{(j+1)}}{\partial
W^{(i)}}-\frac{\partial H^{(j+1)}}{\partial W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim$
$\displaystyle\left\|\phi^{\prime}_{j+1}\big{(}{\tilde{E}_{j+1}}\big{)}-\phi^{\prime}_{j+1}\big{(}{E_{j+1}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{E}_{j+1}}{\partial
W^{(i)}}-\frac{\partial E_{j+1}}{\partial W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim$ $\displaystyle\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta
W^{(t)}\right\|}\Big{)}+\left\|\rho^{\prime}_{j+1}\big{(}{\tilde{F}_{j+1}}\big{)}-\rho^{\prime}_{j+1}\big{(}{F_{j+1}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{F}_{j+1}}{\partial
W^{(i)}}-\frac{\partial F_{j+1}}{W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim$ $\displaystyle\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta
W^{(t)}\right\|}\Big{)}+\left\|\psi^{\prime}_{j+1}\big{(}{\tilde{H}^{(j)}}\big{)}-\psi^{\prime}_{j+1}\big{(}{H^{(j)}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{H}^{(j)}}{\partial
W^{(i)}}-\frac{\partial H^{(j)}}{\partial W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim$ $\displaystyle\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (by Claim A.5 and the
induction step)

The above steps all use Claim A.5. The last step additionally uses the
induction hypothesis. From repeatedly applying the above beginning with $j=i$
along with the base case of equation (17), we conclude that equation (17)
holds.

Next, we consider the base case for equation (18). For the base case $j=i$,
from the chain rule, by Claim A.5, we get:

$\displaystyle\left\|\frac{\partial\tilde{H}^{(i)}}{\partial
U^{(i)}}-\frac{\partial H^{(i)}}{\partial
U^{(i)}}\right\|_{{F}}\lesssim\left\|\phi^{\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{E}_{i}}{\partial
U^{(i)}}-\frac{\partial E_{i}}{\partial
U^{(i)}}\right\|_{{F}}\lesssim{\left\|\Delta W^{(i)}\right\|+\left\|\Delta
U^{(i)}\right\|}.$

Hence, we know that equation (18) will be correct when $j=i$. Assuming that
equation (18) will be correct for any $j$ up to $j\geq i$, we have

$\displaystyle\left\|\frac{\partial\tilde{H}^{(j+1)}}{\partial
U^{(i)}}-\frac{\partial H^{(j+1)}}{\partial
U^{(i)}}\right\|_{{F}}\lesssim\left\|\phi^{\prime}_{j+1}\big{(}{\tilde{E}_{j+1}}\big{)}-\phi^{\prime}_{j+1}\big{(}{E_{j+1}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{E}_{j+1}}{\partial
U^{(i)}}-\frac{\partial E_{j+1}}{\partial U^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta
W^{(t)}\right\|}\Big{)}+\left\|\rho^{\prime}_{j+1}\big{(}{\tilde{F}_{j+1}}\big{)}-\rho^{\prime}_{j+1}\big{(}{F_{j+1}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{F}_{j+1}}{\partial
U^{(i)}}-\frac{\partial F_{j+1}}{\partial U^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta
W^{(t)}\right\|}\Big{)}+\left\|\psi^{\prime}_{j+1}\big{(}{\tilde{H}^{(j)}}\big{)}-\psi^{\prime}_{j+1}\big{(}{H^{(j)}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{H}^{(j)}}{\partial
U^{(i)}}-\frac{\partial H^{(j)}}{\partial U^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
U^{(t)}\right\|+\left\|\Delta W^{(t)}\right\|}\Big{)}.$ (by Claim A.5 and the
induction step)

The second and third steps are based on Claim A.5. From repeatedly applying
the above beginning with $j=i$ along with the base case of equation (18), we
conclude that equation (18) holds. The proof of claim A.6 is complete. ∎

###### Proof of Proposition A.4.

We will consider a fixed $i=1,\dots,l-1$ and take induction over
$j=i,\dots,l-1$. We focus on the proof of equation (12), while the proof of
equation (13) will be similar. To simplify the derivation, we use two
notations for brevity. Let

$\displaystyle F_{j}=P_{{}_{G}}\psi_{j}\big{(}{H^{(j-1)}}\big{)}W^{(j)}\text{
and }E_{j}=XU^{(j)}+\rho_{j}\big{(}{F_{j}}\big{)}.$

First, we consider the base case when $j=i$. By the chain rule, we have: We
use the chain rule to get:

$\displaystyle\frac{\partial^{2}H^{(i)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}=\phi^{\prime\prime}_{i}(E_{i})\odot\frac{\partial
E_{i}}{\partial W_{p,q}^{(i)}}\odot\frac{\partial E_{i}}{\partial
W_{p,q}^{(i)}}+\phi^{\prime}_{i}(E_{i})\odot\rho^{\prime\prime}_{i}(F_{i})\odot{\frac{\partial
F_{i}}{\partial W_{p,q}^{(i)}}}\odot{\frac{\partial F_{i}}{\partial
W^{(i)}_{p,q}}}.$

Hence, the Frobenius norm of the Hessian of $H^{(i)}$ with respect to $W_{i}$
under perturbation on $W$ and $U$ turns to

$\displaystyle\left\|\mathbf{H}_{\bm{W}}^{(i)}[\tilde{H}^{(i)}]-\mathbf{H}_{\bm{W}}^{(i)}[H^{(i)}]\right\|_{{F}}$
$\displaystyle\lesssim\left\|\phi^{\prime\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{E}_{i}}{\partial
W^{(i)}}-\frac{\partial E_{i}}{\partial
W^{(i)}}\right\|_{{F}}+\left\|\phi^{\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}$
$\displaystyle+\left\|\rho^{\prime\prime}_{i}\big{(}{\tilde{F}_{i}}\big{)}-\rho^{\prime\prime}_{i}\big{(}{F_{i}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{F}_{i}}{\partial
W^{(i)}}-\frac{\partial F_{i}}{\partial W^{(i)}}\right\|_{{F}}.$

From Claim A.5, we know

$\displaystyle\left\|\phi^{\prime\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}$
$\displaystyle\leq\kappa_{2}\left\|\tilde{E}_{i}-E_{i}\right\|_{{F}}\lesssim\left\|\Delta
W^{(i)}\right\|+\left\|\Delta U^{(i)}\right\|,$
$\displaystyle\left\|\phi^{\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}$
$\displaystyle\leq\kappa_{1}\left\|\tilde{E}_{i}-E_{i}\right\|_{{F}}\lesssim\left\|\Delta
W^{(i)}\right\|+\left\|\Delta U^{(i)}\right\|,$
$\displaystyle\left\|\rho^{\prime\prime}_{i}\big{(}{\tilde{F}_{i}}\big{)}-\rho^{\prime\prime}_{i}\big{(}{F_{i}}\big{)}\right\|_{{F}}$
$\displaystyle\leq\kappa_{2}\left\|\tilde{F}_{i}-F_{i}\right\|_{{F}}\lesssim\left\|\Delta
W^{(i)}\right\|+\left\|\Delta U^{(i)}\right\|.$

From Claim A.6, we have

$\displaystyle\left\|\frac{\partial\tilde{E}_{i}}{\partial
W^{(i)}}-\frac{\partial E_{i}}{\partial W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\left\|\Delta W^{(i)}\right\|+\left\|\Delta
U^{(i)}\right\|,$ $\displaystyle\left\|\frac{\partial\tilde{F}_{i}}{\partial
W^{(i)}}-\frac{\partial F_{i}}{\partial W^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim\left\|\Delta W^{(i)}\right\|+\left\|\Delta
U^{(i)}\right\|.$

Hence, we know that equation (12) will be correct when $j=i$. Assuming that
equation (12) will be correct for any $j$ up to $j\geq i$, we can get the
following steps, by taking another derivative of the first-order derivative,
we can get the following steps:

$\displaystyle\frac{\partial^{2}H^{(j+1)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}=$
$\displaystyle\,\phi^{\prime\prime}_{j+1}(E_{j+1})\odot\frac{\partial
E_{j+1}}{\partial W_{p,q}^{(i)}}\odot\frac{\partial E_{j+1}}{\partial
W_{p,q}^{(i)}}+\phi^{\prime}_{j+1}(E_{j+1})\odot\rho^{\prime\prime}_{j+1}(F_{j+1})\odot\frac{\partial
F_{j+1}}{\partial W_{p,q}^{(i)}}\odot\frac{\partial F_{j+1}}{\partial
W_{p,q}^{(i)}}$
$\displaystyle+\phi^{\prime}_{j+1}(E_{j+1})\odot\rho^{\prime}_{j+1}(F_{j+1})\odot{P_{{}_{G}}\Big{(}{\psi^{\prime\prime}_{j+1}(H^{(j)})\odot\frac{\partial
H^{(j)}}{\partial W_{p,q}^{(i)}}\odot\frac{\partial H^{(j)}}{\partial
W_{p,q}^{(i)}}+\psi^{\prime}_{j+1}(H^{(j)})\odot\frac{\partial^{2}H^{(j)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}}\Big{)}W^{(j+1)}}.$

Thus, the Frobenius norm of the Hessian of $H^{(j+1)}$ with respect to
$W^{(i)}$ satisfies:

$\displaystyle\left\|\mathbf{H}_{\bm{W}}^{(i)}[\tilde{H}^{(j+1)}]-\mathbf{H}_{\bm{W}}^{(i)}[H^{(j+1)}]\right\|_{{F}}\lesssim\underbrace{\left\|\phi^{\prime\prime}_{j+1}\big{(}{\tilde{E}_{j+1}}\big{)}-\phi^{\prime\prime}_{j+1}\big{(}{E_{j+1}}\big{)}\right\|_{{F}}}_{A_{1}}+\underbrace{\left\|\frac{\partial\tilde{E}_{j+1}}{\partial
W^{(i)}}-\frac{\partial E_{j+1}}{\partial W^{(i)}}\right\|_{{F}}}_{B_{1}}$
$\displaystyle+$
$\displaystyle\underbrace{\left\|\phi^{\prime}_{j+1}\big{(}{\tilde{E}_{j+1}}\big{)}-\phi^{\prime}_{j+1}\big{(}{E_{j+1}}\big{)}\right\|_{{F}}}_{A_{2}}+\underbrace{\left\|\rho^{\prime\prime}_{j+1}\big{(}{\tilde{F}_{j+1}}\big{)}-\rho^{\prime\prime}_{j+1}\big{(}{F_{j+1}}\big{)}\right\|_{{F}}}_{A_{3}}$
$\displaystyle+$
$\displaystyle\underbrace{\left\|\frac{\partial\tilde{F}_{j+1}}{\partial
W^{(i)}}-\frac{\partial F_{j+1}}{\partial
W^{(i)}}\right\|_{{F}}}_{B_{2}}+\underbrace{\left\|\rho^{\prime}_{j+1}\big{(}{\tilde{F}_{j+1}}\big{)}-\rho^{\prime}_{j+1}\big{(}{F_{j+1}}\big{)}\right\|_{{F}}}_{A_{4}}$
$\displaystyle+$
$\displaystyle\underbrace{\left\|\psi^{\prime\prime}_{j+1}\big{(}{\tilde{H}^{(j)}}\big{)}-\psi^{\prime\prime}_{j+1}\big{(}{H^{(j)}}\big{)}\right\|_{{F}}}_{A_{5}}+\underbrace{\left\|\frac{\partial\tilde{H}^{(j)}}{\partial
W^{(i)}}-\frac{\partial H^{(j)}}{\partial W^{(i)}}\right\|_{{F}}}_{B_{3}}$
$\displaystyle+$
$\displaystyle\underbrace{\left\|\psi^{\prime}_{j+1}\big{(}{\tilde{H}^{(j)}}\big{)}-\psi^{\prime}_{j+1}\big{(}{H^{(j)}}\big{)}\right\|_{{F}}}_{A_{6}}+\underbrace{\left\|\mathbf{H}_{\bm{W}}^{(i)}[\tilde{H}^{(j)}]-\mathbf{H}_{\bm{W}}^{(i)}[H^{(j)}]\right\|_{{F}}}_{C_{1}}.$

Similarly, by Claim A.5, we get

$\displaystyle A_{i}\lesssim\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
W^{(t)}\right\|+\left\|\Delta U^{(t)}\right\|}\Big{)},\text{ for }1\leq i\leq
6.$

By Claim A.6, we get

$\displaystyle B_{i}\lesssim\sum_{t=1}^{j+1}\Big{(}{\left\|\Delta
W^{(t)}\right\|+\left\|\Delta U^{(t)}\right\|}\Big{)},\text{ for }1\leq i\leq
3.$

By the induction hypothesis, $C_{1}$ is also less than the above quantity.
From repeatedly applying the above beginning with $j=i$ along with the base
case of equation (12), we conclude that equation (12) holds.

Next, we consider the base case for equation (13). For the base case $j=i$,
from the chain rule, we get:

$\displaystyle\left\|\mathbf{H}_{\bm{U}}^{(i)}[\tilde{H}^{(i)}]-\mathbf{H}_{\bm{U}}^{(i)}[H^{(i)}]\right\|_{{F}}\lesssim$
$\displaystyle\left\|\phi^{\prime\prime}_{i}\big{(}{\tilde{E}_{i}}\big{)}-\phi^{\prime\prime}_{i}\big{(}{E_{i}}\big{)}\right\|_{{F}}+\left\|\frac{\partial\tilde{E}_{i}}{\partial
U^{(i)}}-\frac{\partial E_{i}}{\partial U^{(i)}}\right\|_{{F}}$
$\displaystyle\lesssim$
$\displaystyle\kappa_{2}\left\|\tilde{E}_{i}-E_{i}\right\|_{{F}}+{\left\|\Delta
W^{(i)}\right\|+\left\|\Delta U^{(i)}\right\|}$ (by Claim A.6)
$\displaystyle\lesssim$ $\displaystyle{\left\|\Delta
W^{(i)}\right\|+\left\|\Delta U^{(i)}\right\|}.$ (by Claim A.5)

Hence, we know that equation (13) will be correct when $j=i$. Assuming that
equation (13) will be correct for any $j$ up to $j\geq i$, we obtain the
induction step similar to the proof of equation (12), by Claim A.5, Claim A.6,
and the induction hypothesis, we conclude that equation (13) holds. ∎

###  A.2 Proof for message passing graph neural networks

Next, we present proof for message-passing graph neural networks. First, in
Appendix A.2.1, we derive the trace bound, which separates the trace of the
Hessian matrix into each entry of the weight matrices. Then in Appendix A.2.2
and A.2.3, we provide bounds on the first-order and second-order derivatives
of the Hessian matrix. Last, in Appendix A.2.4, building on these results, we
finish the proof of Theorem 3.1.

####  A.2.1 Proof of Lemma 4.3

###### Proof of Lemma 4.3.

Notice that $f(X,G)=H^{(l)}$. Recall that in each layer for $1\leq i\leq l-1$,
there are two weight matrices, a $d_{i-1}$ by $d_{i}$ matrix denoted as
$W^{(i)}$, and a $d_{0}$ by $U^{(i)}$ matrix denoted as $U^{(i)}$. To deal
with the trace of the Hessian $\mathbf{H}^{(i)}$, we first notice that there
are two parts in the trace:

$\displaystyle{\left|{\operatorname{Tr}\big{[}\mathbf{H}^{(i)}[\ell(H^{(l)},y)]\big{]}}\right|}\leq\underbrace{{\left|{\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\frac{\partial^{2}\ell(H^{(l)},y)}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}}\right|}}_{T_{1}}+\underbrace{{\left|{\sum_{p=1}^{d_{0}}\sum_{q=1}^{d_{i}}\frac{\partial^{2}\ell(H^{(l)},y)}{\partial\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}}\right|}}_{T_{2}}.$

We can inspect $T_{1}$ and $T_{2}$ in the above step separately. First, we
expand out the second-order derivatives in $T_{1}$. This will involve two
terms by the chain rule.

$\displaystyle T_{1}=\,$
$\displaystyle{\left|{\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\Bigg{\langle}\frac{\partial\ell(H^{(l)},y)}{\partial
H^{(l)}},\frac{\partial^{2}H^{(l)}}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}\Bigg{\rangle}}\right|}+{\left|{\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\Bigg{\langle}\frac{\partial^{2}\ell(H^{(l)},y)}{\partial\big{(}H^{(l)}\big{)}^{2}}\frac{\partial
H^{(l)}}{\partial W^{(i)}_{p,q}},\frac{\partial H^{(l)}}{\partial
W^{(i)}_{p,q}}\Bigg{\rangle}}\right|}$ $\displaystyle\leq\,$
$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial\ell(H^{(l)},y)}{\partial
H^{(l)}}\right\|\left\|\frac{\partial^{2}H^{(l)}}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}\right\|+\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}\ell(H^{(l)},y)}{\partial\big{(}H^{(l)}\big{)}^{2}}\right\|\left\|\frac{\partial
H^{(l)}}{\partial W^{(i)}_{p,q}}\right\|^{2}$ $\displaystyle\leq\,$
$\displaystyle\kappa_{0}\sqrt{k}\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(l)}}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}\right\|+\kappa_{1}k\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(l)}}{\partial W^{(i)}_{p,q}}\right\|^{2}.$ (19)

The last step is because $\ell(\cdot)$ is $\kappa_{0}$-Lipschitz continuous
and $\ell^{\prime}(\cdot)$ is $\kappa_{1}$-Lipschitz continuous, under
Assumption A.2. Thus, the Euclidean norm of
${\frac{\partial\ell(H^{(l)},y)}{\partial H^{(l)}}}$ is at most
$\kappa_{0}\sqrt{k}$, since $H^{(l)}$ is a $k$-dimensional vector. Recall from
step (2) that $H^{(l)}=\frac{1}{n}\bm{1}_{n}^{\top}H^{(l-1)}W^{(l)}$. Hence,
we have

$\displaystyle\left\|\frac{\partial H^{(l)}}{\partial W^{(i)}_{p,q}}\right\|$
$\displaystyle=\left\|\frac{1}{n}\bm{1}_{n}^{\top}\frac{\partial
H^{(l-1)}}{\partial W^{(i)}_{p,q}}W^{(l)}\right\|$
$\displaystyle\leq\left\|\frac{1}{n}\bm{1}_{n}^{\top}\right\|\left\|\frac{\partial
H^{(l-1)}}{\partial
W^{(i)}_{p,q}}W^{(l)}\right\|\leq\frac{1}{\sqrt{n}}\left\|\frac{\partial
H^{(l-1)}}{\partial W^{(i)}_{p,q}}\right\|\left\|W^{(l)}\right\|.$ (20)

In a similar vein, the Euclidean norm of
${\frac{\partial^{2}\ell(H^{(l)},y)}{\partial(H^{(l)})^{2}}}$ is at most
$\kappa_{1}k$, since the second-order derivatives become a $k$ by $k$ matrix.
Then, we get

$\displaystyle\left\|\frac{\partial^{2}H^{(l)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}\right\|$
$\displaystyle=\left\|\frac{1}{n}\bm{1}_{n}^{\top}\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}W^{(l)}\right\|$
$\displaystyle\leq\left\|\frac{1}{n}\bm{1}_{n}^{\top}\right\|\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}W^{(l)}\right\|\leq\frac{1}{\sqrt{n}}\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}\right\|\left\|W^{(l)}\right\|.$
(21)

After substituting equations (20) and (21) into equation (19), we get:

$\displaystyle T_{1}$
$\displaystyle\leq\frac{\kappa_{0}\sqrt{k}}{\sqrt{n}}\left\|W^{(l)}\right\|\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}\right\|+\frac{\kappa_{1}k}{n}\left\|W^{(l)}\right\|^{2}\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(l-1)}}{\partial W^{(i)}_{p,q}}\right\|^{2}.$
$\displaystyle\leq\frac{\kappa_{0}\sqrt{k}}{\sqrt{n}}\left\|W^{(l)}\right\|\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}W^{(i)}_{p,q}\big{)}^{2}}\right\|_{{F}}+\frac{\kappa_{1}k}{n}\left\|W^{(l)}\right\|^{2}\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(l-1)}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}.$

The proof for the case of $T_{2}$ concerning $U^{(i)}$ follows the same steps
as above. Without belaboring all the details, one can get that

$\displaystyle
T_{2}\leq\frac{\kappa_{0}\sqrt{k}}{\sqrt{n}}\left\|W^{(l)}\right\|\sum_{p=1}^{d_{0}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(l-1)}}{\partial\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}+\frac{\kappa_{1}k}{n}\left\|W^{(l)}\right\|^{2}\sum_{p=1}^{d_{0}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(l-1)}}{\partial U_{p,q}^{(i)}}\right\|_{{F}}^{2}.$ (22)

This completes the proof of Lemma 4.3. ∎

####  A.2.2 Dealing with first-order derivatives

Based on Lemma 4.3, the analysis involves two parts, one on the first-order
derivatives of $H^{(j)}$ for all layers $j$, and the other on the second-order
derivatives of $H^{(j)}$ for all layers $j$.

######  Proposition A.7.

In the setting of Theorem 3.1, the first-order derivative of $H^{(j)}$ with
respect to $W^{(i)}$ and $U^{(i)}$ satisfies the following, for any
$i=1,\dots,l-1$ and $j\geq i$:

$\displaystyle\left\|\frac{\partial H^{(j)}}{\partial
W^{(i)}}\right\|_{{F}}\leq\,$
$\displaystyle\kappa_{0}^{3(j-i+1)}\sqrt{d_{i}}\left\|P_{{}_{G}}\right\|^{j-i+1}\left\|H^{(i-1)}\right\|_{{F}}\prod_{t=i+1}^{j}\left\|W^{(t)}\right\|,$
(23) $\displaystyle\left\|\frac{\partial H^{(j)}}{\partial
U^{(i)}}\right\|_{{F}}\leq\,$
$\displaystyle\kappa_{0}^{3(j-i)+1}\sqrt{d_{i}}\left\|P_{{}_{G}}\right\|^{j-i+1}\left\|X\right\|_{{F}}\prod_{t=i+1}^{j}\left\|W^{(t)}\right\|.$
(24)

###### Proof.

We will consider a fixed $i=1,\dots,l-1$ and take induction over
$j=i,\dots,l-1$. We focus on the proof of equation (23), while the proof of
equation (24) will be similar. First, we consider the base case when $j=i$.
Let $W^{(i)}_{p,q}$ be the $(p,q)$-th entry of $W^{(i)}$, for any valid
indices $p$ and $q$. Recall that $\phi_{i}(\cdot)$ is $\kappa_{0}$-Lipschitz
continuous from Assumption A.2, for any $i=1,\dots,l-1$. Therefore,

$\displaystyle\left\|\phi^{\prime}_{i}(x)\right\|_{\infty}\leq\kappa_{0},~{}~{}\left\|\psi^{\prime}_{i}(x)\right\|_{\infty}\leq\kappa_{0},~{}~{}\text{and}~{}~{}\left\|\rho^{\prime}_{i}(x)\right\|_{\infty}\leq\kappa_{0}.$
(25)

For each $(p,q)$-entry of $W^{(i)}$, by the chain rule, we have:

$\displaystyle\left\|\frac{\partial H^{(i)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$
$\displaystyle=\left\|\phi^{\prime}_{i}\big{(}{XU^{(i)}+\rho_{i}\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}\big{)}\odot\frac{\partial\big{(}{XU^{(i)}+\rho_{i}\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}\big{)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$ (26)
$\displaystyle\leq\kappa_{0}\left\|\frac{\partial\rho_{i}\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$ (by equation (25))
$\displaystyle=\kappa_{0}\left\|\rho^{\prime}_{i}\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}\odot\frac{\partial\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}{W^{(i)}_{p,q}}\right\|_{{F}}$
$\displaystyle\leq\kappa_{0}^{2}\left\|\frac{\partial\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}.$ (again by equation (25))

Notice that only the $q$-th column of the derivative
$P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}$ is nonzero, which is equal to the
$p$’th column of $P_{{}_{G}}\psi_{i}(H^{(i-1)})$. Thus, the Jacobian of
$H^{(i)}$ over $W^{(i)}$ satisfies:

$\displaystyle\left\|\frac{\partial H^{(i)}}{\partial W^{(i)}}\right\|_{{F}}=$
$\displaystyle\sqrt{\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(i)}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}}$ $\displaystyle\leq$
$\displaystyle\kappa_{0}^{2}\sqrt{\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}{W_{p,q}^{(i)}}\right\|_{{F}}^{2}}=\kappa_{0}^{2}\sqrt{d_{i}}\left\|P_{{}_{G}}\psi_{i}(H^{(i-1)})\right\|_{{F}}.$
(27)

Therefore, the above equation (27) implies that equation (23) holds in the
base case. Next, we consider the induction step from layer $j$ to layer $j+1$.
The derivative of $H^{(j+1)}$ with respect to $W^{(i)}_{p,q}$ satisfies:

$\displaystyle\left\|\frac{\partial H^{(j+1)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$ $\displaystyle=$
$\displaystyle\left\|\phi^{\prime}_{j+1}\big{(}{XU^{(j+1)}+\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}\big{)}\odot\frac{\partial\big{(}{XU^{(j+1)}+\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}\big{)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$ $\displaystyle\leq$
$\displaystyle\kappa_{0}\left\|\frac{\partial\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$ (by equation (25)) $\displaystyle\leq$
$\displaystyle\kappa_{0}\left\|\rho^{\prime}_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}\odot\frac{\partial\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}$ $\displaystyle\leq$
$\displaystyle\kappa_{0}^{2}\left\|P_{{}_{G}}\frac{\partial\psi_{j+1}(H^{(j)})}{\partial
W^{(i)}_{p,q}}W^{(j+1)}\right\|_{{F}}$ (again by equation (25))

By applying equation (25) w.r.t. $\psi^{\prime}_{j+1}$, The above is less
than:

$\displaystyle\kappa_{0}^{2}\left\|P_{{}_{G}}\right\|\left\|\psi_{j+1}^{\prime}(H^{(j)})\odot\frac{\partial
H^{(j)}}{\partial
W^{(i)}_{p,q}}\right\|_{{F}}\left\|W^{(j+1)}\right\|\leq\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\left\|\frac{\partial
H^{(j)}}{\partial W^{(i)}_{p,q}}\right\|_{{F}}.$

Hence, the Jacobian of $H^{(j+1)}$ with respect to $W^{(i)}$ satisfies:

$\displaystyle\left\|\frac{\partial H^{(j+1)}}{\partial
W^{(i)}}\right\|_{{F}}\leq\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\left\|\frac{\partial
H^{(j)}}{\partial W^{(i)}}\right\|_{{F}}.$

From repeatedly applying the above beginning with $j=i$ along with the base
case of equation (27), we conclude that equation (23) holds.

Next, we consider the base case for equation (24). For each $(p,q)$-th entry
of $U^{(i)}$, from the chain rule we get:

$\displaystyle\left\|\frac{\partial H^{(i)}}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}$
$\displaystyle=\left\|\phi^{\prime}_{i}\Big{(}{XU^{(i)}+\rho_{i}\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}\Big{)}\odot\frac{\partial\Big{(}{XU^{(i)}+\rho_{i}\big{(}{P_{{}_{G}}\psi_{i}(H^{(i-1)})W^{(i)}}\big{)}}\Big{)}}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}$
$\displaystyle\leq\kappa_{0}\left\|\frac{\partial(XU^{(i)})}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}.$ (by equation (25))

Therefore, by summing over $p=1,\dots,d_{0}$ and $q=1,\dots,d_{i}$, we get:

$\displaystyle\left\|\frac{\partial H^{(i)}}{\partial U^{(i)}}\right\|_{{F}}$
$\displaystyle=\sqrt{\sum_{p=1}^{d_{0}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(i)}}{\partial U^{(i)}_{p,q}}\right\|_{{F}}^{2}}$
$\displaystyle\leq\kappa_{0}\sqrt{\sum_{p=1}^{d_{0}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial(XU^{(i)})}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}^{2}}=\kappa_{0}\sqrt{d_{i}}\left\|X\right\|_{{F}}.$
(28)

Going from layer $i$ to layer $j+1$, the derivative of $H^{(j+1)}$ with
respect to $U^{(i)}_{p,q}$ satisfies:

$\displaystyle\left\|\frac{\partial H^{(j+1)}}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}$
$\displaystyle=\left\|\phi^{\prime}_{j+1}\big{(}{XU^{(j+1)}+\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}\big{)}\odot\frac{\partial\big{(}{XU^{(j+1)}+\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}\big{)}}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}$
$\displaystyle\leq\kappa_{0}\left\|\frac{\partial\rho_{j+1}\big{(}{P_{{}_{G}}\psi_{j+1}(H^{(j)})W^{(j+1)}}\big{)}}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}$ (by equation (25) w.r.t. $\phi^{\prime}_{j+1}$)
$\displaystyle\leq\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\left\|\frac{\partial
H^{(j)}}{\partial U^{(i)}_{p,q}}\right\|_{{F}}.$ (by equation (25) w.r.t.
$\rho^{\prime}_{j+1},\psi^{\prime}_{j+1}$)

Hence, the Jacobian of $H^{(j+1)}$ with respect to $U^{(i)}$ satisfies:

$\displaystyle\left\|\frac{\partial H^{(j+1)}}{\partial
U^{(i)}}\right\|_{{F}}\leq\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\left\|\frac{\partial
H^{(j)}}{\partial U^{(i)}}\right\|_{{F}}.$

By repeatedly applying the above step beginning with the base case of equation
(28), we have proved that equation (24) holds. The proof of Proposition A.7 is
complete. ∎

####  A.2.3 Deal with second-order derivatives

In the second part towards showing Theorem 3.1 for MPNNs, we look at second-
order derivatives of the embeddings. This will appear later when we deal with
the trace of the Hessian. A fact that we will use throughout the proof is

$\displaystyle\left\|\phi^{\prime\prime}_{i}(x)\right\|_{\infty}\leq\kappa_{1},~{}~{}\left\|\psi^{\prime\prime}_{i}(x)\right\|_{\infty}\leq\kappa_{1},~{}~{}\text{and}~{}~{}\left\|\rho^{\prime\prime}_{i}(x)\right\|_{\infty}\leq\kappa_{1},$
(29)

for any $x$ and $i=1,\dots,l-1$. This is because
$\phi^{\prime}_{i},\psi_{i}^{\prime},$ and $\rho_{i}^{\prime}$ are all
$\kappa_{1}$-Lipschitz continuous from Assumption A.2.

######  Proposition A.8.

In the setting of Theorem 3.1, the second-order derivative of $H^{(l)}$ with
respect to $W^{(i)}$ and $U^{(i)}$ satisfies the following, for any
$i=1,\dots,l-1$ and any $j=i,\dots,l-1$:

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(j)}}{\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}$
$\displaystyle\leq
C_{i,j}\kappa_{1}d_{i}\max(\left\|P_{{}_{G}}\right\|^{j-i+2},\left\|P_{{}_{G}}\right\|^{2(j-i+1)})\left\|H^{(i-1)}\right\|_{{F}}^{2}\prod_{t=i+1}^{j}s_{t}^{2},$
(30)
$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(j)}}{\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}$
$\displaystyle\leq\hat{C}_{i,j}\kappa_{1}d_{i}\max(\left\|P_{{}_{G}}\right\|^{j-i},\left\|P_{{}_{G}}\right\|^{2(j-i)})\left\|X\right\|_{{F}}^{2}\prod_{t=i+1}^{j}s_{t}^{2},$
(31)

where $C_{i,j}$

$C_{i,j}=\left\\{\begin{aligned}
&\kappa_{0}^{3(j-i+1)}\frac{\kappa_{0}^{3(j-i)+2}-1}{\kappa_{0}-1},&\kappa_{0}\neq
1,\\\ &3(j-i)+2,&\kappa_{0}=1,\end{aligned}\right.$

and $\hat{C}_{i,j}$

$\hat{C}_{i,j}=\left\\{\begin{aligned}
&\kappa_{0}^{3(j-i)}\frac{\kappa_{0}^{3(j-i)+1}-1}{\kappa_{0}-1},&\kappa_{0}\neq
1,\\\ &3(j-i)+1,&\kappa_{0}=1.\end{aligned}\right.$

are fixed constants that depend on the Lipschitz-continuity of the activation
mappings.

###### Proof.

First, we will consider equation (30). To simplify the derivation, we
introduce two notations for brevity. Let

$\displaystyle F_{j}=P_{{}_{G}}\psi_{j}\big{(}{H^{(j-1)}}\big{)}W^{(j)}\text{
and }E_{j}=XU^{(j)}+\rho_{j}\big{(}{F_{j}}\big{)}.$

In the base case when $j=i$, from the first-order derivative in equation (26),
we use the chain rule to get:

$\displaystyle\frac{\partial^{2}H^{(i)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}=$
$\displaystyle\phi^{\prime\prime}_{i}(E_{i})\odot\frac{\partial
E_{i}}{\partial W_{p,q}^{(i)}}\odot\frac{\partial E_{i}}{\partial
W_{p,q}^{(i)}}+\phi^{\prime}_{i}(E_{i})\odot\rho^{\prime\prime}_{i}(F_{i})\odot{\frac{\partial
F_{i}}{\partial W_{p,q}^{(i)}}}\odot{\frac{\partial F_{i}}{\partial
W^{(i)}_{p,q}}}.$ (32)

Using equation (29), the maximum entries of
$\phi^{\prime\prime}_{i}(\cdot),\rho^{\prime\prime}_{i}(\cdot)$ are at most
$\kappa_{1}$. Using equation (25), the maximum entry of
$\phi^{\prime}_{i}(\cdot)$ is at most $\kappa_{0}$. Notice that the derivative
of $E_{i}$ can be reduced to the derivative of $F_{i}$ as follows:

$\displaystyle\left\|\frac{\partial E_{i}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}=\left\|\rho^{\prime}_{i}(F_{i})\odot\frac{\partial
F_{1}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}\leq\kappa_{0}^{2}\left\|\frac{\partial
F_{i}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}.$ (33)

Therefore, based on the conditions for first- and second-order derivatives
(cf. (25) and (29)), the Frobenius norm of the above equation (32) is at most:

$\displaystyle\left\|\frac{\partial^{2}H^{(i)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq\kappa_{1}\left\|\frac{\partial
E_{i}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}\kappa_{1}\left\|\frac{\partial
F_{i}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}\leq(\kappa_{0}+1)\kappa_{0}\kappa_{1}\left\|\frac{\partial
F_{i}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}.$

Notice that the derivative of $F_{i}$ with respect to $W^{(i)}_{p,q}$ is
nonzero only in the $q$-th column of $F_{i}$, and is equal to the $p$-th
column of $P_{{}_{G}}g_{i}(H^{(i-1)})$. Therefore, by summing over
$p=1,\dots,d_{i-1}$ and $q=1,\dots,d_{i}$, we get:

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
F_{i}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}^{2}\leq
d_{i}\left\|P_{{}_{G}}\psi_{i}(H^{(i-1)})\right\|_{{F}}^{2}.$

Therefore, we have derived the base case when $j=i$ as:

$\displaystyle\left\|\frac{\partial^{2}H^{(i)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq(\kappa_{0}+1)\kappa_{0}^{3}\kappa_{1}d_{i}\left\|P_{{}_{G}}\right\|^{2}\left\|H^{(i-1)}\right\|_{{F}}^{2}.$
(34)

Next, we consider the induction step from layer $j$ to layer $j+1$. This step
is similar to the base case but also differs since $H^{(j)}$ is now dependent
on $W^{(i)}$. Recall that the second-order derivatives satisfy equation (29).
Based on the Lipschitz-continuity conditions, the Frobenius norm of the
second-order derivatives satisfies:

$\displaystyle\left\|\frac{\partial^{2}H^{(j+1)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}$
(35) $\displaystyle\leq\kappa_{1}\left\|\frac{\partial E_{j+1}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}\kappa_{1}\left\|\frac{\partial
F_{j+1}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}^{2}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\Big{(}{\kappa_{1}\left\|\frac{\partial
H^{(j)}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}\left\|\frac{\partial^{2}H^{(j)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}}\Big{)}$
$\displaystyle\leq(\kappa_{0}+1)\kappa_{0}\kappa_{1}\left\|\frac{\partial
F_{j+1}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}^{2}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\Big{(}{\kappa_{1}\left\|\frac{H^{(j)}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}\left\|\frac{\partial^{2}H^{(j)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}}\Big{)}.$
(36)

The last step follows similarly as equation (33). For the derivative of
$F_{j+1}$, using the chain rule, we get:

$\displaystyle\left\|\frac{\partial F_{j+1}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}$
$\displaystyle=\left\|P_{{}_{G}}\frac{\partial\psi_{j+1}(H^{(j)})}{\partial
W_{p,q}^{(i)}}W^{(j+1)}\right\|_{{F}}^{2}$
$\displaystyle\leq\left\|P_{{}_{G}}\right\|^{2}\left\|W^{(j+1)}\right\|^{2}\left\|\frac{\partial\psi_{j+1}(H^{(j)})}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}$
$\displaystyle\leq\left\|P_{{}_{G}}\right\|^{2}\left\|W^{(j+1)}\right\|^{2}\left\|\psi^{\prime}_{j+1}(H^{(j)})\odot\frac{\partial
H^{(j)}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}$
$\displaystyle\leq\kappa_{0}^{2}\left\|P_{{}_{G}}\right\|^{2}\left\|W^{(j+1)}\right\|^{2}\left\|\frac{\partial
H^{(j)}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}.$

Therefore, combining the above with equations (36) together, we get the
following result:

$\displaystyle\left\|\frac{\partial^{2}H^{(j+1)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq\,$
$\displaystyle\Big{(}{(\kappa_{0}+1)\kappa_{0}^{3}\kappa_{1}\left\|P_{{}_{G}}\right\|^{2}\left\|W^{(j+1)}\right\|^{2}+\kappa_{0}^{2}\kappa_{1}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|}\Big{)}\left\|\frac{\partial
H^{(j)}}{\partial W_{p,q}^{(i)}}\right\|_{{F}}^{2}$
$\displaystyle+\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|\left\|W^{(j+1)}\right\|\left\|\frac{\partial^{2}H^{(j)}}{\partial\big{(}{W^{(i)}_{p,q}}\big{)}^{2}}\right\|_{{F}}$
$\displaystyle\leq$
$\displaystyle\max\Big{(}{\left\|P_{{}_{G}}\right\|,\left\|P_{{}_{G}}\right\|^{2}}\Big{)}s_{j+1}^{2}\Big{(}{(\kappa_{0}^{2}+\kappa_{0}+1)\kappa_{0}^{2}\kappa_{1}\left\|\frac{\partial
H^{(j)}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}+\kappa_{0}^{3}\left\|\frac{\partial^{2}H^{(j)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}}\Big{)}.$

Based on equation (23) of Proposition (A.7), the first-order derivative of
$H^{(j)}$ satisfies:

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial
H^{(j)}}{\partial
W_{p,q}^{(i)}}\right\|_{{F}}^{2}\leq\kappa_{0}^{6(j-i+1)}{d_{i}}\left\|P_{{}_{G}}\right\|^{2(j-i+1)}\left\|H^{(i-1)}\right\|^{2}\prod_{t=i+1}^{j}s_{t}^{2}.$
(37)

Applying equation (37) to the above (and summing over $p=1,\dots,d_{i-1}$ and
$q=1,\dots,d_{i}$) forms the induction step for showing equation (30):

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(j+1)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq\,$
$\displaystyle\frac{\kappa_{0}^{3}-1}{\kappa_{0}-1}\kappa_{0}^{6(j-i+1)+2}\kappa_{1}d_{i}\max\Big{(}{\left\|P_{{}_{G}}\right\|^{2(j-i)+3},\left\|P_{{}_{G}}\right\|^{2(j-i)+4}}\Big{)}\left\|H^{(i-1)}\right\|_{{F}}^{2}\prod_{t=i+1}^{j+1}s_{t}^{2}$
$\displaystyle+\kappa_{0}^{3}\max\big{(}{\left\|P_{{}_{G}}\right\|,\left\|P_{{}_{G}}\right\|^{2}}\big{)}s_{j+1}^{2}\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(j)}}{\partial\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}.$

By repeatedly applying the induction step along with the base case in equation
(34), we have shown that equation (30) holds:

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(j)}}{\big{(}{W_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq
C_{i,j}\kappa_{1}d_{i}\max(\left\|P_{{}_{G}}\right\|^{j-i+2},\left\|P_{{}_{G}}\right\|^{2(j-i+1)})\left\|H^{(i-1)}\right\|_{{F}}^{2}\prod_{t=i+1}^{j}s_{t}^{2},$
(38)

where $C_{i,j}$ satisfies the following equation:

$C_{i,j}=\left\\{\begin{aligned}
&\kappa_{0}^{3(j-i+1)}\frac{\kappa_{0}^{3(j-i)+2}-1}{\kappa_{0}-1},&\kappa_{0}\neq
1,\\\ &3(j-i)+2,&\kappa_{0}=1.\end{aligned}\right.$

In the second part of the proof, we consider equation (31) similar to the
first part. However, the analysis will be significantly simpler. We first
consider the base case. Similar to equation (32), the second-order derivative
of $H^{(i)}$ over $W^{(i)}_{p,q}$ satisfies, for any $p=1,\dots,d_{0}$ and
$q=1,\dots,d_{i}$:

$\displaystyle\left\|\frac{\partial^{2}H^{(i)}}{\partial\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}=\left\|\phi_{i}^{\prime\prime}(E_{i})\odot\frac{\partial
E_{i}}{\partial U^{(i)}_{p,q}}\odot\frac{\partial E_{i}}{\partial
U^{(i)}_{p,q}}\right\|\leq\kappa_{1}\left\|\frac{\partial(XU^{(i)})}{\partial
U^{(i)}_{p,q}}\right\|_{{F}}^{2}.$

Therefore, by summing up the above over all $p$ and $q$, we get the base case
result:

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(i)}}{\partial\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq\kappa_{1}d_{i}\left\|X\right\|_{{F}}^{2}.$
(39)

Next, we consider the induction step from layer $j$ to layer $j+1$. This step
follows the same analysis until equation (38), from which we can similarly
derive that:

$\displaystyle\sum_{p=1}^{d_{i-1}}\sum_{q=1}^{d_{i}}\left\|\frac{\partial^{2}H^{(j)}}{\big{(}{U_{p,q}^{(i)}}\big{)}^{2}}\right\|_{{F}}\leq\hat{C}_{i,j}\kappa_{1}d_{i}\max(\left\|P_{{}_{G}}\right\|^{j-i},\left\|P_{{}_{G}}\right\|^{2(j-i)})\left\|X\right\|_{{F}}^{2}\prod_{t=i+1}^{j}s_{t}^{2}.$
(40)

where $\hat{C}_{i,j}$ satisfies the following equation:

$\hat{C}_{i,j}=\left\\{\begin{aligned}
&\kappa_{0}^{3(j-i)}\frac{\kappa_{0}^{3(j-i)+1}-1}{\kappa_{0}-1},&\kappa_{0}\neq
1,\\\ &3(j-i)+1,&\kappa_{0}=1.\end{aligned}\right.$

∎

####  A.2.4 Proof of Theorem 3.1

Based on Propositions A.7 and A.8, we are ready to present the proof of
Theorem 3.1 for message passing GNNs. First, we will apply the bounds on the
derivatives back in Lemma 4.3. After getting the trace of the Hessians, we
then use the PAC-Bayes bound from Lemma 4.1 to complete the proof.

###### Proof of Theorem 3.1.

By applying equations (23) and (30) into Lemma 4.3’s result, we get that the
trace of $\mathbf{H}^{(l)}$ with respect to $W^{(i)}$ is less than:

$\displaystyle\left(\frac{\kappa_{0}\sqrt{k}}{\sqrt{n}}C_{i,l-1}\kappa_{1}d_{i}\max\big{(}\left\|P_{{}_{G}}\right\|^{l-i+1},\left\|P_{{}_{G}}\right\|^{2(l-i)}\big{)}+\frac{\kappa_{1}k}{n}\kappa_{0}^{6(l-i)}\kappa_{1}d_{i}\left\|P_{{}_{G}}\right\|^{2(l-i)}\right)\left\|H^{(i-1)}\right\|_{{F}}^{2}\prod_{t=i+1}^{l}s_{t}^{2}$
$\displaystyle\leq(\kappa_{0}C_{i,l-1}+\kappa_{0}^{6(l-i)})\sqrt{\frac{{k}}{{n}}}\kappa_{1}d_{i}\max\Big{(}{\left\|P_{{}_{G}}\right\|^{l-i+1},\left\|P_{{}_{G}}\right\|^{2(l-i)}}\Big{)}\left\|H^{(i-1)}\right\|_{{F}}^{2}\prod_{t=i+1}^{l}s_{t}^{2},$
(41)

for any $i=1,2,\cdots,l-1$. Here we have

$\kappa_{0}C_{i,l-1}+\kappa_{0}^{6(l-i)}=\left\\{\begin{aligned}
&\kappa_{0}^{3(l-i)+1}\frac{\kappa_{0}^{3(l-i)}-1}{\kappa_{0}-1},&\kappa_{0}\neq
1,\\\ &3(l-i)-1,&\kappa_{0}=1.\end{aligned}\right.$

It remains to consider the Frobenius norm of $H^{(i-1)}$. Notice that this
satisfies the following:

$\displaystyle\left\|H^{(i-1)}\right\|_{{F}}\leq\kappa_{0}\left\|XU^{(i-1)}+\rho_{i-1}(P_{{}_{G}}\psi_{i-1}(H^{(i-2)}))W^{(i-1)}\right\|_{{F}}$
$\displaystyle\leq\kappa_{0}\left\|U^{(i-1)}\right\|\left\|X\right\|_{{F}}+\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|\left\|W^{(i-1)}\right\|\left\|H^{(i-2)}\right\|_{{F}}\leq\kappa_{0}s_{i}\left\|X\right\|_{{F}}+\kappa_{0}^{3}\left\|P_{{}_{G}}\right\|s_{i}\left\|H^{(i-2)}\right\|_{{F}}.$

By induction over $i$ for the above step, we get that the Frobenius norm of
$H^{(i-1)}$ must be less than:

$\displaystyle\big{(}{\kappa_{0}^{3(i-1)}+\sum_{j=0}^{i-2}\kappa_{0}^{3j+1}}\big{)}{\sqrt{k}}\max_{(X,G,y)\sim\mathcal{D}}\left\|X\right\|\max\big{(}{1,\left\|P_{{}_{G}}\right\|^{i-1}}\big{)}\prod_{j=1}^{i-1}s_{j}.$
(42)

By applying the above (42) back in (41), we have shown that the trace of
$\mathbf{H}^{(l)}$ with respect to $W^{(i)}$ is less than:

$\displaystyle
C^{\prime}\max_{(X,G,y)\sim\mathcal{D}}\left\|X\right\|^{2}\kappa_{1}d_{i}k\max(1,\|P_{{}_{G}}\|^{2(l-1)})\prod_{t=1:t\neq
i}^{l}s_{t}^{2},$ (43)

where $C^{\prime}$ satisfies the following equation:

$C^{\prime}=\left\\{\begin{aligned}
&\frac{(\kappa_{0}^{3l}-1)(\kappa_{0}^{3(l-1)/2}-1)^{2}}{(\kappa_{0}-1)^{3}},&\kappa_{0}\neq
1,\\\ &\frac{4}{9}l^{3},&\kappa_{0}=1.\end{aligned}\right.$

To be specific, when $\kappa_{0}=1$, $(3(l-i)-1)i^{2}\leq\frac{4}{9}l^{3}$ .
If $\kappa_{0}\neq 1$ and $i\geq 2$, we have

$\displaystyle\Big{(}{\kappa_{0}^{3(l-i)+1}\frac{\kappa_{0}^{3(l-i)}-1}{\kappa_{0}-1}}\Big{)}\big{(}{\kappa_{0}^{3(i-1)}+\sum_{j=0}^{i-2}\kappa_{0}^{3j+1}}\big{)}^{2}\leq\kappa_{0}^{3(l-i)+3}\frac{\kappa_{0}^{3(l-i)}-1}{\kappa_{0}-1}\frac{(\kappa_{0}^{3(i-1)}-1)^{2}}{(\kappa_{0}-1)^{2}}$
$\displaystyle=\frac{\kappa_{0}^{3l}-\kappa_{0}^{3(l-i+1)}}{(\kappa_{0}-1)^{3}}\Big{(}{(\kappa_{0}^{3(l-i)}-1)(\kappa_{0}^{3(i-1)}-1)}\Big{)}\leq\frac{(\kappa_{0}^{3l}-1)(\kappa_{0}^{3(l-1)/2}-1)^{2}}{(\kappa_{0}-1)^{3}}.$

If $\kappa_{0}\neq 1$ and $i=1$, we obtain

$\displaystyle\Big{(}{\kappa_{0}^{3(l-i)+1}\frac{\kappa_{0}^{3(l-i)}-1}{\kappa_{0}-1}}\Big{)}\Big{(}{\kappa_{0}^{3(i-1)}+\sum_{j=0}^{i-2}\kappa_{0}^{3j+1}}\Big{)}^{2}=\kappa_{0}^{3l-2}\frac{\kappa_{0}^{3(l-1)}-1}{\kappa_{0}-1}\leq\frac{(\kappa_{0}^{3l}-1)(\kappa_{0}^{3(l-1)/2}-1)^{2}}{(\kappa_{0}-1)^{3}}.$

The above works for the layers from the beginning until layer $l-1$. Last, we
consider the trace of $\mathbf{H}^{(l)}$ with respect to $W^{(l)}$ (notice
that $\bm{U}$ is not needed in the readout layer). Similar to equation (19),
one can prove that the trace of the Hessian with respect to $W^{(l)}$
satisfies:

$\displaystyle{\left|{\operatorname{Tr}\big{[}\mathbf{H}^{(l)}[\ell(H^{(l)},y)]\big{]}}\right|}\leq\,$
$\displaystyle\kappa_{0}\sqrt{k}\sum_{p=1}^{d_{l-1}}\sum_{q=1}^{d_{l}}\left\|\frac{\partial^{2}H^{(l)}}{\partial\big{(}W^{(l)}_{p,q}\big{)}^{2}}\right\|+\kappa_{1}k\sum_{p=1}^{d_{l-1}}\sum_{q=1}^{d_{l}}\left\|\frac{\partial
H^{(l)}}{\partial W^{(l)}_{p,q}}\right\|^{2}$ $\displaystyle\leq\,$
$\displaystyle\kappa_{0}\sqrt{k}\sum_{p=1}^{d_{l-1}}\sum_{q=1}^{d_{l}}\left\|\frac{1}{n}\bm{1}_{n}^{\top}H^{(l-1)}\frac{\partial^{2}W^{(l)}}{\partial\big{(}W^{(l)}_{p,q}\big{)}^{2}}\right\|+\kappa_{1}k\sum_{p=1}^{d_{l-1}}\sum_{q=1}^{d_{l}}\left\|\frac{1}{n}\bm{1}_{n}^{\top}H^{(l-1)}\frac{\partial
W^{(l)}}{\partial W^{(l)}_{p,q}}\right\|^{2}$ $\displaystyle\leq\,$
$\displaystyle\kappa_{1}k\sum_{p=1}^{d_{l-1}}\sum_{q=1}^{d_{l}}\left\|\frac{1}{n}\bm{1}_{n}\right\|^{2}\left\|H^{(l-1)}\frac{\partial
W^{(l)}}{\partial W^{(l)}_{p,q}}\right\|^{2}$ $\displaystyle=\,$
$\displaystyle\kappa_{1}\frac{k}{n}{d_{l}\left\|H^{(l-1)}\right\|_{{F}}^{2}}$

By equation (42), the above is bounded by

$\displaystyle\kappa_{1}\frac{k}{n}d_{l}\big{(}{\kappa_{0}^{3(l-1)}+\sum_{j=0}^{l-2}\kappa_{0}^{3j+1}}\big{)}^{2}\max_{(X,G,y)\sim\mathcal{D}}\left\|X\right\|_{{F}}^{2}\max\big{(}{1,\left\|P_{{}_{G}}\right\|^{2(l-1)}}\big{)}\prod_{j=1}^{l-1}s_{j}^{2}$
$\displaystyle\leq~{}$ $\displaystyle
C_{l}\max_{(X,G,y)\sim\mathcal{D}}\left\|X\right\|^{2}\kappa_{1}d_{l}k\max\Big{(}{1,\left\|P_{{}_{G}}\right\|^{2(l-1)}}\Big{)}\prod_{t=1:t\neq
l}^{l}s_{t}^{2},$

since $\frac{\left\|X\right\|_{{F}}^{2}}{n}\leq\left\|X\right\|^{2}$, where
$C_{l}$ satisfies the following equation:

$C_{l}=\left\\{\begin{aligned}
&\kappa_{0}^{2}\frac{(\kappa_{0}^{3(l-1)}-1)^{2}}{(\kappa_{0}-1)^{2}},&\kappa_{0}\neq
1,\\\ &l^{2},&\kappa_{0}=1.\end{aligned}\right.$

Finally, let

$\displaystyle\tilde{C}$ $\displaystyle=\max(C^{\prime},C_{l}).$ (44)

From the value of $C^{\prime}$ above and the value of $C_{l}$, we get that
$\tilde{C}$ is equal to

$\displaystyle\tilde{C}$ $\displaystyle=\left\\{\begin{aligned}
&\frac{(\kappa_{0}^{3l}-1)(\kappa_{0}^{3(l-1)/2}-1)^{2}}{(\kappa_{0}-1)^{3}},&\kappa_{0}\neq
1,\\\ &\frac{1}{2}l^{3},&\kappa_{0}=1.\end{aligned}\right.$

Similarly by applying equations (24) and (31) into Lemma 4.3, the trace of