# A PAC-Bayesian Approach to Generalization Bounds for Graph Neural Networks

Renjie Liao1,2,3, Raquel Urtasun1,2,3,4, Richard Zemel1,3,4
University of Toronto1, Uber ATG2, Vector Institute3, Canadian Institute for
Advanced Research4
{rjliao, urtasun<EMAIL_ADDRESS>
###### Abstract

In this paper, we derive generalization bounds for the two primary classes of
graph neural networks (GNNs), namely graph convolutional networks (GCNs) and
message passing GNNs (MPGNNs), via a PAC-Bayesian approach. Our result reveals
that the maximum node degree and spectral norm of the weights govern the
generalization bounds of both models. We also show that our bound for GCNs is
a natural generalization of the results developed in (Neyshabur et al., 2017)
for fully-connected and convolutional neural networks. For message passing
GNNs, our PAC-Bayes bound improves over the Rademacher complexity based bound
in (Garg et al., 2020), showing a tighter dependency on the maximum node
degree and the maximum hidden dimension. The key ingredients of our proofs are
a perturbation analysis of GNNs and the generalization of PAC-Bayes analysis
to non-homogeneous GNNs. We perform an empirical study on several real-world
graph datasets and verify that our PAC-Bayes bound is tighter than others.

##  1 Introduction

Graph neural networks (GNNs) (Gori et al., 2005; Scarselli et al., 2008;
Bronstein et al., 2017; Battaglia et al., 2018) have become very popular
recently due to their ability to learn powerful representations from graph-
structured data, and have achieved state-of-the-art results in a variety of
application domains such as social networks (Hamilton et al., 2017), quantum
chemistry (Gilmer et al., 2017), computer vision (Monti et al., 2017),
reinforcement learning (Sanchez-Gonzalez et al., 2018), robotics (Casas et
al., 2019), and physics (Henrion et al., 2017). Given a graph along with
node/edge features, GNNs learn node/edge representations by propagating
information on the graph via local computations shared across the nodes/edges.
Based on the specific form of local computation employed, GNNs can be divided
into two categories: graph convolution based GNNs (Bruna et al., 2013;
Duvenaud et al., 2015; Kipf & Welling, 2016) and message passing based GNNs
(Li et al., 2015; Dai et al., 2016; Gilmer et al., 2017). The former
generalizes the convolution operator from regular graphs (_e.g_., grids) to
ones with arbitrary topology, whereas the latter mimics message passing
algorithms and parameterizes the shared functions via neural networks.

Due to the tremendous empirical success of GNNs, there is increasing interest
in understanding their theoretical properties. For example, some recent works
study their expressiveness (Maron et al., 2018; Xu et al., 2018; Chen et al.,
2019), that is, what class of functions can be represented by GNNs. However,
only few works investigate why GNNs generalize so well to unseen graphs. They
are either restricted to a specific model variant (Verma & Zhang, 2019; Du et
al., 2019; Garg et al., 2020) or have loose dependencies on graph statistics
(Scarselli et al., 2018).

On the other hand, GNNs have close ties to standard feedforward neural
networks, _e.g_., multi-layer perceptrons (MLPs) and convolutional neural
networks (CNNs). In particular, if each i.i.d. sample is viewed as a node,
then the whole dataset becomes a graph without edges. Therefore, GNNs can be
seen as generalizations of MLPs/CNNs since they model not only the
regularities within a sample but also the dependencies among samples as
defined in the graph. It is therefore natural to ask if we can generalize the
recent advancements on generalization bounds for MLPs/CNNs (Harvey et al.,
2017; Neyshabur et al., 2017; Bartlett et al., 2017; Dziugaite & Roy, 2017;
Arora et al., 2018, 2019) to GNNs, and how would graph structures affect the
generalization bounds?

In this paper, we answer the above questions by proving generalization bounds
for the two primary classes of GNNs, _i.e_., graph convolutional networks
(GCNs) (Kipf & Welling, 2016) and message-passing GNNs (MPGNNs) (Dai et al.,
2016; Jin et al., 2018).

Our generalization bound for GCNs shows an intimate relationship with the
bounds for MLPs/CNNs with ReLU activations (Neyshabur et al., 2017; Bartlett
et al., 2017). In particular, they share the same term, _i.e_., the product of
the spectral norms of the learned weights at each layer multiplied by a factor
that is additive across layers. The bound for GCNs has an additional
multiplicative factor $d^{{(l-1)}/{2}}$ where $d-1$ is the maximum node degree
and $l$ is the network depth. Since MLPs/CNNs are special GNNs operating on
graphs without edges (_i.e_., $d-1=0$), the bound for GCNs coincides with the
ones for MLPs/CNNs with ReLU activations (Neyshabur et al., 2017) on such
degenerated graphs. Therefore, our result is a natural generalization of the
existing results for MLPs/CNNs.

Our generalization bound for message passing GNNs reveals that the governing
terms of the bound are similar to the ones of GCNs, _i.e_., the geometric
series of the learned weights and the multiplicative factor $d^{l-1}$. The
geometric series appears due to the weight sharing across message passing
steps, thus corresponding to the product term across layers in GCNs. The term
$d^{l-1}$ encodes the key graph statistics. Our bound improves the dependency
on the maximum node degree and the maximum hidden dimension compared to the
recent Rademacher complexity based bound (Garg et al., 2020). Moreover, we
compute the bound values on four real-world graph datasets (_e.g_., social
networks and protein structures) and verify that our bounds are tighter.

In terms of the proof techniques, our analysis follows the PAC-Bayes framework
in the seminal work of (Neyshabur et al., 2017) for MLPs/CNNs with ReLU
activations. However, we make two distinctive contributions which are
customized for GNNs. First, a naive adaptation of the perturbation analysis in
(Neyshabur et al., 2017) does not work for GNNs since ReLU is not 1-Lipschitz
under the spectral norm, _i.e_., $\|\text{ReLU}(X)\|_{2}\leq\|X\|_{2}$ does
not hold for some real matrix $X$. Instead, we construct the recursion on
certain node representations of GNNs like the one with maximum $\ell_{2}$
norm, so that we can perform perturbation analysis with vector 2-norm. Second,
in contrast to (Neyshabur et al., 2017) which only handles the homogeneous
networks, _i.e_., $f(ax)=af(x)$ when $a\geq 0$, we properly construct a
quantity of the learned weights which 1) provides a way to satisfy the
constraints of the previous perturbation analysis and 2) induces a finite
covering on the range of the quantity so that the PAC-Bayes bound holds for
all possible weights. This generalizes the analysis to non-homogeneous GNNs
like typical MPGNNs.

The rest of the paper is organized as follows. In Section 2, we introduce
background material necessary for our analysis. We then present our
generalization bounds and the comparison to existing results in Section 3. We
also provide an empirical study to support our theoretical arguments in
Section 4. At last, we discuss the extensions, limitations and some open
problems.

##  2 Background

In this section, we first explain our analysis setup including notation and
assumptions. We then describe the two representative GNN models in detail.
Finally, we review the PAC-Bayes analysis.

###  2.1 Analysis Setup

In the following analysis, we consider the $K$-class graph classification
problem which is common in the GNN literature, where given a graph sample $z$,
we would like to classify it into one of the predefined $K$ classes. We will
discuss extensions to other problems like graph regression in Section 5. Each
graph sample $z$ is a triplet of an adjacency matrix $A$, node features
$X\in\mathbb{R}^{n\times h_{0}}$ and output label $y\in\mathbb{R}^{1\times
K}$, _i.e_. $z=(A,X,y)$, where $n$ is the number of nodes and $h_{0}$ is the
input feature dimension. We start our discussion by defining our notations.
Let $\mathbb{N}_{k}^{+}$ be the first $k$ positive integers, _i.e_.,
$\mathbb{N}_{k}^{+}=\\{1,2,\dots,k\\}$, $|\cdot|_{p}$ the vector $p$-norm and
$\|\cdot\|_{p}$ the operator norm induced by the vector $p$-norm. Further,
$\|\cdot\|_{F}$ denotes the Frobenius norm of a matrix, $e$ the base of the
natural logarithm function $\log$, $A[i,j]$ the $(i,j)$-th element of matrix
$A$ and $A[i,:]$ the $i$-th row. We use parenthesis to avoid the ambiguity,
_e.g_., $(AB)[i,j]$ means the $(i,j)$-th element of the product matrix $AB$.
We then introduce some terminologies from statistical learning theory and
define the sample space as $\mathcal{Z}$, $z=(A,X,y)\in\mathcal{Z}$ where
$X\in\mathcal{X}$ (node feature space) and $A\in\mathcal{G}$ (graph space),
data distribution $\mathcal{D}$, $z\overset{iid}{\sim}\mathcal{D}$, hypothesis
(or model) $f_{w}$ where $f_{w}\in\mathcal{H}$ (hypothesis class), and
training set $S$ with size $m$, $S=\\{z_{1},\dots,z_{m}\\}$. We make the
following assumptions which also appear in the literature:

  1. A1

Data, _i.e_., triplets $(A,X,y)$, are i.i.d. samples drawn from some unknown
distribution $\mathcal{D}$.

  2. A2

The maximum hidden dimension across all layers is $h$.

  3. A3

Node feature of any graph is contained in a $\ell_{2}$-ball with radius $B$.
Specifically, we have $\forall i\in\mathbb{N}^{+}_{n}$, the $i$-th node
feature
$X[i,:]\in\mathcal{X}_{B,h_{0}}=\\{x\in\mathbb{R}^{h_{0}}|\sum_{j=1}^{h_{0}}x_{j}^{2}\leq
B^{2}\\}$.

  4. A4

We only consider simple graphs (_i.e_., undirected, no loops111Here loop means
an edge that connects a vertex to itself, a.k.a., self-loop., and no multi-
edges) with maximum node degree as $d-1$.

Note that it is straightforward to estimate $B$ and $d$ empirically on real-
world graph data.

###  2.2 Graph Neural Networks (GNNs)

In this part, we describe the details of the GNN models and the loss function
we used for the graph classification problem. The essential idea of GNNs is to
propagate information over the graph so that the learned representations
capture the dependencies among nodes/edges. We now review two classes of GNNs,
GCNs and MPGNNs, which have different mechanisms for propagating information.
We choose them since they are the most popular variants and represent two
common types of neural networks, _i.e_., feedforward (GCNs) and recurrent
(MPGNNs) neural networks. We discuss the extension of our analysis to other
GNN variants in Section 5. For ease of notation, we define the model to be
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$
where $w$ is the vectorization of all model parameters.

##### GCNs:

Graph convolutional networks (GCNs) (Kipf & Welling, 2016) for the $K$-class
graph classification problem can be defined as follows,

$\displaystyle H_{k}$
$\displaystyle=\sigma_{k}\left(\tilde{L}H_{k-1}W_{k}\right)\quad$
$\displaystyle(k\text{-th Graph Convolution Layer})$ $\displaystyle H_{l}$
$\displaystyle=\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}\quad$
$\displaystyle(\text{Readout Layer}),$ (1)

where $k\in\mathbb{N}^{+}_{l-1}$, $H_{k}\in\mathbb{R}^{n\times h_{k}}$ are the
node representations/states, $\bm{1}_{n}\in\mathbb{R}^{1\times n}$ is a all-
one vector, $l$ is the number of layers.222We count the readout function as a
layer to be consistent with the existing analysis of MLPs/CNNs. and $W_{j}$ is
the weight matrix of the $j$-th layer. The initial node state is the observed
node feature $H_{0}=X$. For both GCNs and MPGNNs, we consider $l>1$ since
otherwise the model degenerates to a linear transformation which does not
leverage the graph and is trivial to analyze. Due to assumption A2, $W_{j}$ is
of size at most $h\times h$, _i.e_., $h_{k}\leq h$, $\forall
k\in\mathbb{N}^{+}_{l-1}$. The graph Laplacian $\tilde{L}$ is defined as,
$\tilde{A}=I+A$, $\tilde{L}=D^{-\frac{1}{2}}\tilde{A}D^{-\frac{1}{2}}$ where
$D$ is the degree matrix of $\tilde{A}$. Note that the maximum eigenvalue of
$\tilde{L}$ is $1$ in this case. We absorb the bias into the weight by
appending constant $1$ to the node feature. Typically, GCNs use ReLU as the
non-linearity, _i.e_., $\sigma_{i}(x)=\max(0,x),\forall i=1,\cdots,l-1$. We
use the common mean-readout to obtain the graph representation where
$H_{l-1}\in\mathbb{R}^{n\times h_{l-1}}$, $W_{l}\in\mathbb{R}^{h_{l-1}\times
K}$, and $H_{l}\in\mathbb{R}^{1\times K}$.

##### MPGNNs:

There are multiple variants of message passing GNNs, _e.g_., (Li et al., 2015;
Dai et al., 2016; Gilmer et al., 2017), which share the same algorithmic
framework but instantiate a few components differently, _e.g_., the node state
update function. We choose the same class of models as in (Garg et al., 2020)
which are popular in the literature (Dai et al., 2016; Jin et al., 2018) in
order to fairly compare bounds. This MPGNN model can be written in matrix
forms as follows,

$\displaystyle M_{k}$
$\displaystyle=g(C_{\operatorname{out}}^{\top}H_{k-1})\quad$
$\displaystyle(k\text{-th step Message Computation})$
$\displaystyle\bar{M}_{k}$ $\displaystyle=C_{\operatorname{in}}M_{k}\quad$
$\displaystyle(k\text{-th step Message Aggregation})$ $\displaystyle H_{k}$
$\displaystyle=\phi\left(XW_{1}+\rho\left(\bar{M}_{k}\right)W_{2}\right)\quad$
$\displaystyle(k\text{-th step Node State Update})$ $\displaystyle H_{l}$
$\displaystyle=\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}\quad$
$\displaystyle(\text{Readout Layer}),$ (2)

where $k\in\mathbb{N}^{+}_{l-1}$, $H_{k}\in\mathbb{R}^{n\times h_{k}}$ are
node representations/states and $H_{l}\in\mathbb{R}^{1\times K}$ is the output
representation. Here we initialize $H_{0}=\bm{0}$. W.l.o.g., we assume
$\forall k\in\mathbb{N}^{+}_{l-1}$, $H_{k}\in\mathbb{R}^{n\times h}$ and
$M_{k}\in\mathbb{R}^{n\times h}$ since $h$ is the maximum hidden dimension.
$C_{\operatorname{in}}\in\mathbb{R}^{n\times c}$ and
$C_{\operatorname{out}}\in\mathbb{R}^{n\times c}$ ($c$ is the number of edges)
are the incidence matrices corresponding to incoming and outgoing nodes333 For
undirected graphs, we convert each edge into two directed edges.
respectively. Specifically, rows and columns of $C_{\operatorname{in}}$ and
$C_{\operatorname{out}}$ correspond to nodes and edges respectively.
$C_{\operatorname{in}}[i,j]=1$ indicates that the incoming node of the $j$-th
edge is the $i$-th node. Similarly, $C_{\operatorname{out}}[i,j]=1$ indicates
that the outgoing node of the $j$-th edge is the $i$-th node. $g,\phi,\rho$
are nonlinear mappings, _e.g_., ReLU and Tanh. Technically speaking,
$g:\mathbb{R}^{h}\rightarrow\mathbb{R}^{h}$,
$\phi:\mathbb{R}^{h}\rightarrow\mathbb{R}^{h}$, and
$\rho:\mathbb{R}^{h}\rightarrow\mathbb{R}^{h}$ operate on vector-states of
individual node/edge. However, since we share these functions across
nodes/edges, we can naturally generalize them to matrix-states, _e.g_.,
$\tilde{\phi}:\mathbb{R}^{n\times h}\rightarrow\mathbb{R}^{n\times h}$ where
$\tilde{\phi}(X)[i,:]=\phi(X[i,:])$. By doing so, the same function could be
applied to matrices with varying size of the first dimension. For simplicity,
we use $g,\phi,\rho$ to denote such generalization to matrices. We denote the
Lipschitz constants of $g,\phi,\rho$ under the vector $2$-norm as
$C_{g},C_{\phi},C_{\rho}$ respectively. We also assume $g(\bm{0})=\bm{0}$,
$\phi(\bm{0})=\bm{0}$, and $\rho(\bm{0})=\bm{0}$ and define the _percolation
complexity_ as $\mathcal{C}=C_{g}C_{\phi}C_{\rho}\|W_{2}\|_{2}$ following
(Garg et al., 2020).

##### Multiclass Margin Loss:

We use the multi-class $\gamma$-margin loss following (Bartlett et al., 2017;
Neyshabur et al., 2017). The generalization error is defined as,

$\displaystyle
L_{\mathcal{D},\gamma}(f_{w})=\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{w}(X,A)[y]\leq\gamma+\max_{j\neq
y}f_{w}(X,A)[j]\right),$ (3)

where $\gamma>0$ and $f_{w}(X,A)$ is the $l$-th layer representations, _i.e_.,
$H_{l}=f_{w}(X,A)$. Accordingly, we can define the empirical error as,

$\displaystyle L_{S,\gamma}(f_{w})=\frac{1}{m}\sum_{z_{i}\in
S}\bm{1}\left(f_{w}(X,A)[y]\leq\gamma+\max_{j\neq y}f_{w}(X,A)[j]\right).$ (4)

###  2.3 Background of PAC-Bayes Analysis

PAC-Bayes (McAllester, 1999, 2003; Langford & Shawe-Taylor, 2003) takes a
Bayesian view of the probably approximately correct (PAC) learning theory
(Valiant, 1984). In particular, it assumes that we have a prior distribution
$P$ over the hypothesis class $\mathcal{H}$ and obtain a posterior
distribution $Q$ over the same support through the learning process on the
training set. Therefore, instead of having a deterministic model/hypothesis as
in common learning formulations, we have a distribution of models. Under this
Bayesian view, we define the generalization error and the empirical error as,

$\displaystyle L_{S,\gamma}(Q)=\mathbb{E}_{w\sim
Q}[L_{S,\gamma}(f_{w})],\qquad\qquad
L_{\mathcal{D},\gamma}(Q)=\mathbb{E}_{w\sim
Q}[L_{\mathcal{D},\gamma}(f_{w})].$

Since many interesting models like neural networks are deterministic and the
exact form of the posterior $Q$ induced by the learning process and the prior
$P$ is typically unknown, it is unclear how one can perform PAC-Bayes
analysis. Fortunately, we can exploit the following result from the PAC-Bayes
theory.

######  Theorem 2.1.

(McAllester, 2003) (Two-sided) Let $P$ be a prior distribution over
$\mathcal{H}$ and let $\delta\in(0,1)$. Then, with probability $1-\delta$ over
the choice of an i.i.d. size-$m$ training set $S$ according to $\mathcal{D}$,
for all distributions $Q$ over $\mathcal{H}$ and any $\gamma>0$, we have

$\displaystyle L_{\mathcal{D},\gamma}(Q)\leq
L_{S,\gamma}(Q)+\sqrt{\frac{D_{\mathrm{KL}}(Q\|P)+\ln\frac{2m}{\delta}}{2(m-1)}}.$

Here $D_{\mathrm{KL}}$ is the KL-divergence. The nice thing about this result
is that the inequality holds for all possible prior $P$ and posterior $Q$
distributions. Hence, we have the freedom to construct specific priors and
posteriors so that we can work out the bound. Moreover, McAllester (2003);
Neyshabur et al. (2017) provide a general recipe to construct the posterior
such that for a large class of models, including deterministic ones, the PAC-
Bayes bound can be computed. Taking a neural network as an example, we can
choose a prior distribution with some known density, _e.g_., a fixed Gaussian,
over the initial weights. After the learning process, we can add random
perturbations to the learned weights from another known distribution as long
as the KL-divergence permits an analytical form. This converts the
deterministic model into a distribution of models while still obtaining a
tractable KL divergence. Leveraging Theorem 2.1 and the above recipe,
Neyshabur et al. (2017) obtained the following result which holds for a large
class of deterministic models.

######  Lemma 2.2.

