# Commute Graph Neural Networks

Wei Zhuo
Shenzhen Campus of Sun Yat-sen University
<EMAIL_ADDRESS>
&Guang Tan
Shenzhen Campus of Sun Yat-sen University
<EMAIL_ADDRESS>
###### Abstract

Graph Neural Networks (GNNs) have shown remarkable success in learning from
graph-structured data. However, their application to directed graphs
(digraphs) presents unique challenges, primarily due to the inherent asymmetry
in node relationships. Traditional GNNs are adept at capturing unidirectional
relations but fall short in encoding the mutual path dependencies between
nodes, such as asymmetrical shortest paths typically found in digraphs.
Recognizing this gap, we introduce Commute Graph Neural Networks (CGNN), an
approach that seamlessly integrates node-wise commute time into the message
passing scheme. The cornerstone of CGNN is an efficient method for computing
commute time using a newly formulated digraph Laplacian. Commute time
information is then integrated into the neighborhood aggregation process, with
neighbor contributions weighted according to their respective commute time to
the central node in each layer. It enables CGNN to directly capture the
mutual, asymmetric relationships in digraphs.

##  1 Introduction

Directed graphs (digraphs) are widely employed to model relational structures
in diverse domains, such as social networks (Cross et al., 2001) and
recommendation systems (Qiu et al., 2020). Recently, the advances of graph
neural networks (GNNs) have inspired various attempts to adopt GNNs for
analyzing digraphs (Tong et al., 2020a, b, 2021; Zhang et al., 2021; Rossi et
al., 2023; Geisler et al., 2023). The essence of GNN-based digraph analysis
lies in utilizing GNNs to learn expressive node representations that encode
edge direction information.

To achieve this, modern digraph neural networks are designed to integrate edge
direction information into the message passing process by distinguishing
between incoming and outgoing edges. This distinction enables the central node
to learn directionally discriminative information from its neighbors. As
illustrated in the digraph of Fig. 1, given a central node $v_{i}$, a 1-layer
digraph neural network can aggregate messages from $v_{i}$’s incoming neighbor
$v_{m}$ and outgoing neighbor $v_{j}$, and simultaneously capture edge
directions by applying direction-specific aggregation functions (Rossi et al.,
2023), or by predefining edge-specific weights (Zhang et al., 2021; Tong et
al., 2020b).

Figure 1: A digraph and its undirected counterpart. Blue arrows indicate
unidirectional paths, together with longer paths in the gray area, forming
commute closed loops between the central node $v_{i}$ and its outgoing
neighbors $v_{j}$ and $v_{k}$. In the undirected graph, shortest path
distances (SPD) between nodes are symmetric. However, in the digraph, the fact
that unidirectional SPDs are equal does not imply that mutual SPDs will also
be equal. For instance, while the SPD from $v_{i}$ to $v_{j}$ and $v_{k}$ is
identical, the reverse SPD from $v_{j}$ and $v_{k}$ back to $v_{i}$ do not
necessarily match this distance.

Despite the advancements, current digraph neural networks primarily capture
unidirectional111‘unidirectional’ refers to relationships in digraphs where
edges have a specific direction from one node to another.  relationships
between nodes, neglecting the complexity arising from path asymmetry. For
instance, a $k$-layer GNN aggregates the neighbors within the shortest path
$k$ for the central node. If the graph is undirected, the shortest path
between any two nodes is symmetric, as shown in the undirected graph of Fig.
1. This symmetry simplifies the representation of node relationships, implying
that if the SPDs from one node to two other nodes are identical, then the SPDs
from these two nodes back to the source node must also be the same.
Conversely, such symmetry is absent in digraphs. Considering the digraph in
Fig. 1, the shortest paths between $v_{i}$ and $v_{j}$ are asymmetric.
Therefore, although $v_{j}$ and $v_{k}$ are both immediate outgoing neighbors
of $v_{i}$, the strength of their relationships with the central node differs
significantly. Existing methods (Rossi et al., 2023; Tong et al., 2020b; Zhang
et al., 2021), by focusing solely on unidirectional shortest paths (blue and
red arrows), fail to capture the asymmetry phenomenon, which conveys valuable
information of node relationships. Take social networks as an example: an
ordinary user can directly follow a celebrity, yielding a short path to the
celebrity, yet the reverse path from the celebrity to the follower might be
much longer. Considering only the short path from the follower to celebrity
could falsely suggest a level of closeness that does not exist. In contrast,
accounting for the mutual paths between users yields a more precise and robust
measure of their relationship, with stronger mutual interactions implying
stronger connections.

To capture the mutual path interactions in GNNs, we adapt the concept of
commute time, the expected number of steps to traverse from a source node to a
target and back, from the Markov chain theory to the domain of graph learning.
To this end, we first generalize the graph Laplacian to the digraph by
defining the divergence of the gradient on the digraph. Utilizing this
digraph-specific Laplacian, we develop an efficient method to compute commute
time, ensuring sparsity and computational feasibility. Then we incorporate the
commute-time-based proximity measure into the message passing process by
assigning aggregation weights to neighbors. The intuition behind is that the
immediate and unidirectional neighboring relationships do not necessarily
imply strong similarity, but the mutual proximity is a more reliable indicator
of relationship closeness. Our experimental results demonstrate the
effectiveness of CGNN.

##  2 Related Work

###  2.1 Digraph Laplacian

While the Laplacian for undirected graphs has been extensively studied, the
area of Laplace operator digraphs remains underexplored. Chung (2005) pioneers
this area by defining a normalized Laplace operator specifically for strongly
connected directed graphs with nonnegative weights. This operator is expressed
as
$\mathbf{I}-\frac{\pi^{1/2}\mathbf{P}\pi^{-1/2}+\pi^{-1/2}\mathbf{P}^{*}\pi^{1/2}}{2}$.
Key to this formulation is the use of the transition probability operator
$\mathbf{P}$ and the Perron vector $\pi$, with the operator being self-
adjoint. Building on the undirected graph Laplacian, Singh et al. (2016) adapt
this concept to accommodate the directed structure, focusing particularly on
the in-degree matrix. They define the directed graph Laplacian as
$\mathbf{D}_{\text{in}}-\mathbf{A}$, where
$\mathbf{D}_{\text{in}}=\operatorname{diag}\left(\left\\{d_{i}^{\text{in}}\right\\}_{i=1}^{N}\right)$
represents the in-degree matrix. Li & Zhang (2012) uses stationary
probabilities of the Markov chain governing random walks on digraphs to define
the Laplacian as $\pi^{\frac{1}{2}}(\mathbf{I}-\mathbf{P})\pi^{-\frac{1}{2}}$,
which underscores the importance of random walks and their stationary
distributions in understanding digraph dynamics. Hermitian Laplacian Furutani
et al. (2020) consider the edge directionality and node connectivity
separately, and encode the edge direction into the argument in the complex
plane. Diverging from existing Laplacians, our proposed DiLap
$\mathbf{\Pi}(\mathbf{D}^{-1}-\mathbf{P})$ is grounded in graph signal
processing principles, conceptualized as the divergence of a signal’s gradient
on the digraph. It encompasses the degree matrix $\mathbf{D}$ to preserve
local connectivity, the transition matrix $\mathbf{P}$ to maintain the graph’s
directed structure, and the diagonalized Perron vector $\mathbf{\Pi}$,
capturing critical global graph attributes such as node structural importance,
global connectivity, and expected reachability (Chung, 1997).

###  2.2 Digraph Neural Networks

To effectively capture the directed structure with GNNs, spectral-based
methods (Zhang et al., 2021; Tong et al., 2020a, b) have been proposed to
preserve the underlying spectral properties of the digraph by performing
spectral analysis based on the digraph Laplacian proposed by (Chung, 2005).
MagNet (Zhang et al., 2021) utilizes magnetic Laplacian to derive a complex-
valued Hermitian matrix to encode the asymmetric nature of digraphs. Spatial
GNNs also offer a natural approach to capturing directed structures. For
instance, GraphSAGE (Hamilton et al., 2017) allows for controlling the
direction of information flow by considering in-neighbors or out-neighbors
separately. DirGNN (Rossi et al., 2023) further extends this framework by
segregating neighbor aggregation according to edge directions, offering a more
refined method to handle the directed nature of graphs.

##  3 Random Walk Distance and GNNs

