\section{Introduction}

Graph neural network (GNN) \citep{scarselli2008graph,bruna2013spectral,defferrard2016convolutional,kipf2016semi,velivckovic2017graph,hamilton2017inductive,gilmer2017neural,xu2018powerful} has become the prominent approach to learn representations for graphs, such as social networks \citep{li2019encoding}, biomedical information networks \citep{yan2019groupinn}, communication networks \citep{suarez2021graph}, n-body systems \citep{kipf2018neural}, etc. Most GNNs rely on the message passing mechanism \citep{gilmer2017neural} to implement the interactions between neighbouring nodes. Despite its huge success, message passing GNNs still faces two fundamental but fatal drawbacks: it fails to generalize to heterophily where neighboring nodes share dissimilar features or labels, and a simple multilayer perceptron can outperform many GNNs \citep{zhu2020beyond}, this limit GNNs extending to many real-world networks with heterophily;
it is also observed the node representations become indistinguishable when stacking multiple layers, and suffers sharply performance drop, resulting in the so-called ``over-smoothing'' problem \citep{li2018deeper}, which prevent GNNs to utilize high-order neighborhood information effectively.


To address these two drawbacks, numerous approaches have been proposed. Most of them concentrate on the aggregation stage of message passing. Some design signed messages to distinguish neighbors belong to different classes \citep{yang2021diverse,bo2021beyond,luan2021heterophily,yan2021two}, allowing GNNs to capture high-frequency signals; \citet{min2020scattering} design specific filters to capture band-pass signals; some apply personalized aggregation with reinforcement learning \citep{lai2020policy} or neural architecture search \citep{wang2022graph}; others attempt to aggregate messages not only from the direct neighbors, but also from the embedding space \citep{pei2020geom} or higher-order neighbors \citep{zhu2020beyond}. These aggregator designs have achieved good performance, however, they primarily focus on the single-round message passing process and ignore the integration of messages from multiple hops. Another line of works focus on the effective utilization of multiple hops information, which is mainly accomplished by designing various skip-connections. \citet{klicpera2018predict,chen2020simple} propose initial connection to prevent ego or local information from being ``washed out'' by stacking multiple GNN layers; inspired by ResNet \citep{he2016deep}, some works \citep{li2019deepgcns,chen2020simple,cong2021provable} explores the application of residual connection on GNNs to improve the gradients; others combined the output of intermediate GNN layers with well-designed components, such as concat \citep{xu2018representation,zhu2020beyond}, learnable weights \citep{zhu2020beyond,abu2019mixhop,liu2020towards}, signed weights \citep{chien2020adaptive}, or RNN-like architectures \citep{xu2018representation,sun2019adagcn}. These works are simple yet effective, however, they can only model the information \textit{within few hops}, but unable to model the information exactly \textit{at some orders}, this lead to a mixing of features at different orders;
besides, many of these approaches \citep{li2019deepgcns,abu2019mixhop,chen2020simple,zhu2020beyond,chien2020adaptive,cong2021provable} are unable to make personalized decisions for each node. These deficiencies result in suboptimal performance. In addition to caring about the model side, other approaches focus on how to modify the graph structure. These methods are called ``graph rewiring'', including randomly removing edges \citep{rong2019dropedge} or nodes \citep{feng2020graph}, or computing a new graph with heuristic algorithms \citep{suresh2021breaking,zeng2021decoupling}. In general, these algorithms are not learnable and thus only applicable to certain graphs.


\begin{wrapfigure}{r}{0.41\textwidth}
\vspace{-0.4cm}
\centering
\includegraphics[width=1.0\linewidth]{nested_tree.pdf}
\caption{Aligning the hierarchy of a rooted-tree $\mathcal{T}_{v}^{(k)}$ underlying the graph with the node embedding of the root node $v$. Neighboring nodes within $k$ hops of edges to $v$ naturally form a depth $k$ subtree. Messages passed to $v$ from nodes within this subtree are restricted to the first $P_v^{(k)}$ neurons in the node embedding of $v$.}
\label{fig.nested_tree}
\vspace{-0.45cm}
\end{wrapfigure}

Unlike the previous works, we address both problems by designing the combine stage of message passing and emphasize the importance of it.
\textit{The key idea is to integrate an inductive bias from rooted-tree hierarchy, let GNN encode the neighborhood information exactly at some orders and avoid feature mixing within hops.}
The combine stage has been rarely focused before, most works simply implement it as a self-loop. This would result in an unreasonable mixing of node features \citep{zhu2020beyond}.
To avoid this ``mixing'', \citet{hamilton2017inductive,xu2018representation,zhu2020beyond} concat the node representation and the aggregated message, which has been identified as an effective design to deal with heterophily \citep{zhu2020beyond}. However, keeping the embedding dimension constant across layers, the local information will be squeezed at an exponential rate. The most related work to ours is Gated GNN \citep{li2015gated}, it applys a GRU in the combine stage and strengthen the expressivity, but fails to prevent feature mixing, limiting the performance.

In this paper, we present the message passing mechanism in an ordered form. That is, the neurons in the node embedding of a certain node is aligned with the hierarchies in the rooted-tree of the node. Here, by rooted-tree of a node we refer to the tree with the node itself as the root and its neighbors as its children. Recursively, for each child, its children are again the neighboring nodes of the child. (c.f. Figure \ref{fig.nested_tree}). We achieve the alignment by proposing a novel ordered gating mechanism, which controls the assignment of neurons to encode subtrees with different depth. Experimental results on an extensive set of datasets show that our model could alleviate the heterophily and over-smoothing problem at the same time. Our model provides the following practical benefits:

\begin{itemize}
\item We design the combine stage guided by the rooted-tree hierarchy, a very general topological inductive bias with least assumption about the neighbors' distribution, this allows a flexible integration of information at different orders, and leading superior performance in both heterophily and homophily.

