On Addressing the Limitations of Graph Neural Networks
This report gives a summary of two problems about graph convolutional
networks (GCNs): over-smoothing and heterophily challenges, and outlines future
directions to explore.
Sitao Luan
# On Addressing the Limitations of Graph Neural Networks

###### Abstract

This report gives a comprehensive summary of two problems about graph convolutional networks (GCNs): over-smoothing and heterophily challenges, and outlines future directions to explore.

## 1 Introduction

Many real-world problems can be modeled as graphs. Recently, neural network based approaches have achieved significant progress for solving large, complex, graph-structured problems [9, 12, 17, 21, 24, 50, 12]. Inspired by the success of Convolutional Neural Networks (CNNs) [31] in computer vision [33], graph convolution defined on the graph Fourier domain stands out as the key operator and one of the most powerful tools for using machine learning to solve graph problems. Although with high expressive power, GCNs still suffer from several difficulties, _e.g._ the over-smoothing problem limits deep GCNs to sufficiently exploit multi-scale information, heterophily problem makes the graph-aware models underperform the graph-agnostic models. This report summarizes the methods we have proposed to address those challenges and puts forward some research problems we will investigate.

To fully explain the above problems, in subsection 1.1, we will first introduce the notation and background knowledge of graph networks. In section 2, we introduce the loss of expressive power of deep graph neural networks (GNNs) and propose snowball and truncated Krylov architecture to address it; in section 3, we analyze heterophily problems for the existing GNNs and propose Adaptive Channel Mixing (ACM) architecture to address it.

Main ContributionIn section 2, we first point out that the output of deep GCN with ReLU activation function will suffer from loss of rank problem under certain conditions and this can cause deep GCN lose expressive power. We then prove that Tanh is better at preserving the rank of the output and verify this claim with numerical tests. Then we find a way to deepen GCN in block Krylov form and propose snowball and truncated Krylov networks which perform better than state-of-the-arts (SOTA) model on semi-supervised node classification tasks on 3 benchmark
datasets. Besides, we point out that finding **a specifically tailored weight initialization scheme for GCNs** can be an promising direction to address over-smoothing efficiently in section 2.2. In section 3, we first illustrate the insufficiency of the current homophily metrics and propose aggregation homophily based on a new similarity matrix. We then show the advantage of the new homophily metric over the existing ones on synthetic graph. Based on the similarity matrix, we define diversification distinguishability of a node and demonstrate why high-pass filters can help to address heterophily problem. To include both low-pass and high-pass in GNNs, we extend filterbank method and propose ACM and ACMII frameworks that can boost the performance of baseline GNNs on heterophilous graphs.

### Notation and Background Knowledge

Suppose we have an undirected graph \(\mathcal{G}=(\mathcal{V},\mathcal{E},A)\), where \(\mathcal{V}\) is the node set with \(|\mathcal{V}|=N\); \(\mathcal{E}\) is the edge set without self-loop; \(A\in\mathbb{R}^{N\times N}\) is the symmetric adjacency matrix with \(A_{ij}=1\) if and only if \(e_{ij}\in\mathcal{E}\), otherwise \(A_{ij}=0\); \(D\) is the diagonal degree matrix, _i.e._\(D_{ii}=\sum_{j}A_{ij}\) and \(\mathcal{N}_{i}=\{j:e_{ij}\in\mathcal{E}\}\) is the neighborhood set of node \(i\). A graph signal is a vector \(\mathbf{x}\in\mathbb{R}^{N}\) defined on \(\mathcal{V}\), where \(x_{i}\) is defined on the node \(i\). We also have a feature matrix \(X\in\mathbb{R}^{N\times F}\) whose columns are graph signals and each node \(i\) has a corresponding feature vector \(X_{i:}\) with dimension \(F\), which is the \(i\)-th row of \(X\). We denote \(Z\in\mathbb{R}^{N\times C}\) as label encoding matrix, where \(Z_{i:}\) is the one hot encoding of the label of node \(i\) and \(C\) is the total number of classes.

The (combinatorial) graph Laplacian is defined as \(L=D-A\), which is a Symmetric Positive Semi-Definite (SPSD) matrix [7]. Its eigendecomposition gives \(L=U\Lambda U^{T}\), where the columns of \(U\in\mathbb{R}^{N\times N}\) are orthonormal eigenvectors, namely the _graph Fourier basis_, \(\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{N})\) with \(\lambda_{1}\leq\cdots\leq\lambda_{N}\), and these eigenvalues are also called _frequencies_. The graph Fourier transform of the graph signal \(\mathbf{x}\) is defined as \(\mathbf{x}_{\mathcal{F}}=U^{-1}\mathbf{x}=U^{T}\mathbf{x}=[\mathbf{u}_{1}^{T}\mathbf{x},\ldots, \mathbf{u}_{N}^{T}\mathbf{x}]^{T}\), where \(\mathbf{u}_{i}^{T}\mathbf{x}\) is the component of \(\mathbf{x}\) in the direction of \(\mathbf{u}_{i}\).

Some graph Laplacian variants are commonly used, _e.g._ the symmetric normalized Laplacian \(L_{\text{sym}}=D^{-1/2}LD^{-1/2}=I-D^{-1/2}AD^{-1/2}\) and the random walk normalized Laplacian \(L_{\text{rw}}=D^{-1}L=I-D^{-1}A\). The eigenvalues of \(L_{\text{rw}}\) and \(L_{\text{sym}}\) are the same and are in \([0,2)\), and their corresponding eigenvectors satisfy \(\mathbf{u}_{\text{rw}}^{i}=D^{-1/2}\mathbf{u}_{\text{sym}}^{i}\).

The affinity (transition) matrices can be derived from the Laplacians, _e.g._\(A_{\text{rw}}=I-L_{\text{rw}}=D^{-1}A\), \(A_{\text{sym}}=I-L_{\text{sym}}=D^{-1/2}AD^{-1/2}\). Then \(\lambda_{i}(A_{\text{rw}})=\lambda_{i}(A_{\text{sym}})=1-\lambda_{i}(A_{\text {sym}})=1-\lambda_{i}(A_{\text{rw}})\in(-1,1]\). Renormalized affinity and Laplacian matrices are introduced in [24] as \(\widehat{A}_{\text{sym}}=\widetilde{D}^{-1/2}\tilde{A}\widetilde{D}^{-1/2}\), \(\hat{L}_{\text{sym}}=I-\hat{A}_{\text{sym}}\), where \(\tilde{A}\equiv A+I,\widetilde{D}\equiv D+I\) and it essentially adds a self-loop and is widely used in Graph Convolutional Network (GCN) as follows:

\[Y=\text{softmax}(\hat{A}_{\text{sym}}\text{ ReLU}(\hat{A}_{\text{sym}}XW_{0})\ W_{1}) \tag{1}\]

where \(W_{0}\in\mathbb{R}^{F\times F_{1}}\) and \(W_{1}\in\mathbb{R}^{F_{1}\times O}\) are parameter matrices. GCN can learn by minimizing the following cross entropy loss

\[\mathcal{L}=-\text{trace}(Z^{T}\log Y). \tag{2}\]

The random walk renormalized matrix \(\hat{A}_{\text{rw}}=\tilde{D}^{-1}\tilde{A}\) can also be applied to GCN and it has the same eigenvalues as \(\hat{A}_{\text{sym}}\). The corresponding Laplacian is defined as \(\hat{L}_{\text{rw}}=I-\hat{A}_{\text{rw}}\) Specifically, the nature of random walk matrix makes \(\hat{A}_{\text{rw}}\) behaves as a mean aggregator
\(\sum_{j\in\{\mathcal{N}_{i}\cup i\}}x_{j}/(D_{ii}+1)\) which is applied in [17] and is important to bridge the gap between spatial- and spectral-based graph convolution methods.

## 2 Loss of Expressive Power of Deep Graph Neural Networks

One major problem of the existing GCNs is the low expressive power limited by their shallow learning mechanisms [66, 61]. There are mainly two reasons why an architecture that is scalable in depth has not been achieved yet. First, this problem is difficult: considering graph convolution as a special form of Laplacian smoothing [32], networks with multiple convolutional layers will suffer from an over-smoothing problem that makes the representation of even distant nodes indistinguishable [66]. Second, some people think it is unnecessary: for example, [4] states that it is not necessary for the label information to totally traverse the entire graph and one can operate on the multi-scale coarsened input graph and obtain the same flow of information as GCNs with more layers. Acknowledging the difficulty, we hold on to the objective of deepening GCNs since the desired compositionality1 will yield easy articulation and consistent performance for problems with different scales.