We start by establishing the notations. We then show that message passing
based GNNs naturally capture the concept of hitting time during information
propagation across the graph, due to the unidirectional nature of the
neighborhood aggregation. Subsequently, we argue for the significance of
commute time, highlighting it as a more compact measure of mutual node-wise
interactions in random walks.

#### Notations

Consider $G=(V,E,\mathbf{X})$ as an unweighted digraph comprising $N$ nodes,
where $V=\\{v_{i}\\}^{N}_{i=1}$ is the node set, $E\subseteq(V\times V)$ is
the edge set, $\mathbf{X}\in\mathbb{R}^{N\times d}$ is the node feature
matrix. $Y=\\{y_{1},\cdots,y_{N}\\}$ is the set of labels for $V$. Let
$\mathbf{A}\in\mathbb{R}^{N\times N}$ be the adjacency matrix and
$\mathbf{D}=\mathrm{diag}(d_{1},\cdots,d_{N})\in\mathbb{R}^{N\times N}$ be the
degree matrix of $\mathbf{A}$, where $d_{i}=\sum_{v_{j}\in V}\mathbf{A}_{ij}$
is the out-degree of $v_{i}$. Let $\tilde{\mathbf{A}}=\mathbf{A}+\mathbf{I}$
and $\tilde{\mathbf{D}}=\mathbf{D}+\mathbf{I}$ denote the adjacency and degree
matrix with self-loops, respectively. The transition probability matrix of the
Markov chain associated with random walks on $G$ is defined as
$\mathbf{P}=\mathbf{D}^{-1}\mathbf{A}$, where
$\mathbf{P}_{ij}=\mathbf{A}_{ij}/\mathrm{deg}(v_{i})$ is the probability of a
1-step random walk starting from $v_{i}$ to $v_{j}$. Graph Laplacian
formalized as $\mathbf{L}=\mathbf{D}-\mathbf{A}$ is defined on the undirected
graph whose adjacency matrix is symmetric. The symmetrically normalized
Laplacian with self-loops (Wu et al., 2019) is defined as
$\hat{\mathbf{L}}=\tilde{\mathbf{D}}^{-\frac{1}{2}}\tilde{\mathbf{L}}\tilde{\mathbf{D}}^{-\frac{1}{2}}$,
where $\tilde{\mathbf{L}}=\tilde{\mathbf{D}}-\tilde{\mathbf{A}}$.

#### Digraph Neural Networks

DirGNN (Rossi et al., 2023) is a general framework that generalizes the
message passing paradigm to digraphs by adapting to the directionality of
edges. It involves separate aggregation processes for incoming and outgoing
neighbors of each node as follows:

$\displaystyle m_{i,\text{in}}^{(\ell)}$
$\displaystyle=\mathrm{Agg}_{\text{in}}^{(\ell)}\left(\left\\{h_{j}^{(\ell-1)}:v_{j}\in\mathcal{N}_{i}^{\text{in}}\right\\}\right)$
(1) $\displaystyle m_{i,\text{out}}^{(\ell)}$
$\displaystyle=\mathrm{Agg}_{\text{out}}^{(\ell)}\left(\left\\{h_{j}^{(\ell-1)}:v_{j}\in\mathcal{N}_{i}^{\text{out}}\right\\}\right)$
$\displaystyle h_{i}^{(l)}$
$\displaystyle=\mathrm{Comb}^{(\ell)}\left(h_{i}^{(\ell-1)},m_{i,\text{in}}^{(\ell)},m_{i,\text{out}}^{(\ell)}\right),$

where $\mathcal{N}_{i}^{\text{in}}$ and $\mathcal{N}_{i}^{\text{out}}$ are
respectively incoming and outgoing neighbors of $v_{i}$.
$\mathrm{Agg}_{\text{in}}^{(\ell)}(\cdot)$ and
$\mathrm{Agg}_{\text{out}}^{(\ell)}(\cdot)$ are specialized aggregation
functions of $\mathcal{N}_{i}^{\text{in}}$ and $\mathcal{N}_{i}^{\text{out}}$
at layer $\ell$, used to encode the directional characteristics of the edges
connected to $v_{i}$.

###  3.1 Can GNNs Capture Random Walk Distance?

In the context of random walks on a digraph, hitting time and commute time,
collectively referred to as random walk distances, serve as key metrics for
assessing node connectivity and interaction strength. Hitting time
$h(v_{i},v_{j})$ is the expected number of steps a random walk takes to reach
a specific target node $v_{j}$ for the first time, starting from a given
source node $v_{i}$. Commute time $c(v_{i},v_{j})$ is the expected number of
steps required for a random walk to start at $v_{i}$, reach $v_{j}$, and come
back. A high hitting (commute) time indicates difficulty in achieving
unidirectional (mutual) visits to each other in a random walk. As illustrated
in the digraph of Fig. 1, commute time $c(v_{i},v_{j})>c(v_{i},v_{k})$, while
the hitting time $h(v_{m},v_{i})=h(v_{i},v_{j})=h(v_{i},v_{k})$.

#### Motivation

Given these definitions, two questions arise: How crucial is it to retain
these measures in graph learning? Also, are message-passing GNNs capable of
preserving these characteristics? Firstly, both hitting time and commute time
are critical in understanding the structural dynamics of graphs. Hitting time,
analogous to the shortest path, measures the cost of reaching one node from
another, reflecting the directed influence or connectivity. Commute time,
encompassing the round-trip journey, offers insights into the mutual
relationships between nodes, which is especially evident in social networks,
as illustrated by celebrity-follower relationships. Secondly, message-passing
GNNs are somewhat effective in capturing hitting time, as they propagate
information across the graph in a manner similar to a random walk, where
quickly reached nodes are preferentially aggregated, and the influence of
nodes exponentially diminishes with increasing distance (Topping et al.,
2022). However, GNNs face challenges in preserving commute time due to their
requirement for comprehending mutual path relations, which are inherently
asymmetric and often involve longer-range interactions especially in digraphs,
which are not naturally captured in the basic message-passing framework.

Taking the digraph in Fig. 1 as an example, a 1-layer DirGNN defined in Eq. 1
can encode $v_{m}$, $v_{j}$ and $v_{k}$ into the representation of $v_{i}$,
while also capturing the directionality of edges from these neighbors by using
distinct aggregation functions for incoming and outgoing neighbors. It shows
that DirGNN can capture the hitting time, as neighbors with lower hitting
times, $h(v_{i},v_{k})$, $h(v_{i},v_{k})$ and $h(v_{m},v_{i})$, are aggregated
preferentially. However, DirGNN inherently focuses on unidirectional
interactions and overlooks mutual path dependencies. Notably, a 1-layer DirGNN
is insufficient for capturing the asymmetric interactions indicated by the
paths from $v_{j}$ and $v_{k}$ returning to the central node $v_{i}$ (gray
areas). One potential approach to address this limitation is to stack
additional message passing layers to encompass the entire commute path between
nodes, thereby capturing mutual path interactions. Nevertheless, this strategy
is non-trivial because the commute paths vary considerably across different
node pairs, complicating the determination of an appropriate number of layers.
Additionally, stacking multiple layers to cover these paths can introduce
irrelevant non-local information and lead to oversmoothing.

#### Goal

We expect to directly encode node-wise commute time into the node
representations to accurately reflect the true interaction strength between
adjacent nodes during neighbor aggregation, accounting for both forward and
backward paths. For instance, even though $h(v_{i},v_{j})=h(v_{i},v_{k})$, a
shorter commute time $c(v_{i},v_{k})<c(v_{i},v_{j})$ suggests a stronger
interaction from $v_{k}$ to $v_{i}$ compared to $v_{j}$ to $v_{i}$.
Consequently, the contribution of neighbor $v_{k}$ to the representation of
$v_{i}$ should be greater than that of $v_{j}$.

###  3.2 Commute Time Computation

Based on the standard Markov chain theory, a useful tool to study random walk
distances is the fundamental matrix (Aldous & Fill, 2002). We first establish
the following assumptions required to support the theorem.

######  Assumption 3.1.

The digraph $G$ is irreducible and aperiodic.

These two properties pertain to the Markov chain’s stationary probability
distribution $\pi$ (i.e., Perron vector) corresponding to the given graph.
Irreducibility ensures that it is possible to reach any node (state) from any
other node, preventing $\pi$ from converging to 0. Aperiodicity ensures that
the Markov chain does not get trapped in cycles of a fixed length, thus
guaranteeing the existence of a unique $\pi$. Existence and uniqueness of
$\pi$ facilitate deterministic analysis and computation. For a more intuitive
understanding of the assumptions, we give the sufficient conditions of digraph
under the irreducibility and aperiodicity assumptions.

