# Residual Connections and Normalization Can Provably Prevent Oversmoothing in
GNNs

Michael Scholkemper
Department of Computer Science
RWTH Aachen University
<EMAIL_ADDRESS>
&Xinyi Wu∗
Institute for Data, Systems, and Society
Massachusetts Institute of Technology
<EMAIL_ADDRESS>
Ali Jadbabaie
Institute for Data, Systems, and Society
Massachusetts Institute of Technology
<EMAIL_ADDRESS>
&Michael T. Schaub
Department of Computer Science
RWTH Aachen University
<EMAIL_ADDRESS>
Both authors contributed equally.

###### Abstract

Residual connections and normalization layers have become standard design
choices for graph neural networks (GNNs), and were proposed as solutions to
the mitigate the oversmoothing problem in GNNs. However, how exactly these
methods help alleviate the oversmoothing problem from a theoretical
perspective is not well understood. In this work, we provide a formal and
precise characterization of (linearized) GNNs with residual connections and
normalization layers. We establish that (a) for residual connections, the
incorporation of the initial features at each layer can prevent the signal
from becoming too smooth, and determines the subspace of possible node
representations; (b) batch normalization prevents a complete collapse of the
output embedding space to a one-dimensional subspace through the individual
rescaling of each column of the feature matrix. This results in the
convergence of node representations to the top-$k$ eigenspace of the message-
passing operator; (c) moreover, we show that the centering step of a
normalization layer — which can be understood as a projection — alters the
graph signal in message-passing in such a way that relevant information can
become harder to extract. We therefore introduce a novel, principled
normalization layer called GraphNormv2 in which the centering step is learned
such that it does not distort the original graph signal in an undesirable way.
Experimental results confirm the effectiveness of our method.

##  1 Introduction

In recent years, graph neural networks (GNNs) have gained significant
popularity due to their ability to process complex graph-structured data and
extract features in an end-to-end trainable fashion [21, 38, 5, 14, 12, 30,
43]. They have shown empirical success in a highly diverse set of problems [4,
45, 48, 34, 27, 8]. Most GNNs follow a message-passing paradigm [20], where
node representations are learned by recursively aggregating and transforming
the representations of the neighboring nodes. Through repeated message-passing
on the graph, the graph information is implicitly incorporated.

A prevalent problem in message-passing GNNs is their tendency to oversmooth
[36, 7, 46, 29], which refers to the observation that node signals (or
representations) tend to contract to a one-dimensional subspace as the number
of layers increases. While a certain amount of smoothing is desirable to
average out noise in the node features and render information extraction from
the graph more reliable [29, 47], excessive smoothing leads to an information
loss as node signals become virtually indistinguishable and can thus not be
exploited for downstream tasks. For this reason, most GNN architectures use
only few message-passing layers [30, 43, 48]. This contrasts with the trend in
deep learning, where much deeper architectures are considered preferable [24,
28]

To mitigate the oversmoothing problem for GNNs, several practical solutions
have been proposed. In particular, both residual connections [9, 31] and
normalization layers [52, 53, 22] have been specifically proposed to address
the oversmoothing problem. Despite empirical observations that the
introduction of these methods can alleviate oversmoothing and enable deeper
architectures, a theoretical analysis on the effects of these methods on
oversmoothing and the resultant expressive power of GNNs is still lacking.
More technically, the analysis of oversmoothing typically relies on a single
node similarity measure as a description of the whole underlying system [36,
6, 46]. Such measures, however, can only identify a complete collapse of the
node signals to one-dimensional subspace and are unable to capture the more
refined geometry of the system beyond a complete collapse. These observations
motivate the following questions:

> How do residual connections and normalization layers affect oversmoothing
> and therefore the practical expressive power of GNNs? How do they compare?

In this work, we answer the above questions with a refined characterization of
the underlying system of linearized GNNs [29, 47] with residual connections
and normalization layers. In particular, we establish that both methods can
alleviate oversmoothing by preventing node signals from a complete collapse to
a one-dimensional subspace, which is the case for standard GNNs [36, 6, 46].

Our contributions can be summarized as follows:

  * •

We characterize the system of linearized GNNs with residual connections. We
show that residual connections prevent complete rank collapse to a one-
dimensional subspace by incorporating the initial features at each step. As a
result, the subspace of possible node representations that a GNN can compute
is determined by the initial features.

  * •

We analyze the system of linearized GNNs with normalization layers and show
that normalization layers prevent oversmoothing of node signals through the
scaling step. Nonetheless, the node representations exponentially converge to
the top-$k$ eigenspace of the message-passing operator.

  * •

We separately identify the role of the centering step in normalization layers.
Our results suggest that the centering step can distort the message-passing in
undesired ways. Consequently, relevant graph information becomes dampened and
can even be lost in message-passing.

  * •

Based on these theoretical findings, we propose a novel normalization layer
named GraphNormv2 that learns a centering operation that does not distort the
subspaces of the message-passing operator in an uncontrolled way, but contains
a learnable projection. Experimental results show the effectiveness of our
proposed method.

##  2 Related Work

#### Oversmoothing in GNNs

Oversmoothing is a known challenge for developing deeper and more powerful
GNNs, and many techniques have been proposed to mitigate the issue in
practice. Among them, residual connections [9] and normalization layers [52,
53, 22] are popular methods that have empirically been shown to mitigate
oversmoothing to certain extent. While many theoretical works on the
underlying mechanism of oversmoothing exist [36, 7, 29, 47, 46], these studies
focus on standard GNNs without these additional modules. It is thus still an
open question how residual connections and normalization can mitigate
oversmoothing and subsequently affect the practical expressive power of GNNs
from a theoretical perspective.

#### Theoretical studies on residual connections and normalization in deep
learning

The empirical success of residual connections and normalization in enhancing
training deep neural networks has inspired research into their underlying
mechanisms [1, 11, 49, 33, 23, 18]. Specifically, it has been shown that batch
normalization avoids rank collapse for randomly initialized deep linear
networks [10] and that residual connections alleviate rank collapse in
transformers [13]. Rank collapse of node representations due to oversmoothing
has also been a notable issue in building deeper GNNs. However, a theoretical
analysis of how residual connections and normalization layers combat
oversmoothing and their additional effects on message-passing in GNNs is still
lacking.

##  3 Notation and Preliminaries

#### Notation

We use the shorthand $[n]:=\\{1,\ldots,n\\}$. We denote the zero-vector by
$\mathbf{0}$, the all-one vector by $\mathbf{1}$ and the identity matrix by
$I$. Let $\left\|\cdot\right\|_{2}$, $\left\|\cdot\right\|_{F}$ be the
$2$-norm and Frobenius norm, respectively. Lastly, for a matrix $M$, we denote
its $i^{th}$ row by $M_{i,:}$ and $j^{th}$ column by $M_{:,j}$.

#### Graph Neural Networks

Most message-passing GNN models — which we will simply refer to as GNNs from
now on — can be described by the following update equation:

$X^{(t+1)}=\sigma(AX^{(t)}W^{(t)})\,,$ (1)

where $X^{(t)},X^{(t+1)}\in\mathbb{R}^{n\times k}$ are the input and output
node representations of the $t^{\text{th}}$ layer, respectively;
$\sigma(\cdot)$ is an element-wise non-linearity such as the ReLU function;
$W^{(t)}\in\mathbb{R}^{k\times k}$ is a learnable linear transformation and
$A\in\mathbb{R}^{n\times n}_{\geq 0}$ represents a message-passing operation
and reflects the graph structure. For example, if $A$ is the graph adjacency
matrix, $A=A_{\operatorname{adj}}$ we recover the Graph Isomorphism Network
(GIN) [48], and for
$A=D^{-\frac{1}{2}}A_{\operatorname{adj}}D^{-\frac{1}{2}}$, where
$D=\text{diag}(A_{\operatorname{adj}}\mathbf{1})$ is the diagonal degree
matrix, we obtain a Graph Convolutional Network (GCN) [30].

Throughout the paper, we assume $A$ to be symmetric and non-negative.
Furthermore, we assume that eigenvalues $\lambda_{i}$ of a matrix are
organized in a non-increasing order in terms of absolute value:
$|\lambda_{1}|\geq|\lambda_{2}|\geq\cdots\geq|\lambda_{n}|$.

#### Normalization Layers

Normalization has been believed to be beneficial to deep neural networks for
nearly two decades [32]. Most normalization layers perform a centering
operation and a scaling operation on the input features. Centering usually
consists of subtracting the mean, centering the features around zero along the
chosen dimension. Scaling usually consists of scaling the features such that
the features along the chosen dimension have unit variance. The two most
popular approaches are batch and layer normalization [26, 2]. In our work, we
focus on batch normalization (BatchNorm), denoted as
$\operatorname{BN}(\cdot)$, which is defined as follows: Let
$X\in\mathbb{R}^{n\times k}$, then

$\operatorname{BN}(X)=\begin{bmatrix}\cdots,\operatorname{bn}(X_{:,i}),\cdots\end{bmatrix},$

where
$\operatorname{bn}(x)=(x-\sum_{i=1}^{n}x_{i}/n)/{\left\|x-\sum_{i=1}^{n}x_{i}/n\right\|_{2}}\,.$
We note that in the case of GNNs, especially graph classification tasks,
batches may contain nodes from different graphs. In our analysis, we consider
only the nodes in the same graph for normalization and nodes in different
graphs do not influence each other. This is sometimes called instance
normalization [42].

#### GraphNorm

In Cai et al. [7], the authors proposed GraphNorm, a normalization layer
specifically designed for GNNs that is like BatchNorm in terms of acting on
the columns of the feature matrix. Nonetheless, compared to BatchNorm, instead
of subtracting the mean, the method learns to subtract an $\tau_{j}$ portion
of the mean for the $j^{th}$ column:

$\operatorname{GraphNorm}(X_{:,j})=\gamma_{j}\frac{X_{:,j}-\tau_{j}\mathbf{1}\mathbf{1}^{\top}X_{:,j}/n}{\sigma_{j}}+\beta_{j}\,,$

where
$\sigma_{j}=\|X_{:,j}-\tau_{j}\mathbf{1}\mathbf{1}^{\top}X_{:,j}/n\|_{2}/\sqrt{n}$,
$\gamma_{j},\beta_{j}$ are the learnable affine parameters as in the
implementation of other normalization methods [2, 26]. Notably, by introducing
$\tau_{j}$ for each feature dimension $j$, [7] claims an advantage of
GraphNorm over the original batch normalization in that GraphNorm is able to
learn how much of the information to keep in the mean rather than always
subtracting the entire mean away.

#### Weisfeiler-Leman and Structural Eigenvectors

In theoretical studies about GNNs, an algorithm that comes up frequently is
the so-called Weisfeiler-Leman (WL) algorithm [44]. This algorithm iteratively
assigns a color $c(v)\in\mathbb{N}$ to each node $v$ starting from a constant
initial coloring $c^{(0)}(v)=1$ for all nodes. In each iteration, an update of
the following form is computed:
$c^{(t+1)}(v)=\text{hash}\left(c^{(t)}(v),\\{\\!\\{c^{(t)}(x)\mid
x\in\mathcal{N}(v)\\}\\!\\}\right)$, where hash is an injective hash function,
$\\{\\!\\{\cdot\\}\\!\\}$ denotes a multiset in which elements can appear more
than once, and $\mathcal{N}(v)$ is the set of neighboring nodes of $v$. The
algorithm returns the final colors $c^{(\infty)}$ when the partition
$\\{(c^{(t)})^{-1}(c^{(t)}(v))|v\\}$ no longer changes for consecutive update
$t$. For GNNs, Morris et al. [34] and Xu et al. [48] have showed that for
graphs with constant initial node features, GNNs cannot compute different
features for nodes that are in the same class in the final coloring
$c^{(\infty)}$.

For this paper, an equivalent algebraic perspective of the WL algorithm111This
algebraic characterisation is also known as the coarsest equitable partition
(cEP) of the graph. will be more useful (see Section A.1 for a detailed
discussion): Given $c^{(\infty)}$ with
$\text{Im}(c^{(\infty)})=\\{c_{1},..,c_{k}\\}$, define $H\in\\{0,1\\}^{n\times
k}$ such that $H_{v,i}=1$ if and only if $c^{(\infty)}(v)=c_{i}$. It holds
that

$AH=HA^{\pi}\,,$ (2)

where $A^{\pi}\in\mathbb{R}^{k\times k}$ is the adjacency matrix of the
_quotient graph_ , which is fixed given $A$ and $H$. In words, a node in the
quotient graph represents a class of nodes in the original graph who share the
same number of neighbors in each final color.

Most relevantly, the adjacency matrix $A$ of the original graph inherits all
the eigenpairs from the quotient graph: If $\nu^{\pi}$ is the eigenvector of
$A^{\pi}$ with eigenvalue $\lambda$, then $H\nu^{\pi}$ is an eigenvector of
$A$ with eigenvalue $\lambda$. We call such eigenvectors of $A$ the
_structural eigenvectors_. These eigenvectors are important for understanding
dynamical systems on graphs [39, 51], and play a role for centrality measures
such as PageRank [37] and others [41]. In fact, from the results of Morris et
al. [34] and Xu et al. [48], we can directly infer that given constant initial
features, GNNs effectively compute node features on this quotient graph,
meaning the features lie in the space spanned by the structural eigenvectors.

##  4 Main Results: Defying Oversmoothing

Both residual connections and normalization methods are commonly used to
mitigate oversmoothing [9, 31, 53, 52, 22]. In this section, we show that
these methods in fact _provably_ prevent oversmoothing. We further analyze
what implications the use of these methods has for the underlying system of
GNNs.

Consider the following node similarity measure [36, 46], which is a common
metric in the literature to detect oversmoothing:

$\mu(X)\vcentcolon=\left\|X-\\!\mathbf{1}\mathbf{1}^{\top}X/n\right\|_{F}\\!=\left\|(I-\mathbf{1}\mathbf{1}^{\top}/n)X\right\|_{F}\,.$
(3)

$\mu(X)$ equals zero if and only if all the columns of $X$ collapse to the
one-dimensional subspace where all nodes have the same value. In Wu et al.
[46], oversmoothing is defined to be the case where $\mu(X^{(t)})\rightarrow
0$ as $t\rightarrow\infty$ and the authors show that oversmoothing happens
exponentially for random walk GCNs and more general attention-based GNNs.
Meanwhile, Oono and Suzuki [36], Cai et al. [7] establish a similar result for
standard GCNs stating that node signals converge to the dominant eigenspace of
the message-passing operator
$D^{-\frac{1}{2}}A_{\operatorname{adj}}D^{-\frac{1}{2}}$ by using a variant of
$\mu(\cdot)$. These results suggest that for all such GNNs with non-diverging
weights, complete collapse of node signals to a one-dimensional subspace
occurs after repeated message-passing, irrespective of the initial features.
However, we will show below that these results do not hold for GNNs with
residual connections or BatchNorm.

