Taming Gradient Oversmoothing and Expansion in Graph Neural Networks

## Abstract

Abstract Oversmoothing has been claimed as a primary bottleneck for multi-layered graph neural networks (GNNs). Multiple analyses have examined how and why oversmoothing occurs. However, none of the prior work addressed how optimization is performed under the oversmoothing regime. In this work, we show the presence of gradient oversmoothing preventing optimization during training. We further analyze that GNNs with residual connections, a well-known solution to help gradient flow in deep architecture, introduce gradient expansion , a phenomenon of the gradient explosion in diverse directions. Therefore, adding residual connections cannot be a solution for making a GNN deep. Our analysis reveals that constraining the Lipschitz bound of each layer can neutralize the gradient expansion. To this end, we provide a simple yet effective normalization method to prevent the gradient expansion. An empirical study shows that the residual GNNs with hundreds of layers can be efficiently trained with the proposed normalization without compromising performance. Additional studies show that the empirical observations corroborate our theoretical analysis.

## 1 Introduction

Increasing depths of layers is recognized as a key reason for the success of many deep neural networks. Despite this widely believed claim, the graph neural network (GNN) does not necessarily benefit from deeper architectures. In fact, increasing the number of layers in GNNs often leads to problems such as oversmoothing, where node representations from different classes become indistinguishable, resulting in a degradation of performance. Although residual connections contribute to the success of deep architectures in many types of neural networks, they are generally ineffective at resolving the oversmoothing in GNNs.

To provide a deeper understanding of the oversmoothing, several studies have addressed the mechanism behind the oversmoothing. Li et al. () explain that the graph convolutional network (GCN) is a special form of Laplacian smoothing. Chamberlain et al. () show that GNN can be interpreted as a discretization of the heat diffusion process, implying the diffusion of the node representation to reach equilibrium. Cai & Wang () and Oono & Suzuki () theoretically show that oversmoothing occurs exponentially when piling the layers in GCN under certain conditions. Wu et al. () show Graph Attention Networks (GATs), a generalization of GCN, also lose their expressive power exponentially with the number of layers. However, none of the previous work on oversmoothing investigates the gradient flow of the deep GNNs.

In this work, we investigate gradient flows under the oversmoothing regime. Understanding the exact characteristics of gradients can help mitigate the oversmoothing problem. To do this, we first analyze the gradients with the graph convolutional network (GCN). The analysis shows the presence of *gradient oversmoothing*, which implies a gradient of node features are similar to each other. We then investigate the effect of the residual connections in GNNs, which is known to make the optimization work with many variants of the neural networks. Interestingly, we show that the opposite effect, named *gradient expansion*, outbreaks as a side effect of the residual connection. When the gradient expansion occurs, the gradients of node features explode in diverse directions.
We find that the gradient expansion is tightly connected with the upper bound of the Lipschitz constants of GNN layers. To relax the gradient expansion, we propose a simple weight normalization that controls the Lipschitz upper bound.

Empirical studies on node classification tasks corroborate our theoretical analysis, showing that gradient smoothing and gradient expansion happen with deep GNNs and deep residual GNNs, respectively. We also show that the similarity between gradients better represents the test performance than the similarity between node representations, a commonly used oversmoothing measure. Finally, we present how the deep GNNs perform with the Lipschitz constraints. With the Lipschitz constraints, we can successfully train the GNNs with residual connections without sacrificing the test accuracy on node classification tasks.

## 2 Related Work

Recently, many studies have investigated the reason for performance degradation in deep GNNs.
Oversmoothing, the phenomenon where the representations of neighboring nodes become similar as the number of GNN layers increases, has been identified as one of the main reasons.
Several studies have theoretically demonstrated that oversmoothing occurs due to the structural characteristics of message-passing GNNs.
Li et al. () reveal that GCN is a special form of Laplacian smoothing and points out that repeated applications lead to convergence to the same value.
Oono & Suzuki () demonstrate that, as the depth of the GCN architecture increases, the distance between representations and invariant spaces decays exponentially, leading to oversmoothing.
Keriven () show that, in linear GNNs whose aggregation matrix takes the form of a stochastic matrix, the output converges to the average value of the training labels as the depth approaches infinity.
Wu et al. () analyze the mixing and denoising effects of GNN layers on graphs sampled from the contextual stochastic block model, providing insight into why oversmoothing can occur even before depth approach infinity.
Wu et al. () prove that oversmoothing occurs exponentially even in the graph attention network architectures with non-linear activations.

Based on these analyses, many studies try to mitigate oversmoothing to stack multiple layers without performance degradation.
PairNorm (Zhao & Akoglu, ) keeps the total pairwise representation distances constant by normalization.
DropEdge (Rong et al., ) randomly drops edges from the graph to slow down oversmoothing.
Energetic graph neural networks  (Zhou et al., ) utilize the Dirichlet energy term, which measures oversmoothing in the loss function to regularize oversmoothing.
Ordered GNN (Song et al., ) preserves aggregated representation within specific hops separately.
Gradient Gating  (Rusch et al., ) proposes a mechanism to stop learning in a node-wise fashion before local oversmoothing occurs. Park et al. () propose to reverse the aggregation process.
Some studies propose to change the dynamics of GNN, based on the findings of Chamberlain et al. (), which shows that classical GNN resembles diffusion-like dynamics.
GRAND++ (Thorpe et al., ) proposes to add a source term to prevent oversmoothing.
PDE-GCN (Eliasof et al., ) proposes wave-like GNN architecture.
Graph-coupled oscillator networks (Rusch et al., ) propose a network based on non-linear oscillator structures.

Some studies propose to improve the optimization process for multi-GNN layer stacking. Li et al. () adopt residual connection to improve the training process.
Zhang et al. () also point out model degradation as a main reason for the performance degradation of deep GNNs and proposes an adaptive initial residual to improve optimization.
However, no studies explain the exact reason for model degradation and residual connection’s effectiveness.

## 3 Preliminary

### 3.1 Graph convolutional networks

