Commute Graph Neural Networks

## Abstract

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 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. Extensive experiments confirm the superior performance of CGNN.

## 1 Introduction

Directed graphs (digraphs) are widely employed to model relational structures in diverse domains, such as social networks (Cross et al., ; Zhuo & Tan, , ), recommendation systems (Qiu et al., ) and fraud detection (Zhuo et al., ). Recently, the advances of graph neural networks (GNNs) have inspired various attempts to adopt GNNs for analyzing digraphs (Tong et al., , , ; Zhang et al., ; Rossi et al., ; Geisler et al., ; Zhuo & Tan, ). 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., ), or by predefining edge-specific weights (Zhang et al., ; Tong et al., ).

Figure: 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 SPDs from $v_{i}$ to $v_{j}$ and $v_{k}$ are identical, the reverse SPD from $v_{j}$ and $v_{k}$ back to $v_{i}$ do not necessarily match these distances.
Refer to caption: x1.png

Despite the advancements, current digraph neural networks primarily capture unidirectional(^1^11‘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., ; Tong et al., ; Zhang et al., ), 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, ) 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., ) 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, ) 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., ) 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, ).

### 2.2 Digraph Neural Networks

To effectively capture the directed structure with GNNs, spectral-based methods (Zhang et al., ; Tong et al., , ) have been proposed to preserve the underlying spectral properties of the digraph by performing spectral analysis based on the digraph Laplacian proposed by (Chung, ). (Koke & Cremers, ) introduce holomorphic filters as spectral filters for digraphs, and investigate their optimal filter-bank. MagNet (Zhang et al., ) and its extensions (Lin & Gao, ; Fiorini et al., ) 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., ) allows for controlling the direction of information flow by considering in-neighbors or out-neighbors separately. DirGNN (Rossi et al., ) further extends this framework by segregating neighbor aggregation according to edge directions, offering a more refined method to handle the directed nature of graphs. Transformer-based methods capture directed structure by specific positional encoding modules, such as directional random walk encoding (Geisler et al., ) and partial order encoding (Luo et al., ). DUPLEX (Ke et al., ) utilizes Hermitian adjacency matrix decomposition for neighbor aggregation and incorporates a dual GAT encoder for modeling directional neighbors.

## 3 Preliminary

#### 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., ) 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., ) 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}$.

## 4 Random Walk Distance and GNNs

Based on the established 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.

### 4.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., ). 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}$.

### 4.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, ). We first establish the following assumptions required to support the theorem.

###### Assumption 4.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 4.2 .

(Li & Zhang, )
Given 4.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, ) 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.

###### Assumption 4.1 .

###### Proposition 1 .

###### Theorem 4.2 .

## 5 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.

### 5.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., ; Hamilton, ): 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., ). Therefore, $\pi$ effectively indicates the centrality and influence over the long term in the graph. Perron-Frobenius Theorem (Horn & Johnson, ) 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}$.

### 5.2 Similarity-based Graph Rewiring

Both the fundamental matrix defined in Eq. 3 and Weighted DiLap requires 4.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., ) introduce a teleporting probability uniformly distributed across all nodes. This method, inspired by PageRank (Page et al., ), 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: 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}$.
Refer to caption: x2.png

This rewiring approach satisfies 4.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}}$.

### 5.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 5.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 5.1 and using Eq. 4, we can further compute the hitting times and commute times in terms of $\widetilde{\mathcal{T}}$ with the following theorem.

###### Theorem 5.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}}}.$ $$

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 4.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., ; Cai et al., ), 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}$.

###### Lemma 5.1 .

###### Theorem 5.2 .

### 5.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 of CGNN in Algorithm 1.

Figure: Algorithm 1 CGNN

#### Complexity Analysis

