On Addressing the Limitations of Graph Neural Networks

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,34,50]. 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 Section 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 Contribution

In 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.

1.1. 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 $\left|\mathcal{V}\right|=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 =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_N)$ with ${\lambda }_1\le \cdots \le {\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}={[{\mathit{u}}_1^T\mathbf{x},\dots ,{\mathit{u}}_N^T\mathbf{x}]}^T$, where ${\mathit{u}}_i^T\mathbf{x}$ is the component of $\mathbf{x}$ in the direction of ${\mathit{u}}_{\mathbf{i}}$.

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

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

$$Y=\mathrm{softmax}\left({\widehat{A}}_{\mathrm{sym}}\mathrm{ReLU}\left({\widehat{A}}_{\mathrm{sym}}X{W}_0\right)W_1\right)$$

(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}=-\mathrm{trace}(Z^{T}logY).$$

(2)

The random walk renormalized matrix ${\widehat{A}}_{\mathrm{rw}}={\tilde{D}}^{-1}\tilde{A}$ can also be applied to GCN and it has the same eigenvalues as ${\widehat{A}}_{\mathrm{sym}}$. The corresponding Laplacian is defined as ${\widehat{L}}_{\mathrm{rw}}=I-{\widehat{A}}_{\mathrm{rw}}$ Specifically, the nature of random walk matrix makes ${\widehat{A}}_{\mathrm{rw}}$ behaves as a mean aggregator ${\left({\widehat{A}}_{\mathrm{rw}}\mathbf{x}\right)}_i={\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 [61,66]. 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.

In Section 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 Section 2.2, we propose to study an over-smoothing problem and give some ideas.

2.1. A Stronger Multi-Scale Deep GNN with Truncated Krylov Architecture

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

$$Y=\mathrm{softmax}\left({\widehat{A}}_{\mathrm{sym}}\mathrm{ReLU}(\cdots {\widehat{A}}_{\mathrm{sym}}\mathrm{ReLU}\left({\widehat{A}}_{\mathrm{sym}}\mathrm{ReLU}\left({\widehat{A}}_{\mathrm{sym}}X{W}_0\right)W_1\right)W_2\cdots )W_n\right)\equiv \mathrm{softmax}\left(Y^{\prime }\right)$$

(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 ${\widehat{A}}_{\mathrm{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 $\parallel W_j{\parallel }_2\le 1$ for $j=0,1,\dots$. If $\mathcal{G}$ has no bipartite components, then in (3), as $n\to \infty$, $\mathrm{rank}\left(Y^{\prime }\right)\le k$.

Theorem 2. 

Suppose the n-dimensional $\mathbf{x}$ and $\mathbf{y}$ are independently sampled from a continuous distribution and the activation function $\mathrm{Tan}\mathrm{h}\left(z\right)=\frac{e^z-e^{-z}}{e^z+e^{-z}}$ is applied to $[\mathbf{x},\mathbf{y}]$ pointwisely, then

$$\mathbb{P}(\mathrm{rank}\left(\mathrm{Tan}\mathrm{h}\left(\right[\mathbf{x},\mathbf{y}\left]\right)\right)=\mathrm{rank}\left([\mathbf{x},\mathbf{y}]\right))=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 1a 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 ${\widehat{A}}_{\mathrm{sym}}\in {\mathbb{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.

We take $\mathbb{S}={\mathbb{R}}^{F\times F}$. Given a set of block vectors ${\left\{X_k\right\}}_{k=1}^m\subset {\mathbb{R}}^{N\times F}$, the $\mathbb{S}$-span of ${\left\{X_k\right\}}_{k=1}^m$ is defined as ${\mathrm{span}}^{\mathbb{S}}\{X_1,\cdots ,X_m\}=\{\sum _{k=1}^m{X}_k{C}_k:C_k\in \mathbb{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 $\mathbb{S}$, and its corresponding block Krylov matrix are respectively defined as

$${\mathcal{K}}_m^{\mathbb{S}}(A,B)\equiv {\mathrm{span}}^{\mathbb{S}}\{B,AB,\cdots ,A^{m-1}B\},K_m(A,B)\equiv [B,AB,\cdots ,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\ge m$, $A^{k}B\in {\mathcal{K}}_m^{\mathbb{S}}(A,B)$, where m depends on A and B.

Let $\rho \left({\widehat{A}}_{\mathrm{sym}}\right)$ denote the spectrum radius of ${\widehat{A}}_{\mathrm{sym}}$ and suppose $\rho \left({\widehat{A}}_{\mathrm{sym}}\right)<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\left({\widehat{A}}_{\mathrm{sym}}\right)X=\sum _{n=0}^{\infty }\frac{g^{\left(n\right)}\left(0\right)}{n!}{\widehat{A}}_{\mathrm{sym}}^{n}X\equiv \left[X,{\widehat{A}}_{\mathrm{sym}}X,\cdots ,{\widehat{A}}_{\mathrm{sym}}^{m-1}X\right]{\left[{\left({{\Gamma }_0}^{\mathbb{S}}\right)}^T,{\left({{\Gamma }_1}^{\mathbb{S}}\right)}^T,\cdots ,{\left({\Gamma }_{m-1}^{\mathbb{S}}\right)}^T\right]}^T\equiv K_m({\widehat{A}}_{\mathrm{sym}},X){\Gamma }^{\mathbb{S}}$$

(4)

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

$$g\left({\widehat{A}}_{\mathrm{sym}}\right)X{W}^{\prime }=K_m({\widehat{A}}_{\mathrm{sym}},X){\Gamma }^{\mathbb{S}}{W}^{\prime }=K_m({\widehat{A}}_{\mathrm{sym}},X)W^{\mathbb{S}}$$

(5)

where $W^{\prime }\in {\mathbb{R}}^{F\times O}$ is a parameter matrix, and $W^{\mathbb{S}}\equiv {\Gamma }^{\mathbb{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].

Table 1. Accuracy without Validation. 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.

Table 1. Accuracy without Validation. 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.

2.2. Future Works on Over-Smoothing

2.2.0.2. Weight Initialization for GNNs

Even 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 $\widehat{A}$. To analyze the variance propagation, we use deep GCN as an example, use $\widehat{A}={\widehat{A}}_{\mathrm{rw}}$ and decompose it as follows,

$$\begin{array}{cc}\hfill & Y_0=X,H_1={\widehat{A}}_{\mathrm{rw}}X{W}_0,Y_1=f\left(H_1\right),H_{l+1}={\widehat{A}}_{\mathrm{rw}}{Y}_l{W}_l,Y_{l+1}=f\left(H_{l+1}\right),l=1,\cdots ,n\hfill \\ \hfill & Y=\mathrm{softmax}\left({\widehat{A}}_{\mathrm{rw}}{Y}_n{W}_n\right)\equiv \mathrm{softmax}\left(H_{n+1}\right),\mathcal{L}=-\mathrm{trace}\left(Z^T\log Y\right)\hfill \end{array}$$

(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,

$$\begin{array}{cc}\hfill & \frac{\partial \mathcal{L}}{\partial H_l}=\frac{\partial \mathcal{L}}{\partial Y_l}\odot f^{\prime }\left(H_l\right),\frac{\partial \mathcal{L}}{\partial W_{l-1}}=Y_{l-1}^T{\widehat{A}}_{\mathrm{rw}}\frac{\partial \mathcal{L}}{\partial H_l},\frac{\partial \mathcal{L}}{\partial Y_{l-1}}={\widehat{A}}_{\mathrm{rw}}\frac{\partial \mathcal{L}}{\partial H_l}{W}_{l-1}^T\hfill \end{array}$$

(7)

2.2.0.3. Variance Analysis: Forward View

Consider element $i,j$ in matrix $H_{l+1}$ during the feed-forward process in (6),

$${\left(H_{l+1}\right)}_{ij}={\left({\widehat{A}}_{\mathrm{rw}}\right)}_{i,:}{Y}_l{\left(W_l\right)}_{:,j}=\sum _{t=1}^{F_l}\sum _{k=1}^N{\left({\widehat{A}}_{\mathrm{rw}}\right)}_{ik}{\left(Y_l\right)}_{kt}{\left(W_l\right)}_{t,j},Y_{l+1}=f\left(H_{l+1}\right),l=1,\cdots ,n$$

(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({\left(W_l\right)}_{ij}\right)=0$; $E\left({\left(Y_l\right)}_{kt}\right)=0$ and all elements in $Y_l$ are independent 2. Then, $\mathrm{Var}\left({\left(Y_{l+1}\right)}_{ij}\right)=\mathrm{Var}\left(\left(Y_{l+1}\right)\right)$ can be written as

$$\begin{array}{cc}\hfill & \mathrm{Var}\left(\sum _{t=1}^{F_l}\sum _{k=1}^N{\left({\widehat{A}}_{\mathrm{rw}}\right)}_{ik}{\left(Y_l\right)}_{kt}{\left(W_l\right)}_{t,j}\right)=\sum _{t=1}^{F_l}\sum _{k=1}^N\mathrm{Var}\left({\left({\widehat{A}}_{\mathrm{rw}}\right)}_{ik}{\left(Y_l\right)}_{kt}{\left(W_l\right)}_{t,j}\right)=\frac{F_l}{d_i+1}\mathrm{Var}\left(Y_l\left)\mathrm{Var}\right(W_l\right)\hfill \end{array}$$

(9)

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

$$\mathrm{Var}\left(W_l\right)=\frac{d_i+1}{F_l}$$

(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

$$\mathrm{Var}\left(W_l\right)\approx \frac{{\displaystyle \sum _{i=1}^N}(d_i+1)}{N{F}_l}=\frac{1+\mathrm{average}\mathrm{node}\mathrm{degree}}{F_l}$$

(11)

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

$$\mathrm{Var}\left(W\right)={\displaystyle \sum _{i=1}^N}\frac{d_i+1}{{\displaystyle \sum _{j=1}^N}{d}_j+1}\frac{d_i+1}{F_l}=\frac{{\displaystyle \sum _{i=1}^N}{(d_i+1)}^2}{({\displaystyle \sum _{i=1}^N}{d}_i+1)F_l}$$

(12)

Variance Analysis: Backward View

Under 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)

$$\begin{array}{c}\hfill \frac{\partial \mathcal{L}}{\partial H_l}=\frac{\partial \mathcal{L}}{\partial Y_l}={\widehat{A}}_{\mathrm{rw}}\frac{\partial \mathcal{L}}{\partial H_{l+1}}{W}_l^T,\frac{\partial \mathcal{L}}{\partial W_{l-1}}=Y_{l-1}^T{\widehat{A}}_{\mathrm{rw}}\frac{\partial \mathcal{L}}{\partial H_l}\end{array}$$

(13)

Then,

$$\begin{array}{cc}\hfill {\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)}_{ij}& =\sum _{t=1}^{F_{l+1}}\sum _{k=1}^N{\left({\widehat{A}}_{\mathrm{rw}}\right)}_{ik}{\left(\frac{\partial \mathcal{L}}{\partial H_{l+1}}\right)}_{kt}{\left(W_l^T\right)}_{t,j}\hfill \\ \hfill {\left(\frac{\partial \mathcal{L}}{\partial W_{l-1}}\right)}_{ij}& ={\left({\widehat{A}}_{\mathrm{rw}}{Y}_{l-1}\right)}_{·i}^T{\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)}_{·j}=\sum _{k=1}^N\left(\sum _{t=1}^N{\left({\widehat{A}}_{\mathrm{rw}}\right)}_{kt}{\left(Y_{l-1}\right)}_{ti}\right){\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)}_{kj},\hfill \end{array}$$

(14)

Thus,

$$\begin{array}{cc}\hfill \mathrm{Var}\left({\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)}_{ij}\right)& =\mathrm{Var}\left(\sum _{t=1}^{F_{l+1}}\sum _{k=1}^N{\widehat{A}}_{{\mathrm{rw}}_{ik}}{\left(\frac{\partial \mathcal{L}}{\partial H_{l+1}}\right)}_{kt}{\left(W_l^T\right)}_{t,j}\right)=\frac{F_{l+1}}{d_i+1}\mathrm{Var}\left(\frac{\partial \mathcal{L}}{\partial H_{l+1}}\right)\mathrm{Var}\left(W_l\right)\hfill \\ \hfill \mathrm{Var}\left({\left(\frac{\partial \mathcal{L}}{\partial W_{l-1}}\right)}_{ij}\right)& =\mathrm{Var}\left(\sum _{k=1}^N\left(\sum _{t=1}^N{\left({\widehat{A}}_{\mathrm{rw}}\right)}_{kt}{\left(Y_{l-1}\right)}_{ti}\right){\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)}_{kj}\right)=\left(\sum _{k=1}^N\frac{1}{d_k+1}\right)\mathrm{Var}\left(Y_{l-1}\right)\mathrm{Var}\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)\hfill \end{array}$$

(15)

$$\mathrm{Var}\left(W_l\right)=\frac{{\displaystyle \sum _{i=1}^N}(d_i+1)}{N{F}_{l+1}}\approx \frac{1+\mathrm{average}\mathrm{node}\mathrm{degree}}{F_{l+1}}$$

(16)

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

$$\begin{array}{cc}\hfill & \mathrm{Var}\left(\frac{\partial \mathcal{L}}{\partial H_l}\right)\approx \frac{N{F}_{l+1}}{\sum _{i=1}^N{d}_i+1}\mathrm{Var}\left(\frac{\partial \mathcal{L}}{\partial H_{l+1}}\right)\mathrm{Var}\left(W_l\right)\approx \mathrm{Var}\left(\frac{\partial \mathcal{L}}{\partial H_{n+1}}\right)\prod _{l^{\prime }=l+1}^{n+1}\frac{N{F}_{l^{\prime }}}{\sum _{i=1}^N{d}_i+1}\mathrm{Var}\left(W_{l^{\prime }-1}\right)\hfill \\ \hfill & \mathrm{Var}\left(Y_{l-1}\right)\approx \frac{N{F}_{l-2}}{{\sum }_k{d}_k+1}\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}-\mathbf{2}}\left)\mathrm{Var}\right(W_{l-2}\right)\approx \mathrm{Var}\left({\mathit{Y}}_{\mathbf{0}}\right)\prod _{l^{\prime }=0}^{l-2}\frac{N{F}_{l^{\prime }}}{{\sum }_k{d}_k+1}\mathrm{Var}\left(W_{l^{\prime }}\right)\hfill \end{array}$$

(17)

From (15), if each $\mathrm{Var}\left(W_{l^{\prime }}\right)$ equals to $\mathrm{Var}\left(W\right)$ and each $F_l$ equals to F, then

$$\begin{array}{cc}\hfill & \mathrm{Var}\left({\left(\frac{\partial \mathcal{L}}{\partial W_{l-1}}\right)}_{ij}\right)\approx \left(\sum _{k=1}^N\frac{1}{d_k+1}\right)\mathrm{Var}\left({\mathit{Y}}_{\mathbf{0}}\right)\mathrm{Var}\left(\frac{\partial \mathcal{L}}{\partial H_{n+1}}\right){\left(\frac{NF}{{\sum }_k{d}_k+1}\mathrm{Var}\left(W\right)\right)}^n\hfill \end{array}$$

(18)

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

$$\mathrm{Var}\left(W_l\right)\approx \frac{2{\displaystyle \sum _{i=1}^N}(d_i+1)}{N(F_l+F_{l+1})}=\frac{2(1+\mathrm{average}\mathrm{node}\mathrm{degree})}{(F_l+F_{l+1})}$$

(19)

Thus, each element in $W_i$ can be drawn from $N(0,\sqrt{\frac{2(1+\mathrm{average}\mathrm{node}\mathrm{degree})}{(F_l+F_{l+1})}})$.

Adaptive ReLU (AdaReLU) Activation Function

To 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\left({\mathbf{x}}_i\right)=\left\{\begin{array}{cc}{\beta }_i{\mathbf{x}}_i,\hfill & \mathrm{if}{\mathbf{x}}_i>0\hfill \\ {\alpha }_i{\mathbf{x}}_i,\hfill & \mathrm{if}{\mathbf{x}}_i\le 0\hfill \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,42,64,69,70].

In this section , we first introduce the most commently used homophily metrics in Section 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 Section 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 Section 3.3. With this fact and knowledge of filterbanks, we propose the Adaptive Channel Mixing (ACM) framework in Section 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 Section 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 Section 3.6

3.1. 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 [1,70], node homophily [51], and class homophily [35] 4 defined as follows:

$$\begin{array}{cc}\hfill & H_{\mathrm{edge}}\left(\mathcal{G}\right)=\frac{\left|\right\{e_{uv}\mid e_{uv}\in \mathcal{E},Z_{u,:}=Z_{v,:}\left\}\right|}{\left|\mathcal{E}\right|},H_{\mathrm{node}}\left(\mathcal{G}\right)=\frac{1}{\left|\mathcal{V}\right|}\sum _{v\in \mathcal{V}}\frac{\left|\right\{u\mid u\in {\mathcal{N}}_v,Z_{u,:}=Z_{v,:}\left\}\right|}{d_v},\hfill \\ \hfill & H_{\mathrm{class}}\left(\mathcal{G}\right)=\frac{1}{C-1}\sum _{k=1}^C{\left[h_k-\frac{\left|\right\{v\mid Z_{v,k}=1\left\}\right|}{N}\right]}_{+},h_k=\frac{{\sum }_{v\in \mathcal{V}}\left|\{u\mid Z_{v,k}=1,u\in {\mathcal{N}}_v,Z_{u,:}=Z_{v,:}\}\right|}{{\sum }_{v\in \left\{v\right|Z_{v,k}=1\}}{d}_v}\hfill \end{array}$$

(20)

where ${\left[a\right]}_{+}=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_{\mathrm{edge}}\left(\mathcal{G}\right)$ measures the proportion of edges that connect two nodes in the same class; $H_{\mathrm{node}}\left(\mathcal{G}\right)$ evaluates the average proportion of edge-label consistency of all nodes; $H_{\mathrm{class}}\left(\mathcal{G}\right)$ tries to avoid the sensitivity to imbalanced class, which can cause $H_{\mathrm{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.

3.2. 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_{\mathrm{node}}$ by examples to show that different graph typologies with the same $H_{\mathrm{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 $\widehat{A}\in {\mathbb{R}}^{N\times N}$ denote a general aggregation operator. Then, Equation (1) can be simplified as,

$$\begin{array}{c}\hfill Y=\mathrm{softmax}\left(\widehat{A}XW\right)=\mathrm{softmax}\left(Y^{\prime }\right)\end{array}$$

(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,

$$\begin{array}{c}\hfill \Delta Y^{\prime }=\widehat{A}X\Delta W=\gamma \widehat{A}X\frac{d\mathcal{L}}{dW}\propto \widehat{A}X\frac{d\mathcal{L}}{dW}=\widehat{A}X{X}^T{\widehat{A}}^T(Z-Y)=S(\widehat{A},X)(Z-Y)\end{array}$$

(22)

where $S(\widehat{A},X)\equiv \widehat{A}X{\left(\widehat{A}X\right)}^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 {\left(Y^{\prime }\right)}_{i,:}={\sum }_{j\in \mathcal{V}}{\left[S(\widehat{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_{\mathrm{agg}}\left(S(\widehat{A},X)\right)=\frac{\left|\left\{v|{\mathrm{Mean}}_u\left(\left\{S{(\widehat{A},X)}_{v,u}\right|Z_{u,:}=Z_{v,:}\}\right)\ge {\mathrm{Mean}}_u\left(\left\{S{(\widehat{A},X)}_{v,u}\right|Z_{u,:}\ne Z_{v,:}\}\right)\right\}\right|}{\left|\mathcal{V}\right|}$$

(23)

where ${\mathrm{Mean}}_u\left(\{·\}\right)$ takes the average over u of a given multiset of values or variables.

$S_{\mathrm{agg}}\left(S(\widehat{A},X)\right)$ 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_{\mathrm{agg}}\left(S(\widehat{A},X)\right)\in [0,1]$. But in practice, we observe that in most datasets, we will have $S_{\mathrm{agg}}\left(S(\widehat{A},X)\right)\ge 0.5$. Based on this observation, we rescale (23) to the following modified aggregation similarity for practical usage,

$$S_{\mathrm{agg}}^M\left(S(\widehat{A},X)\right)={\left[2S_{\mathrm{agg}}\left(S(\widehat{A},X)\right)-1\right]}_{+}$$

(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 ($\widehat{A}$) homophily and its modified version as

$$H_{\mathrm{agg}}\left(\mathcal{G}\right)=S_{\mathrm{agg}}\left(S(\widehat{A},Z)\right),H_{\mathrm{agg}}^M\left(\mathcal{G}\right)=S_{\mathrm{agg}}^M\left(S(\widehat{A},Z)\right)$$

(25)

In practice, we will only check $H_{\mathrm{agg}}\left(\mathcal{G}\right)$ when $H_{\mathrm{agg}}^M\left(\mathcal{G}\right)=0$. As Figure 3 shows, when $\widehat{A}={\widehat{A}}_{\mathrm{rw}}$, $H_{\mathrm{agg}}\left(\mathcal{G}\right)=H_{\mathrm{agg}}^M\left(\mathcal{G}\right)=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_{\mathrm{agg}}\left(\mathcal{G}\right)$ and $H_{\mathrm{agg}}^M\left(\mathcal{G}\right)$ by additionally considering information from aggregation operator $\widehat{A}$ and the similarity matrix.

Comparison of Homophily Metrics on Synthetic Graphs

To comprehensively compare $H_{\mathrm{agg}}^M\left(\mathcal{G}\right)$ 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 4a–c show that the performance curves under $H_{\mathrm{edge}}\left(\mathcal{G}\right),H_{\mathrm{node}}\left(\mathcal{G}\right)$ and $H_{\mathrm{class}}\left(\mathcal{G}\right)$ are U-shaped 5, while Figure 2d reveals a nearly monotonic curve only with a little numerical perturbation around 1. This indicates that $H_{\mathrm{agg}}^M\left(\mathcal{G}\right)$ can describe how the graph structure affects the performance of SGC-1 and GCN more appropriately and adequately than the existing metrics.

3.3. How Diversification Operation Helps with Harmful Heterophily

We first consider the example shown in Figure 5. From $S(\widehat{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-\widehat{A}$[10] (will be introduced in the next subsection) to extract the information of neighborhood differences and address harmful heterophily. As $S(I-\widehat{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-\widehat{A},X)$.

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

$$\begin{array}{c}\hfill \mathbf{1}.{\mathrm{Mean}}_u\left(\left\{S{(I-\widehat{A},X)}_{v,u}\right|u\in \mathcal{V}\wedge Z_{u,:}=Z_{v,:}\}\right)>0;\\ \hfill \mathbf{2}.{\mathrm{Mean}}_u\left(\left\{S{(I-\widehat{A},X)}_{v,u}\right|u\in \mathcal{V}\wedge Z_{u,:}\ne Z_{v,:}\}\right)\le 0\end{array}$$

(26)

Then, graph diversification distinguishability value is defined as

$${\mathrm{DD}}_{\widehat{A},X}\left(\mathcal{G}\right)=\frac{1}{\left|\mathcal{V}\right|}\left|\left\{v\right|visdiversificationdistinguishable\}\right|$$

(27)

We can see that ${\mathrm{DD}}_{\widehat{A},X}\left(\mathcal{G}\right)\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,\widehat{A}={\widehat{A}}_{\mathrm{rw}}$. Then, for a binary classification problem, i.e. $C=2$, all nodes are diversification distinguishable, i.e.${\mathrm{DD}}_{\widehat{A},Z}\left(\mathcal{G}\right)=1$. 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.

3.4. Filterbank and Adaptive Channel Mixing(ACM) GNN Framework

Filterbank

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

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

Filterbank in Spatial Form

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

$${\left(H_{\mathrm{LP}}\mathbf{x}\right)}_i=\sum _{j\in \{{\mathcal{N}}_i\cup i\}}{\widehat{A}}_{i,j}{\mathbf{x}}_{\mathbf{j}},{\left(H_{\mathrm{HP}}\mathbf{x}\right)}_i={\mathbf{x}}_{\mathbf{i}}-\sum _{j\in \{{\mathcal{N}}_i\cup i\}}{\widehat{A}}_{i,j}{\mathbf{x}}_{\mathbf{j}}$$

(28)

where ${\widehat{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_{\mathrm{LP}}$ and $H_{\mathrm{HP}}$ to represent general LP and HP filters. The ACM framework includes 3 steps as follows,

$$\begin{array}{cc}\hfill & \mathbf{Step}\mathbf{1}.\mathbf{Feature}\mathbf{Extraction}\mathbf{for}\mathbf{Each}\mathbf{Channel}:\hfill \\ \hfill & \mathrm{Option}1\text{:}{H}_L^l=\mathrm{ReLU}\left(H_{\mathrm{LP}}{H}^{l-1}{W}_L^{l-1}\right),H_H^l=\mathrm{ReLU}\left(H_{\mathrm{HP}}{H}^{l-1}{W}_H^{l-1}\right),H_I^l=\mathrm{ReLU}\left(I{H}^{l-1}{W}_I^{l-1}\right);\hfill \\ \hfill & \mathrm{Option}2\text{:}{H}_L^l=H_{\mathrm{LP}}\mathrm{ReLU}\left(H^{l-1}{W}_L^{l-1}\right),H_H^l=H_{\mathrm{HP}}\mathrm{ReLU}\left(H^{l-1}{W}_H^{l-1}\right),H_I^l=I\mathrm{ReLU}\left(H^{l-1}{W}_I^{l-1}\right);\hfill \\ \hfill & W_L^{l-1},W_H^{l-1},W_I^{l-1}\in {\mathbb{R}}^{F_{l-1}\times F_l};\hfill \\ \hfill & \mathbf{Step}\mathbf{2}.\mathbf{Feature}-\mathbf{based}\mathbf{Weight}\mathbf{Learning}\hfill \\ \hfill & {\tilde{\alpha }}_L^l=\sigma \left(H_L^l{\tilde{W}}_L^l\right),{\tilde{\alpha }}_H^l=\sigma \left(H_H^l{\tilde{W}}_H^l\right),{\tilde{\alpha }}_I^l=\sigma \left(H_I^l{\tilde{W}}_I^l\right),{\tilde{W}}_L^{l-1},{\tilde{W}}_H^{l-1},{\tilde{W}}_I^{l-1}\in {\mathbb{R}}^{F_l\times 1}\hfill \\ \hfill & \left[{\alpha }_L^l,{\alpha }_H^l,{\alpha }_I^l\right]=\mathrm{Softmax}\left(\left[{\tilde{\alpha }}_L^l,{\tilde{\alpha }}_H^l,{\tilde{\alpha }}_I^l\right]W_{\mathrm{Mix}}^l/T,\right),W_{\mathrm{Mix}}^l\in {\mathbb{R}}^{3\times 3},T\in \mathbb{R}\mathrm{is}\mathrm{the}\mathrm{temperature};\hfill \\ \hfill & \mathbf{Step}\mathbf{3}.\mathbf{Node}-\mathbf{wise}\mathbf{Channel}\mathbf{Mixing}:\hfill \\ \hfill & H^l=\left(\mathrm{diag}\left({\alpha }_L^l\right)H_L^l+\mathrm{diag}\left({\alpha }_H^l\right)H_H^l+\mathrm{diag}\left({\alpha }_I^l\right)H_I^l\right).\hfill \end{array}$$

(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_{\mathrm{LP}}$ and $H_{\mathrm{HP}}$ by aggregation and diversification operations as shown in (28).

Complexity

Number 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 $N{F}_l(8+6F_{l-1})+2F_l(\mathrm{nnz}\left(H_{\mathrm{LP}}\right)+\mathrm{nnz}\left(H_{\mathrm{HP}}\right))+18N$ flops, while GCN layer takes $2N{F}_{l-1}{F}_l+2F_l\left(\mathrm{nnz}\left(H_{\mathrm{LP}}\right)\right)$ flops, where $\mathrm{nnz}(·)$ is the number of non-zero elements.

Performance Comparison

We 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 ${\widehat{A}}_{\mathrm{rw}}$ as LP filter and the corresponding HP filter is $I-{\widehat{A}}_{\mathrm{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.

3.5. 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 aggregation 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 GRPGNN to adapt heterophily structure of graph and is able to handle both high- and low-frequency parts of the graph signals.

3.6. Future Work

Limitation of Diversification Operation

Diversification 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.

3.6.0.12. Exploring Different Ways for Adjacency Candidate Selection

Some 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 pre-defined 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.

4. Graph Representation Learning for Reinforcement Learning

4.1. 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}\left(s_t\right)$ according to a policy $\pi (·|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 }\left(s\right)=E_{\pi }\left[\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}\right|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 }$$

(30)

where ${\mathbf{V}}_{\pi }={\left[V_{\pi }\left(s\right)\right]}_{s\in \mathcal{S}}^T\in {\mathbb{R}}^{\left|\mathcal{S}\right|},{\mathbf{r}}_{\pi }={\left[r_{\pi }\left(s\right)\right]}_{s\in \mathcal{S}}^T\in {\mathbb{R}}^{\left|\mathcal{S}\right|},P_{\pi }={\left[P_{\pi }\left(s^{\prime }\right|s)\right]}_{s^{\prime },s\in \mathcal{S}}\in {\mathbb{R}}^{\left|\mathcal{S}\right|\times \left|\mathcal{S}\right|}$. 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].

4.2. 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.

4.3. 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 $\left|S\right|\left|A\right|\times \left|S\right|$ environment transition matrix as P, where

  $$P_{\left(sa,s^{\prime }\right)}=\sum _{r}p\left(s^{\prime },r|s,a\right)$$

  (31)

  and $P_{\left(sa,s^{\prime }\right)}\ge 0$, ${\sum }_{s^{\prime }}{P}_{\left(sa,s^{\prime }\right)}=1$, for all $s,a$. Note that P is not a square matrix.
• We rewrite the policy $\pi$ by an $\left|S\right|\times \left|S\right|\left|A\right|$ matrix $\Pi$, where ${\Pi }_{\left(s,s^{\prime }a\right)}=\pi \left(a\right|s)\mathrm{if}{s}^{\prime }=s$, otherwise 0:

  $$\Pi =\mathrm{diag}(\pi {(·|s_1)}^T,\cdots ,\pi {(·|s_{\left|S\right|})}^T)$$

  (32)

  where $\pi {(·|s_i)}^T$ is an $\left|A\right|$-dimensional row vector. From this definition, one can quickly verify that the matrix product $\Pi P$ gives the $\left|S\right|\times \left|S\right|$ state-to-state transition matrix $P_{\pi }$ (asymmetric) induced by the policy $\pi$ in the environment P, and the $\left|S\right|\left|A\right|\times \left|S\right|\left|A\right|$ matrix product $P\Pi$ gives the state-action-to-state-action transition matrix $P_{\pi }^{\prime }$ (asymmetric) induced by policy $\pi$ in the environment P.
• We denote the $\left|S\right|\left|A\right|\times 1$ reward vector as $\mathbf{r}$, whose entry $r_{\left(sa\right)}$ specifies the reward obtained when taking action a in state s, i.e.

  $${\mathbf{r}}_{\left(sa\right)}=E\left[r\right|s,a]=\sum _{s^{\prime }\in \mathcal{S}}{P}_{sa,s^{\prime }}·r(s,a,s^{\prime }).$$

  (33)
• The state value function and state-action value function can be represented by

  $${\mathbf{V}}_{\mathbf{\pi }}=\sum _{i=0}^{\infty }{\gamma }^i{\left(\Pi P\right)}^i\Pi \mathbf{r}=\Pi \mathbf{r}+\gamma \Pi P{\mathbf{V}}_{\mathbf{\pi }}\in {\mathbb{R}}^{\left|S\right|\times 1},{\mathbf{Q}}_{\mathbf{\pi }}=\sum _{i=0}^{\infty }{\gamma }^i{\left(P\Pi \right)}^i\mathbf{r}=\mathbf{r}+\gamma P\Pi {\mathbf{Q}}_{\mathbf{\pi }}\in {\mathbb{R}}^{\left|S\right|\left|A\right|\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.

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 }\left(s^{\prime }\right|s)$ for a fixed policy $\pi$, we can learn the state-action transition probability $P\left(s^{\prime }\right|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 \left(a\right|s)$ and $P\left(s^{\prime }\right|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 $\left(\mathrm{head},\mathrm{relation},\mathrm{tail}\right)$ 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.

Appendix A.1. Background

We first decompose the deep GCN architecture as follows

$$\begin{array}{cc}\hfill & {\mathit{Y}}_{\mathbf{0}}=\mathit{X},{\mathit{H}}_{\mathbf{1}}=\widehat{A}\mathit{X}{W}_0,{\mathit{Y}}_{\mathbf{1}}=f\left({\mathit{H}}_{\mathbf{1}}\right)\hfill \\ \hfill & {\mathit{H}}_{\mathbf{l}+\mathbf{1}}=\widehat{A}{\mathit{Y}}_{\mathbf{l}}{W}_l,{\mathit{Y}}_{\mathbf{l}+\mathbf{1}}=f\left({\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right),l=1,\cdots ,n\hfill \\ \hfill & \mathit{Y}=\mathrm{softmax}\left(\widehat{A}{\mathit{Y}}_{\mathbf{n}}{W}_n\right)\equiv \mathrm{softmax}\left({\mathit{H}}_{\mathbf{n}+\mathbf{1}}\right)\hfill \\ \hfill & L=-\mathrm{trace}\left({\mathbf{Z}}^T\log \mathit{Y}\right)\hfill \end{array}$$

(A1)

where ${\mathit{H}}_{\mathbf{l}},{\mathit{Y}}_{\mathbf{l}}\in {\mathbb{R}}^{N\times F_l}$, $W_l\in {\mathbb{R}}^{F_l\times F_{l+1}}$; $\mathbf{Z}\in {\mathbb{R}}^{N\times C}$ is the ground truth matrix with one-hot label vector ${\mathbf{Z}}_{i,:}$ in each row, C is number of classes; L is the scalar loss. Then the gradient propagates in the following way

$$\begin{array}{cc}\hfill \mathrm{Output}& \frac{\partial L}{\partial {\mathit{H}}_{\mathbf{n}+\mathbf{1}}}=\mathrm{softmax}\left({\mathit{H}}_{\mathbf{n}+\mathbf{1}}\right)-\mathbf{Z}\hfill \\ \hfill & \frac{\partial L}{\partial W_n}={{\mathit{Y}}_{\mathbf{n}}}^T\widehat{A}\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{n}+\mathbf{1}}},\frac{\partial L}{\partial {\mathit{Y}}_{\mathbf{n}}}=\widehat{A}\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{n}+\mathbf{1}}}{W}_n^T\hfill \\ \hfill \mathrm{Hidden}& \frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}=\frac{\partial L}{\partial {\mathit{Y}}_{\mathbf{l}}}\odot f^{\prime }\left({\mathit{H}}_{\mathbf{l}}\right),\frac{\partial L}{\partial W_{l-1}}={{\mathit{Y}}_{\mathbf{l}-\mathbf{1}}}^T\widehat{A}\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}},\hfill \\ \hfill & \frac{\partial L}{\partial {\mathit{Y}}_{\mathbf{l}-\mathbf{1}}}=\widehat{A}\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}{W}_{l-1}^T\hfill \end{array}$$

(A2)

where ⊙ is the Hadamard product. The gradient propagation of GCN differs from that of multi-layer perceptron (MLP) by an extra multiplication of $\widehat{A}$ when the gradient signal flows through ${\mathit{Y}}_{\mathbf{l}}$.

Appendix A.2. Variance Analysis

Appendix A.2.1. Forward View

$$\begin{array}{cc}\hfill & {\mathit{H}}_{\mathbf{l}+\mathbf{1}}=\widehat{A}{\mathit{Y}}_{\mathbf{l}}{W}_l,{\left({\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right)}_{ij}={\widehat{A}}_{i,:}{\mathit{Y}}_{\mathbf{l}}{\left(W_l\right)}_{:,j}=\sum _{t=1}^{F_l}\sum _{k=1}^N{\widehat{A}}_{ik}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_{kt}{\left(W_l\right)}_{t,j}\hfill \\ \hfill & {\mathit{Y}}_{\mathbf{l}+\mathbf{1}}=f\left({\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right),l=1,\cdots ,n\hfill \end{array}$$

(A3)

Suppose the activation function is identity function such as [60], all element in W share the same variance and each element in $\mathit{X}$ are independent and share the same variance, $E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_{kt}\right)=0$ and $E\left({\left(W_l\right)}_{ij}\right)=0$, $\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}+\mathbf{1}}\right)}_{ij}\right)$ can be written as

$$\begin{array}{cc}\hfill & \mathrm{Var}\left(\sum _{t=1}^{F_l}\sum _{k=1}^N{\widehat{A}}_{ik}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_{kt}{W}_{t,j}\right)=\sum _{t=1}^{F_l}\sum _{k=1}^N\mathrm{Var}\left({\widehat{A}}_{ik}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_{kt}{W}_{t,j}\right)\hfill \\ \hfill & =F_l(d_i+1)·\frac{1}{{(d_i+1)}^2}\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}}\right)\mathrm{Var}\left(W\right)\hfill \\ \hfill & =\frac{F_l}{d_i+1}\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}}\left)\mathrm{Var}\right(W\right)=\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}+\mathbf{1}}\right)\hfill \end{array}$$