For the following theoretical analysis, we investigate a linearized GNN,
meaning that $\sigma(\cdot)$ is the identity map. For simplicity, we assume,
if not specified otherwise, that the weights
$(W_{1}^{(t)})_{i,j},(W_{2}^{(t)})_{i,j},W^{(t)}_{i,j}$ are randomly sampled
i.i.d. Gaussians, which is typical for GNNs before training. Such a setting is
relevant in practice as oversmoothing is present in GNNs before training has
started, making the gradients used for back propagation almost vanish and
training of the network becomes much harder. Yet, most results hold under more
general conditions. The complete proofs of all the results in the main text
and the results under general conditions are provided in Appendix B.

###  4.1 Residual Connections Prevent Complete Collapse

For our analysis of residual connections, we focus on the commonly used
initial residual connections, which are deployed, e.g., in GCNII [9]. Such
residual connections are closely related to the Personalized PageRank
propagation [31]. For generality, we write the unified layer-wise update rule
for both methods as follows:

$X^{(t+1)}=(1-\alpha)AX^{(t)}W_{1}^{(t)}+\alpha X^{(0)}W_{2}^{(t)},$ (4)

where $\alpha X^{(0)}$ corresponds to the initial residual connection, and
$\alpha\in(0,1)$ can be seen as the strength of residual connections or
alternatively as the teleportation probability in the Personalized PageRank
Propagation. Note that for the Personalized PageRank Propagation method
(APPNP) proposed in [31],
$A=D^{-\frac{1}{2}}A_{\operatorname{adj}}D^{-\frac{1}{2}}$,
$W_{1}^{(t)}=I_{k}$ for all $t\geq 0$.

Intuitively, compared to the case for standard message-passing GNNs in (1), at
each update step, a linear combination of the residual signal $\alpha X^{(0)}$
is now added to the features. As long as the weight matrices are not chosen
such that they annihilate the residual signal, this will prevent the features
from collapsing to a smaller subspace. This implies that $\mu(X^{(t)})$ would
be strictly greater than zero.

######  Proposition 4.1.

If $\mu(X^{(0)})>0$, then $\mu(X^{(t)})>0$ with probability $1$.

The above result suggests that, with proper initialization of node features,
initial residual connections will alleviate oversmoothing almost surely,
meaning the features will not be smooth after iterative message-passing at
initialization. Nonetheless, it is worth noting that the node similarity
measure $\mu(X)$ can only identify a complete collapse to a one-dimensional
subspace, where all nodes share the same representation vector, in which case
$\mu(X)$ equals zero. Even if $\mu(X)$ can remain strictly positive with
residual connections, this does not eliminate the possibility that there may
still be partial collapse of the signal to a lower-dimensional subspace.
However, as we will show with the following result, with residual connections
(4), no such partial collapse occurs.

######  Proposition 4.2.

Let $x_{i}=X^{(0)}_{:,i}$ and let $\left\|x_{i}\right\|_{2}=1$ for $i\in[k]$.
Let each
$(W^{(t)}_{l})_{y,z}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{N}(\eta,s^{2})$.
Then for any $\epsilon>0$,
$\left\|x_{i}^{\top}X^{(t)}\right\|_{2}\geq\epsilon$ with probability at least
$p=1-\exp\left(-\frac{\epsilon^{2}}{2\alpha^{2}s^{2}}\right)$.

This result essentially says that a part of the initial signal is maintained
after each layer with high probability. We can further give a more refined
characterization of what node features the system in (4) can compute after
repeated rounds of message-passing.

######  Proposition 4.3.

Let
$\operatorname{Kr}(A,X^{(0)})=\operatorname{Span}(\\{A^{i-1}X^{(0)}_{:,j}\\}_{i\in[n],j\in[k]})$
be the Krylov subspace. Let $Y\in\mathbb{R}^{n\times k}$. Then there exist a
$T\in\mathbb{N}$ and sequence of weights
$W_{1}^{(0)},W_{2}^{(0)},\cdots,W_{1}^{(T-1)},W_{2}^{(T-1)}$ such that
$X^{(T)}=Y$ if and only if $Y\in\operatorname{Kr}(A,X^{(0)})$.

The result shows that for a message-passing GNN with residual connections, the
subspace that the embedding $X^{(t)}$ lies in is governed by the initial
features $X^{(0)}$ together with the message-passing operator $A$ and that any
such signal is reachable by a sequence of weights. This is in contrast with
the behavior of standard message-passing GNNs, for which node representations
eventually becomes “memoryless”, i.e., independent of initial features. In
particular, for standard GCNs, the subspace the system converges to is
completely governed by the message-passing operator $A$ [36].

######  Remark 4.4.

Our results also suggest that whether initial residual connections would work
for oversmoothing heavily depends on the initialization of features. If chosen
poorly, they may not be able to prevent oversmoothing. In particular, a
constant initialization of features is unhelpful to prevent oversmoothing with
initial residual connections.

###  4.2 BatchNorm Prevents Complete Collapse

Having discussed the effect of residual connections, in this section, we
switch gears and analyze how BatchNorm affects GNNs. We consider the following
combination of GNNs with BatchNorm:

$X^{(t+1)}=\operatorname{BN}(Y^{(t+1)}),\quad
Y^{(t+1)}=\sigma(AX^{(t)}W^{(t)})\,.$ (5)

Similar to the analysis on GNNs with initial residual connections, we will
show that BatchNorm prevents complete collapse of the output embedding space
to the all-ones subspace. We again consider the case where $\sigma(\cdot)$ is
the identity map. Note that even though the GNN computation is now linearized,
the overall system remains non-linear because of the BatchNorm operation. Let
us consider now $V_{\neq 0}\in\mathbb{R}^{n\times k^{\prime}}$ to be the
matrix of eigenvectors associated to non-zero eigenvalues of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$.

######  Proposition 4.5.

Let $\operatorname{Rank}(V_{\neq 0}^{\top}X^{(0)})>1$, then
$\mu(X^{(t)})\geq\sqrt{2}$ for all $t\geq 1$ with probability $1$ .

The above result suggests that with BatchNorm, $\mu(X^{(t)})$ can be
maintained strictly greater than a nontrivial constant at each layer,
indicating that complete collapse to the all-ones subspace does not happen.
Notably, a generalized result also holds for the case where $W^{(t)}$ is the
identity matrix for all $t\geq 0$, meaning that the corresponding system does
not oversmooth in terms of $\mu(\cdot)$ with BatchNorm (see Section B.6). This
again, stands in contrast to the case in standard GNNs without BatchNorm, in
which having $W^{(t)}=I_{k}$ has been proven to oversmooth in GCNs and more
general attention-based GNNs [36, 46].

######  Remark 4.6.

A more general version of the above result holds for both linear and non-
linear GNNs, when assuming $\sigma(\cdot)$ is injective:
$\mu(X^{(t)})=\sqrt{k}$ for all $t\geq 1$, where $k$ is the hidden dimension
of features. See Section B.7 for a detailed discussion. If nonlinearity
$\sigma(\cdot)$ is applied after BatchNorm, then it also holds that
$\mu(X^{(t)})>0$ for all $t\geq 1$ under proper assumptions. See Section B.8
for a detailed discussion.