We begin by considering undirected graphs $G=({\mathcal{V}},{\mathcal{E}},{\mathbf{X}})$, where ${\mathcal{V}}$ represents the set of $N\in{\mathbb{N}}$ *nodes*, ${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$ forms the *edge* set, a symmetric relation, and ${\mathbf{X}}\in{\mathbb{R}}^{N\times d}$ is a collection of $d$-dimensional feature matrix for each node. We use ${\mathbf{X}}_{i}\in{\mathbb{R}}^{d}$ to represent the feature vector of node $i$. The connectivity between nodes can also be represented through an *adjacency matrix* ${\mathbf{A}}\in\{0,1\}^{N\times N}$.
We focus on a graph convolutional network (GCN) (Kipf & Welling, ) and its variants.
With an initial input feature matrix ${\mathbf{X}}^{(0)}\coloneqq{\mathbf{X}}$, at each layer $0\leq\ell<L$, the *representation of nodes* ${\mathbf{X}}^{(\ell)}$ is updated as follows:

$$ $\displaystyle{\mathbf{X}}^{(\ell+1)}\coloneqq\sigma(\hat{{\mathbf{A}}}{\mathbf {X}}^{(\ell)}{\mathbf{W}}^{(\ell)}),$ (1) $$

where $\hat{{\mathbf{A}}}=\tilde{{\mathbf{D}}}^{-\frac{1}{2}}\tilde{{\mathbf{A}}}
\tilde{{\mathbf{D}}}^{-\frac{1}{2}}$ is an normalized adjacency matrix with $\tilde{{\mathbf{A}}}={\mathbf{A}}+{\mathbf{I}}$, diagonal matrix $\tilde{D}_{ii}=\sum_{j}\tilde{A}_{ij}$, $\sigma$ is non-linear activation function, and ${\mathbf{W}}^{(\ell)}$ a layer specific learnable weight matrix. Linear Graph Neural networks (LGN) represent the GCN without an activation function.

There are several variations of GCN through different definitions of the normalized adjacency matrix. For example, Graph Attention Networks (GATs)  use a learnable adjacency matrix through an attention layer (Veličković et al., ). Regardless of the precise formulation of the normalized matrix, most variations share similar characteristics; the node representations are updated by aggregating the neighborhood representations.

The residual connections can be incorporated into the GNN to help the gradient flow with many layers. The update rule of the node representation with the residual connection is as follows:

$$ $\displaystyle{\mathbf{X}}^{(\ell+1)}\coloneqq\sigma(\hat{{\mathbf{A}}}{\mathbf {X}}^{(\ell)}{\mathbf{W}}^{(\ell)})+{\mathbf{X}}^{(\ell)}.$ (2) $$

### 3.2 Similarity measure

We follow the definition of node similarity measure of node feature matrix ${\mathbf{X}}$ defined in Wu et al. () as

$$ $\mu\left({{\mathbf{X}}}\right)\coloneqq\lVert{{\mathbf{X}}-\mathbf{1}\gamma_{X }}\rVert_{F},\text{ where }\gamma_{X}=\frac{\mathbf{1}^{\top}{\mathbf{X}}}{N}\;,$ (3) $$

and $\lVert{\cdot}\rVert_{F}$ indicates Frobenius norm.
Considering ${\mathbf{B}}\in\mathbb{R}^{(N-1)\times(N-1)}$, the orthogonal projection into the space perpendicular to $\mathsf{span}{\mathbf{1}}$, the definition of $\mu$ satisfies $\mu\left({{\mathbf{X}}}\right)=\lVert{{\mathbf{B}}{\mathbf{X}}}\rVert_{F}$.
This definition of node similarity satisfies the following two axioms (Rusch et al., ):

- 1.
$\exists{\bm{c}}\in\mathbb{R}^{d}$ such that ${\mathbf{X}}_{i}={\bm{c}}$ for all node $i$ if and only if $\mu({\mathbf{X}})=0$, for ${\mathbf{X}}\in{\mathbb{R}}^{N\times d}$;
- 2.
$\mu({\mathbf{X}}+{\mathbf{Y}})\leq\mu({\mathbf{X}})+\mu({\mathbf{Y}})$, for all ${\mathbf{X}},{\mathbf{Y}}\in{\mathbb{R}}^{N\times d}$.

*Representation oversmoothing* with respect to $\mu$ is then characterized as the layer-wise convergence of the node similarity measure $\mu$ to zero, i.e.,

$$ $\displaystyle\lim_{\ell\rightarrow\infty}\mu\left({{\mathbf{X}}^{(\ell)}} \right)=0.$ $$

The node similarity measure can further be used to compute the similarity between the gradients of node features. Let $\frac{\partial{{\mathcal{L}}}}{\partial{{\mathbf{X}}^{(\ell)}}}\in{\mathbb{R}}
^{N\times d}$ be the partial derivatives w.r.t some loss function ${\mathcal{L}}$ given input feature matrix and output labels. With the same measure $\mu$, the gradient similarity at layer $\ell$ can be computed by $\mu\left({\frac{\partial{{\mathcal{L}}}}{\partial{{\mathbf{X}}^{(\ell)}}}}\right)$.
*Gradient smoothing* can happen when the gradient similarity converges to zero.
We further characterize the *gradient expansion*, which describes the cases where the gradient similarity of the first layer diverges, i.e.,

$$ $\displaystyle\lim_{L\to\infty,\ell\to 1}\mu\left({\frac{\partial{{\mathcal{L}} }}{\partial{{\mathbf{X}}^{(L-\ell)}}}}\right)=\infty.$ $$

Since $\lVert{{\mathbf{X}}}\rVert_{F}^{2}=\mu\left({{\mathbf{X}}}\right)^{2}+\lVert{{
\mathbf{X}}-{\mathbf{B}}{\mathbf{X}}}\rVert_{F}^{2}$, gradient expansion implies a gradient explosion. We note that, however, gradient expansion and gradient explosion are not the same concept.
The gradient expansion occurs when there is a *gradient explosion in diverse directions* and differs from the gradient explosion where we only measure the magnitude of the gradients.

## 4 Gradient Oversmoothing and Expansion

It is widely recognized that increasing the depth of graph neural networks (GNNs) typically results in performance degradation. While oversmoothing has been identified as a key cause, questions remain about why this issue cannot be resolved through optimization. One possible explanation is that optimization is inherently tricky in GNNs, such as gradient vanishing and exploding, commonly observed in the earlier stages of neural networks. However, it is also observed that the performance of deep GNNs has limitations with the more advanced techniques known to improve optimization, such as the residual connections (He et al., ).

To answer these questions, we thoroughly analyze the gradient of the linear GCNs with and without residual connections first. The analysis shows the presence of gradient oversmoothing that is similar to the representation oversmoothing but a fundamental reason why the training is incomplete with deep GCNs. Our analysis further identifies that gradient expansion is a more common problem with the residual connections, explaining the failure of the residual connections in deep GCNs. We then extend our analysis of the GNNs with a non-linear activation through empirical experiments to show that gradient oversmoothing and expansion can also be observed in practice.
We conduct our experiment also on GATs, showing that gradient oversmoothing and expansion is not the only problem in GCN, although our theoretical analyses are limited to GCN.

### 4.1 Theoretical Analysis

To understand how the optimization is performed withGCNs, we first analyze the gradients of the GCN during training. First, let us show the exact gradients of the parameter matrix at each layer during training.

###### Lemma 1 .

(Blakely et al., ) Let $\mathsf{LGN}$ be $L$-layered linear GCN without an activation function, and $\mathcal{L}_{\mathsf{LGN}}({\mathbf{W}},{\mathbf{X}},{\mathbf{y}})$ be a loss function of the linear GCN with a set of parameters ${\mathbf{W}}=({\mathbf{W}}^{(0)},\cdots,{\mathbf{W}}^{(L-1)})$ and input feature matrix ${\mathbf{X}}$ and output label ${\mathbf{y}}$. The gradient of the parameter at layer $\ell$, i.e., ${\mathbf{W}}^{(\ell)}$, with respect to the loss function is

$$ $\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{W}}}^{({\ell})}}=( \hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})})^{\top}\frac{\partial{\mathcal{L}_ {\mathsf{LGN}}}}{\partial{{{\mathbf{X}}}^{({\ell+1})}}}\;,$ $$

where

$$ $\frac{\partial{\mathcal{L}_{\mathsf{LGN}}}}{\partial{{{\mathbf{X}}}^{({\ell})} }}=(\hat{{\mathbf{A}}}^{\top})^{L-\ell}\frac{\partial\mathcal{L}_{\mathsf{LGN} }}{\partial{{\mathbf{X}}}^{({L})}}\prod_{i=1}^{L-\ell}({\mathbf{W}}^{(L-i)})^{ \top}\;.$ $$

Lemma 1 shows that the formulation of gradient $\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{\mathbf{W}}^{(\ell)}}$ resembles the forward propagation steps of the LGN. Specifically, $\frac{\partial{\mathcal{L}_{\mathsf{LGN}}}}{\partial{{{\mathbf{X}}}^{({\ell})}}}$ can be interpreted as the update equation of $(L-\ell)$-layered LGN with $\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{X}}}^{({L})}}$ as an input feature matrix, and $\hat{{\mathbf{A}}}^{\top}$ and ${\mathbf{W}}^{(\ell)\top}$ as a normalized adjacency and weight matrix at layer $\ell$, respectively. Consequently, following the general analysis of the oversmoothing in node representations, we can show the oversmoothing in gradients of the GNN. The following theorem characterized the upper bound on the similarity between the gradient of different nodes in a deeper layer.

###### Theorem 1 (Gradient oversmoothing in LGN) .

With the same assumptions used in Lemma 1 and a additional assumption that the graph $G$ is connected and non-bipartite, there exists $0<q<1$ and constant $C_{q}>0$ such that

$$ $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{LGN}}}}{\partial{{{\mathbf{X}}}^ {({\ell})}}}}\right)\leq C_{q}\left(q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2} \right)^{(L-\ell)},\quad 0\leq\ell<L\;,$ (4) $$

where $\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}$ is the maximum spectral norm, i.e., $*=\operatorname*{arg\,max}_{i}\lVert{{\mathbf{W}}^{(i)}}\rVert_{2}$ for $1\leq i\leq L-\ell$.

Theorem 1 implies that the similarity measure of the gradient decreases exponentially while approaching to the first layer if $q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}<1$. In other words, as the number of depths increases, the gradient information required to update the layers closer to the input becomes oversmoothed regardless of the input feature. When $q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}>1$, a gradient expansion occurs. As a result, optimization becomes extremely difficult in GNNs with deeper layers due to the structural characteristics.