(A4)

Then, if we want $\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}+\mathbf{1}}\right)=\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}}\right)$ we will have

$$\mathrm{Var}\left(W\right)=\frac{d_i+1}{F_l}$$

(A5)

Since the parameter matrix is shared by all nodes, we cannot design a node specified initialization scheme for the W. Thus, we initialize each element in W by the average values as follows

$$\mathrm{Var}\left(W\right)=\frac{\sum _{i=1}^N(d_i+1)}{N{F}_l}=\frac{1+\mathrm{average}\mathrm{node}\mathrm{degree}}{F_l}$$

(A6)

If we relax the assumption $E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_{kt}\right)=0$ and assume it nonzero, as shown in [18], if we use ReLU activation function, we will have

$$\mathrm{Var}\left(W\right)=\frac{2\sum _{i=1}^N(d_i+1)}{N{F}_l}=\frac{2(1+\mathrm{average}\mathrm{node}\mathrm{degree})}{F_l}$$

(A7)

If we consider the correlation between nodes and keep the assumption that each dimension of the feature (columns in $\mathit{X}$) is independent and $E\left(X\right)\ne 0$, we will have

$$\begin{array}{cc}\hfill & \mathrm{Var}\left({\left({\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right)}_{ij}\right)=\mathrm{Var}\left(\sum _{t=1}^{F_l}\sum _{k=1}^N{\widehat{A}}_{ik}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_{kt}{W}_{t,j}\right)\hfill \\ \hfill & =F_l·\mathrm{Var}\left(\left(\sum _{k=1}^N{A}_{ik}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)W\right)=F_l·E{\left(\sum _{k=1}^N{A}_{ik}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)}^2\mathrm{Var}\left(W\right)\hfill \\ \hfill & =\frac{F_l}{{(d_i+1)}^2}E{\left(\sum _{k\in {\mathcal{N}}_i}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)}^2\mathrm{Var}\left(W\right)\hfill \\ \hfill & =\frac{F_l}{{(d_i+1)}^2}E\left(\sum _{k\in {\mathcal{N}}_i}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2+\sum _{k,j\in {\mathcal{N}}_i,k\ne j}{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)\mathrm{Var}\left(W\right)\hfill \end{array}$$

(A8)

Since

$$E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)=\mathrm{Cov}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k,{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)+E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)$$