Footnote 1: The expressive power of a sound deep Neural Network (NN) architecture should be expected to grow with the increment of network depth [19, 30].

In subsection 2.1, we first analyze the limits of deep GCNs brought by over-smoothing and the activation functions. Then, we show that any graph convolution with a well-defined analytic spectral filter can be written as a product of a block Krylov matrix and a learnable parameter matrix in a special form. Based on this, we propose two GCN architectures that leverage multi-scale information in different ways and are scalable in depth, with stronger expressive powers and abilities to extract richer representations of graph-structured data. For empirical validation, we test different instances of the proposed architectures on multiple node classification tasks. The results show that even the simplest instance of the architectures achieves state-of-the-art performance, and the complex ones achieve surprisingly higher performance. In subsection 2.2, we propose to study an over-smoothing problem and give some ideas.

### A Stronger Multi-scale Deep GNN with Truncated Krylov Architecture

Suppose we deepen GCN in the same way as [32, 24], we have

\[Y=\text{softmax}(\hat{A}_{\text{sym}}\text{ ReLU}(\cdots\hat{A}_{\text{sym}}\text{ ReLU}(\hat{A}_{\text{sym}}\text{ ReLU}(\hat{A}_{\text{sym}}XW_{0})\text{ }W_{1})\text{ }W_{2}\cdots)\text{ }W_{n})\equiv\text{softmax}(Y^{ \prime}) \tag{3}\]

For this architecture, without considering the ReLU activation function, [32] shows that \(Y^{\prime}\) will converge to a space spanned by the eigenvectors of \(\hat{A}_{\text{sym}}\) with eigenvalue 1. Taking activation function into consideration, our analyses on (3) can be summarized in the following theorems (see proof in the appendix of [44]).

**Theorem 1**.: Suppose that \(\mathcal{G}\) has \(k\) connected components. Let \(X\in\mathbb{R}^{N\times F}\) be any feature matrix and let \(W_{j}\) be any non-negative parameter matrix with \(\|W_{j}\|_{2}\leq 1\) for \(j=0,1,\ldots\). If \(\mathcal{G}\) has no bipartite components, then in (3), as \(n\rightarrow\infty\), \(\text{rank}(Y^{\prime})\leq k\).
**Theorem 2**.: Suppose the \(n\)-dimensional \(\mathbf{x}\) and \(\mathbf{y}\) are independently sampled from a continuous distribution and the activation function \(\mathrm{Tanh}(z)=\frac{e^{z}-e^{-z}}{e^{z}+e^{-z}}\) is applied to \([\mathbf{x},\mathbf{y}]\) pointwisely, then

\[\mathbbm{P}(\mathrm{rank}\left(\mathrm{Tanh}([\mathbf{x},\mathbf{y}])\right)=\mathrm{ rank}([\mathbf{x},\mathbf{y}]))=1\]

Theorem 1 shows that, even considering ReLU, if we simply deepen GCN as (3), the extracted features will degrade under certain conditions, _i.e._\(Y^{\prime}\) only contains the stationary information of the graph structure and loses all the local information in node for being smoothed. In addition, from the proof we see that the pointwise ReLU transformation is a conspirator. Theorem 2 tells us that Tanh is better at keeping linear independence among column features. We design a numerical experiment on synthetic data to test, under a 100-layer GCN architecture, how activation functions affect the rank of the output in each hidden layer during the feedforward process. As Figure 1(a) shows, the rank of hidden features decreases rapidly with ReLU, while having little fluctuation under Tanh, and even the identity function performs better than ReLU. So we propose to replace ReLU by Tanh.

Besides activation function, to find a way to deepen GCN, we first show that any graph convolution with well-defined analytic spectral filter defined on \(\hat{A}_{\mathrm{sym}}\in\mathbbm{R}^{N\times N}\) can be written as a product of a block Krylov matrix with a learnable parameter matrix in a specific form. Based on this, we propose snowball network and truncated Krylov network.

Figure 1: Changes in the number of independent features with the increment of network depth

Figure 2: Snowball and Truncated Krylov Architectures
We take \(\mathsf{S}=\mathbb{R}^{F\times F}\). Given a set of block vectors \(\{X_{k}\}_{k=1}^{m}\subset\mathbb{R}^{N\times F}\), the \(\mathsf{S}\)-span of \(\{X_{k}\}_{k=1}^{m}\) is defined as \(\text{span}^{\mathsf{S}}\{X_{1},\ldots,X_{m}\}:=\{\sum\limits_{k=1}^{m}X_{k}C_{ k}:C_{k}\in\mathsf{S}\}\). Then, the order-\(m\) block Krylov subspace with respect to the matrix \(A\in\mathbb{R}^{N\times N}\), the block vector \(B\in\mathbb{R}^{N\times F}\) and the vector space \(\mathsf{S}\), and its corresponding block Krylov matrix are respectively defined as

\[\mathcal{K}_{m}^{\mathsf{S}}(A,B)\equiv\text{span}^{\mathsf{S}}\{B,AB,\ldots, A^{m-1}B\},\;\;K_{m}(A,B)\equiv[B,AB,\ldots,A^{m-1}B]\in\mathbb{R}^{N\times mF}.\]

It is shown in [11, 15] that there exists a smallest \(m\) such that for any \(k\geq m\), \(A^{k}B\in\mathcal{K}_{m}^{\mathsf{S}}(A,B)\), where \(m\) depends on \(A\) and \(B\).

Let \(\rho(\hat{A}_{\text{sym}})\) denote the spectrum radius of \(\hat{A}_{\text{sym}}\) and suppose \(\rho(\hat{A}_{\text{sym}})<R\) where \(R\) is the radius of convergence for a real analytic scalar function \(g\). Based on the above definitions and conclusions, the graph convolution can be written as

\[g(\hat{A}_{\text{sym}})X=\sum\limits_{n=0}^{\infty}\frac{g^{(n)}(0)}{n!}\hat{A }_{\text{sym}}^{n}X\equiv\left[X,\hat{A}_{\text{sym}}X,\ldots,\hat{A}_{\text{ sym}}^{m-1}X\right]\left[(\Gamma_{0}{}^{\mathsf{S}})^{T},(\Gamma_{1}{}^{\mathsf{S}})^{T}, \cdots,(\Gamma_{m-1}^{\mathsf{S}})^{T}\right]^{T}\equiv K_{m}(\hat{A}_{\text{ sym}},X)\Gamma^{\mathsf{S}} \tag{4}\]

where \(\Gamma_{i}^{\mathsf{S}}\in\mathbb{R}^{F\times F}\) for \(i=1,\ldots,m-1\) are parameter matrix blocks and \(\Gamma^{\mathsf{S}}\in\mathbb{R}^{mF\times F}\). Then, a graph convolutional layer can generally be written as

\[g(\hat{A}_{\text{sym}})XW^{\prime}=K_{m}(\hat{A}_{\text{sym}},X)\Gamma^{ \mathsf{S}}W^{\prime}=K_{m}(\hat{A}_{\text{sym}},X)W^{\mathsf{S}} \tag{5}\]

where \(W^{\prime}\in\mathbb{R}^{F\times O}\) is a parameter matrix, and \(W^{\mathsf{S}}\equiv\Gamma^{\mathsf{S}}W^{\prime}\in\mathbb{R}^{mF\times O}\). The essential number of learnable parameters is \(mF\times O\).

The block Krylov form provides an insight about why an architecture that concatenates multi-scale features in each layer will boost the expressive power of GCN. Based on this idea, we propose the snowball and truncated Block Krylov architectures [44] shown in Figure 2, where we stack multi-scale information in each layer. From the performance comparison on semi-supervised node classification tasks with different label percentage in table 1, we can see that the proposed models consistently perform better than the state-of-the-art models, especially when there are less labeled nodes. See detailed experimental results in [44].

### Future Works on Over-smoothing

Weight Initialization for GNNsEven without aggregation in each hidden layer, an NN with deep architecture still suffers from vanishing activation variances and back-propagated gradients variance problem [13], which make the training of deep NN hard. In last decade, designing new parameter initialization methods is proved to be effective [13, 18] to address the variance reduction problem during feedforward and backpropagation process. This motivates us to investigate the variance propagation in GNNs and analyze if the current weight initialization methods are suitable for GNNs or not. To this end, we can show that the vanishing variance caused by aggregation operation in GNNs is more serious than NN. Designing a new parameter initialization scheme for GNNs is potentially a feasible way to address this problem and empirically achieves promising performance [45]. we will propose a new method in this subsection.