Although Theorem 1 shows the challenges in GCN training, similar challenges are also prevalent in any neural network architectures.
The skip connection or residual connection is the off-the-shelf solution to improve the gradient flows in deep architectures (He et al., ).
The residual connections are also suggested in previous studies to improve the performance of GNNs with multiple layers (Li et al., ).
Although they successfully stacked $56$ layers, this is relatively shallow compared to ResNet, which stacks more than 152 layers.
To better understand the gradient flow with the residual connections in GCN, we first state the exact gradients of the parameter matrix at each layer.

###### Lemma 2 .

Let $\mathsf{resLGN}$ be a linear graph neural network with residual connections at each layer, and $\mathcal{L}_{\mathsf{resLGN}}({\mathbf{W}},{\mathbf{X}},{\mathbf{y}})$ be a loss function of $\mathsf{resLGN}$ with a set of parameters ${\mathbf{W}}=({\mathbf{W}}^{(0)},\cdots,{\mathbf{W}}^{(L-1)})$ and input feature matrix ${\mathbf{X}}$ and output label ${\mathbf{y}}$. The gradient of the parameter at layer $\ell$, i.e., ${\mathbf{W}}^{(\ell)}$, with respect to the loss function is

$$ $\frac{\partial\mathcal{L}_{\mathsf{resLGN}}}{\partial{{\mathbf{W}}}^{({\ell})} }=\left(\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})}\right)^{\top}\frac{ \partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({\ell+1})}} }\;,$ $$

where

$$ $\displaystyle\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf {X}}}^{({\ell})}}}$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}$ $\displaystyle+\sum_{p=1}^{L-\ell}\left(\underbrace{\sum_{\{i_{1}:i_{1}\geq\ell \}}\sum_{\{i_{2}:i_{2}>i_{1}\}}\cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}}_{\text{ all possible length-$p$ paths in back-propagation}}\left(\hat{\mathbf{A}}^{ \top}\right)^{p}\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{ \top}\right)\;.$ (5) $$

Compared to Lemma 1, the gradient of resLGN is expressed as the summation of $L-\ell+1$ terms, each of which considers all possible length $p$ back-propagation paths. Each path resembles the forward propagation rule of the LGN similar to Lemma 1.
As a result, the gradients are unlikely to be oversmoothed, as terms with shorter path lengths that are not oversmoothed still remain.
However, the issue of gradient explosion still exists as the number of layers $L$ increases.

The following theorem states the upper bound of the gradient similarity w.r.t. node representation at layer $\ell$.

###### Theorem 2 (Gradient expansion in resLGN) .

With the same assumptions used in Lemma 2, there exists $0<q<1$ and constant $C_{q}>0$, such that

$$ $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}({\mathbf{W}})}}{ \partial{{{\mathbf{X}}}^{({\ell})}}}}\right)\leq\mu\left({\frac{\partial{ \mathcal{L}_{\mathsf{resGNN}}({\mathbf{W}})}}{\partial{{{\mathbf{X}}}^{({L})}} }}\right)+C_{q}\left(\left(1+q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}\right)^{L- \ell}-1\right),\quad 0\leq\ell<L\;,$ (6) $$

where $*=\operatorname*{arg\,max}_{i}\lVert{{\mathbf{W}}^{(i)}}\rVert_{2}$ for $1\leq i\leq L-\ell$.

Theorem 2 implies that the gradients can *expand* with the residual connections when $q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}>0$.
However, unlike in LGNs, the upper bound on gradient similarity cannot converge to zero unless $q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}$ is zero. Therefore, with the residual connection, the gradient oversmoothing is less of a concern than the gradient expansion.

###### Lemma 1 .

###### Theorem 1 (Gradient oversmoothing in LGN) .

###### Lemma 2 .

###### Theorem 2 (Gradient expansion in resLGN) .

#### Lipschitz upper bound constraint.

The Lipschitz upper bound of a GCN layer is $\lVert{{\mathbf{W}}}\rVert_{2}$ since renormalization trick makes $\lVert{\tilde{{\mathbf{A}}}}\rVert_{2}=1$ (Kipf & Welling, ). Therefore, Theorem 2 implies that we can stabilize the training procedure and mitigate the gradient explosion by controlling the Lipschitz constant of the GCN layer. The exact $\lVert{{\mathbf{W}}}\rVert_{2}$ is the largest singular value of the parameter matrix, which is often difficult to compute.
In this work, we use a simple alternative to control the upper bound on the Lipschitz constant by normalizing the weight through its Frobenius norm, i.e., ${\mathbf{W}}\xleftarrow{}c{\mathbf{W}}/\lVert{{\mathbf{W}}}\rVert_{F}$, where $c$ is a hyperparameter controlling the upper bound of the Lipschitz constant as shown in Park et al. ().
We note that since the largest singular value of the right stochastic matrix corresponds to one, we can apply the same approach in GAT.
In our experiments, we normalize the weight matrices at each training iteration. The efficiency of the normalization is shown in Section 6.

Figure: Figure 1: Gradient similarity measure $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{GNN}}({\mathbf{W}})}}{\partial{{ {\mathbf{X}}}^{({\ell})}}}}\right)$ over different layers and activation functions of $128$-layer GCN and GAT in three datasets: Cora, CiteSeer, and Chameleon.
Refer to caption: x1.png

### 4.2 Empirical Studies

We conduct an empirical study to show that the theoretical findings of linear GNNs can also be observed in non-linear GNNs.
For the experiments, we use GCN and GAT and their residual counterparts, resGCN, and resGAT, respectively, as baseline models. To train the models, we use two homophilic datasets, Cora and CiteSeer, and one heterophilic dataset, Chameleon. For each dataset, we train 64- and 128-layer GCN, GAT, resGCN, and resGAT.
We consider three activation functions: ReLU, LeakyReLU with 0.8 negative slope values, and GeLU. We measure the gradient similarity when the test performance is measured through validation.
All experiments are repeated five times, and their average values are reported.
The band indicates min-max values.

#### GNNs without residual connections

Figure 1 illustrates that the gradients of node representations are getting exponentially close to each other over the back-propagation paths for both 128-layer GCN and GAT, independent of the activation functions and datasets. In other words, all nodes are likely to have the same gradient signals near the input layers.
Among three activation functions, GELU leads to the most severe oversmoothing in gradients, followed by ReLU and LeakyReLU.
The experimental results align with the one with node representations but in the other direction (Wu et al., ).
Additionally, we report the results with 64-layer GCN and GAT in Appendix C.

Figure: Figure 2: Gradient similarity measure $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}({\mathbf{W}})}}{ \partial{{{\mathbf{X}}}^{({\ell})}}}}\right)$ over different layers of $64$-layer GCN and GAT with residual connections in three datasets: Cora, CiteSeer, and Chameleon. The similarity measures with ”NaN” value are not indicated in the plot.
Refer to caption: x2.png

#### GNNs with residual connections

We visualize the gradient similarity of 64-layer resGCN and resGAT in Figure 2 with varying activation functions.
The results show that, with the residual connection, the gradients expand regardless of the activations, datasets, and models, while the expansion rate differs depending on the choice of activation, model, and dataset. Note that gradient expansion happens when the magnitude of the gradient is large, and the gradient diverges in directions.
However, among all experiments, we could not find any gradient smoothing, showing that gradient expansion is a more severe problem with residual connections.
In conclusion, these results corroborate the theoretical analysis in Theorem 2.
Additional results with 128-layer resGCN and resGAT are provided in Appendix C.

To explore the training behavior of GNN models, we trained $4$-, $16$-, and $64$-layer GCN, GAT, resGCN, and resGAT models on three datasets: Cora, Citeseer, and Chameleon, with a fixed learning rate of $0.001$ for $1,000$ epochs. Visualizations of the training accuracy and training loss are provided in Appendix D. We observe that both GCN and GAT, with and without residual connections, suffer from more severe underfitting as the number of layers increases. We infer that this phenomenon is related to the presence of gradient oversmoothing and gradient expansion.

## 5 Analysis on Gradient Similarity