The randomized truncated SVD to compute $\mathcal{R}^{\dagger}$ is $\mathcal{O}(q|E|)$ where $q$ is the required rank, and the message passing iteration has the same time complexity as DirGNN with $\mathcal{O}(L|E|d^{2})$. Therefore, the overall time complexity of CGNN is $\mathcal{O}((Ld^{2}+q)|E|)$. In practice, $q$ is typically set to 5, rendering the time complexity effectively linear with respect to the number of edges $|E|$. In GNN domain, models with a complexity less than $\mathcal{O}(N^{2})$ are generally considered feasible by researchers (Wu et al., ). Given that real-world networks are often extremely sparse, i.e., $|E|\ll N^{2}$, CGNN demonstrates its feasibility as a model within the GNN family.

## 6 Experiments

### 6.1 Experimental Settings

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

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., ). For the remaining datasets, we adopt the same splits as (Tong et al., , ), 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.

### 6.2 Data Statistics

The datasets used in Section 6 are Squirrel, Chameleon (Rozemberczki et al., ), citepseer (Sen et al., ), CoraML (Bojchevski & Günnemann, ), AM-Photo (Shchur et al., ), Snap-Patents, Roman-Empire, and Arxiv-Year (Rossi et al., ). We summarize their statistics in Table 1. $\mathrm{homo_ratio}$ represents the homophily ratio, a metric proposed by (Zhu et al., ). 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 1: 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 |
| citepseer | 3,312 | 4,715 | 3,703 | 6 | 0.74 |
| AM-Photo | 7,650 | 238,162 | 745 | 8 | 0.83 |
| Snap-Patents | 2,923,922 | 13,975,791 | 269 | 5 | 0.22 |
| Roman-Empire | 22,662 | 44,363 | 300 | 18 | 0.05 |
| Arxiv-Year | 169,343 | 1,166,243 | 128 | 40 | 0.22 |

Table: Table 2: 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. OOM indicates out-of-memory.

### 6.3 Overall Results and Analysis

Table 2 reports the node classification results across eight digraph datasets. Our method CGNN achieves new state-of-the-art results on 6 out of 8 datasets, and comparable results on the remaining ones, 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 (e.g., Squirrel and Chameleon), supporting the findings of (Rossi et al., ) 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 3 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: Table 3: The distance between $\mathcal{M}$ and $\mathbf{A}$, and between $\mathcal{M}$ and $\widetilde{\mathcal{C}}$.

#### Application scope analysis.

Can commute times always enhance message passing on directed graphs? To answer this, we analyze the scope of use for CGNN based on Table 2. For example, on the Snap-Patents and CoraML dataset, we observed that adding commute time-based weights during message passing did not significantly enhance performance. Now we can analyze the reason from the perspective of dataset. CoraML is a directed citation network where nodes predominantly link to other nodes within the same research area. However, in such networks, reciprocal citations between two papers are impossible due to their chronological sequence. Consequently, mutual path dependencies do not exist, and thus, incorporating commute times to adjust neighbor weights might might (slightly) hurt performance. A similar situation exists with the Snap-Patents dataset, where each directed edge represents a citation from one patent to another, again indicating the absence of mutual path dependencies.

In conclusion, in networks like citation networks where mutual relationships inherently do not exist, applying commute times to enforce these relationships is unnecessary. Conversely, our model is particularly effective in networks like webpage networks and social networks—examples being Squirrel and AM-Photo—where mutual relationships are prevalent. For instance, in a social network, an ordinary user may follow a celebrity, creating a short path to the celebrity. However, the reverse path from the celebrity back to the user might be considerably longer. Therefore, in such networks, considering mutual relationships based on commute times can provide a more accurate description of node relationships.

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

### 6.4 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.

### 6.5 Component Analysis

Table: Table 5: Changes in commute times before and after rewiring.

#### Comparison between graph rewiring and PPR.