######  Proposition 1.

A strongly connected digraph, in which a directed path exists between every
pair of vertices, is irreducible. A digraph with self-loops in each node is
aperiodic.

Given the above assumption, the fundamental matrix $\mathbf{Z}$ is defined as
the sum of an infinite matrix series:

$\mathbf{Z}=\sum_{t=0}^{\infty}\left(\mathbf{P}^{t}-\mathbf{J}\mathbf{\Pi}\right)=\sum_{t=0}^{\infty}\left(\mathbf{P}^{t}-e\pi^{\top}\right),$
(2)

where $e$ is the all-one column vector, then we have $\mathbf{J}=e\cdot
e^{\top}$ is the all-one matrix, and $\mathbf{\Pi}=\mathrm{diag}(\pi)$ is the
diagonal matrix of $\pi$.

######  Theorem 3.2.

(Li & Zhang, 2012) Given 3.1, the fundamental matrix $\mathbf{Z}$ defined in
Eq. 2 converges to:

$\mathbf{Z}=(\mathbf{I}-\mathbf{P}+\mathbf{J}\mathbf{\Pi})^{-1}-\mathbf{J}\mathbf{\Pi},$
(3)

where $\mathbf{I}$ is an identity matrix.

The hitting time and commute time on $G$ can then be expressed as $\mathbf{Z}$
(Aldous & Fill, 2002) as follows:

$h(v_{i},v_{j})=\frac{\mathbf{Z}_{jj}-\mathbf{Z}_{ij}}{\pi_{j}},\quad
c(v_{i},v_{j})=h(v_{i},v_{j})+h(v_{j},v_{i}).$ (4)

However, directly calculating the complete fundamental matrix $\mathbf{Z}$ and
the commute times for all node pairs is computationally expensive and yields a
dense matrix. Moreover, integrating the random walk distances computation,
defined in Eq. 3 and Eq. 4, into the message passing framework is non-trivial,
which concerns the scalability of the model.

##  4 Commute Graph Neural Networks

In this section, we present Commute Graph Neural Networks (CGNN) to encode the
commute time information into message passing. We first establish the
relationship between random walk distances and the digraph Laplacian.

###  4.1 Digraph Laplacian (DiLap)

Contrary to the traditional graph Laplacian, typically defined as a symmetric
positive semi-definite matrix derived from the symmetric adjacency matrix, our
proposed DiLap is built upon the transition matrix to preserve the directed
structure. Specifically, the classical graph Laplacian
$\mathbf{L}=\mathbf{D}-\mathbf{A}$ is interpreted as the divergence of the
gradient of a signal on an undirected graph (Shuman et al., 2013; Hamilton,
2020): given a graph signal $s\in\mathbb{R}^{N}$,
$(\mathbf{L}s)_{i}=\sum_{j\in\mathcal{N}_{i}}\mathbf{A}_{ij}(s_{i}-s_{j})$.
Intuitively, graph Laplacian corresponds to the difference operator on the
signal $s$, and acts as a node-wise measure of local smoothness. In line with
this conceptual foundation, we generalize the graph Laplacian to digraphs by
defining DiLap $\mathbf{T}$:

$\displaystyle(\mathbf{T}s)_{i}=\sum_{v_{j}\in\mathcal{N}^{\text{out}}_{i}}\mathbf{P}_{ij}(s_{i}-s_{j})=((\mathbf{D}^{-1}-\mathbf{P})s)_{i},\quad\text{we
have}\quad\mathbf{T}=\mathbf{D}^{-1}-\mathbf{P},$ (5)

where $\mathcal{N}^{\text{out}}_{i}$ is the set of $v_{i}$’s out-neighbors,
and $\mathbf{P}_{ij}=\frac{1}{d^{\text{out}}_{i}}$ with $d^{\text{out}}_{i}$
being the out-degree of $v_{i}$. The detailed derivation of $\mathbf{T}$ is
included in Section A.1. Eq. 5 illustrates that DiLap $\mathbf{T}$ acts as the
divergence of the gradient of a given signal on a digraph, and can be used to
measure the local smoothness in a digraph.

Considering the Laplacian operator’s function in measuring the smoothness of a
signal across the entire graph, it is crucial to assign greater weights to
nodes of higher structural importance. It is because the smoothness at nodes
central to the graph structure should have a more pronounced impact on the
global smoothness measurement. Thus, we further define the Weighted DiLap
$\mathcal{T}$:

$\displaystyle(\mathcal{T}s)_{i}=\sum_{v_{j}\in\mathcal{N}^{\text{out}}_{i}}\pi_{i}\mathbf{P}_{ij}(s_{i}-s_{j})=(\mathbf{\Pi}(\mathbf{D}^{-1}-\mathbf{P})s)_{i},\quad\text{we
have}\quad\mathcal{T}=\mathbf{\Pi}(\mathbf{D}^{-1}-\mathbf{P}).$ (6)

Here we utilize the $i$-th element of the Perron vector $\pi$ to quantify the
structural importance of $v_{i}$, reflecting its eigenvector centrality. This
is based on the principle that a node’s reachability is directly proportional
to its corresponding value in the Perron vector (Xu et al., 2018). Therefore,
$\pi$ effectively indicates the centrality and influence over the long term in
the graph. Perron-Frobenius Theorem (Horn & Johnson, 2012) establishes that
$\pi$ satisfies $\sum_{i}\pi(i)=1$, is strictly positive, and converges to the
left eigenvector of the dominant eigenvalue of $\mathbf{P}$.

###  4.2 Similarity-based Graph Rewiring

Both the fundamental matrix defined in Eq. 3 and Weighted DiLap requires 3.1
to ensure the existence and uniqueness of the Perron vector $\pi$, conditions
that are not universally met in general graphs. To fulfill the irreducibility
and aperiodicity assumptions, Tong et al. (2020a) introduce a teleporting
probability uniformly distributed across all nodes. This method, inspired by
PageRank (Page et al., 1999), amends the transition matrix to
$\mathbf{P}_{pr}=\gamma\mathbf{P}+(1-\gamma)\frac{ee^{\top}}{N}$, where
$\gamma\in(0,1)$. $\mathbf{P}_{pr}$ allows for the possibility that a random
walker may choose a non-neighbor node for the next step with a probability of
$\frac{1-\gamma}{N}$. This adjustment ensures that $\mathbf{P}_{pr}$ is
irreducible and aperiodic, so it has a unique $\pi$. However, this approach
leads to a complete graph represented by a dense matrix $\mathbf{P}_{pr}$,
posing significant challenges for subsequent computational processes.

Rather than employing $\mathbf{P}_{pr}$ as the transition matrix, we introduce
a graph rewiring method based on feature similarity to make a given graph
irreducible, while maintaining the sparsity. As outlined in Proposition 1, to
transform the digraph into a strongly connected structure, it is essential
that each node possesses a directed path to every other node. To this end, we
initially construct a simple and irreducible graph $G^{\prime}$ with all $N$
nodes, then add all edges from $G^{\prime}$ the original digraph $G$, thereby
ensuring $G$’s irreducibility. The construction of $G^{\prime}$ begins with
the calculation of the mean of node features as the anchor vector $\bm{a}$.
Then we determine the similarity between each node and the anchor, sort the
similarity values, and return the sorted node indices, denoted as
$S\in\mathbb{R}^{N}$:

$\bm{a}=\frac{\sum_{i}\mathbf{X}_{i}}{N},\quad
s_{i}=\cos(\bm{a},\mathbf{X}_{i}),\quad
S=\arg\mathrm{sort}(\\{s_{i}\\}_{i=1}^{N})$ (7)

where $\cos(\bm{a},\mathbf{X}_{i})$ is the cosine similarity between node
features of $v_{i}$ and $\bm{a}$, and $\mathrm{argsort}(\cdot)$ yields the
indices of nodes that sort similarity values $\\{s_{i}\\}_{i=1}^{N}$. We then
connect the nodes one by one with undirected (bidirectional) edges following
the order in $S$ to construct $G^{\prime}$, as shown in Fig. 2. Given that
$G^{\prime}$ is strongly connected, adding all its edges into $G$ results in a
strongly connected digraph $\widetilde{G}$, which is irreducible. To achieve
aperiodicity, self-loops are further added to $\widetilde{G}$.

