# On Performance Discrepancies Across Local Homophily Levels in Graph Neural
Networks

Donald Loveland
University of Michigan, Ann Arbor
<EMAIL_ADDRESS>Zhu
University of Michigan, Ann Arbor
<EMAIL_ADDRESS>Heimann
Lawrence Livermore National Lab
<EMAIL_ADDRESS>Fish
University of Michigan, Ann Arbor
<EMAIL_ADDRESS>T. Schaub
RWTH Aachen University
<EMAIL_ADDRESS>Koutra
University of Michigan, Ann Arbor
<EMAIL_ADDRESS>
###### Abstract

Graph Neural Network (GNN) research has highlighted a relationship between
high homophily (i.e., the tendency of nodes of the same class to connect) and
strong predictive performance in node classification. However, recent work has
found the relationship to be more nuanced, demonstrating that simple GNNs can
learn in certain heterophilous settings. To resolve these conflicting findings
and align closer to real-world datasets, we go beyond the assumption of a
global graph homophily level and study the performance of GNNs when the local
homophily level of a node deviates from the global homophily level. Through
theoretical and empirical analysis, we systematically demonstrate how shifts
in local homophily can introduce performance degradation, leading to
performance discrepancies across local homophily levels. We ground the
practical implications of this work through granular analysis on five real-
world datasets with varying global homophily levels, demonstrating that (a)
GNNs can fail to generalize to test nodes that deviate from the global
homophily of a graph, and (b) high local homophily does not necessarily confer
high performance for a node. We further show that GNNs designed for globally
heterophilous graphs can alleviate performance discrepancy by improving
performance across local homophily levels, offering a new perspective on how
these GNNs achieve stronger global performance.

##  1 Introduction

Deep learning with Graph Neural Networks (GNNs) has become common in many
learning tasks over collaboration networks [1], social networks [2], financial
networks [3], and more [4, 5, 6]. However, given the relative infancy of GNNs,
retrospectives on GNN performance are limited. Understanding the conditions
that will cause a GNN’s performance to degrade promotes a proactive approach
to GNN development, rectifying any issues that may arise once deployed to the
public. One previously studied degradation mechanism of GNN performance is the
presence of heterophilous connections [7, 8]. Heterophily, or the tendency for
nodes of different classes to connect, has been found in a variety of graph
applications where sensitive factors may influence the connective patterns,
necessitating its study [9, 10, 11, 12]. However, more recent retrospectives
have argued that performance does not necessarily degrade with heterophily,
and in fact simple GNN architectures, such as GCN, can perform well in certain
settings [13, 14, 15]. These seemingly conflicting results indicate a gap in
understanding, demanding further research on the influence of heterophily on
GNNs.

To better understand the factors which govern GNN performance, we begin by
exploring the assumptions made in previous works. Surprisingly, many assume
that constituent nodes of a graph possess local homophily levels similar to
the global homophily of the entire graph, leading to a disregard for the
impact of a node’s local homophily level on performance [7, 13, 14]. For the
few works that consider local homophily, it is often assumed that homophilous
nodes should perform better, creating unclear conclusions due to biased
interpretation of the localized results [16, 17]. Moreover, local homophily
has yet to be the focal point of previous works, leading to limited discussion
on why certain patterns emerge. Practically, assuming a constant local
homophily makes it difficult to determine if new models are improving
performance across all nodes, or simply increasing performance for certain
node subsets. Furthermore, by myopically assuming higher homophily is
indicative of higher performance, analysis on discrepancies that arise from
model choice, global homophily, and local homophily is limited, slowing the
development of new GNNs that could address these concerns.

This work. We investigate how shifts in local homophily can impact GNNs,
extending beyond current assumptions and aligning closer to real-world
settings. Our analysis considers a GNN trained on nodes biased towards a
graph’s global homophily, and then applied to test nodes of varying local
homophily levels. We theoretically analyze the scenario by obtaining a closed-
form solution for a GNN’s weight matrix and demonstrate, through perturbation
analysis, that a GNN’s performance can degrade when a node’s local homophily
level is shifted relative to the global homophily of the graph. We also
theoretically demonstrate that when adopting aggregation mechanisms tailored
for heterophilous graphs, predictions are less likely to degrade. We show that
these findings generalize to a variety of settings through a broad empirical
analysis facilitated by our proposed synthetic graph generator that enables
control over the local homophily levels of nodes. We also show the practical
repercussions of our theoretical and empirical analyses on a representative
set of five real-world datasets with varying global homophily levels. Across
our synthetic and real-world datasets, we additionally study nine different
GNN architectures, demonstrating that those tailored to handle heterophily
often maintain more uniform performance, minimizing discrepancies. Together,
our theoretical and empirical analysis describes a new failure point of GNNs –
an expected distribution of labels over a node’s neighbors, stemming from an
over-reliance on the global homophily of a graph – presenting a challenge for
nodes with underrepresented homophily levels to be correctly predicted. While
previous works have noted fairness issues in GNNs based on sensitive
attributes, e.g. race or gender, determined exogenously [18, 19], our results
point to a novel inequality rooted in a network’s structure that could lead to
unfairness in human-centric settings. Our main contributions are:

  * •

Theoretical Analysis: We show how a GNN’s predictions change under a shift in
local homophily level, providing intuition on how GNN performance can degrade
for nodes with local homophily levels which differ from the global homophily
level. We analyze this premise for specific homophilous and heterophilous
designs, showing that heterophilous GNNs will generally be less susceptible to
performance degradation as a node’s local homophily level varies.

  * •

Synthetic Experiments and Model Comparison: We perform empirical analysis by
modifying the preferential attachment model to allow for more granular control
over the distribution of local homophily levels. This capability facilitates
empirical verification of our theory under more general graph structures and
GNN architectures. Additionally, we perform the first node-level analysis that
directly compares GNNs that assume homophily and GNNs that are adjusted for
heterophily, demonstrating different levels of performance discrepancy across
GNN designs.

  * •

Real-world Experiments: We provide the first granular analysis of GNN
performance as local homophily levels are varied across a set of five real-
world datasets. We find that our theoretical performance degradation trends
hold more generally, confirming GNNs designed for heterophily can aid in
minimizing performance discrepancy across nodes with varying local homophily
patterns.

##  2 Related Work

In this section, we begin by discussing GNN architectures designed to improve
learning under heterophily. We then detail previous approaches towards local
property and discrepancy analysis.

Learning GNNs in Diverse Neighborhoods. GNNs adopt an aggregation function to
combine the ego-node’s (the node being updated) features and the neighboring
nodes’ features. Depending on the neighborhood of a node, a particular
aggregation mechanism may be insufficient to learn representations [8]. For
example, GCN [20], GAT [21], and SGC [22] were built to learn over homophilous
neighborhoods through their weighted average of the ego-node and neighboring
nodes’ features. To remedy this issue, models such as GraphSAGE [23], GPR-GNN
[24], FA-GCN [25], GCNII [26], and LINKX [27] separate the ego and neighbor
embeddings, either through a residual connection or concatenation. We note
that LINKX is MLP-based, directly embedding the adjacency matrix through an
MLP rather than traditional message passing. GPR-GNN and FA-GCN additionally
follow a predict-before-propagate paradigm to help alleviate the harm that can
come from mixing representation learning and aggregation [28, 29], while GCNII
utilizes identity mapping to mitigate oversmoothing, a known problem for
homophily [8]. H2GCN adopts further decoupling across higher order
neighborhoods, aggregating each $k$-hop neighborhood separately [7]. In
previous works, there is limited analysis demonstrating the impact of GNN
architectures on the performance of nodes with varying local homophily. In
this work, we provide the first granular analysis of these models, showing a
new perspective on how they perform and their (in)ability to mitigate
performance discrepancies.

Local Property Analysis. Studies on GNN performance relative to an input
graph’s structural properties have gained traction; however, the adjustment to
considering a per-node local perspective is still under-explored. For
instance, many studies have argued the conditions in which a node is able to
benefit from message passing with respect to homophily, but only consider a
constant local homophily level for all nodes [13, 14, 8]. Du et al. offers the
first local analysis, however the results are contradictory across datasets
and are only performed for a single model [16]. More recent work has developed
other homophily-inspired metrics to contextualize local performance, however
the proposed metric can still fail to explain performance depending on the
dataset [17]. Both works assume that higher local homophily should always
improve performance, ultimately guiding their development of new architectures
and metrics that reinforce this assumption. However, the conflicting results
across datasets seen in both works indicate that this assumption may
oversimplify the behavior of a GNN and fail to consider other drivers for
performance degradation. Closely related to our work, Ma et al. analyzes the
disparate treatment of individual nodes defined by their shortest path
distance to the training dataset, showing a degradation in performance as
distance increases [30]. We build upon the idea of structural property
subgroup analysis, but instead consider variations in local homophily rather
than distance to the training set, creating a shift in how performance is
analyzed in the context of homophily. Thus, we analyze the performance of
GNNs, breaking the assumption that the local homophily levels are constant,
and demonstrate how node predictions systematically degrade as the local
homophily levels deviate from the global homophily of the training graph.

##  3 Preliminaries

In this section, we provide key notations and definitions, the notation is
summarized in App. A.1.

###  3.1 Graphs

Let $G=(V,E,\mathbf{X},\mathbf{Y})$ denote a simple graph with node set $V$
and edge set $E$, where $\mathbf{X}\in\mathbb{R}^{|V|\times f}$ represents the
node feature matrix with $f$ features per node and
$\mathbf{Y}\in\\{0,1\\}^{|V|\times c}$ represents the one-hot encoded node
label matrix with $c$ classes. A specific node $i\in G$ has feature vector
$\mathbf{x_{i}}$, class $y_{i}\in\\{1,...,c\\}$, and one-hot encoded class
label vector $\mathbf{y_{i}}$. The edge set can also be represented as an
adjacency matrix, $\mathbf{A}\in\\{0,1\\}^{|V|\times|V|}$, where a value of
$1$ at index $(i,j)$ denotes an edge between nodes $i$ and $j$ in $G$,
otherwise the index is set to $0$. We use both $E$ and $\mathbf{A}$ throughout
the paper, opting for $E$ when explicitly discussing the edges of $G$ and
$\mathbf{A}$ when describing matrix computations on the edge set. A k-hop
neighborhood of node $i\in V$, $N_{k}(i)$, denotes the subgraph induced by the
nodes that are reachable within $k$-steps of $i$.

###  3.2 Node Classification with GNNs

We focus on node classification through a GNN, where the goal is to learn a
mapping between $\mathbf{X}$ and $\mathbf{Y}$. This mapping is estimated
through a subset of $V$, referred to as the set of training nodes $n_{train}$.
For a $k$-layer GNN, learning is facilitated through message passing over
$k$-hop neighborhoods of a graph. The steps, at a high level, include (1)
embedding $\mathbf{X}$ through a non-linear transformation parameterized by a
weight matrix $\mathbf{W}$ and (2) aggregating the embedded features across
neighborhoods of each node. Message passing over all nodes in the graph can be
computed through matrix multiplication, where the most basic formulation
updates node representations through
$\mathbf{R}_{l}=\sigma(\mathbf{(A+I)R}_{l-1}\mathbf{W}_{l})$ for a layer
$l\in\\{1,2,...,k\\}$ of the GNN, where $\mathbf{R}_{0}=\mathbf{X}$ and
$\sigma$ is an activation function. The update is applied $k$ times, resulting
in final representations for each node that can be used for classification.

###  3.3 Homophily and Heterophily

In this work, we focus on edge homophily and present the following definitions
to describe our homophily-based analysis. We begin with the global homophily
ratio of a graph, $h$.
Definition 1 - Global Homophily Ratio. The global homophily ratio $h$ over a
graph’s edge set $E$ is the fraction of edges in $E$ that connect two nodes,
$u$ and $v$, with the same label, $y_{u}$ and $y_{v}$:

$h=\frac{|\\{(u,v):(u,v)\in E\land y_{u}=y_{v}\\}|}{|E|}.$ (1)

The global homophily ratio is used to describe the overall homophily level in
graphs; $h=0$ indicates a fully heterophilous graph and $h=1$ indicates a
fully homophilous graph [10]. Other global homophily metrics are discussed in
Appendix A.3. Additionally, the empirical class compatibility matrix of a
graph, $[\mathbf{H}_{L}]$ describes the probability of two nodes with certain
labels connecting, where the $(u,v)$-th entry is the fraction of edges between
a node in class $u$ and a node in class $v$:
$[\mathbf{H}_{L}]_{u,v}=\dfrac{|\\{(i,j):(i,j)\in E\land y_{i}=u\land
y_{j}=v\\}|}{|\\{(i,j):(i,j)\in E\land y_{i}=u\\}|}$. However, both the global
homophily ratio and compatibility matrix oversimplify the mixing patterns in a
graph when there are varying neighborhood compositions. To perform more
granular analysis on a per-node basis, we also define the local homophily
ratio of a node $t$, $h_{t}$.