we can have several reasonable assumptions over $\mathrm{Cov}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k,{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)$ to get different results

• The adjacency matrix with self-loop can be considered as a prior covariance matrix and thus a reasonable assumption is $\mathrm{Cov}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k,{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)=\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)=\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)$.
• Consider symmetric normalized $\widehat{A}$ as a prior covariance matrix and we have $\mathrm{Cov}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k,{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)=\sqrt{\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)}=\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)$

These assumptions all lead us to

$$E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)=\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)+E^2\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)=E{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2$$

Thus we have

$$\begin{array}{cc}\hfill & \mathrm{Var}\left({\left({\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right)}_{ij}\right)=F_l·\frac{1}{{(d_i+1)}^2}E\left({(d_i+1)}^2{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2\right)\mathrm{Var}\left(W\right)\hfill \\ \hfill & =F_l·E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2\right)\mathrm{Var}\left(W\right)=F_l·\frac{1}{2}\mathrm{Var}\left({\mathit{H}}_l\right)\mathrm{Var}\left(W\right)\hfill \end{array}$$

(A9)

Thus

$$\begin{array}{c}\hfill \mathrm{Var}\left(W\right)=\frac{2}{F_l}\end{array}$$

(A10)

Suppose $E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)\ge 0$ and $Cov({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k,{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j)\ge 0$, since $Cov({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k,{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j)\le \sqrt{\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)}=\mathrm{Var}\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k\right)$, we have $E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)\le E{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2$ and we assume $E\left({\left({\mathit{Y}}_{\mathbf{l}}\right)}_k{\left({\mathit{Y}}_{\mathbf{l}}\right)}_j\right)=\alpha E{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2,\alpha \in [0,1]$. Then,

$$\begin{array}{cc}\hfill & \mathrm{Var}\left({\left({\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right)}_{ij}\right)=F_l·\frac{1}{{(d_i+1)}^2}(d_i+1)(\alpha d_i+1)E{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2\mathrm{Var}\left(W\right)\hfill \\ \hfill & =F_l·\frac{\alpha d_i+1}{d_i+1}E{\left({\mathit{Y}}_{\mathbf{l}}\right)}_k^2\mathrm{Var}\left(W\right)\hfill \\ \hfill & =F_l·\frac{\alpha d_i+1}{2(d_i+1)}\mathrm{Var}\left({\mathit{H}}_{\mathbf{l}}\right)\mathrm{Var}\left(W\right)\hfill \end{array}$$

(A11)

Thus,

$$\begin{array}{cc}\hfill & \mathrm{Var}\left(W\right)=\frac{2(d_i+1)}{F_l(\alpha d_i+1)}\hfill \end{array}$$

(A12)

$(d_i+1)/(\alpha d_i+1)$ can be considered as the effective degree of node i and an estimation is

$$\begin{array}{cc}\hfill \widehat{\mathrm{Var}}\left(W\right)& =\frac{{\sum }_{i}2(d_i+1)/(\alpha d_i+1)}{N{F}_l}\hfill \\ \hfill & =\frac{2\times \mathrm{average}\mathrm{effective}\mathrm{node}\mathrm{degree}}{F_l}\hfill \end{array}$$

(A13)

Appendix A.2.2. Backward View

Under linear assumption, we have ${\mathit{H}}_{\mathbf{l}}={\mathit{Y}}_{\mathbf{l}}$ and

$$\begin{array}{c}\hfill \frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}=\frac{\partial L}{\partial {\mathit{Y}}_{\mathbf{l}}}=\widehat{A}\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}+\mathbf{1}}}{W}_l^T,\frac{\partial L}{\partial W_{l-1}}={{\mathit{Y}}_{\mathbf{l}-\mathbf{1}}}^T\widehat{A}\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}},\end{array}$$