In Section 5.2, we construct a rewired graph $\widetilde{G}$ based on feature similarity to guarantee the irreducibility and aperiodicity. This approach introduces at most two additional edges per node, specifically targeting those with the highest feature similarity, while minimally altering the original graph structure to preserve semantic information. To investigate changes in commute times before and after rewiring, for the original graph, we use its largest connected component, removing absorbing nodes (i.e., nodes without outgoing edges) to ensure that we can compute meaningful and deterministic commute times. We denote the average normalized commute time for the original graph as $c_{orig}$; For the rewired graph, we directly compute the commute time and denote this average normalized value as $c_{rew}$. We then use $\delta=\frac{\|c_{orig}-c_{rew}\|_{2}}{\|c_{orig}\|}$ to quantify the changes, which can be intepretated as the proportion of commute information changed in the original graph. As shown in Table 5, the graph rewiring method can effectively preserve the original commute times of the graph.

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 4. 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: Table 6: Impact of directed structure.

#### 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’. Table 6 shows the accuracy of CGNN and CGNNsym on three datasets. We find that edge direction can significantly influence the prediction accuracy for our model.

## 7 Conclusion

Identifying and encoding asymmetric mutual path dependencies in directed graphs is essential for accurately representing real-world relationships between entities. In this work, we utilize the concept of commute time to assess the strength of relationships in directed graphs and introduce the Commute Graph Neural Network (CGNN) to incorporate node-wise commute time into node representations. To achieve this, we propose DiLap, a novel Laplacian formulation derived from the divergence of the gradient of signals on directed graphs, along with an efficient computational method for deterministic commute times. By integrating commute times into GNN message passing through neighbor weighting, CGNN effectively harnesses path asymmetry in directed graphs, thereby enhancing node representation learning. Our extensive experiments across eight directed graph datasets demonstrate that CGNN significantly outperforms existing methods.

## 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 5.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.
∎

###### Proof.

### A.3 SVD for ℛ † superscript ℛ † \mathcal{R}^{\dagger} caligraphic_R start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT

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., ; Cai et al., ) 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 5.1.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: 1] The derivation of your DiLap operator appears to be flawed. In particular, in the second line of Equation (12) when you pull $s_i$ out of the sum, the term $s_i$ arises $d_i^{out}$ times and therefore, Line 2 should be $s_i - \frac{1}{d_i^{out}} \sum s_j.$ The correction of this error means that you should be working with the random walk Laplacian (see e.g. [1]), which would be far more intuitive. To me it seems that for it to be possible to accept this paper at ICLR, the derivation of your operator needs to be corrected and the subsequent experiments should be adjusted. 

2] I am unsure what adjacency matrix you use in Line 314 to sparsify the matrix $\tilde{\mathcal{C}}.$ Are you using the adjacency matrix corresponding to the graph in which you have added the node feature similarity edges? And if not, could you hypothesise how severe the impact may be of calculating the commute times on a rewired graph and to then message pass with the original graph. 

3] Your method appears to be relatively memory intensive. In particular, you seem to require the evaluation of the exponential function fo the dense matrix $\tilde{\mathcal{C}}$. Empirical evaluation of, not only the time, but also memory complexity of your method in comparison to your baseline methods would be very valuable. 

4] The ablation study in Table 2 is very interesting! I think it should be extended in scope to also extend to homophilic datasets and to also include other commonly used message passing operators, such as the symmetrically normalised adjacency matrix used in the GCN and the PageRank matrix used in the PPrGo model [2].

5] It does not seem sensible to me to compare your sparisfied commute time based CGNN to the CGNN$_{ppr}$ using the dense PageRank matrix. It seems to be fairer to me to either compare dense versions of both matrices or sparse versions of both matrices. In particular, since you sparify your commute time matrix with the adjacency matrix, it would be interesting to compare your model to a PageRank-based scheme, where the PageRank matrix is also sparsified with the adjacency matrix. 

6] Minor comments:

6.1] The abbreviation "SPD" is used in Line 45 before its definition. 

6.2] The contributions of your paper are not explicitly listed. 


[1] Von Luxburg, U., 2007. A tutorial on spectral clustering. Statistics and computing, 17, pp.395-416.

[2] Bojchevski, A., Gasteiger, J., Perozzi, B., Kapoor, A., Blais, M., Rózemberczki, B., Lukasik, M. and Günnemann, S., 2020, August. Scaling graph neural networks with approximate pagerank. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (pp. 2464-2473).
soundness: 2
strengths: - The rewiring scheme that you propose is simple, but rather nice in my opinion. It would be interesting to see further study of its impact on the overall graph structure. 