We have seen gradients exhibit oversmoothing or expansion based on the model architecture. In this section, we further analyze the relationship between representation similarity and gradient similarity and the changes in the gradient similarity over the training.

### 5.1 Which similarity measure better represents test performance?

Figure: Figure 3: Scatter plots between (a) representation similarity vs gradient similarity (b) test accuracy and representation similarity, (c) test accuracy and gradient similarity on the Cora dataset.
Refer to caption: x3.png

To obtain a deeper insight into the representation and gradient similarities and their implications on model performances, we train GCN, GAT, resGCN, and resGAT with varying numbers of layers on the Cora and Chameleon datasets and observe empirical similarities and model performances. The node similarity is measured from the representation of the last layer, and the gradient similarity is measured from the first layer of the model. All similarities are measured when the test performances are measured through validation.
We plot the empirical relation between 1) representation similarity, 2) gradient similarity, and 3) accuracy on the test set in Figure 3 through scatter plots.

Figure 3(a) shows the relationship between the representation similarity and the gradient similarity. We observe that, in general, there are positive correlations between these two similarities, but they cannot fully explain each other in all regions. For example, when the number of layers is relatively low, the gradient similarity ranges around one while the node similarity ranges between $10^{2}$ and $10^{11}$. When both similarities are close to zero, then the correlation is weak, and we cannot estimate one from the other. It is worth noting that the high gradient similarity can only be observed from the models with residual connections, corroborating the findings in Theorem 2. Similarly, the node representation expansion occurs with the residual connections.

If two similarities are not the same, which better represents the model performance on the test set? In Figure 3(b,c), we plot the relationship between each measure and test accuracy to answer this question.
In Figure 3(b), we observe that when the model achieves relatively good test accuracy, the representation similarity ranges between $10^{2}$ and $10^{11}$.
Since this range is quite broad, cases with low accuracy are also observed in this area. In addition, test accuracy is generally relatively low when the node similarity is too high or too low.
To compare, we find that the range of gradient similarity is narrow when the model performs well in Figure 3(c). Specifically, the gradient similarity ranges between $10^{-2}$ and $1$ when the test performance is relatively high. More importantly, there is no case with relatively low performance in this range, showing that the gradient similarity is a more reliable estimate of the test performance. We provide the results on the CiteSeer and Chameleon dataset in Appendix E, where a similar trend was observed.

### 5.2 How does the gradient similarity change over training?

Figure: Figure 4: Gradient similarity measure $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{GNN}}({\mathbf{W}})}}{\partial{{ {\mathbf{X}}}^{({\ell})}}}}\right)$ over training with 4-, 16-, 64- layered GCN and GAT. We report the similarity measures over two datasets: Cora and Chameleon. The dashed and solid lines represent the similarity measured at the start and end of the training, respectively. The shaded area represents the maximum and minimum similarities over training.
Refer to caption: x4.png

Figure: Figure 5: Gradient similarity measure $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}({\mathbf{W}})}}{ \partial{{{\mathbf{X}}}^{({\ell})}}}}\right)$ over training with 4-, 16-, 64- layered resGCN and resGAT. The dashed and solid lines represent the similarity measured at the start and end of the training, respectively. The shaded area represents the maximum and minimum similarities over training. Residual connections with deep layers introduce gradient expansion.
Refer to caption: x5.png

We investigate how the gradient similarity changes over the training to check whether the optimization can mitigate the gradient expansion over time. With the similar experimental settings used in Section 5.1, we plot how the gradient similarity changes with different models and layers in Figures 4 and 5. We apply early stopping for all experiments.

Figure 4 shows the changes in gradient similarity without residual connections. For all experiments, the similarity increases over training, showing that some level of smoothing happens at the early stage of training. The results show that mild oversmoothing with 16 layers can be overcome through training, but hard oversmoothing with 64 layers cannot be overcome for both models. Figure 5 shows the changes in gradient similarity with the residual connections. The gradient expansion can be observed with 64-layered models. Similar to the oversmoothing, once the expansion happens, the training process cannot be done properly (cf., Figure 3).

## 6 Effectiveness of Liptschitz upper bound on deep GNNs

We propose a Lipschitz upper bound constraint on each layer of GNN through the weight normalization in Section 4.1.
In this section, we show the effectiveness of a Lipschitz upper bound constraint in improving test performance and training behavior in deep GNNs.

To evaluate test performance, we train resGCN and resGAT on Cora, CiteSeer, Chameleon, and Squirrel datasets with varying numbers of layers from $1$ to $512$. For the weight normalization, we test three different Lipschitz upper bound $c$: $4$, $1$, $0.25$. All experiments are conducted with early stopping and learning rates of $0.001,0.005,0.01$. We report average values of five repeated runs.

Figure: Figure 6: Performance of resGCN and resGAT on Cora, CiteSeer, Chameleon and Squirrel datasets with varying coefficient values. The solid curve indicates the average over five runs and the shaded area represents one standard deviation around the average.
Refer to caption: x6.png

The test accuracy of resGCN is reported in Figure 6 with and without weight normalization.
For all datasets, we can observe a significant performance drop when the depth exceeds $32$ without the weight normalization.
Similar patterns are observed when a loose Lipschitz upper bound of $c=4$ is applied.
With the tight upper bounds of $c=1,0.25$, the test accuracy remains relatively consistent over $512$ layers.
Specifically, when $c=0.25$, the accuracy of resGCN consistently improves over a hundred layers, reaching the best performance at a depth of $512$ layers on Squirrel datasets.
However, small upper bounds can limit the expressive power of each layer. As a consequence, test accuracy with a small upper bound requires more layers to reach a similar level of accuracy that can be achieved by a small number of layers without weight normalization.

In the case of homophilic datasets, despite the multiple layers, the highest accuracy is similar regardless of the upper bound constraints. For example, resGCN without normalization achieves the best performance of $85.43$ when the number of layers is $16$. In contrast, resGCN with an upper bound of $0.25$ achieves the best performance of $85.14$ when the number of layers is $128$.
However, in heterophilic datasets, deep models with weight normalization outperform the models without the upper bound constraints.
For instance, in the Chameleon dataset, without normalization, the best accuracy of $43.75$ is achieved at a depth of $32$.
In contrast, the best accuracy of $47.04$ is achieved at a depth of $64$ with an upper bound of $4$.
It is claimed that the heterophilic node classification requires capturing long-range dependencies via multiple layers (Rusch et al., ). Our results also support this claim with carefully designed multi-layered GNNs.

Figure: Figure 7: Training curves with Lipschitz upper bound constraint applied 64-layer resGCN over three datasets: Cora, Citeseer and Chameleon with learning rate $0.001$. The solid line indicates the average over five runs and the shaded area represents maximum and minimum value.
Refer to caption: x7.png

To demonstrate the training process, we train $64$—and $128$—layer resGCN and resGAT on Cora, CiteSeer, and Chameleon datasets with $1,000$ epochs and a fixed learning rate of $0.001$.
For the weight normalization, we test three different Lipschitz upper bound $c$: $4$, $1$, $0.25$.
The visualized training accuracy and loss on $64$-layered resGCN are provided in Figure 7.
We observe that applying the Lipschitz upper bound mitigates the underfitting problem of resGCN in all datasets.
Specifically, resGCN with $c=0.25$ achieves $100\%$ training accuracy in all datasets with the fastest convergence rate.
Our results show that applying the Lipschitz upper bound effectively resolves the gradient expansion problem and helps optimization.
We provide the results for $128$—layer resGCN, and $64$—and $128$—layer resGAT in Appendix F.

## 7 Conclusion

In this work, we explored the challenges associated with optimizing deep GNNs. Until now, the oversmoothing has focused on the representation space.
However, we demonstrate that oversmoothing also occurs in the gradient space, a major cause making the training process complicated.
We further analyze the effect of the residual connections in the gradients of deepGNNs, revealing that the residual connection causes gradient expansion.
To alleviate the gradient expansion during training, we constrain the Lipschitz upper bound on each layer through the Frobenius normalization. Experimental results show that the training procedure can be stabilized with the normalization over several hundred layers.

#### Limitation.