(A14)

We have

$$\begin{array}{cc}\hfill {\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)}_{ij}& =\sum _{t=1}^{F_{l+1}}\sum _{k=1}^N{\widehat{A}}_{ik}{\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}+\mathbf{1}}}\right)}_{kt}{\left(W_l^T\right)}_{t,j}\hfill \\ \hfill {\left(\frac{\partial L}{\partial W_{l-1}}\right)}_{ij}& ={\left(\widehat{A}{\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)}_{·i}^T{\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)}_{·j}=\sum _{k=1}^N\left(\sum _{t=1}^N{\widehat{A}}_{kt}{\left({\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)}_{ti}\right){\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)}_{kj},\hfill \end{array}$$

(A15)

And

$$\begin{array}{cc}\hfill \mathrm{Var}\left({\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)}_{ij}\right)& =\mathrm{Var}\left(\sum _{t=1}^{F_{l+1}}\sum _{k=1}^N{\widehat{A}}_{ik}{\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}+\mathbf{1}}}\right)}_{kt}{\left(W_l^T\right)}_{t,j}\right)\hfill \\ \hfill & =\frac{F_{l+1}}{d_i+1}\mathrm{Var}\left({\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}+\mathbf{1}}}\right)}_{kt}{\left(W_l^T\right)}_{t,j}\right)\hfill \\ \hfill & =\frac{F_{l+1}}{d_i+1}\mathrm{Var}\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}+\mathbf{1}}}\right)\mathrm{Var}\left(W_l\right)\hfill \\ \hfill & \approx \frac{N{F}_{l+1}}{{\sum }_i(d_i+1)}\mathrm{Var}\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}+\mathbf{1}}}\right)\mathrm{Var}\left(W_l\right)\hfill \\ \hfill & \approx \mathrm{Var}\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{n}+\mathbf{1}}}\right)\prod _{l^{\prime }=l}^n\frac{N{F}_{l+1}}{{\sum }_i(d_i+1)}\mathrm{Var}\left(W_l\right)\hfill \\ \hfill \mathrm{Var}\left({\left(\frac{\partial L}{\partial W_{l-1}}\right)}_{ij}\right)& =\mathrm{Var}\left(\sum _{k=1}^N\left(\sum _{t=1}^N{\widehat{A}}_{kt}{\left({\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)}_{ti}\right){\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)}_{kj}\right)\hfill \\ \hfill & =\sum _{k=1}^N\mathrm{Var}\left(\sum _{t=1}^N{\widehat{A}}_{kt}{\left({\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)}_{ti}\right)\mathrm{Var}\left({\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)}_{kj}\right)\hfill \\ \hfill & =\left(\sum _{k=1}^N\frac{1}{d_k+1}\right)\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)\mathrm{Var}\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{l}}}\right)\hfill \end{array}$$