Definition 2 - Local Homophily Ratio. The local homophily ratio of a node $t$,
$h_{t}$, is the fraction of edges in the neighborhood of $t$ that connect $t$
to a neighbor $u$ with the same class:

$h_{t}=\frac{|\\{(u,t):(u,t)\in N_{1}(t)\land y_{u}=y_{t}\\}|}{|N_{1}(t)|}.$
(2)

Given GNNs are often shallow and only depend on a small $k$-hop neighborhood
for a single node prediction, it is natural to analyze GNNs through the local,
rather than global, homophily ratio. Moreover, many real-world graphs display
a wide range of local homophily ratios across the constituent nodes, as seen
in App. A.5.2, necessitating local analysis.

##  4 Relationship between a Node’s Local Homophily Level and Performance

In this section, we aim to characterize the impact of local homophily on the
accuracy of node-level predictions by considering shifts in local homophily
levels, at test time, relative to the global homophily level the GNN was
trained on. We begin by revealing the drivers for performance discrepancies
through theoretical analysis and discuss their implications on node-level
performance. Leveraging these insights, we relax our assumptions in Section 9
and show that our theory holds in more general settings via extensive
empirical analysis on synthetic data. Additionally, we consider even more
general real-world graphs (without any constraints) in Section 6.

Setup. Following previous theoretical GNN work [7, 16, 31, 14] and popular
models such as SGC and LightGCN [22, 32], we study two different GNNs with
different aggregation mechanisms. Our homophilous GNN, $F$, is formulated as
$\mathbf{(A+I)XW}$, where $\mathbf{A+I}$ is $G$’s adjacency matrix with self-
loops, directly mixing the features of the ego-node and neighbor nodes [20,
22, 32]. Our heterophilous GNN, $F^{\prime}$ is formulated as
$\mathbf{(X\mathbin{\|}AX)W}$, where $\mathbf{(X\mathbin{\|}AX)}$ is the
concatenation of the ego-node features and aggregated neighboring features
[23, 7, 27]. Both designs have yet to be analyzed through a localized
perspective to characterize discrepancy across predictions.

Similar to the setup in [33], we consider a graph $G$ with a subset of
training nodes $n_{train}$, each of which has an associated node feature
vector $\mathbf{x_{i}}$, one-hot encoded class label vector $\mathbf{y_{i}}$,
1-hop homophily ratio $h$, and degree $d$. For brevity, we focus on binary
classification (though we consider multi-class settings in our experiments)
and represent $\mathbf{y_{i}}$=
$onehot(y_{i})=\begin{bmatrix}1&0\end{bmatrix}$ when node i’s class is
$y_{i}=0$ and $\begin{bmatrix}0&1\end{bmatrix}$ when $y_{i}=1$. We consider
node feature vectors from a uniform distribution and biased towards a
particular class: when
$y_{i}=0,\mathbf{x_{i}}=\begin{bmatrix}(0.5+p)&(0.5-p)\end{bmatrix}$ and when
$y_{i}=1,\mathbf{x_{i}}=\begin{bmatrix}(0.5-p)&(0.5+p)\end{bmatrix}$, where
parameter $p\in[0,0.5]$ controls the ‘agreement’ between the node features and
its class label, i.e. as $p$ approach $0.5$ the features become more similar
to the class labels. The final prediction for node $i$ is
$\mathbf{argmax}\,\mathbf{{z}_{i}}$ where $\mathbf{{z}_{i}}$ is the output
logit vector of the GNNs. We begin by solving for both GNN’s optimal weight
matrix W, and then apply $F$ and $F^{\prime}$ to a test node $t$. We specify
the local homophily ratio for $t$ as $h+\alpha_{t}=h_{t}$, where
$\alpha_{t}\in[-h,1-h]$ is $t$’s shift in local homophily level compared to
the global homophily level. Under this setup, we analyze how $t$’s prediction
is impacted as its local homophily ratio $h_{t}$ shifts relative to $h$, the
global homophily ratio used to train $F$. In Theorem 1, without loss of
generality, we consider the impacts of $h_{t}$ when $t$’s class label is
$y_{t}=0$ for GNN $F$. In Theorem 2, we consider the same setup for GNN
$F^{\prime}$.

Analysis & Implications. Our first theorem provides a direct relationship
between $t$’s performance under a homophilous GNN and its local homophily
ratio as it deviates from the global homophily.

###### Theorem 1

Consider a test node $t$ with local homophily ratio $h+\alpha_{t}$, label
$\mathbf{y_{t}}=\begin{bmatrix}1&0\end{bmatrix}$, and node features
$\mathbf{x_{t}}=\begin{bmatrix}(0.5+p)&(0.5-p)\end{bmatrix}$. The class
prediction from $F$ for node $t$ is a function of the global homophily level
and the shift of the local homophily level, given by
$\mathbf{argmax}\,\mathbf{z_{t}}$, where
$\mathbf{z_{t}}=\mathbf{y_{t}}+b_{1}\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}$
and $b_{1}=d/(1+d(2h-1))$.

Proof.  The proof can be found in App. A.2.1. We provide additional analysis
on a 2-layer variant of $F$, formulated as $\mathbf{(A+I)^{2}XW}$ in App.
A.2.4.

Intuitively, this theorem implies that homophilous GNNs are susceptible to
performance discrepancies. Specifically, we can expect the performance to
degrade when a test node either becomes more homophilous relative to a
heterophilous graph or more heterophilous relative to a homophilous graph. To
further understand the implications of $\alpha_{t}$, we analyze three settings
that naturally arise for the global homophily: (1) $0\leq h<0.5$, (2) $h=0.5$,
and (3) $0.5<h\leq 1$. We note that while the node degrees can influence the
conditions that cause $\alpha_{t}$ to impact $\mathbf{z_{t}}$ (through
$b_{1}$), we show in App. A.2.1 that this mostly occurs for extremely low-
degree nodes. Previous work corroborates our findings regarding the difficulty
with low-degree nodes under heterophily [7, 8, 34]; however, our analysis
extends this observation to demonstrate a significantly more complex interplay
between node degree, global homophily, and shift in local homophily.

Setting 1: Heterophilous $(\mathbf{0\leq\textit{h}<0.5}$): In this scenario,
when $d(2h-1)<-1$, $b_{1}<0$, leading to
$\mathbf{z_{t}}=\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}-\alpha_{t}&\alpha_{t}\end{bmatrix}$,
where $b_{1}$’s sign has been distributed into the vector. Thus,

$\mathbf{z_{t}}=\begin{cases}\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}|\alpha_{t}|&-|\alpha_{t}|\end{bmatrix},&\text{if}\
h_{t}\leq h\\\
\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}-\alpha_{t}&\alpha_{t}\end{bmatrix},&\text{if}\
h_{t}>h,\end{cases}$ (3)

where the sign of $\alpha_{t}$ has been integrated into the vectors. We can
then deduce that (globally) heterophilous graphs, when $b_{1}<0$ is satisfied,
will cause $F$ to degrade in performance as the test node’s local homophily
increases, denoted by the score increase of the wrong class in the second case
of Equation (3). Additionally, when local homophily decreases, the predictions
will improve given the increase in score for the correct class of the first
case in Equation (3), however as $h<0.5$, $\alpha_{t}$ has a smaller possible
range of values, minimizing the impact on $F$’s predictions.

Setting 2: Mixed Homophily ($\mathbf{\textit{h}=0.5}$): When the graph is not
strongly homophilous nor strongly heterophilous (i.e., $h=0.5$), $b_{1}=d$,
leading to:

$\mathbf{z_{t}}=\begin{cases}\mathbf{y_{t}}+d\begin{bmatrix}-|\alpha_{t}|&|\alpha_{t}|\end{bmatrix},&\text{if}\
h_{t}\leq 0.5\\\
\mathbf{y_{t}}+d\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix},&\text{if}\
h_{t}>0.5.\end{cases}$ (4)

In this case, we find that $F$ will have improved performance when the local
homophily of a test node is increased. Conversely, decreased local homophily
for a test node will decrease performance. This is the only case that agrees
with previous work regarding high homophily as a direct indicator of
performance. Notably, the prediction directly depends on $d$, potentially
leading to performance variations that are dominated by degree, rather than
local homophily.

Setting 3: Homophilous ($\mathbf{0.5<\textit{h}\leq 1}$): In this scenario,
$b_{1}>0$, leading to:

$\mathbf{z_{t}}=\begin{cases}\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}-|\alpha_{t}|&|\alpha_{t}|\end{bmatrix},&\text{if}\
h_{t}\leq h\\\
\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix},&\text{if}\
h_{t}>h.\end{cases}$ (5)

We can then deduce that (globally) homophilous graphs will cause $F$ to
degrade in performance as a test node’s local homophily decreases, denoted by
the score increase of the first case of Equation (5). When local homophily
increases, the predictions will improve, however as $h>0.5$, $\alpha_{t}$ has
a smaller range of values, minimizing impact on the predictions.

Going beyond homophilous GNNs, Theorem 2 provides a direct relationship
between $t$’s performance under a heterophilous GNN and its local homophily
ratio as it deviates from the global homophily ratio. At a high level, it
implies that heterophilous GNNs help alleviate performance discrepancies.

###### Theorem 2

Following the same setup from Theorem 1, the class prediction from
$F^{\prime}$ for node $t$ is a function of the global homophily level and the
shift of the local homophily level, given by
$\mathbf{argmax}\,\mathbf{z_{t}}$, where
$\mathbf{z_{t}}=\mathbf{y_{t}}+b^{\prime}_{1}\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}$
and $b^{\prime}_{1}=d^{2}(2h-1)/(1+d^{2}(2h-1)^{2})$.

Proof.  The proof can be found in App. A.2.2.

Comparing to Theorem 1, we can still expect performance degradation, although
at a different rate due to $b^{\prime}_{1}$. Our goal is to determine which of
the two models produces the largest perturbations to the output logit vector,
$\mathbf{z_{t}}$, introducing more significant discrepancies. In Figure 3 in
App. A.2.2, we compare $b^{\prime}_{1}$ and $b_{1}$, the coefficients on the
vector $\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}$ for the
heterophilous and homophilous GNNs, respectively. We find that the magnitude
of $b^{\prime}_{1}$ is often significantly smaller than the magnitude of
$b_{1}$ when varying values of $d$ and $h$. Thus, we can expect
$\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}$ to have less impact on
$\mathbf{z_{t}}$ for heterophilous models, creating less discrepancy as
compared to homophilous models. The full analysis comparing $b_{1}$ and
$b^{\prime}_{1}$ can be found in App. A.2.2.

##  5 Generalization of Theoretical Results via Synthetic Data Analysis

To display how the theoretical relationship between local homophily and
classification performance generalizes, we introduce a graph generator that
enables control over the local homophily ratios across a graph and conduct an
extensive empirical analysis to study the following research questions: (Q1)
What performance discrepancies arise across the range of local homophily
values as the global homophily is varied? and (Q2) Do GNNs built specifically
for heterophily display different discrepancy patterns across local homophily
ranges as compared to homophilous architectures?

###  5.1 Synthetic Data Generation

Building on the preferential attachment model where a compatibility matrix
governs edge likelihood [12, 7], we modify the generator to allow a node’s
homophily level to be either randomly assigned or defined by the compatibility
matrix. We explain the generation process below, and provide the explicit
steps in Algorithm 1, found in App. A.4.2. At a high level, the steps to add a
node to a synthetic graph are: (1) Sample a class label, (2) Generate node
features, and (3) Add connections based on assigned homophily level. The code
is available in an anonymized
repository111https://anonymous.4open.science/r/HeterophilyDiscrepancyGNN-85FB.
We provide related work for graph generation in App. A.4.1, and property
analysis in App. A.4.3.

Class and Feature Generation. For a node $i$, label $y_{i}$ is sampled from
probability distribution $P(\\{0,...,c\\})$ with $c$ possible classes.
Features $\mathbf{x_{i}}$ are generated from a 2D Gaussian, where each
dimension has a mean of $\epsilon y_{i}$ and standard deviation of 1, where
$\epsilon\in[0,1]$ introduces noise into the features. When $\epsilon=1$, the
label $y_{i}$ is explicitly encoded, when $\epsilon=0$, $y_{i}$ is
unrecoverable.