Figure 2: The sorted node indices in $S$ are connected one by one with
undirected edges to construct $G^{\prime}$, then adding all edges from
$G^{\prime}$ to $G$ to generate $\widetilde{G}$.

This rewiring approach satisfies 3.1 and maintains graph sparsity.
Additionally, adding edges between nodes with similar features only minimally
alters the overall semantics of the original graph. Based on $\widetilde{G}$
and its corresponding $\widetilde{\mathbf{P}}$ and $\widetilde{\mathbf{\Pi}}$,
we have the deterministic Weighted DiLap $\widetilde{\mathcal{T}}$.

###  4.3 From DiLap to Commute Time

Given the Weighted DiLap $\widetilde{\mathcal{T}}$, we can unify the commute
time information into the message passing by building the connection between
$\widetilde{\mathcal{T}}$ and the fundamental matrix $\mathbf{Z}$:

######  Lemma 4.1.

Given a rewired graph $\widetilde{G}$, the Weighted DiLap is defined as
$\widetilde{\mathcal{T}}=\widetilde{\mathbf{\Pi}}(\widetilde{\mathbf{D}}^{-1}-\widetilde{\mathbf{P}})$.
Then the fundamental matrix $\mathbf{Z}$ of $\widetilde{G}$ can be solved by:

$\mathbf{Z}=\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}\mathcal{R}^{\dagger}\widetilde{\mathbf{\Pi}}^{\frac{1}{2}},$
(8)

where
$\mathcal{R}=\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}\widetilde{\mathcal{T}}\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I}$
and the superscript $\dagger$ means Moore–Penrose pseudoinverse of the matrix.

The proof is given in Section A.2. Leveraging Lemma 4.1, we can further
compute the hitting times and commute times in terms of
$\widetilde{\mathcal{T}}$ with the following theorem.

######  Theorem 4.2.

Given $\widetilde{G}$, the hitting time and commute time from $v_{i}$ to
$v_{j}$ on $\widetilde{G}$ can be computed as follows:

$\displaystyle h(v_{i},v_{j})$
$\displaystyle=\frac{\mathcal{R}^{\dagger}_{jj}}{\pi_{j}}-\frac{\mathcal{R}^{\dagger}_{ij}}{\sqrt{\pi_{i}\pi_{j}}},$
(9) $\displaystyle c(v_{i},v_{j})$
$\displaystyle=h(v_{i},v_{j})+h(v_{j},v_{i})=\frac{\mathcal{R}^{\dagger}_{jj}}{\pi_{j}}+\frac{\mathcal{R}^{\dagger}_{ii}}{\pi_{i}}-\frac{\mathcal{R}^{\dagger}_{ij}}{\sqrt{\pi_{i}\pi_{j}}}-\frac{\mathcal{R}^{\dagger}_{ji}}{\sqrt{\pi_{i}\pi_{j}}}.$

For a detailed proof of this theorem, refer to Theorem 3 in (Li & Zhang,
2012), where the corresponding term should be replaced with $\mathcal{R}$.
Then we can derive the matrix forms of the hitting time $\mathcal{H}$ and
commute time $\mathcal{C}$ as per Eq. 9:

$\displaystyle\mathcal{H}=(e\otimes\pi^{-1})(\mathcal{R}^{\dagger}\odot\mathbf{I})-\mathcal{R}^{\dagger}\odot(\pi^{-\frac{1}{2}}\otimes\pi^{-\frac{1}{2}}),\quad\mathcal{C}=\mathcal{H}+\mathcal{H}^{\top}$
(10)

where $\odot$ denotes Hadamard product, and $\otimes$ is outer product.
$\mathcal{H}_{ij}$ and $\mathcal{C}_{ij}$ correspond to the hitting and
commute time from $v_{i}$ to $v_{j}$ respectively. The computation of commute
times via DiLap, in contrast to the method delineated in Theorem 3.2, is
primarily motivated by efficiency concerns. Specifically, Eq. 3 necessitates
the inversion of a dense matrix with complexity $\mathcal{O}(N^{3})$, whereas
our DiLap-based method hinges on computing the pseudoinverse of a sparse
matrix $\mathcal{R}$. The pseudoinverse of $\mathcal{R}$ can be efficiently
determined using SVD. Given the sparse nature of $\mathcal{R}$, we can employ
well-established techniques such as the randomized truncated SVD algorithm
(Halko et al., 2011; Cai et al., 2023), which takes advantage of sparsity, to
reduce the time complexity to $\mathcal{O}(q|E|)$, where $|E|$ denotes the
number of edges reflecting the sparsity (See Section A.3). Next, we present
Commute Graph Neural Networks (CGNN) based on $\mathcal{C}$.

###  4.4 CGNN

$\mathcal{C}\in\mathbb{R}^{N\times N}$ quantifies the strength of mutual
relations between node pairs in the random walk context. Notably, smaller
values in $\mathcal{C}$ correspond to stronger mutual reachability, indicating
stronger relations between node pairs. Thus, $\mathcal{C}$ is a positive
symmetric matrix, and the commute-time-based node proximity matrix can be
expressed as $\widetilde{\mathcal{C}}=\exp(-\mathcal{C})$. Since the directed
adjacency matrix $\mathbf{A}$ represents the outgoing edges of each node,
$\mathbf{A}^{\top}$ therefore accounts for all incoming edges. Then we have
$\widetilde{\mathcal{C}}^{\text{out}}=\mathbf{A}\odot\widetilde{\mathcal{C}}$
and
$\widetilde{\mathcal{C}}^{\text{in}}=\mathbf{A}^{\top}\odot\widetilde{\mathcal{C}}$
represent the proximity between adjacent nodes under outgoing and incoming
edges, respectively. We further perform row-wise max-normalization on
$\widetilde{\mathcal{C}}^{\text{out}}$ and
$\widetilde{\mathcal{C}}^{\text{in}}$ to rescale the maximum value in each row
to 1. Given the original graph $G$ as input, we define the $\ell$-th layer of
CGNN as:

$\displaystyle m_{i,\text{in}}^{(\ell)}$
$\displaystyle=\mathrm{Agg}_{\text{in}}^{(\ell)}\left(\left\\{\widetilde{\mathcal{C}}^{\text{in}}_{ij}\cdot
h_{j}^{(\ell-1)}:v_{j}\in\mathcal{N}_{i}^{\text{in}}\right\\}\right)$ (11)
$\displaystyle m_{i,\text{out}}^{(\ell)}$
$\displaystyle=\mathrm{Agg}_{\text{out}}^{(\ell)}\left(\left\\{\widetilde{\mathcal{C}}^{\text{out}}_{ij}\cdot
h_{j}^{(\ell-1)}:v_{j}\in\mathcal{N}_{i}^{\text{out}}\right\\}\right)$
$\displaystyle h_{i}^{(l)}$
$\displaystyle=\mathrm{Comb}^{(\ell)}\left(h_{i}^{(\ell-1)},m_{i,\text{in}}^{(\ell)},m_{i,\text{out}}^{(\ell)}\right),$

where $\mathrm{Agg}_{\text{in}}^{(\ell)}(\cdot)$ and
$\mathrm{Agg}_{\text{out}}^{(\ell)}(\cdot)$ are mean aggregation functions
with different feature transformation weights, and
$\mathrm{Comb}^{(\ell)}(\cdot)$ is a mean operator. Within each layer, the
influence of $v_{j}$ on the central node $v_{i}$ is modulated by the commute-
time-based proximity $\widetilde{\mathcal{C}}$ based on the edge
directionality. We present the pseudocode and complexity of CGNN in Algorithm
1.

##  5 Experiments

We conduct node classification experiments on five digraph datasets.
Experimental details and data statistics are provided in Section C.1 and
Section C.2. We provide a performance comparison with 12 baselines including
1) General GNNs: GCN (Kipf & Welling, 2017), GAT (Veličković et al., 2018),
and GraphSAGE (Hamilton et al., 2017); 2) Non-local GNNs: APPNP (Klicpera et
al., 2019), MixHop (Abu-El-Haija et al., 2019), GPRGNN (Chien et al., 2021),
and GCNII (Ming Chen et al., 2020); 3) Digraph NNs: DGCN (Tong et al., 2020b),
DiGCN (Tong et al., 2020a), MagNet (Zhang et al., 2021), DiGCL (Tong et al.,
2021), and DirGNN (Rossi et al., 2023). For all baselines, we apply both the
symmetrized and asymmetric adjacency matrix for experiments. The results
reported are the better of the two results.