(Neyshabur et al., 2017)444The constants slightly differ from the original
paper since we use a two-sided version of Theorem 2.1. Let
$f_{w}(x):\mathcal{X}\rightarrow\mathbb{R}^{K}$ be any model with parameters
$w$, and let $P$ be any distribution on the parameters that is independent of
the training data. For any $w$, we construct a posterior $Q(w+u)$ by adding
any random perturbation $u$ to $w$, _s.t_.,
$\mathbb{P}(\max_{x\in\mathcal{X}}|f_{w+u}(x)-f_{w}(x)|_{\infty}<\frac{\gamma}{4})>\frac{1}{2}$.
Then, for any $\gamma,\delta>0$, with probability at least $1-\delta$ over an
i.i.d. size-$m$ training set $S$ according to $\mathcal{D}$, for any $w$, we
have:

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\sqrt{\frac{2D_{\mathrm{KL}}(Q(w+u)\|P)+\log\frac{8m}{\delta}}{2(m-1)}}.$

This lemma guarantees that, as long as the change of the output brought by the
perturbations is small with a large probability, one can obtain the
corresponding generalization bound.

##  3 Generalization Bounds

In this section, we present the main results: generalization bounds of GCNs
and MPGNNs using a PAC-Bayesian approach. We then relate them to existing
generalization bounds of GNNs and draw connections to the bounds of MLPs/CNNs.
We summarize the key ideas of the proof in the main text and defer the details
to the appendix.

###  3.1 PAC-Bayes Bounds of GCNs

As discussed above, in order to apply Lemma 2.2, we must ensure that the
change of the output brought by the weight perturbations is small with a large
probability. In the following lemma, we bound this change using the product of
the spectral norms of learned weights at each layer and a term depending on
some statistics of the graph.

######  Lemma 3.1.

(GCN Perturbation Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-layer GCN. Then for any $w$, and $x\in\mathcal{X}_{B,h_{0}}$, and any
perturbation $u=\text{vec}(\\{U_{i}\\}_{i=1}^{l})$ such that $\forall
i\in\mathbb{N}^{+}_{l}$, $\|U_{i}\|_{2}\leq\frac{1}{l}\|W_{i}\|_{2}$, the
change in the output of GCN is bounded as,

$\displaystyle\left|f_{w+u}(X,A)-f_{w}(X,A)\right|_{2}\leq
eBd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\sum_{k=1}^{l}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}.$

The key idea of the proof is to decompose the change of the network output
into two terms which depend on two quantities of GNNs respectively: the
maximum change of node representations
$\max_{i}\left|H_{l-1}^{\prime}[i,:]-H_{l-1}[i,:]\right|_{2}$ and the maximum
node representation $\max_{i}\left|H_{l-1}[i,:]\right|_{2}$. Here the
superscript prime denotes the perturbed model. These two terms can be bounded
by an induction on the layer index. From this lemma, we can see that the most
important graph statistic for the stability of GCNs is the maximum node
degree, _i.e_., $d-1$. Armed with Lemma 3.1 and Lemma 2.2, we now present the
PAC-Bayes generalization bound of GCNs as Theorem 3.2.

######  Theorem 3.2.

(GCN Generalization Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$ layer GCN. Then for any $\delta,\gamma>0$, with probability at least
$1-\delta$ over the choice of an i.i.d. size-$m$ training set $S$ according to
$\mathcal{D}$, for any $w$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}d^{l-1}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}({\|W_{i}\|_{F}^{2}}/{\|W_{i}\|_{2}^{2}})+\log\frac{ml}{\delta}}{\gamma^{2}m}}\right).$

Since it is easy to show GCNs are homogeneous, the proof of Theorem 3.2
follows the one for MLPs/CNNs with ReLU activations in (Neyshabur et al.,
2017). In particular, we choose the prior distribution $P$ and the
perturbation distribution to be zero-mean Gaussians with the same diagonal
variance $\sigma$. The key steps of the proof are: 1) constructing a quantity
of learned weights $\beta=(\prod_{i=1}^{l}\|W_{i}\|_{2})^{1/l}$; 2) fixing any
$\tilde{\beta}$, considering all $\beta$ that are in the range
$|\beta-\tilde{\beta}|\leq\beta/l$ and choosing $\sigma$ which depends on
$\tilde{\beta}$ so that one can apply Lemma 3.1 and 2.2 to obtain the PAC-
Bayes bound; 3) taking a union bound of the result in the 2nd step by
considering multiple choices of $\tilde{\beta}$ so that all possible values of
$\beta$ (corresponding to all possible weight $w$) are covered. Although Lemma
2.2 and 3.1 have their own constraints on the random perturbation, above steps
provide a way to set the variance $\sigma$ which satisfies these constraints
and the independence w.r.t. learned weights. The latter is important since
$\sigma$ is also the variance of the prior $P$ which should not depend on
data.

###  3.2 PAC-Bayes Bounds of MPGNNs

For MPGNNs, we again need to perform a perturbation analysis to make sure that
the change of the network output brought by the perturbations on weights is
small with a large probability. Following the same strategy adopted in proving
Lemma 3.1, we prove the following Lemma.

######  Lemma 3.3.

(MPGNN Perturbation Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-step MPGNN. Then for any $w$, and $x\in\mathcal{X}_{B,h_{0}}$, and any
perturbation $u=\text{vec}(\\{U_{1},U_{2},U_{l}\\})$ such that
$\eta=\max\left(\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}},\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}},\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\right)\leq\frac{1}{l}$,
the change in the output of MPGNN is bounded as,

$\displaystyle|f_{w+u}($ $\displaystyle X,A)-f_{w}(X,A)|_{2}\leq
eBl\eta\|W_{1}\|_{2}\|W_{l}\|_{2}C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1},$

where $\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$.

The proof again involves decomposing the change into two terms which depend on
two quantities respectively: the maximum change of node representations
$\max_{i}\left|H_{l-1}^{\prime}[i,:]-H_{l-1}[i,:]\right|_{2}$ and the maximum
node representation $\max_{i}\left|H_{l-1}[i,:]\right|_{2}$. Then we perform
an induction on the layer index to obtain their bounds individually. Due to
the weight sharing across steps, we have a form of geometric series
$((d\mathcal{C})^{l-1}-1)/(d\mathcal{C}-1)$ rather than the product of
spectral norms of each layer as in GCNs. Technically speaking, the above lemma
only works with $d\mathcal{C}\neq 1$. We refer the reader to the appendix for
the special case of $d\mathcal{C}=1$. We now provide the generalization bound
for MPGNNs.

######  Theorem 3.4.

(MPGNN Generalization Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-step MPGNN. Then for any $\delta,\gamma>0$, with probability at least
$1-\delta$ over the choice of an i.i.d. size-$m$ training set $S$ according to
$\mathcal{D}$, for any $w$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\left(\max\left(\zeta^{-(l+1)},(\lambda\xi)^{(l+1)/l}\right)\right)^{2}l^{2}h\log(lh)|w|_{2}^{2}+\log\frac{m(l+1)}{\delta}}{\gamma^{2}m}}\right),$

where $\zeta=\min\left(\|W_{1}\|_{2},\|W_{2}\|_{2},\|W_{l}\|_{2}\right)$,
$|w|_{2}^{2}=\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}$,
$\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$,
$\lambda=\|W_{1}\|_{2}\|W_{l}\|_{2}$, and
$\xi=C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1}$.

The proof also contains three steps: 1) since MPGNNs are typically non-
homogeneous, _e.g_., when any of $\phi$, $\rho$, and $g$ is a bounded non-
linearity like Sigmoid or Tanh, we design a special quantity of learned
weights $\beta=\max(\zeta^{-1},(\lambda\xi)^{{1}/{l}})$. 2) fixing any
$\tilde{\beta}$, considering all $\beta$ that are in the range
$|\beta-\tilde{\beta}|\leq\beta/(l+1)$ and choosing $\sigma$ which depends on
$\tilde{\beta}$ so that one can apply Lemma 3.3 and 2.2 to work out the PAC-
Bayes bound; 3) taking a union bound of the previous result by considering
multiple choices of $\tilde{\beta}$ so that all possible values of $\beta$ are
covered. The case with $d\mathcal{C}=1$ is again included in the appendix. The
first step is non-trivial since we do not have the nice construction as in the
homogeneous case, _i.e_., normalizing the weights so that the spectral norms
of weights across layers are the same while the network output is unchanged.
Moreover, the quantity is vital to the whole proof framework since it
determines whether one can 1) satisfy the constraints on the random
perturbation (so that Lemma 2.2 and 3.3 are applicable) and 2) simultaneously
induce a finite covering on its range (so that the bound holds for any $w$).
Since it highly depends on the form of the perturbation bound and the network
architecture, there seems to be no general recipe on how to construct such a
quantity.

###  3.3 Comparison with Other Bounds

Statistics |  | Max Node Degree
---
$d-1$
| Max Hidden Dim
---
$h$
| Spectral Norm of
---
Learned Weights
| VC-Dimension
---
(Scarselli et al., 2018)
- | $\mathcal{O}\left(h^{4}\right)$ | -
| Rademacher
---
Complexity
(Garg et al., 2020)
$\mathcal{O}\left(d^{l-1}\sqrt{\log(d^{2l-3})}\right)$ | $\mathcal{O}\left(h\sqrt{\log h}\right)$ | $\mathcal{O}\left(\lambda\mathcal{C}\xi\sqrt{\log\left(\|W_{2}\|_{2}\lambda\xi^{2}\right)}\right)$
Ours | $\mathcal{O}\left(d^{l-1}\right)$ | $\mathcal{O}\left(\sqrt{h\log h}\right)$ | $\mathcal{O}\left(\lambda^{1+\frac{1}{l}}\xi^{1+\frac{1}{l}}\sqrt{\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}}\right)$

Table 1: Comparison of generalization bounds for GNNs. “-” means inapplicable.
$l$ is the network depth. Here
$\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$,
$\xi=C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1}$,
$\zeta=\min\left(\|W_{1}\|_{2},\|W_{2}\|_{2},\|W_{l}\|_{2}\right)$, and
$\lambda=\|W_{1}\|_{2}\|W_{l}\|_{2}$. More details about the comparison can be
found in Appendix A.5.

In this section, we compare our generalization bounds with the ones in the GNN
literature and draw connections with existing MLPs/CNNs bounds.

####  3.3.1 Comparison with Existing GNN Generalization Bounds

We compare against the VC-dimension based bound in (Scarselli et al., 2018)
and the most recent Rademacher complexity based bound in (Garg et al., 2020).
Our results are not directly comparable to (Du et al., 2019) since they
consider a “infinite-wide” class of GNNs constructed based on the neural
tangent kernel (Jacot et al., 2018), whereas we focus on commonly-used GNNs.
Comparisons to (Verma & Zhang, 2019) are also difficult since: 1) they only
show the bound for one graph convolutional layer, _i.e_., it does not depend
on the network depth $l$; and 2) their bound scales as
$\mathcal{O}\left(\lambda_{\max}^{2T}/m\right)$, where $T$ is the number of
SGD steps and $\lambda_{\max}$ is the maximum absolute eigenvalue of Laplacian
$L=D-A$. Therefore, for certain graphs555Since
$\lambda_{\max}=\max_{v\neq\bm{0}}({v^{\top}(D-A)v})/({v^{\top}v})$, we have
$\lambda_{\max}\geq(D-A)[i,i]$ by choosing $v=\bm{e}_{i}$, _i.e_., the $i$-th
standard basis. We can pick any node $i$ which has more than 1 neighbor to
make $\lambda_{\max}>1$., the generalization gap is monotonically increasing
with $T$, which cannot explain the generalization phenomenon. We compare
different bounds by examining their dependency on three terms: the maximum
node degree, the spectral norm of the learned weights, and the maximum hidden
dimension. We summarize the overall comparison in Table 1 and leave the
details such as how we convert bounds into our context to Appendix A.5.

##### Max Node Degree $(d-1)$:

The Rademacher complexity bound scales as
$\mathcal{O}\left(d^{l-1}\sqrt{\log(d^{2l-3})}\right)$ whereas ours scales as
$\mathcal{O}(d^{l-1})$666Our bound actually scales as
$\mathcal{O}\left(d^{{(l+1)(l-2)}/{l}}\right)$ which is upper bounded by
$\mathcal{O}\left(d^{l-1}\right)$.. Many real-world graphs such as social
networks tend to have large hubs (Barabási et al., 2016), which lead to very
large node degrees. Thus, our bound would be significantly better in these
scenarios.

##### Max Hidden Dimension $h$:

Our bound scales as $\mathcal{O}(\sqrt{h\log h})$ which is tighter than the
Rademacher complexity bound $\mathcal{O}\left(h\sqrt{\log h}\right)$ and the
VC-dimension bound $\mathcal{O}(h^{4})$.

##### Spectral Norm of Learned Weights:

As shown in Table 1, we cannot compare the dependencies on the spectral norm
of learned weights without knowing the actual values of the learned weights.
Therefore, we perform an empirical study in Section 4.

####  3.3.2 Connections with Existing Bounds of MLPs/CNNs

As described above, MLPs/CNNs can be viewed as special cases of GNNs. In
particular, we have two ways to show the inclusion relationship. First, we can
treat each i.i.d. sample as a node and the whole dataset as a graph without
edges. Then conventional tasks (_e.g_., classification) become node-level
tasks (_e.g_., node classification) on this graph. Second, we can treat each
i.i.d. sample as a single-node graph. Then conventional tasks (_e.g_.,
classification) becomes graph-level tasks (_e.g_., graph classification).
Since we focus on the graph classification, we adopt the second view. In
particular, MLPs/CNNs with ReLU activations are equivalent to GCNs with the
graph Laplacian $\tilde{L}=I$ (hence $d=1$). We leave the details of this
conversion to Appendix A.6. We restate the PAC-Bayes bound for MLPs/CNNs with
ReLU activations in (Neyshabur et al., 2017) as follows,

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\left({B^{2}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}({\|W_{i}\|_{F}^{2}}/{\|W_{i}\|_{2}^{2}})+\log\frac{ml}{\delta}}\right)/{\gamma^{2}m}}\right).$

Comparing it with our bound for GCNs in Theorem 3.2, it is clear that we only
add a factor $d^{l-1}$ to the first term inside the square root which is due
to the underlying graph structure of the data. If we apply GCNs to single-node
graphs, the two bounds coincide since $d=1$. Therefore, our Theorem 3.2
directly generalizes the result in (Neyshabur et al., 2017) to GCNs, which is
a strictly larger class of models than MLPs/CNNs with ReLU activations.

##  4 Experiments

(a) MPGNNs with $l=2$.

(b) MPGNNs with $l=4$.

Figure 1: Bound evaluations on real-world datasets. The maximum node degrees
(_i.e_., $d-1$) of four datasets from left to right are: $25$ (PROTEINS), $88$
(IMDB-M), $135$ (IMDB-B), and $491$ (COLLAB).

In this section, we perform an empirical comparison between our bound and the
Rademacher complexity bound for MPGNNs. We experiment on 6 synthetic datasets
of random graphs (corresponding to 6 random graph models), 3 social network
datasets (COLLAB, IMDB-BINARY, IMDB-MULTI), and a bioinformatics dataset
PROTEINS from (Yanardag & Vishwanathan, 2015). In particular, we create
synthetic datasets by generating random graphs from the Erdős–Rényi model and
the stochastic block model with different settings (_i.e_., number of blocks
and edge probabilities). All datesets focus on graph classifications. We
repeat all experiments 3 times with different random initializations and
report the means and the standard deviations. Constants are considered in the
bound computation. More details of the experimental setup, dataset statistics,
and the bound computation are provided in Appendix A.7.

As shown in Fig. 1 and Fig. 2, our bound is mostly tighter than the Rademacher
complexity bound with varying message passing steps $l$ on both synthetic and
real-world datasets. Generally, the larger the maximum node degree is, the
more our bound improves777Note that it may not be obvious from the figure as
the y axis is in log domain. Please refer to the appendix where the actual
bound values are listed in the table. over the Rademacher complexity bound
(_c.f_., PROTEINS vs. COLLAB). This could be attributed to the better
dependency on $d$ of our bound. For graphs with large node degrees (_e.g_.,
social networks like Twitter have influential users with lots of followers),
the gap could be more significant. Moreover, with the number of steps/layers
increasing, our bound also improves more in most cases. It may not be clear to
read from the figures since the y-axis is in the log domain and its range
differ from figure to figure. We also provide the numerical values of the
bound evaluations in the appendix for an exact comparison. The number of steps
is chosen to be no larger than 10 as GNNs are generally shown to perform well
with just a few steps/layers (Kipf & Welling, 2016; Jin et al., 2018). We
found $d\mathcal{C}>1$ and the geometric series
$((d\mathcal{C})^{l-1}-1)/(d\mathcal{C}-1)\gg 1$ on all datasets which imply
learned GNNs are not contraction mappings (_i.e_., $d\mathcal{C}<1$). This
also explains why both bounds becomes larger with more steps. At last, we can
see that bound values are much larger than $1$ which indicates both bounds are
still vacuous, similarly to the cases for regular neural networks in (Bartlett
et al., 2017; Neyshabur et al., 2017).

(a) MPGNNs with $l=2$.

(b) MPGNNs with $l=4$.

(c) MPGNNs with $l=6$.

(d) MPGNNs with $l=8$.

Figure 2: Bound evaluations on synthetic datasets. The maximum node degrees
(_i.e_., $d-1$) of datasets from left to right are: $25$ (ER-1), $48$ (ER-2),
$69$ (ER-3), $87$ (ER-4), $25$ (SBM-1), and $36$ (SBM-2). ‘ER-X’ and ‘SBM-X’
denote the Erdős–Rényi model and the stochastic block model with the ‘X’-th
setting respectively. Please refer to the appendix for more details.

##  5 Discussion

In this paper, we present generalization bounds for two primary classes of
GNNs, _i.e_., GCNs and MPGNNs. We show that the maximum node degree and the
spectral norms of learned weights govern the bound for both models. Our
results for GCNs generalize the bounds for MLPs/CNNs in (Neyshabur et al.,
2017), while our results for MPGNNs improve over the state-of-the-art
Rademacher complexity bound in (Garg et al., 2020). Our PAC-Bayes analysis can
be generalized to other graph problems such as node classification and link
prediction since our perturbation analysis bounds the maximum change of any
node representation. Other loss functions (_e.g_., ones for regression) could
also work in our analysis as long as they are bounded.

However, we are far from being able to explain the practical behaviors of
GNNs. Our bound values are still vacuous as shown in the experiments. Our
perturbation analysis is in the worst-case sense which may be loose for most
cases. We introduce Gaussian posterior in the PAC-Bayes framework to obtain an
analytical form of the KL divergence. Nevertheless, the actual posterior
induced by the prior and the learning process may likely to be non-Gaussian.
We also do not explicitly consider the optimization algorithm in the analysis
which clearly has an impact on the learned weights.

This work leads to a few interesting open problems for future work: (1) Is the
maximum node degree the only graph statistic that has an impact on the
generalization ability of GNNs? Investigating other graph statistics may
provide more insights on the behavior of GNNs and inspire the development of
novel models and algorithms. (2) Would the analysis still work for other
interesting GNN architectures, such as those with attention (Veličković et
al., 2017) and learnable spectral filters (Liao et al., 2019)? (3) Can recent
advancements for MLPs/CNNs, _e.g_., the compression technique in (Arora et
al., 2018) and data-dependent prior of (Parrado-Hernández et al., 2012), help
further improve the bounds for GNNs? (4) What is the impact of the
optimization algorithms like SGD on the generalization ability of GNNs? Would
graph structures play a role in the analysis of optimization?

#### Acknowledgments

We would like to thank Vikas K. Garg and Stefanie Jegelka for the valuable
feedback and discussion.