Structure Generation. Similar to previous heterophily analyses, we define a
class compatibility matrix $\mathbf{H}_{L}$ with diagonal elements $h$,
denoting the probability of connecting nodes with similar classes
(homophilous), and off-diagonals elements $(1-h)/c$, denoting the probability
of connecting nodes with different classes (heterophilous) [12, 7]. During
generation, a new node $u$ is attached to an existing node $v$ as
$P\left((u,v)\in E\right)\propto[\mathbf{H}_{L}]_{y_{u},y_{v}}$. To control
the local homophily ratios, we introduce a uniformity parameter $\rho$ such
that with probability $\rho$, a node $i$’s local homophily ratio $h_{i}$ is
sampled at random from a uniform distribution, $U(0,1)$. As $\rho$ increases,
more local homophily ratios follow the random distribution, rather than the
compatibility matrix. Since the preferential attachment model adds nodes
sequentially, it is possible that the local homophily of nodes early in the
generation process drift from their original values. For example, if a node
$i$ is initially generated with two homophilous connections and three
heterophilous connections ($h_{i}=0.4$), later nodes attached to $i$ can
result in a final $h_{i}\neq 0.4$. To correct this, we keep track of how far a
node $i$ has drifted from its original $h_{i}$ through $drift_{i}$, and
prioritize connections to high-drift nodes (i.e., a node $drift_{i}>\delta$,
where $\delta$ is a hyperparameter that defines the drift threshold) that
would return the node’s local homophily ratio back to its original value.

###  5.2 Synthetic Evaluation Setup

We begin by detailing the GNNs and the synthetic graph generation process used
in our experiments.

GNN Models. For our experiments, we consider a diverse set of GNN models
including SGC [22], GCN [20], GAT [21], GraphSAGE [25], GCNII [26], H2GCN [7],
GPR-GNN [24], FA-GCN [25], and LINKX [27]. SGC, GCN and GAT act as homophilous
baselines, while GCNII, H2GCN, GPR-GNN, FA-GCN, and LINKX adopt mechanisms to
improve learning in heterophilous settings. While GraphSAGE was not built for
heterophily, later works [7] have noted its design can improve learning over
heterophily; we include it to represent a midpoint between the models. We also
include results for a graph-agnostic MLP. Experiments across SGC [22], GPR-GNN
[24], and FA-GCN [25] can be found in Appendix A.4.4. All models are
hyperparameter tuned (parameters detailed in App A.6) and the model which
maximizes performance on the validation set is applied to the test set.

Data and Evaluation. Using the proposed generator, we generate graphs by
varying the global homophily ratio $h\in\\{0.1,0.3,0.5,0.7,0.9\\}$ with
$\rho=0.5$, $\epsilon=0.5$, and $\delta=5$, allowing us to study GNNs while
jointly varying global and local homophily levels. An additional study for
$\rho=0.75$ is provided in Figure 9 of App. A.4.4, demonstrating how
discrepancy can be alleviated when the entire range of local homophily ratios
has a sufficient number of nodes. For each combination of $h$ and $\rho$, we
generate 10 graph with 5000 nodes and 100k edges (i.e. n = 5000, m = 20), and
split the nodes into a 50-25-25% split for train, validation, and test.
Additional discussion and experiments on varying the training set size are
provided in Appendix A.4.4. To match our theoretical analysis, we focus on
binary classification. For evaluation, we compare the performance across
models and global homophily ratios as the local homophily ratio is varied. To
compute localized performance, we split the test nodes into four groups based
on their local homophily ratios and calculate an F1 score per group. We choose
F1 score as the node subsets can become imbalanced with respect to their class
after splitting. Per bin, we report the average and standard deviation for F1
over the 10 generated graphs. Results are presented in Figures 1 of the main
text, and Figures 7, 8 in App. A.4.4

Figure 1: Performance for different GNN models (columns) on synthetic data
generated with global homophily ratios $h\in\\{0.1,0.5,0.9\\}$ (rows) and
uniformity $\rho=0.5$, reported for different local homophily ranges. As the
local homophily deviates from the global homophily, homophilous models (blue
bars) fail to generalize, resulting in performance discrepancy, while
heterophilous models (red bars) retain higher relative performance by up to
0.5 F1-score, alleviating discrepancy. Exact differences can be found in
Figure 6. MLP achieves constant performance.

###  5.3 Synthetic Results

(Q1) Performance discrepancies across local and global homophily ratios. In
Figure 1, we visualize the performance of each model by binning nodes of
similar homophily for each of the global homophily ratios (analysis on
$h\in\\{0.3,0.7\\}$ and additional GNNs in App A.4.4). We first find that the
MLP, as expected, performs significantly worse with an F1 score of 0.75 as
compared to the GNN’s which achieve an F1 score of 0.9 and higher,
highlighting the positive contribution of the graph structure on predictions.
Next, we highlight the clear trend between performance and local homophily
depending on the global homophily ratio: when $h<0.5$ (lower global
homophily), performance often degrades as local homophily increases, while
when $h>0.5$ (higher global homophily), performance often degrades as local
homophily decreases. The results align well with our theoretical analysis,
showing strong generalization of our findings to more complex graph structures
and GNN architectures. Together, our theoretical and empirical results
indicate that assuming high homophily always indicates high performance may
oversimplify the GNN’s behavior, leading retrospective analyses astray when
diagnosing performance degradation.

(Q2) Performance discrepancies for homophilous vs. heterophilous GNNs. To
understand how different GNN designs amplify or reduce discrepancy, we analyze
the trend of local performances across models. While all models achieve
similar global performance, we observe that the models perform differently
depending on the local homophily range. As shown in Figure 1 (and Figure 6 in
App. A.4.4), homophilous models often have higher performance for test nodes
with local homophily levels that are close to the global homophily, as
expected from Theorem 1, whereas heterophily-based models perform better than
homophilous models for nodes with local homophily levels far from the global
homophily, agreeing with Theorem 2. These insights highlight the different
behaviors of the two designs, indicating that models designed for heterophily
are able to alleviate performance discrepancies across nodes, while
homophilous models exacerbate discrepancy in strongly
homophilous/heterophilous graphs,  Interestingly, LINKX, despite being MLP
based, is competitive as compared to the other heterophilous GNNs in
minimizing discrepancy. We attribute this to the separation of ego-node and
neighbor-node embeddings, a key design in the other heterophilous models [7,
35]. These results verify that heterophilous models offer a better performance
trade-off between nodes with over- and under-represented local homophily
levels, displaying minor degradation on the over-represented nodes and
significant improvement on the under-represented nodes.

##  6 Real-world Empirical Evaluation

We now demonstrate how our results extend to real-world datasets.

Data and Setup. We choose five real-world graphs: two homophilous (Cora [36],
Coauthor-CS [37]) and three heterophilous (Wisconsin [38], Squirrel [39], and
Arxiv-Year [27]). Each heterophilous dataset varies greatly in size, providing
an opportunity to measure the impact of the number of nodes and density on
discrepancy. Additionally, we choose Wisconsin and Squirrel due to their
historically inconclusive performance when comparing GNNs to non-graph-based
deep learning baselines [7, 13]. Our analysis aims to provide insight into why
prior results have been inconclusive and explain their poor performance from a
local perspective. For Wisconsin and Squirrel, we perform 30 random 50-25-25%
splits of the nodes to obtain the train, validation, and test sets. For Arxiv-
Year, Cora, and Coauthor-CS we perform 5 random splits over the same ratio as
the graph are much larger. We train the same architectures as in the synthetic
experiments, seen in Section 5.2.

We report results in Figure 2 of the main text, and Figure 11 in App. 11.
While we group all test nodes into four local homophily ranges, we limit the
ranges of Wisconsin in Figure 2 due to having less than three test nodes in
local homophily ranges above 0.6. Additionally, features may be more
informative in certain local homophily ranges, obscuring when performance
degrades due to homophily or uninformative features. Thus, we report the
difference in F1 score between the different GNNs and the graph-agnostic MLP,
$\Delta$F1 $\equiv$ (F1-score GNN - F1-score MLP), to disentangle how
performance relies on the node features as compared to the graph structure.

###  6.1 Performance Discrepancies in Homophilous vs. Heterophilous Real-World
Datasets

Homophilous Graphs. The homophilous graphs are shown in the two right-most
plots of Figure 2. First, we highlight the nearly 0.6 drop in $\Delta$F1
across both datasets for test nodes with local homophily ratios far from the
global homophily ratio, performing worse than an MLP as denoted by the
negative F1 score difference. Furthermore, this degradation is consistent
across all of the GNN architectures, with H2GCN, which leverages heterophilous
designs, maintaining generally higher performance as the nodes become more
heterophilous. This result demonstrates the practical implications of our
theoretical analysis, allowing us to additionally demonstrate that degradation
can occur under rich feature sets. Notably, the MLP outperforms all of the GNN
architectures in the heterophilous parts of the graph, implying the
neighborhood information is actively corrupting performance in these regions
due to the GNN’s reliance on the global homophily.

Heterophilous Graphs. The heterophilous datasets are shown on the three left-
most plots of Figure 2. Compared to the homophilous datasets, Wisconsin and
Squirrel have poor performance relative to the MLP for nearly all local
homophily ranges. Our local analysis identifies a mechanism for how this
arises: exacerbated performance degradation for nodes with local homophily
ratios far from the global homophily ratio. Despite the MLP model
outperforming the GNN models across Wisconsin and Squirrel, there is still a
notable degradation in performance as the local homophily ranges increase,
causing a 0.2-0.4 drop in $\Delta$F1 between the bins furthest from the global
homophily ratio. We hypothesize that degree can influence whether GNNs degrade
in performance on heterophilous graphs, as seen in our theoretical analysis in
App. A.2.1, causing this drop to be lower as compared to the homophilous
graphs. Aligning closer to the homophilous datasets, each model’s best
performance on Arxiv-Year is within the bin closest to the global homophily
ratio, demonstrating a scenario where the heterophilous neighborhood is
beneficial to learning. Following the trends seen in the synthetic data, the
heterophilous models also display higher uniformity across the homophily
levels for each datasets. Surprisingly, for Squirrel, we observe a significant
difference in MLP performance across the various local homophily ranges,
indicating that the features are intrinsically more informative in certain
regions. When accounting for this through $\Delta$F1, we see similar
degradation to that in the synthetic and homophilous datasets, again with the
heterophilous GNNs tending to perform best.

Figure 2: Real-words graphs: Difference in performance, $\Delta$F1=F1(GNN) -
F1(MLP), for GNN models across ranges of local homophily ratios (more GNNs in
App. 11), averaged over multiple splits. Our results elucidate how models
achieve different global performance, where heterophilous models (bars in
reds) better tend to combat the systematic performance discrepancy seen in
homophilous models (bars in blues). Gray indicates the range that the global
homophily ratio falls in; negative bars indicate worse performance than MLP.

###  6.2 Discussion

Two key insights emerge that have yet to be established in other analyses: (a)
nodes with higher local homophily are not inherently easier to classify, as
seen by the relative drops in performance on heterophilous datasets, and (b)
heterophilous designs improve learning across nearly all local homophily
ranges, not just one particular range, alleviating performance discrepancies
across nodes. We note that there might be additional factors which give rise
to performance variations across local homophily ratios. For instance,
previous works have identified that, under certain settings, nodes with
extreme homophily levels can be easier to classify [13, 16, 17]. Moreover, the
ease of classification has been tied to the interaction of degree and
homophily, pointing towards high degree nodes as the easiest setting [8].
Performance variations due to raw structural properties (e.g. degree or local
homophily ratio) and local homophily shift are not necessarily mutually
exclusive, nor do they conflict with our results. Instead, we hypothesize both
factors can interact, amplifying performance disparities across homophily
ranges, further necessitating future work which studies their interplay.

##  7 Conclusion

In this work, we take a local perspective focused on discerning how nodes of
varying local homophily levels can experience performance discrepancies. We
first theoretically demonstrated that classification performance degrades as
the local homophily ratio of a node deviates from the global homophily ratio.
To demonstrate the generalizability of our findings, we proposed a new
parsimonious synthetic graph generator that allows generating graphs with
varying global and local homophily. We demonstrated that our theoretical
insights still hold in more general settings, finding that performance
degradation can occur in either highly homophilous or heterophilous settings.
Furthermore, we showed that this disparity in performance can be reduced by
using GNN models which adopt explicit mechanisms to support heterophily,
providing insight into how GNNs with heterophilous designs improve performance
globally. Our experiments on real-world datasets of varying global homophily
ratios confirm the practical implications of our insights, exhibiting similar
disparity patterns. The discovery and characterization of GNN degradation
through shifts in local homophily relative to a graph’s global homophily
necessitates the development of new GNNs that are able to explicitly handle
such data shifts. Additionally, our findings highlight how, in human-facing
applications of GNNs, individuals might experience disparate treatment under a
GNN due to structural properties of the underlying graph, opening new research
directions in algorithmic fairness.