- Your proposed CGNNs are compared to a comprehensive set of baseline models, which great to see.  

- The analysis of the application scope of your proposed model is very strong indeed and something that is generally not done enough in our literature.
summary: In the submitted manuscript, the authors propose a novel digraph Laplacian, which is later used to more efficiently calculate the commute time of pairs of nodes. They furthermore propose a simple node feature based rewiring scheme, which allows them to ensure that the resulting graph gives rise to an aperiodic, irreducible Markov Chain, which has a unique steady state. The authors then propose to calculate commute times on this rewired graph, to transform these commute times by taking the exponential function of this matrix and subsequently sparsifying it with the adjacency matrix. This then allows the authors to propose a variant of the DirGNN, called CGNN, in which edges are reweighted by their transformed commute times. The authors finally evaluate the empirical performance of their CGNNs against a large variety of baseline models on a large number of datasets and find consistently good, although sometimes marginal performance improvements. They furthermore analyse these results and provide several insightful further experiments on runtimes and ablation studies of different model components.
weaknesses: Please find further details on my listed weaknesses in my questions below.
- The proposed model boils down a weighting of the DirGNN by an efficiently calculated function of the commute times, which is a rather trivial change. 
- The derivation of the DiLap matrix appears to be flawed. 
-  Your ablation studies could be improved and extended.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: As outline in the weaknesses above, can the authors provide a principled reason for why features aggregated from a node's neighbors should be weighed by the commute distance of the node to its neighbors? Can the authors provide any theoretical justification for using commute distance for such weighing?

The authors provide one figure where they show for Chameleon and Squirrel data sets that an adjacency matrix constructed from nearest neighbors weighted by their commute distance more closely resembles an adjacency matrix constrained to connecting nodes within the same class. Can the authors include other empirical measures of how the weighting of neighbors can improve performance of GNNs? Can the authors look at properties of longer hops (going beyond one-hop neighbors) and how information is aggregated across nodes within or outside the same class? Would it be possible to construct a synthetic data set that would shed light on the mechanism behind why they see an improvement?

How does the proposed model's performance change with depth? The authors claim that the weighted neighbors avoids the problem of aggregation of irrelevant information as depth is increased? It would be informative to see an empirical demonstration of this. The authors should plot the performance of their model as a function of model depth and compare with existing models.

In addition, can the authors apply their approach a real world problem on directed graphs that requires long range information transmission, such as power grids or traffic flow data? This would bolster the empirical support for their method.

As noted in the weaknesses above, can the authors try alternative methods for graph rewiring that also produce sparse graphs, such as constructing a kNN graph using the node features? The authors should include an empirical comparison or provide theoretical justification for their method. Is the proposed similarity-based rewiring mode optimal?

How much of the improvement is coming from the fact that the Laplacian proposed by the authors is weighted by the stationary probability of random walks on the graph (or what the authors call the importance of a node)? Can the authors do an ablation study to disentangle this from the commute distance?

This reviewer was confused by the comment that general GNNs can outperform models tailored for directed graphs with hyper-parameter tuning. In Table 1, DirGNN, a model tailored for directed graphs, outperforms GCN. Can the authors clarify what they mean by this comment?
soundness: 2
strengths: This is an interesting paper with a clever insight. It is sensible to say that not all nearest neighbors on a directed graph are created equal. Weighing the features of some neighbors more during the aggregation step of a GNN based on the shorter commute distance of those neighbors to the original node is an intriguing idea. Weighing based on the commute distance certainly sounds reasonable. The proposed method for rewiring a graph to make it irreducible and aperiodic while still retaining sparsity is also clever as is the weighted Laplacian that can be used to efficiently compute the commute distance leveraging its sparsity using methods such as randomized truncated singular value decomposition.