The theoretical analysis is based on linear GCNs. Although we have shown that similar patterns are observable with non-linear GCNs and GATs through empirical studies, the non-linearity can introduce subtle differences compared with our analysis.
We propose a weight normalization to restrict the Lipschitz constant. There are, however, several normalization approaches, such as the batch normalization (Ioffe & Szegedy, ) and the layer normalization (Ba, ). In a previous study, these normalizations are not helpful in deepGNNs (Chen et al., ). To investigate the reason further, we need to study the gradients of these normalizations in GNNs. We leave this future work.

## Appendix A Missing Proofs in Section 4.1

In this section, we provide missing proofs in Section 4.1.

### A.1 Proof of Lemma 1

###### Proof.

We recap the derivation of the gradient of $\mathsf{LGN}$, which is originally shown in Blakely et al. ().
Using the chain rule, we can compute:

$$ $\displaystyle\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{W}}}^ {({\ell})}}$ $\displaystyle=\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{X}}} ^{({\ell+1})}}\frac{\partial{{{\mathbf{X}}}^{({\ell+1})}}}{\partial{{{\mathbf{ W}}}^{({\ell})}}}$ $\displaystyle=\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{X}}} ^{({\ell+1})}}\frac{\partial{}}{\partial{{{\mathbf{W}}}^{({\ell})}}}\left(\hat {{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})}{{\mathbf{W}}}^{({\ell})}\right)$ $\displaystyle=(\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})})^{\top}\frac{ \partial{\mathcal{L}_{\mathsf{LGN}}}}{\partial{{{\mathbf{X}}}^{({\ell+1})}}}\;.$ $$

Again, by using the chain rule, we have:

$$ $\displaystyle\frac{\partial{\mathcal{L}_{\mathsf{LGN}}}}{\partial{{{\mathbf{X} }}^{({\ell})}}}$ $\displaystyle=\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{X}}} ^{({L})}}\frac{\partial{{{\mathbf{X}}}^{({L})}}}{\partial{{{\mathbf{X}}}^{({L- 1})}}}\cdots\frac{\partial{{{\mathbf{X}}}^{({\ell+2})}}}{\partial{{{\mathbf{X} }}^{({\ell+1})}}}\frac{\partial{{{\mathbf{X}}}^{({\ell+1})}}}{\partial{{{ \mathbf{X}}}^{({\ell})}}}$ $\displaystyle=\left(\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{ \mathbf{X}}}^{({L})}}\frac{\partial{{{\mathbf{X}}}^{({L})}}}{\partial{{{ \mathbf{X}}}^{({L-1})}}}\cdots\frac{\partial{{{\mathbf{X}}}^{({\ell+2})}}}{ \partial{{{\mathbf{X}}}^{({\ell+1})}}}\right)\frac{\partial{}}{\partial{{{ \mathbf{X}}}^{({\ell})}}}\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})}{{\mathbf{ W}}}^{({\ell})}$ $\displaystyle=\hat{{\mathbf{A}}}^{\top}\left(\frac{\partial\mathcal{L}_{ \mathsf{LGN}}}{\partial{{\mathbf{X}}}^{({L})}}\frac{\partial{{{\mathbf{X}}}^{( {L})}}}{\partial{{{\mathbf{X}}}^{({L-1})}}}\cdots\frac{\partial{{{\mathbf{X}}} ^{({\ell+2})}}}{\partial{{{\mathbf{X}}}^{({\ell+1})}}}\right)({{\mathbf{W}}}^{ ({\ell})})^{\top}$ $\displaystyle=(\hat{{\mathbf{A}}}^{\top})^{2}\left(\frac{\partial\mathcal{L}_{ \mathsf{LGN}}}{\partial{{\mathbf{X}}}^{({L})}}\frac{\partial{{{\mathbf{X}}}^{( {L})}}}{\partial{{{\mathbf{X}}}^{({L-1})}}}\cdots\frac{\partial{{{\mathbf{X}}} ^{({\ell+3})}}}{\partial{{{\mathbf{X}}}^{({\ell+2})}}}\right)({{\mathbf{W}}}^{ ({\ell+1})})^{\top}({{\mathbf{W}}}^{({\ell})})^{\top}$ $\displaystyle=\cdots$ $\displaystyle=(\hat{{\mathbf{A}}}^{\top})^{L-\ell}\frac{\partial\mathcal{L}_{ \mathsf{LGN}}}{\partial{{\mathbf{X}}}^{({L})}}\prod_{i=1}^{L-\ell}({\mathbf{W} }^{(L-i)})^{\top}\;.$ $$

∎

###### Proof.

### A.2 Proof of Theorem 1

###### Proof.

To prove Theorem 1, we recap and use the following definitions and lemmas from Wu et al. ().

###### Definition 1 (Ergodicity) .

Let ${\mathbf{B}}\in\mathbb{R}^{(N-1)\times N}$ be the orthogonal projection onto the space orthogonal to $\mathsf{span}\{\mathbf{1}\}$. A sequence of matrices $\{M^{(n)}\}_{n=1}^{\infty}$ is ergodic if

$$ $\lim_{t\to\infty}{\mathbf{B}}\prod_{n=0}^{t}{\mathbf{M}}^{(n)}=0.$ $$

By using the projection matrix ${\mathbf{B}}$, we can express $\mu({\mathbf{X}})$ as $\lVert{{\mathbf{B}}{\mathbf{X}}}\rVert_{F}$.
We want to show that $\{\hat{{\mathbf{A}}}\}_{n=1}^{\infty}$ is ergodic.
Since $\hat{{\mathbf{A}}}=\tilde{{\mathbf{D}}}^{-\frac{1}{2}}\tilde{{\mathbf{A}}}
\tilde{{\mathbf{D}}}^{-\frac{1}{2}}$ and $\tilde{{\mathbf{D}}}^{-1}\tilde{{\mathbf{A}}}$ have a same spectrum, $\{\hat{{\mathbf{A}}}\}_{n=1}^{\infty}$ is ergodic if $\{\tilde{{\mathbf{D}}}^{-1}\tilde{{\mathbf{A}}}\}_{n=1}^{\infty}$ is ergodic.
Using Lemma 3 in (Wu et al., ), $\{\tilde{{\mathbf{D}}}^{-1}\tilde{{\mathbf{A}}}\}_{n=1}^{\infty}$ is ergodic.

Next, we need notion of joint spectral radius to show the convergence rate.

###### Definition 2 (Joint Spectral Radius) .

For a collection of matrices $\mathcal{A}$, the joint spectral radius $\operatorname{JSR}(\mathcal{A})$ is defined to be

$$ $\operatorname{JSR}(\mathcal{A})=\limsup_{k\rightarrow\infty}\sup_{{\mathbf{A}} _{1},{\mathbf{A}}_{2},\ldots,{\mathbf{A}}_{k}\in\mathcal{A}}\left\|{\mathbf{A} }_{1}{\mathbf{A}}_{2}\cdots{\mathbf{A}}_{k}\right\|^{\frac{1}{k}}$ $$

and it is independent of the norm used.

By the definition, it is straight forward that if $\operatorname{JSR}(\mathcal{A})<1$, for any $\operatorname{JSR}(\mathcal{A})<q<1$, there exists a $C$ for which satisfies $||{\mathbf{A}}_{1}{\mathbf{A}}_{2}\cdots{\mathbf{A}}_{k}\mathbf{y}||\leq Cq^{k
}||\mathbf{y}||$ for ${\mathbf{A}}_{1},{\mathbf{A}}_{2},\cdots,{\mathbf{A}}_{k}\in\mathcal{A}$.
Let $\hat{\mathcal{A}}=\{\hat{{\mathbf{A}}}\}$ and $\tilde{\mathcal{A}}=\{\tilde{{\mathbf{D}}}^{-1}\tilde{{\mathbf{A}}}\}$.
Similar to Ergodicity, we can get $\operatorname{JSR}(\hat{\mathcal{A}})<1$ if $\operatorname{JSR}(\tilde{\mathcal{A}})<1$.
Using Lemma 6 in (Wu et al., ), $\operatorname{JSR}(\tilde{\mathcal{A}})<1$.

Now, we can derive:

$$ $\displaystyle\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{LGN}}}}{\partial{{{ \mathbf{X}}}^{({\ell})}}}}\right)$ $\displaystyle=\lVert{{\mathbf{B}}\frac{\partial{\mathcal{L}_{\mathsf{LGN}}}}{ \partial{{{\mathbf{X}}}^{({\ell})}}}}\rVert_{F}$ $\displaystyle=\lVert{{\mathbf{B}}(\hat{{\mathbf{A}}}^{\top})^{L-\ell}\frac{ \partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{\mathbf{X}}}^{({L})}}\prod_{i=1} ^{L-\ell}({\mathbf{W}}^{(L-i)})^{\top}}\rVert_{F}$ $\displaystyle\leq\lVert{\frac{\partial\mathcal{L}_{\mathsf{LGN}}}{\partial{{ \mathbf{X}}}^{({L})}}}\rVert_{F}\lVert{{\mathbf{B}}(\hat{{\mathbf{A}}}^{\top}) ^{L-\ell}}\rVert_{2}\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}^{L-\ell}\;,$ $$

where $\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}$ is the maximum spectral norm, i.e., $*=\operatorname*{arg\,max}_{i}\lVert{{\mathbf{W}}^{(i)}}\rVert_{2}$ for $1\leq i\leq L-\ell$.
By using $\operatorname{JSR}(\hat{\mathcal{A}})<1$, we can result in to Theorem 1.

∎

###### Proof.

###### Definition 1 (Ergodicity) .

###### Definition 2 (Joint Spectral Radius) .

### A.3 Proof of Lemma 2

###### Proof.

By applying the chain rule, we can obtain:

$$ $\displaystyle\frac{\partial\mathcal{L}_{\mathsf{resLGN}}}{\partial{{\mathbf{W} }}^{({\ell})}}$ $\displaystyle=\frac{\partial\mathcal{L}_{\mathsf{resLGN}}}{\partial{{\mathbf{X }}}^{({\ell+1})}}\frac{\partial{{{\mathbf{X}}}^{({\ell+1})}}}{\partial{{{ \mathbf{W}}}^{({\ell})}}}$ $\displaystyle=\frac{\partial\mathcal{L}_{\mathsf{resLGN}}}{\partial{{\mathbf{X }}}^{({\ell+1})}}\frac{\partial{}}{\partial{{{\mathbf{W}}}^{({\ell})}}}\left({ {\mathbf{X}}}^{({\ell})}+\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})}{{\mathbf{ W}}}^{({\ell})}\right)$ $\displaystyle=\left(\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})}\right)^{\top} \frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({\ell +1})}}}\;.$ $$

We have:

$$ $\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({\ell })}}}=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{ ({L})}}}\frac{\partial{{{\mathbf{X}}}^{({L})}}}{\partial{{{\mathbf{X}}}^{({L-1 })}}}\frac{\partial{{{\mathbf{X}}}^{({L-1})}}}{\partial{{{\mathbf{X}}}^{({L-2} )}}}\cdots\frac{\partial{{{\mathbf{X}}}^{({\ell+1})}}}{\partial{{{\mathbf{X}}} ^{({\ell})}}}\;.$ (7) $$

We will show Lemma 2 holds inductively.
First two components of Equation 7 can be derived as follows:

$$ $\displaystyle\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf {X}}}^{({L})}}}\frac{\partial{{{\mathbf{X}}}^{({L})}}}{\partial{{{\mathbf{X}}} ^{({L-1})}}}$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}\frac{\partial{}}{\partial{{{\mathbf{X}}}^{({L-1})}}}({{ \mathbf{X}}}^{({L-1})}+\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({L-1})}{{\mathbf{W}} }^{({L-1})})$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}+(\hat{{\mathbf{A}}})^{\top}\frac{\partial{\mathcal{L}_{ \mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}({{\mathbf{W}}}^{({L-1})}) ^{\top}\;,$ (8) $$

which satisfies Lemma 2 with $\ell=L-1$.
Suppose we have $\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({\ell
+1})}}}$ which satisfies Lemma 2.
We can derive:

$$ $\displaystyle\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf {X}}}^{({\ell})}}}$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({\ell+1})}}}\frac{\partial{{{\mathbf{X}}}^{({\ell+1})}}}{ \partial{{{\mathbf{X}}}^{({\ell})}}}$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({\ell+1})}}}\frac{\partial{}}{\partial{{{\mathbf{X}}}^{({\ell}) }}}({{\mathbf{X}}}^{({\ell})}+\hat{{\mathbf{A}}}{{\mathbf{X}}}^{({\ell})}{{ \mathbf{W}}}^{({\ell})})$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({\ell+1})}}}+\hat{{\mathbf{A}}}^{\top}\frac{\partial{\mathcal{L }_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({\ell+1})}}}({{\mathbf{W}}}^{( {\ell})})^{\top}$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}+\hat{{\mathbf{A}}}^{\top}\frac{\partial{\mathcal{L}_{ \mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}({{\mathbf{W}}}^{({\ell})} )^{\top}$ $\displaystyle+\sum_{p=1}^{L-\ell-1}\left(\sum_{\{i_{1}:i_{1}\geq\ell+1\}} \cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p} \frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})} }}\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right)$ $\displaystyle+\sum_{p=1}^{L-\ell-1}\left(\sum_{\{i_{1}:i_{1}\geq\ell+1\}} \cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p+ 1}\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L} )}}}\left(\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right) ({{\mathbf{W}}}^{({\ell})})^{\top}\right)$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}+\hat{{\mathbf{A}}}^{\top}\frac{\partial{\mathcal{L}_{ \mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}({{\mathbf{W}}}^{({\ell})} )^{\top}+\sum_{\ell+1\leq i_{1}<L}\hat{{\mathbf{A}}}^{\top}\frac{\partial{ \mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}({{\mathbf{W} }}^{({i_{1}})})^{\top}$ $\displaystyle+\sum_{p=2}^{L-\ell-1}\left(\sum_{\{i_{1}:i_{1}\geq\ell+1\}} \cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p} \frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})} }}\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right)$ $\displaystyle+\sum_{p=1}^{L-\ell-2}\left(\sum_{\{i_{1}:i_{1}\geq\ell+1\}} \cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p+ 1}\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L} )}}}\left(\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right) ({{\mathbf{W}}}^{({\ell})})^{\top}\right)$ $\displaystyle+\left(\hat{\mathbf{A}}^{\top}\right)^{L-\ell}\frac{\partial{ \mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}\left(\prod_{ i=1}^{L-\ell}({\mathbf{W}}^{(L-i)})^{\top}\right)\;.$ $$

The fourth term and the fifth term can be merged into:

$$ $\sum_{p=2}^{L-\ell-1}\left(\sum_{\{i_{1}:i_{1}\geq\ell\}}\cdots\sum_{\{i_{p}:i _{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p}\frac{\partial{ \mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}\prod_{k=0}^{ p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right)\;,$ $$

since $\binom{L-\ell-1}{p+1}+\binom{L-\ell-1}{p}=\binom{L-\ell}{p+1}$.
Therefore,

$$ $\displaystyle\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf {X}}}^{({\ell})}}}$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}+\sum_{\ell\leq i_{1}<L}\hat{{\mathbf{A}}}^{\top}\frac{ \partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}({{ \mathbf{W}}}^{({i_{1}})})^{\top}$ $\displaystyle+\sum_{p=2}^{L-\ell-1}\left(\sum_{\{i_{1}:i_{1}\geq\ell\}}\cdots \sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p}\frac{ \partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}} \prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right)$ $\displaystyle+\left(\hat{\mathbf{A}}^{\top}\right)^{L-\ell}\frac{\partial{ \mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}\left(\prod_{ i=1}^{L-\ell}({\mathbf{W}}^{(L-i)})^{\top}\right)$ $\displaystyle=\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{ \mathbf{X}}}^{({L})}}}+\sum_{p=1}^{L-\ell}\left(\sum_{\{i_{1}:i_{1}\geq\ell\}} \cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L\}}\left(\hat{\mathbf{A}}^{\top}\right)^{p} \frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})} }}\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}\right)\;,$ $$

which satisfies Lemma 2.
∎

###### Proof.

### A.4 Proof of Theorem 2

###### Proof.

We can obtain:

$$ $\displaystyle\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}({\mathbf{W }})}}{\partial{{{\mathbf{X}}}^{({\ell})}}}}\right)$ $\displaystyle=\lVert{{\mathbf{B}}\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}( {\mathbf{W}})}}{\partial{{{\mathbf{X}}}^{({\ell})}}}}\rVert_{F}$ $\displaystyle\leq\lVert{{\mathbf{B}}\frac{\partial{\mathcal{L}_{\mathsf{resLGN }}}}{\partial{{{\mathbf{X}}}^{({L})}}}}\rVert_{F}$ $\displaystyle+\sum_{p=1}^{L-\ell}\sum_{\{i_{1}:i_{1}\geq\ell\}}\cdots\sum_{\{i _{p}:i_{p-1}<i_{p}<L\}}\lVert{{\mathbf{B}}\left(\hat{\mathbf{A}}^{\top}\right) ^{p}\frac{\partial{\mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({ L})}}}\prod_{k=0}^{p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}}\rVert_{F}\;.$ $$

Since there are $\binom{L-\ell}{p}$ sets of $(i_{1},...i_{p})$ which satisfies $\ell\leq i_{1}<i_{2}<\cdots<i_{p}<L$,

$$ $\displaystyle\sum_{\{i_{1}:i_{1}\geq\ell\}}\cdots\sum_{\{i_{p}:i_{p-1}<i_{p}<L \}}\lVert{{\mathbf{B}}\left(\hat{\mathbf{A}}^{\top}\right)^{p}\frac{\partial{ \mathcal{L}_{\mathsf{resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}\prod_{k=0}^{ p-1}\left({\mathbf{W}}^{(i_{p-k})}\right)^{\top}}\rVert_{F}$ $\displaystyle\leq\binom{L-\ell}{p}\lVert{\frac{\partial\mathcal{L}_{\mathsf{ LGN}}}{\partial{{\mathbf{X}}}^{({L})}}}\rVert_{F}\lVert{{\mathbf{B}}(\hat{{ \mathbf{A}}}^{\top})^{L-\ell}}\rVert_{2}\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}^{ L-\ell}\;,$ $$

where $\lVert{{\mathbf{W}}^{(*)}}\rVert_{2}$ is the maximum spectral norm, i.e., $*=\operatorname*{arg\,max}_{i}\lVert{{\mathbf{W}}^{(i)}}\rVert_{2}$ for $1\leq i\leq L-\ell$.
By using $\operatorname{JSR}(\hat{\mathcal{A}})<1$, there exists $0<q<1$ and $C_{q}>0$ such that

$$ $\displaystyle\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}({\mathbf{W }})}}{\partial{{{\mathbf{X}}}^{({\ell})}}}}\right)$ $\displaystyle\leq\mu\left({{\mathbf{B}}\frac{\partial{\mathcal{L}_{\mathsf{ resLGN}}}}{\partial{{{\mathbf{X}}}^{({L})}}}}\right)+\sum_{p=1}^{L-\ell}C_{q} \binom{L-\ell}{p}(q\lVert{{\mathbf{W}}^{(*)}}\rVert_{2})^{L-\ell}$ $\displaystyle=\mu\left({{\mathbf{B}}\frac{\partial{\mathcal{L}_{\mathsf{resLGN }}}}{\partial{{{\mathbf{X}}}^{({L})}}}}\right)+C_{q}\left(\left(1+q\lVert{{ \mathbf{W}}^{(*)}}\rVert_{2}\right)^{L-\ell}-1\right)\;.$ $$

∎

###### Proof.

## Appendix B Dataset Statistics

**Table 1: Statistics of the dataset utilized in the experiments.**
| Dataset | Type | # nodes | # edges | # classes | avg degree |
| --- | --- | --- | --- | --- | --- |
| Cora | homophilic | 2,708 | 5,278 | 7 | 3.90 |
| CiteSeer | homophilic | 3,327 | 4,552 | 6 | 2.74 |
| Squirrel | heterophilic | 2,223 | 46,998 | 5 | 42.28 |
| Chameleon | heterophilic | 890 | 8,854 | 5 | 19.90 |

Statistics of datasets used in our experiments can be found in Table 1.
We note that the filtering process proposed by Platonov et al. () was adopted in the case of Chameleon and Squirrel datasets.

## Appendix C Supplementary Visualizations of Gradient Similarity Evaluations

Supplementary visualizations for the analysis in Section 4.2 are provided below. We did not plot gradient similarity measure of $128$-layer GAT, as the measures yielded ”NaN” values in all layers except on last layer because of expansion.

Figure: Figure 8: Gradient similarity measure $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{GNN}}({\mathbf{W}})}}{\partial{{ {\mathbf{X}}}^{({\ell})}}}}\right)$ over different layers and activation functions. We report the oversmoothing measures of $64$-layer GCN and GAT over three datasets: Cora, CiteSeer, and Chameleon.
Refer to caption: x8.png

Figure: Figure 9: Gradient similarity measure $\mu\left({\frac{\partial{\mathcal{L}_{\mathsf{resGNN}}({\mathbf{W}})}}{ \partial{{{\mathbf{X}}}^{({\ell})}}}}\right)$ over different layers of $128$-layer GCN with residual connections.
Refer to caption: x9.png

## Appendix D Visualizations of Training Curves

Supplementary visualizations for the $4$-, $16$-, $64$- layer GCN, GAT, resGCN and resGAT are provided below.

Figure: Figure 10: Training curves $4$-, $16$-, $64$- layer GCN and resGCN over three datasets: Cora, Citeseer, and Chameleon with learning rate $0.001$. The dashed lines represent the existence of residual connection while solid lines represent the vanilla model.
Refer to caption: x10.png

Figure: Figure 11: Training curves $4$-, $16$-, $64$- layer GAT and resGAT over three datasets: Cora, Citeseer, and Chameleon with learning rate $0.001$. The dashed lines represent the existence of residual connection, while solid lines represent the vanilla model.
Refer to caption: x11.png

## Appendix E Supplementary Scatter Plots

Scatter plots for the analysis in Section 5.1 are provided below.

Figure: Figure 12: Scatter plots between (a) representation similarity vs gradient similarity (b) test accuracy and representation similarity, (c) test accuracy and gradient similarity on the CiteSeer dataset.
Refer to caption: x12.png

Figure: Figure 13: Scatter plots between (a) representation similarity vs gradient similarity (b) test accuracy and representation similarity, (c) test accuracy and gradient similarity on the Chameleon dataset.
Refer to caption: x13.png

## Appendix F Supplementary Visualizations of Training Curves with Lipschitz Upper Bound Constraint

Supplementary visualizations of training curves with Lipschitz upper bound constraint for Section 6 are provided below.

Figure: Figure 14: Training curves with Lipschitz upper bound constraint applied $128$-layer resGCN over three datasets: Cora, Citeseer and Chameleon with learning rate $0.001$.
Refer to caption: x14.png

Figure: Figure 15: Training curves with Lipschitz upper bound constraint applied $64$-layer resGAT over three datasets: Cora, Citeseer and Chameleon with learning rate $0.001$.
Refer to caption: x15.png

Figure: Figure 16: Training curves with Lipschitz upper bound constraint applied $128$-layer resGAT over three datasets: Cora, Citeseer and Chameleon with learning rate $0.001$.
Refer to caption: x16.png

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 3
presentation: 1
questions: Can you try some of the experiments with GNN architectures that are known to not suffer from gradient vanishing/explosion or oversmoothing? If your gradient similarity measure would go to zero or explode for these cases, it would strengthen your claim.
soundness: 2
strengths: The paper studies gradient behavior during training of deep GNNs. While this is a widely studied field, mainly linked to oversquashing these days, its full qualitative behavior remains to be studied in more detail. The paper is further a good combination of theoretical and empirical work.
summary: The paper argues that gradient oversmoothing and gradient expansion pose challenges in training deep Graph Neural Networks (GNNs). Unlike previous approaches that focus on node features, the authors apply a previously defined similarity measure to gradients. They also establish an asymptotic upper bound for this measure, describing the phenomenon where it approaches zero as "gradient oversmoothing," and linking its growth to "gradient expansion." To address these issues, the authors propose a novel normalization technique. Empirical results on various graph datasets and architectures demonstrate that this measure decreases as layer depth increases, and the new normalization approach enables successful training of deep GNNs.
weaknesses: My biggest concern is that the statement on gradient oversmoothing in the asymptotic case is likely to be equivalent to gradient vanishing, i.e., that gradients vanish whenever the proposed similarity measure vanishes and vice versa. In fact, very similar statements to Theorem 1 (i.e., bounds depending on the weight matrix W) can be made for the gradients directly (see Di Giovanni et al, 2024 for instance). Of course it can happen that in the pre-asymptotic case this gradient similarity measure is small or even zero at some point, but this does not denote an issue. 