The current initialization scheme of GNNs still follows the Xavier initialization [13], _i.e._\(W_{i}\sim U\left[-\frac{\sqrt{6}}{\sqrt{n_{j}+n_{j+1}}},\frac{\sqrt{6}}{\sqrt{n_{j} +n_{j+1}}}\right]\), or He (or Kaiming) initialization [18], _i.e._\(W_{i}\sim N\left(0,\sqrt{2/n_{i}}\right)\)
which is designed for traditional multilayer perceptron (MLP), where \(W_{i}\) is the parameter matrix of layer \(i\) and \(n_{i}\) is the number of hidden units of layer \(i\). These two initialization methods are derived by studying the variance propagation between layers during feedforward and backpropagation process. These two processes are different in GNNs by an extra multiplication of aggregation operator \(\hat{A}\). To analyze the variance propagation, we use deep GCN as an example, use \(\hat{A}=\hat{A}_{\text{rw}}\) and decompose it as follows,

\[\begin{split}& Y_{0}=X,\ H_{1}=\hat{A}_{\text{rw}}XW_{0},\ Y_{1}=f(H_{1}),\ H_{l+1}=\hat{A}_{\text{rw}}Y_{l}W_{l},\ Y_{l+1}=f(H_{l+1}),\ l=1,\ldots,n\\ & Y=\text{softmax}(\hat{A}_{\text{rw}}Y_{n}W_{n})\equiv\text{ softmax}(H_{n+1}),\ \mathcal{L}=-\text{trace}(Z^{T}\text{log}Y)\end{split} \tag{6}\]

where \(H_{l},Y_{l}\in\mathbb{R}^{N\times F_{l}}\), \(W_{l}\in\mathbb{R}^{F_{l}\times F_{l+1}}\); \(Z\in\mathbb{R}^{N\times C}\) is the ground truth matrix with one-hot label vector. Then the gradient propagates in the following way,

\[\frac{\partial\mathcal{L}}{\partial H_{l}}=\frac{\partial\mathcal{L}}{ \partial Y_{l}}\odot f^{\prime}(H_{l}),\ \frac{\partial\mathcal{L}}{\partial W_{l-1}}=Y_{l-1}^{T}\hat{A}_{\text{rw}} \frac{\partial\mathcal{L}}{\partial H_{l}},\ \frac{\partial\mathcal{L}}{ \partial Y_{l-1}}=\hat{A}_{\text{rw}}\frac{\partial\mathcal{L}}{\partial H_{l }}W_{l-1}^{T} \tag{7}\]

Variance Analysis: Forward ViewConsider element \(i,j\) in matrix \(H_{l+1}\) during the feed-forward process in (6),

\[(H_{l+1})_{ij}=(\hat{A}_{\text{rw}})_{i,z}Y_{l}(W_{l})_{:,j}=\sum_{t=1}^{F_{l }}\sum_{k=1}^{N}(\hat{A}_{\text{rw}})_{ik}\left(Y_{l}\right)_{kt}\left(W_{l} \right)_{t,j},\ Y_{l+1}=f(H_{l+1}),\ l=1,\ldots,n \tag{8}\]

Suppose we have linear activation function such as that proposed in [60]; each element in \(W_{l}\) is _i.i.d._ initialized with \(E\left((W_{l})_{ij}\right)=0\); \(E\left((Y_{l})_{kt}\right)=0\) and all elements in \(Y_{l}\) are independent 2. Then, \(\text{Var}\left((Y_{l+1})_{ij}\right)=\text{Var}\left((Y_{l+1})\right)\) can be written as

\begin{table}
\begin{tabular}{c|c c c c|c c c|c c c c} \multirow{2}{*}{Algorithms} & \multicolumn{4}{c|}{Cora} & \multicolumn{4}{c}{CiteSeer} & \multicolumn{4}{c}{PubMedMed} \\  & 0.5\% & 1\% & 2\% & 3\% & 4\% & 5\% & 0.5\% & 1\% & 2\% & 3\% & 4\% & 5\% & 0.03\% & 0.05\% & 0.1\% & 0.3\% \\ \hline LP & 56.4 & 62.3 & 65.4 & 67.5 & 69.0 & 70.2 & 34.8 & 40.2 & 45.6 & 45.3 & 46.4 & 47.3 & 61.4 & 66.4 & 65.4 & 66.8 \\ Cheby & 38.0 & 52.0 & 62.4 & 70.8 & 74.1 & 77.6 & 31.7 & 42.8 & 59.9 & 66.2 & 68.3 & 69.3 & **40.4** & **73.5** & 51.2 & 72.8 \\ Co-training & 56.6 & 66.4 & 73.5 & 75.9 & 78.9 & 80.8 & 47.5 & 57.6 & 62.1 & 62.5 & 64.5 & 65.5 & 62.2 & 68.3 & 72.7 & 78.2 \\ Self-training & 53.7 & 66.1 & 73.8 & 77.2 & 79.4 & 80.0 & 43.3 & 58.1 & 68.2 & 69.8 & 70.4 & 71.0 & 51.9 & 87.6 & 66.8 & 77.0 \\ Union & 58.5 & 69.9 & 75.9 & 78.5 & 80.4 & 81.7 & 46.3 & 59.1 & 66.7 & 67.6 & 67.6 & 68.2 & 58.4 & 64.0 & 70.7 & 79.2 \\ Intersection & 49.7 & 65.0 & 72.9 & 77.1 & 79.4 & 80.2 & 42.9 & 59.1 & 68.6 & 70.1 & 70.8 & 71.2 & 52.0 & 59.3 & 69.7 & 77.6 \\ MultiStage & 61.1 & 63.7 & 74.4 & 76.1 & 77.2 & 53.0 & 78.3 & 68.3 & 68.0 & 69.0 & 57.4 & 64.3 & 70.2 \\ M3S & 61.5 & 67.2 & 76.5 & 77.8 & 78.0 & 56.1 & 62.1 & 66.4 & 70.3 & 70.5 & 59.2 & 64.4 & 70.6 \\ GCN & 42.6 & 56.9 & 67.8 & 49.7 & 77.6 & 79.3 & 38.4 & 46.5 & 62.6 & 66.9 & 68.7 & 69.6 & 46.4 & 49.7 & 56.3 & 76.6 \\ GCN-SVAT & 43.6 & 53.9 & 71.4 & 75.6 & 78.3 & 78.5 & 47.0 & 52.4 & 65.8 & 68.6 & 69.5 & 70.7 & 52.1 & 56.9 & 63.5 & 77.2 \\ GCN-DVAT & 49 & 61.8 & 71.9 & 75.9 & 78.4 & 78.6 & 51.5 & 58.5 & 67.4 & 69.2 & 70.8 & 71.3 & 53.3 & 58.6 & 66.3 & 77.3 \\ \hline \multicolumn{2}{c}{_linear Snowball_} & _67.6_ & _74.6_ & _78.9_ & _80.9_ & _82.3_ & _82.9_ & _56.0_ & _63.4_ & _69.3_ & _70.6_ & **72.5** & **72.6** & **65.5** & _68.5_ & _73.6_ & _79.7_ \\ \multicolumn{2}{c}{_Snowball_} & _68.4_ & _73.2_ & _78.4_ & _80.8_ & _82.3_ & _83.0_ & _56.4_ & _63.9_ & _68.7_ & _70.5_ & _71.8_ & _72.8_ & _66.5_ & _68.6_ & _73.2_ & _80.1_ \\ \multicolumn{2}{c}{_truncated Krylov_} & _71.8_ & _76.5_ & _80.0_ & _82.0_ & _83.0_ & _94.1_ & _59.9_ & _66.1_ & _69.8_ & _71.3_ & _72.3_ & _73.7_ & **68.7** & _71.4_ & _75.5_ & _80.4_ \\ \hline \hline \end{tabular}

* For each (column), the greener the cell, the better the performance. The redder, the worse. If our methods achieve better performance than all others, the corresponding cell will be in bold.

\end{table}
Table 1: Accuracy without Validation
\[\text{Var}\left(\sum_{t=1}^{F_{l}}\sum_{k=1}^{N}(\hat{A}_{\text{rw}})_{jk}\left(Y_ {l}\right)_{kt}\left(W_{l}\right)_{t,j}\right)=\sum_{t=1}^{F_{l}}\sum_{k=1}^{N} \text{Var}\left((\hat{A}_{\text{rw}})_{ik}\left(Y_{l}\right)_{kt}\left(W_{l} \right)_{t,j}\right)=\frac{F_{l}}{d_{i}+1}\text{Var}\left(Y_{l}\right)\text{ Var}(W_{l}) \tag{9}\]