Table 1: Node classification results. Accuracy (%) with standard deviation for
10 runs. We highlight/underline the best/second best method. For general GNN
and non-local GNN baselines, we conduct experiments on both symmetrized
versions and their directed counterparts, reporting better results from these
two settings.

Method Squirrel Chameleon Citeseer CoraML AM-Photo GCN 52.43$\pm$2.01
67.96$\pm$1.82 66.03$\pm$1.88 70.92$\pm$0.39 88.52$\pm$0.47 GAT 40.72$\pm$1.55
60.69$\pm$1.95 65.58$\pm$1.39 72.22$\pm$0.57 88.36$\pm$1.25 GraphSAGE
41.61$\pm$0.74 62.01$\pm$1.06 66.81$\pm$1.38 74.16$\pm$1.55 89.71$\pm$0.57
APPNP 51.91$\pm$0.56 45.37$\pm$1.62 66.90$\pm$1.82 70.31$\pm$0.67
87.43$\pm$0.98 MixHop 43.80$\pm$1.48 60.50$\pm$2.53 56.09$\pm$2.08
65.89$\pm$1.50 87.17$\pm$1.34 GPRGNN 50.56$\pm$1.51 66.31$\pm$2.05
61.74$\pm$1.87 73.31$\pm$1.37 90.23$\pm$0.34 GCNII 38.47$\pm$1.58
63.86$\pm$3.04 58.32$\pm$1.93 64.84$\pm$0.71 83.40$\pm$0.79 DGCN
37.16$\pm$1.72 50.77$\pm$3.31 66.37$\pm$1.93 75.02$\pm$0.50 87.74$\pm$1.02
DiGCN 33.44$\pm$2.07 50.37$\pm$4.31 64.99$\pm$1.72 77.03$\pm$0.70
88.66$\pm$0.51 MagNet 39.01$\pm$1.93 58.22$\pm$2.87 65.04$\pm$0.47
76.32$\pm$0.10 86.80$\pm$0.65 DiGCL 35.82$\pm$1.73 56.45$\pm$2.77
67.42$\pm$0.14 77.53$\pm$0.14 89.41$\pm$0.11 DirGNN 75.19$\pm$1.26
79.11$\pm$2.28 66.57$\pm$0.74 75.33$\pm$0.32 88.09$\pm$0.46 CGNN
77.20$\pm$2.42 79.54$\pm$1.82 70.27$\pm$0.66 77.02$\pm$0.67 90.41$\pm$0.32

###  5.1 Overall Results

Table 1 reports the node classification results across five digraph datasets.
Our method CGNN achieves new state-of-the-art results on 4 out of 5 datasets,
and comparable results on the remaining one, validating the superiority of
CGNN. We provide more observations as follows. Firstly, while non-local GNNs
have the potential to cover the commute paths between adjacent nodes by
stacking multiple layers, they do not consistently outperform general, shallow
GNN models. It suggests that coarsely aggregating all nodes in commute paths
is ineffective. The reason is that deeper models may aggregate excessive
irrelevant information for the central node. Our goal is to encode mutual
relationships between adjacent nodes by considering their commute times.
Aggregating all nodes along the entire path introduces excessive information
about other nodes unrelated to the direct relationship between the target
nodes. Secondly, GNNs tailored for digraphs do not seem to bring substantial
gains. Our results show that with careful hyper-parameter tuning, general GNNs
can achieve results comparable to, or even better than, those tailored for
digraphs, as evidenced in the Squirrel, Chameleon, and AM-Photo datasets.
Thirdly, CGNN achieves state-of-the-art results on both homophilic and
heterophilic digraph benchmarks. Notably, DirGNN also performs comparably on
heterophilic graphs (Squirrel and Chameleon), supporting the findings of Rossi
et al. (2023) that distinguishing edge directionality during message passing
enables the central node to adaptively balance information flows from both
heterophilic and homophilic neighbors, effectively mitigating the impact of
heterophily. Moreover, CGNN, an enhanced version of DirGNN, further improves
performance on these graphs by effectively incorporating commute times to
refine the strength of relationships between nodes, enhancing model robustness
under heterophily.

To illustrate this, we further examine the relations between commute-time-
based proximity and label similarity along edges. As shown in Eq. 11, we use
commute-time-based proximity $\widetilde{\mathcal{C}}$ to weigh the neighbors
during neighbor aggregation. Then we define a label similarity matrix
$\mathcal{M}$ where $\mathcal{M}_{ij}=1$ if $v_{j}\in\mathcal{N}_{i}$ and
$y_{i}=y_{j}$; otherwise $\mathcal{M}_{ij}=0$. Essentially, $\mathcal{M}$
extracts the edges connecting nodes with the same classes from the adjacency
matrix $\mathbf{A}$. Thus a higher value of
$\|\mathcal{M}-(\mathbf{A}+\mathbf{A}^{\top})\|_{2}^{2}$ indicates a more
pronounced negative impact of heterophily on the model’s performance. On the
other hand, we compute
$\|\mathcal{M}-(\widetilde{\mathcal{C}}^{\text{in}}+\widetilde{\mathcal{C}}^{\text{out}})\|_{2}^{2}$
to evaluate the efficacy of $\widetilde{\mathcal{C}}$ in filtering
heterophilic information. The closer
$(\widetilde{\mathcal{C}}^{\text{in}}+\widetilde{\mathcal{C}}^{\text{out}})$
is to $\mathcal{M}$, the more effectively it aids the model in discarding
irrelevant heterophilic information. Table 2 visually demonstrates these
relationships. We observe that in heterophilic datasets, the commute-time-
based proximity matrix
$(\widetilde{\mathcal{C}}^{\text{in}}+\widetilde{\mathcal{C}}^{\text{out}})$,
aligns more closely with the label similarity matrix $\mathcal{M}$ than
$(\mathbf{A}+\mathbf{A}^{\top})$. It indicates that $\widetilde{\mathcal{C}}$
effectively filters out irrelevant information during message passing by
appropriately weighting neighbors, which explains the exceptional performance
of CGNN on heterophilic datasets.

Table 2: The distance between $\mathcal{M}$ and $\mathbf{A}$, and between
$\mathcal{M}$ and $\widetilde{\mathcal{C}}$.

(a) Squirrel

(b) AM-Photo

Figure 3: Accuracy vs. running time on Squirrel and AM-Photo.

###  5.2 Efficiency Comparsion

Fig. 3 compares the accuracy of CGNN and other baseline models with their
running times. Despite the additional computational load of calculating
commute-time-based proximity, the results show that CGNN provides the best
trade-off between effectiveness and efficiency. In particular, on the Squirrel
dataset, CGNN has the third-fastest calculation speed while yielding accuracy
nearly double that of all other methods. On AM-Photo, CGNN achieves the
highest accuracy while maintaining moderate efficiency.

###  5.3 Component Analysis

#### Comparison between graph rewiring and PPR.

In Section 4.2, we construct a rewired graph $\widetilde{G}$ based on feature
similarity to guarantee the irreducibility and aperiodicity. On the other
hand, In contrast, the classic PageRank transition matrix, defined as
$\mathbf{P}_{pr}=\gamma\mathbf{P}+(1-\gamma)\frac{ee^{\top}}{N}$, achieves a
similar objective but results in a completely connected graph $G_{pr}$.
However, this approach tends to overlook the sparse structure of the original
graph, which may alter the semantic information in the graph. Additionally,
computing commute times using a dense transition matrix incurs a high
computational cost. To validate the effectiveness of the rewiring approach
over the PPR method, we conduct an experiment where $\widetilde{G}$ is
replaced with $G_{pr}$ in the computation of commute-time-based proximity. We
denote this variant as ‘CGNNppr’ and the results of accuracy and efficiency
are reported in Table 3. The findings reveal that the PPR approach is
suboptimal in terms of both accuracy and efficiency, thereby underscoring the
effectiveness of our rewiring-based approach.

Table 3: Accuracy and running time (s) of CGNN and CGNNppr.

Squirrel Chameleon Citeseer CoraML AM-Photo Method Acc. Time Acc. Time Acc.
Time Acc. Time Acc. Time CGNN 77.20 102.93 79.54 112.24 80.27 92.26 77.02
127.15 90.41 121.80 CGNNppr 68.37 257.84 71.69 253.05 78.59 137.82 76.23
192.09 88.52 203.04