(A16)

From (A4), $\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)$ can be approximately written as

$$\begin{array}{cc}\hfill & \mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}-\mathbf{1}}\right)=\mathrm{Var}\left({\mathit{H}}_{\mathbf{l}-\mathbf{1}}\right)\approx \frac{N{F}_{l-2}}{{\sum }_k{d}_k+1}\mathrm{Var}\left({\mathit{Y}}_{\mathbf{l}-\mathbf{2}}\left)\mathrm{Var}\right(W\right)\hfill \\ \hfill & \approx \mathrm{Var}\left({\mathit{Y}}_{\mathbf{0}}\right)\prod _{l^{\prime }=0}^{l-2}\frac{N{F}_{l^{\prime }}}{{\sum }_k{d}_k+1}\mathrm{Var}\left(W\right)\hfill \end{array}$$

(A17)

Then

$$\begin{array}{cc}\hfill & \mathrm{Var}\left({\left(\frac{\partial L}{\partial W_{l-1}}\right)}_{ij}\right)\hfill \\ & \approx \left(\sum _{k=1}^N\frac{1}{d_k+1}\right)\mathrm{Var}\left({\mathit{Y}}_{\mathbf{0}}\right)\mathrm{Var}\left(\frac{\partial L}{\partial {\mathit{H}}_{\mathbf{n}+\mathbf{1}}}\right)\prod _{l^{\prime }=0,l^{\prime }\ne l-1}^n\frac{N{F}_{l^{\prime }}}{{\sum }_k{d}_k+1}\mathrm{Var}\left(W\right)\hfill \end{array}$$