## References

  * Mccallum et al. [2000] Andrew Mccallum, Kamal Nigam, and Jason Rennie.  Automating the construction of internet portals.  03 2000.
  * Fan et al. [2019] Wenqi Fan, Yao Ma, Qing Li, Yuan He, Eric Zhao, Jiliang Tang, and Dawei Yin.  Graph neural networks for social recommendation.  In _WWW_ , 2019.
  * [3] Daixin Wang, Zhiqiang Zhang, Jun Zhou, Peng Cui, Jingli Fang, Quanhui Jia, Yanming Fang, and Yuan Qi.  Temporal-Aware Graph Neural Network for Credit Risk Prediction.  In _SDM_.
  * Ciotti et al. [2015] Valerio Ciotti, Moreno Bonaventura, Vincenzo Nicosia, Pietro Panzarasa, and Vito Latora.  Homophily and missing links in citation networks.  _EPJ Data Science_ , 5, 11 2015.  doi: 10.1140/epjds/s13688-016-0068-2.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang.  Connecting the dots: Multivariate time series forecasting with graph neural networks.  _In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2020.
  * Vashishth et al. [2020] Shikhar Vashishth, Soumya Sanyal, Vikram Nitin, and Partha Talukdar.  Composition-based multi-relational graph convolutional networks, 2020\.
  * Zhu et al. [2020] Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra.  Beyond homophily in graph neural networks: Current limitations and effective designs.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Yan et al. [2022] Yujun Yan, Milad Hashemi, Kevin Swersky, Yaoqing Yang, and Danai Koutra.  Two sides of the same coin: Heterophily and oversmoothing in graph convolutional neural networks.  _IEEE ICDM_ , 2022.
  * Hofstra et al. [2017] Bas Hofstra, Rense Corten, Frank Van Tubergen, and Nicole B Ellison.  Sources of segregation in social networks: A novel approach using facebook.  _American Sociological Review_ , 82(3):625–656, 2017.
  * [10] Miller McPherson, Lynn Smith-Lovin, and James M Cook.  Birds of a feather: Homophily in social networks.  _Annual review of sociology_.
  * Brea et al. [2018] Jorge Brea, Javier Burroni, and Carlos Sarraute.  Inference of users demographic attributes based on homophily in communication networks.  _NetMob_ , 2018.
  * Karimi et al. [2018] Fariba Karimi, Mathieu Génois, Claudia Wagner, Philipp Singer, and Markus Strohmaier.  Homophily influences ranking of minorities in social networks.  _Scientific reports_ , 2018.
  * Ma et al. [2022] Yao Ma, Xiaorui Liu, Neil Shah, and Jiliang Tang.  Is homophily a necessity for graph neural networks?  _ICLR_ , 2022.
  * Luan et al. [2021] Sitao Luan, Chenqing Hua, Qincheng Lu, Jiaqi Zhu, Mingde Zhao, Shuyuan Zhang, Xiao-Wen Chang, and Doina Precup.  Is heterophily a real nightmare for graph neural networks to do node classification?  _ICLR 2022_ , 2021.
  * [15] Jiong Zhu, Yujun Yan, Mark Heimann, Lingxiao Zhao, Leman Akoglu, and Danai Koutra.  Heterophily and graph neural networks: Past, present and future.  _IEEE Data Engineering Bulletin_.  URL https://par.nsf.gov/biblio/10435541.
  * Du et al. [2022] Lun Du, Xiaozhou Shi, Qiang Fu, Xiaojun Ma, Hengyu Liu, Shi Han, and Dongmei Zhang.  Gbk-gnn: Gated bi-kernel graph neural networks for modeling both homophily and heterophily.  In _Proceedings of the ACM Web Conference 2022_ , pages 1550–1558, 2022.
  * Cavallo et al. [2022] Andrea Cavallo, Claas Grohnfeldt, Michele Russo, Giulio Lovisotto, and Luca Vassio.  2-hop neighbor class similarity (2ncs): A graph structural metric indicative of graph neural network performance.  _arXiv preprint arXiv:2212.13202_ , 2022.
  * Dong et al. [2022] Yushun Dong, Jing Ma, Chen Chen, and Jundong Li.  Fairness in graph mining: A survey.  _arXiv preprint arXiv:2204.09888_ , 2022.
  * Zhang et al. [2022] Wenbin Zhang, Jeremy C Weiss, Shuigeng Zhou, and Toby Walsh.  Fairness amidst non-iid graph data: A literature review.  _arXiv preprint arXiv:2202.07170_ , 2022.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _ICLR_ , 2017.
  * Veličković et al. [2018] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio.  Graph attention networks.  In _ICLR_ , 2018.
  * Wu et al. [2019] Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger.  Simplifying graph convolutional networks.  In _In ICML_ , 2019.
  * Hamilton et al. [2017] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In _In NeurIPS_ , 2017.
  * Chien et al. [2021] Eli Chien, Jianhao Peng, Pan Li, and Olgica Milenkovic.  Adaptive universal generalized pagerank graph neural network.  _In ICLR_ , 2021.
  * Bo et al. [2021] Deyu Bo, Xiao Wang, Chuan Shi, and Huawei Shen.  Beyond low-frequency information in graph convolutional networks.  In _In AAAI_ , 2021.
  * Chen et al. [2020] Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li.  Simple and deep graph convolutional networks.  In _ICML_ , 2020.
  * Lim et al. [2021] Derek Lim, Felix Hohne, Xiuyu Li, Sijia Linda Huang, Vaishnavi Gupta, Omkar Bhalerao, and Ser Nam Lim.  Large scale learning on non-homophilous graphs: New benchmarks and strong simple methods.  _Advances in Neural Information Processing Systems_ , 34:20887–20902, 2021.
  * Cui et al. [2020] Ganqu Cui, Jie Zhou, Cheng Yang, and Zhiyuan Liu.  Adaptive graph encoder for attributed graph embedding.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2020.
  * Gasteiger et al. [2019] Johannes Gasteiger, Aleksandar Bojchevski, and Stephan Günnemann.  Predict then propagate: Graph neural networks meet personalized pagerank.  _In ICLR_ , 2019.
  * Ma et al. [2021] Jiaqi Ma, Junwei Deng, and Qiaozhu Mei.  Subgroup generalization and fairness of graph neural networks.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Wang et al. [2022] Yu Wang, Yuying Zhao, Yushun Dong, Huiyuan Chen, Jundong Li, and Tyler Derr.  Improving fairness in graph neural networks via mitigating sensitive attribute leakage.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , KDD ’22, page 1938–1948, New York, NY, USA, 2022\. Association for Computing Machinery.  ISBN 9781450393850.  doi: 10.1145/3534678.3539404.  URL https://doi.org/10.1145/3534678.3539404.
  * He et al. [2020] Xiangnan He, Kuan Deng, Xiang Wang, Yan Li, Yongdong Zhang, and Meng Wang.  Lightgcn: Simplifying and powering graph convolution network for recommendation, 2020.
  * Zhu et al. [2022] Jiong Zhu, Junchen Jin, Donald Loveland, Michael T Schaub, and Danai Koutra.  On the relationship between heterophily and robustness of graph neural networks.  _KDD_ , 2022.
  * Lin et al. [2023] Bei Lin, You Li, Ning Gui, Zhuopeng Xu, and Zhiwu Yu.  Multi-view graph representation learning beyond homophily.  _ACM Trans. Knowl. Discov. Data_ , 17(8), jun 2023.  ISSN 1556-4681.  URL https://doi.org/10.1145/3592858.
  * Platonov et al. [2023] Oleg Platonov, Denis Kuznedelev, Michael Diskin, Artem Babenko, and Liudmila Prokhorenkova.  A critical look at the evaluation of gnns under heterophily: are we really making progress?  _ICLR_ , 2023.
  * Sen et al. [2008] Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad.  Collective classification in network data.  _AI Magazine_ , 29(3):93, Sep. 2008.  doi: 10.1609/aimag.v29i3.2157.  URL https://ojs.aaai.org/aimagazine/index.php/aimagazine/article/view/2157.
  * Bojchevski and Günnemann [2017] Aleksandar Bojchevski and Stephan Günnemann.  Deep gaussian embedding of graphs: Unsupervised inductive learning via ranking.  _arXiv preprint arXiv:1707.03815_ , 2017.
  * Pei et al. [2020a] Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang.  Geom-gcn: Geometric graph convolutional networks.  In _International Conference on Learning Representations_ , 2020a.
  * Rozemberczki et al. [2021] Benedek Rozemberczki, Carl Allen, and Rik Sarkar.  Multi-scale attributed node embedding.  _Journal of Complex Networks_ , 9(2):cnab014, 2021\.
  * Palowitch et al. [2022] John Palowitch, Anton Tsitsulin, Brandon Mayer, and Bryan Perozzi.  Graphworld: Fake graphs bring real insights for gnns.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , KDD ’22, page 3691–3701, New York, NY, USA, 2022\. Association for Computing Machinery.  ISBN 9781450393850.  doi: 10.1145/3534678.3539203.  URL https://doi.org/10.1145/3534678.3539203.
  * Wang et al. [2021] Chaokun Wang, Binbin Wang, Bingyang Huang, Shaoxu Song, and Zai Li.  Fastsgg: Efficient social graph generation using a degree distribution generation model.  In _2021 IEEE 37th International Conference on Data Engineering (ICDE)_ , pages 564–575, 2021.  doi: 10.1109/ICDE51399.2021.00055.
  * Maekawa et al. [2023] Seiji Maekawa, Yuya Sasaki, George Fletcher, and Makoto Onizuka.  Gencat: Generating attributed graphs with controlled relationships between classes, attributes, and topology.  _Information Systems_ , 115:102195, 2023.  ISSN 0306-4379.  doi: https://doi.org/10.1016/j.is.2023.102195.  URL https://www.sciencedirect.com/science/article/pii/S0306437923000315.
  * Barabasi and Albert [1999] Albert-Laszlo Barabasi and Reka Albert.  Emergence of scaling in random networks.  _Science_ , 286(5439):509–512, Oct. 1999.
  * Pei et al. [2020b] Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang.  Geom-gcn: Geometric graph convolutional networks.  2020b.  URL https://openreview.net/forum?id=S1e2agrFvS.
  * Hu et al. [2020] Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  _Advances in neural information processing systems_ , 33:22118–22133, 2020.

##  Appendix A Appendix

###  A.1 Notation Table

In this section, we give a summary of the notation used throughout the main
text.

Table 1: Graph- and Node-level Notations and Definitions Notation |  Definitions
---|---
$G$ |  Graph $G$
$V$ |  Node set $V$ for $G$
$E$ |  Edge set $E$ for $G$
$\mathbf{I}$ |  Identity matrix
$\mathbf{A}$ |  Adjacency matrix for $G$. Shape $|V|\times|V|$
$\mathbf{X}$ |  Feature matrix for $G$. Shape $|V|\times f$ where $f$ is the number of features.
$\mathbf{Y}$ |  One-hot encoded class matrix for $G$. Shape $|V|\times c$ where $c$ is the number of classes.
$\mathbf{H}_{L}$ |  Class compatibility matrix..
$h$ |  Global homophily ratio.
$i$ |  A node $i\in G$
$\mathbf{x_{i}}$ |  Feature vector for node $i$
$y_{i}$ |  Class label for node $i$
$\mathbf{y_{i}}$ |  One-hot encoded class label vector for node $i$
$N_{k}(i)$ |  $k$-hop neighborhood for a node $i$
$h_{i}$ |  Local homophily ratio for a node $i$
$\alpha_{i}$ |  Shift in local homophily ratio for a node $i$, relative to the global homophily ratio of its graph, $h+\alpha_{i}=h_{i}.$
Table 2: Graph Neural Network Notations and Definitions Notation |  Definitions
---|---
$n_{train}$ |  Training nodes of a GNN
$\mathbf{W_{l}}$ |  The weight matrix for a layer $l$ of a GNN
$\mathbf{R_{l}}$ |  The hidden representations for nodes at layer $l$ of a GNN. $R_{0}$ denotes the input features.
$\mathbf{r_{i}^{l}}$ |  The hidden representation vector for a node $i$ at layer $l$
$\sigma$ |  Activation function of a GNN.
$\mathbf{z_{i}}$ |  Output logit vector of a GNN for node $i$.
Table 3: Synthetic Graph/Evaluation Notations and Definitions Notation |  Definitions
---|---
$n$ |  Number of nodes
$m$ |  Number of edges added per generation step
$\epsilon$ |  Noise strength for generated features.
$\rho$ |  Probability to sample a new homophily ratio.
$\delta$ |  Hyperparameter to determine if a node has drifted in homophily level.
$\Delta F1$ |  F1 Score of GNN - F1 Score of MLP (graph agnostic).

###  A.2 Theoretical Relationship between a Node’s Local Homophily Level and
Performance

In this section, we provide the proof for Theorem 1 in the main text. We
additionally perform theoretical analysis on the influence of degree on
Theorem 1, as well as the impact of higher order aggregation. We provide a
high-level takeaway for each additional analysis.

####  A.2.1 Proof for Theorem 1: Homophilous Model $\mathbf{(A+I)XW}$