The state-of-the-art performance achieved using the author's method on some of the most commonly used directed graph data sets, such as Squirrel, Chameleon, etc. is impressive and provides reasonable empirical proof of the validity of the proposed approach. The authors also include solid empirical evidence on running times of their algorithm and convincing comparisons to PageRank for graph rewiring and calculation of commute distances.
summary: In this paper, the authors propose a method for weighing the features of the neighbors of a node during aggregation step of GNN based on the commute distance of the neighbor to the node. The commute distance between nodes A and B is the average number of steps that a random walk takes traversing from node A to node B and back to node A. This distance is particularly relevant for directed graphs because although all nearest neighbors are one hop away from a node, their commute distances might vary because of the constraints imposed by the directions of the edges. One neighbor might require a circuitous path along many other nodes before returning to the original node (have longer commute distance) whereas another neighbor could be closer. The authors' key idea is to weight the importance of the features of the neighbors of a node during the aggregation step of GNN based on the commute distance of the node to those neighbors. Besides this weighing, the aggregation and update scheme that they use is based on that of Rossi et al. where the features of the incoming and outgoing neighbors are aggregated separately and used alongside the node's own features during the update step. The authors also propose an efficient way of computing commute distance. To do so, they introduce a weighted Laplacian for directed graphs that accounts both for the directional connectivity of the nodes and their importance (computed as the stationary probability at each node of a random walk on the graph). The authors also introduce a way to rewire a graph to ensure that it is irreducible and aperiodic while keeping the graph sparse (unlike alternative methods such as PageRank). The commute distance can then be efficiently computed using the sparse weighted Laplacian. Finally, the authors show empirically that their proposed approach improves on existing methods when applied to many standard directed graph data sets such as Squirrel and Chameleon.
weaknesses: To this reviewer, the biggest weakness of the paper was that although weighing neighbors by commute time is sensible, it is not necessarily principled. Is there any reason to a priori expect that neighbors of a node that have shorter commute times to that node somehow contain more relevant features for learning on graphs? This seems to depend on the nature of the learning problem and the data set. Now, the authors can argue that their empirical evidence is sufficient to motivate their approach. However, more should be done here to support the author's proposal. Some evidence is provided in that an adjacency matrix constructed by weighing the neighbors by their commute distance more closely resembles an adjacency matrix constrained to edges that connect nodes within the same class. The authors should expand on this. What does this look like for other data sets? How does aggregation of information using these weighted neighbors look across multiple hops and longer distances across the graph? The authors should have come up with synthetic data sets that can elucidate the mechanism behind the improvement that they are seeing.

If empirical evidence is the main motivation behind the proposed schemes, the authors could have dome more to build a stronger case. An argument is made in the paper that with the weighing of the neighbors less irrelevant information is aggregated as the GNN models go deeper. The authors should empirically demonstrate this by showing how the performance of their model changes with depth and contrast with existing models. In general, it would have been very interesting to see the impact of the weights proposed in this paper on multi-hop GNN models, such MixHop, Shortest Path Networks, or DRew. In addition, it would have been more convincing if the authors had applied their approach to real world problems of directed graphs in addition to standard benchmark data sets used in Table 1 such as temporal web traffic data, power grids, traffic flow, etc.

The method proposed in the paper to rewire the graph to ensure irreducible and aperiodic graphs is also ad hoc and not very principled. The method certainly produces a sparse graph unlike PageRank, however, other approaches can also be used to generate irreducible graphs that are sparse such as generating a kNN graph based on the node features. It is not clear that the proposed method is optimal in any way other than that it outperforms the amended probabilities used in PageRank. See Questions below for suggestions on how the authors can evaluate this empirically.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: I have no further questions, though I would recommend that the authors address the previously noted weaknesses.
soundness: 3
strengths: This paper is engaging and well-written, with a thorough background review that enhances accessibility and readability. Key strengths I noted include:

1. **Novel Approach with Significant Potential:** The proposed method, particularly the newly formulated digraph Laplacian, offers a fresh perspective with substantial potential for future research and applications.

2. **Comprehensive Component Analysis (Section 5.3):** The inclusion of a component analysis strengthens the paper by providing an effective ablation study.