Figure 4: Impact of directed structure.

Squirrel CoraML AM-Photo CGNN 77.20 77.02 90.41 CGNNsym 71.31 70.25 88.42

#### Directed vs. Undirected.

To validate the critical role of directed structures in our model, we
transform all directed edges into undirected ones by adding their reverse
counterparts. This process results in a symmetric adjacency matrix, denoted as
$\mathbf{A}_{\mathrm{sym}}$. Subsequently, the commute time is calculated
based on the transition matrix derived from $\mathbf{A}_{\mathrm{sym}}$. We
refer this variant as ‘CGNNsym’. Fig. 4 shows the accuracy of CGNN and CGNNsym
on three datasets. We find that edge direction can significantly influence the
prediction accuracy for our model.

##  6 Conclusion

In this work, we introduce the Commute Graph Neural Network (CGNN) to
integrate commute time into GNNs for digraphs. To this end, we propose DiLap,
a new Laplacian formulation for digraphs, along with a rapid computational
method for determining commute times. By integrating commute times into GNN
message passing through neighbor weighting, CGNN effectively leverages the
path asymmetry in digraphs, thus improving node representation learning. Our
comprehensive experiments demonstrate the superior performance of CGNN over
existing methods.

## References

  * Abu-El-Haija et al. (2019) Abu-El-Haija, S., Perozzi, B., Kapoor, A., Alipourfard, N., Lerman, K., Harutyunyan, H., Ver Steeg, G., and Galstyan, A.  Mixhop: Higher-order graph convolutional architectures via sparsified neighborhood mixing.  In _international conference on machine learning_ , pp. 21–29. PMLR, 2019.
  * Aldous & Fill (2002) Aldous, D. and Fill, J.  Reversible markov chains and random walks on graphs, 2002.
  * Bojchevski & Günnemann (2017) Bojchevski, A. and Günnemann, S.  Deep gaussian embedding of attributed graphs: Unsupervised inductive learning via ranking.  _arXiv preprint arXiv:1707.03815_ , 2017.
  * Cai et al. (2023) Cai, X., Huang, C., Xia, L., and Ren, X.  Lightgcl: Simple yet effective graph contrastive learning for recommendation.  _arXiv preprint arXiv:2302.08191_ , 2023.
  * Chien et al. (2021) Chien, E., Peng, J., Li, P., and Milenkovic, O.  Adaptive universal generalized pagerank graph neural network.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=n6jl7fLxrP.
  * Chung (2005) Chung, F.  Laplacians and the cheeger inequality for directed graphs.  _Annals of Combinatorics_ , 9(1):1–19, 2005.
  * Chung (1997) Chung, F. R.  _Spectral graph theory_ , volume 92.  American Mathematical Soc., 1997.
  * Cross et al. (2001) Cross, R., Borgatti, S. P., and Parker, A.  Beyond answers: Dimensions of the advice network.  _Social networks_ , 23(3):215–235, 2001.
  * Furutani et al. (2020) Furutani, S., Shibahara, T., Akiyama, M., Hato, K., and Aida, M.  Graph signal processing for directed graphs based on the hermitian laplacian.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2019, Würzburg, Germany, September 16–20, 2019, Proceedings, Part I_ , pp. 447–463. Springer, 2020.
  * Geisler et al. (2023) Geisler, S., Li, Y., Mankowitz, D. J., Cemgil, A. T., Günnemann, S., and Paduraru, C.  Transformers meet directed graphs.  In _International Conference on Machine Learning_ , pp. 11144–11172. PMLR, 2023.
  * Halko et al. (2011) Halko, N., Martinsson, P.-G., and Tropp, J. A.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM review_ , 53(2):217–288, 2011.
  * Hamilton (2020) Hamilton, W. L.  _Graph representation learning_.  Morgan & Claypool Publishers, 2020.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, R., and Leskovec, J.  Inductive representation learning on large graphs.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , pp. 1025–1035, 2017.
  * Horn & Johnson (2012) Horn, R. A. and Johnson, C. R.  _Matrix analysis_.  Cambridge university press, 2012.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations (ICLR)_ , 2017.
  * Klicpera et al. (2019) Klicpera, J., Bojchevski, A., and Günnemann, S.  Predict then propagate: Graph neural networks meet personalized pagerank.  In _International Conference on Learning Representations (ICLR)_ , 2019.
  * Li & Zhang (2012) Li, Y. and Zhang, Z.-L.  Digraph laplacian and the degree of asymmetry.  _Internet Mathematics_ , 8(4):381–401, 2012.
  * Ming Chen et al. (2020) Ming Chen, Z. W., Zengfeng Huang, B. D., and Li, Y.  Simple and deep graph convolutional networks.  2020\.
  * Page et al. (1999) Page, L., Brin, S., Motwani, R., and Winograd, T.  The pagerank citation ranking: Bringing order to the web.  Technical report, Stanford InfoLab, 1999.
  * Pei et al. (2019) Pei, H., Wei, B., Chang, K. C.-C., Lei, Y., and Yang, B.  Geom-gcn: Geometric graph convolutional networks.  In _International Conference on Learning Representations_ , 2019.
  * Qiu et al. (2020) Qiu, R., Yin, H., Huang, Z., and Chen, T.  Gag: Global attributed graph neural network for streaming session-based recommendation.  In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pp. 669–678, 2020.
  * Rossi et al. (2023) Rossi, E., Charpentier, B., Giovanni, F. D., Frasca, F., Günnemann, S., and Bronstein, M. M.  Edge directionality improves learning on heterophilic graphs.  In _The Second Learning on Graphs Conference_ , 2023.  URL https://openreview.net/forum?id=T4LRbAMWFn.
  * Rozemberczki et al. (2021) Rozemberczki, B., Allen, C., and Sarkar, R.  Multi-scale attributed node embedding.  _Journal of Complex Networks_ , 9(2):cnab014, 2021.
  * Sen et al. (2008) Sen, P., Namata, G., Bilgic, M., Getoor, L., Galligher, B., and Eliassi-Rad, T.  Collective classification in network data.  _AI magazine_ , 29(3):93–93, 2008.
  * Shchur et al. (2018) Shchur, O., Mumme, M., Bojchevski, A., and Günnemann, S.  Pitfalls of graph neural network evaluation.  _arXiv preprint arXiv:1811.05868_ , 2018.
  * Shuman et al. (2013) Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., and Vandergheynst, P.  The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains.  _IEEE signal processing magazine_ , 30(3):83–98, 2013.
  * Singh et al. (2016) Singh, R., Chakraborty, A., and Manoj, B.  Graph fourier transform based on directed laplacian.  In _2016 International Conference on Signal Processing and Communications (SPCOM)_ , pp. 1–5. IEEE, 2016.
  * Tong et al. (2020a) Tong, Z., Liang, Y., Sun, C., Li, X., Rosenblum, D., and Lim, A.  Digraph inception convolutional networks.  _Advances in neural information processing systems_ , 33, 2020a.
  * Tong et al. (2020b) Tong, Z., Liang, Y., Sun, C., Rosenblum, D. S., and Lim, A.  Directed graph convolutional network.  _arXiv preprint arXiv:2004.13970_ , 2020b.
  * Tong et al. (2021) Tong, Z., Liang, Y., Ding, H., Dai, Y., Li, X., and Wang, C.  Directed graph contrastive learning.  _Advances in neural information processing systems_ , 34:19580–19593, 2021.
  * Topping et al. (2022) Topping, J., Di Giovanni, F., Chamberlain, B. P., Dong, X., and Bronstein, M. M.  Understanding over-squashing and bottlenecks on graphs via curvature.  2022\.
  * Veličković et al. (2018) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y.  Graph Attention Networks.  _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=rJXMpikCZ.  accepted as poster.
  * Wu et al. (2019) Wu, F., Souza, A., Zhang, T., Fifty, C., Yu, T., and Weinberger, K.  Simplifying graph convolutional networks.  In _International conference on machine learning_ , pp. 6861–6871. PMLR, 2019.
  * Xu et al. (2018) Xu, K., Li, C., Tian, Y., Sonobe, T., Kawarabayashi, K.-i., and Jegelka, S.  Representation learning on graphs with jumping knowledge networks.  In _International conference on machine learning_ , pp. 5453–5462. PMLR, 2018.
  * Zhang et al. (2021) Zhang, X., He, Y., Brugnone, N., Perlmutter, M., and Hirn, M.  Magnet: A neural network for directed graphs.  _Advances in neural information processing systems_ , 34:27003–27015, 2021.
  * Zhu et al. (2020) Zhu, J., Yan, Y., Zhao, L., Heimann, M., Akoglu, L., and Koutra, D.  Beyond homophily in graph neural networks: Current limitations and effective designs.  _Advances in Neural Information Processing Systems_ , 33, 2020.