As described in the main text, we consider a graph $G$ with a subset of nodes
$n_{train}$. Each node $i\in n_{train}$ has an associated node feature vector
$\mathbf{x_{i}}$, one-hot encoded class label vector $\mathbf{y_{i}}$, 1-hop
homophily ratio $h$, and degree $d$. For brevity, we focus on binary
classification, and represent $\mathbf{y_{i}}$ as
$\begin{bmatrix}1&0\end{bmatrix}$ when $y=0$ and
$\begin{bmatrix}0&1\end{bmatrix}$ when $y=1$. The node feature vectors are
initialized as $\mathbf{x_{i}}=\begin{bmatrix}(0.5+p)&(0.5-p)\end{bmatrix}$
when $y_{i}=0$ and
$\mathbf{x_{i}}=\begin{bmatrix}(0.5-p)&(0.5+p)\end{bmatrix}$ when $y_{i}=1$.
The GNN, $F$, is formulated as $\mathbf{(A+I)XW}$, where $\mathbf{A+I}$ is
$G$’s adjacency matrix with self-loops and $\mathbf{W}$ is $F$’s weight
matrix, and trained through $n_{train}$. The final prediction for a node $i$
is $\mathbf{argmax}\,\mathbf{{z}_{i}}$ where $\mathbf{{z}_{i}}$ is the output
logit vector of $F$. Focusing first on the (A + I)X term, a row (node) $i$ in
this matrix with a class 0 will be equal to:

$\displaystyle
y_{i}=0:\mathbf{r_{i}^{1}}=\mathbf{x_{i}}+\dfrac{hd}{2}\begin{bmatrix}1+p&1-p\end{bmatrix}+\dfrac{(1-h)d}{2}\begin{bmatrix}1-p&1+p\end{bmatrix}$
(6)

similarly, a row (node) $i$ in this matrix of class 1 will be equal to:

$\displaystyle
y_{i}=1:\mathbf{r_{i}^{1}}=\mathbf{x_{i}}+\dfrac{hd}{2}\begin{bmatrix}1-p&1+p\end{bmatrix}+\dfrac{(1-h)d}{2}\begin{bmatrix}1+p&1-p\end{bmatrix}$
(7)

The computations for each $\mathbf{r_{i}^{1}}$ represent the sum of the node
vector for node $i$, the feature vectors from the homophilous connections of
node $i$, and the feature vectors of the heterophilous connections for node
$i$. Condensing these terms, the simplified expressions for each
$\mathbf{r_{i}^{1}}$ is shown in Equations 8 and 9

$\displaystyle
y_{i}=0:\mathbf{r_{i}^{1}}=\dfrac{1}{2}\begin{bmatrix}(1+d+p+dp(2h-1))&(1-p+d(1+p-2hp))\end{bmatrix}$
(8)

$\displaystyle
y_{i}=1:\mathbf{r_{i}^{1}}=\dfrac{1}{2}\begin{bmatrix}(1-p+d(1+p-2hp))&(1+d+p+dp(2h-1))\end{bmatrix}$
(9)

We can then represent (A+I)XW as RW as seen in Equation 10 and set it equal to
the label matrix $\mathbf{Y}$ after a transformation by the weight matrix W.

$\displaystyle\textbf{R}=\dfrac{1}{2}\left[{\begin{array}[]{cc}\vdots&\vdots\\\
(1+d+p+dp(2h-1))&(1-p+d(1+p-2hp))\\\ \vdots&\vdots\\\
(1-p+d(1+p-2hp))&(1+d+p+dp(2h-1))\\\ \vdots&\vdots\\\
\end{array}}\right]_{|T|\times 2}$ (10)

(11) $\displaystyle\textbf{R}_{|T|\times 2}\textbf{W}_{2\times
2}=\textbf{Y}_{|T|\times 2}=\begin{bmatrix}\vdots&\vdots\\\ 1&0\\\
\vdots&\vdots\\\ 0&1\\\ \vdots&\vdots\\\ \end{bmatrix}$

Our goal is to now solve for the optimal W in this system of equations. While
the system of equations is overdetermined, each row of the same class share
the same solution. Thus, we can sample the unique data points, leaving us with
the optimal W in Equation 12 from solving the reduced system of equations,
where $c_{1}=\dfrac{1}{2(p-d^{2}p+2dhp+2d^{2}hp)}$.

$\displaystyle\textbf{W}=c_{1}\left[{\begin{array}[]{cc}1+d+p+dp(2h-1)&-1+p-d(1+p-2hp)\\\
-1+p-d(1+p-2hp)&1+d+p+dp(2h-1)\\\ \end{array}}\right]$ (12)

We now consider a new test point, $t$ with label $y_{t}=0$ and subsequently
features $\mathbf{x_{t}}=\dfrac{1}{2}\begin{bmatrix}(1+p)&(1-p)\end{bmatrix}$.
In addition, we assume that $t$ has a homophily ratio
$h+\alpha_{t}=h_{t}\in[0,1]$ where $\alpha_{t}$ represents a shift away from
$h$. When $\alpha_{t}$ is positive, we can interpret this as a test node with
a higher homophily ratio, and when $\alpha_{t}$ is negative, a lower homophily
ratio. Under this new homophily ratio, $\mathbf{r_{t}^{1}}$ can be calculated
as:

$\displaystyle\mathbf{r_{t}^{1}}=\dfrac{h_{t}d+1}{2}\begin{bmatrix}1+p&1-p\end{bmatrix}+\dfrac{(1-h_{t})d}{2}\begin{bmatrix}1-p&1+p\end{bmatrix}$
(13)

Similar to equation 8, the simplified expressions for $\mathbf{r_{t}^{1}}$ is
given as:

$\displaystyle\mathbf{r_{t}^{1}}=\dfrac{1}{2}\begin{bmatrix}(1+d+p+dp(2h_{t}-1))&(1-p+d(1+p-2h_{t}p))\end{bmatrix}$
(14)

We can now compute $\mathbf{r_{t}^{1}}\textbf{W}$ and analyze the associated
predictions. The output of $\mathbf{r_{t}^{1}}\textbf{W}$ is:

$\mathbf{z_{t}}=\dfrac{1}{1+d(2h-1)}\begin{bmatrix}(1+d(h+h_{t}-1))&d(h-h_{t})\end{bmatrix}$
(15)

It is easy to check that when $h_{t}=h$, we recover the correct prediction of
$[1,0]$ as expected. We now consider what happens when $h_{t}=h+\alpha_{t}$
where $\alpha_{t}\in[-h,1-h]$. Since we are interested in understanding how
the predictions would change as a byproduct of $\alpha_{t}$, we look at the
difference between predictions when $h_{t}=h$. Below we compute
$\Delta\mathbf{z_{t}}$, the change in the logit vector output for a test point
with local homophily $h_{t}$, where $b_{1}=\dfrac{d}{1+d(2h-1)}$.

$\displaystyle\Delta\mathbf{z_{t}}$
$\displaystyle=\dfrac{b_{1}}{d}\begin{bmatrix}(1+d(h+(h+\alpha_{t})-1))&d(h-(h+\alpha_{t}))\end{bmatrix}$
$\displaystyle-\dfrac{b_{1}}{d}\begin{bmatrix}(1+d(h+h-1))&d(h-h)\end{bmatrix}$
(16) $\displaystyle\Delta\mathbf{z_{t}}$
$\displaystyle=b_{1}\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}\quad\square$
(17)

####  A.2.2 Proof for Theorem 2: Heterophilous Model
$\mathbf{(X\mathbin{\|}AX)W}$

Following the same set up previously described for the GNN parameterized as
$\mathbf{(A+I)XW}$, we also consider a GNN, $F^{\prime}$, parameterized as
$\mathbf{(X\mathbin{\|}AX)W}$, concatenating the features from the ego-node
and the aggregated neighboring features. This separation through concatenation
is found in a variety of GNNs, and is often denoted as a core design to
promote heterophilous learning [23, 7].

Focusing first on the aggregation step of $F^{\prime}$, $\mathbf{AX}$, a node
$i$ of class $0$ will have aggregated features from its neighbors,
$\mathbf{r_{i}^{1}}$, defined as:

$\displaystyle
y_{i}=0:\mathbf{r_{i}^{1}}=\dfrac{hd}{2}\begin{bmatrix}1+p&1-p\end{bmatrix}+\dfrac{(1-h)d}{2}\begin{bmatrix}1-p&1+p\end{bmatrix}$
(18)

similarly, a node $i$ of class 1 will have aggregated features from its
neighbors:

$\displaystyle
y_{i}=1:\mathbf{r_{i}^{1}}=\dfrac{hd}{2}\begin{bmatrix}1-p&1+p\end{bmatrix}+\dfrac{(1-h)d}{2}\begin{bmatrix}1+p&1-p\end{bmatrix}$
(19)

Then, the concatenated vectors $\mathbf{(X\mathbin{\|}AX)W}$ for nodes of
class 0 and 1 are expressed as:

$\displaystyle
y_{i}=0:\mathbf{x_{i}\mathbin{\|}r_{i}^{1}}=\dfrac{1}{2}\begin{bmatrix}1+p&1-p&d+dp(2h-1)&d(1+p-2hp)\end{bmatrix}$
(20) $\displaystyle
y_{i}=1:\mathbf{x_{i}\mathbin{\|}r_{i}^{1}}=\dfrac{1}{2}\begin{bmatrix}1-p&1+p&d(1+p-2hp)&d+dp(2h-1)\end{bmatrix}$
(21)

We refer to the collection of concatenated vectors as
$\mathbf{R}\in\mathcal{R}^{|V|\times 4}$. Then, the goal, similar to the case
of the GNN $F$, is to solve for $W$ such that:

(22) $\displaystyle\textbf{R}_{|T|\times 4}\textbf{W}_{4\times
2}=\textbf{Y}_{|T|\times 2}=\begin{bmatrix}\vdots&\vdots\\\ 1&0\\\
\vdots&\vdots\\\ 0&1\\\ \vdots&\vdots\\\ \end{bmatrix}$

Following a similar procedure to the analysis of $F$, we reduce the system of
equations by sampling the unique points and solve for the optimal
$\mathbf{W}$. However, unlike in $F$, $\mathbf{R}$ is not a square matrix
after sampling, and thus we use a right pseudo-inverse through
$\mathbf{R^{-1}=R^{T}(RR^{T})^{-1}}$ to attain a solution to the system of
equations.

We now consider a new test point, $t$ with label $y_{t}=0$, similar to the set
up for $F$. Then:

$\displaystyle\mathbf{x_{t}\mathbin{\|}r_{t}^{1}}=\dfrac{1}{2}\begin{bmatrix}1+p&1-p&d+dp(2h_{t}-1)&d(1+p-2h_{t}p)\end{bmatrix}$
(23)

We can now compute $\mathbf{r_{t}^{1}}\mathbf{W}=\mathbf{z_{t}}$:

$\displaystyle\mathbf{z_{t}}=\dfrac{1}{1+d^{2}(2h-1)^{2}}\begin{bmatrix}1+d^{2}(2h-1)(h+h_{t}-1))&(d^{2}(2h-1)(h-h_{t}))\end{bmatrix}$
(24)

As before, when $h_{t}=h$, it is easy to see that the label is recovered. We
now consider when $h_{t}=h+\alpha_{t}$ and look at the difference in
predictions relative to when $h_{t}=h$. $\Delta\mathbf{z_{t}}$ for a test
point with local homophily $h_{t}$ predicted by the heterophilous model, where
$b^{\prime}_{1}=\dfrac{d^{2}(2h-1)}{1+d^{2}(2h-1)^{2}}$ is:

$\displaystyle\Delta\mathbf{z_{t}}$
$\displaystyle=b^{\prime}_{1}\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}\quad\square$
(25)

Thus, we can still expect heterophilous models which adopt the ego-neighbor
separation design to degrade in performance as the local homophily of a node
changes, just as studied in the Section 4 of the main text. However, the
change in coefficient, $b_{1}$ vs. $b^{\prime}_{1}$, demonstrates a
significantly different degradation pattern.

Figure 3: The difference in magnitudes for $b^{\prime}_{1}$ and $b_{1}$ for
parameter combinations of degree and $h$. Red indicates parameter combinations
where $b_{1}$, coming from the homophilous model, is smaller, and blue
indicates parameter combinations where $b^{\prime}_{1}$, coming from the
heterophilous model, is smaller. Smaller coefficients are preferred to
minimize the discrepancies on predictions induced by $\alpha_{t}$.
Predominantly negative values indicates that the heterophilous model
formulated as $\mathbf{(X\mathbin{\|}AX)W}$ is less sensitive to $\alpha_{t}$
and thus produces smaller discrepancies.

The impact of $b_{1}$ and $b^{\prime}_{1}$. We analyze how the different
coefficients influence the degradation of performance. To do this, we compare
the difference in magnitude between $b^{\prime}_{1}$ and $b_{1}$, where a
positive difference indicates a smaller magnitude for $b_{1}$ (indicating that
$b^{\prime}_{1}$ has a stronger influence on $\Delta\mathbf{z_{t}}$ as
compared to $b_{1}$), while a negative difference indicates a smaller
magnitude for $b^{\prime}_{1}$ (indicating that $b_{1}$ has a stronger
influence on $\Delta\mathbf{z_{t}}$ as compared to $b^{\prime}_{1}$). As the
goal is to identify which GNN formulation is likely to produce discrepancies,
we prefer models with smaller magnitudes given this will minimize the
influence of $\alpha_{t}$.