However, as what we have discussed for the case of residual connections, the
node similarity measure $\mu(X)$ can only capture complete collapse to a one-
dimensional subspace. In what follows, we will provide a more precise
characterization of the convergence behaviors of GNNs with BatchNorm. Notably,
since the scaling operation of BatchNorm guarantees that the system will not
diverge or collapse, we can give an exact asymptotic characterization. Let
$V_{k}\in\mathbb{R}^{n\times k}$ be the matrix of the top-$k$ eigenvectors of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$. We will show that the resulting
linearized GNN converges to a rank-$k$ subspace spanned by the top-$k$
eigenvectors of $(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$.

######  Proposition 4.7.

Suppose $V_{k}^{\top}X^{(0)}$ has rank $k$, then for all weights $W^{(t)}$,
the GNN with BatchNorm given in (5) exponentially converges to the column
space of $V_{k}$.

In the proof of Proposition 4.7, it becomes clear that the centering operation
in BatchNorm is the reason for which we use the _centered_ message-passing
operator $(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$. If we left out the
centering step and only keep the scaling step of BatchNorm, the proof would
work in the same way (switching all occurrences of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$ with $A$). This implies that the
column-wise scaling operation is responsible for the preservation of the rank
of the features. Furthermore, there are no requirements for the weights as in
the case without BatchNorm [36, 46]: even extremely large or random weights
can be chosen, as the scaling ensures that the system neither diverges nor
collapses.

Moreover, the convergence of the linearized system in (5) to a $k$-dimensional
subspace can be shown to be tight: we can choose weights such that the top-$k$
eigenvectors are exactly recovered. This is of course only possible if these
eigenvectors do have a nonzero eigenvalue, which we assume for the result
below, in that, none of the top-$k$ eigenvectors of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$ has eigenvalue $\hat{\lambda}_{i}$
zero.

######  Proposition 4.8.

Suppose $|\hat{\lambda}_{k}|>0$ and $V_{k}^{\top}X^{(0)}$ has rank $k$. For
any $\epsilon>0$, there exists $T>0$ and a sequence of weights
$W^{(0)},W^{(1)},...,W^{(T)}$ such that for all $t\geq T$ and $i\in[k]$,

$\left\|\nu_{i}^{\top}X^{(t)}_{:,i}\right\|_{2}\geq 1/\sqrt{1+\epsilon}\,,$

where $\nu_{i}$ denotes the $i$-th eigenvector of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$.

In fact, the above result ties in nicely with recent results showing that GNNs
are strengthened through _positional encodings_ [15, 16], where the features
are augmented by the top-$k$ eigenvectors of the graph. This can, in some
sense, be seen as emulating a deep GNN, which would converge to the top-$k$
eigenspace using BatchNorm. Yet it is worth noting that while BatchNorm
improves the practical expressive power of GNNs by converging to a larger
subspace, under the type of convergence given in Proposition 4.7, the
information in the eigenvectors associated with small magnitude eigenvalues is
still dampened after repeated message-passing.

######  Remark 4.9.

If no initial node features $X^{(0)}$ are given by the dataset, common choices
in practice are to initialize the features randomly or identically for each
node. In the former case, the prerequisites of 4.7 and Proposition 4.8 are
satisfied. In the latter case, the conditions are not met, as $X^{(0)}$ has
rank one. In that case, the system still converges. In fact, it retains its
rank and converges to the dominant eigenvector of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$.

###  4.3 Comparison between normalization and residual connections

From previous sections, we have seen that both residual connections and batch
normalization are able to prevent a complete collapse of the node embeddings
to a one-dimensional subspace. In both cases, the embeddings converge to a
larger subspace and thus oversmoothing is alleviated. However, it is clear
that different mechanisms are at play to mitigate oversmoothing. With residual
connections, the system is able to keep the dimensions of the initial input
features by incorporating the initial features $X^{(0)}$ at each layer; while
with normalization, the system converges to the subspace of the top
eigenvectors of the message-passing operator through the scaling step.

##  5 Centering Distorts the Graph

So far, we have analyzed the effects of residual connections and normalization
layers on oversmoothing. Specifically, we have shown that the incorporation of
initial features of residual connections and the scaling effect of
normalization help alleviate complete rank collapse of node features. However,
there are two steps in normalization layers: centering and scaling. If already
scaling helps preventing a complete collapse, a natural question is _what is
the role of centering in the process?_ In this section, we will show that the
current centering operation used in normalization layers can in fact have an
undesirable effect altering the graph signal that message-passing can extract,
as if message-passing happens on a different graph.

###  5.1 Centering Interferes the Structural Eigenvectors

The centering operation in normalization layers takes away the (scaled) mean
across all rows in each column, and thus can be written as applying the
operator $I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n$ to the input, where $\tau$
indicates how much mean is taken away in centering. Specifically, for $\tau=1$
we recover BatchNorm, whereas for $\tau\in\mathbb{R}$, we recover GraphNorm
[7]. To analyze how this step would alter the graph signal message-pasing can
extract, we make use of the concepts of quotient graph and structural
eigenvectors as introduced in Section 3.

Given a symmetric, non-negative adjacency matrix $A\in\mathbb{R}^{n\times
n}_{\geq 0}$, let $H\in\\{0,1\\}^{n\times m}$ be the indicator matrix of its
final WL coloring $c^{(\infty)}$. Consider the eigenpairs
$\mathcal{V}=\\{...,(\lambda,\nu),...\\}$ of $A$ and divide them into the set
of structural eigenpairs
$\mathcal{V}_{\text{struc}}=\\{(\lambda,\nu)\in\mathcal{V}\mid\nu=H\nu^{\pi}\\}$,
and the remaining eigenpairs
$\mathcal{V}_{\text{rest}}=\mathcal{V}\backslash\mathcal{V}_{\text{struc}}$.
Similarly, let $\hat{\mathcal{V}}=\\{...,(\hat{\lambda},\hat{\nu}),...\\}$ be
the eigenpairs of $(I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)A$. We now analyze
what happens to these distinct sets of eigenvectors when applying the
centering operation. Notice that the parameter $\tau$ controls how much of the
mean is taken away and thus how much the centering influences the input graph.
However, as long as $\tau$ is not zero, there is always an effect altering the
graph operator used in message-passing:

######  Proposition 5.1.

Assuming $\tau\neq 0$,

  1. 1.

$\mathcal{V}_{\text{rest}}\subset\hat{\mathcal{V}}$.

  2. 2.

Assume that $A$ is not regular, then the dominant eigenvector $\nu$ of $A$ is
not an eigenvector in $\hat{\mathcal{V}}$ for any eigenvalue.

  3. 3.

$\sum_{(\lambda,\nu)\in\mathcal{V}_{\text{struc}}}\lambda>\sum_{(\hat{\lambda},\hat{\nu})\in\hat{V}\backslash\mathcal{V}_{\text{rest}}}\hat{\lambda}$
.

Here, we consider a general case where the graph is not regular, meaning there
is more than one color in the final WL coloring $c^{(\infty)}$. Assuming
constant initialization of node features, $\mathcal{V}_{\text{struc}}$ spans
the space of all possible node features that a GNN can compute and Proposition
5.1 states that it is exactly the space that the centering acts on.
Specifically, while leaving the eigenvectors and eigenvalues in
$\mathcal{V}_{\text{rest}}$ untouched, centering changes the eigenvector basis
of the space spanned by $\mathcal{V}_{\text{struc}}$ in two ways: some
vectors, such as the dominant eigenvector, are affected by this transformation
and thus no longer convey the same information. At the same time, the
centering transformation may change the magnitude of eigenvalues — that is,
the dominant eigenvector may not be dominant anymore. Meanwhile, the whole
space is pushed downward in the spectrum, meaning that after the centering
transformation, the signal components within the structural eigenvectors are
dampened and thus become less pronounced in the node representations given by
the GNN.

Notably, such an effect altering the graph signal not only applies to
BatchNorm with $\tau=1$, but also GraphNorm [7], which was proposed
specifically for GNNs. In their paper, the authors address the problem that
BatchNorm’s centering operation completely nullifies the graph signal on
regular graphs. Their remedy is to only subtract an $\tau$ portion of the mean
instead. However, Proposition 5.1 shows that a similar underlying problem
altering the graph signal would persist for general graphs even switching from
BatchNorm to GraphNorm.

#### Comparison with standard neural networks

The use of normalization techniques in GNNs was inspired by use of
normalization methods in standard feed-forward neural networks [26]. Here, we
want to emphasize that the issue centering causes in GNNs as described above
is not a problem for standard neural networks, as standard neural networks do
not need to incorporate the graph information in the forward pass. As a
result, in standard neural networks, normalization transforms the input in a
way that it does not affect the classification performance: for each neural
network, there exists a neural network that yields an equivalent
classification after normalization. However, this is not the case for GNNs. In
GNNs, the graph information is added in through message-passing, which can be
heavily altered by normalization as shown above. Such normalization can lead
to information loss, negatively impacting the model’s performance.

###  5.2 Our Method: GraphNormv2

Based on our theoretical analysis, we propose a new normalization layer for
GNNs which has a similar motivation as the original GraphNorm but improves the
centering operation to not affect the graph information in message-passing.
Specifically, instead of naive centering, which can be thought of as
subtracting a projection to the all-ones space, we use a learned projection.
Learning a completely general projection may have certain downsides, however.
As graphs can have different sizes, we either would need to learn different
projections for different graph sizes or use only part of the learned
projection for smaller graphs. We thus opt for learning a centering that
transforms the features by subtracting an $(\tau_{j})_{i}$ portion of the
$i$-th eigenvector from the $j$-th column. Our proposed graph normalization is
thus:

$\operatorname{\text{GraphNormv2}}(X_{:,j})=\gamma_{j}\frac{X_{:,j}-(V_{k+}\tau_{j}\tau_{j}^{\top}V_{k+}^{\top})X_{:,j}}{\sigma_{j}}+\beta_{j}\,,$

where $\tau_{j}\in\mathbb{R}^{k+1}$ is a learnable parameter,
$\sigma_{j}=\|X_{:,j}-(V_{k+}\tau_{j}\tau_{j}^{\top}V_{k+}^{\top})X_{:,j}\|_{2}$,
and $\gamma_{j},\beta_{j}$ are the learnable affine parameters. Instead of
just using the top-$k$ eigenvectors $V_{k}$ of the message-passing operator,
we use $V_{k}$ and one additional vector
$r=\mathbf{1}-V_{k}V_{k}^{\dagger}\mathbf{1}$ that guarantee one to retrieve
the all-ones vector $\mathbf{1}$. We denote the set
$V_{k+}\vcentcolon=V_{k}\cup\\{r\\}$. This ensures backward compatibility, in
that, GraphNormv2 can emulate GraphNorm and BatchNorm.

##  6 Numerical Experiments

Figure 1: Long Term behavior of GCN. Progression of $\mu(X^{(t)})$ and
$\operatorname{Rank}(X^{(t)})$ over $256$ iterations of message-passing in
both linear and non-linear GCN. In the linear case, $\mu(X^{(t)})$ remains
constant for all methods, indicating that complete collapse to the all-ones
space does not happen. However, PairNorm and the non-normalized variant do
collapse in terms of rank, while the other methods maintain a rank greater
than $2$. All the phenomena are explained by our theory. In the non-linear
case, the models behave similarly. Notably, centering seems to prevent rank
collapse in the non-linear case as PairNorm no longer collapses in rank.

In this section, we investigate the benefits that can be derived from the
proposed graph normalization method GraphNormv2. We examine the long term
behavior of linear and non-linear GNNs by conducting an ablation study on
randomly initialized, untrained GNNs. We then go on to inspect the practical
relevance of our proposed method. More details about the experiments are
provided in Appendix C.

Table 1: Performance under different normalization layers. Performance of GIN,
GCN, and GAT with different normalization layers on graph classification task
(MUTAG, PROTEINS, and PTC-MR) and node classification task (Cora, CiteSeer,
and ogbn-arxiv). Results are reported as the mean accuracy (in %) $\pm$ std.
over 10 independent trials and 5 folds. Best results are highlighted in blue;
second best results are highlighted in gray.

c c| c c c | c c c | c c [colortbl-like] Graph Classification Node Classification Node Classification ($\\#$ layers=20)
MUTAG PROTEINS PTC-MR Cora CiteSeer ogbn-arxiv Cora CiteSeer
GCN no norm $78.2\pm 7.8$$70.5\pm 4.0$ $\pagecolor{gray!20}57.7\pm 3.0$
$82.6\pm 4.6$ $89.4\pm 0.6$ $67.4\pm 0.2$ $46.0\pm 11.5$ $65.6\pm 4.3$
batch $\pagecolor{gray!20}81.5\pm 2.0$$70.4\pm 1.8$$56.0\pm 3.5$ $85.2\pm 0.6$
$89.4\pm 0.7$ $\pagecolor{best}{68.3\pm 0.3}$ $\pagecolor{gray!20}84.0\pm 1.6$
$82.2\pm 4.4$
graph $81.1\pm 4.7$$\pagecolor{gray!20}71.4\pm 3.6$$\pagecolor{best}{58.1\pm
4.8}$ $85.1\pm 0.8$ $89.1\pm 0.7$ $68.3\pm 0.2$ $80.5\pm 6.3$
$\pagecolor{gray!20}82.7\pm 2.1$
pair $61.3\pm 10.1$ $59.6\pm 0.0$ $55.8\pm 0.0$ $84.2\pm 0.7$ $88.2\pm 0.7$
$63.0\pm 1.2$ $82.3\pm 0.5$ $74.0\pm 14.1$
graphv2$\pagecolor{best}{82.6\pm 4.6}$$\pagecolor{best}{72.6\pm 2.6}$$56.5\pm
4.8$ $\pagecolor{best}{85.8\pm 0.5}$ $\pagecolor{best}{89.5\pm 0.6}$
$\pagecolor{best}{68.3\pm 0.3}$ $\pagecolor{best}{84.8\pm 0.8}$
$\pagecolor{best}{83.4\pm 2.0}$
GAT no norm $78.5\pm 5.3$$\pagecolor{gray!20}71.2\pm 2.0$
$\pagecolor{best}{61.3\pm 4.2}$ $85.9\pm 0.9$ $93.7\pm 0.5$ $68.3\pm 0.9$
$79.2\pm 0.7$ $89.1\pm 0.5$
batch $\pagecolor{best}{82.5\pm 4.1}$$68.4\pm 5.2$ $55.7\pm 0.7$ $86.2\pm 0.7$
$\pagecolor{best}{94.9\pm 0.5}$ $\pagecolor{best}{71.0\pm 0.1}$ $85.8\pm 0.7$
$\pagecolor{gray!20}94.2\pm 0.6$
graph $81.2\pm 4.8$$\pagecolor{best}{71.9\pm 2.4}$$\pagecolor{gray!20}60.3\pm
5.4$ $\pagecolor{gray!20}86.3\pm 0.9$ $\pagecolor{gray!20}94.8\pm 0.5$
$\pagecolor{best}{71.0\pm 0.1}$ $\pagecolor{best}{85.7\pm 0.8}$ $94.0\pm 0.6$
pair $59.7\pm 11.1$ $60.1\pm 0.9$ $55.7\pm 1.1$ $85.4\pm 1.0$ $93.4\pm 1.2$
$70.0\pm 0.2$ $84.9\pm 0.7$ $92.5\pm 1.0$
graphv2$\pagecolor{gray!20}81.6\pm 3.8$$71.0\pm 2.6$$59.5\pm 4.4$
$\pagecolor{best}{86.3\pm 0.8}$ $\pagecolor{gray!20}94.8\pm 0.5$ $70.9\pm 0.2$
$\pagecolor{best}{85.7\pm 0.8}$ $\pagecolor{best}{94.4\pm 0.5}$
GIN no norm $79.7\pm 5.9$$70.7\pm 3.6$$\pagecolor{gray!20}59.2\pm 3.9$
$33.4\pm 45.4$ $47.0\pm 3.1$ $\pagecolor{gray!20}21.6\pm 0.0$ $28.2\pm 6.2$
$23.4\pm 5.1$
batch $82.5\pm 5.0$$69.2\pm 5.3$$52.9\pm 5.3$ $\pagecolor{gray!20}68.4\pm 4.0$
$49.7\pm 6.8$ $21.1\pm 7.9$ $\pagecolor{gray!20}31.4\pm 5.4$ $23.6\pm 6.1$
graph $\pagecolor{gray!20}83.7\pm 4.2$$\pagecolor{best}{72.6\pm 2.4}$$59.1\pm
5.1$ $33.2\pm 3.4$ $\pagecolor{gray!20}63.0\pm 11.0$
$\pagecolor{gray!20}21.6\pm 0.0$ $26.9\pm 5.2$ $\pagecolor{gray!20}27.1\pm
4.6$
pair $65.2\pm 3.2$ $64.5\pm 4.3$ $55.5\pm 1.6$ $37.2\pm 3.5$ $51.1\pm 9.2$
$20.7\pm 9.6$ $28.6\pm 6.2$ $26.9\pm 4.7$
graphv2$\pagecolor{best}{84.9\pm 3.6}$$\pagecolor{gray!20}71.0\pm
3.6$$\pagecolor{best}{60.1\pm 5.5}$ $\pagecolor{best}{86.5\pm 0.9}$
$\pagecolor{best}{94.9\pm 0.5}$ $\pagecolor{best}{68.9\pm 0.3}$
$\pagecolor{best}{83.8\pm 2.0}$ $\pagecolor{best}{93.8\pm 0.7}$

#### Rank collapse in linear and non-linear GNNs

We investigate the effect of normalization in deep (linear) GNNs on the Cora
dataset [50]. We employ seven different architectures: an architecture using
residual connections, BatchNorm [26], BatchNorm without the centering
operation, PairNorm [52], GraphNorm [7], GraphNormv2, and finally no
normalization as a baseline. Each of these methods is used in an untrained,
randomly initialized GCN [30] (without biases) with 256 layers.

We compare these models using the following measures of convergence:
$\mu(X^{(t)})$ as defined in (3), the eigenvector space projection:
$d_{\text{ev}}(X)\vcentcolon=\frac{1}{n}\left\|X-VV^{\top}X\right\|_{F},$
where $V\in\mathbb{R}^{n\times n}$ is the set of normalized eigenvectors of
$A$, and the numerical rank of the features $\operatorname{Rank}(X^{(t)})$.
The results are shown in Figure 1. The same experiment with both GIN and GAT
can be found in the appendix together with additional measures of convergence.

The two left panels of Figure 1 show that the commonly considered metric for
measuring oversmoothing $\mu(X)$ indeed does not detect any collapse of the
feature space — neither in the linear nor non-linear case. However, the right
two panels show that there is more going on as PairNorm and no norm do
collapse to a rank one subspace in the linear case. Specifically, they
converge to the dominant eigenvector of
$(I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)D^{-\frac{1}{2}}AD^{-\frac{1}{2}})$
and $D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$, respectively, as evident from the
middle panel. The other methods are able to preserve a rank greater than two
over all iterations. However, they do converge to a low-dimensional subspace,
e.g. GraphNormv2 converges to the top-$k$ eigenvector space as can be seen in
the middle panel. All of these phenomena are explained by our theoretical
analysis. As for the non-linear case, the models behave similarly apart from
two things: (a) there is no convergence to the linear subspace due to the non-
linearity $\sigma(\cdot)$, although the rank can still be preserved. (b) The
centering operation seems to prevent rank collapse in the non-linear case as
PairNorm no longer collapses in rank.

#### Classification performance

We evaluate the effectiveness of our method GraphNormv2for real graph learning
tasks. We perform graph classification tasks on the standard benchmark
datasets MUTAG [40], PROTEINS [35] and PTC-MR [3] as well as node
classification tasks on Cora, Citeseer [50] and large-scale ogbn-arxiv [25].
Following the general set-up of [17], we investigate the performance of GIN,
GCN and GAT in a 5-fold cross-validation setting. Details on hyperparameter
tuning and other specifics can be found in Appendix C. The final test scores
are obtained as the mean scores across the 5 folds and 10 independent trials
with the selected hyperparameters. We then also repeat the same experiment on
Cora and Citeseer where we fix the depth of the models to $20$ layers. The
results are reported in Section 6.

Section 6 shows improvements on most benchmarks for our proposed normalization
technique GraphNormv2. Although our method yields weak performance
improvements in certain cases, this trend is apparent across datasets and
tasks and even seems to be independent of the architecture. It is however
worth mentioning that our method seems not to perform as well with GAT. This
may be the result of GAT having both asymmetric and time-varying message-
passing operators due to the attention mechanism [46] — both aspects are
outside the scope of our current analysis.

##  7 Discussion

In this paper, we have analyzed the effect of both residual connections and
normalization layers in GNNs. We show that both methods provably alleviate
oversmoothing through the incorporation of the initial features and the
scaling operation, respectively. In addition, by identifying that the
centering operation in a normalization layer alters the graph information in
message-passing, we proposed GraphNormv2, a novel normalization layer which
does not distort the graph and empirically verified its effectiveness.

Our experiments showed that the trends described by our theoretical analysis
are visible even in the non-linear case. Future work may concern itself with
closing the gap and explaining how centering together with non-linearity can
prevent node representations from collapsing to a low-dimensional subspace.

## References

  * Arora et al. [2018] S. Arora, Z. Li, and K. Lyu.  Theoretical analysis of auto rate-tuning by batch normalization.  In _ICLR_ , 2018.
  * Ba et al. [2016] J. Ba, J. R. Kiros, and G. E. Hinton.  Layer normalization.  _ArXiv_ , abs/1607.06450, 2016.
  * Bai et al. [2019] Y. Bai, H. Ding, Y. Qiao, A. Marinovic, K. Gu, T. Chen, Y. Sun, and W. Wang.  Unsupervised inductive graph-level representation learning via graph-graph proximity.  In _IJCAI_ , 2019.
  * Battaglia et al. [2016] P. Battaglia, R. Pascanu, M. Lai, D. Jimenez Rezende, and k. kavukcuoglu.  Interaction networks for learning about objects, relations and physics.  In _NeurIPS_ , 2016.
  * Bruna et al. [2014] J. Bruna, W. Zaremba, A. D. Szlam, and Y. LeCun.  Spectral networks and locally connected networks on graphs.  In _ICLR_ , 2014.
  * Cai and Wang [2020] C. Cai and Y. Wang.  A note on over-smoothing for graph neural networks.  In _ICML Graph Representation Learning and Beyond (GRL+) Workshop_ , 2020.
  * Cai et al. [2021] T. Cai, S. Luo, K. Xu, D. He, T.-y. Liu, and L. Wang.  Graphnorm: A principled approach to accelerating graph neural network training.  In _ICML_ , 2021.
  * Cappart et al. [2021] Q. Cappart, D. Chételat, E. B. Khalil, A. Lodi, C. Morris, and P. Velickovic.  Combinatorial optimization and reasoning with graph neural networks.  _Journal of Machine Learning Research_ , 2021.
  * Chen et al. [2020] M. Chen, Z. Wei, Z. Huang, B. Ding, and Y. Li.  Simple and deep graph convolutional networks.  In _ICML_ , 2020.
  * Daneshmand et al. [2020] H. Daneshmand, J. M. Kohler, F. R. Bach, T. Hofmann, and A. Lucchi.  Batch normalization provably avoids ranks collapse for randomly initialised deep networks.  In _NeurIPS_ , 2020.
  * Daneshmand et al. [2021] H. Daneshmand, A. Joudaki, and F. R. Bach.  Batch normalization orthogonalizes representations in deep random networks.  In _NeurIPS_ , 2021.
  * Defferrard et al. [2016] M. Defferrard, X. Bresson, and P. Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  In _NeurIPS_ , 2016.
  * Dong et al. [2021] Y. Dong, J.-B. Cordonnier, and A. Loukas.  Attention is not all you need: Pure attention loses rank doubly exponentially with depth.  In _ICML_ , 2021.
  * Duvenaud et al. [2015] D. K. Duvenaud, D. Maclaurin, J. Aguilera-Iparraguirre, R. Gómez-Bombarelli, T. D. Hirzel, A. Aspuru-Guzik, and R. P. Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In _NeurIPS_ , 2015.
  * Dwivedi et al. [2021] V. P. Dwivedi, A. T. Luu, T. Laurent, Y. Bengio, and X. Bresson.  Graph neural networks with learnable structural and positional representations.  _arXiv preprint arXiv:2110.07875_ , 2021.
  * Dwivedi et al. [2023] V. P. Dwivedi, C. K. Joshi, A. T. Luu, T. Laurent, Y. Bengio, and X. Bresson.  Benchmarking graph neural networks.  _Journal of Machine Learning Research_ , 24(43):1–48, 2023.
  * Errica et al. [2019] F. Errica, M. Podda, D. Bacciu, and A. Micheli.  A fair comparison of graph neural networks for graph classification.  _arXiv preprint arXiv:1912.09893_ , 2019.
  * Esteve et al. [2020] C. Esteve, B. Geshkovski, D. Pighin, and E. Zuazua.  Large-time asymptotics in deep learning.  _ArXiv_ , abs/2008.02491, 2020.
  * Fey and Lenssen [2019] M. Fey and J. E. Lenssen.  Fast graph representation learning with PyTorch Geometric.  In _ICLR Workshop on Representation Learning on Graphs and Manifolds_ , 2019.
  * Gilmer et al. [2017] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl.  Neural message passing for quantum chemistry.  In _ICML_ , 2017.
  * Gori et al. [2005] M. Gori, G. Monfardini, and F. Scarselli.  A new model for learning in graph domains.  In _IJCNN_ , 2005.
  * Guo et al. [2023] X. Guo, Y. Wang, T. Du, and Y. Wang.  Contranorm: A contrastive learning perspective on oversmoothing and beyond.  In _ICLR_ , 2023.
  * Hardt and Ma [2017] M. Hardt and T. Ma.  Identity matters in deep learning.  In _ICLR_ , 2017.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _CVPR_ , 2016.
  * Hu et al. [2020] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  _arXiv preprint arXiv:2005.00687_ , 2020.
  * Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _ICML_ , 2015.
  * Jha et al. [2022] K. Jha, S. Saha, and H. Singh.  Prediction of protein–protein interaction using graph neural networks.  _Scientific Reports_ , 2022.
  * Kaplan et al. [2020] J. Kaplan, S. McCandlish, T. J. Henighan, T. B. Brown, B. Chess, R. Child, S. Gray, A. Radford, J. Wu, and D. Amodei.  Scaling laws for neural language models.  _ArXiv_ , abs/2001.08361, 2020.
  * Keriven [2022] N. Keriven.  Not too little, not too much: a theoretical analysis of graph (over)smoothing.  In _NeurIPS_ , 2022.
  * Kipf and Welling [2017] T. N. Kipf and M. Welling.  Semi-supervised classification with graph convolutional networks.  In _ICLR_ , 2017.
  * Klicpera et al. [2019] J. Klicpera, A. Bojchevski, and S. Günnemann.  Predict then propagate: Graph neural networks meet personalized pagerank.  In _ICLR_ , 2019.
  * LeCun et al. [2002] Y. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller.  Efficient backprop.  In _Neural networks: Tricks of the trade_ , pages 9–50. Springer, 2002.
  * Liu et al. [2019] T. Liu, M. Chen, M. Zhou, S. S. Du, E. Zhou, and T. Zhao.  Towards understanding the importance of shortcut connections in residual networks.  In _NeurIPS_ , 2019.
  * Morris et al. [2019] C. Morris, M. Ritzert, M. Fey, W. L. Hamilton, J. E. Lenssen, G. Rattan, and M. Grohe.  Weisfeiler and leman go neural: Higher-order graph neural networks.  In _AAAI_ , 2019.
  * Morris et al. [2020] C. Morris, N. M. Kriege, F. Bause, K. Kersting, P. Mutzel, and M. Neumann.  Tudataset: A collection of benchmark datasets for learning with graphs.  In _ICML Graph Representation Learning and Beyond (GRL+) Workshop_ , 2020.
  * Oono and Suzuki [2020] K. Oono and T. Suzuki.  Graph neural networks exponentially lose expressive power for node classification.  In _ICLR_ , 2020.
  * Sánchez-García [2020] R. J. Sánchez-García.  Exploiting symmetry in network analysis.  _Communications Physics_ , 2020.
  * Scarselli et al. [2009] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini.  The graph neural network model.  _IEEE Transactions on Neural Networks_ , 20:61–80, 2009\.
  * Schaub et al. [2016] M. T. Schaub, N. O’Clery, Y. N. Billeh, J.-C. Delvenne, R. Lambiotte, and M. Barahona.  Graph partitions and cluster synchronization in networks of oscillators.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 26(9), 2016.
  * Schlichtkrull et al. [2017] M. Schlichtkrull, T. Kipf, P. Bloem, R. van den Berg, I. Titov, and M. Welling.  Modeling relational data with graph convolutional networks.  In _Extended Semantic Web Conference_ , 2017.
  * Stamm et al. [2023] F. I. Stamm, M. Scholkemper, M. Strohmaier, and M. T. Schaub.  Neighborhood structure configuration models.  In _WWW_ , 2023.
  * Ulyanov et al. [2016] D. Ulyanov, A. Vedaldi, and V. S. Lempitsky.  Instance normalization: The missing ingredient for fast stylization.  _ArXiv_ , abs/1607.08022, 2016.
  * Veličković et al. [2018] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio.  Graph attention networks.  In _ICLR_ , 2018.
  * Weisfeiler and Leman [1968] B. Weisfeiler and A. Leman.  The reduction of a graph to canonical form and the algebra which appears therein.  _NTI, Series_ , 2(9):12–16, 1968.
  * Wu et al. [2020] S. Wu, W. Zhang, F. Sun, and B. Cui.  Graph neural networks in recommender systems: A survey.  _ACM Computing Surveys_ , 55:1 – 37, 2020.
  * Wu et al. [2023a] X. Wu, A. Ajorlou, Z. Wu, and A. Jadbabaie.  Demystifying oversmoothing in attention-based graph neural networks.  In _NeurIPS_ , 2023a.
  * Wu et al. [2023b] X. Wu, Z. Chen, W. Wang, and A. Jadbabaie.  A non-asymptotic analysis of oversmoothing in graph neural networks.  In _ICLR_ , 2023b.
  * Xu et al. [2019] K. Xu, W. Hu, J. Leskovec, and S. Jegelka.  How powerful are graph neural networks?  In _ICLR_ , 2019.
  * Yang et al. [2019] G. Yang, J. Pennington, V. Rao, J. N. Sohl-Dickstein, and S. S. Schoenholz.  A mean field theory of batch normalization.  In _ICLR_ , 2019.
  * Yang et al. [2016] Z. Yang, W. W. Cohen, and R. Salakhutdinov.  Revisiting semi-supervised learning with graph embeddings.  In _ICML_ , 2016.
  * Yuan et al. [2013] Y. Yuan, G.-B. Stan, L. Shi, M. Barahona, and J. Goncalves.  Decentralised minimum-time consensus.  _Automatica_ , 2013.
  * Zhao and Akoglu [2020] L. Zhao and L. Akoglu.  Pairnorm: Tackling oversmoothing in gnns.  In _ICLR_ , 2020.
  * Zhou et al. [2021] K. Zhou, Y. Dong, K. Wang, W. S. Lee, B. Hooi, H. Xu, and J. Feng.  Understanding and resolving performance degradation in deep graph convolutional networks.  In _CIKM_ , 2021.

##  Appendix A Further Background Materials

###  A.1 An algebraic perspective on the Weisfeiler-Leman coloring

We briefly restate part of what has already been stated in the main text. So
that notation is clear again. The Weisfeiler-Leman (WL) algorithm iteratively
assigns a color $c(v)\in\mathbb{N}$ to each node $v\in V$ starting from a
constant initial coloring $c^{(0)}(v)=1\ \ \forall v\in V$. In each iteration,
an update of the following form is computed:

$c^{(t+1)}(v)=\text{hash}\left(c^{(t)}(v),\\{\\!\\{c^{(t)}(x)\mid x\in
N(v)\\}\\!\\}\right)$ (6)

where hash is an injective hash function, and $\\{\\!\\{\cdot\\}\\!\\}$
denotes a multiset in which elements can appear more than once. The algorithm
returns the final colors $c^{(\infty)}$ when the partition
$\\{(c^{(t)})^{-1}(c^{(t)}(v))|v\in V\\}$ no longer changes for consecutive
$t$. Morris et al. [34] and Xu et al. [48] showed that GNNs cannot compute
different features for nodes that are in the same class in the final coloring
$c^{(\infty)}$.

For this paper, the equivalent algebraic perspective of the WL algorithm will
be more useful: Given $c^{(\infty)}$ with
$\text{Im}(c^{(\infty)})=\\{c_{1},..,c_{k}\\}$, define $H\in\\{0,1\\}^{n\times
k}$ such that $H_{v,i}=1$ if and only if $c^{(\infty)}(v)=c_{i}$. It holds
that

$AH=H(H^{\top}H)^{-1}H^{\top}AH=HA^{\pi}\,,$ (7)

where $A^{\pi}:=(H^{\top}H)^{-1}H^{\top}AH\in\mathbb{R}^{k\times k}$ is the
adjacency matrix of the quotient graph. Looking at this equation, there are
two things of interest here.

Firstly, $AH=HA^{\pi}$. Considering this equation more closely, the left side
$AH$ counts, for each node, the number of neighbors of color $c_{i}$. More
formally,

$(AH)_{v,i}=\sum_{x\in(c^{(\infty)})^{-1}(c_{i})}[x\in N(v)]=\sum_{x\in
N(v)}[c(x)=c_{i}]$

where the Iverson bracket $[\cdot]$ returns $1$ if the statement is satisfied
and $0$ if it is not. The right hand side of the equation ($HA^{\pi}$) states
that nodes in the same class have the same rows. It is not hard to verify that
if $c^{(\infty)}(v)=c_{i}$

$(HA^{\pi})_{v,:}=(A^{\pi})_{i,:}$

Now, combining both sides of the equation, the statement $AH=HA^{\pi}$ reads,
that the number of neighbors of any color that a node has is the same for all
nodes of the same class meaning $\\{\\!\\{c^{(t)}(x)\mid x\in N(v)\\}\\!\\}$
is the same for all nodes of the class. It becomes clear that at this point we
have a fixed point of the WL update equation Equation 6 and conversely
whenever we have such a fixed point, the equation $AH=HA^{\pi}$ holds. The
difference between the WL algorithm and Equation 7 is that the latter equation
holds for any so called _equitable partition_ while the WL algorithm converges
to the _coarsest_ partition - meaning the partition with fewest distinct
colors. For instance, on regular graphs, the WL algorithm returns the
partition with $1$ color i.e. $H=\mathbb{1}$. However, the trivial partition
with $n$ colors $H=I$ also fulfills Equation 7.

Secondly $A^{\pi}:=(H^{\top}H)^{-1}H^{\top}AH\in\mathbb{R}^{k\times k}$ is the
adjacency matrix of the so-called _quotient graph_. Noticing that
$H^{\top}H=\operatorname{diag}(\\{|(c^{(\infty)})^{-1}(c_{i})|\\}_{i\in[k]})$
it quickly becomes clear that $A^{\pi}$ is the graph of mean connectivity
between colors. In other words, a supernode in the quotient graph represent a
class of nodes in the original graph who share the same number of neighbors in
each final coloring and an edge connecting two such supernodes is weighted by
the number of edges there were going from nodes of the one color to nodes of
the other color. In this sense, the quotient graph is a compression of the
graph structure. It now only has $k$ supernodes compared to the $n$ nodes that
there were in the original graph. Still the quotient graph is an adequate
depiction of the structure of the graph. Most relevantly, the adjacency matrix
$A$ of the original graph inherits all eigenpairs from the quotient graph: Let
$(\lambda,\nu^{\pi})$ be an eigenpair of $A^{\pi}$, then

$AH\nu^{\pi}=HA^{\pi}\nu^{\pi}=\lambda H\nu^{\pi}\,.$

We call such eigenvectors of $A$ the structural eigenvectors. They are
profoundly important and may even completely determine processes that move
over the edges of a network. The structural eigenvectors span a linear
subspace that is invariant to multiplication with $A$. This means that once
inside this subspace a graph dynamical system cannot leave it. This holds true
for GNNs as Morris et al. [34] and Xu et al. [48] showed.

##  Appendix B Proofs

###  B.1 Proof for Proposition 4.1

See 4.1

###### Proof.

If $\mu(X^{(0)})>0$, then, there must be a column in $X^{(0)}$ that is not a
scaled version of the all-ones vector - otherwise $\mu(X^{(0)})=0$ follows
directly. Take this column $X^{(0)}_{:,i}\neq c\mathbb{1}$. Then there exists
indices $j,l$ s.t. $X^{(0)}_{j,i}\neq X^{(0)}_{l,i}$. We now prove by
induction that with probability $1$, $X^{(t)}_{j,i}\neq X^{(t)}_{l,i}$ \-
which will directly yield the statement.

The base case for the initial features is obviously true.

Consider w.l.o.g the first column in the next iteration:

$\displaystyle X^{(t+1)}_{:,1}$
$\displaystyle=(1-\alpha)AX^{(t)}(W^{(t)}_{1})_{:,1}+\alpha
X^{(0)}(W^{(t)}_{2})_{:,1}$
$\displaystyle=(1-\alpha)\sum_{q=1}^{k}(AX^{(t)})_{:,q}(W^{(t)}_{1})_{q,1}+\alpha\sum_{q=1}^{k}X^{(0)}_{:,q}(W^{(t)}_{2})_{q,1}$
$\displaystyle=\sum_{q=1}^{k}(1-\alpha)(AX^{(t)})_{:,q}(W^{(t)}_{1})_{q,1}+\alpha
X^{(0)}_{:,q}(W^{(t)}_{2})_{q,1}$

Now, looking at the entries $j,l$, which were distinct for column $i$ in the
initial features, we get that:

$\displaystyle X^{(t+1)}_{j,1}$
$\displaystyle=\sum_{q=1}^{k}(1-\alpha)(AX^{(t)})_{j,q}(W^{(t)}_{1})_{q,1}+\alpha
X^{(0)}_{j,q}(W^{(t)}_{2})_{q,1}$ $\displaystyle=\sum_{q=1,q\neq
i}^{k}(1-\alpha)(AX^{(t)})_{j,q}(W^{(t)}_{1})_{q,1}+\alpha
X^{(0)}_{j,q}(W^{(t)}_{2})_{q,1}$
$\displaystyle\quad\quad+(1-\alpha)(AX^{(t)})_{j,i}(W^{(t)}_{1})_{i,1}+\alpha
X^{(0)}_{j,i}(W^{(t)}_{2})_{i,1}$ $\displaystyle=\varphi^{(t)}_{j}+\alpha
X^{(0)}_{j,i}(W^{(t)}_{2})_{i,1}$

By the same reasoning, $X^{(t+1)}_{l,1}=\varphi^{(t)}_{l}+\alpha
X^{(0)}_{l,i}(W^{(t)}_{2})_{i,1}$. Note that both $\varphi^{(t)}_{j}$ and
$\varphi^{(t)}_{l}$ are not influenced by $(W^{(t)}_{2})_{i,1}$. To wrap up
the proof, we consider the probability that the opposite of what we want to
show happens, and show that this has probability $0$. Consider:

$\displaystyle\mathcal{A}$
$\displaystyle=\\{W^{(t)}_{1},W^{(t)}_{2}\in\mathbb{R}^{k\times k}\mid
X^{(t+1)}_{j,1}=X^{(t+1)}_{l,1}\\}$
$\displaystyle=\\{W^{(t)}_{1},W^{(t)}_{2}\in\mathbb{R}^{k\times
k}\mid\varphi^{(t)}_{j}+\alpha
X^{(0)}_{j,i}(W^{(t)}_{2})_{i,1}-\varphi^{(t)}_{l}-\alpha
X^{(0)}_{l,i}(W^{(t)}_{2})_{i,1}=0\\}$

$\mathcal{A}$ defines a proper hyperplane in the space of the randomly sampled
weights. As such it has Lebesgue measure $0$ \- which in turn means it has
probability $0$. Thus, the probability of the opposite event

$\mathcal{B}=\\{W^{(t)}_{1},W^{(t)}_{2}\in\mathbb{R}^{k\times k}\mid
X^{(t)}_{j,i}\neq
X^{(t)}_{l,i}\\}=\\{W^{(t)}_{1},W^{(t)}_{2}\in\mathbb{R}^{k\times
k}\mid\mu(X^{(t)})>0\\}$

is $1$. This concludes the induction and the proof. ∎

###  B.2 Proposition 4.1: deterministic case

For the deterministic version, we adopt the following regularity conditions on
weight matrices:

######  Assumption B.1.

For the system described in (4), assume: there exists $\epsilon>0$ such that

  1. 1.

$\sum_{m=0}^{t}\alpha(1-\alpha)^{m}\lambda_{i}^{m}W_{2}^{(t-m)}W_{1}^{(t-m+1)}\cdots
W_{1}^{(t)}+(1-\alpha)^{t+1}\lambda_{i}^{t+1}W_{1}^{(0)}\cdots W_{1}^{(t)}$
has smallest singular value $\sigma_{\min}\geq\epsilon$ for all $i\in[n]$.

  2. 2.

$\sum_{m=0}^{t}\alpha(1-\alpha)^{m}\lambda_{i}^{m}W_{2}^{(t-m)}W_{1}^{(t-m+1)}\cdots
W_{1}^{(t)}+(1-\alpha)^{t+1}\lambda_{i}^{t+1}W_{1}^{(0)}\cdots W_{1}^{(t)}$
converges as $t\to\infty$ and has smallest singular value
$\sigma_{\min}\geq\epsilon$ for all $i\in[n]$.

Suppose $A$ is full-rank. If weights $W_{1}^{(t)},W_{2}^{(t)}$ are orthogonal,
then B.1.1 holds. On the other hand, B.1.2 is an asymptotic technical
condition to ensure that weights are non-collapsing and non-diverging in the
limit. Some ways to satisfy the assumptions is to have the spectral radius of
$A$, $\rho(A)\leq 1$ and $W_{1}^{(t)},W_{2}^{(t)}=I_{k}$ for any $t\geq 0$.

We restate Proposition 4.1 with full conditions: let $V\in\mathbb{R}^{n\times
n}$ be the matrix of eigenvectors for $A$.

######  Proposition B.2.

Under B.1, let $\nu_{q}\in V$ be such that $v_{q}^{\top}\mathbf{1}=0$. If
$X^{(0)}$ is properly initialized, such as if $X^{(0)}$ is not the zero matrix
and $\|\nu_{q}^{\top}X^{(0)}\|_{2}=c>0$, then $\mu(X^{(t)})\geq
c\epsilon/\sqrt{k}$ for all $t\geq 0$ and $\lim_{t\to\infty}\mu(X^{(t)})\geq
c\epsilon/\sqrt{k}$.

###### Proof.

Writing (4) recursively, we get that

$\displaystyle X^{(t+1)}$
$\displaystyle=\alpha\sum_{m=0}^{t}(1-\alpha)^{m}A^{m}X^{(0)}W_{2}^{(t-m)}W_{1}^{(t-m+1)}...W_{1}^{(t)}$
$\displaystyle+(1-\alpha)^{t+1}A^{t+1}X^{(0)}W_{1}^{(0)}...W_{1}^{(t)}\,.$

For each column $X^{(t+1)}_{:,i}$, similarly, one can prove by induction that

$X^{(t+1)}_{:,i}=\sum_{l=1}^{n}\sum_{m=0}^{t}\sigma_{l,i}^{(t,m)}\lambda_{l}^{m}\nu_{l}\,,$

where

  * •

$\Sigma^{(0,0)}=V^{\top}X^{(0)}$,

  * •

$\Sigma^{(t,0)}=\alpha\Sigma^{(0,0)}W_{2}^{(t-1)}$ for all $t\geq 0$,

  * •

$\Sigma^{(t,m)}=\alpha(1-\alpha)^{m}\Sigma^{(0,0)}W_{2}^{(t-m-1)}W_{1}^{(t-m)}...W_{1}^{(t-1)}$,
for all $1\leq m\leq t-1$,

  * •

$\Sigma^{(t,t)}=(1-\alpha)^{t}\Sigma^{(0,0)}W_{1}^{(0)}...W_{1}^{(t-1)}$.

Then $\nu_{q}^{\top}X^{(t)}$

$=\nu_{q}^{\top}X^{(0)}\left(\sum_{m=0}^{t}\alpha(1-\alpha)^{m}\lambda_{i}^{m}W_{2}^{(t-m)}W_{1}^{(t-m+1)}\cdots
W_{1}^{(t)}+(1-\alpha)^{t+1}\lambda_{i}^{t+1}W_{1}^{(0)}\cdots
W_{1}^{(t)}\right).$ (8)

Since by construction, $\|\nu_{q}^{\top}X^{(0)}\|_{2}=c$, it follows from the
regularity conditions on weights that

$\|\nu_{q}^{\top}X^{(t)}\|_{2}\geq c\epsilon\,.$

This implies that

$\|\nu_{q}^{\top}X^{(t)}\|_{\infty}\geq c\epsilon/\sqrt{k}\,,$

which means that there exists $i\in[k]$ such that

$\left|\nu_{q}^{\top}X_{:,i}^{(t)}\right|=c\epsilon/\sqrt{k}\,.$

Note that since $\nu_{q}\in V$ and $\nu_{q}^{\top}\mathbf{1}=0$, we get that

$\displaystyle\mu(X^{(t)})$
$\displaystyle=\left\|X^{(t)}-\frac{\mathbf{1}\mathbf{1}^{\top}}{n}X^{(t)}\right\|_{F}=\sqrt{\sum_{l=1}^{k}\left\|X_{:,l}^{(t)}-\frac{\mathbf{1}\mathbf{1}^{\top}}{n}X_{:,l}^{(t)}\right\|^{2}_{2}}$
$\displaystyle\geq\sqrt{\left\|X_{:,i}^{(t)}-\frac{\mathbf{1}\mathbf{1}^{\top}}{n}X_{:,i}^{(t)}\right\|^{2}_{2}}$
$\displaystyle\geq\left|\nu_{q}^{\top}X_{:,i}^{(t)}\right|$

which means that $\mu(X^{(t)})\geq c\epsilon/\sqrt{k}$.

Similarly, we can show that $\lim_{t\to\infty}\mu(X^{(t)})\geq
c\epsilon/\sqrt{k}$. ∎

###  B.3 Proof for Proposition 4.2

See 4.2

###### Proof.

We start by deconstructing $X^{(t)}$ as

$X^{(t)}=(1-\alpha)AX^{(t-1)}W_{1}^{(t-1)}+\alpha X^{(0)}W_{2}^{(t-1)}$

This means that:

$\displaystyle x_{i}^{\top}X^{(t)}$
$\displaystyle=(1-\alpha)x_{i}^{\top}AX^{(t-1)}W_{1}^{(t-1)}+\alpha
x_{i}^{\top}X^{(0)}W_{2}^{(t-1)}$ (9) $\displaystyle=\varphi+\alpha
x_{i}^{\top}X^{(0)}W_{2}^{(t-1)}$

Resulting in:

$\displaystyle\left\|x_{i}^{\top}X^{(t)}\right\|_{2}\geq\left\|x_{i}^{\top}X^{(t)}\right\|_{\infty}$
$\displaystyle=\left\|\varphi+\alpha
x_{i}^{\top}X^{(0)}W_{2}^{(t-1)}\right\|_{\infty}$ (10)
$\displaystyle=\max_{j}|\varphi_{j}+\alpha(x_{i}^{\top}X^{(0)}W_{2}^{(t-1)})_{j}|$
$\displaystyle\geq|\varphi_{j}+\alpha(x_{i}^{\top}X^{(0)}W_{2}^{(t-1)})_{j}|$
$\displaystyle=|\varphi_{j}+\alpha x_{i}^{\top}(X^{(0)}W_{2}^{(t-1)})_{:,j}|$
$\displaystyle=|\varphi_{j}+\alpha x_{i}^{\top}X^{(0)}(W_{2}^{(t-1)})_{:,j}|$
$\displaystyle=|\varphi_{j}+\alpha\sum_{a}(x_{i}^{\top})_{a}(X^{(0)})_{a,:}(W_{2}^{(t-1)})_{:,j}|$
$\displaystyle=|\varphi_{j}+\alpha\sum_{b}\sum_{a}(x_{i}^{\top})_{a}(X^{(0)})_{a,b}(W_{2}^{(t-1)})_{b,j}|=|\hat{Z}|$

$\hat{Z}$ is a weighted sum of Gaussian variables and as such, is Gaussian
itself ($\mathcal{N}(\hat{\eta},\hat{s}^{2})$) with mean
$\hat{\eta}=\eta(\varphi_{j})+\eta\sum_{b}\sum_{a}(x_{i}^{\top})_{a}(X^{(0)})_{a,b}$
and variance
$\hat{s}^{2}=s^{2}(\varphi_{j})+\alpha^{2}s^{2}\sum_{b}(\sum_{a}(x_{i}^{\top})_{a}(X^{(0)})_{a,b})^{2}$.
Because $x_{i}$ is normalized, we have that
$\sum_{a}(x_{i}^{\top})_{a}(X^{(0)})_{a,i}=1$ and as such
$\hat{s}^{2}\geq\alpha^{2}s^{2}$.

Define the helper variable $\hat{Z}/\alpha$, which is Gaussian with
$\hat{Z}/\alpha\sim\mathcal{N}(\frac{\hat{\eta}}{\alpha},\frac{\hat{s^{2}}}{\alpha^{2}})$
and define the variable $Z\sim\mathcal{N}(0,s^{2})$. Notice that
$\hat{Z}/\alpha$ has higher variance than $Z$. To finish the proof, notice
that

$Pr(|\hat{Z}|\geq\alpha\epsilon)=Pr(|\hat{Z}|/\alpha\geq\epsilon)\geq
Pr(|Z|\geq\epsilon)=1-Pr(|Z|<\epsilon)\geq
1-\exp\left(-\frac{\epsilon^{2}}{2s^{2}}\right)\,,$

where the last inequality is based on the Chernoff Bound. ∎

###  B.4 Proposition 4.2: deterministic case

We complement 4.2 with the following result. Let $V\in\mathbb{R}^{n\times n}$
be the matrix of eigenvectors of $A$ and define

$V^{\star}\vcentcolon=\\{\nu_{q}\in
V:\left\|\nu_{q}^{\top}X^{(0)}\right\|_{2}=c_{q}\\}$
$V^{0}\vcentcolon=\\{\nu_{p}\in
V:\left\|\nu_{p}^{\top}X^{(0)}\right\|_{2}=0\\}\,,$

where $c_{q}>0$ for all $q$. In words, $V^{\star}$ is the set of eigenspaces
of $A$ onto which $X^{(0)}$ has a non-trivial projection, and $V^{0}$ is the
set of eigenspaces of $A$ onto which $X^{(0)}$ has no projection.

######  Proposition B.3.

Under B.1.1, for all $t\geq 1$,

$\left\|\nu_{q}^{\top}X^{(t)}\right\|_{2}=c_{q}\epsilon,\forall\nu_{q}\in
V^{\star},\quad\left\|\nu_{p}^{\top}X^{(t)}\right\|_{2}=0,\forall\nu_{p}\in
V^{0}.$

Under B.1.2,

$\lim_{t\to\infty}\left\|\nu_{q}^{\top}X^{(t)}\right\|_{2}>c_{q}\epsilon,\forall\nu_{q}\in
V^{\star}\,,\quad\lim_{t\to\infty}\left\|\nu_{p}^{\top}X^{(t)}\right\|_{2}=0,\forall\nu_{p}\in
V^{0}\,.$

###### Proof.

The proof follows directly from the form (8). ∎

The above result states that the signal excited in the original graph input
$X^{(0)}$ is precisely what stays and the signal that is not excited in
$X^{(0)}$ can never be created through message-passing.

We give the following concrete example of the above result:

#### Example

Suppose $W_{1}^{(t)},W_{2}^{(t)}=I$, then

$X^{(t+1)}=\left(\alpha\sum_{k=0}^{t}(1-\alpha)^{k}A^{k}+(1-\alpha)^{t+1}A^{t+1}\right)X\,.$

Note that when $\rho(A)<1/(1-\alpha)$ such as
$A=D^{-1/2}A_{\text{adj}}D^{-1/2}$, as $t\to\infty$,

$\underset{t\to\infty}{\lim}X^{(t)}=\alpha(I_{n}-(1-\alpha)A)^{-1}X\,.$

Let $(\lambda_{i},\nu_{i})$ be the $i$-th eigenpair of $A$ and
$\sigma_{l,i}=\langle v_{l},X_{:,i}\rangle$ , then

$X^{(t+1)}_{:,i}=\sum_{l=1}^{n}\left(\sum_{k=1}^{t}\alpha(1-\alpha)^{k}\lambda_{l}^{k}+(1-\alpha)^{t+1}\lambda_{l}^{t+1}\right)\sigma_{l,i}\nu_{l}\,,$

and when $\rho(A)<1/(1-\alpha)$ such as $A=D^{-1/2}A_{\text{adj}}D^{-1/2}$,

$\underset{t\to\infty}{\lim}X^{(t+1)}_{:,i}=\sum_{l=1}^{n}\frac{\alpha}{1-(1-\alpha)\lambda_{l}}\sigma_{l,i}\nu_{l}\,.$

This implies that for all $\nu_{q}$,

$\displaystyle\|\nu_{q}^{\top}X^{(t)}\|_{2}$
$\displaystyle=\sqrt{\sum_{i=1}^{k}\left(\left(\sum_{m=1}^{t-1}\alpha(1-\alpha)^{m}\lambda_{q}^{m}+(1-\alpha)^{t}\lambda_{q}^{t}\right)\sigma_{q,i}\right)^{2}}$
$\displaystyle=\left(\sum_{m=1}^{t-1}\alpha(1-\alpha)^{m}\lambda_{q}^{m}+(1-\alpha)^{t}\lambda_{q}^{t}\right)\|\sigma_{q,:}\|_{2}$

and as $t\to\infty$,

$\lim_{t\to\infty}\|\nu_{q}^{\top}X^{(t)}\|_{2}=\frac{\alpha}{1-(1-\alpha)\lambda_{q}}\|\sigma_{q,:}\|_{2}\,.$

###  B.5 Proof of Proposition 4.3

See 4.3

###### Proof.

Set $T=n$, $\alpha=0.5$ and $W^{(t)}_{1}=I$ for $t>0$ and
$W^{(0)}_{1}=\mathbb{0}$. Begin by unrolling the recursive equation Equation
4:

$X^{(1)}=(1-\alpha)AX^{(0)}W_{1}^{(0)}+\alpha X^{(0)}W^{(0)}_{2}$

And in turn:

$\displaystyle X^{(2)}$ $\displaystyle=(1-\alpha)AX^{(1)}W_{1}^{(1)}+\alpha
X^{(0)}W^{(1)}_{2}$
$\displaystyle=(1-\alpha)A((1-\alpha)AX^{(0)}W_{1}^{(0)}+\alpha
X^{(0)}W^{(0)}_{2})I+\alpha X^{(0)}W^{(1)}_{2}$
$\displaystyle=0.5^{2}AX^{(0)}W^{(0)}_{2}+0.5X^{(0)}W^{(1)}_{2}$

Iterating this, yields:

$X^{(n)}=\sum_{i=1}^{n}0.5^{i}A^{i-1}X^{(0)}W^{(i-1)}_{2}$

Now consider a single column of $X^{(n)}$:

$\displaystyle(X^{(n)})_{:,j}$
$\displaystyle=\sum_{i=1}^{n}0.5^{i}A^{i-1}X^{(0)}(W^{(i-1)}_{2})_{:,j}$
$\displaystyle=\sum_{l=1}^{k}\sum_{i=1}^{n}0.5^{i}(W^{(i-1)}_{2})_{l,j}A^{i-1}X^{(0)}_{:,l}$

Now let $Y_{:,j}\in\operatorname{Kr}(A,X^{(0)})$ be in the Krylov subspace.
Then $Y_{:,j}=\sum_{l=1}^{k}\sum_{i=1}^{n}w_{l,i}A^{i-1}X^{(0)}_{:,l}$.
Setting $(W_{2}^{(i-1)})_{l,j}=\frac{w_{l,i}}{0.5^{i}}$ yields
$X^{(n)}_{:,j}=Y_{:,j}$.

For the other direction, begin similiarly by unrolling the recursive equation:

$\displaystyle X^{(2)}$ $\displaystyle=(1-\alpha)AX^{(1)}W_{1}^{(1)}+\alpha
X^{(0)}W^{(1)}_{2}$
$\displaystyle=(1-\alpha)A((1-\alpha)AX^{(0)}W_{1}^{(0)}+\alpha
X^{(0)}W^{(0)}_{2})W_{1}^{(1)}+\alpha X^{(0)}W^{(1)}_{2}$
$\displaystyle=(1-\alpha)^{2}A^{2}X^{(0)}W_{1}^{(0)}W_{1}^{(1)}+(1-\alpha)\alpha
AX^{(0)}W^{(0)}_{2}W_{1}^{(1)}+\alpha X^{(0)}W^{(1)}_{2}$

Iterating this, yields:

$X^{(n)}=\sum_{i=1}^{n+1}A^{i-1}X^{(0)}\mathcal{W}^{(i-1)}$

Now consider a single column of $X^{(n)}$:

$\displaystyle(X^{(n)})_{:,j}$
$\displaystyle=\sum_{l=1}^{k}\sum_{i=1}^{n+1}A^{i-1}X^{(0)}_{:,l}\mathcal{W}^{(i-1)}_{l,j}$

Now, setting $w_{l,i}=\mathcal{W}^{(i-1)}_{l,j}$ and
$Y_{:,j}=\sum_{l=1}^{k}\sum_{i=1}^{n+1}w_{l,i}A^{i-1}X^{(0)}_{:,l}$ verifies
that $X^{(n)}_{:,j}=Y_{:,j}\in\operatorname{Kr}(A,X^{(0)})$. ∎

###  B.6 Proof for Proposition 4.5

See 4.5

###### Proof.

We prove this by induction. The base case for $0$ holds. Assume
$\operatorname{Rank}(V_{\neq 0}^{\top}X^{(t)})$ has rank at least $2$. Then,
there exist at least $2$ columns $X^{(t)}_{:,i}$ and $X^{(t)}_{:,j}$ such that
$\operatorname{Rank}(V_{\neq
0}^{\top}\begin{bmatrix}X^{(t)}_{:,i},X^{(t)}_{:,j}\end{bmatrix})=2$. Consider
their eigenvector decomposition in terms of eigenvectors of
$(I-\mathbf{1}\mathbf{1}^{\top}/n)A$:

$\displaystyle X^{(t)}_{:,i}=\sum_{l=1}^{n}\sigma^{(t)}_{l,i}v_{l},\quad
X^{(t)}_{:,j}=\sum_{l=1}^{n}\sigma^{(t)}_{l,j}v_{l}\,.$

Consider the action of $(I-\mathbf{1}\mathbf{1}^{\top}/n)A$:

$\displaystyle\tilde{X}^{(t)}_{:,i}=(I-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)}_{:,i}=\sum_{l=1}^{n}\lambda_{l}\sigma^{(t)}_{l,i}v_{l},\quad\tilde{X}^{(t)}_{:,j}=(I-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)}_{:,j}=\sum_{l=1}^{n}\lambda_{l}\sigma^{(t)}_{l,j}v_{l}\,.$