## References

  * Arora et al. (2018) Sanjeev Arora, Rong Ge, Behnam Neyshabur, and Yi Zhang.  Stronger generalization bounds for deep nets via a compression approach.  _arXiv preprint arXiv:1802.05296_ , 2018.
  * Arora et al. (2019) Sanjeev Arora, Simon S Du, Wei Hu, Zhiyuan Li, and Ruosong Wang.  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  _arXiv preprint arXiv:1901.08584_ , 2019.
  * Barabási et al. (2016) Albert-László Barabási et al.  _Network science_.  Cambridge university press, 2016.
  * Bartlett et al. (2017) Peter L Bartlett, Dylan J Foster, and Matus J Telgarsky.  Spectrally-normalized margin bounds for neural networks.  In _NIPS_ , pp. 6240–6249, 2017.
  * Battaglia et al. (2018) Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al.  Relational inductive biases, deep learning, and graph networks.  _arXiv preprint arXiv:1806.01261_ , 2018.
  * Bronstein et al. (2017) Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst.  Geometric deep learning: going beyond euclidean data.  _IEEE Signal Processing Magazine_ , 34(4):18–42, 2017.
  * Bruna et al. (2013) Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  _arXiv preprint arXiv:1312.6203_ , 2013.
  * Casas et al. (2019) Sergio Casas, Cole Gulino, Renjie Liao, and Raquel Urtasun.  Spatially-aware graph neural networks for relational behavior forecasting from sensor data.  _arXiv preprint arXiv:1910.08233_ , 2019.
  * Chen et al. (2019) Zhengdao Chen, Soledad Villar, Lei Chen, and Joan Bruna.  On the equivalence between graph isomorphism testing and function approximation with gnns.  In _NeurIPS_ , pp. 15894–15902, 2019.
  * Dai et al. (2016) Hanjun Dai, Bo Dai, and Le Song.  Discriminative embeddings of latent variable models for structured data.  In _International conference on machine learning_ , pp. 2702–2711, 2016.
  * Du et al. (2019) 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_ , pp. 5723–5733, 2019.
  * Duvenaud et al. (2015) David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In _NIPS_ , pp. 2224–2232, 2015.
  * Dziugaite & Roy (2017) Gintare Karolina Dziugaite and Daniel M Roy.  Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data.  _arXiv preprint arXiv:1703.11008_ , 2017.
  * Garg et al. (2020) Vikas K Garg, Stefanie Jegelka, and Tommi Jaakkola.  Generalization and representational limits of graph neural networks.  _arXiv preprint arXiv:2002.06157_ , 2020.
  * Gilmer et al. (2017) Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  _arXiv preprint arXiv:1704.01212_ , 2017.
  * Gori et al. (2005) Marco Gori, Gabriele Monfardini, and Franco Scarselli.  A new model for learning in graph domains.  In _Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005._ , volume 2, pp. 729–734. IEEE, 2005.
  * Hamilton et al. (2017) Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In _NIPS_ , pp. 1024–1034, 2017.
  * Harvey et al. (2017) Nick Harvey, Christopher Liaw, and Abbas Mehrabian.  Nearly-tight vc-dimension bounds for piecewise linear neural networks.  In _Conference on Learning Theory_ , pp. 1064–1068, 2017.
  * Henrion et al. (2017) Isaac Henrion, Johann Brehmer, Joan Bruna, Kyunghyun Cho, Kyle Cranmer, Gilles Louppe, and Gaspar Rochette.  Neural message passing for jet physics.  _NIPS Workshop on Deep Learning for Physical Sciences_ , 2017.
  * Jacot et al. (2018) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In _NeurIPS_ , pp. 8571–8580, 2018.
  * Jin et al. (2018) Wengong Jin, Kevin Yang, Regina Barzilay, and Tommi Jaakkola.  Learning multimodal graph-to-graph translation for molecular optimization.  _arXiv preprint arXiv:1812.01070_ , 2018.
  * Kipf & Welling (2016) Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Langford & Shawe-Taylor (2003) John Langford and John Shawe-Taylor.  Pac-bayes & margins.  In _NIPS_ , pp. 439–446, 2003.
  * Li et al. (2015) Yujia Li, Daniel Tarlow, Marc Brockschmidt, and Richard Zemel.  Gated graph sequence neural networks.  _arXiv preprint arXiv:1511.05493_ , 2015.
  * Liao et al. (2019) Renjie Liao, Zhizhen Zhao, Raquel Urtasun, and Richard Zemel.  Lanczosnet: Multi-scale deep graph convolutional networks.  In _ICLR_ , 2019.
  * Maron et al. (2018) Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman.  Invariant and equivariant graph networks.  _arXiv preprint arXiv:1812.09902_ , 2018.
  * McAllester (2003) David McAllester.  Simplified pac-bayesian margin bounds.  In _Learning theory and Kernel machines_ , pp. 203–215. Springer, 2003.
  * McAllester (1999) David A McAllester.  Pac-bayesian model averaging.  In _Proceedings of the twelfth annual conference on Computational learning theory_ , pp. 164–170, 1999.
  * Monti et al. (2017) Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda, and Michael M Bronstein.  Geometric deep learning on graphs and manifolds using mixture model cnns.  In _CVPR_ , pp. 5115–5124, 2017.
  * Neyshabur et al. (2017) Behnam Neyshabur, Srinadh Bhojanapalli, and Nathan Srebro.  A pac-bayesian approach to spectrally-normalized margin bounds for neural networks.  _arXiv preprint arXiv:1707.09564_ , 2017.
  * Parrado-Hernández et al. (2012) Emilio Parrado-Hernández, Amiran Ambroladze, John Shawe-Taylor, and Shiliang Sun.  Pac-bayes bounds with data dependent priors.  _The Journal of Machine Learning Research_ , 13(1):3507–3531, 2012.
  * Sanchez-Gonzalez et al. (2018) Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller, Raia Hadsell, and Peter Battaglia.  Graph networks as learnable physics engines for inference and control.  _arXiv preprint arXiv:1806.01242_ , 2018.
  * Scarselli et al. (2008) Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  _IEEE Transactions on Neural Networks_ , 20(1):61–80, 2008.
  * Scarselli et al. (2018) Franco Scarselli, Ah Chung Tsoi, and Markus Hagenbuchner.  The vapnik–chervonenkis dimension of graph and recursive neural networks.  _Neural Networks_ , 108:248–259, 2018.
  * Tropp (2012) Joel A Tropp.  User-friendly tail bounds for sums of random matrices.  _Foundations of computational mathematics_ , 12(4):389–434, 2012.
  * Valiant (1984) Leslie G Valiant.  A theory of the learnable.  _Communications of the ACM_ , 27(11):1134–1142, 1984.
  * Veličković et al. (2017) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ , 2017.
  * Verma & Zhang (2019) Saurabh Verma and Zhi-Li Zhang.  Stability and generalization of graph convolutional neural networks.  In _KDD_ , pp. 1539–1548, 2019.
  * Xu et al. (2018) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  _arXiv preprint arXiv:1810.00826_ , 2018.
  * Yanardag & Vishwanathan (2015) Pinar Yanardag and SVN Vishwanathan.  Deep graph kernels.  In _KDD_ , pp. 1365–1374, 2015.

##  Appendix A Appendix

We summarize the notations used throughout the paper in Table 2. In the
following sections, we provide proofs of all results in the main text and the
additional details.

Symbol | Meaning
---|---
$\mathbb{R}$ | the set of real numbers
$\mathbb{R}^{m\times n}$ | the set of $m\times n$ real matrices
$\mathbb{N}^{+}_{k}$ | the set of first $k$ positive numbers
$|\cdot|_{p}$ | the vector p-norm
$\|\cdot\|_{p}$ | the operator norm induced by the vector p-norm
$\|\cdot\|_{F}$ | the Frobenius norm
$X[i,j]$ | the ($i,j$)-th element of matrix $X$
$X[i,:]$ | the $i$-th row of matrix $X$
$X[:,i]$ | the $i$-th column of matrix $X$
$\bm{1}_{n}$ | a all-one vector with size $n$
$A$ | an adjacency matrix
$\tilde{A}$ | an adjacency matrix plus the identity matrix
$B$ | the radius of the $\ell_{2}$-ball where an input node feature lies
$C_{\text{in}}$ | an incidence matrix of incoming nodes
$C_{\text{out}}$ | an incidence matrix of outgoing nodes
$\phi$, $\rho$, $g$ | non-linearities in MPGNN
$C_{\phi}$, $C_{\rho}$, $C_{g}$ | Lipschitz constants of $\phi$, $\rho$, $g$ under the vector 2-norm
$\mathcal{C}$ | the percolation complexity
$D$ | the degree matrix
$\mathcal{D}$ | the unknown data distribution
$d$ | the maximum node degree plus one
$e$ | the Euler’s number
$f_{w}$ | a model parameterized by vector $w$
$\mathcal{G}$ | the space of graph
$h$ | the maximum hidden dimension
$H$ | a node representation matrix
$\mathcal{H}$ | the hypothesis/model class
$I$ | the identity matrix
$l$ | the number of graph convolution layers / message passing steps
$L$ | the loss function
$\tilde{L}$ | the Laplacian matrix
$m$ | the number of training samples
$\mathbb{P}$, $\mathbb{E}$ | probability and expectation of a random variable
$P$ | the prior distribution over hypothesis class
$Q$ | the posterior distribution over hypothesis class
$S$ | a set of training samples
$W$ | a weight matrix
$X$ | a node feature matrix where each row corresponds to a node
$\mathcal{X}$ | the space of node feature
$y$ | the graph class label
$\gamma$ | the margin parameter
$z$ | a data triplet $(A,X,y)$
$\mathcal{Z}$ | the space of data triplet
$\log$ | the natural logarithm
Table 2: Summary of important notations.

###  A.1 PAC Bayes Results

For completeness, we provide the proofs of the standard PAC-Bayes results as
below.

######  Lemma A.1.

For non-negative continuous random variables $X$, we have

$\displaystyle\mathbb{E}[X]=\int_{0}^{\infty}\mathbb{P}(X\geq\nu)\mathrm{d}\nu.$

###### Proof.

$\displaystyle\mathbb{E}[X]$
$\displaystyle=\int_{0}^{\infty}X\mathbb{P}(X)\mathrm{d}X$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{X}\bm{1}\mathrm{d}\nu\mathbb{P}(X)\mathrm{d}X$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{X}\mathbb{P}(X)\mathrm{d}\nu\mathrm{d}X$
$\displaystyle=\int_{0}^{\infty}\int_{\nu}^{\infty}\mathbb{P}(X)\mathrm{d}X\mathrm{d}\nu\qquad(\text{region
of the integral is the same})$
$\displaystyle=\int_{0}^{\infty}\mathbb{P}(X\geq\nu)\mathrm{d}\nu$

∎

######  Lemma A.2.

[2-side] Let $X$ be a random variable satisfying
$\mathbb{P}(X\geq\epsilon)\leq e^{-2m\epsilon^{2}}$ and
$\mathbb{P}(X\leq-\epsilon)\leq e^{-2m\epsilon^{2}}$ where $m\geq 1$ and
$\epsilon>0$, we have

$\displaystyle\mathbb{E}[e^{2(m-1)X^{2}}]\leq 2m.$

###### Proof.

If $m=1$, the inequality holds trivially. Let us now consider $m>1$.

$\displaystyle\mathbb{E}[e^{2(m-1)X^{2}}]$
$\displaystyle=\int_{0}^{\infty}\mathbb{P}\left(e^{2(m-1)X^{2}}\geq\nu\right)\mathrm{d}\nu\qquad(\text{Lemma
\ref{lemma:expect_nonneg_RV}})$
$\displaystyle=\int_{0}^{\infty}\mathbb{P}\left(X^{2}\geq\frac{\log\nu}{2(m-1)}\right)\mathrm{d}\nu\qquad$
$\displaystyle=\int_{0}^{\infty}\mathbb{P}\left(X\geq\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu+\int_{0}^{\infty}\mathbb{P}\left(X\leq-\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu$
(5)

$\displaystyle\int_{0}^{\infty}\mathbb{P}\left(X\geq\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu$
$\displaystyle=\int_{0}^{1}\mathbb{P}\left(X\geq\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu+\int_{1}^{\infty}\mathbb{P}\left(X\geq\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu$
$\displaystyle\leq
1+\int_{1}^{\infty}\mathbb{P}\left(X\geq\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu$
$\displaystyle\leq
1+\int_{1}^{\infty}e^{-2m\frac{\log\nu}{2(m-1)}}\mathrm{d}\nu$
$\displaystyle=1+\left(\left.-(m-1)\nu^{-\frac{1}{m-1}}\right\rvert_{1}^{\infty}\right)$
$\displaystyle=m$ (6)

Similarly, we can show that

$\displaystyle\int_{0}^{\infty}\mathbb{P}\left(X\leq-\sqrt{\frac{\log\nu}{2(m-1)}}\right)\mathrm{d}\nu\leq
m$ (7)

Combining Eq. (A.1) and Eq. (A.1), we finish the proof. ∎

######  Theorem 2.1.

(Two-side) Let $P$ be a prior distribution over $\mathcal{H}$ and let
$\delta\in(0,1)$. Then, with probability $1-\delta$ over the choice of an
i.i.d. training set $S$ according to $\mathcal{D}$, for all distributions $Q$
over $\mathcal{H}$ and any $\gamma>0$, we have

$\displaystyle L_{\mathcal{D},\gamma}(Q)\leq
L_{S,\gamma}(Q)+\sqrt{\frac{D_{\mathrm{KL}}(Q\|P)+\log\frac{2m}{\delta}}{2(m-1)}}$

###### Proof.

Let $\Delta(h)=L_{\mathcal{D},\gamma}(h)-L_{S,\gamma}(h)$. For any function
$f(h)$, we have

$\displaystyle\mathbb{E}_{h\sim Q}[f(h)]$ $\displaystyle=\mathbb{E}_{h\sim
Q}[\log e^{f(h)}]$ $\displaystyle=\mathbb{E}_{h\sim Q}[\log
e^{f(h)}+\log\frac{Q}{P}+\log\frac{P}{Q}]$
$\displaystyle=D_{\mathrm{KL}}(Q\|P)+\mathbb{E}_{h\sim
Q}\left[\log\left(\frac{P}{Q}e^{f(h)}\right)\right]$ $\displaystyle\leq
D_{\mathrm{KL}}(Q\|P)+\log\mathbb{E}_{h\sim
Q}\left[\frac{P}{Q}e^{f(h)}\right]\qquad(\text{Jensen's inequality})$
$\displaystyle=D_{\mathrm{KL}}(Q\|P)+\log\mathbb{E}_{h\sim
P}\left[e^{f(h)}\right].$ (8)

Let $f(h)=2(m-1)\Delta(h)^{2}$. We have

$\displaystyle 2(m-1)\mathbb{E}_{h\sim Q}[\Delta(h)]^{2}$ $\displaystyle\leq
2(m-1)\mathbb{E}_{h\sim Q}[\Delta(h)^{2}]\qquad(\text{Jensen's inequality})$
$\displaystyle\leq D_{\mathrm{KL}}(Q\|P)+\log\mathbb{E}_{h\sim
P}\left[e^{2(m-1)\Delta(h)^{2}}\right].$ (9)

Since $L_{\mathcal{D}}(h)\in[0,1]$, based on Hoeffding’s inequality, for any
$\epsilon>0$, we have

$\displaystyle\mathbb{P}(\Delta(h)\geq\epsilon)\leq e^{-2m\epsilon^{2}}$
$\displaystyle\mathbb{P}(\Delta(h)\leq-\epsilon)\leq e^{-2m\epsilon^{2}}$

Hence, based on Lemma A.2, we have

$\displaystyle\mathbb{E}_{S}\left[e^{2(m-1)\Delta(h)^{2}}\right]\leq 2m~{}$
$\displaystyle\Rightarrow~{}\mathbb{E}_{h\sim
P}\left[\mathbb{E}_{S}\left[e^{2(m-1)\Delta(h)^{2}}\right]\right]\leq 2m$
$\displaystyle\Leftrightarrow~{}\mathbb{E}_{S}\left[\mathbb{E}_{h\sim
P}\left[e^{2(m-1)\Delta(h)^{2}}\right]\right]\leq 2m$

Based on Markov’s inequality, we have

$\displaystyle\mathbb{P}\left(\mathbb{E}_{h\sim
P}\left[e^{2(m-1)\Delta(h)^{2}}\right]\geq\frac{2m}{\delta}\right)\leq\frac{\delta\mathbb{E}_{S}\left[\mathbb{E}_{h\sim
P}\left[e^{2(m-1)\Delta(h)^{2}}\right]\right]}{2m}\leq\delta.$ (10)

Combining Eq. (A.1) and Eq. (10), with probability $1-\delta$, we have

$\displaystyle\mathbb{E}_{h\sim
Q}[\Delta(h)]^{2}\leq\frac{D_{\mathrm{KL}}(Q\|P)+\log\left(\frac{2m}{\delta}\right)}{2(m-1)}$
(11)

which proves the theorem. ∎

######  Lemma 2.2.

Let $f_{w}(x):\mathcal{X}\rightarrow\mathbb{R}^{k}$ be any model with
parameters $w$, and $P$ be any distribution on the parameters that is
independent of the training data. For any $w$, we construct a posterior
$Q(w+u)$ by adding any random perturbation $u$ to $w$, _s.t_.,
$\mathbb{P}(\max_{x\in\mathcal{X}}|f_{w+u}(x)-f_{w}(x)|_{\infty}<\frac{\gamma}{4})>\frac{1}{2}$.
Then, for any $\gamma,\delta>0$, with probability at least $1-\delta$ over an
i.i.d. size-$m$ training set $S$ according to $\mathcal{D}$, for any $w$, we
have:

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\sqrt{\frac{2D_{\mathrm{KL}}(Q(w+u)\|P)+\log\frac{8m}{\delta}}{2(m-1)}}$

###### Proof.

Let $\tilde{w}=w+u$. Let $\mathcal{C}$ be the set of perturbation with the
following property,

$\displaystyle\mathcal{C}=\left\\{w^{\prime}\middle|\max_{x\in\mathcal{X}}|f_{w^{\prime}}(x)-f_{w}(x)|_{\infty}<\frac{\gamma}{4}\right\\}.$
(12)

$\tilde{w}=w+u$ ($w$ is deterministic and $u$ is stochastic) is distributed
according to $Q(\tilde{w})$. We now construct a new posterior $\tilde{Q}$ as
follows,

$\displaystyle\tilde{Q}(\tilde{w})=\begin{cases}\frac{1}{Z}Q(\tilde{w})&\tilde{w}\in\mathcal{C}\\\
0&\tilde{w}\in\bar{\mathcal{C}}.\end{cases}$ (13)

Here
$Z=\int_{\tilde{w}\in\mathcal{C}}\mathrm{d}Q(\tilde{w})=\underset{\tilde{w}\sim
Q}{\mathbb{P}}(\tilde{w}\in\mathcal{C})$ and $\bar{\mathcal{C}}$ is the
complement set of $\mathcal{C}$. We know from the assumption that
$Z>\frac{1}{2}$. Therefore, for any $\tilde{w}\sim\tilde{Q}$, we have

$\displaystyle\left.\max_{i\in\mathbb{N}_{k}^{+},j\in\mathbb{N}_{k}^{+},x\in\mathcal{X}}\middle|\left|f_{\tilde{w}}(x)[i]-f_{\tilde{w}}(x)[j]\right|-|f_{w}(x)[i]-f_{w}(x)[j]|\right|$
$\displaystyle\leq$
$\displaystyle\left.\max_{i\in\mathbb{N}_{k}^{+},j\in\mathbb{N}_{k}^{+},x\in\mathcal{X}}\middle|f_{\tilde{w}}(x)[i]-f_{\tilde{w}}(x)[j]-f_{w}(x)[i]+f_{w}(x)[j]\right|$
$\displaystyle\leq$
$\displaystyle\left.\max_{i\in\mathbb{N}_{k}^{+},j\in\mathbb{N}_{k}^{+},x\in\mathcal{X}}\middle|f_{\tilde{w}}(x)[i]-f_{w}(x)[i]\middle|+\middle|f_{\tilde{w}}(x)[j]-f_{w}(x)[j]\right|$
$\displaystyle\leq$
$\displaystyle\left.\max_{i\in\mathbb{N}_{k}^{+},x\in\mathcal{X}}\middle|f_{\tilde{w}}(x)[i]-f_{w}(x)[i]\middle|+\max_{j\in\mathbb{N}_{k}^{+},x\in\mathcal{X}}\middle|f_{\tilde{w}}(x)[j]-f_{w}(x)[j]\right|$
$\displaystyle<\frac{\gamma}{4}+\frac{\gamma}{4}=\frac{\gamma}{2}$ (14)

Recall that

$\displaystyle L_{\mathcal{D}}(f_{w},0)$
$\displaystyle=\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{w}(x)[y]\leq\max_{j\neq
y}f_{w}(x)[j]\right)$ $\displaystyle
L_{\mathcal{D}}(f_{\tilde{w}},\frac{\gamma}{2})$
$\displaystyle=\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{\tilde{w}}(x)[y]\leq\frac{\gamma}{2}+\max_{j\neq
y}f_{\tilde{w}}(x)[j]\right).$

Denoting $j_{1}^{\ast}=\operatorname*{arg\,max}_{j\neq y}f_{\tilde{w}}(x)[j]$
and $j_{2}^{\ast}=\operatorname*{arg\,max}_{j\neq y}f_{w}(x)[j]$, from Eq.
(A.1), we have

$\displaystyle\left|f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{2}^{\ast}]-f_{w}(x)[y]+f_{w}(x)[j_{2}^{\ast}]\middle|<\frac{\gamma}{2}\right.$
$\displaystyle\Rightarrow~{}$ $\displaystyle
f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{2}^{\ast}]<f_{w}(x)[y]-f_{w}(x)[j_{2}^{\ast}]+\frac{\gamma}{2}$
(15)

Note that since $f_{\tilde{w}}(x)[j_{1}^{\ast}]\geq
f_{\tilde{w}}(x)[j_{2}^{\ast}]$, we have

$\displaystyle f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{1}^{\ast}]$
$\displaystyle\leq f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{2}^{\ast}]$
$\displaystyle\leq
f_{w}(x)[y]-f_{w}(x)[j_{2}^{\ast}]+\frac{\gamma}{2}\qquad(\text{Eq.
(\ref{eq:pac_bayes_deterministic_tmp_2})})$

Therefore, we have

$\displaystyle f_{w}(x)[y]-f_{w}(x)[j_{2}^{\ast}]\leq
0~{}\Rightarrow~{}f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{1}^{\ast}]\leq\frac{\gamma}{2},$

which indicates $\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{w}(x)[y]\leq
f_{w}(x)[j_{2}^{\ast}]\right)\leq\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{\tilde{w}}(x)[y]\leq
f_{\tilde{w}}(x)[j_{1}^{\ast}]+\frac{\gamma}{2}\right)$, or equivalently

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{\mathcal{D},\frac{\gamma}{2}}(f_{\tilde{w}}).$ (16)

Note that this holds for any perturbation $\tilde{w}\sim\tilde{Q}$.

Again, recall that

$\displaystyle L_{\mathcal{D},\frac{\gamma}{2}}(f_{\tilde{w}})$
$\displaystyle=\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{\tilde{w}}(x)[y]\leq\frac{\gamma}{2}+\max_{j\neq
y}f_{\tilde{w}}(x)[j]\right)$ $\displaystyle L_{\mathcal{D},\gamma}(f_{w})$
$\displaystyle=\underset{z\sim\mathcal{D}}{\mathbb{P}}\left(f_{w}(x)[y]\leq\gamma+\max_{j\neq
y}f_{w}(x)[j]\right)$