The comparison of $b^{\prime}_{1}$ and $b_{1}$ is performed in Figure 3, where
we sample points for the degree from $\\{1,5,10,15,20,25\\}$ and $h$ from $0$
to $1$ in intervals of $0.001$. We can see that the difference in the
coefficients is predominantly negative, indicating that the GNN $F^{\prime}$,
formulated as $\mathbf{(X\mathbin{\|}AX)W}$, is less impacted by shifts in
local homophily. Thus, although both the homophilous and heterophilous
formulations are impacted by discrepancies induced by $\alpha_{t}$, we can
expect the predictions of $F^{\prime}$ to experience less discrepancy as the
magnitude of $\alpha_{t}$ becomes larger.

####  A.2.3 The Impact of Degree on $\Delta\mathbf{z_{t}}$ on Homophilous
Models ($\mathbf{AXW}$)

As in the main text, we provide an outline of the cases below, with commentary
on the impact of degree $d$ on the cases.

Setting 1: Heterophilous $(\mathbf{0\leq\textit{h}<0.5}$): In this scenario,
when $d(2h-1)<-1$, $b_{1}<0$, leading to
$\mathbf{z_{t}}=\mathbf{y_{t,0}}+|b_{1}|\begin{bmatrix}-\alpha_{t}&\alpha_{t}\end{bmatrix}$
where the sign of $b_{1}$ has been distributed into the vector. There are then
two cases:

$\mathbf{z_{t}}=\begin{cases}\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}|\alpha_{t}|&-|\alpha_{t}|\end{bmatrix},&\text{if}\
h_{t}\leq h\\\
\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}-\alpha_{t}&\alpha_{t}\end{bmatrix},&\text{if}\
h_{t}>h.\end{cases}$ (26)

where the sign of $\alpha_{t}$ has been integrated into the vectors. To better
understand the conditions that lead to this scenario, we note that $d$ must be
greater than $\dfrac{1}{1-2h}$ to lead to $b_{1}<0$. However, even for
moderate values of $h$, the required degree to elicit a change is small, thus
the heterophily shift is often what will dominate the prediction, not the
degree. Furthermore, as $d(2h-1)$ becomes more negative, $b_{1}$ will tend
towards $\dfrac{1}{1+(2h-1)}$, demonstrating that the influence of $d$ on
predictions decreases as $d$ gets larger. As such, we can expect changes in
the local homophily to be the dominant factor in prediction changes.

Setting 2: Mixed Homophily ($\mathbf{\textit{h}=0.5}$): When the graph is not
strongly homophilous nor strongly heterophilous, $h=0.5\rightarrow b_{1}=d$,
leading to the following changes for $\mathbf{z_{t}}$:

$\mathbf{z_{t}}=\begin{cases}\mathbf{y_{t}}+d\begin{bmatrix}-|\alpha_{t}|&|\alpha_{t}|\end{bmatrix},&\text{if}\
h_{t}\leq h\\\
\mathbf{y_{t}}+d\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix},&\text{if}\
h_{t}>h.\end{cases}$ (27)

As the exception from other cases, the change in performance directly depends
on the degree $d$ and thus performance variations may be dominated by degree,
rather than local homophily, as $d$ becomes large. Intuitively, if the global
homophily level is 0.5, we can expect the model to not make use of this
signal.

Setting 3: Homophilous ($\mathbf{0.5<\textit{h}\leq 1}$): In this scenario,
$b_{1}>0$, thus:

$\mathbf{z_{t}}=\begin{cases}\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}-|\alpha_{t}|&|\alpha_{t}|\end{bmatrix},&\text{if}\
h_{t}\leq h\\\
\mathbf{y_{t}}+|b_{1}|\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix},&\text{if}\
h_{t}>h.\end{cases}$ (28)

Similar to Case 1, as $d(2h-1)$ becomes more positive, $b_{1}$ will approach
$\dfrac{1}{1+(2h-1)}$, leading to $d$ becoming less influential on the
prediction. As such, we can expect changes in the local homophily to be the
dominant factor in prediction changes.

High-level Takeaway: Notably, the heterophilous and homophilous settings are
not necessarily symmetric – the heterophilous setting requires the degree
criteria to be satisfied in order to attain the predictive pattern, while the
homophilous setting is always true. We suspect this plays a role in our (as
well as previous) analysis as our theoretical analysis suggests homophilous
datasets will always degrade in performance as test nodes become
heterophilous, while heterophilous datasets will only degrade when the degree
criteria is met. This phenomenon is seen in Figure 2 of our real-world
analysis, where the performance disparities are generally more drastic in the
homophilous datasets.

####  A.2.4 Additional Theoretical Analysis on 2-layer GNN

We follow a similar set up as in the 1-layer case for the graph $G$, but now
consider a 2-layer GNN, $F$, formulated as $\mathbf{(A+I)^{2}XW}$. During the
first aggregation step of $F$, a node $i$ of class $0$ will have updated
features as specified in Equation 6, and a node $i$ of class $1$ will have
updated features as specified in Equation 7.

During the second aggregation step of $F$, a node $i$ of class 0 will have
updated features $\mathbf{r_{i}^{2}}$ defined as:

$\displaystyle
y_{i}=0:\mathbf{r_{i}^{2}}=\dfrac{hd+1}{2}\begin{bmatrix}(1+d+p+dp(2h-1))&(1-p+d(1+p-2hp))\end{bmatrix}+$
(29)
$\displaystyle\dfrac{(1-h)d}{2}\begin{bmatrix}(1-p+d(1+p-2hp))&(1+d+p+dp(2h-1))\end{bmatrix}$

similarly, a node $i$ of class 1 will have updated features:

$\displaystyle
y_{i}=1:\mathbf{r_{i}^{2}}=\dfrac{hd+1}{2}\begin{bmatrix}(1-p+d(1+p-2hp))&(1+d+p+dp(2h-1))\end{bmatrix}+$
(30)
$\displaystyle\dfrac{(1-h)d}{2}\begin{bmatrix}(1+d+p+dp(2h-1))&(1-p+d(1+p-2hp))\end{bmatrix}$

Each $\mathbf{r_{i}^{2}}$ is the sum of the updated node vector for node $i$
after one aggregation step ($\mathbf{r_{i}^{1}}$), the updated feature vectors
from the homophilous connections of node $i$, and the updated feature vectors
of the heterophilous connections for node $i$. Condensing these terms, the
simplified expressions for each $\mathbf{r_{i}^{2}}$ is shown in Equations 31
and 32

$\displaystyle
y_{i}=0:\mathbf{r_{i}^{2}}=\dfrac{1}{2}\left[\begin{matrix}1+p+d^{2}(1+(1-2h)^{2}p)+d(2+(4h-2)p),\\\
\end{matrix}\right.$ (31)
$\displaystyle\qquad\qquad\left.\begin{matrix}1-p+d(2+(2-4h)p)-d^{2}(-1+(1-2h)^{2}p)\end{matrix}\right]$

$\displaystyle
y_{i}=1:\mathbf{r_{i}^{2}}=\dfrac{1}{2}\left[\begin{matrix}1-p+d(2+(2-4h)p)-d^{2}(-1+(1-2h)^{2}p),\\\
\end{matrix}\right.$ (32)
$\displaystyle\qquad\qquad\left.\begin{matrix}1+p+d^{2}(1+(1-2h)^{2}p)+d(2+(4h-2)p)\end{matrix}\right]$

Similar to the 1-layer case, $\mathbf{(A+I)^{2}XW}$ can be represented as RW,
where $\mathbf{R}$ is the matrix of updated features after two layers of
feature aggregation. $\mathbf{RW}$ is then set equal to the label matrix
$\mathbf{Y}$ after a transformation by the weight matrix W. Again, we sample
the unique data points, leaving us with the optimal W after solving the
reduced system of equations.

We now consider a new test point, $t$ with label $y_{t}=0$ and features
$\mathbf{x_{t}}=\dfrac{1}{2}\begin{bmatrix}(1+p)&(1-p)\end{bmatrix}$. In
addition, we assume that $t$ has a homophily ratio
$h+\alpha_{t}=h_{t}\in[0,1]$ where $\alpha_{t}$ represents a shift away from
$h$. When $\alpha_{t}$ is positive, we can interpret this as a test node with
a higher homophily ratio, and when $\alpha_{t}$ is negative, a lower homophily
ratio. Under this new homophily ratio, $\mathbf{r_{t}^{2}}$ can be calculated
as:

$\displaystyle\mathbf{r_{t}^{2}}=\dfrac{1}{2}\left[\begin{matrix}1+p+d^{2}(1+(1-2h_{t})^{2}p)+d(2+(4h_{t}-2)p),\\\
\end{matrix}\right.$ (33)
$\displaystyle\qquad\qquad\left.\begin{matrix}1-p+d(2+(2-4h_{t})p)-d^{2}(-1+(1-2h_{t})^{2}p)\end{matrix}\right]$

We can now compute $\mathbf{r_{t}^{2}}\textbf{W}=\mathbf{z_{t}^{2}}$, the
output of the 2-layer GNN, and analyze the associated predictions. The output
of $\mathbf{r_{t}^{2}}\textbf{W}$ is shown in Equation 34, where
$b_{2}=\dfrac{d}{(1+d(2h-1))^{2}}$

$\displaystyle\mathbf{z_{t}^{2}}=\dfrac{b_{2}}{d}\left[\begin{matrix}1+2d(h+h_{t}-1)+d^{2}(1-2h+2h^{2}-2h_{t}+2h_{t}^{2}),\\\
\end{matrix}\right.$ (34)
$\displaystyle\qquad\qquad\left.\begin{matrix}2d(h-h_{t})(1+d(-1+h+h_{t}))\end{matrix}\right]$

We now consider what happens when $h_{t}=h+\alpha_{t}$ where
$\alpha_{t}\in[-h,1-h]$. Since we are interested in understanding how the
predictions would change as a byproduct of $\alpha_{t}$, we look at the
difference between predictions when $h_{t}=h$. Below we compute
$\Delta\mathbf{z_{t}^{2}}$, the change in the logit vector output for the
2-layer GNN on the test point with local homophily $h_{t}$.

$\displaystyle\mathbf{z_{t}^{2}}=\dfrac{b_{2}}{d}\left[\begin{matrix}1+2d(h+(h+\alpha_{t})-1)+d^{2}(1-2h+2h^{2}-2(h+\alpha_{t})+2(h+\alpha_{t})^{2}),\\\
\end{matrix}\right.$ (35)
$\displaystyle\qquad\qquad\left.\begin{matrix}2d(h-(h+\alpha_{t}))(1+d(-1+h+(h+\alpha_{t})))\end{matrix}\right]$

$\displaystyle\mathbf{\Delta
z_{t}^{2}}=2b_{2}\begin{bmatrix}\alpha_{t}(1+d(2h-1+\alpha_{t}))&-\alpha_{t}(1+d(2h-1+\alpha_{t}))\end{bmatrix}$
(36) $\displaystyle\mathbf{\Delta
z_{t}^{2}}=2(1+d(2h-1+\alpha_{t}))b_{2}\begin{bmatrix}\alpha_{t}&-\alpha_{t}\end{bmatrix}\quad\square$