\item The ordered gating mechanism prevent the mixing of node features within hops, enable us to model information at different orders. This open a door to extract similar neighborhood patterns for each node under heterophily; this also make it easy to preserve ego and local information, then effectively alleviating over-smoothing.

\item Our model aligns neighboring structures with blocks in node embeddings through an explicit gating mechanism, thus the gating mechanism could provide visualizations to reveal the connecting type of the data and offer explainability.
\end{itemize}


\section{Our Approach}
In this section, we will first present a general formalization of GNNs. Then we show how to align the hierarchy of the rooted-tree to the node embeddings, the intuitions behind it, as well as the two necessary components to imeplement it. Finally, our model is presented in Section \ref{puttogether}, with few cases to explain it's practical benefits.

\subsection{Graph Neural Networks}
Consider an unweighted and self-loop-free graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, with $\mathcal{V}$ and $\mathcal{E}$ denoting its nodes and edges, respectively. Let $v \in \mathcal{V}$ being one of the nodes with its label $y_{v} \in \mathcal{Y}$, and assume that there are $N$ nodes in the graph. The adjacency matrix for graph $\mathcal{G}$ is denoted as $A \in \mathbb{R}^{N \times N}$.
Assume that the input features for each node form a matrix $X \in \mathbb{R}^{N \times F}$, where $F$ is the dimension of input features, with $X_v \in \mathbb{R}^N$ being the features of node $v$. Before feeding it into the GNN, an input projection function $f_{\theta}(\cdot)$ is optionally incorporated to nonlinearly project the features into a space mathing the node embedding dimensions, i.e., 
\begin{equation}
Z_v=f_{\theta}(X_v)
\end{equation}
where $D$ is the dimension of node embeddings, and $Z_v \in \mathbb{R}^{D}$ is feed to the first layer of GNN, i.e., $h_{v}^{(0)}=Z_{v}$. In our work, we simply parameterize $f_{\theta}(\cdot)$ with an MLP. 


Generally, a GNN layer with message passing can be divided into two stages: the \textit{aggregate} stage and the \textit{combine} stage. For the node $v$ at the $k$-th layer, it first aggregates information from its neighbors, forming the message vector $m_{v}^{(k)}$. Then in the second stage $m_{v}^{(k)}$ is combined with $v$'s ego representation $h_{v}^{(k-1)}$ from the last layer. Formally,
\begin{equation}
m_{v}^{(k)}=\operatorname{AGGREGATE}^{(k)}\left(\left\{h_{u}^{(k-1)}: u \in \mathcal{N}(v)\right\}\right), h_{v}^{(k)}=\operatorname{COMBINE}^{(k)}\left(h_{v}^{(k-1)}, m_{v}^{(k)}\right)
\end{equation}
where $\mathcal{N}(v)$ is the set of node $v$'s direct neighbors. The $\operatorname{AGGREGATE}(\cdot)$ operator denotes a permutation-invariant function on all the neighboring nodes of $v$, such as addition, mean pooling or max pooling, optionally followed by a transformation function. The transformation function could be an MLP \citep{xu2018powerful} or transformation matrix \citep{kipf2016semi}. The $\operatorname{COMBINE}(\cdot)$ operator denotes a combine function that fuses the ego information and the aggrated message. 

Our model only modifies the \textit{combine} stage, which means it could be seamlessly extended to many popular GNNs, such as GCN \citep{kipf2016semi}, GAT \citep{velivckovic2017graph}, GIN \citep{xu2018powerful}, etc., by changing the \textit{aggregate} stage correspondingly. In this paper, we simply implement the \textit{aggregate} stage as computing the mean value of neighboring vectors.


\subsection{Aligning Rooted-Tree with Node Embedding}
\label{align_rooted_tree}
From the viewpoint of a node, the neighboring nodes within multiple hops around it can be naturally presented as a rooted-tree \citep{xu2018powerful,xu2018representation,liu2020towards}. For node $v$, its $k$-th order rooted-tree $\mathcal{T}_{v}^{(k)}$ is composed by putting its first order neighbors as its children, and then recursively for each child, putting the neighboring nodes of the child as its children (except the child that is also the ancester of the node), unitl its $k$-th order neighbors (Figure \ref{fig.nested_tree}). 

Since there is only one round of message passing in each GNN layer, stacking multiple layers of GNNs corresponds to doing multiple rounds of message passing. And since for each single round of message passing, only direct neighbors exchange their information, the node representation of $v$ in the $k$-th layer $h_v^{(k)}$ is only depending on the nodes within its $k$-th order rooted-tree $\mathcal{T}_{v}^{(k)}$. Our goal is to properly encode the information within the $k$-th order rooted-tree into the fixed size vector $h_v^{(k)}$. 

Note that there is a nested structure between rooted-trees of $v$, and we are going to utilize this structure as anchors to align the rooted-trees with the node embedding of $v$. It is obvious that the $k-1$-th rooted-tree $\mathcal{T}_{v}^{(k-1)}$ is a subtree of the $k$-th rooted-tree  $\mathcal{T}_{v}^{(k)}$, with the same root $v$ and the same tree structure up to depth $k-1$. Thus we have (assuming the model has $K$ layers)
\begin{equation} \label{relativemagnitude}
    \mathcal{T}_{v}^{(0)} \subseteq \mathcal{T}_{v}^{(1)} \subseteq \cdots \subseteq \mathcal{T}_{v}^{(k)} \subseteq \cdots \subseteq \mathcal{T}_{v}^{(K)}.