This has to be addressed, otherwise this paper cannot be published. In particular, the authors have to show that there exists an asymptotic regime, where the similarity measure in Theorem 1 goes to zero, but the gradient norms do not. I suspect there are no such cases, but I am happy to be convinced otherwise. The same has to be done for Theorem 2.

That being said, this has also be demonstrated empirically. The provided plots in the paper only show the effect of varying the depth on the similarity measure. The same has to be done for the gradient norms for the exact same setup, i.e., same architecture with same weights. This would demonstrate that there are in fact cases, where gradients oversmooth but do not vanish.

Other issues:
* Some mathematical statements are wrong. For instance, the sentence after equation (3). Here, if B would be perpendicular to span($\bf 1$), it should be in $\mathbb{R}^{N-1,N}$, not in $\mathbb{R}^{N-1,N-1}$. Moreover, this sentence does not add anything to the context and is simply taken from Wu et al. which used it to demonstrate the similarity between this measure and another oversmoothing measure in the literature. 
* most of the theory seems to be taken from other papers. It would be good if the authors could state explicitly what is their own theoretical contribution
* it would be good to extend to nonlinear layers, since in practice linear GNNs are not very common
* Only four small-scale graph datasets are considered. It would strengthen the paper if this would be extended.

Minor issues:
* The writing needs to be revised. There are many typos.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: Please see my suggestions in the weaknesses section.
soundness: 3
strengths: In terms of strengths, the paper is nicely written and easy to follow. It also sheds more light on the oversmoothing problem. The experiments conducted to understand the theoretical findings are in-order.
summary: This paper studies the behavior of gradients in deep GNNs by considering the oversmoothing problem. The authors show that when GNNs are very deep, the representations become similar (which is a known thing), and that the gradients in first layers become very similar (makes sense by inspecting the Jacobian of the GNN).

The authors perform several experiments, mostly on simple node classification datasets, and show the training/test accuracies alongside the gradient and feature similarities.
weaknesses: In terms of relevance and novelty, as well as related works, I think that the main issue is that most of the presented results are already quite known within the community, and can be found in the literature. For example see "A survey on oversmoothing in graph neural networks"  and "Simplifying the Theory on Over-Smoothing". Also, the authors lack a discussion of "Revisiting Graph Neural Networks: All We Have is Low-Pass Filters".

In terms of using the Lipschitz constant and bounding it, it was shown in "On the Robustness of Graph Neural Diffusion to Topology Perturbations" and "Contractive Systems Improve Graph Neural Networks Against Adversarial Attacks" that using it can help to address robustness problems, so it would be interesting to understand the connection between oversmoothing and such approaches. 

In terms of experiments, the authors used mostly simple datasets, and I think that it would be beneficial to also study the performance on other tasks and datasets. I think it would also be interesting to see a report of the Dirichlet energy, that is usually reported in oversmoothing studies. Also, I feel that the experiments lack comparisons with other methods that can address oversmoothing (there are already many such methods).

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 4
questions: - As the graph size (number of nodes specifically) diverges, what happens to the proposed normalization?  

- Do you not have any theoretical guarantees for the proposed procedure since you could not obtain matching lower bounds?  I assume that would be required to argue that the normalization does the trick.

- Can you add a few more words explaining what the trouble is with bipartide?
soundness: 3
strengths: The numerical experitments are very nicely executed, clear, convincing, and reproducible (but I still have questions as they do not seem to match the maths).  The bounds are reasonable (even if the proof is extremely straightforward....) and interesting.
summary: The authors analyze the effects of the addition of  residual connections in **linear** GNNs **with no non-linear activation function** (specifically GCNNs) as an effort to compact oversmoothing (in extremely deep GNNs).  In a pair of two results (Theorems 1 and 2) the authors analyize the decay of a (reasonable) orthogonal projection of the gradient of the loss with respect to pertrubations of entires of the feature matrix.  

The use their theoretical upper bounds to suggest that spectral weight normalization, in the spirit of [1], would remedy oversmoothing.  Only empirical studies are then used to support at (presumably) theoretical claim.



[1] Neyshabur, Behnam, Russ R. Salakhutdinov, and Nati Srebro. "Path-sgd: Path-normalized optimization in deep neural networks." Advances in neural information processing systems 28 (2015).
weaknesses: - **Identity Activation Only**: Only identity activation are considered (it doesn't seem to be difficult to obtain a result for non-linear activation fixing the origin from your proofs....

- **Connected Graphs** Why not consider disconnected graphs?  It seems to me that with only a bit more work, you can generalize your result by noting that the normalized adjacency matrix has a block-diagonal form in the general (non-bipartide case).


- Results are for non-activated GNNs (identity activation function) but the numerics have non-linear activations (and result with identity activation are not plotted).  This latter point, namely that there are no numerical experiments matching the theoretical setting, makes me wonder if there is a gap...Does the theory really come through in practice.  One **needs** illustrations matching the architectural setup in the theorem.

- Theory is for very basic GCNN model (Which is definitely not universal by any means) but many other GNN models are considered in the experiments section.  What is the relevance?

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Refer to weaknesses
soundness: 3
strengths: 1. The paper presents a comprehensive theoretical framework for analyzing the over-smoothing phenomenon in GNNs and gradient expansion in residual GNNs, introducing the concepts of representation similarity and gradient similarity to further explain these phenomena.

2. The mathematical derivations are rigorous yet easy to comprehend.

3. The author extensively investigates gradient similarity and representation similarity through numerous experiments, examining their variation throughout the training process.

4. The paper proposes a Lipschitz constraint method that effectively alleviates the effects of over-smoothing and gradient expansion during the training process.
summary: The paper analyzed the over-smoothing phenomenon in Graph Neural Networks (GNNs) from a unique perspective, uncovering the issue of gradient expansion in residual connections. To address the problem of gradient explosion, the author introduced a layer-wise Lipschitz constraint, which facilitates efficient training of residual GNNs and enhances their performance.
weaknesses: 1. In line 147, the equation $$ \||\mathbf{X}\||_F^2 = \mu(\mathbf{X})^2 + \||\mathbf{X} - \mathbf{B}\mathbf{X}\||_F^2 $$ is presented. a step-by-step derivation or proof of this equation in an appendix, which will greatly enhance the clarity of the methodology presented.

2. The authors provided a mathematical derivation for representation similarity and gradient similarity, validating their correctness in Figures 1 and 2. However, it is unclear which metric better measures over-smoothing, as the relationships among representation similarity, gradient similarity, and test accuracy appear ambiguous in Figure 3. It seems the authors primarily used gradient similarity to assess over-smoothing in Section 5. In references [1] and [2], Dirichlet energy appears to be a more effective metric compared to gradient similarity.  It would be helpful if  authors could provide a more detailed comparison between gradient similarity and other metrics like Dirichlet energy, including quantitative results, which would help clarify the advantages of their chosen metric.

3. In my view, large datasets can mitigate oversmoothing compared to smaller ones. Intuitively, when aggregating features, each node has more distant neighbors, making it harder for features to become overly similar. As shown in [2], Pubmed and Ogbn-arxiv are less affected, as seen in Table 1. Could authors extend their theoretical analysis to consider the impact of dataset size on oversmoothing to enhance the theoretical persuasiveness of the paper?

4. The Lipschitz constraint method is an effective approach to alleviating over-smoothing and gradient expansion. It would enhance the paper's strength to incorporate 2-3 specific baseline methods from the related work section and utilize 1-2 large datasets from references [1] and [2] that are particularly relevant for comparison.

5. It appears that the theoretical derivation in the paper is based on the GCN model, yet the authors have also conducted some experiments on GAT, such as in Figure 4. Perhaps the authors should provide a corresponding theoretical analysis to enhance the persuasiveness of the paper. Additionally, if the model is changed to another one, such as GraphSage, would the theory still be useful?

[1] A Survey on Oversmoothing in Graph Neural Networks  
[2] Dirichlet Energy Constrained Learning for Deep Graph Neural Networks