Since $X^{(t)}_{:,i}$ and $X^{(t)}_{:,j}$ are linearly independent and the
eigenvectors of a symmetric matrix are orthogonal, there exists $q$ such that
$\sigma^{(t)}_{q,i}\neq\sigma^{(t)}_{q,j}$ with $\lambda_{q}\neq 0$. This
exists because $\operatorname{Rank}(V_{\neq
0}^{\top}\begin{bmatrix}X^{(t)}_{i,:},X^{(t)}_{j,:}\end{bmatrix})=2$. It
follows that $\sigma^{(t)}_{q,i}\lambda_{q}\neq\sigma^{(t)}_{q,j}\lambda_{q}$,
and thus the centered features
$\tilde{X}^{(t)}_{:,i}\neq\tilde{X}^{(t)}_{:,j}$ and neither
$\tilde{X}^{(t)}_{:,i}=\mathbb{0}$ nor $\tilde{X}^{(t)}_{:,j}=\mathbb{0}$
(otherwise $X^{(t)}_{:,i}$ would be a $0$ eigenvector and be orthogonal to
$V_{\neq 0}$). Thus, they are linearly independent. Furthermore,
$\tilde{X}^{(t)}_{:,i},\tilde{X}^{(t)}_{:,j},\mathbf{1}$ are linearly
independent, since $(I-\mathbf{1}\mathbf{1}^{\top}/n)$ projects to the space
orthogonal to $\operatorname{Span}\\{\mathbf{1}\\}$. Write:

$\displaystyle X^{(t+1)}_{:,i}$
$\displaystyle=\frac{1}{\Gamma_{i}}(I-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)}W^{(t)}_{:,i}$
$\displaystyle=\frac{1}{\Gamma_{i}}\sum_{a=1}^{k}W^{(t)}_{a,i}\tilde{X}^{(t)}_{:,a}$
$\displaystyle=\frac{1}{\Gamma_{i}}\sum_{a=1,a\neq i,a\neq
j}^{k}W^{(t)}_{a,i}\tilde{X}^{(t)}_{:,a}+W^{(t)}_{i,i}\tilde{X}^{(t)}_{:,i}+W^{(t)}_{j,i}\tilde{X}^{(t)}_{:,j}$
$\displaystyle=\frac{1}{\Gamma_{i}}(\varphi^{(i)}+W^{(t)}_{i,i}\tilde{X}^{(t)}_{:,i}+W^{(t)}_{j,i}\tilde{X}^{(t)}_{:,j})\,.$

We now consider the event that column $i$ collapses to the all-ones space.
Notice that dividing the whole column by $\Gamma_{i}$ does not change whether
or not the column has converged to the all-ones space or not. Thus,

$\displaystyle\mathcal{A}$ $\displaystyle=\\{W^{(t)}\in\mathbb{R}^{k\times
k}\mid X^{(t+1)}_{:,i}=\beta 1\\}$
$\displaystyle=\\{W^{(t)}\in\mathbb{R}^{k\times
k}\mid\frac{1}{\Gamma_{i}}(\varphi^{(i)}+W^{(t)}_{i,i}\tilde{X}^{(t)}_{:,i}+W^{(t)}_{j,i}\tilde{X}^{(t)}_{:,j})=\beta\mathbf{1}\\}$
$\displaystyle=\\{W^{(t)}\in\mathbb{R}^{k\times
k}\mid\varphi^{(i)}+W^{(t)}_{i,i}\tilde{X}^{(t)}_{:,i}+W^{(t)}_{j,i}\tilde{X}^{(t)}_{:,j}=\beta^{\prime}\mathbf{1}\\}$
$\displaystyle=\\{W^{(t)}\in\mathbb{R}^{k\times k}\mid
W^{(t)}_{i,i}\tilde{X}^{(t)}_{:,i}+W^{(t)}_{j,i}\tilde{X}^{(t)}_{:,j}-\beta^{\prime}\mathbf{1}=-\varphi^{(i)}\\}$