The change in the logit vector, $\mathbf{\Delta z_{t}^{2}}$, for the 2-layer
GNN, is similar to the change in the logit vector for the 1-layer GNN, except
with an additional multiplicative factor. Specifically, if we refer to the
change in the logit vector for the 1-layer GNN as $\mathbf{\Delta z_{t}^{1}}$,
$\mathbf{\Delta
z_{t}^{2}}=(2(1+d(2h-1+\alpha_{t}))\dfrac{b_{1}}{d}\mathbf{\Delta z_{t}^{1}}$.
We focus on analyzing this multiplicative factor, given we already understand
the behavior of $\mathbf{\Delta z_{t}^{1}}$ as studied in Section 4. Our
analysis is facilitated by letting the element-wise division of
$\mathbf{\Delta z_{t}^{2}}$ and $\mathbf{\Delta z_{t}^{1}}$ be equal to
$\mathbf{q}=\dfrac{\mathbf{\Delta z_{t}^{2}}}{\mathbf{\Delta
z_{t}^{1}}}=\dfrac{2b_{1}(1+d(2h-1+\alpha_{t}))}{d}\begin{bmatrix}1&1\end{bmatrix}=\dfrac{2(1+d(2h-1+\alpha_{t}))}{1+d(2h-1)}\begin{bmatrix}1&1\end{bmatrix}=\dfrac{2(1+d(2h-1)+d\alpha_{t})}{1+d(2h-1)}\begin{bmatrix}1&1\end{bmatrix}$
and studying the coefficient of $\mathbf{q}$. Specifically, we analyze two
extreme cases where $t$ is either part of a strongly heterophilous or strongly
homophilous graph to see what patterns emerge.

Setting 1: Strongly Heterophilous $(\mathbf{\textit{h}\approx 0}$): We study
this scenario by analyzing the limit of q as $h$ approaches 0, with $d>1$.

$\lim_{h\to 0}\textbf{q}=\dfrac{2(1+d(\alpha_{t}-1))}{1-d}$ (37)

In this scenario $0\leq a_{t}<1$. When $d(\alpha_{t}-1)<-1$, we see that the
coefficient is positive, thus following the same degradation trend as in the
1-layer GNN. Conversely, when $d(\alpha_{t}-1)>-1$ the coefficient is
negative, indicating that the performance degradation does not occur. However,
this requires that $\alpha_{t}>\dfrac{d-1}{d}$ which quickly becomes unlikely
to occur as $d$ grows, as even low values of $d$ require extremely large
$\alpha_{t}$ values to satisfy the inequality, e.g. $d=4$ requires
$\alpha_{t}>\dfrac{3}{4}$. We additionally analyze the magnitude of this
coefficient, finding that $\dfrac{2(1+d(\alpha_{t}-1))}{1-d}>1$ when
$\alpha_{t}<\dfrac{d-1}{2d}$, suggesting that the 2-layer GNN more rapidly
degrades in performance as compared to the 1-layer GNN when for moderate
values of $\alpha_{t}$, and behaves similarly to the 1-layer when $\alpha_{t}$
is large and $d$ is small.

Setting 2: Strongly Heterophilous $(\mathbf{\textit{h}\approx 1}$): We study
this scenario by analyzing the limit of q as $h$ approaches 1, with $d>1$.

$\lim_{h\to 1}\textbf{q}=\dfrac{2(1+d(\alpha_{t}+1))}{1+d}$ (38)

In this scenario $-1<a_{t}\leq 0$. As in the 1-layer case, the coefficient is
strictly positive, thus following the same degradation trend as in the 1-layer
GNN. Similar to the heterophilous case, we find that
$\dfrac{2(1+d(\alpha_{t}+1))}{1+d}>1$ when $\alpha_{t}>-\dfrac{d-1}{2d}$,
mirroring the behavior of the 2-layer GNN applied to heterophilous graphs with
increased degradation.

High-level Takeaway: When moving from a 1-layer to 2-layer GNN, we find that
the trend of performance degradation still holds, with performance degrading
as the local homophily ratio of a node deviates from the global homophily used
to train the GNN. Additionally, we show that the rate of performance
degradation for local homophily ratios far from the global homophily ratio is
generally increased when a 2-layer GNN is applied, relative to a 1-layer GNN.

###  A.3 Evaluation through Edge Homophily

In this section we discuss other edge-based homophily metrics proposed in the
literature. Additionally, we discuss how our results generalize to higher-
order neighborhoods.

Other Edge-based Homophily Metrics:  In this work, we focus on edge homophily
as a means of characterizing performance. Recent work has pointed out that the
global edge homophily metric, provided in Equation 1, can become susceptible
to imbalance in multi-class settings [27]. To remedy this issue, the authors
propose a class homophily metric which accounts for the size of each class
within the global homophily calculation. As our analysis focuses on local
homophily analysis, the issues in regards to the global homophily metric do
not directly apply when stratifying across local homophily ranges.
Additionally, the class homophily metric has no direct application to local
analysis, given it depends on graph-level properties. While it is possible
that the new class homophily metric can change the global homophily
characterization, we find that the heterophilous datasets in this work have
similar metrics as computed by both the class homophily and traditional global
edge homophily [27], suggesting either metric would lead to similar
conclusions. On heterophilous datasets with larger class imbalances, it may be
prudent to use the class homophily metric to better contextualize discrepancy.

Higher-Order Behavior:  Two natural question arise as to whether the results
hold when the GNN operates on higher-order information, as well as if the
results can be contextualized by higher-order information. We first note that
through the theoretical analysis in Appendix A.2.4, we can see that a GNN
operating on higher order information (2-hop neighborhoods) will experience
similar degradation patterns. Given the relationship under the 2-hop setting
is more complex, we prioritize the 1-hop theorem in the main text to provide a
more intuitive understanding of how performance degrades. Moreover, all of the
GNNs trained are cross validated on depths between 2 and 4, indicating that
empirically the higher order information again clearly trends with our
homophily measure. From the perspective of characterization, operating on
higher-order neighborhoods often follows similar patterns as the one-hop case
[27, 17], up to a certain point, as these metrics will approximate the global
homophily when the size of the neighborhoods approach the entire graph.
Moreover, there are practical challenges such as extremely large computation
times, and ambiguity in metric design, that make it difficult to apply higher-
order metrics. For instance, it is unclear if the metric should compare (a)
the ego-node label with the labels for nodes in the $k$-hop neighborhood, (b)
the labels between neighbors within the $k$-hop neighborhood, or (c) the ego-
node label as compared with the labels of nodes exactly $k$-hops from the ego-
node. In this work, we focus on the 1-hop case given its ubiquity across
research on homophily/heterophily, as well as its intuitive definition.

###  A.4 Synthetic Datasets

####  A.4.1 Graph Generation Related Work

Karimi et al. proposed the homophily-based preferential attachment model that
is used as the foundation for our synthetic generation method [12]. More
recently, GraphWorld has been proposed to allow for the study of GNN
performance across diverse graphs with varying global homophily, size, and
class imbalance [40]. Additionally, FastSGG has been proposed to generate
graphs for social network settings with an emphasis on control of the graph’s
degree distribution [41]. Despite these varying approaches, GraphWorld,
FastSGG, and Karimi et al’s preferential attachment model only consider
homophily through a fixed global parameter, providing no insight into, or
control of, the local homophily distribution in each graph. GenCAT remedies
this issue by introducing a node-class matrix to control the homophily level
for each node [42]. However, while this matrix provides the flexibility to
generate a graph with an arbitrary local homophily distribution, in practice
there are no mechanisms proposed to enforce a particular distribution over it.
Instead, each node’s homophily is assumed to be completely independent, which
is both unrealistic and cumbersome to instantiate. As such, systematically
controlling this distribution would require the introduction of a
parameterization over the node-class matrix, which is not proposed in the
work. Our approach solves this problem by including a new single parameter
into the preferential attachment model that enforces nodes to deviate from the
global homophily ratio, introducing nodes with a local homophily that has been
shifted relative to the global homophily ratio. This capability allows us to
easily learn over graphs with a wide range of local homophily ratios in a
highly controlled manner, demonstrating the practical implications of our
theory across different GNNs.

####  A.4.2 Synthetic Generation Algorithm - Additional Details

Initializing the Generated Graph: We follow the same strategy to initialize
the graph as is done in previous homophily-based preferential attachment
models, as well as in the classic Barabasi-Albert (BA) algorithm [43, 12, 7].
We begin by generating a node with label assigned from $P({0,...,c})$ and add
it to the empty graph $G$. We then sequentially add $m-1$ nodes to $G$, where
each node is first assigned a class label and then attached to one other node
in the graph. The attachment process for these nodes follows the GetNeighbors
function as outlined in the main text, with only one edge added, rather than
$m$ edges. The result of this process is a connected graph of $m$ nodes and
$m$ edges, allowing any subsequent nodes to have $m$ connections without
needing to create multiple edges between the same nodes.

Input: Total nodes $n$, $m$ edges to add per step, label probability
distribution $P(\\{0,...,c\\})$, uniformity probability $\rho$, class
compatibility matrix $\mathbf{H}_{L}$, drift change cutoff $\delta$

Initialize $G$ with $m$ nodes and $m$ edges according to $\mathbf{H}_{L}$ ;

// Details in App. A.4.2

1 Initialize vector $drift$ to hold change in homophily per node
$drift_{0},...,drift_{m}=0$

2 for  _$u=m$ to $n$_ do

3 Add node $u$ to $G$

4 Sample label $y_{u}\sim P(\\{0,...,c\\})$

5 $N_{1}(u)$ = GetNeighbors($G$, $m$, $u$, $drift$, $\rho$, $\mathbf{H}_{L}$,
$\delta$)

6 for  _$v$ in $N_{1}(u)$_ do

7 Add edge $(u,v)$ to $G$

8 if  _$y_{u}$ == $y_{v}$_ then

$drift[v]$ += 1 ;

// Drifted more homophilous

9

10 else

$drift[v]$ -= 1 ;

// Drifted more heterophilous

11

12

13 $drift[u]=0$

Result: Graph $G$

Algorithm 1 Synthetic Graph Generation

Input: $G$, $m$ edges to add, node $u$, $drift$ vector, uniformity probability
$\rho$, class compatibility matrix $\mathbf{H}_{L}$, drift change cutoff
$\delta$

1

2if  _$q\sim U(0,1)\leq\rho$_ then // Connect based on random local homophily
$h^{\prime}$

3 Sample new local homophily ratio $h^{\prime}\sim U(0,1)$

$hom\\_drift=\\{q:drift_{q}>\delta\\}$ ;

// Get nodes with large homophilous drift

$het\\_drift=\\{s:drift_{s}<-\delta\\}$ ;

// Get nodes with large heterophilous drift

4 $N_{1}(u)$ = {Sample up to $round(mh^{\prime})$ nodes with label $y_{u}$
from $het\\_drift$}

5 if  _$|N_{1}(u)|$ < $round(mh^{\prime})$_ then

6 $N_{1}(u)\cup$ {Sample $round(mh^{\prime})$ \- $|(N_{1}(u))|$ nodes with
label $y_{u}$ from $G\setminus het\\_drift$}

7 $N_{1}(u)\cup$ {Sample up to $round(m(1-h^{\prime}))$ nodes with label $\neq
y_{u}$ from $hom\\_drift$}

8 if  _$|N_{1}(u)|$ < $m$_ then

9 $N_{1}(u)\cup$ {Sample $m-|N_{1}(u)|$ nodes with label $\neq y_{u}$ from
$G\setminus hom\\_drift$ }

10

11else // Connect based on original compatibility matrix $\mathbf{H}_{L}$

12 $N_{1}(u)$ = {Sample $m$ nodes with probability based on $\mathbf{H_{L}}$}

Result: $N_{1}(u)$

Algorithm 2 GetNeighbors

####  A.4.3 Synthetic Dataset Properties

The synthetic dataset generation process is described in Section 5 of the main
text. Here, we detail some properties of the synthetic data, including the
distribution of local homophily ratios as the uniformity parameter is varied
and degree information.

Figure 4: Distribution of local homophily ratios generated as the global class
homophily is varied (rows) and the uniformity $\rho$ is varied (columns). We
see the nodes spread out to occupy the full range of homophily values around
$\rho=0.5$ while still maintaining a clear peak. Figure 5: Distribution of
degrees generated as the global class homophily is varied (rows) and the
uniformity $\rho$ is varied (columns). We see that $\rho$ does not
significantly impact the degree distribution other than for a small percentage
of nodes.

In Figure 4, we show how the distribution of local homophily changes as a
function of the uniformity parameter. On the left-most plot for each row, we
show the distribution of local homophily ratios on the unmodified preferential
attachment model. Despite the minor deviation from the set global homophily
ratio, the values are highly concentrated and do not provide samples across
the entire range of local homophily values. This is similarly seen for the
plot second to the left, which denotes $\rho=0.25$, where the extreme
homophily and heterophily values do not have enough nodes for analysis. The
plot second from the right, $\rho=0.5$, demonstrates the first scenario where
there are ample data points across the the full range of local homophily and
thus what we focus on in our analysis. Finally, the plot farthest to the
right, $\rho=0.75$, demonstrates a scenario where the distribution approaches
uniformity. As we are interested in the problem of distribution shift at test
time, this scenario is of less interest.

In Figure 5, we show the degree distribution for each of the graphs of varying
uniformity level. As compared to the unmodified preferential attachment model,
the uniformity parameter introduces a small amount of high degree nodes seen
to the right of the distribution in each plot. This is due to the correction
phase of the generation process, as described in the main text. However, as
the uniformity parameter increases, the amount of high degree nodes becomes
less severe. Moreover, even in the case of $h=0.5,\rho=0.25$, the high degree
nodes only account for around $1\%$ of the nodes with the other graphs only
containing even less.

####  A.4.4 Additional Synthetic Results

In this section, we include a series of additional experiments to supplement
the findings in the main text. For our synthetic data, we include results for
$h\in\\{0.3,0.7\\}$, as well as additional models across all synthetic graph
setups.

Pair-wise Differences between Homophilous and Heterophilous Models. In this
section we re-plot Figure 1, directly subtracting the performance of the
homophilous models from the heterophilous models. We show the difference in
performance for all comparisons in Figure 6, where the homophilous model is
subtracted from the heterophilous model. Thus, positive values are indicative
of the homophilous model performing better, while negative values are
indicative of the heterophilous model performing better. We show this for our
two extreme settings, $h=0.1,0.9$. These plots visually reinforce our findings
as the local homophily ranges near the global homophily ratio are always
positive (homophilous models performing better), while the local homophily
ranges far from the global homophily ratio are always negative (heterophilous
models performing better). Additionally, the relative performance differences
between the homophilous and heterophilous models in the positive and negative
sections of Figure 6 are significant – heterophilous models achieve upwards of
0.45 F1 improvement in under-represented local homophily ranges while often
losing less than 0.05 F1 on the over-represented local homophily ranges.

Figure 6: Performance differences between homophilous and heterophilous models
(columns) on synthetic data generated with global homophily ratios (columns)
$h\in\\{0.1,0.9\\}$ and uniformity $\rho=0.5$. Each bar represents the average
difference in F1 score between the two models specified in the column header
for nodes with a local homophily ratio between the values specified on the
x-axis.

Additional Experiments for $\mathbf{\textit{h}\in\\{0.3,0.7\\}}$. In this
section we highlight the performance of GCN, GAT, SAGE, GCNII, H2GCN, LINKX,
and MLP models on synthetic settings where $h\in\\{0.3,0.7\\}$. The results
are shown in 7. While the discrepancy is not as significant as before, this is
not unexpected as the theory indicates models will experience the most
degradation at extreme local homophily shifts. However, almost every model,
except for GCNII, degrades with a drop in F1 around 0.2 as the local homophily
level shifts relative to the global homophily. Regarding model design, GCNII
maintains a strong performance across all ranges in both settings, as does
H2GCN in the $h=0.7$ setting, showing the importance of heterophilous design
to minimize disparity.

Figure 7: Performance across models (columns) on synthetic data generated with
global homophily ratios (rows) $h\in\\{0.3,0.7\\}$ and uniformity $\rho=0.5$.
Each bar represents the F1 score for nodes with a local homophily ratio
between the values specified on the x-axis. Error bars denote the standard
deviation in performance. Blue bars indicate models with homophilous designs,
red bar indicates models with heterophilous designs.

Additional Models for All $h$. In addition to the models used in the main
text, i.e. GCN, GAT, SAGE, GCNII, and H2GCN, we include results on SGC, FAGCN
and GPRGNN. SGC is a linearized GNN architecture with a simple weighted
aggregation tailored towards homophily settings. FAGCN and GPRGNN are both
architectures designed to improve learning under heterophilous graphs, with
the ability to differentiate homophilic and heterophilic nodes through gating
mechanisms and adaptable weights, respectively. In Figure 8, we show the
results for these three models on each set of synthetic graph parameters.
Similar to the main text, we find that SGC significantly degrades in
performance as the local homophily ratio deviates from the global homophily
ratio, even more so than GCN or GAT. Additionally, we see that FAGCN and
GPRGNN both retain relatively strong performance, albeit slightly less than
H2GCN and GCNII. We find that FAGCN and GPRGNN are much more sensitive to
hyperparameters, as compared to H2GCN and GCNII, where the wrong choice can
cause their performance to devolve close to that of GCN. While the general
trend of heterophilous GNNs performing better on the full range of local
homophily values remains true, H2GCN’s concatenation based mechanisms provides
significant practical value due to not relying on hyperparameters for strong
performance.

Figure 8: Performance across models (rows) on synthetic data generated with
global homophily ratios (columns) $h\in\\{0.1,0.3,0.5,0.7,0.9\\}$ and
uniformity $\rho=0.5$. Each bar represents the F1 score for nodes with a local
homophily ratio between the values specified on the x-axis. Error bars denote
the standard deviation in performance. Blue bars indicate models with
homophilous designs, red bar indicates models with heterophilous designs.
Figure 9: Performance across models (columns) on synthetic data generated with
global homophily ratios (rows) $h\in\\{0.1,0.9\\}$ and uniformity $\rho=0.75$.
Each bar represents the F1 score for nodes with a local homophily ratio
between the values specified on the x-axis. Error bars denote the standard
deviation in performance. Blue bars indicate models with homophilous designs,
red bar indicates models with heterophilous designs. By introducing more
support across the range of local homophily ratios, each GNN is able to
achieve much lower discrepancy across groups.

Additional Performance for $\rho=0.75$. To further explain our choice for
$\rho=0.5$ in Section 9, we train all of our models (GNNs and MLP) on
additional synthetic datasets where $\rho=0.75$. The distribution of local
homophily ratios for $\rho=0.75$ can be seen in Figure 4. We do not consider
$\rho=0.25$ as there are not enough data points across the full local
homophily range to achieve a reasonable signal, as seen in Figure 4. For our
experiments, we follow the same setup as in Section 9 by training each model
on 10 different graphs of a particular parameter combination and analyze the
performance trends across groups of nodes with varying local homophily level.
We focus on $h=0.1$ and $h=0.9$ to maximize the possible shifts in local
homophily, relative to $h$. The results are displayed in Figure 9 for all
models. As expected, the performance discrepancy is significantly less than in
Figure 1, across all model. Moreover, the difference in performance between
the homophilous and heterophilous GNNs is less significant. As each GNN now
has ample support across the full range local homophily range during training,
each model is able to generalize, irrespective of the homophily level. While
$\rho=0.5$ is a useful tool to demonstrate how discrepancy can arise in GNN
models, Figure 9 is valuable in showing a set of conditions that can allow
GNNs to perform well, irrespective of design.

Additional Discussion and Experiments on Training Set Size. For our synthetic
analysis, we use a training/validation/test split ratio of 50/25/25%. Based on
previous works for GNNs applied to heterophilous settings, most use similar
(or even higher) ratio to study their models [7, 44, 45]. While studies have
shown that GNNs applied to homophilous settings can attain strong performance
with low training node ratios [20], it is common for heterophilous graphs to
require more training data as the local homophily patterns are more complex.
Despite this known behavior, the impact of training size on disparity is not
well studied. To address this gap, we include a set of experiments on our
synthetic data where the amount of training data is varied. All of the
parameters, aside from the amount of training data, is kept the same as
described in Section of the main text. As the main text uses 50% of the nodes
for training, we study the scenarios of 30%, and 10%.

Through the experiments depicted in Figure 10, we find that the amount of
training data does not significantly alter the discrepancy patterns. This
further verifies the findings of the main text that discrepancy is a
fundamental issue embedded within the models and training distribution, rather
than an artifact of training set size. We do note that many of the models have
amplified degradation as the homophily level deviates from the global
homophily ratio, as compared to the original results with a 50% training set
ratio. For example, GCN, GAT, and H2GCN nearly drop to an F1 score of 0 for
the furthest bins from the global homophily ratio as the training set ratio
drops. This is not unexpected, as the smaller training sets can cause less
nodes to fall within these bins during training. Thus, training set size can
influence the distribution of local homophily levels, which impacts the
discrepancy patterns, but it is not a direct issue. Some models retain
comparable, or even slightly better, performance irrespective of training set
size, however this is likely due to variations in training and are often
accompanied by larger variance bars.

Figure 10: Performance across models (columns) on synthetic data generated
with global homophily ratios (rows) $h\in\\{0.1,0.9\\}$ and uniformity
$\rho=0.5$ with training set ratios of 50% (top, the original results), 30%
(middle) and 10% (bottom). Performance degradation patterns generally remain
the same, highlighting discrepancy as a fundamental model and data
distribution problem, rather than a direct artifact of training set size.

###  A.5 Real-world Datasets

####  A.5.1 Additional Models for Real-world Datasets

Figure 11: F1 performance difference between GNNs and MLP ($\Delta$F1) for
varying local homophily ratios in real-world graphs. For a homophily range
specified on the x-axis, the 5 bars denote $\Delta$F1 for each GNN across test
nodes that fall within the range (ranges adjusted per dataset such that each
range has at least three nodes). Gray regions indicate the range that the
global homophily ratio falls in. More negative bars indicate worse
performance, while more positive bars indicate better performance. Error bars
are the standard deviation in performance across experiment runs. Blue bars
denote GNNs with homophilous designs, while red bars denote GNNs with
heterophilous designs. All results align with previous findings where best
performance is generally within the gray regions, and systematic performance
disparity emerges as the ranges move further from the global homophily.

Similar to our expanded analysis on the synthetic data, we include SGC, FAGCN,
and GPRGNN results for the real-world datasets from the main text. As seen in
Figure 11, the performance trends across the varying local homophily ranges
generally stay consistent with the findings in the main text. In particular,
we see that for Wisconsion, Arxiv-Year, Cora, and Coauthor-CS, the highest
performing subgroups for each model is in, or near, the gray region,
indicating that the global homophily ratio is influencing the performance.
Squirrel demonstrates a slightly different behavior by not degrading at the
furthest group, however does follow the trend across the first three
groupings. In the main text, we note that there is noise in this regime due to
the lack of nodes. Additionally, previous work has pointed out that the
squirrel dataset can become easier to predict on due to duplicates across the
structure [35]. As such, it may be possible that the data leakage caused by
these duplicates can influence performance, making certain homophily ranges
easier.

####  A.5.2 Local Homophily Distribution for Real-world Datasets

For each the dataset analyzed in this work, we include a histogram of their
local homophily distribution. While Cornell and Wisconsin are very similar and
predominately heterophilous, the other datasets provide a wide spread of nodes
across the varying local homophily ratios.

Figure 12: Histograms depicting the distribution of local homophily ratios for
the real-world datasets analyzed in the main text. Notably, all datasets,
despite their clear peak, possess nodes across the entire range of possible
local homophily ratios

###  A.6 Models and Hyperparameter Tuning

In this section we detail the models and hyperparameters tuned for each
experiment.

For synthetic experiments:

  1. 1.

MLP: Implementation directly from PyTorch

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Dropout: 0.3, 0.5

  2. 2.

GCN: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.gcn_conv)

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Dropout: 0.3, 0.5

  3. 3.