From Eq. (A.1), we have

$\displaystyle\left|f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{1}^{\ast}]-f_{w}(x)[y]+f_{w}(x)[j_{1}^{\ast}]\middle|<\frac{\gamma}{2}\right.$
$\displaystyle\Rightarrow~{}$ $\displaystyle
f_{w}(x)[y]-f_{w}(x)[j_{1}^{\ast}]<f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{1}^{\ast}]+\frac{\gamma}{2}$
(17)

Note that since $f_{{w}}(x)[j_{2}^{\ast}]\geq f_{{w}}(x)[j_{1}^{\ast}]$, we
have

$\displaystyle f_{{w}}(x)[y]-f_{{w}}(x)[j_{2}^{\ast}]$ $\displaystyle\leq
f_{{w}}(x)[y]-f_{{w}}(x)[j_{1}^{\ast}]$ $\displaystyle\leq
f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{1}^{\ast}]+\frac{\gamma}{2}\qquad(\text{Eq.
(\ref{eq:pac_bayes_deterministic_tmp_4})})$

Therefore, we have

$\displaystyle
f_{\tilde{w}}(x)[y]-f_{\tilde{w}}(x)[j_{1}^{\ast}]\leq\frac{\gamma}{2}~{}\Rightarrow~{}f_{{w}}(x)[y]-f_{{w}}(x)[j_{2}^{\ast}]\leq\gamma,$

which indicates $L_{\mathcal{D},\frac{\gamma}{2}}(f_{\tilde{w}})\leq
L_{\mathcal{D},\gamma}(f_{{w}})$. Therefore, from the perspective of the
empirical estimation of the probability, for any $\tilde{w}\sim\tilde{Q}$, we
almost surely have

$\displaystyle L_{S,\frac{\gamma}{2}}(f_{\tilde{w}})\leq
L_{S,\gamma}(f_{{w}}).$ (18)

Now with probability at least $1-\delta$, we have

$\displaystyle L_{\mathcal{D},0}(f_{w})$
$\displaystyle\leq\mathbb{E}_{\tilde{w}\sim\tilde{Q}}\left[L_{\mathcal{D},\frac{\gamma}{2}}(f_{\tilde{w}})\right]\qquad(\text{Eq.
(\ref{eq:pac_bayes_deterministic_tmp_3})})$
$\displaystyle\leq\mathbb{E}_{\tilde{w}\sim\tilde{Q}}\left[L_{S,\frac{\gamma}{2}}(f_{\tilde{w}})\right]+\sqrt{\frac{D_{\mathrm{KL}}(\tilde{Q}\|P)+\log\frac{2m}{\delta}}{2(m-1)}}\qquad(\text{Theorem
\ref{thm:pac_bayes}})$ $\displaystyle\leq
L_{S,\gamma}(f_{w})+\sqrt{\frac{D_{\mathrm{KL}}(\tilde{Q}\|P)+\log\frac{2m}{\delta}}{2(m-1)}}\qquad(\text{Eq.
(\ref{eq:pac_bayes_deterministic_tmp_5})})$ (19)

Note that

$\displaystyle D_{\mathrm{KL}}(Q\|P)$
$\displaystyle=\underset{{\tilde{w}\in\mathcal{C}}}{\int}Q\log\frac{Q}{P}\mathrm{d}\tilde{w}+\underset{\tilde{w}\in\bar{\mathcal{C}}}{\int}Q\log\frac{Q}{P}\mathrm{d}\tilde{w}$
$\displaystyle=\underset{\tilde{w}\in\mathcal{C}}{\int}\frac{QZ}{Z}\log\frac{Q}{ZP}\mathrm{d}\tilde{w}+\underset{\tilde{w}\in\mathcal{C}}{\int}Q\log
Z\mathrm{d}\tilde{w}$
$\displaystyle\quad+\underset{\tilde{w}\in\bar{\mathcal{C}}}{\int}\frac{Q(1-Z)}{1-Z}\log\frac{Q}{(1-Z)P}\mathrm{d}\tilde{w}+\underset{\tilde{w}\in\bar{\mathcal{C}}}{\int}Q\log(1-Z)\mathrm{d}\tilde{w}$
$\displaystyle=ZD_{\mathrm{KL}}(\tilde{Q}\|P)+(1-Z)D_{\mathrm{KL}}(\bar{Q}\|P)-H(Z),$
(20)

where $\bar{Q}$ denotes the normalized density of $Q$ restricted to
$\bar{\mathcal{C}}$. $H(Z)$ is the entropy of a Bernoulli random variable with
parameter $Z$. Since we know $\frac{1}{2}\leq Z\leq 1$ from the beginning,
$0\leq H(Z)\leq\log 2$, and $D_{\mathrm{KL}}$ is nonnegative, from Eq. (A.1),
we have

$\displaystyle D_{\mathrm{KL}}(\tilde{Q}\|P)$
$\displaystyle=\frac{1}{Z}\left[D_{\mathrm{KL}}(Q\|P)+H(Z)-(1-Z)D_{\mathrm{KL}}(\bar{Q}\|P)\right]$
$\displaystyle\leq\frac{1}{Z}\left[D_{\mathrm{KL}}(Q\|P)+H(Z)\right]$
$\displaystyle\leq 2D_{\mathrm{KL}}(Q\|P)+2\log 2.$ (21)

Combining Eq. (A.1) and Eq. (A.1), we have

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\sqrt{\frac{D_{\mathrm{KL}}(Q\|P)+\frac{1}{2}\log\frac{8m}{\delta}}{m-1}},$
(22)

which finishes the proof. ∎

###  A.2 Graph Results

In this part, we provide a result on the graph Laplacian used by GCNs along
with the proof. It is used in the perturbation analysis of GCNs.

######  Lemma A.3.

Let $A$ be the binary adjacency matrix of an arbitrary simple graph $G=(V,E)$
and $\tilde{A}=A+I$. We define the graph Laplacian
$L=D^{-\frac{1}{2}}\tilde{A}D^{-\frac{1}{2}}$ where $D$ is the degree matrix
of $\tilde{A}$. Then we have $\|L\|_{1}=\|L\|_{\infty}\leq\sqrt{d}$,
$\|L\|_{2}\leq 1$, and $\|L\|_{F}\leq\sqrt{r}$ where $r$ is the rank of $L$
and $d-1$ is the maximum node degree of $G$.

###### Proof.

First, $\tilde{A}$ is symmetric and element-wise nonnegative. Denoting
$n=|V|$, we have $\tilde{A}\in\mathbb{R}^{n\times n}$,
$D_{i}=\sum_{j=1}^{n}\tilde{A}[i,j]$, and $1\leq D_{i}\leq d,\forall
i\in\mathbb{N}^{+}_{n}$. It is easy to show that
$L[i,j]={\tilde{A}[i,j]}/{\sqrt{D_{i}D_{j}}}$.

For the infinity norm and 1-norm, we have
$\|L\|_{1}=\|L^{\top}\|_{\infty}=\|L\|_{\infty}$. Moreover,

$\displaystyle\|L\|_{\infty}$
$\displaystyle=\max_{i\in\mathbb{N}^{+}_{n}}\sum_{j=1}^{n}\left|L[i,j]\right|$
$\displaystyle=\max_{i\in\mathbb{N}^{+}_{n}}\sum_{j=1}^{n}\frac{\tilde{A}[i,j]}{\sqrt{D_{i}D_{j}}}$
$\displaystyle\leq\max_{i\in\mathbb{N}^{+}_{n}}\frac{1}{\sqrt{D_{i}}}\sum_{j=1}^{n}\tilde{A}[i,j]$
$\displaystyle=\max_{i\in\mathbb{N}^{+}_{n}}\sqrt{D_{i}}$
$\displaystyle\leq\sqrt{d}$ (23)

For the spectral norm, based on the definition, we have

$\displaystyle\|L\|_{2}=\sup_{x\neq 0}\frac{|Lx|_{2}}{|x|_{2}}=\sigma_{\max},$
(24)

where $\sigma_{\max}$ is the maximum singular value of $L$. Since $L$ is
symmetric, we have $\sigma_{i}=|\lambda_{i}|$ where $\lambda_{i}$ is the
$i$-th eigenvalue of $L$. Hence, $\sigma_{\max}=\max_{i}|\lambda_{i}|$. From
Raylaigh quotient and Courant–Fischer minimax theorem, we have

$\displaystyle\|L\|_{2}$ $\displaystyle=\max_{i}|\lambda_{i}|=\max_{x\neq
0}\left|\frac{x^{\top}Lx}{x^{\top}x}\right|$ $\displaystyle=\max_{x\neq
0}\left|\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}L[i,j]x_{i}x_{j}}{\sum_{i=1}^{n}x_{i}^{2}}\right|$
$\displaystyle=\max_{x\neq
0}\left|\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}{\tilde{A}[i,j]x_{i}x_{j}}/{\sqrt{D_{i}D_{j}}}}{\sum_{i=1}^{n}x_{i}^{2}}\right|$
$\displaystyle=\max_{x\neq
0}\left|\frac{\sum_{(i,j)\in\tilde{E}}{x_{i}x_{j}}/{\sqrt{D_{i}D_{j}}}}{\sum_{i=1}^{n}x_{i}^{2}}\right|$
$\displaystyle\leq\max_{x\neq
0}\left|\frac{\frac{1}{2}\sum_{(i,j)\in\tilde{E}}\left({x_{i}^{2}}/{D_{i}}+{x_{j}^{2}}/{D_{j}}\right)}{\sum_{i=1}^{n}x_{i}^{2}}\right|$
$\displaystyle=\max_{x\neq
0}\left|\frac{\sum_{(i,j)\in\tilde{E}}{x_{i}^{2}}/{D_{i}}}{\sum_{i=1}^{n}x_{i}^{2}}\right|=\max_{x\neq
0}\left|\frac{\sum_{i=1}^{n}{x_{i}^{2}}}{\sum_{i=1}^{n}x_{i}^{2}}\right|=1,$
(25)

where $\tilde{E}$ is the union of the set of edges $E$ in the original graph
and the set of self-loops. For Frobenius norm, we have
$\|L\|_{F}\leq\sqrt{r}\|L\|_{2}\leq\sqrt{r}$ where $r$ is the rank of $L$. ∎

###  A.3 GCN Results

In this part, we provide the proofs of the main results regarding GCNs.

######  Lemma 3.1.

(GCN Perturbation Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-layer GCN. Then for any $w$, and $x\in\mathcal{X}_{B,h_{0}}$, and any
perturbation $u=\text{vec}(\\{U_{i}\\}_{i=1}^{l})$ such that $\forall
i\in\mathbb{N}^{+}_{l}$, $\|U_{i}\|_{2}\leq\frac{1}{l}\|W_{i}\|_{2}$, the
change in the output of GCN is bounded as,

$\displaystyle\left|f_{w+u}(X,A)-f_{w}(X,A)\right|_{2}\leq
eBd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\sum_{k=1}^{l}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}$

###### Proof.

We first perform the recursive perturbation analysis on node representations
of all layers except the last one, _i.e_., the readout layer. Then we derive
the bound for the graph representation of the last readout layer.

##### Perturbation Analysis on Node Representations.

In GCN, for any layer $j<l$ besides the last readout one, the node
representations are,

$\displaystyle
f_{w}^{j}(X,A)=H_{j}=\sigma_{j}\left(\tilde{L}H_{j-1}W_{j}\right).$ (26)

We add perturbation $u$ to the weights $w$, _i.e_., for the $j$-th layer, the
perturbed weights are $W_{j}+U_{j}$. For the ease of notation, we use the
superscript of prime to denote the perturbed node representations, _e.g_.,
$H_{j}^{\prime}=f_{w+u}^{j}(X,A)$. Let
$\Delta_{j}=f_{w+u}^{j}(X,A)-f_{w}^{j}(X,A)=H_{j}^{\prime}-H_{j}$. Note that
$\Delta_{j}\in\mathbb{R}^{n\times h_{j}}$. Let
$\Psi_{j}=\max\limits_{i}\left|\Delta_{j}[i,:]\right|_{2}=\max\limits_{i}\left|H_{j}^{\prime}[i,:]-H_{j}[i,:]\right|_{2}$
and $\Phi_{j}=\max\limits_{i}\left|H_{j}[i,:]\right|_{2}$. We denote the
$u_{j}^{\ast}=\operatorname*{arg\,max}\limits_{i}\left|\Delta_{j}[i,:]\right|_{2}$
and
$v_{j}^{\ast}=\operatorname*{arg\,max}\limits_{i}\left|H_{j}[i,:]\right|_{2}$.

##### Upper Bound on the Max Node Representation

For any layer $j<l$, we can derive an upper bound on the maximum (w.r.t.
$\ell_{2}$ norm) node representation as follows,

$\displaystyle\Phi_{j}$
$\displaystyle=\max_{i}\left|H_{j}[i,:]\right|_{2}=\left|\left(\sigma_{j}\left(\tilde{L}H_{j-1}W_{j}\right)\right)[v_{j}^{\ast},:]\right|_{2}=\left|\sigma_{j}\left(\left(\tilde{L}H_{j-1}W_{j}\right)[v_{j}^{\ast},:]\right)\right|_{2}$
$\displaystyle\leq\left|\left(\tilde{L}H_{j-1}W_{j}\right)[v_{j}^{\ast},:]\right|_{2}\qquad(\text{Lipschitz
property of ReLU under vector 2-norm})$
$\displaystyle=\left|\left(\tilde{L}H_{j-1}\right)[v_{j}^{\ast},:]W_{j}\right|_{2}$
$\displaystyle\leq\left|\left(\tilde{L}H_{j-1}\right)[v_{j}^{\ast},:]\right|_{2}\left\|W_{j}\right\|_{2}=\left|\sum\nolimits_{k\in\mathcal{N}_{v_{j}^{\ast}}}\tilde{L}[v_{j}^{\ast},k]H_{j-1}[k,:]\right|_{2}\left\|W_{j}\right\|_{2}$
$\displaystyle\leq\sum\nolimits_{k\in\mathcal{N}_{v_{j}^{\ast}}}\tilde{L}[v_{j}^{\ast},k]\left|H_{j-1}[k,:]\right|_{2}\left\|W_{j}\right\|_{2}$
$\displaystyle\leq\sum\nolimits_{k\in\mathcal{N}_{v_{j}^{\ast}}}\tilde{L}[v_{j}^{\ast},k]\Phi_{j-1}\left\|W_{j}\right\|_{2}\qquad\left(\text{since
}\forall i,\left|H_{j-1}[i,:]\right|_{2}\leq\Phi_{j-1}\right)$
$\displaystyle\leq d^{\frac{1}{2}}\Phi_{j-1}\left\|W_{j}\right\|_{2}$
$\displaystyle\leq
d^{\frac{j}{2}}\Phi_{0}\prod_{i=1}^{j}\|W_{i}\|_{2}\qquad(\text{unroll the
recursion})$ $\displaystyle\leq d^{\frac{j}{2}}B\prod_{i=1}^{j}\|W_{i}\|_{2},$
(27)

where in the last inequality we use the fact
$\Phi_{0}=\max_{i}\left|X[i,:]\right|_{2}\leq B$ based on the assumption A3.
$\mathcal{N}_{v_{j}^{\ast}}$ is the set of neighboring nodes (including
itself) of node $v_{j}^{\ast}$. In the third from the last inequality, we use
the Lemma A.3 to derive the following fact that $\forall i$,

$\displaystyle\sum_{k\in\mathcal{N}_{i}}\tilde{L}[i,k]=\sum_{k\in\mathcal{N}_{i}}\left|\tilde{L}[i,k]\right|\leq\left\|\tilde{L}\right\|_{\infty}\leq\sqrt{d}.$
(28)

##### Upper Bound on the Max Change of Node Representation.

For any layer $j<l$, we can derive an upper bound on the maximum (w.r.t.
$\ell_{2}$ norm) change between the representations with and without the
weight perturbation for any node as follows,