Since $\tilde{X}^{(t)}_{:,i},\tilde{X}^{(t)}_{:,j},\mathbf{1}$ are linearly
independent, given $\varphi^{(i)}$, there is only 1 solution to this equation.
$\mathcal{A}$ is a proper hyperplane in $\mathbb{R^{k\times k}}$ and as such
has Lebesgue measure $0$. The event $\mathcal{A}$ thus has probability $0$ and
the opposite event, that column $i$ does not collapse to the all-ones space,
has probability $1$.

The same holds for $X^{(t+1)}_{:,j}$ and by the same reasoning,
$X^{(t+1)}_{:,i}$ and $X^{(t+1)}_{:,j}$ are linearly independent with
probability $1$. Notice, that thus $\operatorname{Rank}(V_{\neq
0}^{\top}\begin{bmatrix}X^{(t+1)}_{:,i},X^{(t+1)}_{:,j}\end{bmatrix})=2$ still
holds. Finally, because both $\tilde{X}^{(t)}_{:,i}$ and
$\tilde{X}^{(t)}_{:,j}$ have variance greater than $0$ and the rescaling will
make it so that they have variance exactly $1$. This means, that
$X^{(t+1)}_{:,i},X^{(t+1)}_{:,j}$ has norm $1$ and all columns are centered.
We conclude that with probability $1$,

$\mu(X^{(t+1)})\vcentcolon=\left\|X^{(t+1)}-\mathbf{1}\mathbf{1}^{\top}X^{(t+1)}/n\right\|_{F}=\left\|X^{(t+1)}\right\|_{F}\geq\sqrt{2}\,.$

∎

###  B.7 Proposition 4.5: general conditions

We adopt the following regularity conditions:

######  Assumption B.4.

For the system described in (5), assume:

  1. 1.

For nonzero $x\in\mathbb{R}^{n}$ such that $x^{\top}\mathbf{1}=0$,
$Ax\notin\operatorname{Span}\\{\mathbf{1}\\}$.

  2. 2.

If $\mu((AX^{(t)})_{:,i})>0$, then $\mu((AX^{(t)}W^{(t)})_{:,i})>0$ for all
$i\in[k]$.

B.4.1 ensure that the adversarial situation where one step of message-passing
automatically leads to oversmoothing does not happen. Note that this is a more
relaxed condition than requiring $\operatorname{Span}\\{1\\}^{\perp}$ to be an
invariant subspace of $A$.

B.4.2 ensures that the case where weights are deliberately chosen for
oversmoothing to happen in one layer does not happen, either. Note that for
the second condition, such weights exist, i.e. let $W^{(t)}=I_{k}$. Moreover,
if weights are randomly initialized, then B.4.2 holds almost surely.

We restate Proposition 4.5 under general conditions, which accounts for both
linear and non-linear systems:

######  Proposition B.5.

For the system in (5), suppose $\sigma(\cdot)$ is injective and B.4 holds.
Without loss of generality, also suppose the initial features $X^{(0)}$ are
centered and all columns are nonzero. If $\mu(X^{(0)})>0$, then
$\mu(X^{(t)})=\sqrt{k}$ for all $t\geq 1$.

###### Proof.

For each column in $X^{(0)}$, since it is centered and nonzero,
$X^{(0)}_{:,i}\in\operatorname{Span}\\{1\\}^{\perp}\backslash\\{\mathbf{0}\\}$
for all $i\in[k]$.

Then given B.4.1, $AX^{(0)}_{:,i}\notin\operatorname{Span}\\{\mathbf{1}\\}$
and thus can be written as

$AX^{(0)}_{:,i}=a\mathbf{1}+w_{0}^{\perp}\,,$

where
$a\in\mathbb{R},w_{0}^{\perp}\in\operatorname{Span}\\{\mathbf{1}\\}^{\perp}$
and $w_{0}^{\perp}\neq\mathbf{0}$.