GAT: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.gat_conv)

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Heads: 1, 2

     * •

Dropout: 0.3, 0.5

  4. 4.

SAGE: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.sage_conv)

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Dropout: 0.3, 0.5

  5. 5.

GCN-II: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.gcn2_conv)

     * •

Hidden Dim: 16, 32

     * •

Depth: 4, 8

     * •

Dropout: 0.3, 0.5

  6. 6.

H2GCN: Open source PyTorch implementation from
https://github.com/CUAI/Non-Homophily-Large-Scale

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Dropout: 0.3, 0.5

  7. 7.

SGC: Implementation directly from PyTorch Geometric
(torch_geometric/nn/conv/sg_conv)

     * •

Depth: 2, 3

  8. 8.

GPR-GNN: Open source PyTorch implementation from
https://github.com/jianhao2016/GPRGNN

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Alpha: 0.1, 0.5, 0.9

     * •

K: 10

     * •

Dropout: 0.3, 0.5

  9. 9.

FA-GNN: Open source PyTorch implementation from
https://github.com/bdy9527/FAGCN

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Epsilon: 0.1, 0.5, 0.9

     * •

Dropout: 0.3, 0.5

  10. 10.

LINKX: Implementation directly from PyTorch Geometric
(torch_geometric/nn/models/LINKX)

     * •

Hidden Dim: 16, 32

     * •

Depth: 2, 3

     * •

Dropout: 0.3, 0.5

For real-world experiments:

  1. 1.

MLP: Implementation directly from PyTorch

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Dropout: 0.3, 0.5

  2. 2.

GCN: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.gcn_conv)

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Dropout: 0.3, 0.5

  3. 3.

GAT: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.gat_conv)

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Heads: 1, 2

     * •

Dropout: 0.3, 0.5

  4. 4.

SAGE: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.sage_conv)

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Dropout: 0.3, 0.5

  5. 5.

GCN-II: Implementation directly from PyTorch Geometric
(torch_geometric.nn.conv.gcn2_conv)

     * •

Hidden Dim: 32, 64

     * •

Depth: 4, 8, 16

     * •

Dropout: 0.3, 0.5

  6. 6.

H2GCN: Open source PyTorch implementation from
https://github.com/CUAI/Non-Homophily-Large-Scale

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Dropout: 0.3, 0.5

  7. 7.

SGC: Implementation directly from PyTorch Geometric
(torch_geometric/nn/conv/sg_conv)

     * •

Depth: 2, 3, 4

  8. 8.

GPR-GNN: Open source PyTorch implementation from
https://github.com/jianhao2016/GPRGNN

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Alpha: 0.1, 0.5, 0.9

     * •

K: 10

     * •

Dropout: 0.3, 0.5

  9. 9.

FA-GNN: Open source PyTorch implementation from
https://github.com/bdy9527/FAGCN

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Epsilon: 0.1, 0.5, 0.9

     * •

Dropout: 0.3, 0.5

  10. 10.

LINKX: Implementation directly from PyTorch Geometric
(torch_geometric/nn/models/LINKX)

     * •

Hidden Dim: 32, 64

     * •

Depth: 2, 3, 4

     * •

Dropout: 0.3, 0.5

All models have a final linear layer at the end of the convolution section to
produce the final predictions. Any unspecified parameters, such as the
additional parameters introduced by GCN-II, H2GCN, and LINKX are left as the
defaults in their respective code bases.