$\displaystyle\Psi_{j}$
$\displaystyle=\max\limits_{i}\left|H_{j}^{\prime}[i,:]-H_{j}[i,:]\right|_{2}=\left|\sigma_{j}\left(\tilde{L}H_{j-1}^{\prime}(W_{j}+U_{j})\right)[u_{j}^{\ast},:]-\sigma_{j}\left(\tilde{L}H_{j-1}W_{j}\right)[u_{j}^{\ast},:]\right|_{2}\qquad$
$\displaystyle\leq\left|\left(\tilde{L}H_{j-1}^{\prime}(W_{j}+U_{j})\right)[u_{j}^{\ast},:]-\left(\tilde{L}H_{j-1}W_{j}\right)[u_{j}^{\ast},:]\right|_{2}\qquad(\text{Lipschitz
property of ReLU})$
$\displaystyle=\left|\left((\tilde{L}H_{j-1}^{\prime})[u_{j}^{\ast},:]\right)(W_{j}+U_{j})-\left((\tilde{L}H_{j-1})[u_{j}^{\ast},:]\right)W_{j}\right|_{2}$
$\displaystyle=\left|\left(\left((\tilde{L}H_{j-1}^{\prime})[u_{j}^{\ast},:]\right)-\left((\tilde{L}H_{j-1})[u_{j}^{\ast},:]\right)\right)(W_{j}+U_{j})+\left((\tilde{L}H_{j-1})[u_{j}^{\ast},:]\right)U_{j}\right|_{2}$
$\displaystyle=\left|\left(\sum_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]\left(H_{j-1}^{\prime}[k,:]-H_{j-1}[k,:]\right)\right)(W_{j}+U_{j})+\left(\sum_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]H_{j-1}[k,:]\right)U_{j}\right|_{2}$
$\displaystyle\leq\left|\sum\nolimits_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]\left(H_{j-1}^{\prime}[k,:]-H_{j-1}[k,:]\right)\right|_{2}\left\|W_{j}+U_{j}\right\|_{2}+\left|\sum\nolimits_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]H_{j-1}[k,:]\right|_{2}\left\|U_{j}\right\|_{2}$
$\displaystyle\leq\sum\nolimits_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]\left|H_{j-1}^{\prime}[k,:]-H_{j-1}[k,:]\right|_{2}\left\|W_{j}+U_{j}\right\|_{2}+\sum\nolimits_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]\left|H_{j-1}[k,:]\right|_{2}\left\|U_{j}\right\|_{2}$
$\displaystyle\leq\sum\nolimits_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]\Psi_{j-1}\left\|W_{j}+U_{j}\right\|_{2}+\sum\nolimits_{k\in\mathcal{N}_{u_{j}^{\ast}}}\tilde{L}[u_{j}^{\ast},k]\Phi_{j-1}\left\|U_{j}\right\|_{2}$
$\displaystyle\leq\sqrt{d}\Psi_{j-1}\left\|W_{j}+U_{j}\right\|_{2}+\sqrt{d}\Phi_{j-1}\left\|U_{j}\right\|_{2},$
(29)

where in the second from the last inequality we use the fact $\forall k$,
$\left|H_{j-1}^{\prime}[k,:]-H_{j-1}[k,:]\right|_{2}\leq\Psi_{j-1}$ and
$\forall k$, $\left|H_{j-1}[k,:]\right|_{2}\leq\Phi_{j-1}$. In the last
inequality, we again use the fact in Eq. (28). We can simplify the notations
in Eq. (A.3) as $\Psi_{j}\leq a_{j-1}\Psi_{j-1}+b_{j-1}$ where
$a_{j-1}=\sqrt{d}\|W_{j}+U_{j}\|_{2}$ and
$b_{j-1}=\sqrt{d}\Phi_{j-1}\|U_{j}\|_{2}$. Since $\Delta_{0}=X-X=\bm{0}$, we
have $\Psi_{0}=0$. It is straightforward to work out the recursion as,

$\displaystyle\Psi_{j}$
$\displaystyle\leq\sum_{k=0}^{j-1}b_{k}\left(\prod_{i=k+1}^{j-1}a_{i}\right)=\sum_{k=0}^{j-1}d^{\frac{1}{2}}\Phi_{k}\|U_{k+1}\|_{2}\left(\prod_{i=k+1}^{j-1}d^{\frac{1}{2}}\|W_{i+1}+U_{i+1}\|_{2}\right)$
$\displaystyle=\sum_{k=0}^{j-1}d^{\frac{j-k}{2}}\Phi_{k}\|U_{k+1}\|_{2}\left(\prod_{i=k+2}^{j}\|W_{i}+U_{i}\|_{2}\right).$
(30)

Based on Eq. (A.3), we can instantiate the bound in Eq. (A.3) as

$\displaystyle\Psi_{j}$
$\displaystyle\leq\sum_{k=0}^{j-1}d^{\frac{j-k}{2}}\Phi_{k}\|U_{k+1}\|_{2}\left(\prod_{i=k+2}^{j}\|W_{i}+U_{i}\|_{2}\right)$
$\displaystyle\leq\sum_{k=0}^{j-1}d^{\frac{j-k}{2}}\left(d^{\frac{k}{2}}B\prod_{i=1}^{k}\|W_{i}\|_{2}\right)\|U_{k+1}\|_{2}\left(\prod_{i=k+2}^{j}\left(\|W_{i}\|_{2}+\|U_{i}\|_{2}\right)\right)$
$\displaystyle\leq
B\sum_{k=0}^{j-1}d^{\frac{j}{2}}\left(\prod_{i=1}^{k}\|W_{i}\|_{2}\right)\|U_{k+1}\|_{2}\left(\prod_{i=k+2}^{j}\left(1+\frac{1}{l}\right)\|W_{i}\|_{2}\right)$
$\displaystyle=B\sum_{k=0}^{j-1}d^{\frac{j}{2}}\left(\prod_{i=1}^{k+1}\|W_{i}\|_{2}\right)\frac{\|U_{k+1}\|_{2}}{\|W_{k+1}\|_{2}}\left(\prod_{i=k+2}^{j}\left(1+\frac{1}{l}\right)\|W_{i}\|_{2}\right)$
$\displaystyle=Bd^{\frac{j}{2}}\left(\prod_{i=1}^{j}\|W_{i}\|_{2}\right)\sum_{k=0}^{j-1}\frac{\|U_{k+1}\|_{2}}{\|W_{k+1}\|_{2}}\left(1+\frac{1}{l}\right)^{j-k-1}$
$\displaystyle\leq
Bd^{\frac{j}{2}}\left(\prod_{i=1}^{j}\|W_{i}\|_{2}\right)\sum_{k=1}^{j}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\left(1+\frac{1}{l}\right)^{j-k}$
(31)

##### Final Bound on the Readout Layer

Now let us consider the average readout function in the last layer, _i.e_.,
the $l$-th layer. Based on Eq. (A.3) and Eq. (A.3), we can bound the change of
GCN’s output with and without the weight perturbation as follows,

$\displaystyle|\Delta_{l}|_{2}=$
$\displaystyle\left|\frac{1}{n}\bm{1}_{n}H_{l-1}^{\prime}(W_{l}+U_{l})-\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}\right|_{2}$
$\displaystyle=$
$\displaystyle\left|\frac{1}{n}\bm{1}_{n}\Delta_{l-1}(W_{l}+U_{l})+\frac{1}{n}\bm{1}_{n}H_{l-1}U_{l}\right|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\left|\bm{1}_{n}\Delta_{l-1}(W_{l}+U_{l})\right|_{2}+\frac{1}{n}\left|\bm{1}_{n}H_{l-1}U_{l}\right|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\|W_{l}+U_{l}\|_{2}|\bm{1}_{n}\Delta_{l-1}|_{2}+\frac{1}{n}\|U_{l}\|_{2}|\bm{1}_{n}H_{l-1}|_{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\left\|W_{l}+U_{l}\right\|_{2}\left|\sum_{i=1}^{n}\Delta_{l-1}[i,:]\right|_{2}+\frac{1}{n}\left\|U_{l}\right\|_{2}\left|\sum_{i=1}^{n}H_{l-1}[i,:]\right|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\left\|W_{l}+U_{l}\right\|_{2}\left(\sum_{i=1}^{n}\left|\Delta_{l-1}[i,:]\right|_{2}\right)+\frac{1}{n}\left\|U_{l}\right\|_{2}\left(\sum_{i=1}^{n}\left|H_{l-1}[i,:]\right|_{2}\right)$
$\displaystyle\leq$
$\displaystyle\left\|W_{l}+U_{l}\right\|_{2}\Psi_{l-1}+\left\|U_{l}\right\|_{2}\Phi_{l-1}$
$\displaystyle\leq$
$\displaystyle\left\|W_{l}+U_{l}\right\|_{2}Bd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l-1}\|W_{i}\|_{2}\right)\sum_{k=1}^{l-1}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\left(1+\frac{1}{l}\right)^{l-1-k}+\left\|U_{l}\right\|_{2}Bd^{\frac{l-1}{2}}\prod_{i=1}^{l-1}\|W_{i}\|_{2}$
$\displaystyle=$ $\displaystyle
Bd^{\frac{l-1}{2}}\left[\left\|W_{l}+U_{l}\right\|_{2}\left(\prod_{i=1}^{l-1}\|W_{i}\|_{2}\right)\sum_{k=1}^{l-1}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\left(1+\frac{1}{l}\right)^{l-1-k}+\left\|U_{l}\right\|_{2}\prod_{i=1}^{l-1}\|W_{i}\|_{2}\right]$
$\displaystyle=$ $\displaystyle
Bd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\left[\frac{\left\|W_{l}+U_{l}\right\|_{2}}{\left\|W_{l}\right\|_{2}}\sum_{k=1}^{l-1}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\left(1+\frac{1}{l}\right)^{l-1-k}+\frac{\left\|U_{l}\right\|_{2}}{\left\|W_{l}\right\|_{2}}\right]$
$\displaystyle\leq$ $\displaystyle
Bd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\left[\left(1+\frac{1}{l}\right)\sum_{k=1}^{l-1}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\left(1+\frac{1}{l}\right)^{l-1-k}+\frac{\left\|U_{l}\right\|_{2}}{\left\|W_{l}\right\|_{2}}\right]$
$\displaystyle=$ $\displaystyle
Bd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\left(1+\frac{1}{l}\right)^{l}\left[\sum_{k=1}^{l-1}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\left(1+\frac{1}{l}\right)^{-k}+\frac{\left\|U_{l}\right\|_{2}}{\left\|W_{l}\right\|_{2}}\left(1+\frac{1}{l}\right)^{-l}\right]$
$\displaystyle\leq$ $\displaystyle
eBd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\left[\sum_{k=1}^{l}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}\right]\qquad\left(\text{Use
}1\leq\left(1+\frac{1}{l}\right)^{l}\leq e\right)$ (32)

which proves the lemma. ∎

######  Theorem 3.2.

(GCN Generalization Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-layer GCN. Then for any $\delta,\gamma>0$, with probability at least
$1-\delta$ over the choice of an i.i.d. size-$m$ training set $S$ according to
$\mathcal{D}$, for any $w$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}d^{l-1}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}+\log\frac{ml}{\delta}}{\gamma^{2}m}}\right)$

###### Proof.

Let $\beta=\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)^{1/l}$. We normalize the
weights as $\tilde{W}_{i}=\frac{\beta}{\|W_{i}\|_{2}}W_{i}$. Due to the
homogeneity of ReLU, _i.e_., $a\phi(x)=\phi(ax)$, $\forall a\geq 0$, we have
$f_{w}=f_{\tilde{w}}$. We can also verify that
$\prod_{i=1}^{l}\|W_{i}\|_{2}=\prod_{i=1}^{l}\|\tilde{W}_{i}\|_{2}$ and
$\|W_{i}\|_{{F}}/\|W_{i}\|_{2}=\|\tilde{W}_{i}\|_{{F}}/\|\tilde{W}_{i}\|_{2}$,
_i.e_., the terms appear in the bound stay the same after applying the
normalization. Therefore, w.l.o.g., we assume that the norm is equal across
layers, _i.e_., $\forall i$, $\|W_{i}\|_{2}=\beta$.

Consider the prior $P=\mathcal{N}(0,\sigma^{2}I)$ and the random perturbation
$u\sim\mathcal{N}(0,\sigma^{2}I)$. Note that the $\sigma$ of the prior and the
perturbation are the same and will be set according to $\beta$. More
precisely, we will set the $\sigma$ based on some approximation
$\tilde{\beta}$ of $\beta$ since the prior $P$ can not depend on any learned
weights directly. The approximation $\tilde{\beta}$ is chosen to be a cover
set which covers the meaningful range of $\beta$. For now, let us assume that
we have a fix $\tilde{\beta}$ and consider $\beta$ which satisfies
$|\beta-\tilde{\beta}|\leq\frac{1}{l}\beta$. Note that this also implies

$\displaystyle|\beta-\tilde{\beta}|\leq\frac{1}{l}\beta$
$\displaystyle~{}\Rightarrow~{}\left(1-\frac{1}{l}\right)\beta\leq\tilde{\beta}\leq\left(1+\frac{1}{l}\right)\beta$
$\displaystyle~{}\Rightarrow~{}\left(1-\frac{1}{l}\right)^{l-1}\beta^{l-1}\leq\tilde{\beta}^{l-1}\leq\left(1+\frac{1}{l}\right)^{l-1}\beta^{l-1}$
$\displaystyle~{}\Rightarrow~{}\left(1-\frac{1}{l}\right)^{l}\beta^{l-1}\leq\tilde{\beta}^{l-1}\leq\left(1+\frac{1}{l}\right)^{l}\beta^{l-1}$
$\displaystyle~{}\Rightarrow~{}\frac{1}{e}\beta^{l-1}\leq\tilde{\beta}^{l-1}\leq
e\beta^{l-1}$ (33)

From Tropp (2012), for $U_{i}\in\mathbb{R}^{h\times h}$ and
$U_{i}\sim\mathcal{N}(\bm{0},\sigma^{2}I)$, we have,

$\displaystyle\mathbb{P}\left(\|U_{i}\|_{2}\geq t\right)\leq
2he^{-t^{2}/2h\sigma^{2}}.$ (34)

Taking a union bound, we have

$\displaystyle\mathbb{P}\left(\|U_{1}\|_{2}<t~{}\&~{}\cdots~{}\&~{}\|U_{l}\|_{2}<t\right)$
$\displaystyle=1-\mathbb{P}\left(\exists i,\|U_{i}\|_{2}\geq t\right)$
$\displaystyle\geq 1-\sum_{i=1}^{l}\mathbb{P}\left(\|U_{i}\|_{2}\geq t\right)$
$\displaystyle\geq 1-2lhe^{-t^{2}/2h\sigma^{2}}.$ (35)

Setting $2lhe^{-t^{2}/2h\sigma^{2}}=\frac{1}{2}$, we have
$t=\sigma\sqrt{2h\log(4lh)}$. This implies that the probability that the
spectral norm of the perturbation of any layer is no larger than
$\sigma\sqrt{2h\log(4lh)}$ holds with probability at least $\frac{1}{2}$.
Plugging this bound into Lemma 3.1, we have with probability at least
$\frac{1}{2}$,

$\displaystyle\left|f_{w+u}(X,A)-f_{w}(X,A)\right|_{2}$ $\displaystyle\leq
eBd^{\frac{l-1}{2}}\left(\prod_{i=1}^{l}\|W_{i}\|_{2}\right)\sum_{k=1}^{l}\frac{\|U_{k}\|_{2}}{\|W_{k}\|_{2}}$
$\displaystyle=eBd^{\frac{l-1}{2}}\beta^{l}\sum_{k=1}^{l}\frac{\|U_{k}\|_{2}}{\beta}$
$\displaystyle\leq eBd^{\frac{l-1}{2}}\beta^{l-1}l\sigma\sqrt{2h\log(4lh)}$
$\displaystyle\leq
e^{2}Bd^{\frac{l-1}{2}}\tilde{\beta}^{l-1}l\sigma\sqrt{2h\log(4lh)}\leq\frac{\gamma}{4},$
(36)

where we can set
$\sigma=\frac{\gamma}{42Bd^{\frac{l-1}{2}}\tilde{\beta}^{l-1}l\sqrt{h\log(4lh)}}$
to get the last inequality. Note that Lemma 3.1 also requires $\forall
i\in\mathbb{N}^{+}_{l}$, $\|U_{i}\|_{2}\leq\frac{1}{l}\|W_{i}\|_{2}$. The
requirement is satisfied if $\sigma\leq\frac{\beta}{l\sqrt{2h\log(4lh)}}$
which in turn can be satisfied if

$\displaystyle\frac{\gamma}{4eBd^{\frac{l-1}{2}}\beta^{l-1}l\sqrt{2h\log(4lh)}}\leq\frac{\beta}{l\sqrt{2h\log(4lh)}},$
(37)

since the chosen value of $\sigma$ satisfies
$\sigma\leq\frac{\gamma}{4eBd^{\frac{l-1}{2}}\beta^{l-1}l\sqrt{2h\log(4lh)}}$.
Note that Eq. (37) is equivalent to
$\frac{\gamma}{4eB}d^{\frac{1-l}{2}}\leq\beta^{l}$. We will see how to satisfy
this condition later.

We now compute the KL term in the PAC-Bayes bound in Lemma 2.2.

$\displaystyle\text{KL}\left(Q\|P\right)$
$\displaystyle=\frac{|w|_{2}^{2}}{2\sigma^{2}}=\frac{42^{2}B^{2}d^{l-1}\tilde{\beta}^{2l-2}l^{2}h\log(4lh)}{2\gamma^{2}}\sum_{i=1}^{l}\|W_{i}\|_{F}^{2}$
$\displaystyle\leq\mathcal{O}\left(\frac{B^{2}d^{l-1}\beta^{2l}l^{2}h\log(lh)}{\gamma^{2}}\sum_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\beta^{2}}\right)$
$\displaystyle\leq\mathcal{O}\left(B^{2}d^{l-1}l^{2}h\log(lh)\frac{\prod_{i=1}^{l}\|W_{i}\|_{2}^{2}}{\gamma^{2}}\sum_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}\right).$
(38)

From Lemma 2.2, fixing any $\tilde{\beta}$, with probability $1-\delta$ and
for all $w$ such that $|\beta-\tilde{\beta}|\leq\frac{1}{l}\beta$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}d^{l-1}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}+\log\frac{m}{\delta}}{\gamma^{2}m}}\right).$
(39)

Finally, we need to consider multiple choices of $\tilde{\beta}$ so that for
any $\beta$, we can bound the generalization error like Eq. (39). First, we
only need to consider values of $\beta$ in the following range,

$\displaystyle\frac{1}{\sqrt{d}}\left(\frac{\gamma\sqrt{d}}{2B}\right)^{1/l}\leq\beta\leq\frac{1}{\sqrt{d}}\left(\frac{\gamma\sqrt{md}}{2B}\right)^{1/l},$
(40)

since otherwise the bound holds trivially as $L_{\mathcal{D},0}(f_{w})\leq 1$
by definition. Note that the lower bound in Eq. (40) ensures that Eq. (37)
holds which in turn justifies the applicability of Lemma 3.1. If
$\beta<\frac{1}{\sqrt{d}}\left(\frac{\gamma\sqrt{d}}{2B}\right)^{1/l}$, then
for any $(X,A)$ and any $j\in\mathbb{N}^{+}_{K}$,
$|f(X,A)[j]|\leq\frac{\gamma}{2}$. To see this, we have,

$\displaystyle\left|f_{w}(X,A)[j]\right|$
$\displaystyle\leq\left|f_{w}(X,A)\right|_{2}=|\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}|_{2}$
$\displaystyle\leq\frac{1}{n}|\bm{1}_{n}H_{l-1}|_{2}\|W_{l}\|_{2}$
$\displaystyle\leq\|W_{l}\|_{2}\max_{i}|H_{l-1}[i,:]|_{2}$ $\displaystyle\leq
Bd^{\frac{l-1}{2}}\prod_{i=1}^{l}\|W_{i}\|_{2}=d^{\frac{l-1}{2}}\beta^{l}B\qquad(\text{Use
Eq. (\ref{eq:gcn_output_bound})})$
$\displaystyle=d^{\frac{l-1}{2}}B\frac{\gamma}{2Bd^{\frac{l-1}{2}}}\leq\frac{\gamma}{2}.$
(41)

Therefore, by the definition in Eq. (4), we always have
$L_{S,\gamma}(f_{w})=1$ when
$\beta<\frac{1}{\sqrt{d}}\left(\frac{\gamma\sqrt{d}}{2B}\right)^{1/l}$.
Alternatively, if
$\beta>\frac{1}{\sqrt{d}}\left(\frac{\gamma\sqrt{md}}{2B}\right)^{1/l}$, the
term inside the big-O notation in Eq. (39) would be,

$\displaystyle\sqrt{\frac{B^{2}d^{l-1}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}+\log\frac{m}{\delta}}{\gamma^{2}m}}$
$\displaystyle\geq\sqrt{\frac{l^{2}h\log(lh)}{4}\sum\limits_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}}$
$\displaystyle\geq\sqrt{\frac{l^{2}h\log(lh)}{4}}\geq 1,$ (42)

where we use the facts that $\|W_{i}\|_{F}\geq\|W_{i}\|_{2}$ and we typically
choose $h\geq 2$ in practice and $l\geq 2$.