##  Appendix A Proofs and Derivations

###  A.1 Derivation of DiLap T

Since the digraph is unweighted, then
$\mathbf{P}_{ij}=\frac{1}{d^{\text{out}}_{i}}$ and
$\frac{1}{d^{\text{out}}_{i}}\mathbf{A}_{i,:}=\mathbf{P}_{i,:}$. We have:

$\displaystyle(\mathbf{T}s)_{i}$
$\displaystyle=\sum_{v_{j}\in\mathcal{N}^{\text{out}}_{i}}\mathbf{P}_{ij}(s_{i}-s_{j})$
(12)
$\displaystyle=\frac{1}{d^{\text{out}}_{i}}s_{i}-\frac{1}{d^{\text{out}}_{i}}\sum_{v_{j}\in\mathcal{N}^{\text{out}}_{i}}s_{j}$
$\displaystyle=(\mathbf{D}^{-1}s)_{i}-\frac{1}{d^{\text{out}}_{i}}\mathbf{A}_{i,:}s$
$\displaystyle=(\mathbf{D}^{-1}s)_{i}-(\mathbf{P}s)_{i}$
$\displaystyle=((\mathbf{D}^{-1}-\mathbf{P})s)_{i}.$

Thus, DiLap $\mathbf{T}$ can be represented as $\mathbf{D}^{-1}-\mathbf{P}$.

###  A.2 Proof of Lemma 4.1

###### Proof.

Since the transition matrix $\widetilde{\mathbf{P}}$ is row-stochastic, it
follows that $\widetilde{\mathbf{P}}^{t}\mathbf{J}=\mathbf{J}$. In light of
Eq. 2 and considering that $\pi$ is stochastic, we have
$\mathbf{Z}\mathbf{J}=\bm{0}_{n\times n}$. Let
$\mathcal{K}=\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}\widetilde{\mathcal{T}}\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}=\widetilde{\mathbf{\Pi}}^{\frac{1}{2}}(\widetilde{\mathbf{D}}^{-1}-\widetilde{\mathbf{P}})\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}$,
$\mathcal{J}=\widetilde{\mathbf{\Pi}}^{\frac{1}{2}}\mathbf{J}\widetilde{\mathbf{\Pi}}^{\frac{1}{2}}$,
and
$\mathcal{Z}=\widetilde{\mathbf{\Pi}}^{\frac{1}{2}}\mathbf{Z}\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}$,
incorporating these into Eq. 3, we have:

$\mathcal{Z}+\mathcal{J}=(\mathcal{K}+\mathcal{J}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})^{-1}.$
(13)

By post-multiplying Eq. 13 from the right by
$(\mathcal{K}+\mathcal{J}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})$, we have:

$\displaystyle\mathbf{I}$
$\displaystyle=(\mathcal{Z}+\mathcal{J})(\mathcal{K}+\mathcal{J}-\tilde{\mathbf{D}}_{c}^{-1}+\mathbf{I})$
(14)
$\displaystyle=\mathcal{Z}\mathcal{K}+\mathcal{Z}\mathcal{J}-\mathcal{Z}\widetilde{\mathbf{D}}^{-1}+\mathcal{Z}+\mathcal{J}\mathcal{K}+\mathcal{J}^{2}-\mathcal{J}\widetilde{\mathbf{D}}^{-1}+\mathcal{J}.$

Since $\mathcal{J}^{2}=\mathcal{J}$,
$\mathcal{J}\mathcal{K}=\mathcal{J}(\widetilde{\mathbf{D}}^{-1}-\mathbf{I})$,
and
$\mathcal{Z}\mathcal{J}=\widetilde{\mathbf{\Pi}}^{\frac{1}{2}}\mathbf{Z}\mathbf{J}\widetilde{\mathbf{\Pi}}^{\frac{1}{2}}=\bm{0}_{n\times
n}$, Eq. 14 can be simplified to:

$\mathcal{Z}(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})=\mathbf{I}-\mathcal{J}.$
(15)

Similarly, by multiplying from the left, we establish that
$(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})\mathcal{Z}=\mathbf{I}-\mathcal{J}$.
As $\pi^{\top}\mathbf{Z}=\bm{0}$, we have $\mathcal{J}\mathcal{Z}=0$. Thus,
$\mathcal{Z}(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})\mathcal{Z}=\mathcal{Z}$.
Furthermore, it can be straightforwardly inferred that
$(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})\mathcal{J}=0$, then have
$(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})\mathcal{Z}(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})=(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})$.
Besides, considering the symmetry of the left part of Eq. 15, we have
$(\mathcal{Z}(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I}))^{\top}=\mathcal{Z}(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})$.
Similarly,
$((\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})\mathcal{Z})^{\top}=(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})\mathcal{Z}$.
These derivations satisfy the sufficient conditions for the Moore–Penrose
pseudoinverse, such that

$\mathcal{Z}=(\mathcal{K}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I})^{\dagger}.$
(16)

Finally, recovering $\mathcal{Z}$ and $\mathcal{K}$, which concludes the
proof. ∎

###  A.3 SVD for $\mathcal{R}^{\dagger}$

Given a matrix $\mathcal{R}\in\mathbb{R}^{N\times N}$, its Moore-Penrose
pseudoinverse can be directly computed with an SVD-based method. Specifically,
we first perform truncated SVD on
$\mathcal{R}=\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}\widetilde{\mathcal{T}}\widetilde{\mathbf{\Pi}}^{-\frac{1}{2}}-\widetilde{\mathbf{D}}^{-1}+\mathbf{I}\approx\mathbf{U}_{q}\Sigma_{q}\mathbf{V}_{q}^{\top}$,
where $\mathbf{U}_{q}\in\mathbb{R}^{N\times q}$ and
$\mathbf{V}_{q}\in\mathbb{R}^{N\times q}$ contains the first $q$ columns of
$\mathbf{U}$ and $\mathbf{V}$. $\Sigma_{q}\in\mathbb{R}^{q\times q}$ is the
diagonal matrix of $q$ largest singular values. It is a $q$-rank approximation
of $\mathcal{R}$, which holds that $\mathrm{rank}(\mathrm{R})=q$. Then the
Moore-Penrose pseudoinverse of $\mathcal{R}$ can be easily computed as
follows:

$\mathcal{R}^{\dagger}=\mathbf{U}_{q}\Sigma_{q}^{-1}\mathbf{V}_{q}^{\top}.$
(17)

To leverage sparsity of $\mathcal{R}$ to avoid $\mathcal{O}(N^{3})$
complexity, we adopt the randomized SVD algorithm proposed by (Halko et al.,
2011; Cai et al., 2023) to first approximate the range of the input matrix
with a low-rank orthonormal matrix, and then perform SVD on this smaller
matrix:

$\hat{\bm{U}}_{q},\hat{\Sigma}_{q},\hat{\bm{V}}_{q}^{\top}=\operatorname{ApproxSVD}(\mathcal{R},q),\quad\hat{\mathcal{R}}_{SVD}=\hat{\bm{U}}_{q}\hat{\Sigma}_{q}\hat{\bm{V}}_{q}^{\top},$
(18)

where $\hat{\bm{U}}_{q}$, $\hat{\Sigma}_{q}$, and $\hat{\bm{V}}_{q}$ are the
approximated versions of $\mathbf{U}_{q}$, $\Sigma_{q}$, and $\mathbf{V}_{q}$.
Then the Moore-Penrose pseudoinverse of $\mathcal{R}$ can be computed by:

$\mathcal{R}^{\dagger}=\hat{\bm{U}}_{q}\hat{\Sigma}^{-1}_{q}\hat{\bm{V}}_{q}^{\top}.$
(19)

The computation cost of randomized truncated SVD takes $\mathcal{O}(qK)$,
where $K$ is the number of non-zero elements in $\mathcal{R}$, so we have
$K=|E|$. Thus, the sparsity degree of $\mathcal{R}$ can determine the time
complexity of its Moore-Penrose pseudoinverse, which demonstrates the
importance of Lemma 4.1.