Suppose each element in \(Y_{l}\) shares the same variance denoted as \(\text{Var}\left(Y_{l}\right)\). To prevent variance vanishing between layers, _i.e._\(\text{Var}\left(Y_{l+1}\right)=\text{Var}\left(Y_{l}\right)\), from (8) we can approximately have (see computation in Appendix A.2.1)

\[\text{Var}(W_{l})=\frac{d_{i}+1}{F_{l}} \tag{10}\]

This tells us that the variance of \(W_{l}\) depends on the degree of a node, but since the parameter matrix is shared by all nodes, we cannot design a node specified initialization scheme. Thus, we make a compromise between nodes as follows

\[\text{Var}(W_{l})\approx\frac{\sum\limits_{i=1}^{N}(d_{i}+1)}{NF_{l}}=\frac{1 +\text{average node degree}}{F_{l}} \tag{11}\]

Another way is to use weighted average by considering the node degree as the weight of each node. Through this way, we have

\[\text{Var}(W)=\sum\limits_{i=1}^{N}\frac{d_{i}+1}{\sum\limits_{j=1}^{N}d_{j}+ 1}\frac{d_{i}+1}{F_{l}}=\frac{\sum\limits_{i=1}^{N}(d_{i}+1)^{2}}{(\sum\limits _{i=1}^{N}d_{i}+1)F_{l}} \tag{12}\]

Variance Analysis: Backward ViewUnder the same assumption as feedforward view and suppose each element in \(\frac{\partial\mathcal{L}}{\partial H_{l}}\) and \(\frac{\partial\mathcal{L}}{\partial W_{l-1}}\) are independent to each other and has zero mean, from (7) we can approximately (see computation in Appendix A.2.2)

\[\frac{\partial\mathcal{L}}{\partial H_{l}}=\frac{\partial\mathcal{L}}{\partial Y _{l}}=\hat{A}_{\text{rw}}\frac{\partial\mathcal{L}}{\partial H_{l+1}}W_{l}^{ T},\ \frac{\partial\mathcal{L}}{\partial W_{l-1}}=Y_{l-1}^{T}\hat{A}_{\text{rw}}\frac{ \partial\mathcal{L}}{\partial H_{l}} \tag{13}\]

Then,

\[\left(\frac{\partial\mathcal{L}}{\partial H_{l}}\right)_{ij}=\sum \limits_{t=1}^{F_{l+1}}\sum\limits_{k=1}^{N}(\hat{A}_{\text{rw}})_{ik}\left( \frac{\partial\mathcal{L}}{\partial H_{l+1}}\right)_{kt}(W_{l}^{T})_{t,j} \tag{14}\] \[\left(\frac{\partial\mathcal{L}}{\partial W_{l-1}}\right)_{ij}=( \hat{A}_{\text{rw}}Y_{l-1})_{i}^{T}\left(\frac{\partial\mathcal{L}}{\partial H _{l}}\right)_{.j}=\sum\limits_{k=1}^{N}(\sum\limits_{t=1}^{N}(\hat{A}_{\text{rw }})_{kt}(Y_{l-1})_{ti})\left(\frac{\partial\mathcal{L}}{\partial H_{l}}\right)_ {kj},\]
Thus,

\[\text{Var}\left(\left(\frac{\partial\mathcal{L}}{\partial H_{l}} \right)_{ij}\right) =\text{Var}\left(\sum\limits_{t=1}^{F_{l+1}}\sum\limits_{k=1}^{N} \hat{A}_{\text{rw}_{k}}\left(\frac{\partial\mathcal{L}}{\partial H_{l+1}} \right)_{kt}(W_{l}^{T})_{t,j}\right)=\frac{F_{l+1}}{d_{i}+1}\text{Var}\left( \frac{\partial\mathcal{L}}{\partial H_{l+1}}\right)\text{Var}\left(W_{l}\right) \tag{15}\] \[\text{Var}\left(\left(\frac{\partial\mathcal{L}}{\partial W_{l-1}} \right)_{ij}\right) =\text{Var}\left(\sum\limits_{k=1}^{N}(\sum\limits_{t=1}^{N}( \hat{A}_{\text{rw}})_{kt}(Y_{l-1})_{ti})\left(\frac{\partial\mathcal{L}}{ \partial H_{l}}\right)_{kj}\right)=\left(\sum\limits_{k=1}^{N}\frac{1}{d_{k}+1 }\right)\text{Var}\left(Y_{l-1}\right)\text{Var}\left(\frac{\partial\mathcal{ L}}{\partial H_{l}}\right)\]

\[\text{Var}(W_{l})=\frac{\sum\limits_{i=1}^{N}(d_{i}+1)}{NF_{l+1}}\approx\frac {1+\text{average node degree}}{F_{l+1}} \tag{16}\]

From (9)(15), \(\text{Var}\left(Y_{l-1}\right)\) can be approximately written as

\[\text{Var}\left(\frac{\partial\mathcal{L}}{\partial H_{l}} \right)\approx\frac{NF_{l+1}}{\sum\limits_{i=1}^{N}d_{i}+1}\text{Var}\left( \frac{\partial\mathcal{L}}{\partial H_{l+1}}\right)\text{Var}\left(W_{l} \right)\approx\text{Var}\left(\frac{\partial\mathcal{L}}{\partial H_{l+1}} \right)\prod\limits_{l=1}^{n+1}\frac{NF_{l^{\prime}}}{\sum\limits_{i=1}^{N}d _{i}+1}\text{Var}\left(W_{l^{\prime}-1}\right) \tag{17}\] \[\text{Var}\left(Y_{l-1}\right)\approx\frac{NF_{l-2}}{\sum_{k}d_ {k}+1}\text{Var}\left(\mathbf{Y_{l-2}}\right)\text{Var}(W_{l-2})\approx\text{ Var}(\mathbf{Y_{0}})\prod\limits_{l^{\prime}=0}^{l-2}\frac{NF_{l^{\prime}}}{\sum_{k}d_{k}+1} \text{Var}(W_{l^{\prime}})\]

From (15), if each \(\text{Var}(W_{l^{\prime}})\) equals to \(\text{Var}(W)\) and each \(F_{l}\) equals to \(F\), then

\[\text{Var}\left(\left(\frac{\partial\mathcal{L}}{\partial W_{l-1}} \right)_{ij}\right)\approx\left(\sum\limits_{k=1}^{N}\frac{1}{d_{k}+1}\right) \text{Var}(\mathbf{Y_{0}})\text{Var}\left(\frac{\partial\mathcal{L}}{\partial H_{ n+1}}\right)\left(\frac{NF}{\sum_{k}d_{k}+1}\text{Var}(W)\right)^{n} \tag{18}\]

Combined with (11), we can set the variance of the parameter matrix as

\[\text{Var}(W_{l})\approx\frac{2\sum\limits_{i=1}^{N}(d_{i}+1)}{N(F_{l}+F_{l+1} )}=\frac{2(1+\text{average node degree})}{(F_{l}+F_{l+1})} \tag{19}\]

Thus, each element in \(W_{i}\) can be drawn from \(N(0,\sqrt{\frac{2(1+\text{average node degree})}{(F_{l}+F_{l+1})}})\).

Adaptive ReLU (AdaReLU) Activation FunctionTo satisfy the assumption that the activation function is linear at the beginning of training process and to still learn a nonlinear function during training, we design the following adaptive ReLU (AdaReLU) activation function