Since we only need to consider $\beta$ in the range of Eq. (40), a sufficient
condition to make $|\beta-\tilde{\beta}|\leq\frac{1}{l}\beta$ hold would be
$|\beta-\tilde{\beta}|\leq\frac{1}{l\sqrt{d}}\left(\frac{\gamma\sqrt{d}}{2B}\right)^{1/l}$.
Therefore, if we can find a covering of the interval in Eq. (40) with radius
$\frac{1}{l\sqrt{d}}\left(\frac{\gamma\sqrt{d}}{2B}\right)^{1/l}$ and make
sure bounds like Eq. (39) holds while $\tilde{\beta}$ takes all possible
values from the covering, then we can get a bound which holds for all $\beta$.
It is clear that we only need to consider a covering $C$ with size
$|C|=\frac{l}{2}\left(m^{\frac{1}{2l}}-1\right)$. Therefore, denoting the
event of Eq. (39) with $\tilde{\beta}$ taking the $i$-th value of the covering
as $E_{i}$, we have

$\displaystyle\mathbb{P}\left(E_{1}~{}\&~{}\cdots~{}\&~{}E_{|C|}\right)$
$\displaystyle=1-\mathbb{P}\left(\exists i,\bar{E}_{i}\right)\geq
1-\sum_{i=1}^{|C|}\mathbb{P}\left(\bar{E}_{i}\right)\geq 1-|C|\delta.$ (43)

Note $\bar{E}_{i}$ denotes the complement of $E_{i}$. Hence, with probability
$1-\delta$ and for all $w$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})$ $\displaystyle\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}d^{l-1}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}+\log\frac{m|C|}{\delta}}{\gamma^{2}m}}\right)$
$\displaystyle=L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}d^{l-1}l^{2}h\log(lh)\prod\limits_{i=1}^{l}\|W_{i}\|_{2}^{2}\sum\limits_{i=1}^{l}\frac{\|W_{i}\|_{F}^{2}}{\|W_{i}\|_{2}^{2}}+\log\frac{ml}{\delta}}{\gamma^{2}m}}\right),$
(44)

which proves the theorem. ∎

###  A.4 MPGNNs Results

In this part, we provide the proofs of the main results regarding MPGNNs.

######  Lemma 3.3.

(MPGNN Perturbation Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-step MPGNN. Then for any $w$, and $x\in\mathcal{X}_{B,h_{0}}$, and any
perturbation $u=\text{vec}(\\{U_{1},U_{2},U_{l}\\})$ such that
$\eta=\max\left(\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}},\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}},\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\right)\leq\frac{1}{l}$,
the change in the output of MPGNN is bounded as,

$\displaystyle|f_{w+u}($ $\displaystyle
X,A)-f_{w}(X,A)|_{2}\leq\begin{cases}eB\left(l+1\right)^{2}\eta\|W_{1}\|_{2}\|W_{l}\|_{2}C_{\phi},&\text{if}\
d\mathcal{C}=1\\\
eBl\eta\|W_{1}\|_{2}\|W_{l}\|_{2}C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1},&\text{otherwise}\end{cases}$

where $\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$.

###### Proof.

We first perform the recursive perturbation analysis on node representations
of all steps except the last one, i.e., the readout step. Then we derive the
bound for the graph representation of the last readout step.

##### Perturbation Analysis on Node Representations.

In message passing GNNs, for any step $j<l$ besides the last readout one, the
node representations are,

$\displaystyle\bar{M}_{j}$
$\displaystyle=C_{\operatorname{in}}g\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)$
$\displaystyle H_{j}$
$\displaystyle=\phi\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right),$
(45)

where the incidence matrices $C_{\operatorname{in}}\in\mathbb{R}^{n\times c}$
and $C_{\operatorname{out}}\in\mathbb{R}^{n\times c}$ (recall $c$ is the
number of edges). Moreover, since each edge only connects one incoming and one
outgoing node, we have,

$\displaystyle\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]$
$\displaystyle\leq\max_{i}\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]=\|C_{\operatorname{in}}\|_{\infty}\leq
d$ $\displaystyle\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]$
$\displaystyle\leq\max_{k}\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]\leq\|C_{\operatorname{out}}\|_{1}\leq
1$ (46)

where $d-1$ is the maximum node degree. Note that one actually has
$\|C_{\operatorname{in}}\|_{\infty}\leq d-1$ for simple graphs. Since some
models in the literature pre-process the graphs by adding self-loops, we thus
relax it to $\|C_{\operatorname{in}}\|_{\infty}\leq d$ which holds in both
cases.

We add perturbation $u$ to the weights $w$, _i.e_., the perturbed weights are
$W_{1}+U_{1}$, $W_{2}+U_{2}$ and $W_{l}+U_{l}$. For the ease of notation, we
use the superscript of prime to denote the perturbed node representations,
_e.g_., $H_{j}^{\prime}=f_{w+u}^{j}(X,A)$. Let
$\Delta_{j}=f_{w+u}^{j}(X,A)-f_{w}^{j}(X,A)=H_{j}^{\prime}-H_{j}$. Note that
$\Delta_{j}\in\mathbb{R}^{n\times h_{j}}$. Let
$\Psi_{j}=\max\limits_{i}\left|\Delta_{j}[i,:]\right|_{2}=\max\limits_{i}\left|H_{j}^{\prime}[i,:]-H_{j}[i,:]\right|_{2}$
and $\Phi_{j}=\max\limits_{i}\left|H_{j}[i,:]\right|_{2}$. We denote the
$u_{j}^{\ast}=\operatorname*{arg\,max}\limits_{i}\left|\Delta_{j}[i,:]\right|_{2}$
and
$v_{j}^{\ast}=\operatorname*{arg\,max}\limits_{i}\left|H_{j}[i,:]\right|_{2}$.
To simplify the derivation, we abbreviate the following statistics
$\kappa=C_{\phi}B\|W_{1}\|_{2}$ and $\tau=d\mathcal{C}$ throughout the proof
where $\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$ is the _percolation
complexity_.

##### Upper Bound on the Max Node Representation.

For any step $j<l$, we can derive an upper bound on the $\ell_{2}$ norm of the
aggregated message of any node $i$ as follows,

$\displaystyle\left|\bar{M}_{j}[i,:]\right|_{2}$
$\displaystyle=\left|\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]\left(g\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)\right)[k,:]\right|_{2}=\left|\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]g\left(C_{\operatorname{out}}^{\top}[k,:]H_{j-1}\right)\right|_{2}$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]\left|g\left(C_{\operatorname{out}}^{\top}[k,:]H_{j-1}\right)\right|_{2}$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left|C_{\operatorname{out}}^{\top}[k,:]H_{j-1}\right|_{2}=\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left|\sum_{t=1}^{n}C_{\operatorname{out}}^{\top}[k,t]H_{j-1}[t,:]\right|_{2}$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left(\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]\left|H_{j-1}[t,:]\right|_{2}\right)$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left(\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]\Phi_{j-1}\right)$
$\displaystyle\leq dC_{g}\Phi_{j-1}.$ (47)

Then we can derive an upper bound on the maximum (w.r.t. $\ell_{2}$ norm) node
representation as follows,

$\displaystyle\Phi_{j}$
$\displaystyle=\max\limits_{i}\left|H_{j}[i,:]\right|_{2}=\left|\phi\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right)[v_{j}^{\ast},:]\right|_{2}$
$\displaystyle=\left|\phi\left(\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right)[v_{j}^{\ast},:]\right)\right|_{2}$
$\displaystyle\leq
C_{\phi}\left|\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right)[v_{j}^{\ast},:]\right|_{2}$
$\displaystyle=C_{\phi}\left|\left(XW_{1}\right)[v_{j}^{\ast},:]+\left(\rho\left(\bar{M}_{j}\right)W_{2}\right)[v_{j}^{\ast},:]\right|_{2}$
$\displaystyle\leq
C_{\phi}\left|\left(XW_{1}\right)[v_{j}^{\ast},:]\right|_{2}+C_{\phi}\left|\left(\rho\left(\bar{M}_{j}\right)W_{2}\right)[v_{j}^{\ast},:]\right|_{2}$
$\displaystyle=C_{\phi}\left|X[v_{j}^{\ast},:]W_{1}\right|_{2}+C_{\phi}\left|\rho\left(\bar{M}_{j}\right)[v_{j}^{\ast},:]W_{2}\right|_{2}$
$\displaystyle\leq
C_{\phi}\left|X[v_{j}^{\ast},:]\right|_{2}\|W_{1}\|_{2}+C_{\phi}\left|\rho\left(\bar{M}_{j}\right)[v_{j}^{\ast},:]\right|_{2}\|W_{2}\|_{2}$
$\displaystyle\leq
C_{\phi}B\|W_{1}\|_{2}+C_{\phi}\left|\rho\left(\bar{M}_{j}[v_{j}^{\ast},:]\right)\right|_{2}\|W_{2}\|_{2}$
$\displaystyle\leq
C_{\phi}B\|W_{1}\|_{2}+C_{\phi}C_{\rho}\left|\bar{M}_{j}[v_{j}^{\ast},:]\right|_{2}\|W_{2}\|_{2}$
$\displaystyle\leq
C_{\phi}B\|W_{1}\|_{2}+dC_{\phi}C_{\rho}C_{g}\Phi_{j-1}\|W_{2}\|_{2}=\kappa+\tau\Phi_{j-1}$
$\displaystyle\leq\tau^{j}\Phi_{0}+\sum_{i=0}^{j-1}\tau^{j-1-i}\kappa\qquad(\text{Unroll
recursion})$ $\displaystyle=\sum_{i=0}^{j-1}\tau^{j-1-i}\kappa\qquad(\text{Use
}\Phi_{0}=0)$ $\displaystyle=\begin{cases}j\kappa,&\text{if}\ \tau=1\\\
\kappa\frac{\tau^{j}-1}{\tau-1},&\text{otherwise}\end{cases}$ (48)

##### Upper Bound on the Max Change of Node Representation.

For any step $j<l$, we can derive an upper bound on the maximum (w.r.t.
$\ell_{2}$ norm) change between the aggregated message with and without the
weight perturbation for any node $i$ as follows,

$\displaystyle\left|\bar{M}_{j}^{\prime}[i,:]-\bar{M}_{j}[i,:]\right|_{2}$
$\displaystyle\leq\left|\left(C_{\operatorname{in}}g\left(C_{\operatorname{out}}^{\top}H_{j-1}^{\prime}\right)\right)[i,:]-\left(C_{\operatorname{in}}g\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)\right)[i,:]\right|_{2}$
$\displaystyle=\left|\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]\left(g\left(C_{\operatorname{out}}^{\top}H_{j-1}^{\prime}\right)-g\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)\right)[k,:]\right|_{2}$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]\left|\left(g\left(C_{\operatorname{out}}^{\top}H_{j-1}^{\prime}\right)-g\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)\right)[k,:]\right|_{2}$
$\displaystyle=\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]\left|g\left(\left(C_{\operatorname{out}}^{\top}H_{j-1}^{\prime}\right)[k,:]\right)-g\left(\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)[k,:]\right)\right|_{2}$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left|\left(C_{\operatorname{out}}^{\top}H_{j-1}^{\prime}\right)[k,:]-\left(C_{\operatorname{out}}^{\top}H_{j-1}\right)[k,:]\right|_{2}$
$\displaystyle=\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left|\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]H_{j-1}^{\prime}[t,:]-\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]H_{j-1}[t,:]\right|_{2}$
$\displaystyle=\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left|\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]\left(H_{j-1}^{\prime}[t,:]-H_{j-1}[t,:]\right)\right|_{2}$
$\displaystyle\leq\sum_{k=1}^{c}C_{\operatorname{in}}[i,k]C_{g}\left(\sum_{t=1}^{n}C_{\operatorname{out}}[t,k]\left|H_{j-1}^{\prime}[t,:]-H_{j-1}[t,:]\right|_{2}\right)$
$\displaystyle\leq dC_{g}\Psi_{j-1}$ (49)

Based on Eq. (A.4), we can derive an upper bound on the maximum (w.r.t.
$\ell_{2}$ norm) change between the representations with and without the
weight perturbation for any node as follows,

$\displaystyle\Psi_{j}=$
$\displaystyle\max\limits_{i}\left|H_{j}^{\prime}[i,:]-H_{j}[i,:]\right|_{2}$
$\displaystyle=$
$\displaystyle\left|\left(\phi\left(X(W_{1}+U_{1})+\rho\left(\bar{M}_{j}^{\prime}\right)(W_{2}+U_{2})\right)\right)[u_{j}^{\ast},:]-\left(\phi\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right)\right)[u_{j}^{\ast},:]\right|_{2}$
$\displaystyle=$
$\displaystyle\left|\phi\left(\left(X(W_{1}+U_{1})+\rho\left(\bar{M}_{j}^{\prime}\right)(W_{2}+U_{2})\right)[u_{j}^{\ast},:]\right)-\phi\left(\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right)[u_{j}^{\ast},:]\right)\right|_{2}$
$\displaystyle\leq$ $\displaystyle
C_{\phi}\left|\left(X(W_{1}+U_{1})+\rho\left(\bar{M}_{j}^{\prime}\right)(W_{2}+U_{2})\right)[u_{j}^{\ast},:]-\left(XW_{1}+\rho\left(\bar{M}_{j}\right)W_{2}\right)[u_{j}^{\ast},:]\right|_{2}$
$\displaystyle\leq$ $\displaystyle
C_{\phi}\left|X[u_{j}^{\ast},:]U_{1}+\left(\rho\left(\bar{M}_{j}^{\prime}\right)\right)[u_{j}^{\ast},:](W_{2}+U_{2})-\left(\rho\left(\bar{M}_{j}\right)\right)[u_{j}^{\ast},:]W_{2}\right|_{2}$
$\displaystyle=$ $\displaystyle
C_{\phi}\left|X[u_{j}^{\ast},:]U_{1}+\left(\rho\left(\bar{M}_{j}^{\prime}\right)-\rho\left(\bar{M}_{j}\right)\right)[u_{j}^{\ast},:](W_{2}+U_{2})+\rho\left(\bar{M}_{j}\right)[u_{j}^{\ast},:]U_{2}\right|_{2}$
$\displaystyle\leq$ $\displaystyle
C_{\phi}B\left\|U_{1}\right\|_{2}+C_{\phi}C_{\rho}\left|\bar{M}_{j}^{\prime}[u_{j}^{\ast},:]-\bar{M}_{j}[u_{j}^{\ast},:]\right|_{2}\left\|W_{2}+U_{2}\right\|_{2}+C_{\phi}C_{\rho}\left|\bar{M}_{j}[u_{j}^{\ast},:]\right|_{2}\left\|U_{2}\right\|_{2}$
$\displaystyle\leq$
$\displaystyle\kappa\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}}+d\mathcal{C}\Psi_{j-1}\frac{\|W_{2}+U_{2}\|_{2}}{\|W_{2}\|_{2}}+d\mathcal{C}\Phi_{j-1}\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}\qquad(\text{Use
Eq. (\ref{eq:mpgnn_msg_bound}) and (\ref{eq:mpgnn_msg_diff_bound})})$
$\displaystyle\leq$
$\displaystyle\tau\left(1+\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}\right)\Psi_{j-1}+\kappa\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}}+\tau\Phi_{j-1}\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}$
(50)

If $\tau=1$, then we have,

$\displaystyle\Psi_{j}\leq$
$\displaystyle\tau\left(1+\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}\right)\Psi_{j-1}+\kappa\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}}+\tau\Phi_{j-1}\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}$
$\displaystyle\leq$
$\displaystyle\left(1+\frac{\left\|U_{2}\right\|_{2}}{\left\|W_{2}\right\|_{2}}\right)\Psi_{j-1}+\kappa\left(\frac{\left\|U_{1}\right\|_{2}}{\left\|W_{1}\right\|_{2}}+\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}(j-1)\right)\qquad(\text{Use
Eq. (\ref{eq:mpgnn_output_bound})})$ $\displaystyle\leq$
$\displaystyle\left(1+\eta\right)\Psi_{j-1}+\kappa\eta\left(1+(j-1)\right)\qquad\left(\text{Use
}\eta=\max\left(\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}},\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}},\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\right)\right)$
$\displaystyle=$ $\displaystyle\left(1+\eta\right)\Psi_{j-1}+\kappa\eta j.$
(51)

If $\tau\neq 1$, then we have,

$\displaystyle\Psi_{j}\leq$
$\displaystyle\tau\left(1+\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}\right)\Psi_{j-1}+\kappa\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}}+\tau\Phi_{j-1}\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}}$
$\displaystyle\leq$
$\displaystyle\tau\left(1+\frac{\left\|U_{2}\right\|_{2}}{\left\|W_{2}\right\|_{2}}\right)\Psi_{j-1}+\kappa\left(\frac{\left\|U_{1}\right\|_{2}}{\left\|W_{1}\right\|_{2}}+\tau\frac{\left\|U_{2}\right\|_{2}}{\left\|W_{2}\right\|_{2}}\frac{\tau^{j-1}-1}{\tau-1}\right)\qquad(\text{Use
Eq. (\ref{eq:mpgnn_output_bound})})$ $\displaystyle\leq$
$\displaystyle\tau\left(1+\eta\right)\Psi_{j-1}+\kappa\eta\left(1+\frac{\tau^{j}-\tau}{\tau-1}\right)\qquad\left(\text{Use
}\eta=\max\left(\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}},\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}},\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\right)\right)$
$\displaystyle=$
$\displaystyle\tau\left(1+\eta\right)\Psi_{j-1}+\kappa\eta\left(\frac{\tau^{j}-1}{\tau-1}\right).$
(52)

Recall from Eq. (A.3), if $\Psi_{j}\leq a_{j-1}\Psi_{j-1}+b_{j-1}$ and
$\Psi_{0}=0$, then
$\Psi_{j}\leq\sum_{k=0}^{j-1}b_{k}\left(\prod_{i=k+1}^{j-1}a_{i}\right)$. If
$\tau=1$, then we have $a_{j-1}=1+\eta$, $b_{j-1}=\kappa\eta j$ in our case
and,

$\displaystyle\Psi_{j}\leq$
$\displaystyle\sum\nolimits_{k=0}^{j-1}b_{k}\left(\prod\nolimits_{i=k+1}^{j-1}a_{i}\right)=\sum\nolimits_{k=0}^{j-1}\kappa\eta(k+1)\left(1+\eta\right)^{j-k-1}$
$\displaystyle\leq$
$\displaystyle\kappa\eta\left(1+\frac{1}{l}\right)^{j}\sum_{k=0}^{j-1}(k+1)\left(1+\frac{1}{l}\right)^{-k-1}\qquad\left(\text{Use
}\eta\leq\frac{1}{l}\right)$ $\displaystyle=$
$\displaystyle\kappa\eta\left(1+\frac{1}{l}\right)^{j}\sum_{k=1}^{j}k\left(1+\frac{1}{l}\right)^{-k}$
$\displaystyle=$
$\displaystyle\kappa\eta\left(1+\frac{1}{l}\right)^{j}\left(1+\frac{1}{l}\right)^{-1}\frac{1-(j+1)\left(1+\frac{1}{l}\right)^{-j}+j\left(1+\frac{1}{l}\right)^{-j-1}}{\left(1-\left(1+\frac{1}{l}\right)^{-1}\right)^{2}}$
$\displaystyle=$
$\displaystyle\kappa\eta\frac{\left(1+\frac{1}{l}\right)^{j+1}-(j+1)\left(1+\frac{1}{l}\right)+j}{\left(\left(1+\frac{1}{l}\right)^{1}-1\right)^{2}}$
$\displaystyle=$ $\displaystyle\kappa\eta
l^{2}\left(\left(1+\frac{1}{l}\right)^{j+1}-(j+1)\left(1+\frac{1}{l}\right)+j\right)$
$\displaystyle\leq$ $\displaystyle\kappa\eta
l^{2}\left(\left(1+\frac{1}{l}\right)^{j+1}-1\right)$ $\displaystyle\leq$
$\displaystyle\kappa\eta l\left(l+1\right)\left(1+\frac{1}{l}\right)^{j}$ (53)

If $\tau\neq 1$, then we have $a_{j-1}=\tau\left(1+\eta\right)$,
$b_{j-1}=\kappa\eta\left(\frac{\tau^{j}-1}{\tau-1}\right)$ in our case and,