Then given B.4.2,
$AX^{(0)}W^{(0)}_{:,i}\notin\operatorname{Span}\\{\mathbf{1}\\}$ and since
$\sigma(\cdot)$ is injective,
$\sigma(AX^{(0)}W^{(0)}_{:,i})\notin\operatorname{Span}\\{\mathbf{1}\\}$ and

$\sigma(AX^{(0)}W^{(0)}_{:,i})=b\mathbf{1}+w_{1}^{\perp}\,,$

where
$b\in\mathbb{R},w_{1}^{\perp}\in\operatorname{Span}\\{\mathbf{1}\\}^{\perp}$
and $w_{1}^{\perp}\neq\mathbf{0}$.

Then after the centering step of batch normalization,

$(I-\mathbf{1}\mathbf{1}^{\top}/n)\sigma(AX^{(0)}W^{(0)}_{:,i})=w_{1}^{\perp}\,,$

and after the scaling step of batch normalization,

$X^{(1)}_{:,i}=\operatorname{BN}(\sigma(AX^{(0)}W^{(0)}_{:,i}))=\frac{w_{1}^{\perp}}{\|w_{1}^{\perp}\|_{2}}\,,$

which means each column of
$X^{(1)}\in\operatorname{Span}\\{\mathbf{1}\\}^{\perp}\backslash\\{\mathbf{0}\\}$
and has norm $1$. This directly implies $\mu(X^{(1)})=\sqrt{k}$.

The above argument applies for all $t\geq 0$ going from $X^{(t)}$ to
$X^{(t+1)}$, which concludes the proof.

∎

###  B.8 If nonlinearity $\sigma(\cdot)$ is applied after batch normalization

Suppose the GNN with batcn normalization defined in (5) is instead

$X^{(t+1)}=\sigma(Y^{(t+1)}),\quad
Y^{(t+1)}=\operatorname{BN}(AX^{(t)}W^{(t)})\,.$ (11)

Suppose $\mathbf{1}$ is an eigenvector of $A$. Then without
$\operatorname{BN}(\cdot)$, [36] showed that $\mu(X^{(t)})$ under appropriate
assumptions on the weight matrices. We show the following claim for $X^{(0)}$
where $\mu(X^{(0)}_{:,i})>0$:

######  Proposition B.6.

Suppose $\mathbf{1}$ is an eigenvector of $A$. If $\sigma(\cdot)$ is
injective, then under B.4.2, $\mu(X^{(t)}_{:,i})>0$ for all $t\geq 1$ and
$i\in[k]$.

###### Proof.

Given $\mu(X^{(0)}_{:,i})>0$,
$X^{(0)}_{:,i}\notin\operatorname{Span}\\{\mathbf{1}\\}$ and hence

$X^{(0)}_{:,i}=a\mathbf{1}+w_{0}^{\perp}\,,$

where
$a\in\mathbb{R},w_{0}^{\perp}\in\operatorname{Span}\\{\mathbf{1}\\}^{\perp}$
and $w_{0}^{\perp}\neq\mathbf{0}$.

Since $\mathbf{1}$ is an eigenvector of $A$,

$AX^{(0)}_{:,i}=b\mathbf{1}+w_{1}^{\perp}\,,$

where $b\in\mathbb{R},w_{1}^{\perp}\in\operatorname{Span}\\{1\\}^{\perp}$ and
$w_{1}^{\perp}\neq\mathbf{0}$, meaning that
$AX^{(0)}_{:,i}\notin\operatorname{Span}\\{\mathbf{1}\\}$. Then given B.4.2,
$AX^{(0)}W^{(0)}_{:,i}\notin\operatorname{Span}\\{\mathbf{1}\\}$ and

$AX^{(0)}W^{(0)}_{:,i}=c\mathbf{1}+w_{2}^{\perp}\,,$

where
$c\in\mathbb{R},w_{2}^{\perp}\in\operatorname{Span}\\{\mathbf{1}\\}^{\perp}$
and $w_{2}^{\perp}\neq\mathbf{0}$.

Then after the centering step of batch normalization,

$(I-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(0)}W^{(0)}_{:,i}=w_{2}^{\perp}\,,$

and after the scaling step of batch normalization,

$Y^{(1)}_{:,i}=\operatorname{BN}(AX^{(0)}W^{(0)}_{:,i})=\frac{w_{2}^{\perp}}{\|w_{2}^{\perp}\|_{2}}\,,$

which means each column of
$Y^{(1)}\in\operatorname{Span}\\{\mathbf{1}\\}^{\perp}\backslash\\{\mathbf{0}\\}$
and has norm $1$.

Then since $\sigma(\cdot)$ is injective, it follows that
$X^{(1)}_{:,i}=\sigma(Y^{(1)}_{:,i})\notin\operatorname{Span}\\{\mathbf{1}\\}$
and we conclude that $\mu(X^{(1)}_{:,i})>0$.

Note that the above argument applies for all $t\geq 0$ going from $X^{(t)}$ to
$X^{(t+1)}$, which concludes the proof.

∎

###  B.9 Proof for Proposition 4.7

See 4.7

###### Proof.

We will prove by induction on $t$ that

$X^{(t)}_{:,i}=\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sigma^{(t)}_{l,i}\lambda_{l}^{t}\nu_{l}$
(12)

where $(\lambda_{i},\nu_{i})$ is the $i$-th eigenpair of
$(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)A$ with
$|\lambda_{1}|\geq|\lambda_{2}|\geq\cdots\geq|\lambda_{n}|$, $\Gamma^{(t)}$ is
the normalization factor in the $t$-th round and
$\sigma_{l,i}^{(t)}\in\mathbb{R}$. The base case follows from the
decomposition of $X^{(0)}$ in the eigenvector basis of
$(I_{n}-\mathbf{1}\mathbf{1}^{T}/n)A$: $X^{(0)}_{:,i}=\sum_{l=1}^{n}\langle
X^{(0)}_{:,i},\nu_{l}\rangle\nu_{l}$ (and the fact that $X^{(0)}$ is
normalized).
For the induction step, the system can be rewritten as

$\displaystyle X^{(t+1)}$
$\displaystyle=(AX^{(t)}W^{(t)}-\mathbf{1}\mathbf{1}^{\top}AX^{(t)}W^{(t)}/n)\operatorname{diag}(...,\frac{1}{\operatorname{var}((AX^{(t)}W^{(t)})_{:,i})},...))$
$\displaystyle=(I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)}W^{(t)}D_{\operatorname{var}}\,.$

Assuming Equation 12,

$((I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)})_{:,i}=\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sigma^{(t)}_{l,i}\lambda_{l}^{t+1}\nu_{l}\,.$

Further the action of $W^{(t)}$ is the following:

$\displaystyle((I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)}W^{(t)})_{:,i}$
$\displaystyle=\frac{1}{\Gamma^{(t)}}\sum_{j=1}^{k}{W^{(t)}_{j,i}}\sum_{l=1}^{n}\sigma^{(t)}_{l,j}\lambda_{l}^{t+1}\nu_{l}$
(13)
$\displaystyle=\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sum_{j=1}^{k}({W^{(t)}_{j,i}}\sigma^{(t)}_{l,j})\lambda_{l}^{t+1}\nu_{l}$
$\displaystyle=\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1}\nu_{l}\,.$

Notice that
$\sigma^{(t+1)}_{l,i}=\sum_{j=1}^{k}({W^{(t)}_{j,i}}\sigma^{(t)}_{l,j})$ or
equivalently, $\Sigma^{(t+1)}=\Sigma^{(t)}W^{(t)}$, where
$\Sigma^{(t)}=[\sigma^{(t)}_{l,i}]$. Thus,
$\Sigma^{(t+1)}=\Sigma^{(0)}W^{(0)}W^{(1)}\cdots W^{(t)}$. Lastly, the action
of $D_{\operatorname{var}}$ is:

$\displaystyle((I_{n}-\mathbf{1}\mathbf{1}^{\top}/n)AX^{(t)}W^{(t)}D_{\operatorname{var}})_{:,i}$
$\displaystyle=\frac{1}{||\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1}\nu_{l}||_{2}}\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1}\nu_{l}$
$\displaystyle=\frac{1}{|\frac{1}{\Gamma^{(t)}}\sqrt{\sum_{l=1}^{n}(\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1})^{2}}|}\frac{1}{\Gamma^{(t)}}\sum_{l=1}^{n}\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1}\nu_{l}$
$\displaystyle=\frac{1}{|\sqrt{\sum_{l=1}^{n}(\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1})^{2}}|}\sum_{l=1}^{n}\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1}\nu_{l}$
$\displaystyle=\frac{1}{\Gamma^{(t+1)}}\sum_{l=1}^{n}\sigma^{(t+1)}_{l,i}\lambda_{l}^{t+1}\nu_{l}\,.$

This concludes the induction. Notice that superscripts in brackets do not
denote exponentiation, but rather an iterate at iteration $t$. Using this
intermediate result (12), we show the following:

######  Lemma B.7.

For all $q>k$,

$||\nu_{q}^{\top}X^{(t)}||_{2}\leq
C_{0}\left(\frac{\lambda_{q}}{\lambda_{k}}\right)^{t}\,.$

Proving this directly yields Proposition 4.7. We have that
$\Sigma^{(0)}=V^{\top}X^{(0)}$ due to the base case of the induction and by
assumption, $\Sigma^{(0)}_{:k,:}=V_{k}^{\top}X^{(0)}\in\mathbb{R}^{k\times k}$
has rank $k$ and is therefore full rank. It is therefore invertible, meaning
that there exists $(\Sigma_{:k,:}^{(0)})^{-1}\in\mathbb{R}^{k\times k}$ such
that $V_{k}^{\top}X^{(0)}(\Sigma_{:k,:}^{(0)})^{-1}=I_{k}$.

We can thus write for the simplicity of notation:

$\Sigma^{(0)}=\Sigma^{(0)}(\Sigma_{:k,:}^{(0)})^{-1}\Sigma_{:k,:}^{(0)}=\Sigma^{(\bot)}W^{(\bot)}=\begin{bmatrix}I_{k}\\\
\Sigma^{(\bot)}_{(k+1):,:}\end{bmatrix}W^{(\bot)}\,.$

This has the nice property that $\sigma^{(\bot)}_{i,i}=1$ for $i\leq k$. Now,
let’s revisit Equation 13 in that we can write
$\Sigma^{(t+1)}=\Sigma^{(0)}W^{(0)}W^{(1)}\cdots W^{(t)}$. Let
$\mathcal{W}^{(t)}=W^{(\bot)}W^{(0)}\cdots W^{(t)}$, meaning that
$\sigma_{i,j}^{(t)}=\sum_{l=1}^{k}\mathcal{W}^{(t)}_{i,l}\sigma^{(\bot)}_{l,j}$.
We now have everything to conclude the proof:

Consider now, the contribution of an eigenvector $\nu_{q}$ with $q>k$ at
iteration $t$.

$\displaystyle||\nu_{q}^{\top}X^{(t)}||_{2}$
$\displaystyle=\sqrt{\sum_{i=1}^{k}(\nu_{q}^{\top}X^{(t)}_{:,i})^{2}}=\sqrt{\sum_{i=1}^{k}(\frac{1}{\Gamma_{i}^{(t)}}\sum_{l=1}^{n}\sigma^{(t)}_{l,i}\lambda_{l}^{t}\nu_{q}^{\top}\nu_{l})^{2}}=\sqrt{\sum_{i=1}^{k}(\frac{1}{\Gamma_{i}^{(t)}}\sigma^{(t)}_{q,i}\lambda_{q}^{t})^{2}}$
$\displaystyle=\sqrt{\sum_{i=1}^{k}\frac{(\sigma^{(t)}_{q,i}\lambda_{q}^{t})^{2}}{\sum_{p=0}^{n}(\sigma^{(t)}_{p,i}\lambda_{p}^{t})^{2}}}=\sqrt{\sum_{i=1}^{k}\frac{(\sum_{l=1}^{k}\mathcal{W}^{(t)}_{l,i}\sigma^{(\bot)}_{q,l}\lambda_{q}^{t})^{2}}{\sum_{p=1}^{n}(\sum_{l=1}^{k}\mathcal{W}^{(t)}_{l,i}\sigma^{(\bot)}_{p,l}\lambda_{p}^{t})^{2}}}$
$\displaystyle\leq\sqrt{\sum_{i=1}^{k}\frac{(\sum_{l=1}^{k}\mathcal{W}^{(t)}_{l,i}\sigma^{(\bot)}_{q,l}\lambda_{q}^{t})^{2}}{\sum_{p=1}^{k}(\mathcal{W}^{(t)}_{p,i}\sigma^{(\bot)}_{p,p}\lambda_{p}^{t})^{2}}}=\sqrt{\sum_{i=1}^{k}\frac{(\sum_{l=1}^{k}\mathcal{W}^{(t)}_{l,i}\sigma^{(\bot)}_{q,l}\lambda_{q}^{t})^{2}}{\sum_{p=1}^{k}(\mathcal{W}^{(t)}_{p,i}\lambda_{p}^{t})^{2}}}$
$\displaystyle\leq\sqrt{\sum_{l=1}^{k}(\sigma_{q,l}^{(\bot)})^{2}\sum_{i=1}^{k}\frac{\sum_{l=1}^{k}(\mathcal{W}^{(t)}_{l,i}\lambda_{q}^{t})^{2}}{\sum_{p=1}^{k}(\mathcal{W}^{(t)}_{p,i}\lambda_{p}^{t})^{2}}}$
$\displaystyle\leq\sqrt{\sum_{l=1}^{k}(\sigma_{q,l}^{(\bot)})^{2}\sum_{i=1}^{k}\frac{\sum_{l=1}^{k}(\mathcal{W}^{(t)}_{l,i}\lambda_{q}^{t})^{2}}{\sum_{p=1}^{k}(\mathcal{W}^{(t)}_{p,i}\lambda_{k}^{t})^{2}}}=\sqrt{\sum_{l=1}^{k}(\sigma_{q,l}^{(\bot)})^{2}k}\frac{\lambda_{q}^{t}}{\lambda_{k}^{t}}$
$\displaystyle\leq C_{0}\frac{\lambda_{q}^{t}}{\lambda_{k}^{t}}\,.$

As we have that $\lambda_{q}<\lambda_{k}$ by construction,
$\|\nu_{q}^{\top}X^{(t)}||_{2}\rightarrow 0$ exponentially as
$t\rightarrow\infty$. ∎

###  B.10 Proof for Proposition 4.8

See 4.8

###### Proof.

The proof idea is simple: we use Gaussian elimination to cancel out all
“top-k" eigenvectors but the one in that row and then use the power iteration
until the smaller eigenvectors are “drowned out" below the $\epsilon$ error
margin. Assuming the columns of $X^{(0)}_{:k,:}$ are linearly independent
(which is the case if it has rank $k$), this leads to the desired output:
choose

$W^{(0)}_{j,i}=\begin{cases}1\quad&\text{if }j=i\\\
-\frac{\sigma^{(0)}_{1,i}}{\sigma^{(0)}_{1,1}}&\text{if }j=1\text{ and }i\neq
1\\\ 0&\text{else}\end{cases}$