(A18)

Appendix A.3. Energy Analysis

Another way to design weight initialization is from the flow of energy perspective. Under the linear assumption and suppose we can do QR factorization of the weight matrices to make them orthogonal and $W^{T}W=W{W}^T=\alpha I$, then we have

$$\begin{array}{cc}\hfill & \mathrm{trace}\left({{\mathit{H}}_{\mathbf{l}+\mathbf{1}}}^T{\mathit{H}}_{\mathbf{l}+\mathbf{1}}\right)=\mathrm{trace}\left({\left(\widehat{A}{\mathit{Y}}_{\mathbf{l}}{W}_l\right)}^T\widehat{A}{\mathit{Y}}_{\mathbf{l}}{W}_l\right)\hfill \\ \hfill & =\mathrm{trace}\left(W_l^T{{\mathit{Y}}_{\mathbf{l}}}^T{\widehat{A}}^T\widehat{A}{\mathit{Y}}_{\mathbf{l}}{W}_l\right)=\mathrm{trace}\left(W_l{W}_l^T{{\mathit{Y}}_{\mathbf{l}}}^T{\widehat{A}}^T\widehat{A}{\mathit{Y}}_{\mathbf{l}}\right)\hfill \\ \hfill & =\mathrm{trace}\left(\alpha \sum _{i=1}^N{\lambda }_i^2{{\mathit{Y}}_{\mathbf{l}}}^T{\mathit{u}}_{\mathbf{i}}{\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{Y}}_{\mathbf{l}}\right)=\mathrm{trace}\left(\alpha \sum _{i=1}^N{\lambda }_i^2{\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{Y}}_{\mathbf{l}}{\left({\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{Y}}_{\mathbf{l}}\right)}^T\right)\hfill \\ \hfill & =\mathrm{trace}\left(\alpha \sum _{i=1}^N{\lambda }_i^2{∥{\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{Y}}_{\mathbf{l}}∥}_2^2\right)=\mathrm{trace}\left(\sum _{i=1}^N{∥{\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{H}}_{\mathbf{l}+\mathbf{1}}∥}_2^2\right)\hfill \end{array}$$

(A19)

Suppose all ${∥{\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{Y}}_{\mathbf{l}}∥}_2^2$ and ${\mathit{u}}_{\mathbf{i}}^{\mathbf{T}}{\mathit{H}}_{\mathbf{l}+\mathbf{1}}$ are equal, then

$$\begin{array}{c}\hfill \alpha =\frac{N}{{\displaystyle \sum _{i=1}^N}{\lambda }_i^2}=\frac{N}{{∥\widehat{A}∥}_F^2}=\frac{N}{{\displaystyle \sum _{i=1}^N}\frac{1}{d_i+1}}\end{array}$$

(A20)

If we use ReLU activation function, we have

$$\begin{array}{c}\hfill \alpha =\frac{2N}{{\displaystyle \sum _{i=1}^N}{\lambda }_i^2}=\frac{2N}{{∥\widehat{A}∥}_F^2}=\frac{2N}{{\displaystyle \sum _{i=1}^N}\frac{1}{d_i+1}}\end{array}$$

(A21)