$\displaystyle\Psi_{j}\leq$
$\displaystyle\sum_{k=0}^{j-1}b_{k}\left(\prod_{i=k+1}^{j-1}a_{i}\right)=\sum_{k=0}^{j-1}\kappa\eta\left(\frac{\tau^{k+1}-1}{\tau-1}\right)\tau^{j-k-1}\left(1+\eta\right)^{j-k-1}$
$\displaystyle\leq$
$\displaystyle\kappa\eta\tau^{j}\left(1+\frac{1}{l}\right)^{j}\sum_{k=0}^{j-1}\left(\frac{\tau^{k+1}-1}{\tau-1}\right)\tau^{-k-1}\left(1+\frac{1}{l}\right)^{-k-1}\qquad\left(\text{Use
}\eta\leq\frac{1}{l}\right)$ $\displaystyle\leq$
$\displaystyle\kappa\eta\tau^{j}\left(1+\frac{1}{l}\right)^{j}\sum_{k=0}^{j-1}\left(\frac{1-\tau^{-k-1}}{\tau-1}\right)\left(1+\frac{1}{l}\right)^{-k-1}$
$\displaystyle\leq$
$\displaystyle\frac{\kappa\eta\tau^{j}}{\tau-1}\left(1+\frac{1}{l}\right)^{j}\sum_{k=0}^{j-1}\left(1-\tau^{-k-1}\right)\left(1+\frac{1}{l}\right)^{-k-1}$
$\displaystyle\leq$
$\displaystyle\frac{\kappa\eta\tau^{j}}{\tau-1}\left(1+\frac{1}{l}\right)^{j}\sum_{k=1}^{j}\left(1-\tau^{-k}\right)$
(54)

##### Final Bound with Readout Function

Now let us consider the readout function. Since the last readout layer
produces a vector in $\mathbb{R}^{1\times C}$, we have,

$\displaystyle|\Delta_{l}|_{2}=$
$\displaystyle\left|\frac{1}{n}\bm{1}_{n}H_{l-1}^{\prime}(W_{l}+U_{l})-\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}\right|_{2}$
$\displaystyle=$
$\displaystyle\left|\frac{1}{n}\bm{1}_{n}\Delta_{l-1}(W_{l}+U_{l})+\frac{1}{n}\bm{1}_{n}H_{l-1}U_{l}\right|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\left|\bm{1}_{n}\Delta_{l-1}(W_{l}+U_{l})\right|_{2}+\frac{1}{n}\left|\bm{1}_{n}H_{l-1}U_{l}\right|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\|W_{l}+U_{l}\|_{2}|\bm{1}_{n}\Delta_{l-1}|_{2}+\frac{1}{n}\|U_{l}\|_{2}|\bm{1}_{n}H_{l-1}|_{2}$
$\displaystyle\leq$
$\displaystyle\|W_{l}+U_{l}\|_{2}\Psi_{l-1}+\|U_{l}\|_{2}\Phi_{l-1}$ (55)

If $\tau=1$, we have,

$\displaystyle|\Delta_{l}|_{2}\leq$
$\displaystyle\|W_{l}\|_{2}\left(1+\frac{1}{l}\right)\kappa\eta
l\left(l+1\right)\left(1+\frac{1}{l}\right)^{l-1}+(l-1)\kappa\|U_{l}\|_{2}\qquad(\text{Use
Eq. (\ref{eq:mpgnn_output_bound}), (\ref{eq:psi_value_tau_is_1})})$
$\displaystyle\leq$
$\displaystyle\|W_{l}\|_{2}\left(1+\frac{1}{l}\right)^{l}\kappa\left(\eta
l\left(l+1\right)+(l-1)\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\left(1+\frac{1}{l}\right)^{-l}\right)$
$\displaystyle\leq$
$\displaystyle\|W_{l}\|_{2}e\kappa\eta\left(l(l+1)+(l-1)\right)$
$\displaystyle=$
$\displaystyle\|W_{l}\|_{2}e\kappa\eta\left(l^{2}+2l-1\right)$
$\displaystyle\leq$
$\displaystyle\|W_{l}\|_{2}e\kappa\eta\left(l+1\right)^{2}$ (56)

Otherwise, we have,

$\displaystyle|\Delta_{l}|_{2}\leq$
$\displaystyle\|W_{l}\|_{2}\frac{\kappa\eta\tau^{l-1}}{\tau-1}\left(1+\frac{1}{l}\right)^{l}\sum_{k=1}^{l-1}\left(1-\tau^{-k}\right)+\kappa\|U_{l}\|_{2}\frac{\tau^{l-1}-1}{\tau-1}\qquad(\text{Use
Eq. (\ref{eq:mpgnn_output_bound}), (\ref{eq:psi_value_tau_is_not_1})})$
$\displaystyle\leq$
$\displaystyle\|W_{l}\|_{2}\frac{\kappa\tau^{l-1}}{\tau-1}\left(1+\frac{1}{l}\right)^{l}\left(\eta\sum_{k=1}^{l-1}\left(1-\tau^{-k}\right)+\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}(1-\tau^{1-l})\right)$
(57)

If $\tau>1$, then $\frac{1-\tau^{-k}}{\tau-1}\leq\frac{1-\tau^{1-l}}{\tau-1}$
when $1\leq k\leq l-1$. If $\tau<1$, we also have
$\frac{1-\tau^{-k}}{\tau-1}\leq\frac{1-\tau^{1-l}}{\tau-1}$ when $1\leq k\leq
l-1$. Therefore, Eq. (A.4) can be further relaxed as,

$\displaystyle|\Delta_{l}|_{2}\leq$
$\displaystyle\|W_{l}\|_{2}\frac{\kappa\tau^{l-1}}{\tau-1}\left(1+\frac{1}{l}\right)^{l}\left(\eta\sum_{k=1}^{l-1}\left(1-\tau^{-k}\right)+\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}(1-\tau^{1-l})\right)$
$\displaystyle=$
$\displaystyle\|W_{l}\|_{2}\kappa\tau^{l-1}\left(1+\frac{1}{l}\right)^{l}\left(\eta\sum_{k=1}^{l-1}\frac{1-\tau^{-k}}{\tau-1}+\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\frac{1-\tau^{1-l}}{\tau-1}\right)$
$\displaystyle\leq$
$\displaystyle\|W_{l}\|_{2}\kappa\tau^{l-1}e\left(\eta(l-1)\frac{\left(1-\tau^{1-l}\right)}{\tau-1}+\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\frac{\left(1-\tau^{1-l}\right)}{\tau-1}\right)$
$\displaystyle\leq$ $\displaystyle\|W_{l}\|_{2}\kappa\tau^{l-1}e\eta
l\frac{\left(1-\tau^{1-l}\right)}{\tau-1}\qquad\left(\text{Use
}\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\leq\eta\right)$ $\displaystyle=$
$\displaystyle e\eta l\kappa\|W_{l}\|_{2}\frac{\tau^{l-1}-1}{\tau-1},$ (58)

Therefore, combining Eq. (A.4) and Eq. (A.4), we have,

$\displaystyle|\Delta_{l}|_{2}\leq\begin{cases}e\eta\kappa\left(l+1\right)^{2}\|W_{l}\|_{2},&\text{if}\
d\mathcal{C}=1\\\ e\eta\kappa
l\|W_{l}\|_{2}\frac{\tau^{l-1}-1}{\tau-1},&\text{otherwise.}\end{cases}$ (59)

which proves the lemma. ∎

######  Theorem 3.4.

(MPGNN Generalization Bound) For any $B>0,l>1$, let
$f_{w}\in\mathcal{H}:\mathcal{X}\times\mathcal{G}\rightarrow\mathbb{R}^{K}$ be
a $l$-step MPGNN. Then for any $\delta,\gamma>0$, with probability at least
$1-\delta$ over the choice of an i.i.d. size-$m$ training set $S$ according to
$\mathcal{D}$, for any $w$, we have,

  1. 1.

If $d\mathcal{C}\neq 1$, then

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\left(\max\left(\zeta^{-(l+1)},(\lambda\xi)^{(l+1)/l}\right)\right)^{2}l^{2}h\log(lh)|w|_{2}^{2}+\log\frac{m(l+1)}{\delta}}{\gamma^{2}m}}\right).$

  2. 2.

If $d\mathcal{C}=1$, then

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\max\left(\zeta^{-6},\lambda^{3}C_{\phi}^{3}\right)(l+1)^{4}h\log(lh)|w|_{2}^{2}+\log\frac{m}{\delta}}{\gamma^{2}m}}\right).$

where $\zeta=\min\left(\|W_{1}\|_{2},\|W_{2}\|_{2},\|W_{l}\|_{2}\right)$,
$|w|_{2}^{2}=\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}$,
$\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$,
$\lambda=\|W_{1}\|_{2}\|W_{l}\|_{2}$, and
$\xi=C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1}$.

###### Proof.

We will derive the results conditioning on the value of $d\mathcal{C}$.

##### General Case $d\mathcal{C}\neq 1$

We first consider the general case $d\mathcal{C}\neq 1$. To derive the
generalization bound, we construct a special statistic of the learned weights
$\beta=\max\left(\frac{1}{\zeta},\left(\lambda\xi\right)^{\frac{1}{l}}\right)$.
It is clear that $\frac{1}{\zeta}\leq\beta$, $\lambda\xi\leq\beta^{l}$, and
$\lambda\xi/\zeta\leq\beta^{l+1}$. Note that
$\frac{1}{\zeta}=\max\left(\frac{1}{\|W_{1}\|_{2}},\frac{1}{\|W_{2}\|_{2}},\frac{1}{\|W_{l}\|_{2}}\right)$.

Consider the prior $P=\mathcal{N}(0,\sigma^{2}I)$ and the random perturbation
$u\sim\mathcal{N}(0,\sigma^{2}I)$. Note that the $\sigma$ of the prior and the
perturbation are the same and will be set according to $\beta$. More
precisely, we will set the $\sigma$ based on some approximation
$\tilde{\beta}$ of $\beta$ since the prior $P$ can not depend on any learned
weights directly. The approximation $\tilde{\beta}$ is chosen to be a cover
set which covers the meaningful range of $\beta$. For now, let us fix any
$\tilde{\beta}$ and consider $\beta$ which satisfies
$|\beta-\tilde{\beta}|\leq\frac{1}{l+1}\beta$. This also implies,

$\displaystyle|\beta-\tilde{\beta}|\leq\frac{1}{l+1}\beta~{}\Rightarrow~{}$
$\displaystyle\left(1-\frac{1}{l+1}\right)\beta\leq\tilde{\beta}\leq\left(1+\frac{1}{l+1}\right)\beta$
$\displaystyle~{}\Rightarrow~{}$
$\displaystyle\left(1-\frac{1}{l+1}\right)^{l+1}\beta^{l+1}\leq\tilde{\beta}^{l+1}\leq\left(1+\frac{1}{l+1}\right)^{l+1}\beta^{l+1}$
$\displaystyle~{}\Rightarrow~{}$
$\displaystyle\frac{1}{e}\beta^{l+1}\leq\tilde{\beta}^{l+1}\leq e\beta^{l+1}$
(60)

From Tropp (2012), for $U_{i}\in\mathbb{R}^{h\times h}$ and
$U_{i}\sim\mathcal{N}(\bm{0},\sigma^{2}I)$, we have,

$\displaystyle\mathbb{P}\left(\|U_{i}\|_{2}\geq t\right)\leq
2he^{-t^{2}/2h\sigma^{2}}.$ (61)

Taking a union bound, we have

$\displaystyle\mathbb{P}\left(\|U_{1}\|_{2}<t~{}\&~{}\cdots~{}\&~{}\|U_{l}\|_{2}<t\right)$
$\displaystyle=1-\mathbb{P}\left(\exists i,\|U_{i}\|_{2}\geq t\right)$
$\displaystyle\geq 1-\sum_{i=1}^{l}\mathbb{P}\left(\|U_{i}\|_{2}\geq t\right)$
$\displaystyle\geq 1-2lhe^{-t^{2}/2h\sigma^{2}}.$ (62)

Setting $2lhe^{-t^{2}/2h\sigma^{2}}=\frac{1}{2}$, we have
$t=\sigma\sqrt{2h\log(4lh)}$. This implies that the probability that the
spectral norm of the perturbation of any layer is no larger than
$\sigma\sqrt{2h\log(4lh)}$ holds with probability at least $\frac{1}{2}$.
Plugging this bound into Lemma 3.3, we have with probability at least
$\frac{1}{2}$,

$\displaystyle|f_{w+u}(X,A)-f_{w}(X,A)|_{2}$ $\displaystyle\leq
e\frac{t}{\zeta}lC_{\phi}B\|W_{1}\|_{2}\|W_{l}\|_{2}\frac{(d\mathcal{C})^{l-1}-1}{d\mathcal{C}-1}$
$\displaystyle=etlB\frac{\lambda\xi}{\zeta}$
$\displaystyle=eBl\beta^{l+1}t\leq
e^{2}Bl\tilde{\beta}^{l+1}\sigma\sqrt{2h\log(4lh)}\leq\frac{\gamma}{4},$ (63)

where we can set
$\sigma=\frac{\gamma}{42Bl\tilde{\beta}^{l+1}\sqrt{h\log(4lh)}}$ to get the
last inequality. Note that Lemma 3.3 also requires
$\max\left(\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}},\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}},\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\right)\leq\frac{1}{l}$.
The requirement is satisfied if $\sigma\leq\frac{\zeta}{l\sqrt{2h\log(4lh)}}$
which in turn can be satisfied if

$\displaystyle\frac{\gamma}{4eBl\beta^{l+1}\sqrt{2h\log(4lh)}}\leq\frac{1}{\beta
l\sqrt{2h\log(4lh)}},$ (64)

since the chosen value of $\sigma$ satisfies
$\sigma\leq\frac{\gamma}{4eBl\beta^{l+1}\sqrt{2h\log(4lh)}}$ and
$\frac{1}{\beta}\leq\zeta$. Therefore, one sufficient condition to make Eq.
(64) hold is $\frac{\gamma}{4eB}\leq\beta^{l}$. We will see how to satisfy
this condition later.

We now compute the KL term in the PAC-Bayes bound in Lemma 2.2.

$\displaystyle\text{KL}\left(Q\|P\right)$
$\displaystyle=\frac{|w|_{2}^{2}}{2\sigma^{2}}$
$\displaystyle=\frac{42^{2}B^{2}\tilde{\beta}^{2l+2}l^{2}h\log(4lh)}{2\gamma^{2}}\left(\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}\right)$
$\displaystyle\leq\mathcal{O}\left(\frac{B^{2}\beta^{2l+2}l^{2}h\log(lh)}{\gamma^{2}}\left(\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}\right)\right)$
(65)

From Lemma 2.2, fixing any $\tilde{\beta}$, with probability $1-\delta$ and
for all $w$ such that $|\beta-\tilde{\beta}|\leq\frac{1}{l+1}\beta$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\beta^{2l+2}l^{2}h\log(lh)|w|_{2}^{2}+\log\frac{m}{\delta}}{\gamma^{2}m}}\right).$
(66)

Finally, we need to consider multiple choices of $\tilde{\beta}$ so that for
any $\beta$, we can bound the generalization error like Eq. (66). In
particular, we only need to consider values of $\beta$ in the following range,

$\displaystyle\left(\frac{\gamma}{2B}\right)^{\frac{1}{l}}\leq\beta\leq\left(\frac{\gamma\sqrt{m}}{2B}\right)^{\frac{1}{l}},$
(67)

since otherwise the bound holds trivially as $L_{\mathcal{D},0}(f_{w})\leq 1$
by definition. To see this, if $\beta^{l}<\frac{\gamma}{2B}$, then for any
$(X,A)$ and any $j\in\mathbb{N}^{+}_{K}$, we have,

$\displaystyle\left|f_{w}(X,A)[j]\right|$
$\displaystyle\leq\left|f_{w}(X,A)\right|_{2}=|\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}|_{2}$
$\displaystyle\leq\frac{1}{n}|\bm{1}_{n}H_{l-1}|_{2}\|W_{l}\|_{2}$
$\displaystyle\leq\|W_{l}\|_{2}\max_{i}|H_{l-1}[i,:]|_{2}$ $\displaystyle\leq
BC_{\phi}\|W_{1}\|_{2}\|W_{l}\|_{2}\frac{(d\mathcal{C})^{l-1}-1}{d\mathcal{C}-1}\qquad(\text{Use
Eq. (\ref{eq:mpgnn_output_bound})})$ $\displaystyle\leq
B\lambda\xi\qquad(\text{Use definition of }\lambda\text{ and }\xi)$
$\displaystyle\leq B\beta^{l}\qquad(\text{Use definition of }\beta)$
$\displaystyle<\frac{\gamma}{2}.$ (68)

Therefore, based on the definition in Eq. (4), we always have
$L_{S,\gamma}(f_{w})=1$ when $\beta^{l}<\frac{\gamma}{2B}$. It is hence
sufficient to consider $\beta^{l}\geq\frac{\gamma}{2B}>\frac{\gamma}{4eB}$
which also makes Eq. (64) hold. Alternatively, if
$\beta^{l}>\frac{\gamma\sqrt{m}}{2B}$, the term inside the big-O notation in
Eq. (66) would be,

$\displaystyle\sqrt{\frac{B^{2}\beta^{2l}l^{2}h\log(lh)(\beta^{2}|w|_{2}^{2})+\log\frac{m}{\delta}}{\gamma^{2}m}}\geq\sqrt{\frac{l^{2}h\log(lh)(|w|_{2}^{2}/\zeta^{2})}{4}}\geq
1,$ (69)

The last inequality uses the fact that we typically choose $h\geq 2$ in
practice, $l\geq 2$ and
$|w|_{2}^{2}\geq\min\left(\|W_{1}\|_{F}^{2},\|W_{2}\|_{F}^{2},\|W_{l}\|_{F}^{2}\right)\geq\zeta^{2}$.
Since we only need to consider $\beta$ in the range of Eq. (67), one
sufficient condition to ensure $|\beta-\tilde{\beta}|\leq\frac{1}{l+1}\beta$
holds would be
$|\beta-\tilde{\beta}|\leq\frac{1}{l+1}\left(\frac{\gamma}{2B}\right)^{\frac{1}{l}}$.
Therefore, if we can find a covering of the interval in Eq. (67) with radius
$\frac{1}{l+1}\left(\frac{\gamma}{2B}\right)^{\frac{1}{l}}$ and make sure
bounds like Eq. (66) holds while $\tilde{\beta}$ takes all possible values
from the covering, then we can get a bound which holds for all $\beta$. It is
clear that we only need to consider a covering $C$ with size
$|C|=\frac{(l+1)}{2}\left(m^{1/2l}-1\right)$. Therefore, denoting the event of
Eq. (66) with $\tilde{\beta}$ taking the $i$-th value of the covering as
$E_{i}$, we have

$\displaystyle\mathbb{P}\left(E_{1}~{}\&~{}\cdots~{}\&~{}E_{|C|}\right)$
$\displaystyle=1-\mathbb{P}\left(\exists i,\bar{E}_{i}\right)\geq
1-\sum_{i=1}^{|C|}\mathbb{P}\left(\bar{E}_{i}\right)\geq 1-|C|\delta,$ (70)

where $\bar{E}_{i}$ denotes the complement of $E_{i}$. Hence, with probability
$1-\delta$ and for all $w$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})$ $\displaystyle\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\beta^{2l+2}l^{2}h\log(lh)|w|_{2}^{2}+\log\frac{m|C|}{\delta}}{\gamma^{2}m}}\right)$
$\displaystyle=L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\beta^{2l+2}l^{2}h\log(lh)|w|_{2}^{2}+\log\frac{m(l+1)}{\delta}+\frac{1}{2l}\log
m}{\gamma^{2}m}}\right)$
$\displaystyle=L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\max\left(\zeta^{-1},\left(\lambda\xi\right)^{\frac{1}{l}}\right)^{2(l+1)}l^{2}h\log(lh)|w|_{2}^{2}+\log\frac{m(l+1)}{\delta}}{\gamma^{2}m}}\right)$
(71)

which proves the theorem for the case of $d\mathcal{C}\neq 1$.

##### Special Case $d\mathcal{C}=1$