Then from Equation 13 for each column $X_{:,i}$, we get that

$X^{(1)}_{:,i}=\frac{1}{\Gamma_{i}^{(1)}}\sum_{l=1}^{n}\sum_{j=1}^{k}({W^{(0)}_{j,i}}\sigma^{(0)}_{l,j})\lambda_{l}^{1}\nu_{l}=\frac{1}{\Gamma_{i}^{(1)}}\sum_{l=1}^{n}(\sigma^{(0)}_{l,i}-\frac{\sigma^{(0)}_{1,i}}{\sigma^{(0)}_{1,1}}\sigma^{(0)}_{l,1})\lambda_{l}^{1}\nu_{l}\,.$

Now for $l=1$, the factor
$\sigma^{(0)}_{l,i}-\frac{\sigma^{(0)}_{1,i}}{\sigma^{(0)}_{1,1}}\sigma^{(0)}_{l,1}=0$
yields

$X^{(1)}_{:,i}=\frac{1}{\Gamma_{i}^{(1)}}\sum_{l=2}^{n}(\sigma^{(0)}_{l,i}-\frac{\sigma^{(0)}_{1,i}}{\sigma^{(0)}_{1,1}}\sigma^{(0)}_{l,1})\lambda_{l}^{1}\,.$

Iterating this idea and choosing

$W^{(k)}_{j,i}=\begin{cases}1\quad&\text{if }j=i\\\
-\frac{\sigma^{(k)}_{k+1,i}}{\sigma^{(k)}_{k+1,k+1}}&\text{if }j=k+1\text{ and
}i\neq k+1\\\ 0&\text{else}\,,\end{cases}$

we arrive at

$X^{(k)}_{:,i}=\frac{1}{\Gamma_{i}^{(k)}}\left(\sigma^{(i-1)}_{i,i}\lambda_{i}^{k}\nu_{i}+\sum_{l=k+1}^{n}\sigma^{(k)}_{l,i}\lambda_{l}^{k}\nu_{l}\right)\,.$

Now, we switch gears and use $W^{(t)}=I_{n}$ for $t\geq k$, it follows that

$\displaystyle||\nu_{i}^{\top}(X^{(t)}_{:,i})||_{2}$
$\displaystyle=\frac{\sigma^{(i-1)}_{i,i}\lambda_{i}^{t}}{\sqrt{(\sigma^{(i-1)}_{i,i}\lambda_{i}^{t})^{2}+\sum_{l=k+1}^{n}(\sigma^{(k)}_{l,i}\lambda_{l}^{t})^{2}}}$
$\displaystyle=\frac{1}{\sqrt{1+\frac{\sum_{l=k+1}^{n}(\sigma^{(k)}_{l,i}\lambda_{l}^{t})^{2}}{(\sigma^{(i-1)}_{i,i}\lambda_{i}^{t})^{2}}}}\,.$

The statement left to prove is thus:

$\displaystyle\quad\epsilon\geq\frac{\sum_{l=k+1}^{n}(\sigma^{(k)}_{l,i}\lambda_{l}^{t})^{2}}{(\sigma^{(i-1)}_{i,i}\lambda_{i}^{t})^{2}}$
$\displaystyle\Longleftarrow\epsilon\geq\frac{(n-k)\max_{l}(\sigma^{(k)}_{l,i})^{2}}{(\sigma_{i,i}^{i-1})^{2}}\frac{\lambda_{k+1}^{2t}}{\lambda_{i}^{2t}}$
$\displaystyle\Longleftarrow\frac{\log\left(\frac{\epsilon\sigma^{(i-1)}_{i,i}}{(n-k)\max_{l}(\sigma^{(k)}_{l,j})^{2}}\right)}{2\log\left(\frac{\lambda_{k+1}}{\lambda_{i}}\right)}\leq
t\,.$

Thus, setting $T$ to be larger than this bound yields the desired claim. ∎

###  B.11 Proof for Proposition 5.1

######  Proposition B.8.

Let $A\in\mathbb{R}^{n\times n}_{\geq 0}$ be a symmetric non-negative matrix.
Let $H\in\\{0,1\\}^{n\times m}$ such that $AH=HA^{\pi}$ is the coarsest EP of
$A$. Divide the eigenpairs $\mathcal{V}=\\{...,(\lambda,\nu),...\\}$ of $A$
into the following two sets:
$\mathcal{V}_{\text{struc}}=\\{(\lambda,\nu)\in\mathcal{V}\mid\nu=H\nu^{\pi}\\},\mathcal{V}_{\text{rest}}=\mathcal{V}\backslash\mathcal{V}_{\text{struc}}$.
Let $\hat{\mathcal{V}}=\\{...,(\hat{\lambda},\hat{\nu}),...\\}$ be eigenpairs
of $(I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)A$, for $\tau\neq 0$. Then

  1. 1.

$\mathcal{V}_{\text{rest}}\subset\hat{\mathcal{V}}$.

  2. 2.

Assume $\mathcal{V}_{\text{struc}}\neq\\{(\lambda,\mathbf{1})\\}$. Let
$(\lambda,\nu)$ be the dominant eigenpair of $A$. Then $\nu$ is _not_ an
eigenvector of $(I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)A$.

  3. 3.

$\sum_{(\lambda,\nu)\in\mathcal{V}_{\text{struc}}}\lambda>\sum_{(\hat{\lambda},\hat{\nu})\in\hat{\mathcal{V}}\backslash\mathcal{V}_{\text{rest}}}\hat{\lambda}$

###### Proof.

Let $H\in\\{0,1\\}^{n\times m}$ indicate the coarsest EP of $A$
($AH=HA^{\pi}$). As each node belongs to exactly $1$ class, it holds that
$H\mathbf{1}_{n}=\mathbf{1}_{m}$. We prove that for any eigenpair
$(\lambda,\nu)\in\mathcal{V}_{\text{rest}}$, it holds that
$\nu^{\top}\mathbf{1}=0$. From this, the first statment follows quickly:

$(I_{n}-\tau\mathbf{1}\mathbf{1}^{T}/n)A\nu=A\nu-\tau\mathbf{1}\mathbf{1}^{T}A\nu/n=\lambda\nu-\tau\lambda\mathbf{1}\mathbf{1}^{T}\nu/n=\lambda\nu\,.$

To prove this, let’s look at $A^{\pi}$. $A^{\pi}$ is not symmetric, but its
eigenpairs are associated to $A$’s eigenpairs in the following way. Let
$(\lambda,\nu^{\pi})$ be an eigenpair of $A^{\pi}$, then
$(\lambda,H\nu^{\pi})$ is an eigenpair of $A$:
$AH\nu^{\pi}=HA^{\pi}=H\lambda\nu^{\pi}$. As $A$ is a symmetric matrix, its
eigenvectors are orthogonal implying
$(H\nu_{i}^{\pi})^{\top}(H\nu_{j}^{\pi})=0$.

######  Lemma B.9.

The eigenvectors $\nu_{1}^{\pi},...,\nu_{m}^{\pi}$ of $A^{\pi}$ are linearly
independent.

For simplicity, choose the eigenvectors $\nu_{i}^{\pi}$ to be normalized in
such a way, that $(H\nu_{i}^{\pi})^{\top}(H\nu_{i}^{\pi})=1$. Suppose for a
contradiction, that they are linearly dependent and without loss of
generality, $\nu_{m}^{\pi}=\sum_{i=1}^{m-1}a_{i}\nu_{i}^{\pi}$. Take $j$ such
that $a_{j}\neq 0$, this must exist otherwise $\nu_{m}=\mathbf{0}$. Now,

$(H\nu_{m}^{\pi})^{\top}(H\nu_{j}^{\pi})=(H\sum_{i=1}^{m-1}a_{i}\nu_{i}^{\pi})^{\top}(H\nu_{j}^{\pi})=\sum_{i=1}^{m-1}(Ha_{i}\nu_{i}^{\pi})^{\top}(H\nu_{j}^{\pi})=a_{j}\neq
0\,,$

which is a contradiction.

Since the eigenvectors of $A^{\pi}$ are linearly independent, there exists a
unique $\beta\in\mathbb{R}^{m}$ s.t.
$\sum_{i=1}^{m}\beta_{i}\nu_{i}^{\pi}=\mathbf{1}_{m}$. Finally let
$\varphi\in\mathcal{V}_{\text{rest}}$,

$\varphi^{\top}\mathbf{1}_{n}=\varphi^{\top}H\mathbf{1}_{m}=\varphi^{\top}H\sum_{i=1}^{m}\beta_{i}\nu_{i}^{\pi}=\sum_{i=1}^{m}\beta_{i}\varphi^{\top}H\nu_{i}^{\pi}=0\,.$

This concludes the proof of the first statement.

To prove the second statement, let $(\lambda,\nu)$ be a dominant eigenpair of
$A$, such that $\nu\geq 0$ is non-negative. This exists as $A$ is a non-
negative matrix. Additionally, $\nu$ is not the all-zeros vector and as such
$\mathbf{1}^{\top}\nu>0$. By assumption $\nu\neq\mathbf{1}$. Then:

$\displaystyle(I-\tau\mathbf{1}\mathbf{1}^{\top}/n)A\nu$
$\displaystyle=\hat{\lambda}\nu\iff\lambda\nu-\tau\mathbf{1}\mathbf{1}^{\top}\lambda\nu/n$
$\displaystyle=\hat{\lambda}\nu\iff(\lambda-\hat{\lambda})\nu-\tau\lambda\mathbf{1}\mathbf{1}^{\top}\nu/n$
$\displaystyle=0\,.$

As $\nu\neq\mathbf{1}$, there exist $i,j$ s.t. $\nu_{i}\neq\nu_{j}$ thus,
$(\lambda-\hat{\lambda})\nu_{i}-\tau\lambda(\mathbf{1}\mathbf{1}^{\top}\nu)_{i}/n=0$
and
$(\lambda-\hat{\lambda})\nu_{j}-\tau\lambda(\mathbf{1}\mathbf{1}^{\top}\nu)_{j}/n=0$
cannot be true at the same time. Thus, this equation has no solution and $\nu$
is not an eigenvector of $(I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)A$.

For the third and final statement, notice that the trace of $A$ is larger than
the trace of $(I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)A$. Since the trace of
a matrix is the sum of its eigenvalues, we have

$\sum_{(\lambda,\nu)\in\mathcal{V}}\lambda=\operatorname{Tr}(A)>\operatorname{Tr}((I_{n}-\tau\mathbf{1}\mathbf{1}^{\top}/n)A)=\sum_{(\hat{\lambda},\hat{\nu})\in\hat{\mathcal{V}}}\hat{\lambda}\,.$

Consequently,

$\sum_{(\lambda,\nu)\in\mathcal{V}}\lambda=\sum_{(\lambda,\nu)\in\mathcal{V}_{\text{struc}}}\lambda+\sum_{(\lambda,\nu)\in\mathcal{V}_{\text{rest}}}\lambda>\sum_{(\hat{\lambda},\hat{\nu})\in\hat{\mathcal{V}}\backslash\mathcal{V}_{\text{rest}}}\hat{\lambda}+\sum_{(\lambda,\nu)\in\mathcal{V}_{\text{rest}}}\lambda=\sum_{(\hat{\lambda},\hat{\nu})\in\hat{\mathcal{V}}}\hat{\lambda}\,.$

Subtracting $\sum_{(\lambda,\hat{\nu})\in\mathcal{V}_{\text{rest}}}\lambda$
from both sides yields the final statement.

∎

##  Appendix C Experiments

We compare these models using the measures of convergence: $\mu(X^{(t)})$ as
defined in (3), the numerical rank of the features
$\operatorname{Rank}(X^{(t)})$, the column distance used in [52]:

$d_{\text{col}}(X)\vcentcolon=\frac{1}{d^{2}}\sum_{i,j}\left\|\frac{X_{:,i}}{\left\|X_{:,i}\right\|_{1}}-\frac{X_{:,j}}{\left\|X_{:,j}\right\|_{1}}\right\|_{2},$

the column projection distance:

$d_{\text{p-col}}(X)\vcentcolon=\frac{1}{d^{2}}\sum_{i,j}1-\frac{X_{:,i}^{\top}}{\left\|X_{:,i}\right\|_{2}}\frac{X_{:,j}}{\left\|X_{:,j}\right\|_{2}},$

and the eigenvector space projection:

$d_{\text{ev}}(X)\vcentcolon=\frac{1}{n}\left\|X-VV^{\top}X\right\|_{F},$

where $V\in\mathbb{R}^{n\times n}$ is the set of normalized eigenvectors of
$A$, and finally, the rank of $X$: $\operatorname{Rank}(X)$. The same
experiment with both GIN and GAT can be found below. The trends are similar to
those described in the main manuscript.

Figure 2: Long Term behavior of GCN. Progression of the above measures of
convergence over $256$ iterations of message-passing in both linear and non-
linear GCN.  Figure 3: Long Term behavior of GIN. Progression of the above
measures of convergence over $256$ iterations of message-passing in both
linear and non-linear GIN.  Figure 4: Long Term behavior of GAT. Progression
of the above measures of convergence over $256$ iterations of message-passing
in both linear and non-linear GAT.

#### Datasets

We provide summary statistics for datasets used in experiments in Section 6 in
Table 2.

Table 2: The summary statistics of the datasets used in Section 6. Dataset | #graphs | #nodes | #edges | #features | # classes
---|---|---|---|---|---
MUTAG | 188 |  $\sim$17.9 |  $\sim$39.6 | 7 | 2
PROTEINS | 1,113 |  $\sim$39.1 |  $\sim$145.6 | 3 | 2
PTC-MR | 344 |  $\sim$14.29 |  $\sim$14.69 | 18 | 2
Cora | 1 | 2,708 | 10,556 | 1,433 | 7
Citeseer | 1 | 3,327 | 9,104 | 3,703 | 6
PubMed | 1 | 19,717 | 88,648 | 500 | 3
ogbn-arxiv | 1 | 169,343 | 1,166,243 | 128 | 40

#### Training details

We perform a within-fold 90%/10% train/validation split for model selection.
We train the models for $200$ epochs using the AdamW optimizer and search the
hyperparameter space over the following parameter combinations:

  * •

learning rate $\in\\{10^{-4},10^{-3},10^{-2},10^{-1}\\}$

  * •

feature size $\in\\{32,64\\}$

  * •

weight decay $\in\\{0,10^{-2},10^{-4}\\}$

  * •

number of layers $\in\\{3,5\\}$

We select the hyperparameters of the model with the best mean validation
accuracy over its 30 best epochs. The code and all non publicly available data
will be made available here.

#### Compute

We ran all of our experiments on a system with two NVIDIA L40 GPUs, two AMD
EPYC 7H12 CPUs and 1TB RAM.

#### Licenses

  * •

PyG [19]: MIT license

  * •

OGB [25]: MIT license

  * •

ogbn-arxiv: ODC-BY