\end{equation}
where $\mathcal{T}_v^{(0)}$ corresponds to the neurons used to represent the node's ego information.
As $k$ grows, $\mathcal{T}_{v}^{(k)}$ becomes larger and more complex, requiring more neurons to encode its information. As a result, it is natural to assume that the neurons used to represent message passing within $\mathcal{T}_{v}^{(k-1)}$ should be a subset of those for $\mathcal{T}_{v}^{(k)}$. Correspondingly, we are going to implement the same nested structure in the node embedding of $v$. We first order all the neurons in $h_v$ into a sequence, thus we can index them with numbers. For information within $\mathcal{T}_{v}^{(k-1)}$, we allow it to be encoded within the first $P_v^{(k-1)}$ neurons. As for the next level rooted-tree $\mathcal{T}_{v}^{(k)}$, more neurons are involved, corresponding to the first $P_v^{(k)}$ neurons. Apparently here we have $P_v^{(k-1)} \leq P_v^{(k)}$.

There is a critical property that comes with this way of ordering the neurons. Consider the neurons between the two split points $P_v^{(k-1)}$ and $P_v^{(k)}$. They reflect the delta between $k-1$ and $k$ rounds of message passing, encoding the neighbor information exactly \textit{at $k$-th order}. Intuitively, it organizes the message passing in an \textit{ordered} form, avoid the mixing of neighbor features at different orders within a hop. We implement the ordered neurons with the gating mechanism (which will be introduced in the following sections in detail), thus it's natural to visualize the gates as a way to probe the message passing of the model.

Now we have the node embedding of $v$ split by a set of spliting points 
\begin{equation} \label{relativemagnitude2}
    P_v^{(0)} \leq P_v^{(1)} \leq \cdots \leq P_v^{(k)} \leq \cdots \leq P_v^{(K)}.
\end{equation}
Each of them is predicted by a GNN layer, with $P_v^{(K)}=D$ pointing to the end of the node embedding. If we can learn the positions of these spliting points and ensure their relative magnitudes, then the alignment could be realized. We are going to show how we predict the spliting points through a gating mechanism, and how we lock their relative magnitudes with the differentiable OR operator. Our model naturally emerges by putting everything together in Section \ref{puttogether}.


\subsection{The Split Points}
Consider the split point $P_{v}^{(k)}$ that split up the ordered node embedding into two blocks, with the left block containing neurons indexed between $[0, P_{v}^{(k)}-1]$ and the right block $[P_{v}^{(k)}, D]$. We can represent the split by a $D$-dimensional gating vector $g_v^{(k)}$, with its left $P_{v}^{(k)}$ entries are $1$s and the rest are $0$s. This gating vector could then be used in the \textit{combine} stage to control the integration between the ego representation as input to the $k$-th layer $h_{v}^{(k-1)}$ and the aggregated context $m_{v}^{(k)}$.
\begin{equation}
h_{v}^{(k)}=g_{v}^{(k)} \circ h_{v}^{(k-1)} + \left(1-g_{v}^{(k)}\right) \circ m_{v}^{(k)}
\end{equation}
Where $\circ$ denotes element-wise multiplication. Intuitively, for the left $P_{v}^{(k)}$ neurons in $h_{v}^{(k)}$, since their corresponding gating vectors in $g_{v}^{(k)}$ are $1$s, they suppress the newly aggregated context $m_{v}^{(k)}$ and only inherit the last layer's output $h_{v}^{(k-1)}$. And vice versa for the right block of neurons with gating vectors of $0$.

Ideally, we would prefer a crispy split, with a clear boundary at $P_{v}^{(k)}$ and binary gating vectors in $g_{v}^{(k)}$. However, this will result in discretized operations, rendering the model non-differentiable. In this work, we resort to ``soften'' the gates by predicting the expectations of them, thus keeping the whole model differentiable. The expectation vector $\hat{g}_{v}^{(k)}$ is presented as the cumulative sum of the probability of each position in the node embedding being the split point $P_{v}^{(k)}$.\footnote{We refer the readers to Appendix \ref{proof} for a proof of $\hat{g}_{v}^{(k)}$ being the expectation of the gating vector $g_{v}^{(k)}$.} Concretely, $\hat{g}_{v}^{(k)}$ is parameterized as
\begin{equation}
\label{eq.ordered_gating}
\hat{g}_{v}^{(k)} = \operatorname{cumax}_{\leftarrow}\left( f_{\xi}^{(k)} \left( h_{v}^{(k-1)}, m_{v}^{(k)} \right) \right) = \operatorname{cumax}_{\leftarrow}\left( W^{(k)} \left[ h_{v}^{(k-1)}; m_{v}^{(k)} \right]+b^{(k)} \right)
\end{equation}
where $\operatorname{cumax}_{\leftarrow}(\cdot)=\operatorname{cumsum}_{\leftarrow}(\operatorname{softmax}(\ldots))$ is the composite operator of cumulative sum and softmax, with the subscript $\leftarrow$ indicating that the direction of cumulation is from right to left. $f_{\xi}^{(k)}(\cdot)$ denotes a function fusing the two vectors $h_{v}^{(k-1)}$ and $m_{v}^{(k)}$. Here we simply put it as a linear projection with bias $b^{(k)}$. The $\left[ h_{v}^{(k-1)}; m_{v}^{(k)} \right]$ indicates the concatenation of the two vectors $h_{v}^{(k-1)}$ and $m_{v}^{(k)}$.


\subsection{The Differentiable $\operatorname{OR}$ Operator}
In the above section we provide a differentiable way to predict the position of the split points, as well as corresponding gates that implement the split operation. However, there's no guarantee that the predicted split points will comply with the restrictions in Equation \ref{relativemagnitude2}. Not ensuring Equation \ref{relativemagnitude2} would break the alignment between the rooted-trees and node embeddings. 