3. **Clear Contribution and Baseline Comparison:** The authors clearly articulate their contributions, outlining the distinctions between their method and existing baselines. They explain where prior approaches fall short and demonstrate how their approach addresses these limitations.

4. **Effective Visual Aids:** Figures 1 and 2 are well-designed and enhance understanding by clarifying details within the method.

5. **Robust Experimental Validation:** he paper validates its approach across a wide variety of datasets and multiple baseline comparisons, highlighting the robustness and generalizability of the proposed method.

6. **Reproducibility:** The authors provide code for reproducing their experiments.
summary: The authors present a novel approach that integrates node-wise commute time into a message-passing framework, introducing the Commute Graph Neural Network (CGNN). The central contribution of CGNN is the development of a new directed graph Laplacian, designed to address path asymmetry in directed graphs. The authors demonstrate that CGNN outperforms existing baselines across most datasets and effectively motivates the significance of the problem they address.

Overall, I found the paper well-executed and recommend it for acceptance.
weaknesses: Overall, this paper is strong in its methodology and results, though I have a few recommendations that could enhance its clarity and depth.

1. **Graph Density in Rewiring Approach:** While I appreciate that the authors provided commute times before and after rewiring in Table 3, it would be beneficial to also examine how rewiring affects graph density. This additional metric could offer deeper insights into structural changes post-rewiring.

2. **Unobserved Edges in Definition of $m_{i,in}^{(l)}$ and $m_{i,out}^{(l)}$:** Given that unobserved edges are introduced to the graph, I suggest adjusting the definitions of $m_{i,in}^{(l)}$ and $m_{i,out}^{(l)}$  to account for these edges, potentially assigning them a lower weight than observed edges. This adjustment could yield a more realistic representation of edge significance.

3. **Model Complexity:** The model’s complexity is relatively high, even though it’s reported to be on par with other GNN models. This complexity, particularly in precomputation, might be a barrier in some cases. However, I do not consider this a critical issue, as future work could address and optimize this aspect.

4. Inclusion of Synthetic Datasets: While the paper impressively covers a range of empirical datasets, the addition of synthetic datasets could improve interpretability. By embedding known patterns, synthetic data could highlight the model's strengths and limitations in detecting specific features.

5. Reordering Related Work: Placing the Related Work section (currently Section 6) closer to the beginning would make the reading experience smoother, giving readers essential context before diving into the methodology and results.

These revisions would, in my opinion, strengthen the paper without diminishing its core contributions.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1) Why is the proposed Laplacian sparse? From Eq. (5), the matrix $P$ seems to be a complete matrix.
2) What is the relationship between Eq. (5) and $D^{-2}L$ and why do you should choose Eq. (5)?
3) Why does the rewiring procedure only minimally alter the overall semantics of the original graph?
soundness: 2
strengths: 1) The idea of considering commute time is novel and reasonable.
2) The topic of directed graph neural networks is important.
3) The source code is provided.
summary: This paper proposes  Commute Graph Neural Networks (CGNN) for directed graphs, which is based on a new digraph Laplacian matrix taking into the commute time on a (possibly rewired) strongly connected graph. Theoretical and empirical analysis is provided.
weaknesses: 1) It is unclear why the proposed Laplacian is sparse. From Eq. (5), the matrix $P$ seems to be a complete matrix.
2) It is unclear what the relationship is between Eq. (5) and $D^{-2}L$ and why you should choose Eq. (5).
3) Being strongly connected is too strong an assumption, and it is not clear why the rewiring procedure only minimally alters the overall semantics of the original graph.
4) [1] mentions flow imbalance in directed graphs and is not discussed. It is also unclear whether the idea in [1] is considered unidirectional by the authors.
5) (minor) Grammar issues: e.g., line 115.5 "notations. We" should be "notations, we"

Reference:
 [1] He, Y., Reinert, G., & Cucuringu, M. (2022, December). DIGRAC: digraph clustering based on flow imbalance. In Learning on Graphs Conference (pp. 21-1). PMLR.