\[f(\mathbf{x}_{i})=\left\{\begin{array}{ll}\beta_{i}\mathbf{x}_{i},&\text{if }\mathbf{x}_{i}>0\\ \alpha_{i}\mathbf{x}_{i},&\text{if }\mathbf{x}_{i}\leq 0\end{array}\right.\]

where \(\alpha_{i}\) and \(\beta_{i}\) are learnable parameters and are initialized to be 1. If \(\alpha_{i}=\alpha\) and \(\beta_{i}=\beta\) for all \(i\), we have channel-shared AdaReLU, otherwise we have channel-wise AdaReLU 3. From the preliminary experimental results, channel-wise works better than channel-shared AdaReLU.
There exist some experimental evidence [45] that controlling variance flow by initialization like (19) can relieve the performance decrease of deep GCN. But more tests and hyperparameter tunning still needs to be done. More theoretical analysis on variance propagation needs to be done.

## 3 GNNs on Heterophily Graphs

GNNs are considered as an extension of basic Neural Networks (NNs) by additionally making use of graph structure based on the relational inductive bias (homophily assumption), rather than treating the nodes as collections of independent and identically distributed (_i.i.d._) samples. Though GNNs are believed to outperform basic NNs in real-world tasks, it is found that in some cases, the graph-aware models have little performance gain or even underperform graph-agnostic models [6, 46, 51, 69, 71]. One of the main reasons for the performance degradation is believed to be heterophily, _i.e._ when the connected nodes tend to have different labels [69, 71]. Heterophily challenge has received attention recently and there are increasing number of models being put forward to analyze [39, 43] and address this problem [6, 36, 41, 69, 70, 42, 64].

In this section, we first introduce the most commonly used homophily metrics in subsection 3.1. Then, we show that not all cases of heterophily are harmful for GNNs and propose new metrics based on a similarity matrix which considers the influence of both graph structure and input features on GNNs in subsection 3.2. The metrics demonstrate advantages over the commonly used homophily metrics by tests on synthetic graphs. From the metrics and the observations, we find some cases of harmful heterophily can be addressed by diversification operation and its effectiveness can be proved in subsection 3.3. With this fact and knowledge of filterbanks, we propose the Adaptive Channel Mixing (ACM) framework in subsection 3.4 to adaptively exploit aggregation, diversification and identity channels in each GNN layer to address harmful heterophily. We validate the ACM-augmented baselines with real-world node classification tasks. They consistently achieve significant performance gain and exceed the state-of-the-art GNNs on most of the tasks without incurring significant computational burden. In subsection 3.5, we introduce some prior work on addressing heterophily problems and explain their differences with ACM framework. The limitation of diversification operation and remaining challenges of heterophily problems are discussed in subsection 3.6

### Metrics of Homophily

The metrics of homophily are defined by considering different relations between node labels and graph structures defined by adjacency matrix. There are three commonly used homophily metrics: edge homophily [70, 1], node homophily [51], and class homophily [35]4 defined as follows:

Footnote 4: The authors in [35] did not name this homophily metric. We name it class homophily based on its definition.

\[H_{\text{edge}}(\mathcal{G}) =\frac{\left|\left[e_{uv}\mid e_{uv}\in\mathcal{E},Z_{u,:}=Z_{v,: }\right]\right|}{\left|\mathcal{E}\right|},\;\;H_{\text{node}}(\mathcal{G})= \frac{1}{\left|\mathcal{V}\right|}\sum_{v\in\mathcal{V}}\frac{\left|\left[u \mid u\in\mathcal{N}_{v},Z_{u,:}=Z_{v,:}\right]\right|}{d_{v}}, \tag{20}\] \[H_{\text{class}}(\mathcal{G}) =\frac{1}{C-1}\sum_{k=1}^{C}\left[h_{k}-\frac{\left|\left\{v \mid Z_{v,k}=1\right\}\right|}{N}\right]_{+},\;\;h_{k}=\frac{\sum_{v\in \mathcal{V}}\left|\left\{u\mid Z_{v,k}=1,u\in\mathcal{N}_{v},Z_{u,:}=Z_{v,:} \right\}\right|}{\sum_{v\in\left[v\mid Z_{v,k}=1\right]d_{v}}}\]
where \([a]_{+}=\max(a,0)\); \(h_{k}\) is the class-wise homophily metric [35]. They are all in the range of \([0,1]\) and a value close to \(1\) corresponds to strong homophily while a value close to \(0\) indicates strong heterophily. \(H_{\text{edge}}(\mathcal{G})\) measures the proportion of edges that connect two nodes in the same class; \(H_{\text{node}}(\mathcal{G})\) evaluates the average proportion of edge-label consistency of all nodes; \(H_{\text{class}}(\mathcal{G})\) tries to avoid the sensitivity to imbalanced class, which can cause \(H_{\text{edge}}\) misleadingly large. The above definitions are all based on the graph-label consistency and imply that the inconsistency will cause harmful effect to the performance of GNNs. With this in mind, we will show a counter example to illustrate the insufficiency of the above metrics and propose new metrics in the following subsection.

### Analysis of Heterophily and Aggregation Homophily Metric

Heterophily is believed to be harmful for message-passing based GNNs [6, 51, 70] because intuitively features of nodes in different classes will be falsely mixed and this will lead nodes indistinguishable [70]. Nevertheless, it is not always the case, _e.g._ the bipartite graph shown in Figure 3 is highly heterophilous according to the homophily metrics in (20), but after mean aggregation, the nodes in classes 1 and 2 only exchange colors and are still distinguishable. Authors in [6] also point out the insufficiency of \(H_{\text{node}}\) by examples to show that different graph typologies with the same \(H_{\text{node}}\) can carry different label information.

To analyze to what extent the graph structure can affect the output of a GNN, we first simplify the GCN by removing its nonlinearity as [60]. Let \(\hat{A}\in\mathbb{R}^{N\times N}\) denote a general aggregation operator. Then, equation (1) can be simplified as,

\[Y=\text{softmax}(\hat{A}XW)=\text{softmax}(Y^{\prime}) \tag{21}\]

After each gradient decent step \(\Delta W=\gamma\frac{d\mathcal{L}}{dW}\), where \(\gamma\) is the learning rate, the update of \(Y^{\prime}\) will be,

\[\Delta Y^{\prime}=\hat{A}X\Delta W=\gamma\hat{A}X\frac{d\mathcal{L}}{dW} \propto\hat{A}X\frac{d\mathcal{L}}{dW}=\hat{A}XX^{T}\hat{A}^{T}(Z-Y)=S(\hat{A},X)(Z-Y) \tag{22}\]

where \(S(\hat{A},X)\equiv\hat{A}X(\hat{A}X)^{T}\) is a post-aggregation node similarity matrix, \(Z-Y\) is the prediction error matrix. The update direction of node \(i\) is essentially a weighted sum of the prediction error, _i.e._\(\Delta(Y^{\prime})_{i,:}=\sum_{j\in\mathcal{V}}\left[S(\hat{A},X)\right]_{i,j }(Z-Y)_{j,:}\).

To study the effect of heterophily, we first define the _aggregation similarity score_ as follows.

**Definition 1**.: _Aggregation similarity score_

\[S_{agg}\left(S(\hat{A},X)\right)=\frac{\left|\left[v\left|\operatorname{Mean }_{u}\left(\{S(\hat{A},X)_{v,u}|Z_{u,:}=Z_{v,:}\}\right)\right.\right.\right) \geq\operatorname{Mean}_{u}\left(\{S(\hat{A},X)_{v,u}|Z_{u,:}\neq Z_{v,:}\} \right)\right|}{\left|\mathcal{V}\right|} \tag{23}\]

_where \(\operatorname{Mean}_{u}\left(\{\cdot\}\right)\) takes the average over \(u\) of a given multiset of values or variables._

Figure 3: Example of harmless heterophily
\(S_{\text{agg}}(S(\hat{A},X))\) measures the proportion of nodes \(v\in\mathcal{V}\) that will put relatively larger similarity weights on nodes in the same class than in other classes after aggregation. It is easy to see that \(S_{\text{agg}}(S(\hat{A},X))\in[0,1]\). But in practice, we observe that in most datasets, we will have \(S_{\text{agg}}(S(\hat{A},X))\geq 0.5\). Based on this observation, we rescale (23) to the following modified aggregation similarity for practical usage,

\[S_{\text{agg}}^{M}\left(S(\hat{A},X)\right)=\left[2S_{\text{agg}}\left(S(\hat {A},X)\right)-1\right]_{+} \tag{24}\]

In order to measure the consistency between labels and graph structures without considering node features and make a fair comparison with the existing homophily metrics in (20), we define the graph (\(\mathcal{G}\)) aggregation (\(\hat{A}\)) homophily and its modified version as

\[H_{\text{agg}}(\mathcal{G})=S_{\text{agg}}\left(S(\hat{A},Z)\right),\ H_{ \text{agg}}^{M}(\mathcal{G})=S_{\text{agg}}^{M}\left(S(\hat{A},Z)\right) \tag{25}\]

In practice, we will only check \(H_{\text{agg}}(\mathcal{G})\) when \(H_{\text{agg}}^{M}(\mathcal{G})=0\). As Figure 3 shows, when \(\hat{A}=\hat{A}_{\text{rw}}\), \(H_{\text{agg}}(\mathcal{G})=H_{\text{agg}}^{M}(\mathcal{G})=1\). Thus, this new metric reflects the fact that nodes in classes 1 and 2 are still highly distinguishable after aggregation, while other metrics mentioned before fail to capture the information and misleadingly give value 0. This shows the advantage of \(H_{\text{agg}}(\mathcal{G})\) and \(H_{\text{agg}}^{M}(\mathcal{G})\) by additionally considering information from aggregation operator \(\hat{A}\) and the similarity matrix.

Comparison of Homophily Metrics on Synthetic GraphsTo comprehensively compare \(H_{\text{agg}}^{M}(\mathcal{G})\) with the metrics in (20) in terms of how they reveal the influence of graph structure on the GNN performance, we generate synthetic graphs (\(d\)-regular graphs with edge homophily varied from 0.005 to 0.95) and evaluate SGC with 1-hop aggregation (SGC-1) [60] and GCN [24] on them.

The performance of SGC-1 and GCN are expected to be monotonically increasing with a proper and informative homophily metric. However, Figure 4(a)(b)(c) show that the performance curves under \(H_{\text{edge}}(\mathcal{G}),H_{\text{node}}(\mathcal{G})\) and \(H_{\text{class}}(\mathcal{G})\) are \(U\)-shaped 5, while Figure 4(d) reveals a nearly monotonic curve only with a little numerical perturbation around 1. This indicates that \(H_{\text{agg}}^{M}(\mathcal{G})\) can describe how the graph structure affects the performance of SGC-1 and GCN more appropriately and adequately than the existing metrics.

Figure 4: Comparison of baseline performance under different homophily metrics.
### How Diversification Operation Helps with Harmful Heterophily

We first consider the example shown in Figure 5. From \(S(\hat{A},X)\), nodes 1,3 assign relatively large positive weights to nodes in class 2 after aggregation, which will make node 1,3 hard to be distinguished from nodes in class 2. Despite the fact, we can still distinguish between nodes 1,3 and 4,5,6,7 by considering their neighborhood difference: nodes 1,3 are different from most of their neighbors while nodes 4,5,6,7 are similar to most of their neighbors. This indicates, in some cases, although some nodes become similar after aggregation, they are still distinguishable via their surrounding dissimilarities. This leads us to use _diversification operation_, _i.e._ high-pass (HP) filter \(I-\hat{A}\)[10] (will be introduced in the next subsection) to extract the information of neighborhood differences and address harmful heterophily. As \(S(I-\hat{A},X)\) in Figure 5 shows, nodes 1,3 will assign negative weights to nodes 4,5,6,7 after diversification operation, _i.e._ nodes 1,3 treat nodes 4,5,6,7 as negative samples and will move away from them during backpropagation. Base on this example, we first propose diversification distinguishability as follows to measures the proportion of nodes that diversification operation is potentially helpful for,

**Definition 2**.: _Diversification Distinguishability (DD) based on \(S(I-\hat{A},X)\)._

_Given \(S(I-\hat{A},X)\), a node \(v\) is diversification distinguishable if the following two conditions are satisfied at the same time,_

\[\begin{split}\textbf{1.}&\ \mathrm{Mean}_{u}\left(\{S(I-\hat{A},X)_{v,u}|u \in\mathcal{V}\wedge Z_{u,z}=Z_{v,z}\}\right)>0;\\ \textbf{2.}&\ \mathrm{Mean}_{u}\left(\{S(I-\hat{A},X)_{v,u }|u\in\mathcal{V}\wedge Z_{u,z}\neq Z_{v,z}\}\right)\leq 0\end{split} \tag{26}\]

_Then, graph diversification distinguishability value is defined as_

\[\mathrm{DD}_{\hat{A},X}(\mathcal{G})=\frac{1}{|\mathcal{V}|}\left|\{v|v\ \text{is diversification distinguishable}\}\right| \tag{27}\]

We can see that \(\mathrm{DD}_{\hat{A},X}(\mathcal{G})\in[0,1]\). The effectiveness of diversification operation can be proved for binary classification problems under certain conditions based on definition 2, leading us to:

**Theorem 3**.: Suppose \(X=Z,\hat{A}=\hat{A}_{\mathrm{rw}}\). Then, for a binary classification problem, _i.e._\(C=2\), all nodes are diversification distinguishable, _i.e._\(\mathrm{DD}_{\hat{A},Z}(\mathcal{G})=1\).

Figure 5: Example of how HP filter addresses harmful heterophily
Theorem 3 theoretically demonstrates the importance of diversification operation to extract high-frequency information of graph signal [10]. Combined with aggregation operation, which is a low-pass filter [10, 48], we can get a filterbank which uses both aggregation and diversification operations to distinctively extract the low- and high-frequency information from graph signals. We will introduce filterbank in the next subsection.

### Filterbank and Adaptive Channel Mixing(ACM) GNN Framework

FilterbankFor the graph signal \(\mathbf{x}\) defined on \(\mathbf{\mathcal{G}}\), a 2-channel linear (analysis) filterbank [10]6 includes a pair of low-pass(LP) and high-pass(HP) filters \(H_{\text{LP}},H_{\text{HP}}\), where \(H_{\text{LP}}\) and \(H_{\text{HP}}\) retain the low-frequency and high-frequency content of \(\mathbf{x}\), respectively. Filterbanks with \(H_{\text{LP}}+H_{\text{HP}}=I\) will not lose any information of the input signal, _i.e._ perfect reconstruction property [10].

Footnote 6: In graph signal processing, an additional synthesis filter [10] is required to form the 2-channel filterbank. But synthesis filter is not needed in our framework, so we do not introduce it in our paper.

However, most existing GNNs are under uni-channel filtering architecture [17, 24, 58] with either \(H_{\text{LP}}\) or \(H_{\text{HP}}\) channel that only partially preserves the input information. Generally, the Laplacian matrices (\(L_{\text{sym}}\), \(L_{\text{rw}}\), \(\hat{L}_{\text{sym}}\), \(\hat{L}_{\text{rw}}\)) can be regarded as HP filters [10] and affinity matrices (\(A_{\text{sym}}\), \(A_{\text{rw}}\), \(\hat{A}_{\text{sym}}\), \(\hat{A}_{\text{rw}}\)) can be treated as LP filters [48, 16]. Moreover, we consider MLPs as owing a special identity filterbank with matrix \(I\) that satisfies \(H_{\text{LP}}\) + \(H_{\text{HP}}\) = \(I+0=I\).

Filterbank in Spatial FormFilterbank methods can also be extended to spatial GNNs. Formally, on the node level, left multiplying \(H_{\text{LP}}\) and \(H_{\text{HP}}\) on \(\mathbf{x}\) performs as aggregation and diversification operations, respectively. For example, suppose \(H_{\text{LP}}\) = \(\hat{A}\) and \(H_{\text{HP}}\) = \(I-\hat{A}\), then for node \(i\) we have

\[(H_{\text{LP}}\mathbf{x})_{i}=\sum_{j\in[\mathcal{N}_{i}\cup i]}\hat{A}_{ij}\mathbf{x _{j}},\ (H_{\text{HP}}\mathbf{x})_{i}=\mathbf{x_{i}}-\sum_{j\in[\mathcal{N}_{i}\cup i]}\hat{A} _{ij}\mathbf{x_{j}} \tag{28}\]

where \(\hat{A}_{i,j}\) is the connection weight between two nodes. To leverage HP and identity channels in GNNs, we propose the Adaptive Channel Mixing (ACM) framework which can be applied to lots of baseline GNN. We use GCN as an example and introduce ACM framework in matrix form. We use \(H_{\text{LP}}\) and \(H_{\text{HP}}\) to represent general LP and HP filters. The ACM framework includes 3 steps as
follows,

**Step 1. Feature Extraction for Each Channel:**

Option 1: \(H_{L}^{l}=\text{ReLU}\left(H_{\text{LP}}H^{l-1}W_{L}^{l-1}\right),\ H_{H}^{l}= \text{ReLU}\left(H_{\text{HP}}H^{l-1}W_{H}^{l-1}\right),\ H_{I}^{l}=\text{ ReLU}\left(H^{l-1}W_{I}^{l-1}\right);\)

Option 2: \(H_{L}^{l}=H_{\text{LP}}\text{ReLU}\left(H^{l-1}W_{L}^{l-1}\right),\ H_{H}^{l}=H_{ \text{HP}}\text{ReLU}\left(H^{l-1}W_{H}^{l-1}\right),\ H_{I}^{l}=I\text{ ReLU}\left(H^{l-1}W_{I}^{l-1}\right);\)

\(W_{L}^{l-1},\ W_{H}^{l-1},\ W_{I}^{l-1}\in\mathbb{R}^{F_{l-1}\times F_{l}}\);

**Step 2. Feature-based Weight Learning**

\[\bar{\alpha}_{L}^{l}=\sigma\left(H_{L}^{l}\bar{W}_{L}^{l}\right),\ \bar{\alpha}_{H}^{l}=\sigma\left(H_{H}^{l}\bar{W}_{I}^{l}\right),\ \bar{\alpha}_{I}^{l}=\sigma\left(H_{I}^{l}\bar{W}_{I}^{l}\right),\ \bar{W}_{L}^{l-1},\ \bar{W}_{H}^{l-1},\ \bar{W}_{I}^{l-1}\in \mathbb{R}^{F_{l}\times 1}\]

\[\left[\alpha_{L}^{l},\alpha_{H}^{l},\alpha_{I}^{l}\right]=\text{Softmax}\left( \left[\bar{\alpha}_{L}^{l},\bar{\alpha}_{H}^{l},\bar{\alpha}_{I}^{l}\right]W _{\text{Mix}}^{l}/T\right),\ W_{\text{Mix}}^{l}\in\mathbb{R}^{3\times 3},T\in \mathbb{R}\text{ is the temperature;}\]

**Step 3. Node-wise Channel Mixing:**

\[H^{l}=\left(\text{diag}(\alpha_{L}^{l})H_{L}^{l}+\text{diag}(\alpha_{H}^{l}) H_{H}^{l}+\text{diag}(\alpha_{I}^{l})H_{I}^{l}\right). \tag{29}\]

The framework with option 1 in step 1 is ACM framework and with option 2 is ACMII framework. ACM(II)-GCN first implement distinct feature extractions for 3 channels, respectively. After processed by a set of filterbanks, 3 filtered components \(H_{L}^{l},H_{H}^{l},H_{I}^{l}\) are obtained. Different nodes may have different needs for the information in the 3 channels, _e.g._ in Figure 5, nodes 1,3 demand high-frequency information while node 2 only needs low-frequency information. To adaptively exploit information from different channels, ACM(II)-GCN learns row-wise (node-wise) feature-conditioned weights to combine the 3 channels. ACM(II) can be easily plugged into spatial GNNs by replacing \(H_{\text{LP}}\) and \(H_{\text{HP}}\) by aggregation and diversification operations as shown in (28).

ComplexityNumber of learnable parameters in layer \(l\) of ACM(II)-GCN is \(3F_{l-1}(F_{l}+1)+9\), while it is \(F_{l-1}F_{l}\) in GCN. The computation of step 1-3 takes \(NF_{l}(8+6F_{l-1})+2F_{l}(\text{nnz}(H_{\text{LP}})+\text{nnz}(H_{\text{HP}}))+18N\) flops, while GCN layer takes \(2NF_{l-1}F_{l}+2F_{l}(\text{nnz}(H_{\text{LP}}))\) flops, where \(\text{nnz}(\cdot)\) is the number of non-zero elements.

Performance ComparisonWe implement SGC [60] with 1 hop and 2 hops (SGC-1, SGC-2), GCNII [5], GCNII*[5], GCN [24] and snowball networks with 2 and 3 layers (snowball-2, snowball-3) and apply them in ACM or ACMII framework: we use \(\hat{A}_{\text{rw}}\) as LP filter and the corresponding HP filter is \(I-\hat{A}_{\text{rw}}\). We compare them with several baseline and SOTA GNN models: MLP with 2 layers (MLP-2), GAT [58], APPNP [25], GPRGNN [6], H\({}_{2}\)GCN [70], MixHop [1], GCN+JK [24, 35, 63], GAT+JK [35, 58, 63], FAGCN [2] GraphSAGE [17] and Geom-GCN [51]. Besides the 9 benchmark datasets _Cornell_, _Wisconsin_, _Texas_, _Film_, _Chameleon_, _Squirrel_, _Cora_, _Citeseer_ and _Pubmed_ used in [51], we further test the above models on a new benchmark dataset, _Deezer-Europe_, that is proposed in [35]. On each dataset used in [51], we test the models 10 times following the same early stopping strategy, the same random data splitting method and Adam [23] optimizer as used in GPRGNN [6]. For _Deezer-Europe_, we test the above models 5 times with the same early stopping strategy, the same fixed splits and AdamW [37] used in [35].

To better visualize the performance boost and the comparison with SOTA models, in Figure 6, we plot the bar charts of the test accuracy of SOTA models, 3 selected baselines (GCN, snowball-2, snowball-3) and their ACM and ACMII augmented models on 6 most commonly used benchmark heterophily datasets (See [40] for the full results and comparison). We can see that after being
applied in ACM or ACMII framework, the performance of the 3 baseline models are significantly boosted on all tasks and can achieve SOTA performance. Especially on _Cornell, Texas, Film_ and _Squirrel_, the augmented models significantly outperform the current SOTA models. Overall, It suggests that ACM or ACMII framework can help GNNs to generalize better on node classification tasks on heterophilous graphs.

### Prior Work

We discuss relevant work of GNNs on addressing heterophily challenge in this part. Authors in [1] acknowledge the difficulty of learning on graphs with weak homophily and propose MixHop to extract features from multi-hop neighborhood to get more information. Geom-GCN [51] precomputes unsupervised node embeddings and uses graph structure defined by geometric relationships in the embedding space to define the bi-level aggregation process. Authors in [20] propose measurements based on feature smoothness and label smoothness that are potentially helpful to guide GNNs on dealing with heterophilous graphs. H\({}_{2}\)GCN [70] combines 3 key designs to address heterophily: (1) ego- and neighbor-embedding separation; (2) higher-order neighborhoods; (3) combination of intermediate representations. CPGNN [69] models label correlations by the compatibility matrix, which is beneficial for heterophily settings, and propagates a prior belief estimation into GNNs by the compatibility matrix. FBGNN [47] first proposes to use filterbank to address heterophily problem, but it does not fully explain the insights behind HP filters and does not contain identity channel and node-wise channel mixing mechanism. FAGCN [2] learns edge-level aggre

Figure 6: Comparison of SOTA models (magenta), selected baseline GNNs (red) and their ACM (green) and ACMII (blue) augmented models on 6 selected datasets. The black line and the error bar indicate the standard deviation. The symbol “\(\uparrow\)” means the amount of improvement of the best ACM-baseline and ACM-baseline over the SOTA models.
gation weights as GAT [58] but allows the weights to be negative which enables the network to capture the high-frequency components in graph signals. GPRGNN [6] uses learnable weights that can be both positive and negative for feature propagation, it allows GPRGNN to adapt heterophily structure of graph and is able to handle both high- and low-frequency parts of the graph signals.

### Future Work

Limitation of diversification operationDiversification operation does not work well in all harmful heterophily cases. For example, consider an imbalanced dataset where several small clusters with distinctive labels are densely connected to a large cluster. In this case, the surrounding differences of nodes in small clusters are similar, _i.e._ the neighborhood differences are mainly from their connection to the same large cluster, and this possibly makes diversification operation fail to discriminate them. Thus, it is obvious that ACM framework is not able to handle all heterophily cases.

From Figure 4, we can see that GNNs consistently perform well in the high homophily area. This reveals the fact that all homophily cases are helpful. This reminds us that instead of using a fixed adjacency matrix, we can learn a new adjacency matrix with different homophily level. With this in mind, we design an architecture with additional adjacency learner as shown in Figure 7: instead of using a fixed predefined adjacency matrix, we will learn an adjacency matrix with edges that can reveal the label similarity between nodes, _i.e._ homophily.. This adjacency learner should ideally be trained end-to end. From some preliminary experimental results (not included in this report) of a GCN with a pretrained adjacency learner, this method is promising although there are some stability issues need to be fixed.

Exploring Different Ways for Adjacency Candidate SelectionSome tricks can be explored when we are selecting the adjacency candidates for the adjacency learner:

* Sample or select (top-\(k_{1}\)) nodes from complementary graph, put them together with the predefined neighborhood set to form adjacency candidate set, then sample or select (top-\(k_{2}\)) adjacency candidates for training. Try to train it end-to-end.
* Consider modeling the candidate selection process as a multi-armed bandit problem. Find an efficient way to learn to select good candidates from complementary graph. Can use pseudo count to prevent selecting the same nodes repeatedly.

Figure 7: GNN with adjacency learner
Graph Representation Learning for Reinforcement Learning

### Markov Decision Process (MDP)

MDP is a framework to model the learning process that the agent learns from the interaction with the environment [56, 67, 68]. The interaction happens in discrete time steps, \(t=0,1,2,3,\cdots\). At step \(t\), given a state \(S_{t}=s_{t}\in\mathcal{S}\), the agent picks an action \(a_{t}\in\mathcal{A}(s_{t})\) according to a policy \(\pi(\cdot|s_{t})\), which is a rule of choosing actions given a state. Then, at time \(t+1\), the environmental dynamics \(p:\mathcal{S}\times\mathcal{R}\times\mathcal{A}\times\mathcal{S}\to[0,1]\) take the agent to a new state \(S_{t+1}=s_{t+1}\in\mathcal{S}\) and provide a numerical reward \(R_{t+1}=r_{t+1}(s_{t},a_{t},s_{t+1})\in\mathbb{R}\). Such a sequence of interaction gives us a trajectory \(\tau=\{S_{0},A_{0},R_{1},S_{1},A_{1},R_{2},S_{2},A_{2},R_{3},\cdots\}\). The objective is to find an optimal policy to maximize the expected long-term discounted cumulative reward \(V_{\pi}(s)=E_{\pi}[\sum\limits_{k=0}^{\infty}\gamma^{k}R_{t+k+1}|S_{t}=s]\) for each state \(s\), where \(\gamma\) is the discount factor.

For a given policy \(\pi\), solving its value function \(\mathbf{V}_{\pi}\) is equivalent to solving the following linear system,

\[\mathbf{V}_{\pi}=\mathbf{r}_{\pi}+\gamma P_{\pi}\mathbf{V}_{\pi} \tag{30}\]

where \(\mathbf{V}_{\pi}=[V_{\pi}(s)]_{s\in\mathcal{S}}^{T}\in\mathbb{R}^{|\mathcal{S}|}, \mathbf{r}_{\pi}=[r_{\pi}(s)]_{s\in\mathcal{S}}^{T}\in\mathbb{R}^{|\mathcal{S}|}, P_{\pi}=[P_{\pi}(s^{\prime}|s)]_{s^{\prime},s\in\mathcal{S}}\in\mathbb{R}^{| \mathcal{S}|\times|\mathcal{S}|}\). The state transition matrix \(P_{\pi}\) essentially defines a graph structure over states and the reward vector \(\mathbf{r}_{\pi}\) is a signal defines on graph. Thus, solving value function can be considered as a (supervised or semi-supervised) node regression tasks over graph. Besides solving \(\mathbf{V}_{\pi}\), the graph structure can also be used for reward propagation and representation learning in Reinforcement Learning (RL) [26, 27, 28].

### Graph Representation Learning for MDP

Treating MDP as a graph is an old but never outdated idea. Traditional methods use graph Laplacian for a fixed policy to estimate \(\mathbf{V}_{\pi}\), _e.g._ proto-value function [49]. In addition to value function estimation, [28] proposes to use GCN to learn potential-based reward shaping, which can accelerate the learning process of the agent.

The above methods both construct the graph from the sampled trajectory data. With modern Graph Representation Learning (GRL) methods _e.g._ node embedding methods [3], link prediction methods [52, 54, 57], we can learn to reconstruct the underlying graph (adjacency matrix) from sampled data more efficiently. And label propagation [53], which is a commonly used algorithm for graph semi-supervised learning, can be helpful for efficient reward propagation. In section 4.3, we will introduce the potential of using GRL for reward propagation and representation learning in reinforcement learning.

### Reinforcement Learning with Graph Representation Learning

In this section, we will draw how to represent Markov Decision Process (MDP) with graph and introduce two possible ways of using graph representation learning to address the problems defined on MDP.

Each state can be treated as a node on a graph, the transition probability between each pair of nodes (an element in state transition matrix) can be represented by the edge (or weight) between
them and value function is a function defined on each node of the graph. The details (for finite MDP) are introduced in matrix form as follows [38, 59]:

* Denote \(|S||A|\times|S|\) environment transition matrix as \(P\), where \[P_{(sa,s^{\prime})}=\sum_{r}p\left(s^{\prime},r|s,a\right)\] (31) and \(P_{(sa,s^{\prime})}\geq 0\), \(\sum_{s^{\prime}}P_{(sa,s^{\prime})}=1\), for all \(s,a\). Note that \(P\) is not a square matrix.
* We rewrite the policy \(\pi\) by an \(|S|\times|S||A|\) matrix \(\Pi\), where \(\Pi_{(s,s^{\prime}a)}=\pi(a|s)\) if \(s^{\prime}=s\), otherwise 0: \[\Pi=\text{diag}(\pi(\cdot|s_{1})^{T},\cdots,\pi(\cdot|s_{|S|})^{T})\] (32) where \(\pi(\cdot|s_{i})^{T}\) is an \(|A|\)-dimensional row vector. From this definition, one can quickly verify that the matrix product \(\Pi P\) gives the \(|S|\times|S|\) state-to-state transition matrix \(P_{\pi}\) (asymmetric) induced by the policy \(\pi\) in the environment \(P\), and the \(|S||A|\times|S||A|\) matrix product \(P\Pi\) gives the state-action-to-state-action transition matrix \(P^{\prime}_{\pi}\) (asymmetric) induced by policy \(\pi\) in the environment \(P\).
* We denote the \(|S||A|\times 1\) reward vector as \(\mathbf{r}\), whose entry \(r_{(sa)}\) specifies the reward obtained when taking action \(a\) in state \(s\), i.e. \[\mathbf{r}_{(sa)}=E[r|s,a]=\sum_{s^{\prime}\in\mathcal{S}}P_{sa,s^{\prime}}\cdot r (s,a,s^{\prime}).\] (33)
* The state value function and state-action value function can be represented by \[\mathbf{V}_{\pi}=\sum_{i=0}^{\infty}\gamma^{i}(\Pi P)^{i}\Pi\mathbf{r}=\Pi\mathbf{r}+ \gamma\Pi P\mathbf{V}_{\pi}\in\mathbb{R}^{|S|\times 1},\ \mathbf{Q}_{\pi}=\sum_{i=0}^{\infty} \gamma^{i}(P\Pi)^{i}\mathbf{r}=\mathbf{r}+\gamma P\Pi\mathbf{Q}_{\pi}\in\mathbb{R}^{|S|A| \times 1}\] (34)

#### 4.3.1 Learn Reward Propagation as Label Propagation

The sampling process from an MDP can be considered as a random walk defined on a graph, because the relation (edge) between each pair of states is essentially a transition probability 7. Discovering the underlying graph of a MDP can help us to leverage the correlation between states to learn value function or do to efficient exploration in sparse reward environment.

Footnote 7: From this perspective, we should not treat the trajectory as sequential data, because we do not necessarily have an ordered relation between states on a graph, even for directed graph. Although the observation seems to have an order in it, we actually only have transition relation.

Usually, the graph is constructed from the trajectory data, _i.e._ the pairwise state transition data. But once we update the policy, we need to reconstruct the graph. With graph embedding methods for link prediction tasks, _e.g._, Deepwalk [52], node2vec [14], Line [57], we are able to learn graph reconstruction by inferring some unobserved transition. To be more specifically, instead of learning \(P_{\pi}(s^{\prime}|s)\) for a fixed policy \(\pi\), we can learn the state-action transition probability \(P(s^{\prime}|s,a)\), which is independent of \(\pi\). In this way, we can take use of trajectory data in all history no matter the policy changes or not. And once we are given a policy, we can infer the graph by combining \(\pi(a|s)\) and \(P(s^{\prime}|s,a)\).
#### 4.3.2 Graph Embedding as Auxiliary Task for Representation Learning

Learning auxiliary tasks is showed to be helpful for state representation learning [22], which is critical to learn a good policy for agents. Among the methods, successor representation is showed to be theoretically and empirically important for learning a good state representations [8, 29]. Modeling the successor triplet \((s,a,s^{\prime})\) for MDP is essentially equivalent to modeling the triplet (head, relation, tail) in knowledge graph. And there exist a lot of algorithms in knowledge graph embedding community to address triplet embedding problem, _e.g._, TransE [3], RotateE [55], QuatE [65] and DihEdral [62]. These methods can be borrowed to learn richer representation for RL tasks.