To ensure such a relative magnitude, we conduct a bitwise $\operatorname{OR}$ operator between the newly computed gating vector $\hat{g}_{v}^{(k)}$ and its predecessor. We change the notation to $\tilde{g}_{v}^{(k)}$ to discriminate it from the originally computed gating vector $\hat{g}_{v}^{(k)}$. Again, the $\operatorname{OR}$ operator is not differentiable in nature. We instead implement the softened version of it by 
\begin{equation}
\label{softor}
\tilde{g}_{v}^{(k)} = \operatorname{SOFTOR}(\tilde{g}_{v}^{(k-1)}, \hat{g}_{v}^{(k)}) = \tilde{g}_{v}^{(k-1)} + (1-\tilde{g}_{v}^{(k-1)}) \circ \hat{g}_{v}^{(k)}
\end{equation}
We call it the $\operatorname{SOFTOR}(\cdot)$ operator.


\begin{figure
\centering
\includegraphics[width=\linewidth]{rooted_tree_cases.pdf}
\caption{Some rooted-tree cases, each rooted-tree is constructed from root node by message passing, node's color indicate it's label. Case 1 and Case 2 show two cases under heterophily; Case 3 is under homophily, node 1 \& 2 are at the center and the border of a community, respectively.}
\label{fig.showing_cases}
\end{figure}


\subsection{Putting it All Together}
\label{puttogether}
With the above discussion, we can get the updating rule for the $k$-th GNN layer by putting everything together, leading to the \textit{Ordered GNN} model:

\begin{equation}
\begin{aligned}
m_{v}^{(k)}&=\operatorname{MEAN}\left(\left\{h_{u}^{(k-1)}: u \in \mathcal{N}(v)\right\}\right) \\
\hat{g}_{v}^{(k)}&=\operatorname{cumax}_{\leftarrow}\left( f_{\xi}^{(k)} \left( h_{v}^{(k-1)}, m_{v}^{(k)} \right) \right) \\
\tilde{g}_{v}^{(k)}&=\operatorname{SOFTOR}\left(\tilde{g}_{v}^{(k-1)}, \hat{g}_{v}^{(k)}\right)\\
h_{v}^{(k)}&=\tilde{g}_{v}^{(k)} \circ h_{v}^{(k-1)} + \left(1-\tilde{g}_{v}^{(k)}\right) \circ m_{v}^{(k)}
\end{aligned}
\end{equation}

The first equation corresponds to the \textit{aggregate} stage, while the latter three correspond to the \textit{combine} stage. In practice, we can decrease the parameters in the gating network (parameters in the second equation) by applying a chunking trick, i.e., we set $D_{m}=D/C$ as the size of the gating vector, where $D$ is the dimension of the GNN's hidden state, and $C$ is the chunking size factor. In this setting, each gate controls $C$ neurons instead, thus drastically reducing the number of gates to be predicted.


In order to explain the practical benefits our model provided for heterophily and over-smoothing problems, we illustrate three cases in the Figure \ref{fig.showing_cases} and analyse them with the insights from related theoretical works. \citet{luan2021heterophily,ma2021homophily} notice that heterophily is not always harmful for GNNs, and \citet{ma2021homophily} points out that a simple GCN can work well if similar neighbor patterns are shared across nodes with the same label, they called it ``good heterophily''.
In Case 1, two root nodes share the same first-order neighborhood pattern, one GNN layer can capture this ``good heterophily'' and obtain good performance, however, a multi-layer GNN will mix in the second-order neighbor features that have different pattern, the captured ``good heterophily'' would be destroyed. Our model can protect the captured ``good heterophily'' and reject that ``bad heteropihly'' by learning same split points for layer 1 \& 2.
In Case 2, two root nodes' first-order neighbor patterns differ, but they share the same second-order neighborhood pattern. Since an ordinary GNN can only capture information \textit{within $k$-hops}, they can not capture this ``good heterophily'' \textit{at second-order}, whereas our model could do so with the help of split points prediction. Furthermore, we can filter out the first-order ``bad heterophily''.

Some works \citep{li2018deeper,chen2020measuring,ma2021homophily} show that neighbor feature smoothing is the key reason why GNN works, but too much smoothing lead to the mixing of nodes from different classes. Case 3 shows that while a GNN allows node 1 to smooth neighbor features within it's class, it also mixing node features belong to different classes for node 2. With the fine-grained ordered gating mechanism, our model can allow node 1 to aggregate similar features while filtering out information beyond node 2's own, thus helpful for homophily and over-smoothing.

The analysis of these cases shows that it is critical to predict a set of \textit{flexible} split points. By analyzing the convergence rate of the gates, we prove that our model only constrain the relative position of the split points, the delta between split points is very flexible, i.e. the number of neurons assigned to encode the information at different orders is flexible. We provide the proof in Appendix \ref{proof_conv_speed}.




\section{Experiments}
\label{part.Experiments}
In this section, we test our model on 11 datasets. We mainly evaluate our model in 3 aspects: performances on homophily data, on heterophily data, and the robustness towards over-smoothing. We further use 2 larger datasets to see if our model could scale well w.r.t. the size of the dataset, these two datasets are from large-scale homophily and heterophily benchmarks, respectively.


\begin{wraptable}{r}{0.6\textwidth}
\vspace{-0.5cm}
\caption{Dataset statistics.}
\label{table.dataset_stats}
\centering
\vspace{0.1cm}
\resizebox{0.6\textwidth}{!}{
\begin{tabular}{c|ccccc}
\toprule
Dataset &  Classes &  Nodes &   Edges &  Features & $\mathcal{H}(G)$ \\
\midrule
Cora &        7 &   2,708 &    5,429 &      1,433 & 0.81\\
PubMed &        3 &  19,717 &   44,338 &       500 & 0.80\\
CiteSeer &        6 &   3,327 &    4,732 &      3,703 & 0.74\\
Cornell &        5 &    183 &     280 &      1,703 & 0.30\\
Chameleon &        5 &   2,277 &   31,421 &       2,325 & 0.23\\
Squirrel &        5 &   5,201 &   198,493 &       2,089 & 0.22\\
Actor &        5 &   7,600 &   26,752 &       931 & 0.22\\
Texas &        5 &    183 &     295 &      1,703 & 0.21\\
Wisconsin &        5 &    251 &     466 &      1,703 & 0.11\\
ogbn-arxiv &       40 & 169,343 & 1,166,243 &       128 & 0.66\\
arXiv-year &        5 & 169,343 & 1,166,243 &       128 & 0.22\\
\bottomrule
\end{tabular}
}
\vspace{-0.2cm}
\end{wraptable}

\paragraph{Datasets \& Experiment Setup}
\label{exp_setup}
For homophily data, we use three citation network datasets, i.e., Cora, CiteSeer, and PubMed \citep{sen2008collective}. While six heterophilous web network datasets \citep{pei2020geom,musae}, i.e., Actor, Texas, Cornell, Wisconsin, Squirrel and Chameleon, are used to evaluate our model on the heterophily setting. In citation networks, nodes and edges correspond to documents and citations, respectively. Node features are bag-of-words representation of the documents. In web networks, nodes and edges represent web pages and hyper-links, with node features being bag-of-words representation of web pages. For the over-smoothing problem, following \citet{rong2019dropedge}, we use Cora, CiteSeer, and PubMed as testbeds. For the two larger datasets, we select ogbn-arxiv from \citet{hu2020open} and arXiv-year from \citet{lim2021new}. These two dataset have same nodes and edges that represent papers and citations. In ogbn-arxiv, node labels represent the papers' academic topics, which tend to be a homophily setting, while in arXiv-year, the node labels represent the publishing year, thus exhibiting more heterophily label patterns. We include the edge homophily score $\mathcal{H}(G)$ \citep{yan2021two} for each dataset, which represents the homophily level, together with the statistics of the datasets summarized in Table \ref{table.dataset_stats}.
We use the Adam optimizer \citep{kingma2014adam}, and apply dropout and $L_2$ regularization for $f_{\theta}$ and $f_{\xi}^{(k)}$. We also insert a LayerNorm \citep{ba2016layer} layer after every other layers. We set the hidden dimension as 256 and chunk size factor as 4. For each dataset, unless otherwise noted, a 8 layer model is used (which is already deeper than most popular GNN models).
We perform a grid search to tune the hyper-parameters for all the models.
More details are listed in Appendix \ref{exp_detail}.







\subsection{Homophily and Heterophily}
\label{setting_homo_hete}
We evaluate our model's performance on aforementioned 9 datasets, following the same data splits as \citet{pei2020geom}. Results are shown in Table \ref{table.results.hete_homo}. We report the mean classification accuracy with the standard deviation on the test nodes over 10 random data splits. For baselines, we include classic GNN models such as vanilla GCN \citep{kipf2016semi}, GAT\citep{velivckovic2017graph}, and GraphSAGE\citep{hamilton2017inductive}; state-of-the-art methods that are specifically designed for heterophily: H2GCN \citep{zhu2020beyond}, GPRGNN \citep{chien2020adaptive}, and GGCN \citep{yan2021two}; state-of-the-art models on three homophily citation network datasets: GCNII \citep{chen2020simple}, Geom-GCN \citep{pei2020geom}. We also include models that utilize the information from different hops: MixHop \citep{abu2019mixhop}, JK-Net \citep{xu2018representation}; and finally a 2-layer MLP. For the baselines that have multiple architectural variants, we choose the best one for each dataset and denote the model with an asteriod ``*''. All the baseline results are from \citet{yan2021two}, except for those on H2GCN, JK-Net and MixHop are from \citet{zhu2020beyond}. All of them were evaluated on the same set of data splits as ours.

\begin{table
\caption{Experimental results on homophily and heterophily datasets, for each dataset, we bold the model with the best performance.}
\label{table.results.hete_homo}
\centering
\resizebox{\textwidth}{!}
{
\begin{tabular}{c|cccccc|ccc}
\toprule
& Texas & Wisconsin & Actor & Squirrel & Chameleon & Cornell & CiteSeer & PubMed & Cora \\
\midrule
\textbf{Ordered GNN} & $\textbf{86.22}{\scriptstyle\pm4.12}$ & $\textbf{88.04}{\scriptstyle\pm3.63}$ & $\textbf{37.99}{\scriptstyle\pm1.00}$ & $\textbf{62.44}{\scriptstyle\pm1.96}$ & $\textbf{72.28}{\scriptstyle\pm2.29}$ & $\textbf{87.03}{\scriptstyle\pm4.73}$ & $77.31{\scriptstyle\pm1.73}$ & $\textbf{90.15}{\scriptstyle\pm0.38}$ & $\textbf{88.37}{\scriptstyle\pm0.75}$ \\
\midrule
H2GCN* & $84.86{\scriptstyle\pm7.23}$ & $87.65{\scriptstyle\pm4.98}$ & $35.70{\scriptstyle\pm1.00}$ & $36.48{\scriptstyle\pm1.86}$ & $60.11{\scriptstyle\pm2.15}$ & $82.70{\scriptstyle\pm5.28}$ &  $77.11{\scriptstyle\pm1.57}$ & $89.49{\scriptstyle\pm0.38}$ & $87.87{\scriptstyle\pm1.20}$ \\
GPRGNN & $78.38{\scriptstyle\pm4.36}$ & $82.94{\scriptstyle\pm4.21}$ & $34.63{\scriptstyle\pm1.22}$ & $31.61{\scriptstyle\pm1.24}$ & $46.58{\scriptstyle\pm1.71}$ & $80.27{\scriptstyle\pm8.11}$ &  $77.13{\scriptstyle\pm1.67}$ & $87.54{\scriptstyle\pm0.38}$ & $87.95{\scriptstyle\pm1.18}$ \\
GGCN & $84.86{\scriptstyle\pm4.55}$ & $86.86{\scriptstyle\pm3.29}$ &
$37.54{\scriptstyle\pm1.56}$ & $55.17{\scriptstyle\pm1.58}$ & $71.14{\scriptstyle\pm1.84}$ & $85.68{\scriptstyle\pm6.63}$ &  $77.14{\scriptstyle\pm1.45}$ & $89.15{\scriptstyle\pm0.37}$ & $87.95{\scriptstyle\pm1.05}$ \\
\midrule
GCNII* & $77.57{\scriptstyle\pm3.83}$ & $80.39{\scriptstyle\pm3.4}$ & $37.44{\scriptstyle\pm1.30}$ & $38.47{\scriptstyle\pm1.58}$ & $63.86{\scriptstyle\pm3.04}$ & $77.86{\scriptstyle\pm3.79}$ & $77.33{\scriptstyle\pm1.48}$ &
$\textbf{90.15}{\scriptstyle\pm0.43}$ &
$\textbf{88.37}{\scriptstyle\pm1.25}$ \\
Geom-GCN* & $66.76{\scriptstyle\pm2.72}$ & $64.51{\scriptstyle\pm3.66}$ & $31.59{\scriptstyle\pm1.15}$ & $38.15{\scriptstyle\pm0.92}$ & $60.00{\scriptstyle\pm2.81}$ & $60.54{\scriptstyle\pm3.67}$ &  $\textbf{78.02}{\scriptstyle\pm1.15}$ & $89.95{\scriptstyle\pm0.47}$ & $85.35{\scriptstyle\pm1.57}$ \\
\midrule
MixHop & $77.84{\scriptstyle\pm7.73}$ & $75.88{\scriptstyle\pm4.90}$ & $32.22{\scriptstyle\pm2.34}$ & $43.80{\scriptstyle\pm1.48}$ & $60.50{\scriptstyle\pm2.53}$ & $73.51{\scriptstyle\pm6.34}$ & $76.26{\scriptstyle\pm1.33}$ & $85.31{\scriptstyle\pm0.61}$ & $87.61{\scriptstyle\pm0.85}$ \\ 
JK-Net* & $83.78{\scriptstyle\pm2.21}$ & $82.55{\scriptstyle\pm4.57}$ & $35.14{\scriptstyle\pm1.37}$ & $45.03{\scriptstyle\pm1.73}$ & $63.79{\scriptstyle\pm2.27}$ & $75.68{\scriptstyle\pm4.03}$ & $76.05{\scriptstyle\pm1.37}$ & $88.41{\scriptstyle\pm0.45}$ & $85.96{\scriptstyle\pm0.83}$ \\
\midrule
GCN & $55.14{\scriptstyle\pm5.16}$ & $51.76{\scriptstyle\pm3.06}$ &  $27.32{\scriptstyle\pm1.10}$ & $53.43{\scriptstyle\pm2.01}$ & $64.82{\scriptstyle\pm2.24}$ &
$60.54{\scriptstyle\pm5.30}$ &  $76.50{\scriptstyle\pm1.36}$ & $88.42{\scriptstyle\pm0.50}$ & $86.98{\scriptstyle\pm1.27}$ \\
GAT & $52.16{\scriptstyle\pm6.63}$ & $49.41{\scriptstyle\pm4.09}$ & $27.44{\scriptstyle\pm0.89}$ & $40.72{\scriptstyle\pm1.55}$ & $60.26{\scriptstyle\pm2.5}$ & $61.89{\scriptstyle\pm5.05}$ &  $76.55{\scriptstyle\pm1.23}$ & $86.33{\scriptstyle\pm0.48}$ & $87.30{\scriptstyle\pm1.10}$ \\
GraphSAGE & $82.43{\scriptstyle\pm6.14}$ & $81.18{\scriptstyle\pm5.56}$ & $34.23{\scriptstyle\pm0.99}$ & $41.61{\scriptstyle\pm0.74}$ & $58.73{\scriptstyle\pm1.68}$ & $75.95{\scriptstyle\pm5.01}$ &  $76.04{\scriptstyle\pm1.30}$ & $88.45{\scriptstyle\pm0.50}$ & $86.90{\scriptstyle\pm1.04}$ \\
\midrule
MLP & $80.81{\scriptstyle\pm4.75}$ & $85.29{\scriptstyle\pm3.31}$ & $36.53{\scriptstyle\pm0.70}$ & $28.77{\scriptstyle\pm1.56}$ & $46.21{\scriptstyle\pm2.99}$ & $81.89{\scriptstyle\pm6.40}$ &  $74.02{\scriptstyle\pm1.90}$ & $87.16{\scriptstyle\pm0.37}$ & $75.69{\scriptstyle\pm2.00}$ \\
\bottomrule
\end{tabular}
}
\end{table}

Our Ordered GNN model consistently achieves SOTA both on homophilous and heterophilous datasets. On homophilous datasets, we match the performance of the SOTA model GCNII. However, GCNII performs worse on heterophilous datasets, leaving a gap of more than 20 percents on Squirrel dataset. The amazing performance on this dataset might related to the ``over-squashing'' problem \citep{alon2020bottleneck,topping2021understanding} (see Appendix \ref{Impact}). We also achieve SOTA on the heterophilous datasets and consistently outperform baselines across all datasets. The superior performances on both homophily and heterophily data show the generality of our model.


\begin{figure
\centering
\includegraphics[width=0.9\linewidth]{results_over_smoothing}
\caption{Node classification accuracy results for over-smoothing problem with various depths.}
\label{fig.results_over_smoothing}
\end{figure}

\subsection{Over-smoothing}
\label{setting_smooth}
We test our model's performance with depth 2, 4, 8, 16, 32 to evaluate the ability in alleviating the over-smoothing problem, we include following models for comparison: classical GNN models GCN \citep{kipf2016semi}, GAT \citep{velivckovic2017graph} and GraphSAGE \citep{hamilton2017inductive}, models that modify the architecture to alleviate over-smoothing: ResGCN \citep{li2019deepgcns}, JK-Net \citep{xu2018representation}, and IncepGCN \citep{rong2019dropedge}. We also include three state-of-the-art methods ResGCN+Drop, JK-Net+Drop, and IncepGCN+Drop that is equipped with the advanced training trick DropEdge \citep{rong2019dropedge} which is developed for alleviating over-smoothing. We reuse the results from \citet{rong2019dropedge} for all baselines except GAT. For GAT, we conduct a separate grid search for each depth to tune it's hyperparameters, with the search space the same to ours for fair comparison. Detailed configurations can be found in Appendix \ref{exp_detail}. We follow the same data split from \citet{rong2019dropedge} as used in baselines.

Performances are shown in Figure \ref{fig.results_over_smoothing}. Compared with the baseline models, our model consistently achieve the best performance across all datasets and all layers. Furthermore, for each dataset, our model's performance doesn't drop significantly as some baselines. Note that this is achieved without using training tricks such as DropEdge \citep{rong2019dropedge}. This shows that our model can well alleviate the over-smoothing problem.



\begin{figure}
\centering
\subfigure[Homophily]{
\label{fig.gate_values.Homophily}
\begin{minipage}[b]{0.26\linewidth}
\label{fig.gate_values.CiteSeer_geom_8}
\includegraphics[width=1\linewidth]{CiteSeer_geom_8.png}\vspace{4pt}
\label{fig.gate_values.PubMed_8}
\includegraphics[width=1\linewidth]{PubMed_geom_8.png}\vspace{4pt}
\end{minipage}
}
\subfigure[Heterophily]{
\label{fig.gate_values.Heterophily}
\begin{minipage}[b]{0.26\linewidth}
\label{fig.gate_values.Actor_8}
\includegraphics[width=1\linewidth]{Actor_8.png}\vspace{4pt}
\label{fig.gate_values.Texas_8}
\includegraphics[width=1\linewidth]{Texas_8.png}\vspace{4pt}
\end{minipage}
}
\subfigure[Over-smoothing]{
\label{fig.gate_values.Over-smoothing}
  \begin{minipage}[b]{0.26\linewidth}
\label{fig.gate_values.Cora_32}
\includegraphics[width=1\linewidth]{Cora_full_32.png}\vspace{4pt}
\end{minipage}
}
\caption{Visualization of gating vectors across all layers in our models on (a) homophily, (b) heterophily and (c) over-smoothing settings. Vertical and horizontal axes are gating vectors and chunk channels, respectively. Each box in the boxplot refect the distribution of the value of the corresponding gate across the whole dataset.}
\label{fig.gate_values}
\end{figure}

\subsection{Visualizing the Gating Vectors}
We checkout the gating vectors across all layers in our model, plotted in Figure \ref{fig.gate_values}. We can find that the gates automatically learns to behave differently for homophilous and heterophilous data, and learns to suppress message passing when needed to prevent over-smoothing. 

On homophilous data the split points distributes more uniformly at first layer, with a gradually increasing trend of gating vectors. As the layer gets deeper, gates start to saturate. After layer 4, most of the positions are saturated, leaving the last very few positions open for update. This coinsides with the fact that 4 layers are the depth when many GNNs starts to show the over-smoothing phenomenon.

On heterophilous data, a large block of gates become saturated even in the first layer, emphasizes the importance of a node's ego information by suppressing most of the incoming messages. This observation is consistent with the amazing performance of a simple MLP model when deal with heterophily (see Table \ref{table.results.hete_homo}). Moreover, on Texas dataset the saturated gates keep the same for layer 1 and 2, however enlarge on layer 3, reflect the flexibility of our gating mechanism and the ability to capture ``good heterophily'' for high-order neighbors as we described in Section \ref{puttogether}.

We also plot the gating vectors of our model on Cora dataset in Figure \ref{fig.gate_values.Over-smoothing}. We can see that almost all the gates are closed when layers become very deep. This phenomenon reveals Ordered GNN's superior performance on deep models: it manages to keep node representations distinguishable by learning to reject incoming messages from farther-away neighbors.





\subsection{Scaling to Larger Datasets}
\label{setting_large}
The larger ogbn-arxiv and arXiv-year datasets contain much more nodes, and present a complex linking pattern. We compared our model with these models: baseline graph neural models SGC \citep{wu2019simplifying}, GCN \citep{kipf2016semi}, label propagation based methods LP \citep{zhou2003learning}, C\&S \citep{huang2020combining}, strong baselines: JK-Net \citep{xu2018representation}, APPNP \citep{klicpera2018predict}, and state-of-the art model UniMP \citep{shi2020masked}, 
GCNII \citep{chen2020simple} for ogbn-arxiv, and LINK \citep{lim2021new}, H2GCN \citep{zhu2020beyond} for arXiv-year. Features, links and labels are inherited from \citet{hu2020open}, all the models follow the same data split from \citet{hu2020open} and \citet{lim2021new}. The baselines results are from \citet{lim2021new}, \citet{li2020deepergcn} and \citet{huang2020combining}.

For homophily dataset ogbn-arxiv, our model beats all the baseline methods except UniMP. This is understandable since UniMP based on the graph Transformer much larger model size, and uses the extra node label information. For heterophily dataset, our method outperform the SOTA models. The performance on these two large dataset demonstrate the scalability and expressivity of our model to capture the complex patterns in large networks.


\begin{table}
\caption{Summary of classification accuracy results on larger benchmarks, we bold the model with the best performance.}
\label{table.results.larger_benchmark}
\centering
\resizebox{\textwidth}{!}
{
\begin{tabular}{l|ll|ll|ll|ll|l}
\toprule
& GCN & SGC & LP & C\&S & JK-Net & APPNP & GCNII & UniMP & \textbf{Ordered GNN} \\
ogbn-arxiv & $71.74{\scriptstyle\pm0.29}$ & $69.39$ & $68.50$ & $72.62$ & $72.19{\scriptstyle\pm0.21}$ & $66.38$ & $72.74{\scriptstyle\pm0.16}$ & $\textbf{73.11}{\scriptstyle\pm0.21}$ & $72.81{\scriptstyle\pm0.21}$ \\
\midrule
& GCN & SGC & LP & C\&S & JK-Net & APPNP & LINK & H2GCN & \textbf{Ordered GNN} \\
arXiv-year & $46.02{\scriptstyle\pm0.26}$ & $32.83{\scriptstyle\pm0.13}$ & $46.07{\scriptstyle\pm0.15}$ & $42.17{\scriptstyle\pm0.27}$ & $46.28{\scriptstyle\pm0.29}$ & $38.15{\scriptstyle\pm0.26}$ & $53.97{\scriptstyle\pm0.18}$ & $49.09{\scriptstyle\pm0.10}$ & $\textbf{54.49}{\scriptstyle\pm0.29}$ \\
\bottomrule
\end{tabular}
}
\end{table}


\subsection{Ablation Study}

We perform an ablation study to evaluate our proposed gating mechanism, Table \ref{table.results.ablation} shows the results. ``Bare GNN'' denotes a simple GNN model with a mean aggregator; ``+simple gating'' denotes adding a gating module where each gate is predicted independently by the gating network, this is equivalent to the Gated GNN \citep{li2015gated} with a mean aggregator; ``++ordered gating'' denotes an ordered gating module as we introduced in Equation \ref{eq.ordered_gating}; ``+++SOFTOR'' indicates our proposed model Ordered GNN. We test on three settings: homophily, heterophily, and over-smoothing. We use CiteSeer for both homophily and over-smoothing settings and Wisconsin for heterophily setting. We set model layers as 8 for the first two settings and 32 for the last one.


\begin{wraptable}{r}{0.55\textwidth}
\vspace{-0.55cm}
\caption{Ablation study}
\label{table.results.ablation}
\centering
\vspace{0.1cm}
\resizebox{0.55\textwidth}{!}{
\begin{tabular}{l|ccc}
\toprule
Model & Homophily & Heterophily & Over-smoothing \\
\midrule
Bare GNN & $73.79{\scriptstyle\pm1.57}$ & $53.92{\scriptstyle\pm5.79}$ & $74.70$ \\
\quad +simple gating & $75.78{\scriptstyle\pm1.49}$ & $86.27{\scriptstyle\pm5.23}$ & $79.30$ \\
\quad \quad ++ordered gating & $76.45{\scriptstyle\pm1.35}$ & $85.69{\scriptstyle\pm4.90}$ & $79.40$ \\
\quad \quad \quad +++SOFTOR & $77.31{\scriptstyle\pm1.73}$ & $88.04{\scriptstyle\pm3.63}$ & $80.80$ \\
\bottomrule
\end{tabular}
}
\vspace{-0.25cm}
\end{wraptable}


As we can see, the ``Bare GNN'' model performs poorly in all three settings, especially in the heterophily setting, revealing the challenge of these settings; the ``+simple gating'' model and the ``++ordered gating'' model perform closely, both providing performance improvements, and let the model comparable to strong baselines; with the help of our proposed differentiable OR module, the ``+++SOFTOR'' model obtains SOTA performance, demonstrating the effectiveness of our proposed methods. Overall, the ablation study verified each component we designed for combine stage of message passing is useful and necessary.


\section{Conclusion}
We propose Ordered GNN, a novel GNN model that aligns the hierarchy of the rooted-tree of a central node with the ordered neurons in its node representation. This alignment implements a more organized message passing mechanism, by ordering the passed information into the node embedding w.r.t. their distance towards the node. This model results to be both more explainable and stronger in performance. It achieves the state-of-the-art on various datasets, ranging across both homophily and heterophily settings. Moreover, we show that the model could effectively prevent over-smoothing, when the model gets very deep. We discuss the connections with related works, the limitation of our model and broader impact on building more powerful GNNs in Appendix \ref{Discussion} due to the space limit.


\section{Ethics Statement}
In this work, we propose a novel GNN model to deal heterophily and over-smoothing problem. Our model integrates inductive bias from graph topology and does not contain explicit assumptions about the labels of nodes on graph, thus avoiding the occurrence of (social) bias. However, the gating mechanism of our model needs to learn from samples with labels, and it is possible to be misleading by targeting corrupted graph structure or features, which can result in negative social impact, because both internet data and social networks can be expressed with graphs. We can develop corresponding ethical guidelines, and constraint the usage of our model to avoid such problems.




\section{Reproducibility Statement}
We provide the introduction for our used datasets and model configuration in Section \ref{exp_setup}, the space of grid search to tune our model and the hardware \& software experimental environment are introduced in Appendix \ref{exp_detail}, we also post every set of hyper-parameters to reproduce our results in Table \ref{table.hyper.hete_homo}, Table \ref{table.hyper.over_smoothing} and Table \ref{table.hyper.larger}, including the baseline model's. The download, load, and pre-process steps are introduced in Appendix \ref{exp_detail}. We've test our model on 3 settings, the setting for heterophily \& homophily are introduced in Section \ref{setting_homo_hete}, the setting for over-smoothing is introduced in Section \ref{setting_smooth}, and the setting for larger benchmarks are introduced in Section \ref{setting_large}. In our released source code, we list the steps needed to reproduce results, and also list all the hyper-parameters with yaml files.

\section*{Acknowledgement}
This work was sponsored by the National Natural Science Foundation of China (NSFC) grant (No.
62106143), and Shanghai Pujiang Program (No. 21PJ1405700).