Now we consider $d\mathcal{C}=1$ of which the proof follows the logic of the
one for $d\mathcal{C}\neq 1$. Note that this case happens rarely in practice.
We only include it for the completeness of the analysis. We again construct a
statistic $\beta=\max\left(\frac{1}{\zeta},\sqrt{\lambda C_{\phi}}\right)$.
For now, let us fix any $\tilde{\beta}$ and consider $\beta$ which satisfies
$|\beta-\tilde{\beta}|\leq\frac{1}{3}\beta$. This also implies
$\frac{1}{e}\beta^{3}\leq\tilde{\beta}^{3}\leq e\beta^{3}$. Based on Lemma
3.3, we have,

$\displaystyle|f_{w+u}(X,A)-f_{w}(X,A)|_{2}$ $\displaystyle\leq
e\frac{t}{\zeta}(l+1)^{2}C_{\phi}B\|W_{1}\|_{2}\|W_{l}\|_{2}$
$\displaystyle=et(l+1)^{2}B\frac{\lambda C_{\phi}}{\zeta}\leq
eB(l+1)^{2}\beta^{3}t$ $\displaystyle\leq
e^{2}B(l+1)^{2}\tilde{\beta}^{3}\sigma\sqrt{2h\log(4lh)}\leq\frac{\gamma}{4},$
(72)

where we can set
$\sigma=\frac{\gamma}{42B(l+1)^{2}\tilde{\beta}^{3}\sqrt{h\log(4lh)}}$ to get
the last inequality. Note that Lemma 3.3 also requires
$\max\left(\frac{\|U_{1}\|_{2}}{\|W_{1}\|_{2}},\frac{\|U_{2}\|_{2}}{\|W_{2}\|_{2}},\frac{\|U_{l}\|_{2}}{\|W_{l}\|_{2}}\right)\leq\frac{1}{l}$.
The requirement is satisfied if $\sigma\leq\frac{\zeta}{l\sqrt{2h\log(4lh)}}$
which in turn can be satisfied if

$\displaystyle\frac{\gamma}{4eB(l+1)^{2}\beta^{3}\sqrt{2h\log(4lh)}}\leq\frac{1}{\beta
l\sqrt{2h\log(4lh)}},$ (73)

since the chosen value of $\sigma$ satisfies
$\sigma\leq\frac{\gamma}{4eB(l+1)^{2}\beta^{3}\sqrt{2h\log(4lh)}}$ and
$\frac{1}{\beta}\leq\zeta$. As shown later, we only need to consider a certain
range of values of $\beta$ which naturally satisfy the condition $\frac{\gamma
l}{4eB(l+1)^{2}}\leq\beta^{2}$, _i.e_., the equivalent form of Eq. (73). This
assures the applicability of Lemma 3.3. Now we compute the KL divergence,

$\displaystyle\text{KL}\left(Q\|P\right)$
$\displaystyle=\frac{|w|_{2}^{2}}{2\sigma^{2}}$
$\displaystyle=\frac{42^{2}B^{2}\tilde{\beta}^{6}(l+1)^{4}h\log(4lh)}{2\gamma^{2}}\left(\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}\right)$
$\displaystyle\leq\mathcal{O}\left(\frac{B^{2}\beta^{6}(l+1)^{4}h\log(lh)}{\gamma^{2}}\left(\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}\right)\right)$
(74)

In particular, we only need to consider values of $\beta$ in the following
range,

$\displaystyle\sqrt{\frac{\gamma}{2Bl}}\leq\beta\leq\sqrt{\frac{\gamma\sqrt{m}}{2Bl}},$
(75)

since otherwise the bound holds trivially as $L_{\mathcal{D},0}(f_{w})\leq 1$
by definition. To see this, if $\beta<\frac{\gamma}{2Bl}$, then for any
$(X,A)$ and any $j\in\mathbb{N}^{+}_{K}$, we have,

$\displaystyle\left|f_{w}(X,A)[j]\right|$
$\displaystyle\leq\left|f_{w}(X,A)\right|_{2}=|\frac{1}{n}\bm{1}_{n}H_{l-1}W_{l}|_{2}$
$\displaystyle\leq\frac{1}{n}|\bm{1}_{n}H_{l-1}|_{2}\|W_{l}\|_{2}$
$\displaystyle\leq\|W_{l}\|_{2}\max_{i}|H_{l-1}[i,:]|_{2}$ $\displaystyle\leq
B(l-1)C_{\phi}\|W_{1}\|_{2}\|W_{l}\|_{2}\qquad(\text{Use Eq.
(\ref{eq:mpgnn_output_bound})})$ $\displaystyle\leq B(l-1)\lambda
C_{\phi}\qquad(\text{Use definition of }\lambda)$ $\displaystyle\leq
Bl\beta^{2}\qquad(\text{Use definition of }\beta)$
$\displaystyle<\frac{\gamma}{2}.$ (76)

Therefore, based on the definition in Eq. (4), we always have
$L_{S,\gamma}(f_{w})=1$ when $\beta<\frac{\gamma}{2Bl}$. It it hence
sufficient to consider $\beta^{2}\geq\frac{\gamma}{2Bl}\geq\frac{\gamma
l}{4eBl^{2}}\geq\frac{\gamma l}{4eB(l+1)^{2}}$ which means the condition in
Eq. (73) is indeed satisfied. Alternatively, if
$\beta>\sqrt{\frac{\gamma\sqrt{m}}{2Bl}}$, the term inside the big-O notation
in Eq. (A.4) would be,

$\displaystyle\sqrt{\frac{B^{2}\beta^{4}(l+1)^{4}h\log(lh)\beta^{2}|w|_{2}^{2}+\log\frac{m}{\delta}}{\gamma^{2}m}}\geq\sqrt{\frac{(l+1)^{4}h\log(lh)\frac{|w|_{2}^{2}}{\zeta^{2}}}{4l^{2}}}\geq
1,$ (77)

where the first inequality hold since $\beta\geq\frac{1}{\zeta}$. The last
inequality uses the fact that we typically choose $h\geq 2$ in practice,
$l\geq 2$, and
$|w|_{2}^{2}\geq\min\left(\|W_{1}\|_{F}^{2},\|W_{2}\|_{F}^{2},\|W_{l}\|_{F}^{2}\right)\geq\zeta^{2}$.
Since we only need to consider $\beta$ in the range of Eq. (75), one
sufficient condition to ensure $|\beta-\tilde{\beta}|\leq\frac{1}{3}\beta$
always holds would be
$|\beta-\tilde{\beta}|\leq\frac{1}{3}\sqrt{\frac{\gamma}{2Bl}}$. Therefore, if
we can find a covering of the interval in Eq. (75) with radius
$\frac{1}{3}\sqrt{\frac{\gamma}{2Bl}}$ and make sure bounds like Eq. (66)
holds while $\tilde{\beta}$ takes all possible values from the covering, then
we can get a bound which holds for all $\beta$. It is clear that we only need
to consider a covering $C$ with size
$|C|=\frac{3}{2}\left(m^{\frac{1}{4}}-1\right)$.

Hence, with probability $1-\delta$ and for all $w$, we have,

$\displaystyle L_{\mathcal{D},0}(f_{w})$ $\displaystyle\leq
L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\beta^{6}(l+1)^{4}h\log(lh)|w|_{2}^{2}+\log\frac{m|C|}{\delta}}{\gamma^{2}m}}\right)$
$\displaystyle=L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\beta^{6}(l+1)^{4}h\log(lh)|w|_{2}^{2}+\log\frac{m}{\delta}+\frac{1}{4}\log{m}}{\gamma^{2}m}}\right)$
$\displaystyle=L_{S,\gamma}(f_{w})+\mathcal{O}\left(\sqrt{\frac{B^{2}\max\left(\zeta^{-6},\lambda^{3}C_{\phi}^{3}\right)(l+1)^{4}h\log(lh)|w|_{2}^{2}+\log\frac{m}{\delta}}{\gamma^{2}m}}\right)$
(78)

which proves the theorem for the case of $d\mathcal{C}=1$.

∎

######  Remark.

Note that our proof applies to both homogeneous and non-homogeneous GNNs.

###  A.5 Bound Comparison

Statistics |  | Max Node Degree
---
$d-1$
| Max Hidden Dim
---
$h$
| Spectral Norm of
---
Learned Weights
| VC-Dimension
---
(Scarselli et al., 2018)
- | $\mathcal{O}\left(h^{4}\right)$ | -
| Rademacher Complexity
---
(Garg et al., 2020) Case A
$\mathcal{O}\left(d^{l-1}\sqrt{\log(d^{l-1})}\right)$ | $\mathcal{O}\left(h\right)$ | $\mathcal{O}\left(\lambda\mathcal{C}\xi\sqrt{\log\left(\lambda\mathcal{C}\xi\right)}\right)$
| Rademacher Complexity
---
(Garg et al., 2020) Case B
$\mathcal{O}\left(d^{l-1}\sqrt{\log(d^{l-2})}\right)$ | $\mathcal{O}\left(h\sqrt{\log\sqrt{h}}\right)$ | $\mathcal{O}\left(\lambda\mathcal{C}\xi\sqrt{\log\left(\lambda\xi\right)}\right)$
| Rademacher Complexity
---
(Garg et al., 2020) Case C
$\mathcal{O}\left(d^{l-1}\sqrt{\log(d^{2l-3})}\right)$ | $\mathcal{O}\left(h\sqrt{\log\sqrt{h}}\right)$ | $\mathcal{O}\left(\lambda\mathcal{C}\xi\sqrt{\log\left(\|W_{2}\|_{2}\lambda\xi^{2}\right)}\right)$
Ours Case A | - | $\mathcal{O}\left(\sqrt{h\log h}\right)$ | $\mathcal{O}\left(\zeta^{-(l+1)}\sqrt{\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}}\right)$
Ours Case B | $\mathcal{O}\left(d^{\frac{(l+1)(l-2)}{l}}\right)$ | $\mathcal{O}\left(\sqrt{h\log h}\right)$ | $\mathcal{O}\left(\lambda^{1+\frac{1}{l}}\xi^{1+\frac{1}{l}}\sqrt{\|W_{1}\|_{F}^{2}+\|W_{2}\|_{F}^{2}+\|W_{l}\|_{F}^{2}}\right)$

Table 3: Detailed comparison of Generalization Bounds for GNNs. “-” means
inapplicable. We only consider the general case $d\mathcal{C}\|W_{2}\|_{2}\neq
1$ and simplify the Rademacher complexity based bounds (w.r.t. spectral norm
of weights) based on the assumption that $C_{\phi}\ll d\mathcal{C}\xi$ which
generally holds in practice. Here
$\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$,
$\xi=C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1}$,
$\zeta=\min\left(\|W_{1}\|_{2},\|W_{2}\|_{2},\|W_{l}\|_{2}\right)$, and
$\lambda=\|W_{1}\|_{2}\|W_{l}\|_{2}$. Note that $d^{\frac{(l+1)(l-2)}{l}}\leq
d^{\frac{l^{2}-l}{l}}=d^{l-1}$.

In this section, we explain the details of the comparison with Rademacher
complexity based generalization bounds of GNNs.

####  A.5.1 Rademacher Complexity based Bound

We first restate the Rademacher complexity bound from (Garg et al., 2020) as
below:

$\displaystyle L_{D,0}(f_{w})$ $\displaystyle\leq L_{S,\gamma}(f_{w})$
$\displaystyle+\mathcal{O}\left(\frac{1}{\gamma
m}+hB_{l}Z\sqrt{\frac{\log\left(B_{l}\sqrt{m}\max\left(Z,M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)\right)\right)}{\gamma^{2}m}}+\sqrt{\frac{\log\frac{1}{\delta}}{m}}\right)$
(79)

where
$M=C_{\phi}\frac{(C_{\phi}C_{\rho}C_{g}dB_{2})^{l-1}-1}{C_{\phi}C_{\rho}C_{g}dB_{2}-1}$,
$Z=C_{\phi}\left(BB_{1}+\bar{R}B_{2}\right)$, $\bar{R}\leq
C_{\rho}C_{g}d\min\left(b\sqrt{h},BB_{1}M\right)$, $b$ is the uniform upper
bound of $\phi$ (_i.e_., $\forall x\in\mathbb{R}^{h}$, $\phi(x)\leq b$), and
$B_{1},B_{2},B_{l}$ are the spectral norms of the weight matrices
$W_{1},W_{2},W_{l}$. Note that the numerator of $M$ has the exponent $l-1$
since we count the readout function in the number of layers/steps, _i.e_.,
there are $l-1$ message passing steps in total.

####  A.5.2 Comparison in Our Context

For typical message passing GNNs presented in the literature, node state
update function $\phi$ could be a neural network like MLP or GRU, a ReLU unit,
etc. This makes the assumption of the uniform upper bound on $\phi$
impractical, _e.g_., $b=\infty$ when $\phi$ is ReLU. Therefore, we dot not
adopt this assumption in our analysis888If we introduce the uniform upper
bound on $\phi$ in our analysis, we can also obtain a similar functional
dependency in our bound like $\min(b\sqrt{h},\cdot)$. But as aforementioned,
it is somewhat impractical and leads to a more cumbersome bound..

##### Rademacher Complexity Based Bound

Based on the above consideration, we have $\bar{R}\leq C_{\rho}C_{g}dBB_{1}M$.
We further convert some notations in the original bound to the ones in our
context.

$\displaystyle M$
$\displaystyle=C_{\phi}\frac{(C_{\phi}C_{\rho}C_{g}dB_{2})^{l-1}-1}{C_{\phi}C_{\rho}C_{g}dB_{2}-1}=\xi$
(80) $\displaystyle\bar{R}$ $\displaystyle\leq
C_{\rho}C_{g}dBB_{1}M=C_{\rho}C_{g}dB\|W_{1}\|_{2}\xi$ (81) $\displaystyle Z$
$\displaystyle=C_{\phi}\left(BB_{1}+\bar{R}B_{2}\right)=B\|W_{1}\|_{2}\left(C_{\phi}+d\mathcal{C}\xi\right),$
(82)

where we use the same abbreviations as in Theorem 3.4,
$\xi=C_{\phi}\frac{\left(d\mathcal{C}\right)^{l-1}-1}{d\mathcal{C}-1}$,
$\lambda=\|W_{1}\|_{2}\|W_{l}\|_{2}$,
$\mathcal{C}=C_{\phi}C_{\rho}C_{g}\|W_{2}\|_{2}$.

We need to consider three cases for the big-O term of the original bound in
Eq. (A.5.1) depending on the outcomes of the two point-wise maximum functions.

##### Case A

If $\max\left(Z,M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)\right)=Z$, then
the generalization bound is,

$\displaystyle\mathcal{O}$
$\displaystyle\left(hB_{l}Z\sqrt{\frac{\log\left(B_{l}\sqrt{m}Z\right)}{m}}\right)$
$\displaystyle=\mathcal{O}\left(h\|W_{l}\|B\|W_{1}\|_{2}\left(C_{\phi}+d\mathcal{C}\xi\right)\sqrt{\frac{\log\left(\|W_{l}\|_{2}\sqrt{m}B\|W_{1}\|_{2}\left(C_{\phi}+d\mathcal{C}\xi\right)\right)}{m}}\right)$
$\displaystyle=\mathcal{O}\left(hB\lambda\left(C_{\phi}+d\mathcal{C}\xi\right)\sqrt{\frac{\log\left(\sqrt{m}B\lambda\left(C_{\phi}+d\mathcal{C}\xi\right)\right)}{m}}\right).$
(83)

##### Case B

If
$\max\left(Z,M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)\right)=M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)$
and $BB_{1}=\max\left(BB_{1},\bar{R}B_{2}\right)$, then the generalization
bound is,

$\displaystyle\mathcal{O}$
$\displaystyle\left(hB_{l}Z\sqrt{\frac{\log\left(B_{l}\sqrt{m}M\sqrt{h}BB_{1}\right)}{m}}\right)$
$\displaystyle=\mathcal{O}\left(h\|W_{l}\|B\|W_{1}\|_{2}\left(C_{\phi}+d\mathcal{C}\xi\right)\sqrt{\frac{\log\left(\|W_{l}\|_{2}\sqrt{m}\xi\sqrt{h}B\|W_{1}\|_{2}\right)}{m}}\right)$
$\displaystyle=\mathcal{O}\left(hB\lambda\left(C_{\phi}+d\mathcal{C}\xi\right)\sqrt{\frac{\log\left(\sqrt{m}\lambda\xi\sqrt{h}B\right)}{m}}\right)$
(84)

##### Case C

If
$\max\left(Z,M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)\right)=M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)$
and $\bar{R}B_{2}=\max\left(BB_{1},\bar{R}B_{2}\right)$, then the
generalization bound is,

$\displaystyle\mathcal{O}$
$\displaystyle\left(hB_{l}Z\sqrt{\frac{\log\left(B_{l}\sqrt{m}M\sqrt{h}\bar{R}B_{2}\right)}{m}}\right)$
$\displaystyle=\mathcal{O}\left(h\|W_{l}\|B\|W_{1}\|_{2}\left(C_{\phi}+d\mathcal{C}\xi\right)\sqrt{\frac{\log\left(\|W_{l}\|_{2}\sqrt{m}\xi\sqrt{h}C_{\rho}C_{g}dB\|W_{1}\|_{2}\xi\|W_{2}\|_{2}\right)}{m}}\right)$
$\displaystyle=\mathcal{O}\left(hB\lambda\left(C_{\phi}+d\mathcal{C}\xi\right)\sqrt{\frac{\log\left(\lambda\sqrt{m}\sqrt{h}C_{\rho}C_{g}dB\xi^{2}\|W_{2}\|_{2}\right)}{m}}\right)$
(85)

We show the detailed dependencies of the Rademarcher complexity based bound
under three cases in Table 3. In practice, we found message passing GNNs
typically do not behave like a contraction mapping. In other words, we have
$d\mathcal{C}>1$ and $\xi\gg 1$ hold for many datasets. Therefore, the case C
happens more often in practice, _i.e_.,
$\max\left(Z,M\sqrt{h}\max\left(BB_{1},\bar{R}B_{2}\right)\right)=M\sqrt{h}\bar{R}B_{2}$.

##### PAC Bayes Bound

For our PAC-Bayes bound in Theorem 3.4, we also need to consider two cases
which correspond to
$\max\left({\zeta}^{-1},(\lambda\xi)^{\frac{1}{l}}\right)={\zeta}^{-1}$ (case
A) and
$\max\left({\zeta}^{-1},(\lambda\xi)^{\frac{1}{l}}\right)=(\lambda\xi)^{\frac{1}{l}}$
(case B) respectively. Here
$\zeta=\min\left(\|W_{1}\|_{2},\|W_{2}\|_{2},\|W_{l}\|_{2}\right)$. We show
the detailed dependencies of our bound under three cases in Table 3. Again, in
practice, we found $d\mathcal{C}>1$, $\xi\gg 1$ and $\zeta\leq 1$. Therefore,
case B occurs more often.

##### VC-dim Bound

(Scarselli et al., 2018) show that the upper bound of the VC-dimension of
general GNNs with Sigmoid or Tanh activations is $\mathcal{O}(p^{4}N^{2})$
where $p$ is the total number of parameters and $N$ is the maximum number of
nodes. Since $p=\mathcal{O}(h^{2})$ in our case, the VC-dim bound is
$\mathcal{O}(h^{8}N^{2})$. Therefore, the corresponding generalization bound
scales as $\mathcal{O}(\frac{h^{4}N}{\sqrt{m}})$. Note that $N$ is at least
$d$ and could be much larger than $d$ for some datasets.

###  A.6 Connections with Existing Bounds of MLPs/CNNs

##### ReLU Networks are Special GCNs

Since regular feedforward neural networks could be viewed as a special case of
GNNs by treating each sample as the node feature of a single-node graph, it is
natural to investigate the connections between these two classes of models. In
particular, we consider the class of ReLU networks studied in Neyshabur et al.
(2017),

$\displaystyle H_{0}$ $\displaystyle=X\quad$ $\displaystyle(\text{Input Node
Feature})$ $\displaystyle H_{k}$
$\displaystyle=\sigma_{k}\left(H_{k-1}W_{k}\right)\quad$