##  Appendix B Pseudo Code for CGNN

Algorithm 1 CGNN

Input: Digraph $G=(V,E,\mathbf{X})$; Depth $L$; Hidden size $d^{\prime}$;
Number of classes $K$
Output: Logits $\hat{Y}\in\mathbb{R}^{N\times K}$

1: Compute the anchor $\bm{a}$ and node-anchor similarities to construct
$G^{\prime}$ with Eq. 7.

2: Add all edges from $G^{\prime}$ to $G$ to generate $\widetilde{G}$.

3: Compute the Weight DiLap $\widetilde{\mathcal{T}}$ for $\widetilde{G}$ with
Eq. 6.

4: Compute $\mathcal{R}$ and its Moore-Penrose pseudoinverse with Eq. 8 and
Eq. 19.

5: Compute the commute time matrix $\mathcal{C}$ with Eq. 10.

6: Compute the normalized proximity matrix $\widetilde{\mathcal{C}}$ with
$\widetilde{\mathcal{C}}^{\text{out}}=\mathbf{A}\odot\widetilde{\mathcal{C}}$
and
$\widetilde{\mathcal{C}}^{\text{in}}=\mathbf{A}^{\top}\odot\widetilde{\mathcal{C}}$.

7: for $\ell\in\\{1,\cdots,L\\}$ do

8: Layer-wise message passing with Eq. 11.

9: end for

10: $\mathbf{H}=\mathrm{MLP}(\mathbf{H}^{(L)})$.

11: $\hat{Y}=\mathrm{Softmax}(\mathbf{H})$.

#### Complexity Analysis

The time complexity of randomized truncated SVD to compute
$\mathcal{R}^{\dagger}$ is $\mathcal{O}(q|E|)$, and the message passing
iteration has the same time complexity as GraphSAGE with $\mathcal{O}(|E|)$.
Therefore, the overall time complexity of CGNN is $\mathcal{O}(q|E|)$.

##  Appendix C Implementation Details

###  C.1 Experimental Settings

We evaluate the performance by node classification accuracy with standard
deviation in the semi-supervised setting. For Squirrel and Chameleon, we use
10 public splits (48%/32%/20% for training/validation/testing) provided by
(Pei et al., 2019). For the remaining datasets, we adopt the same splits as
(Tong et al., 2020a, 2021), which chooses 20 nodes per class for the training
set, 500 for the validation set, and allocates the rest to the test set. We
conduct our experiments on 2 Intel Xeon Gold 5215 CPUs and 1 NVIDIA GeForce
RTX 3090 GPU.

###  C.2 Data Statistics

The dataset used in Section 5 are Squirrel, Chameleon (Rozemberczki et al.,
2021), Citeseer (Sen et al., 2008), CoraML (Bojchevski & Günnemann, 2017), and
AM-Photo (Shchur et al., 2018). We summarize their statistics in Table 4.
$\mathrm{homo\\_ratio}$ represents the homophily ratio, a metric proposed by
Zhu et al. (2020). which is employed to gauge the degree of homophily within
the graph. A lower $\mathrm{homo\\_ratio}$ signifies a greater degree of
heterophily, indicating a higher prevalence of edges that connect nodes of
differing classes.

Table 4: Statistics of the datasets. Dataset | $N$ | $|E|$ | # Feat. | # Classes | $\mathrm{homo\\_ratio}$
---|---|---|---|---|---
Squirrel | 5,201 | 217,073 | 2,089 | 5 | 0.22
Chameleon | 2,277 | 36,101 | 2,325 | 5 | 0.23
Cora-ML | 2,995 | 8,416 | 2,879 | 7 | 0.79
Citeseer | 3,312 | 4,715 | 3,703 | 6 | 0.74
AM-Photo | 7,650 | 238,162 | 745 | 8 | 0.83

###  C.3 Hyperparameter Settings

For our model, we tune the hyperparameters based on the highest average
validation accuracy. We utilize the randomized truncated SVD algorithm for
computing the Moore-Penrose pseudoinverse of matrix $\mathcal{R}$, setting the
required rank $q$ to 5 for all datasets. The learning rate $lr$ is selected
from $\\{0.01,0.005\\}$, and the weight decay $wd$ from $\\{0,5e-5,5e-4\\}$.
In the model architecture, the number of layers $L$ vary among
$\\{1,2,3,4,5\\}$ and the dimension $d^{\prime}$ is selected from
$\\{32,64,128,256,512\\}$. The comprehensive hyperparameter configurations for
CGNN are detailed in Table 5.

Table 5: Hyperparameters specifications. Dataset | $lr$ | $wd$ | $L$ | $d^{\prime}$
---|---|---|---|---
Squirrel | 0.005 | 0 | 5 | 128
Chameleon | 0.01 | 0 | 4 | 128
CoraML | 0.01 | 0 | 2 | 64
Citeseer | 0.01 | 0 | 2 | 128
AM-Photo | 0.005 | 0 | 2 | 512

##  Appendix D Additional Experiments

Table 6: Link prediction results.

Direction prediction Existence prediction Method Citeseer CoraML Citeseer
CoraML GCN 77.69$\pm$0.73 75.68$\pm$0.97 70.27$\pm$2.12 80.67$\pm$1.08 GAT
76.09$\pm$0.09 50.02$\pm$0.05 50.82$\pm$0.36 50.00$\pm$0.01 GraphSAGE
75.16$\pm$0.91 70.07$\pm$1.43 70.59$\pm$1.74 69.32$\pm$0.58 APPNP
79.23$\pm$0.35 78.64$\pm$0.70 81.14$\pm$1.21 85.46$\pm$0.80 MixHop
79.40$\pm$0.42 78.62$\pm$1.25 80.11$\pm$1.93 85.21$\pm$0.90 GPRGNN
79.15$\pm$0.35 77.88$\pm$1.30 78.15$\pm$0.24 83.83$\pm$0.74 GCNII
78.60$\pm$0.69 75.51$\pm$1.90 79.48$\pm$1.74 75.51$\pm$1.90 DGCN
82.38$\pm$0.55 74.13$\pm$1.07 81.11$\pm$2.28 80.30$\pm$1.25 DiGCN
83.06$\pm$0.70 74.92$\pm$0.88 85.08$\pm$1.26 81.06$\pm$1.16 MagNet
79.36$\pm$0.44 77.71$\pm$1.04 88.20$\pm$1.05 87.21$\pm$0.72 DiGCL
81.89$\pm$0.76 76.05$\pm$0.97 82.96$\pm$1.32 84.23$\pm$1.26 DirGNN
82.09$\pm$0.88 77.93$\pm$0.76 85.18$\pm$1.02 83.37$\pm$1.83 CGNN
84.97$\pm$0.71 80.33$\pm$0.90 85.27$\pm$1.30 82.08$\pm$1.22

###  D.1 Link Prediction

For link prediction, we compare our method with all baselines on Citeseer and
CoraML datasets. Our experimental setup follows the protocol established by
(Zhang et al., 2021), considering two distinct types of link prediction tasks
for performance evaluation. The first task involves predicting the direction
of edges between pairs of nodes, while the second focuses on predicting the
existence of edges. For both tasks, we partition the data by removing 15% of
the edges for testing and 5% for validation, with the remaining edges used for
training. The results, detailed in Table 6, reveal that CGNN excels in the
direction prediction task. Meanwhile, in the edge existence prediction task,
MagNet yields the best results, with CGNN also achieving comparable
performance. These findings show CGNN’s effectiveness in accurately capturing
directed structural patterns in graphs.

###  D.2 Sensitivity analysis

We investigate the sensitivity of CGNN to key hyperparameters that influence
its performance, specifically focusing on the number of layers $L$ and the
dimension of the hidden layer $d^{\prime}$. We explore a range of values for
$L$, considering $\\{1,2,3,4,5\\}$, and for $d^{\prime}$, considering
$\\{32,64,128,256,512\\}$. From Fig. 5, we observe that we observe that CGNN
achieves optimal performance with $L=5$ and $d^{\prime}=128$ on Squirrel, and
with $L=2$ and $d^{\prime}=64$ on CoraML. This suggests that deeper models are
necessary to effectively aggregate valuable information in heterophilic
graphs, whereas in homophilic graphs, leveraging local neighborhood
information is generally adequate.

(a)

(b)

Figure 5: Sensivisity analysis on Squirrel and CoraML.