Differentiable Cluster Graph Neural Network

## Abstract

Abstract Graph Neural Networks often struggle with long-range information propagation and in the presence of heterophilous neighborhoods.
We address both challenges with a unified framework that incorporates a clustering inductive bias into the message passing mechanism, using additional cluster-nodes.
Central to our approach is the formulation of an optimal transport based implicit clustering objective function. However, the algorithm for solving the implicit objective function needs to be differentiable to enable end-to-end learning of the GNN. To facilitate this, we adopt an entropy regularized objective function and propose an iterative optimization process, alternating between solving for the cluster assignments and updating the node/cluster-node embeddings.
Notably, our derived closed-form optimization steps are themselves simple yet elegant message passing steps operating seamlessly on a bipartite graph of nodes and cluster-nodes.
Our clustering-based approach can effectively capture both local and global information,
demonstrated by extensive experiments on both heterophilous and homophilous datasets.

## 1 Introduction

Graph Neural Networks (GNNs) have emerged as prominent models for learning node representations on graph-structured data. Their architectures predominantly adhere to the message passing paradigm, where node embeddings are iteratively refined using features from its adjacent neighbors [, , ].
While this message passing paradigm has proven effective in numerous applications [, ], two prominent challenges have been observed.
First, long-range information propagation over a sparse graph can be challenging [, , ].
Expanding the network’s reach by increasing the number of layers is often suboptimal as it could encounter issues such as
over-squashing [, ], where valuable long-range information gets diluted as it passes through the graph’s bottlenecks, diminishing its impact on the target nodes.
Second, some graphs exhibit heterophily, where connected nodes are likely to be dissimilar. In such cases, aggregating information from the dissimilar neighbors could potentially be harmful and hinder the graph representation learning performance [, ].

In this paper, we focus on the task of supervised node classification and explore clustering as an inductive bias to address both challenges. We observe that cluster patterns appear globally when nodes that are far apart in the graph have similar features (see Fig. 1).
By clustering nodes into multiple global clusters based on their features, we can establish connections between these distant nodes.
Additionally, clusters may also be present locally among a node and its neighbours in the graph. Particularly for heterophilic graphs, there are likely more than one cluster in the local neighbourhood.
Leveraging these local clusters for information aggregation, rather than across the entire neighborhood at once, can potentially help to avoid the problems caused by the heterophilic nature of the graph.

To embed the clustering inductive bias, we propose the use of a differentiable clustering-based objective function consisting of a weighted sum of global and local clustering terms. The global objective encourages clustering on all the nodes in the graph, while the local objective does so on subsets of nodes, where each subset consists of a node in the graph and its immediate neighbors. Optimizing the objective function would produce a set of node embeddings that can then be used with a readout function for classifying the nodes in the graph.

Figure: Figure 1: On the left is an instance where distant nodes are similar to each other. On the right is the heterophilous ego-neighborhood of node A where cluster patterns appear. Square boxes indicate conceptual cluster centroids.
Refer to caption: extracted/5620147/arxiv_paper/Figures/intro6.png

To optimize the implicit clustering-based objective function,
we employ a block coordinate descent approach by iteratively, (1) solving for a soft cluster assignment matrix that probabilistically assigns a node to a cluster, followed by, (2) updating the node and cluster embeddings given the cluster assignment matrix. We view the first step of solving for the assignment matrix as an Optimal Transport (OT) problem [], which aims to find the most cost-effective way to move from one probability distribution to another. In our context, the cost is defined as the distance between nodes and cluster-centroids, and the goal is to transport nodes to cluster centroids with minimum cost. We adopt the entropic regularized Sinkhorn distance [] approximation for solving the
OT problem. This approximation utilizes differentiable operations, ensuring end-to-end learning.

For the second step of updating the node and cluster embeddings, we derive a closed-form solution that minimizes the objective function given the assignment matrix.
The derived block coordinate descent algorithm can be viewed as a message passing algorithm
on a bipartite graph with an additional set of cluster-nodes. The cluster-nodes represent cluster centroids both globally of the entire graph and locally of each neighborhood. Importantly, our message-passing steps act as optimization steps that enforce the clustering inductive bias.

We name our method Differentiable Cluster Graph Neural Network (DC-GNN) and show its convergence. In DC-GNN, both local and global cluster-nodes serve important but different roles. The local cluster-node connections help to preserve the original graph structure information and enable clustered aggregation.
Therefore, a single update via local cluster-nodes is analogous to one layer of message passing in conventional GNNs , with the cluster structure assisting in handling heterophilic local neighbourhoods.
Simultaneously, each update via global cluster-nodes allows transfer of long-range information from relevant distant nodes. Additionally, both types of cluster-nodes can be viewed as providing new pathways among the original nodes. This helps to lower the graph’s total effective resistance, which is useful in mitigating oversquashing [].

DC-GNN is efficient and has a linear complexity with respect to the graph size. We carry out comprehensive experiments on a variety of datasets, including both heterophilous and homophilous types. The results underscore the efficacy of our method, with our approach attaining state-of-the-art performance in multiple instances.

## 2 Related Work

Prominent GNN models typically follow a message passing paradigm that iteratively aggregates information in a node’s neighborhood
[, , , , , ].
This local message passing, however, requires the stacking of multiple layers to pursue long-range information and can encounter issues such as over-smoothing [, , ] and over-squashing [, , , ].
To tackle oversquashing, existing works design graph rewiring techniques that change graph topology [, , , ].
Similar to some existing methods [], our approach of adding cluster-nodes also changes graph topology and is shown to reduce effective resistance, thereby helpful in mitigating oversquashing [].

Additionally, some graphs contain heterophilous neighborhoods, in which traditional aggregation promoting similarity among neighbors is suboptimal. [, ].
There are three types of approaches to address this, including design of high-pass filters in message passing [, , ], exploring global neighborhoods [, , , ], and use of auxiliary graph structures [, , , ].
Many of these approaches are often computationally expensive and may struggle on homophilous graphs. Our method seeks to explore global neighborhoods by introducing cluster-nodes as auxiliary structures and conduct clustered aggregation in local neighborhoods, with linear complexity.

Our work focuses on using a clustering inductive bias to enhance the supervised node classification task.
This differs from tasks like graph clustering [, ] and graph pooling [, ], which are constrained by graph typology and aim to partition a graph into substructures.
Instead, our approach primarily utilizes feature information for global clustering. Unlike previous methods, we also explore clustering within local neighborhoods, which has not been explored before. Additionally, we integrate the clustering bias into the message passing mechanism in an end-to-end differentiable framework.

Clustering with a differentiable pipeline has been explored in some works including in unsupervised and self-supervised settings [, , , ]. However, these works do not explore application on graph structured data.
Central to our approach is an Optimal Transport based clustering objective.
OT has been recently applied to graph learning for tasks like
graph classification [, ], regularizing node representations [] and finetuning [].
However, most of these methods leverage OT as a separate component and do not integrate it within message passing. Distinct from these approaches, we implicitly optimize an OT-based clustering objective function via message passing.

Our approach of message passing with additional cluster-nodes can also be viewed as the construction of factor graph neural network with the cluster-nodes representing factors, which has been effective in learning representation on structured data
[, , , , ]. Our approach provides further evidence for the usefulness of factor graph based approaches.

## 3 Methodology

In this section, we introduce our unified clustering-based GNN message passing method to address both long-range interactions and heterophilous neighborhood aggregation. This involves transforming the input graph into a bipartite graph by introducing cluster-nodes, defining an optimal transport (OT) based clustering objective function, and optimizing it through our derived message passing steps within a differentiable coordinate descent framework. Notations are introduced where needed throughout the paper. A complete list of these notations is available in Appendix A.

### 3.1 DC-GNN formulation

We begin by constructing a bipartite graph, denoted as $\mathcal{G}=(\mathcal{V},\mathcal{C},\mathcal{E})$. This bipartite graph is derived from the original graph $G=(V,E)$ and comprises two distinct sets of nodes. The first set, $\mathcal{V}$, is a direct copy of the nodes $V$ from the original graph. The second set, $\mathcal{C}$, consists of cluster-nodes divided into two categories: global clusters ($\Omega$) and local clusters ($\Gamma$). $\mathcal{E}$ and $E$ represent the set of edges in the bipartitie graph and original graph respectively.

In this bipartite graph, each global cluster-node from $\Omega$ connects to all nodes in $\mathcal{V}$, thereby facilitating long-range interactions across distant nodes. Meanwhile, each local cluster-node from $\Gamma$ is associated with a specific node $i$ in $\mathcal{V}$, and connected to its ego-neighborhood. The ego-neighborhood of node $i$, denoted as $\mathcal{N}_{i}^{+}$, includes $i$ itself and its one-hop neighbors, formally defined as $\mathcal{N}_{i}^{+}=\mathcal{N}_{i}\cup\{i\}$, where $\mathcal{N}_{i}=\{j:(i,j)\in E\}$. Here, $\Gamma_{i}$ represents the set of local clusters associated with node $i$ in $\mathcal{V}$. The total number of nodes in $\mathcal{C}$ equals the nodes in $\Omega$ plus those in all local clusters, i.e., $|\mathcal{C}|=|\Omega|+\sum_{i\in\mathcal{V}}|\Gamma_{i}|$. An illustration of the bipartite graph construction is provided in Appendix C.

Figure: Figure 2: Overview of DC-GNN. DC-MsgPassing is an iterative algorithm that minimizes $\mathcal{L}_{\rm cluster}^{\lambda}$ in each step of message passing, where $\mathcal{L}_{\rm cluster}^{\lambda}$ is an optimal transport based clustering objective function.
Refer to caption: extracted/5620147/arxiv_paper/Figures/dcgnn.png

Following the bipartite graph construction, we propose our Differentiable Cluster Graph Neural Network (DC-GNN) to learn the node embeddings. An illustration of DC-GNN architecture is presented in Fig. 2.
Given $X_{\rm input}$ as input node features, we first transform it with an MLP to produce $X$. The transformed features together with the bipartitie graph $\mathcal{G}$ are then put through an iterative DC-MsgPassing algorithm to produce the embeddings of nodes $Z$ and cluster-nodes $C$. A learnable readout function such as a simple MLP, is then used on $Z$ to produce the class probabilities $Y$ in our node classification task.
Then, DC-GNN is trained end-to-end with a loss function $\mathcal{L}_{\rm train}$ via backpropagation. At the core of our design is the iterative DC-MsgPassing procedure, an implicit optimization algorithm that optimizes a clustering-based objective function $\mathcal{L}_{\rm cluster}^{\lambda}$.

### 3.2 Clustering-based objective function ℒ cluster subscript ℒ cluster \mathcal{L}_{\rm cluster} caligraphic_L start_POSTSUBSCRIPT roman_cluster end_POSTSUBSCRIPT

In this work, we aim to address over-squashing and heterophily by embedding a clustering inductive bias into the design of the GNN.
One way to achieve this is by optimizing a clustering-based objective function.
Motivated by the theory of optimal transport, we propose a novel soft clustering-based objective function.

We conceptualize the cluster assignment problem as an optimal transport problem, where the cost is defined as the distance between the embeddings of nodes and their cluster centroids. Therefore, we propose to minimize the overall cost, weighted by the soft cluster assignment matrix $P$ which indicates the amount of assignment from a node to a cluster. Specifically, we have a single global soft cluster assignment matrix $P^{\Omega}\in\mathbb{R}^{|\mathcal{V}|\times|\Omega|}$ and local soft cluster assignment matrices $P^{\Gamma_{i}}\in\mathbb{R}^{|\mathcal{N}_{i}^{+}|\times|\Gamma_{i}|}$ for each node $i$.
Let $d(u,v)$ represent the distance between two vectors $u$ and $v$, $z_{i}$ be the node embeddings to be learnt for each node $i\in\mathcal{V}$, $x_{i}$ be the node features after an initial transformation by a multilayer perceptron,
and $c_{j}^{\Omega}$, $c_{j}^{\Gamma_{i}}$ be the embeddings of the $j^{th}$ global and local cluster-node respectively.
Then we define our objective function $\mathcal{L}_{\rm cluster}$ as

$$ $\begin{split}\mathcal{L}_{\rm cluster}=\;&\alpha\underbrace{\sum_{i\in\mathcal {V}}\sum_{j\in\Omega}P_{ij}^{\Omega}d(z_{i},c_{j}^{\Omega})}_{\text{global clustering }}+(1-\alpha)\sum_{i\in\mathcal{V}}\underbrace{\sum_{u\in\mathcal{N }_{i}^{+}}\sum_{j\in\Gamma_{i}}P_{uj}^{\Gamma_{i}}d(z_{u},c_{j}^{\Gamma_{i}})} _{\text{local clustering }}+\;\beta\underbrace{\sum_{i\in\mathcal{V}}d(z_{i},x _{i})}_{\text{node fidelity }}.\end{split}$ (1) $$

The global clustering part optimizes the OT distance between node embeddings and global cluster-node embeddings. The local clustering term optimizes the OT distance between the embeddings of nodes and local cluster-nodes within each ego-neighborhood. The scalar parameter $\alpha\in[0,1]$ is a balancing factor between the two terms. Furthermore, the additional node fidelity term encourages the node embeddings to retain some information from the original node features [, ].

### 3.3 DC-MsgPassing : Optimize ℒ cluster subscript ℒ cluster \mathcal{L}_{\rm cluster} caligraphic_L start_POSTSUBSCRIPT roman_cluster end_POSTSUBSCRIPT with entropic regularization via message passing

Since conventional OT solvers can be computationally prohibitive [], we
adopt the entropy regularized version of the OT distance that is designed for efficiency and offers a good approximation of OT distance [], with details in Section 3.3.1.
let $h(P)$ be the entropy of the assignment matrix $P$, we propose the following refined objective function

$$ $\mathcal{L}_{\rm cluster}^{\lambda}\coloneqq\mathcal{L}_{\rm cluster}-\frac{ \alpha}{\lambda}h(P^{\Omega})-\frac{(1-\alpha)}{\lambda}\sum_{i\in\mathcal{V}} h(P^{\Gamma_{i}}).$ (2) $$

Direct optimization of $\mathcal{L}_{\rm cluster}^{\lambda}$ as a training objective is intractable due to the presence of unobserved assignment matrices $P^{\Omega}$ and $P^{\Gamma_{i}}$, since we cannot compute Eq. (2) without estimating the cluster assignment values.
To overcome this, we propose an iterative block coordinate descent algorithm called DC-MsgPassing.
In each iteration of DC-MsgPassing, we alternatively optimize $\mathcal{L}_{\rm cluster}^{\lambda}$ with respect to one block of variables at a time, while all other variables are held constant.
Specifically, we alternatively update the assignment matrices $P^{\Omega}$, $P^{\Gamma_{i}}$ and embeddings $Z,C$ in each iteration.

#### 3.3.1 Assignment update

##### Update for assignment matrices P Ω , P Γ i superscript 𝑃 Ω superscript 𝑃 subscript Γ 𝑖 P^{\Omega},P^{\Gamma_{i}} italic_P start_POSTSUPERSCRIPT roman_Ω end_POSTSUPERSCRIPT , italic_P start_POSTSUPERSCRIPT roman_Γ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUPERSCRIPT :

With node and cluster-node embeddings kept constant, we aim to update the $P^{\Omega},P^{\Gamma_{i}}$ minimizing $\mathcal{L}_{\rm cluster}$. We approach this cluster assignment problem from the perspective of optimal transport theory, and use $P^{\Omega}$ as an example without loss of generality.

Our objective is to determine the optimal assignment matrix $P^{\Omega*}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Omega|}$ for a given cost matrix $M\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Omega|}$ , aligning with a target clustering distribution. Assuming a uniform distribution of data points across all clusters, we aim to find a mapping from $\mathbf{u}^{\top}=\begin{bmatrix}|\mathcal{V}|^{-1}&\cdots&|\mathcal{V}|^{-1}
\end{bmatrix}_{1\times|\mathcal{V}|}$ to $\mathbf{v}^{\top}=\begin{bmatrix}|\Omega|^{-1}&\cdots&|\Omega|^{-1}\end{
bmatrix}_{1\times|\Omega|}$, minimizing the overall cost. To formalize, optimizing $\mathcal{L}_{\rm cluster}$ with respect to $P^{\Omega}$ is tantamount to solving:

$$ $\min_{P^{\Omega}\in U(\mathbf{u},\mathbf{v})}\langle P^{\Omega},M\rangle,$ (3) $$

where $U(\mathbf{u},\mathbf{v})=\{P\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Omega|}:P
\mathbf{1}_{|\Omega|}=\mathbf{u},P^{\top}\mathbf{1}_{|\mathcal{V}|}=\mathbf{v}\}$. $\langle\cdot,\cdot\rangle$ is the Frobenius dot-product, and $\langle P^{\Omega},M\rangle=\sum_{i\in\mathcal{V}}\sum_{j\in\Omega}P_{ij}^{
\Omega}d(z_{i},c_{j}^{\Omega})$. Importantly, $M_{ij}=d(z_{i},c^{\Omega}_{j})$ can be seen as the cost of assigning node $i$ to cluster $j$, and $P^{\Omega}_{ij}$ indicates the amount of assignment from node $i$ to cluster $j$.

This is a classical optimal transport problem, where Eq. (3) represents the optimal transport distance, also known as the earth mover’s distance.
While such problems are typically solved via linear programming techniques, these approaches are computationally expensive [].
To overcome it, we opt for the Sinkhorn distance [] instead, which offers a good approximation to the optimal transport distance with additional entropic regularization, weighted by scalar $1/\lambda$, where $\lambda>0$. Formally,

$$ $\begin{gathered}\langle P_{\lambda}^{\Omega},M\rangle,\text{ where }P_{\lambda }^{\Omega}=\operatorname*{arg\,min}_{P^{\Omega}\in U(\mathbf{u},\mathbf{v})} \langle P^{\Omega},M\rangle-\frac{1}{\lambda}h(P^{\Omega}).\end{gathered}$ (4) $$

The benefit of having this entropic regularization term $h(P^{\Omega})$ is that the solution $P_{\lambda}^{\Omega*}$ now has the form $P_{\lambda}^{\Omega*}=UBV$ [],
where $B=e^{-\lambda M}$, and $U$ and $V$ are diagonal matrices.
Now the OT problem reduces to the classical matrix scaling problem, for which the objective is to determine if there exist diagonal matrices $U$ and $V$ such that the $i^{th}$ row of the matrix $UBV$ sums to $\mathbf{u}_{i}$ and the $j^{th}$ column of $UBV$ sums to $\mathbf{v}_{j}$.
Since $e^{-\lambda M}$ is strictly positive, there exists a unique $P_{\lambda}^{\Omega*}$ that belongs to $U(\mathbf{u},\mathbf{v})$ [, ], which can be obtained by the well-known Sinkhorn–Knopp algorithm [, ].

To obtain $P_{\lambda}^{\Omega*}$, we run the Sinkhorn–Knopp algorithm which iteratively updates the matrix $B$ by scaling each row of $B$ by the respective row-sum, and each column of $B$ by the respective column-sum. Formally, we have

$$ $B_{ij}^{(t)}=\frac{B^{(t-1)}_{ij}}{\sum_{j}B^{(t-1)}_{ij}}\mathbf{u}_{i},\quad B _{ij}^{(t+1)}=\frac{B^{(t)}_{ij}}{\sum_{i}B^{(t)}_{ij}}\mathbf{v}_{j}.$ (5) $$

After $T$ steps, we update $P^{\Omega}_{ij}$ by the value of $B_{ij}$. Similarly, we update $P^{\Gamma_{i}}$ in the local clustering term. The scaling operations in Sinkhorn–Knopp are fully differentiable, enabling end-to-end learning.

#### 3.3.2 Embeddings update

With the updated assignment matrices $P^{\Omega},P^{\Gamma_{i}}$, we now derive the message passing equations to update the cluster-node and node embeddings.

##### Update for cluster-node embeddings C 𝐶 C italic_C :

To minimize $\mathcal{L}_{\rm cluster}^{\lambda}$with respect to a specific cluster node $c_{j}$, we differentiate $\mathcal{L}_{\rm cluster}^{\lambda}$  in terms of $c_{j}$ with other variables fixed and set the derivative to zero.
If we choose the distance function $d(\cdot)$ to be the squared Euclidean norm, $\mathcal{L}_{\rm cluster}^{\lambda}$ is a quadratic model of node embeddings $c_{j}$. Then we can derive the following closed-form solution:

$$ $C=\text{diag}(\mathbf{k})P^{\top}Z,$ (6) $$

where $\mathbf{k^{\top}}=\begin{bmatrix}|\Omega|\mathbf{1}_{|\Omega|};|\Gamma_{i}|
\mathbf{1}_{|\Gamma_{i}|};\cdots;|\Gamma_{|\mathcal{V}|}|\mathbf{1}_{|\Gamma_{
i}|}\end{bmatrix}_{1\times|\mathcal{C}|}$, ; denotes concatenation. $P\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\mathcal{C}|}$ is the overall assignment matrix that indicates the amount of assignment from all nodes $Z$ to all cluster-nodes $C$. See Appendix B.2 and B.3 for notation and full derivation.

###### Remark 3.1 .

Since $P_{ij}$ can be viewed as edge weight between a node $i\in\mathcal{V}$ and a cluster-node $j\in\mathcal{C}$,
the updating mechanism serves as the message passing function from nodes to cluster-nodes.

###### Remark 3.1 .

##### Update for node embeddings Z 𝑍 Z italic_Z :

Similar to cluster-node embeddings update, we derive the node embeddings update by differentiating $\mathcal{L}_{\rm cluster}^{\lambda}$ with respect to $z_{i}$. Consequently, we obtain

$$ $Z=\gamma\Big{[}\beta X+\alpha P^{\Omega}C^{\Omega}+(1-\alpha)\sum_{u\in \mathcal{N}_{i}^{+}}\hat{P}^{\Gamma_{u}}C^{\Gamma_{u}}\Big{]},$ (7) $$

where $\gamma=({\alpha}{|\mathcal{V}|}^{-1}+\beta+1-\alpha)^{-1}$ is a constant, and $\hat{P}^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Gamma_{u}|}$ is the broadcasted local assignment matrix from $P^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{N}_{u}^{+}|\times|\Gamma_{u}|}$ (see Appendix B.2 and B.4 for notation and full derivation).
Intriguingly, this closed-form solution reveals that $Z$ is updated by a linear combination of global and local cluster-node embeddings, weighted by cluster assignment probabilities.
The hyperparameter $\alpha$ balances the influence of local and long-range interactions while $\beta$ scales the original node features that provide the initial residual [, ].

###### Remark 3.2 .

Viewing the cluster assignment probabilities as edge weights, Eq. (7) represents the message passing from cluster-nodes back to original nodes. Therefore, Eq. (6) and Eq. (7) function as message passing on the bipartite graph, enforcing the clustering inductive bias.

Figure: Algorithm 1 DC-MsgPassing

In summary, DC-MsgPassing is the key to our method.
Each iteration of DC-MsgPassing consists of two alternative steps. First, with fixed embeddings $Z$ and $C$, optimal clustering assignment matrices are calculated via Sinkhorn-Knopp algorithm (Eq. (5)). Then, the cluster-node and node embeddings are refined through message passing with the updated assignment matrices via Eq. (6)
and Eq. (7). In practice, we could also add learnable components such as linear transformation matrices or MLPs for the messages in Eq. (6) and Eq. (7) to allow the network to fit the data distribution better. Each iteration of DC-MsgPassing is one optimization step towards minimizing $\mathcal{L}_{\rm cluster}^{\lambda}$. Thus our message passing mechanism provides the needed inductive bias in learning the local and global clusters present in the data, while simultaneously learning node embeddings using the clustering information. We present the details of DC-MsgPassing in Algorithm 1.

###### Remark 3.2 .

##### Convergence analysis

The DC-MsgPassing algorithm is generally well behaved, converging when enough iterations are performed.
The following theorem presents the convergence analysis of minimizing the objective function $\mathcal{L}^{\lambda}_{cluster}$ with the DC-MsgPassing algorithm.

###### Theorem 3.3 (Convergence of DC-MsgPassing ) .

Assuming the Sinkhorn–Knopp algorithm is run to convergence in each iteration, for any $\lambda>0$, the value of $\mathcal{L}^{\lambda}_{cluster}$ produced by DC-MsgPassing algorithm (Algorithm 1) is guaranteed to converge.

Proof can be found in the Appendix B.1.

###### Theorem 3.3 (Convergence of DC-MsgPassing ) .

##### Complexity analysis

The complexity of the global assignment update step is $O(T|\mathcal{V}||\Omega|)$ = $O(|\mathcal{V}|)$ as $|\Omega|\ll|\mathcal{V}|$ and $T$ are small constants. Both the complexity of the local assignment update step $O(T|\mathcal{E}|)$ = $O(|E|)$ and the embeddings update step $O(|\mathcal{E}|)$ = $O(|E|)$ are linear w.r.t. the number of edges. Thus the overall computational complexity is linear w.r.t. the size of the original graph $O(|E|)$, same order as standard GNNs. Runtime of DC-MsgPassing is measured in Appendix D.3.

### 3.4 Training objective function ℒ train subscript ℒ train \mathcal{L}_{\rm train} caligraphic_L start_POSTSUBSCRIPT roman_train end_POSTSUBSCRIPT

Our task is supervised node-classification. We use cross entropy loss $\mathcal{L}_{ce}$ to train DC-GNN along with two regularizing loss functions to facilitate the clustering process.

$$ $\mathcal{L}_{\rm train}=\mathcal{L}_{\rm ce}+\omega_{1}\mathcal{L}_{\rm ortho} +\omega_{2}\mathcal{L}_{\rm sim},$ (8) $$

where $\mathcal{L}_{\rm ortho}$ and $\mathcal{L}_{\rm sim}$ are orthogonality and similarity losses, weighted by hyperparameters $\omega_{1}$ and $\omega_{2}$, as described below.

##### Orthogonality loss ( ℒ ortho subscript ℒ ortho \mathcal{L}_{\rm ortho} caligraphic_L start_POSTSUBSCRIPT roman_ortho end_POSTSUBSCRIPT ):

To encourage the clusters to be distinct,
we adopt a regularizing orthogonality loss function [] $\mathcal{L}_{\rm ortho}=\bigg{\lVert}\frac{C^{\top}C}{\|C^{\top}C\|_{F}}-\frac
{I_{|\Omega|}}{\sqrt{|\Omega|}}\bigg{\rVert}_{F}$, where $\|\cdot\|_{F}$ is the Frobenius norm. This pushes the cluster-nodes to be orthogonal to each other.

##### Similarity loss ( ℒ sim subscript ℒ sim \mathcal{L}_{\rm sim} caligraphic_L start_POSTSUBSCRIPT roman_sim end_POSTSUBSCRIPT ):

To further enhance the clustering process, we introduce $\mathcal{L}_{\rm sim}$ that encourages node similarity to only a single cluster.
To achieve this, we set $|\Omega|$ to be multiple of the number of classes and associate a set of cluster-nodes $\Omega^{\tau}$ with each class $\tau$.
We then compute distances between the node embedding and the cluster-node embeddings associated with its labelled class,
select the cluster-node embedding that is most similar to the node with a $max$ operator, and push them closer.
If a training node $i$ belongs to class $\tau$, $c_{j}^{\tau}$ is the $j^{th}$ cluster embedding associated with class $\tau$. Let $\Lambda$ be a similarity function, and $\mathcal{V}_{+}$ be the set of training nodes,
then $\mathcal{L}_{\rm sim}$ is defined as

$$ $\displaystyle\mathcal{L}_{\rm sim}=\frac{1}{|\Omega||\mathcal{V}_{+}|}\text{ \Large$\sum_{i\in\mathcal{V}_{+}}$}\;\Biggl{[}s_{i}^{\tau}+\log\sum_{\tau^{ \prime}\neq\tau}\exp{\big{(}-s_{i}^{\tau^{\prime}}\big{)}}\Biggl{]},\text{ where }s_{i}^{\tau}=\max_{j\in|\Omega^{\tau}|}\;\Lambda\Bigl{(}z_{i}^{L},c_{j} ^{\tau}\Bigl{)}.$ (9) $$

## 4 Experiments

In this section, we empirically validate the capabilities of our proposed solution through extensive experiments and ablation studies.

##### Datasets.

We conduct experiments on fourteen datasets, a mix of small-scale and large-scale datasets. Eleven of them are non-homophilous, including:
(1) Roman-empire, Amazon-ratings, Minesweeper, Tolokers, Questions [];
(2) Penn94, Genius, Cornell5, Amherst41
[];
(3) Wisconsin [];
(4) a US election dataset [].
Three are homophilous citation networks: Cora, Citeseer and Pubmed []. We use the original train/validation/test splits when they exist. Otherwise we follow the splits specified in [, , ].
Descriptions and statistics of the datasets are in Appendix. D.4.

Baselines, implementation and training details
can be found in Appendix. D.5.

**Table 1: Classification performance comparison across various heterophilous and homophilous datasets. We report ROC AUC for Genius [] and accuracy for other datasets. OOM refers to out-of-memory errors.**
|  | Heterophilous | Homophilous |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | Penn94 | Genius | Cornell5 | Amherst41 | US-election | Wisconsin | Cora | Citeseer | Pubmed |
| MLP | 73.61 (0.40) | 86.68 (0.09) | 68.86 (1.83) | 60.43 (1.26) | 81.92 (1.01) | 85.29 (3.31) | 75.69 (2.00) | 74.02 (1.90) | 87.16 (0.37) |
| GCN | 82.47 (0.27) | 87.42 (0.37) | 80.15 (0.37) | 81.41 (1.70) | 82.07 (1.65) | 51.76 (3.06) | 86.98 (1.27) | 76.50 (1.36) | 88.42 (0.50) |
| GAT | 81.53 (0.55) | 55.80 (0.87) | 78.96 (1.57) | 79.33 (2.09) | 84.17 (0.98) | 49.41 (4.09) | 87.30 (1.10) | 76.55 (1.23) | 86.33 (0.48) |
| MixHop | 83.47 (0.71) | 90.58 (0.16) | 78.52 (1.22) | 76.26 (2.56) | 85.90 (1.55) | 75.88 (4.90) | 87.61 (0.85) | 76.26 (1.33) | 85.31 (0.61) |
| GCNII | 82.92 (0.59) | 90.24 (0.09) | 78.85 (0.78) | 76.02 (1.38) | 82.90 (0.29) | 80.39 (3.40) | 88.37 (1.25) | 77.33 (1.48) | 90.15 (0.43) |
| H2GCN | 81.31 (0.60) | OOM | 78.46 (0.75) | 79.64 (1.63) | 85.53 (0.77) | 87.65 (4.98) | 87.87 (1.20) | 77.11 (1.57) | 89.49 (0.38) |
| WRGAT | 74.32 (0.53) | OOM | 71.11 (0.48) | 62.59 (2.46) | 84.45 (0.56) | 86.98 (3.78) | 88.20 (2.26) | 76.81 (1.89) | 88.52 (0.92) |
| GPR-GNN | 81.38 (0.16) | 90.05 (0.31) | 73.30 (1.87) | 67.00 (1.92) | 84.49 (1.09) | 82.94 (4.21) | 87.95 (1.18) | 77.13 (1.67) | 87.54 (0.38) |
| GGCN | 73.62 (0.61) | OOM | 71.35 (0.81) | 66.53 (1.61) | 84.71 (2.60) | 86.86 (3.29) | 87.95 (1.05) | 77.14 (1.45) | 89.15 (0.37) |
| ACM-GCN | 82.52 (0.96) | 80.33 (3.91) | 78.17 (1.42) | 70.11 (2.10) | 85.14 (1.33) | 88.43 (3.22) | 87.91 (0.95) | 77.32 (1.70) | 90.00 (0.52) |
| LINKX | 84.71 (0.52) | 90.77 (0.27) | 83.46 (0.61) | 81.73 (1.94) | 84.08 (0.67) | 75.49 (5.72) | 84.64 (1.13) | 73.19 (0.99) | 87.86 (0.77) |
| GloGNN++ | 85.74 (0.42) | 90.91 (0.13) | 83.96 (0.46) | 81.81 (1.50) | 85.48 (1.19) | 88.04 (3.22) | 88.33 (1.09) | 77.22 (1.78) | 89.24 (0.39) |
| DC-GNN | 86.69 (0.22) | 91.70 (0.08) | 84.68 (0.24) | 82.94 (1.59) | 89.59 (1.60) | 91.67 (1.95) | 89.13 (1.18) | 77.93 (1.82) | 91.00 (1.28) |

**Table 2: Classification performance comparison on more recent heterophilous datasets []. Following [], we report accuracy for Roman-empire and Amazon-ratings, and ROC AUC for the rest.**
|  | Roman-empire | Amazon-ratings | Minesweeper | Tolokers | Questions |
| --- | --- | --- | --- | --- | --- |
| ResNet | 65.88 (0.38) | 45.90 (0.52) | 50.89 (1.39) | 72.95 (1.06) | 70.34 (0.76) |
| ResNet+SGC | 73.90 (0.51) | 50.66 (0.48) | 70.88 (0.90) | 80.70 (0.97) | 75.81 (0.96) |
| ResNet+adj | 52.25 (0.40) | 51.83 (0.57) | 50.42 (0.83) | 78.78 (1.11) | 75.77 (1.24) |
| GCN | 73.69 (0.74) | 48.70 (0.63) | 89.75 (0.52) | 83.64 (0.67) | 76.09 (1.27) |
| SAGE | 85.74 (0.67) | 53.63 (0.39) | 93.51 (0.57) | 82.43 (0.44) | 76.44 (0.62) |
| GAT | 80.87 (0.30) | 49.09 (0.63) | 92.01 (0.68) | 83.70 (0.47) | 77.43 (1.20) |
| GAT-sep | 88.75 (0.41) | 52.70 (0.62) | 93.91 (0.35) | 83.78 (0.43) | 76.79 (0.71) |
| GT | 86.51 (0.73) | 51.17 (0.66) | 91.85 (0.76) | 83.23 (0.64) | 77.95 (0.68) |
| GT-sep | 87.32 (0.39) | 52.18 (0.80) | 92.29 (0.47) | 82.52 (0.92) | 78.05 (0.93) |
| GPS^GAT+Performer | 87.04 (0.58) | 49.92 (0.68) | 91.08 (0.58) | 84.38 (0.91) | 77.14 (1.49) |
| H2GCN | 60.11 (0.52) | 36.47 (0.23) | 89.71 (0.31) | 73.35 (1.01) | 63.59 (1.46) |
| CPGNN | 63.96 (0.62) | 39.79 (0.77) | 52.03 (5.46) | 73.36 (1.01) | 65.96 (1.95) |
| GPR-GNN | 64.85 (0.27) | 44.88 (0.34) | 86.24 (0.61) | 72.94 (0.97) | 55.48 (0.91) |
| FSGNN | 79.92 (0.56) | 52.74 (0.83) | 90.08 (0.70) | 82.76 (0.61) | 78.86 (0.92) |
| GloGNN | 59.63 (0.69) | 36.89 (0.14) | 51.08 (1.23) | 73.39 (1.17) | 65.74 (1.19) |
| FAGCN | 65.22 (0.56) | 44.12 (0.30) | 88.17 (0.73) | 77.75 (1.05) | 77.24 (1.26) |
| GBK-GNN | 74.57 (0.47) | 45.98 (0.71) | 90.85 (0.58) | 81.01 (0.67) | 74.47 (0.86) |
| JacobiConv | 71.14 (0.42) | 43.55 (0.48) | 89.66 (0.40) | 68.66 (0.65) | 73.88 (1.16) |
| DC-GNN | 89.96 (0.35) | 51.11 (0.47) | 98.50 (0.21) | 85.88 (0.81) | 78.96 (0.60) |

### 4.1 Comparison with baselines on heterophilous graphs

We first conduct experiments on heterophilous datasets where long-range information is beneficial and neighborhood aggregation needs special attention. We achieve state-of-the-art on all six heterophilous datasets in Tab. 1.
The performance uplift is especially pronounced on US-election and Wisconsin, where our method outperforms baselines by more than 3%.
Furthermore, our method achieves state-of-the-art on four out of five heterophilous datasets proposed by [], as shown in Tab. 2. Our performance is especially strong on Minesweeper, where we surpass existing baselines by 5%.
Baseline results are from [], [] and [] except for Cornell5, Amherst41 and US-election, which we reproduced following the code in []. Details are in Appendix. D.5.

### 4.2 Comparison with baselines on homophilous graphs

For a more comprehensive evaluation, we also run DC-GNN on three well-known homophilous citation network datasets Cora, Citeseer and Pubmed [].
As shown in Tab. 1, DC-GNN achieves the best performance on all three datasets, outperforming both general-purpose GNN baselines and those proposed specifically for heterophilous graphs. Our strong performance on both homophilous and heterophilous graphs shows that our network with the built-in clustering inductive bias is flexible and adaptive, capable of effective information aggregation in graphs with various homophily levels.

### 4.3 Benefit of capturing long range information in sparse label settings

Our method introduces shortcuts between distant nodes via the cluster-nodes.
To empirically evaluate the effects of shortcut construction, we consider a generalized scenario on homophilous graphs where information from labeled nodes needs to propagate over a long distance to reach most unlabeled nodes.
Specifically, we evaluate our method with only a small number of training labels per class.

**Table 3: Classification accuracy with 5 labels per class.**
|  | Cora | Citeseer | Pubmed |
| --- | --- | --- | --- |
| GCN | 69.23 (3.39) | 63.03 (4.48) | 68.00 (3.75) |
| DC-GNN | 72.17 (1.76) | 62.14 (2.65) | 75.07 (0.82) |

Our hypothesis is that information propagation becomes more challenging when useful training information is more scarce, making the effects of shortcut construction more pronounced in sparsely-annotated graph datasets. The hypothesis is supported by experiment results in Tab. 3, where our method outperforms vanilla GCN on Cora and Pubmed with substantial improvements.

### 4.4 Alleviating oversquashing

Figure: (a) Amherst41
Refer to caption: extracted/5620147/arxiv_paper/Figures/Amherst_16_v4.png

Total effective resistance ($R_{\text{tot}}$) is established as an indicator of oversquashing [], prompting the development of various graph rewiring strategies to diminish $R_{\text{tot}}$ within the underlying graph and thus address oversquashing. Our approach contributes to this endeavor by introducing cluster-nodes. This effectively creates new pathways among the original nodes, thereby reducing the graph’s $R_{\text{tot}}$ and aiding in mitigating the oversquashing issue [].
To validate this, we conduct an empirical analysis of the total pairwise effective resistance among the original nodes in our bipartite graph, with varying number of global and local clusters. Fig. 3(a) and Fig. 3(b) display a heatmap of $R_{\text{tot}}$, with darker shades representing higher $R_{\text{tot}}$ values. The results indicate that $R_{\text{tot}}$ decreases sharply as we increase the number of both global and local clusters.
In addition, we also perform experiments on synthetic random graphs with varying degrees of sparsity and show that the bipartite graph construction reduces $R_{\text{tot}}$ on these graphs. Details can be found in Appendix 5.

To further validate the oversquashing mitigation capability, we conduct experiments on Tree-NeighborsMatch dataset [], which requires long-range interaction between leaf nodes and the root node of tree graphs with varying depths. In Fig. 3(c), DC-GNN achieves perfect performance along with GT [] on all depth settings, significantly outperforming other message passing GNNs.

### 4.5 Ablation studies

#### 4.5.1 Effects of each term in ℒ cluster λ superscript subscript ℒ cluster 𝜆 \mathcal{L}_{\rm cluster}^{\lambda} caligraphic_L start_POSTSUBSCRIPT roman_cluster end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_λ end_POSTSUPERSCRIPT

To validate the effectiveness of our DC-MsgPassing algorithm, we conduct an ablation study on the individual components of our objective function $\mathcal{L}_{\rm cluster}^{\lambda}$. Specifically, we vary the parameters by setting (1) $\alpha$ to 0, (2) $\alpha$ to 1 and (3) $\beta$ to 0, aiming to ablate the contributions of the global clustering term, local clustering term and the node fidelity term respectively.

**Table 4: Effects of each term in $\mathcal{L}_{\rm cluster}^{\lambda}$.**
|  | Genius | US-election | Penn94 | Amherst41 |
| --- | --- | --- | --- | --- |
| DC-GNN | 91.70 (0.08) | 89.59 (1.60) | 86.69(0.22) | 82.94 (1.59) |
| (-)Global | 91.62 (0.07) | 88.77 (2.21) | 84.61 (0.42) | 81.43 (1.53) |
| (-)Local | 87.05 (0.09) | 83.26 (1.77) | 86.69 (0.22) | 80.77 (2.04) |
| (-)Fidelity | 91.08 (0.04) | 87.84 (2.67) | 86.69 (0.22) | 82.28 (1.32) |

Results from Tab. 4 indicate that the contributions of global and local clustering vary on different datasets. Specifically, the contribution of local clustering is dominant on Genius, US-election and Amherst41. This is expected as local clustering facilitates message passing via adjacent nodes and embeds graph structure information into the model. The contribution of global clustering is most pronounced on Penn94, indicating the usefulness of long-range information in Penn94 captured by global clustering term. Additionally, the results show that all three terms—local clustering, global clustering, and node fidelity—contribute to the overall efficacy of our model.

#### 4.5.2 Effects of ℒ ortho subscript ℒ ortho \mathcal{L}_{\rm ortho} caligraphic_L start_POSTSUBSCRIPT roman_ortho end_POSTSUBSCRIPT and ℒ sim subscript ℒ sim \mathcal{L}_{\rm sim} caligraphic_L start_POSTSUBSCRIPT roman_sim end_POSTSUBSCRIPT

**Table 5: Effects of $\mathcal{L}_{\rm ortho}$ and $\mathcal{L}_{\rm sim}$.**
|  | Genius | US-election | Penn94 | Amherst41 |
| --- | --- | --- | --- | --- |
| DC-GNN | 91.70 (0.08) | 89.59 (1.60) | 86.69 (0.22) | 82.94 (1.59) |
| (-) $\mathcal{L}_{\rm sim}$ | 91.70 (0.08) | 89.08 (1.46) | 86.65 (0.15) | 82.22 (1.03) |
| (-) $\mathcal{L}_{\rm ortho}$ | 91.68 (0.08) | 88.72 (1.20) | 86.64 (0.27) | 82.35 (1.25) |
| (-)$\mathcal{L}_{\rm sim}$, $\mathcal{L}_{\rm ortho}$ | 91.68 (0.08) | 88.61 (1.57) | 86.40 (0.25) | 81.75 (1.07) |

We introduced orthogonality ($\mathcal{L}_{\rm ortho}$) and similarity ($\mathcal{L}_{\rm sim}$) losses to regularize and assist the clustering process. In this ablation, we evaluate the effects of these losses on model performance. As shown in Tab. 5, $\mathcal{L}_{\rm ortho}$ and $\mathcal{L}_{\rm sim}$ generally help to improve the scores across datasets. The results provide some evidence that the two losses indeed facilitate the clustering process, plausibly by encouraging distinct cluster representations and node-cluster alignment as hypothesized.

## 5 Conclusion

This paper tackles the dual challenges in Graph Neural Networks: capturing global long-range information and preserving performance in heterophilous local neighborhoods. We proposed a novel differentiable framework that seamlessly embeds a clustering inductive bias into the message passing mechanism, facilitated by the introduction of cluster-nodes.
At the heart of our approach is an optimal transport based implicit clustering objective function whose optimization presents a considerable challenge. We addressed this through an iterative optimization strategy, alternating between the computation of cluster assignments and the refinement of node/cluster-node embeddings. Importantly, the derived optimization steps effectively function as message passing steps on the bipartite graph.
The message passing algorithm is efficient and we show that it is guaranteed to converge.
The efficacy of our clustering-centric method in capturing both local nuances and global structures within graphs is supported by extensive experiments on both heterophilous and homophilous datasets.

## Appendix A Notations

All notations are listed in Tab. 6.

**Table 6: Table for notations.**
| Variable | Definition |
| --- | --- |
| $\mathcal{G}$ | bipartite graph, denoted as ($\mathcal{V}$, $\mathcal{C}$, $\mathcal{E}$) |
| $G$ | original graph, denoted as ($V$, $E$) |
| $V$ | set of nodes from the original graph $G$ |
| $\mathcal{V}$ | vertices of $\mathcal{G}$, direct copy of $V$ |
| $E$ | set of edges from the original graph $G$ |
| $\mathcal{E}$ | set of edges in the bipartite graph $\mathcal{G}$ |
| $\mathcal{N}_{i}$ | set of one-hop neighbors of node $i$ |
| $\mathcal{N}_{i}^{+}$ | node $i$ and its one-hop neighbors (ego-neighborhood of $i$) |
| $\mathcal{C}$ | set of cluster-nodes in the bipartite graph $\mathcal{G}$ |
| $\Omega$ | set of global cluster-nodes |
| $\Gamma$ | set of local cluster-nodes |
| $\Gamma_{i}$ | set of local cluster-nodes associated with $\mathcal{N}_{i}^{+}$ |
| $C$ | set of cluster-node embeddings |
| $Z$ | set of node embeddings |
| $Y$ | predicted class probabilities |
| $X_{\rm input}$ | input features |
| $X$ | transformed input features |
| $P$ | overall cluster assignment matrix. $P\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\mathcal{C}|}$ |
| $P^{\Omega}$ | global cluster assignment matrix. $P^{\Omega}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Omega|}$ |
| $P^{\Gamma}$ | local cluster assignment matrix. $P^{\Gamma_{i}}\in\mathbb{R}_{+}^{|\mathcal{N}_{i}^{+}|\times|\Gamma_{i}|}$ for each node $i$ |
| $d(\cdot)$ | distance function |
| $z_{i}$ | node embeddings of node $i$ |
| $x_{i}$ | node features of node $i$ after initial transformation |
| $c_{j}^{\Omega}$ | node embeddings of $j^{th}$ global cluster-node |
| $c_{j}^{\Gamma_{i}}$ | node embeddings of $j^{th}$ local cluster-node for $\mathcal{N}_{i}^{+}$ |
| $\alpha$ | scalar parameter that balances global and local clustering objectives. $\alpha\in[0,1]$ |
| $\beta$ | scalar parameter for node fidelity term |
| $h(\cdot)$ | entropy function |
| $P^{\Omega*}$ | optimal global soft-assignment matrix. $P^{\Omega*}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Omega|}$ |
| $M_{ij}$ | cost of assigning node $i$ to cluster $j$ |
| $M$ | cost matrix |
| $\mathbf{u}^{\top}$ | $\mathbf{u}^{\top}=\begin{bmatrix}|\mathcal{V}|^{-1}&\cdots&|\mathcal{V}|^{-1}<br>\end{bmatrix}_{1\times|\mathcal{V}|}$ |
| $\mathbf{v}^{\top}$ | $\mathbf{v}^{\top}=\begin{bmatrix}|\Omega|^{-1}&\cdots&|\Omega|^{-1}\end{<br>bmatrix}_{1\times|\Omega|}$ |
| $U(\mathbf{u},\mathbf{v})$ | $U(\mathbf{u},\mathbf{v})=\{P\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Omega|}:P<br>\mathbf{1}_{|\Omega|}=\mathbf{u},P^{\top}\mathbf{1}_{|\mathcal{V}|}=\mathbf{v}\}$ |
| $\langle\cdot,\cdot\rangle$ | Frobenius dot product |
| $P^{\Omega}_{ij}$ | the amount of assignment from node $i\in\mathcal{V}$ to cluster $j\in\Omega$ |
| $\lambda$ | scalar for entropy regularization |
| $B$ | initial value $e^{-\lambda M}$ |
| $U$ | diagonal matrix |
| $V$ | diagonal matrix |
| $T$ | number of iterations for Sinkhorn-Knopp algorithm |
| $\mathbf{k}$ | $\mathbf{k^{\top}}=\begin{bmatrix}|\Omega|\mathbf{1}_{|\Omega|};|\Gamma_{i}|<br>\mathbf{1}_{|\Gamma_{i}|};\cdots;|\Gamma_{|\mathcal{V}|}|\mathbf{1}_{|\Gamma_{<br>i}|}\end{bmatrix}_{1\times|\mathcal{C}|}$ |
| $\tau$ | class $\tau$ for our node classification task |
| $\tau^{\prime}$ | class that is not class $\tau$ |
| $\Omega^{\tau}$ | set of global cluster-nodes associated with class $\tau$ |
| $L$ | number of DC-MsgPassing iterations |
| $c_{j}^{\tau}$ | the $j^{th}$ cluster embedding associated with class $\tau$ |
| $z_{i}^{L}$ | embeddings for node $i$ after $L$ iterations |
| $\gamma$ | a constant. $\gamma=({\alpha}{N}^{-1}+\beta+1-\alpha)^{-1}$ |
| $\hat{P}^{\Gamma_{u}}$ | $\hat{P}^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Gamma_{u}|}$ broadcasted local cluster assignment matrix from $P^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{N}_{i}^{+}|\times|\Gamma_{i}|}$, defined in Eq.(16) |
| $\|\cdot\|_{F}$ | Frobenius norm |
| $\Lambda(\cdot)$ | similarity function |
| $\mathcal{V}_{+}$ | set of training nodes |

## Appendix B Proof and derivation

### B.1 Proof of Theorem 3.3

To prove Theorem 3.3, we first introduce the following Lemma.

###### Lemma B.1 .

Let $\mathcal{L}^{\lambda}_{\text{cluster}}$ be the loss function optimized by the DC-MsgPassing algorithm when the optimal transport components are replaced by the entropic regularized versions. Then, $\mathcal{L}^{\lambda}_{\text{cluster}}$ is lower-bounded by:

$$ $\mathcal{L}_{\rm cluster}^{\lambda}\geq\frac{\alpha}{\lambda}\log\frac{1}{| \mathcal{V}||\Omega|}+\frac{1-\alpha}{\lambda}\sum_{i\in\mathcal{V}}\log\frac{ 1}{|\mathcal{N}_{i}^{+}||\Gamma_{i}|}.$ $$

###### Proof.

Recall that $\Omega$ is the set of global cluster-nodes, and $\Gamma_{i}$ refers to the set of local cluster-nodes associated with a node $i$. $\mathcal{V}$ is a direct copy of the nodes $V$ from the original graph. $P^{\Omega}\in\mathbb{R}^{|\mathcal{V}|\times|\Omega|}$ and $P^{\Gamma_{i}}\in\mathbb{R}^{|\mathcal{N}_{i}^{+}|\times|\Gamma_{i}|}$ are the global and local soft cluster assignment matrices respectively. $\mathcal{N}_{i}^{+}$ refers to the ego neighborhood of node $i$.

Let $p(P^{\Omega})$ be the probability of assignment matrix $P^{\Omega}$ and $h(P^{\Omega})$ be its entropy, we can upper bound its entropy by:

$$ $\displaystyle h(P^{\Omega})={}$ $\displaystyle\mathbb{E}\left[\log\frac{1}{p(P^{\Omega})}\right]$ $\displaystyle\leq{}$ $\displaystyle\log\mathbb{E}\left[\frac{1}{p(P^{\Omega})}\right]$ (by Jensen’s inequality) $\displaystyle={}$ $\displaystyle\log\sum_{i\in\mathcal{V},j\in\Omega}p(P^{\Omega}_{ij})\frac{1}{p (P_{ij}^{\Omega})}$ (10) $\displaystyle={}$ $\displaystyle\log\big{(}|\mathcal{V}|\times|\Omega|\big{)}.$ (11) $$

The inequality is due to uniform distribution having the maximum entropy.

Similarly, for local assignment matrices $P^{\Gamma_{i}}$, we have:

$$ $h(P^{\Gamma_{i}})\leq\log\big{(}|\mathcal{N}_{i}^{+}|\times|\Gamma_{i}|\big{)}.$ (12) $$

Note that distance $d(\cdot,\cdot)\geq 0$, with the above results we obtain:

$$ $\displaystyle\mathcal{L}_{\rm cluster}^{\lambda}={}$ $\displaystyle\alpha\left(\sum_{i\in\mathcal{V}}\sum_{j\in\Omega}P_{ij}^{\Omega }d(z_{i},c_{j}^{\Omega})-\frac{1}{\lambda}h(P^{\Omega})\right)+\beta\sum_{i\in \mathcal{V}}d(z_{i},x_{i})$ (13) $\displaystyle+(1-\alpha)\sum_{i\in\mathcal{V}}\left(\sum_{k\in\mathcal{N}_{i}^ {+}}\sum_{j\in\Gamma_{i}}P_{kj}^{\Gamma_{i}}d(z_{k},c_{j}^{\Gamma_{i}})-\frac{ 1}{\lambda}h(P^{\Gamma_{i}})\right)$ (14) $\displaystyle\geq{}$ $\displaystyle-\frac{\alpha}{\lambda}h(P^{\Omega})-\frac{(1-\alpha)}{\lambda} \sum_{i\in\mathcal{V}}h(P^{\Gamma_{i}})$ (by $d(\cdot,\cdot)\geq 0$ ) $\displaystyle\geq{}$ $\displaystyle-\frac{\alpha}{\lambda}\log\big{(}|\mathcal{V}|\times|\Omega|\big {)}-\frac{(1-\alpha)}{\lambda}\sum_{i\in\mathcal{V}}\log\big{(}|\mathcal{N}_{i }^{+}|\times|\Gamma_{i}|\big{)}$ (by ( 11 ) and ( 12 )) $\displaystyle={}$ $\displaystyle\frac{\alpha}{\lambda}\log\frac{1}{|\mathcal{V}|\times|\Omega|}+ \frac{(1-\alpha)}{\lambda}\sum_{i\in\mathcal{V}}\log\frac{1}{|\mathcal{N}_{i}^ {+}|\times|\Gamma_{i}|},$ (15) $$

which concludes the proof.
∎

By Lemma B.1, there exists a lower bound of $\mathcal{L}^{\lambda}_{cluster}$.
Therefore, to prove the convergence of our algorithm, we only need to show that the loss function is guaranteed to decrease monotonically in each iteration until convergence for the assignment update step and for the embeddings update step.

For the assignment update step, let $P^{\Omega}$ be the current assignment from the previous iteration and $P^{\Omega*}$ be the new assignment obtained as

$$ $P^{\Omega*}\in\operatorname*{arg\,min}_{P^{\Omega}\in U(\mathbf{u},\mathbf{v}) }\langle P^{\Omega*},M\rangle-\frac{1}{\lambda}h(P^{\Omega*}).$ $$

The change in the loss function after this assignment step is then given by

$$ $\mathcal{L}_{\rm cluster}^{\lambda}(P^{\Omega*})-\mathcal{L}_{\rm cluster}^{ \lambda}(P^{\Omega})\leq 0,$ $$

where the inequality holds by the convergence of Sinkhorn–Knopp algorithm [, , , , ]. Similarly, we have

$$ $\mathcal{L}_{\rm cluster}^{\lambda}(P^{\Gamma_{i}*})-\mathcal{L}_{\rm cluster} ^{\lambda}(P^{\Gamma_{i}})\leq 0.$ $$

For the embeddings update step, let $Z$ and $C$ be the current embeddings from the previous iteration, and $Z^{*}$ and $C^{*}$ be the new embeddings obtained by Eq. (6) and Eq. (7). Let the $d(\cdot,\cdot)$ be the squared Euclidean distance, then since $P$ is a positive constant in this step, the loss function is convex. Since Eq. (6) and Eq. (7) are closed form solutions to the loss function, we have

$$ $\mathcal{L}_{\rm cluster}^{\lambda}(Z^{*},C^{*})-\mathcal{L}_{\rm cluster}^{ \lambda}(Z,C)\leq 0.$ $$

Since $\mathcal{L}^{\lambda}_{\rm cluster}$ has a lower bound and it decreases monotonically in each iteration, the value of $\mathcal{L}^{\lambda}_{\rm cluster}$ produced by the message passing algorithm is guaranteed to converge.

###### Lemma B.1 .

###### Proof.

### B.2 Broadcasted assignment matrices

Let $f\colon\mathcal{Z}\times\mathcal{Z}\to\mathcal{Z}$ be a mapping from the index $k$ of a node in the ego-neighborhood of the node $u$ to the index $i$ of the same node in the node set $\mathcal{V}$. Let $P^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{N}_{u}^{+}|\times|\Gamma_{u}|}$ be the local assignment matrix for the ego-neighborhood of node $u\in\mathcal{V}$.
For any $u\in\mathcal{V}$, we define the broadcasted local assignment matrix $\hat{P}^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Gamma_{u}|}$ as

$$ $\hat{P}^{\Gamma_{u}}_{ij}=\begin{cases}P^{\Gamma_{u}}_{kj},&\text{if }i=f(u,k) \\ 0,&\text{otherwise }\end{cases}.$ (16) $$

Then, we can define the the overall assignment matrix $P\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\mathcal{C}|}$, where $|\mathcal{C}|=|\Omega|+|\Gamma|$ and $|\Gamma|=\sum_{i\in\mathcal{V}}|\Gamma_{i}|$, as

$$ $P=\left[P^{\Omega}\quad\hat{P}^{\Gamma_{1}}\quad\cdots\quad\hat{P}^{\Gamma_{| \mathcal{V}|}}\right].$ (17) $$

Note that each element $P_{ij}$ can be viewed as the edge weights of a node $i\in\mathcal{V}$ and a specific cluster-node $j\in\mathcal{C}$. In simple words, P includes all the global and local cluster assignment matrices, collated in a single matrix. This allows us to unify the message passing update in Eq.( 7) for both local and global clustering terms.

### B.3 Update for cluster-node embeddings C 𝐶 C italic_C : derivation of Eq. ( 6 )

Without loss of generality, we provide the full derivation for global cluster-node embeddings update function.
With squared Euclidean norm as the distance function, i.e., $d(u,v)=\|u-v\|^{2}$, we can derive that

$$ $\frac{\partial\mathcal{L}_{\rm cluster}^{\lambda}}{\partial c^{\Omega}_{j}}= \frac{\partial}{\partial c^{\Omega}_{j}}\sum_{i\in\mathcal{V}}P^{\Omega}_{ij} \|z_{i}-c_{j}^{\Omega}\|^{2}=0$ $$

Then we have

$$ $2\sum_{i\in\mathcal{V}}P^{\Omega}_{ij}(c_{j}^{\Omega}-z_{i})=0$ $$

By rearranging the terms

$$ $\sum_{i\in\mathcal{V}}P^{\Omega}_{ij}c_{j}^{\Omega}=\sum_{i\in\mathcal{V}}P^{ \Omega}_{ij}z_{i}$ $$

Then one has

$$ $c_{j}^{\Omega}=\frac{\sum_{i\in\mathcal{V}}P^{\Omega}_{ij}z_{i}}{\sum_{i\in \mathcal{V}}P^{\Omega}_{ij}}$ $$

Since $P^{\Omega}\in U(\mathbf{u},\mathbf{v})$, we have $\sum_{i}P^{\Omega}_{ij}=\frac{1}{|\Omega|}$ where $|\Omega|$ is the number of global clusters. Therefore,

$$ $c_{j}^{\Omega}=|\Omega|\sum_{i\in\mathcal{V}}P^{\Omega}_{ij}z_{i}$ (18) $$

Similarly, with $|\Gamma_{i}|$ denoting the number of local clusters within the ego-neighborhood of node $i$, we have

$$ $c_{j}^{\Gamma_{i}}=|\Gamma_{i}|\sum_{u\in\mathcal{N}_{i}^{+}}P^{\Gamma_{i}}_{ ij}z_{u}$ (19) $$

Let $\mathbf{k^{\top}}=\begin{bmatrix}|\Omega|\mathbf{1}_{|\Omega|};|\Gamma_{i}|
\mathbf{1}_{|\Gamma_{i}|};\cdots;|\Gamma_{|\mathcal{V}|}|\mathbf{1}_{|\Gamma_{
i}|}\end{bmatrix}_{1\times|\mathcal{C}|}$, where ; denotes concatenation.
Then we can combine Eq. (18) and Eq. (19) together in one matrix equation,

$$ $C=\text{diag}(\mathbf{k})P^{\top}Z,$ (20) $$

where the overall assignment matrix $P$ is defined in Eq. (17),

### B.4 Update for node embeddings Z 𝑍 Z italic_Z : derivation of Eq. ( 7 )

We provide the full derivation of the node embeddings update function. This manifests as message passing from cluster-nodes to nodes. With squared Euclidean norm as the distance function, i.e., $d(u,v)=\|u-v\|^{2}$, we can derive that

$$ $\frac{\partial\mathcal{L}_{\rm cluster}^{\lambda}}{\partial z_{i}}=\alpha\sum_ {j\in\Omega}2P_{ij}^{\Omega}(z_{i}-c_{j})+2\beta(z_{i}-x_{i})+(1-\alpha)\sum_{ u\in\mathcal{N}_{i}^{+}}\sum_{j\in\Gamma_{u}}2P_{ij}^{\Gamma_{u}}(z_{i}-c_{j})=0$ (21) $$

From Eq. (21), we obtain

$$ $\alpha\sum_{j\in\Omega}P_{ij}^{\Omega}z_{i}-\alpha\sum_{j\in\Omega}P_{ij}^{ \Omega}c_{j}+\beta z_{i}-\beta x_{i}+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+}} \sum_{j\in\Gamma_{u}}P_{ij}^{\Gamma_{u}}z_{i}-(1-\alpha)\sum_{u\in\mathcal{N}_ {i}^{+}}\sum_{j\in\Gamma_{u}}P_{ij}^{\Gamma_{u}}c_{j}=0$ $$

By rearranging the terms, we have

$$ $(\alpha\sum_{j\in\Omega}P_{ij}^{\Omega}+\beta+(1-\alpha)\sum_{u\in\mathcal{N}_ {i}^{+}}\sum_{j\in\Gamma_{u}}P_{ij}^{\Gamma_{u}})z_{i}=\alpha\sum_{j\in\Omega} P_{ij}^{\Omega}c_{j}+\beta x_{i}+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+}}\sum_ {j\in\Gamma_{u}}P_{ij}^{\Gamma_{u}}c_{j}$ $$

Since $P^{\Omega}\in U(\mathbf{u},\mathbf{v})$, we have $\sum_{j}P^{\Omega}_{ij}=\frac{1}{|\mathcal{V}|}$. Similarly, $\sum_{j\in\Gamma_{u}}P_{ij}^{\Gamma_{u}}=\frac{1}{|\mathcal{N}_{i}^{+}|}$. Therefore, we can deduce that

$$ $(\frac{\alpha}{|\mathcal{V}|}+\beta+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+}} \frac{1}{|\mathcal{N}_{i}^{+}|})z_{i}=\alpha\sum_{j\in\Omega}P_{ij}^{\Omega}c_ {j}+\beta x_{i}+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+}}\sum_{j\in\Gamma_{u}}P _{ij}^{\Gamma_{u}}c_{j}$ $$

With $\sum_{u\in\mathcal{N}_{i}^{+}}\frac{1}{|\mathcal{N}_{i}^{+}|}=1$, we have

$$ $(\frac{\alpha}{|\mathcal{V}|}+\beta+1-\alpha)z_{i}=\alpha\sum_{j\in\Omega}P_{ ij}^{\Omega}c_{j}+\beta x_{i}+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+}}\sum_{j \in\Gamma_{u}}P_{ij}^{\Gamma_{u}}c_{j}$ $$

This leads to

$$ $z_{i}=\frac{1}{\frac{\alpha}{|\mathcal{V}|}+\beta+1-\alpha}[\alpha\sum_{j\in \Omega}P_{ij}^{\Omega}c_{j}+\beta x_{i}+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+ }}\sum_{j\in\Gamma_{u}}P_{ij}^{\Gamma_{u}}c_{j}]$ $$

Finally, expressing in matrix form admits the following closed-form solution

$$ $Z=\frac{1}{\frac{\alpha}{|\mathcal{V}|}+\beta+1-\alpha}[\alpha P^{\Omega}C^{ \Omega}+\beta X+(1-\alpha)\sum_{u\in\mathcal{N}_{i}^{+}}\hat{P}^{\Gamma_{u}}C^ {\Gamma_{u}}]$ $$

where $\hat{P}^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{V}|\times|\Gamma_{u}|}$ is the broadcasted local cluster assignment matrix from $\hat{P}^{\Gamma_{u}}\in\mathbb{R}_{+}^{|\mathcal{N}_{u}^{+}|\times|\Gamma_{u}|}$ as defined in Eq.(16).

## Appendix C Illustration of bipartite graph formulation

We construct a bipartite graph, denoted as $\mathcal{G}=(\mathcal{V},\mathcal{C},\mathcal{E})$. The graph is derived from the original graph $G=(V,E)$ and comprises two distinct sets of nodes. The first set, $\mathcal{V}$, is a direct copy of the nodes $V$ from the original graph. The second set, $\mathcal{C}$, consists of cluster nodes divided into two categories: global clusters ($\Omega$) and local clusters ($\Gamma$).

In this bipartite graph, each node from the global clusters $\Omega$ connects to all nodes in $\mathcal{V}$. Meanwhile, each node from the local clusters $\Gamma$ is associated with a specific node $i$ in $\mathcal{V}$, and connected to its ego-neighborhood, which includes node $i$ and its one-hop neighbors. For a node $i$ in $\mathcal{V}$, $\Gamma_{i}$ represents the set of local clusters associated with it. The total number of nodes in $\mathcal{C}$ is the sum of nodes in $\Omega$ and the nodes in all local clusters i.e., $|\mathcal{C}|=|\Omega|+\sum_{i\in\mathcal{V}}|\Gamma_{i}|$. An illustration of the bipartite graph is provided in Fig. 4.

Figure: Figure 4: Based on the original graph on the left, we construct a bipartite graph on the right by adding local and global cluster-nodes. For each node in the original graph, a set of local cluster-nodes, represented by the blue boxes at the top, is connected to its ego-neighborhood. For example, the ego-neighborhood of node a includes itself and its one-hop neighbor node b. Therefore the local cluster-nodes for node a are connected to a and b. Meanwhile, a set of global cluster-nodes are added and connected to all nodes in the original graph, as represented by the blue boxes at the bottom.
Refer to caption: extracted/5620147/arxiv_paper/Figures/bipartite.png

## Appendix D Experiments

### D.1 Effective resistance on random graphs

Figure: (a)
Refer to caption: extracted/5620147/arxiv_paper/Figures/Connected_erods_num_10_0.2.png

Let $u$ and $v$ be vertices of $G$. The effective resistance between $u$ and $v$ is defined as

$$ $R_{u,v}=(1_{u}-1_{v})^{T}L^{+}(1_{u}-1_{v}),$ $$

where $1_{v}$ is the indicator vector of the vertex $v$ []. Let $A$ be the adjacency matrix and $D$ be the degree matrix. The Laplacian is $L=D-A$ and $L^{+}$ is the pseudoinverse of $L$.
The total effective resistance ($R_{\text{tot}}$) of a graph is therefore the total sum of effective resistance between every pair of nodes.

We measure effective resistance ($R_{\text{tot}}$) in synthetic random graphs with different degrees of sparsity. Results in Fig. 5 show that both global and local cluster-nodes contribute to reducing effective resistance, as demonstrated by decreasing $R_{\text{tot}}$ values in both row and column directions. Additionally, the more drastic $R_{\text{tot}}$ decrease from the first to last column in heatmap (a) compared to heatmap (d) show that global cluster-nodes play a more pronounced role in reducing effective resistance at a higher edge sparsity setting.

### D.2 More ablation studies

#### D.2.1 Aggregation operation in similarity loss

**Table 7: Effects of aggregator function in $\mathcal{L}_{\rm sim}$.**
| AGG | Penn94 | Amherst41 |
| --- | --- | --- |
| MEAN | 86.13 (0.12) | 81.23 (1.34) |
| SUM | 85.84 (0.26) | 81.26 (1.58) |
| MAX | 86.69 (0.22) | 82.94 (1.59) |

After computing the similarity between a node and each of the multiple clusters from the same class, the choice of aggregation method is crucial. We evaluate the effectiveness of using the aggregation operator on Amherst41 and Penn94 datasets. Tab. 7 shows the effects of replacing the $\max$ aggregator with $\mathrm{mean}$ and $\mathrm{sum}$ in computing similarity loss. On both datasets, $\max$ outperforms both $\mathrm{sum}$ and $\mathrm{mean}$, indicating the effectiveness of using $\mathrm{max}$ as the aggregation operation.
Intuitively, taking the average of all similarity scores ($\mathrm{mean}$) is sub-optimal. $\mathrm{mean}$ tends to make the node embeddings closer to the average of all clusters belonging to a same class, undermining the purpose of using multiple clusters. Similar to $\mathrm{mean}$, summing up all similarity scores ($\mathrm{sum}$) is more powerful yet requires more data to learn. $\mathrm{max}$ selects the maximum similarity score to compute similarity loss and guides the node embeddings closer to one of the clusters, thus preserving the power of diversity in representation.

#### D.2.2 Additional ablation on the orthogonality loss and similarity loss

**Table 8: Effects of Similarity and Orthogonality losses for datasets in Tab. 1.**
|  | Penn94 | Genius | Cornell5 | Amherst41 | US-election | Wisconsin | Cora | Citeseer | Pubmed |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| DC-GNN | 86.69 (0.22) | 91.70 (0.08) | 84.68 (0.24) | 82.94 (1.59) | 89.59 (1.60) | 91.67 (1.95) | 89.13 (1.18) | 77.93 (1.82) | 91.00 (1.28) |
| (-)$\mathcal{L}_{sim}$, $\mathcal{L}_{ortho}$ | 86.40 (0.25) | 91.68 (0.08) | 84.34 (0.22) | 81.75 (1.07) | 88.61 (1.57) | 89.06 (2.55) | 88.87 (1.23) | 77.27 (1.86) | 91.00 (1.61) |

**Table 9: Effects of Similarity and Orthogonality losses on recent heterophilous datasets [].**
|  | Roman-empire | Amazon-ratings | Minesweeper | Tolokers | Questions |
| --- | --- | --- | --- | --- | --- |
| DC-GNN | 89.96 (0.35) | 51.11 (0.47) | 98.50 (0.21) | 85.88 (0.81) | 78.96 (0.60) |
| (-)$\mathcal{L}_{\rm sim}$, $\mathcal{L}_{\rm ortho}$ | 89.44 (0.72) | 50.32 (0.46) | 98.04 (0.19) | 84.94 (0.59) | 77.30 (0.98) |

We conduct ablation studies on all fourteen datasets we have used. As shown in Tab. 8 and Tab. 9, $\mathcal{L}_{\rm ortho}$ and $\mathcal{L}_{\rm sim}$ generally help to improve the scores across datasets by facilitating the clustering process.

#### D.2.3 Effects of node fidelity term

**Table 10: Dirichlet Energy (DE) with different $\beta$ values. Higher DE indicates increased node distinctiveness.**
|  | $\beta=0.0$ | $\beta=0.5$ | $\beta=1.0$ |
| --- | --- | --- | --- |
| Wisconsin | 0.741288 | 0.963261 | 0.978465 |
| Citeseer | 0.110083 | 0.151699 | 0.206022 |
| Cora | 0.218658 | 0.294024 | 0.330234 |

The node fidelity term encourages the node embeddings to retain some information from the original node features, which serve as initial residual. This technique can also potentially help to alleviate oversmoothing as shown in Klicpera et al. [], Chen et al. []. To validate this, we conduct additional experiments to measure the normalized Dirichlet Energy (DE) [] for DC-GNN on Wisconsin, Cora and Citeseer, using the implementation from  [].

We set $\beta$ to 0, 0.5 and 1 for each dataset to measure how increased weightage of the node fidelity term influences DE, while keeping $\alpha$ constant at 0.5. As observed in Tab. 10, DE positively correlates with $\beta$ on all datasets.

### D.3 Runtime experiments

We measure the average runtime of a DC-MsgPassing layer on the three largest datasets used in our experiments against Pytorch Geometric [] implementation of GATConv and GCNConv. DC-MsgPassing takes less than 4x times GCN and is faster than GAT on these datasets. The results show that DC-MsgPassing is competitive in terms of runtime, confirming our complexity analysis.

**Table 11: Average runtime and dataset statistics comparison on three large-scale datasets. Runtime results are in seconds.**
|  | Penn94 | Cornell5 | Genius |
| --- | --- | --- | --- |
| # Nodes | 41,554 | 18,660 | 421,961 |
| # Edges | 1,362,229 | 790,777 | 984,979 |
| DC-MsgPassing | 0.00852 | 0.00820 | 0.00708 |
| GATConv | 0.01242 | 0.01685 | 0.01636 |
| GCNConv | 0.00220 | 0.00242 | 0.00207 |
| Multiples of GAT | 0.69x | 0.49x | 0.43x |
| Multiples of GCN | 3.88x | 3.38x | 3.42x |

### D.4 Dataset details

#### D.4.1 Dataset description

Roman-empire, Amazon-ratings, Minesweeper, Tolokers and Questions are five datasets proposed in Platonov et al. [] to better evaluate the performance of GNNs under heterophilous settings. The description of each dataset is as follows.

Roman-empire is based on the Roman Empire article from English Wikipedia. Each node in the graph represents one word in the text, and each edge between two words represents either one word following another word or if the two words are connected in the dependency tree. Node features is its FastText word embeddings. The task is to predict a node’s syntactic role.

Amazon-ratings is based on the Amazon product co-purchasing network metadata. Nodes represent products and edges connect products frequently purchased together. Node features are the mean of FastText embeddings for words in product description. The task is to predict the class of products’ ratings.

Minesweeper is a synthetic dataset inspired by the Minesweeper game. The graph is a regular 100x100 grid where each node is connected its eight neighboring nodes. 20% of the nodes are randomly assgined as mines. The node features are one-hot-encoded numbers of the neighboring mines. The task is to predict if the nodes are mines.

Tolokers is based on data from the Toloka crowdsourcing platform. Nodes represent workers who have participated in the selected projects, while edges connect two workers who work on the same task. Node features are based on worker’s profile information and task performance. The task is to predict which workers have been banned.

Questions is based on question-answering data from website Yandex Q. Nodes represent users and edges connect an answer provider to a question provider. Node features are the mean of FastText word embeddings of user profile description, with an additional binary feature indicating users with no descriptions. The task is to predict if the users remain active on the website.

Penn94, Cornell5 and Amherst41 [] are friendship network datasets extracted from Facebook of students from selected universities from 2005. Each node in the datasets represent a student, while node label represents the reported gender of the student. Node features include major, second major/minor, dorm/house, year, and high school.

Wisconsin [] is a web page dataset collected
from the computer science department of Wisconsin Madison. In this dataset, nodes represent web pages and edges are hyperlinks between them. Feature vectors of nodes are bag-of-words representations. The task is to classify the web pages into one of the five categories including student, project, course, staff and faculty.

Genius [] is a sub-network from website genius.com, a crowd-sourced website of song lyrics annotations. Nodes represent users while edges connect users that follow each other. Node features include expertise scores expertise scores, counts of contributions and roles held by users. Around 20% of the users are marked with a "gone" label, indicating that they are more likely to be spam users. The task is to predict which users are marked.

US-election [] is a geographical dataset extracted from statistics of Unite States election of year 2012. Nodes represent US counties, while edges connect bordering counties. Node features include income, education, population etc. The task is a binary classification to predict election outcome.

Cora, Citeseer and Pubmed [] are citation graphs, where each node represents a scientific paper and two papers are connected when a paper cites the other. Each node is labeled with the research field and the task is to predict which field the paper belongs to. All three datasets are homophilous.

#### D.4.2 Dataset statistics

Tab. 12 covers statistics of datasets in Tab. 1. Tab. 13 covers statistics of datasets in Tab. 2.

**Table 12: Dataset statistics for Tab. 1.**
|  | Penn94 | Cornell5 | Amherst41 | Genius | US-election | Wisconsin | Cora | Citeseer | Pubmed |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Edge Hom. | 0.47 | 0.47 | 0.46 | 0.61 | 0.83 | 0.21 | 0.81 | 0.74 | 0.80 |
| Improved Edge Hom. [] | 0.046 | 0.09 | 0.05 | 0.08 | 0.54 | 0.094 | 0.766 | 0.627 | 0.664 |
| # Nodes | 41,554 | 18,660 | 2,235 | 421,961 | 3,234 | 251 | 2,708 | 3,327 | 19,717 |
| # Edges | 1,362,229 | 790,777 | 90,954 | 984,979 | 11,100 | 466 | 5,278 | 4,676 | 44,327 |
| # Node Features | 4814 | 4735 | 1193 | 12 | 6 | 1,703 | 1,433 | 3,703 | 500 |
| # Classes | 2 | 2 | 2 | 2 | 2 | 5 | 6 | 7 | 3 |

**Table 13: Dataset statistics for Tab. 2.**
|  | Roman-Empire | Amazon-Ratings | Minesweeper | Tolokers | Questions |
| --- | --- | --- | --- | --- | --- |
| Edge Hom. | 0.05 | 0.38 | 0.68 | 0.59 | 0.84 |
| Improved Edge Hom. [] | 0.01 | 0.12 | 0.009 | 0.17 | 0.08 |
| # Nodes | 22,662 | 24,492 | 10,000 | 11,758 | 48,921 |
| # Edges | 32,927 | 93,050 | 39,402 | 519,000 | 153,540 |
| # Node Features | 300 | 300 | 7 | 10 | 301 |
| # Classes | 18 | 5 | 2 | 2 | 2 |

##### Homophily matrix

Homophily refers to the degree of similarity between connected neighboring nodes in terms of their features or labels.
There are many types of homophily measures proposed, including edge homophily [], node homophily [], and improved edge homophily []. Homophily matrix proposed in Lim et al. [] is an important metric, since it can better reflect class-wise homophily.
The homophiliy matrix is defined as:

$$ $H_{c_{1},c_{2}}=\frac{|(u,v)\in E:c_{u}=c_{1},\;c_{v}=c_{2}|}{|(u,v)\in E:c_{u }=c_{1}|},$ (22) $$

for classes $c_{1}$ and $c_{2}$, $H_{{c_{1}}{c_{2}}}$ denotes the proportion of edges between from nodes of class $c_{1}$ to nodes of class $c_{2}$. A homophilous graph has high values on the diagonal entries of $H$.

Fig. 6 are the homophily matrices for three well-known homophilous datasets: Cora, Citeseer and Pubmed []. High homophily is signified by the high numbers in diagonal cells, whereas values of non-diagonal cells are mostly less than 0.1. This is different from the homophily matrices of heterophilous datasets, where values of non-diagonal cells are similar or even higher than diagonal cells.

Figure: (a) Cora
Refer to caption: extracted/5620147/arxiv_paper/Figures/cora.png

We show in Fig. 7 the homophily matrices for some heterophilous datasets for comparison.

Figure: (a) Penn94
Refer to caption: extracted/5620147/arxiv_paper/Figures/penn_94_metric.png

### D.5 More experiment details

##### Baselines

To comprehensively evaluate the effectiveness of our method, we compare it against various strong baselines following Platonov et al. [] and  Li et al. []. This includes (1) graph-agnostic model MLP and ResNet [], with two modified versions ResNet+SGC [] and ResNet+adj [];
(2) general-purpose GNN architectures: GCN[], Graph-SAGE [], MixHop [] and GCNII [];
(3) GNN models that leverage attention-based aggregation: GAT [] and Graph Transformer(GT) []. Following Platonov et al. [], we also include two modified architectures GAT-sep and GT-sep, where ego- and neighbor-embeddings are aggregated separately. Following [], we also include GPS^GAT+Performer that achieves the best performance among GPS variants on most datasets as a baseline [, , ].(4) GNN models designed for heterophilous graphs: H2GCN [], CPGNN [], GPR-GNN [], FSGNN [], FFAGCN [], GBK-GNN [], Jacobi-Conv [], WRGAT [],
GPR-GNN [], GGCN [],
ACM-GCN [], LINKX [] and GloGNN/GloGNN++ [].

##### Implementation details.

For global cluster-nodes, we use trainable lookup embeddings to initialize the embeddings. For local cluster-nodes, we fix the number of local clusters to be 2 for every ego-neighborhoods, and initialize the two cluster-node embeddings by the central node features and the average of neighboring node features respectively.
In practice, for local clustering cost matrices M, we rescale and normalize the distance before running the Sinkhorn-Knopp algorithm for numerical stability. We apply non-linear activation function tanh to the messages.

##### Training Settings

We conduct each experiment of DC-GNN using three distinct data splits and present the corresponding mean and standard deviation of the performance metrics. The experiments are executed on a single GPU. The GPUs are from various types—specifically, the V100, A100, GeForce RTX 2080, or 3090—based on their availability at the time the experiments are conducted. For optimization, we employ the Adam optimizer and undertake a grid search of hyperparameters, the specifics of which are in Tab. 14 and Tab. 15. Should the baseline results be publicly accessible, we directly incorporate them into our report. For the datasets where baseline results are missing (Cornell5, Amherst41, US-election), we reproduced them following the code in Li et al. [].

##### Hyperparameters for DC-GNN

For experiments in Tab. 1 and Tab. 2, we fix some hyperparameters and perform grid search for other hyperparameters. To facilitate reproducibility, we document the details of the hyperparameters and search space in Tab. 14 and Tab. 15 respectively. $T^{\Omega}$ refer to the number of iterations of Sinkhorn–Knopp algorithm when solving $P^{\Omega}$, and $T^{\Gamma}$ is the number of Sinkhorn–Knopp iterations for solving $P^{\Gamma}$. We set $|\Omega|$ to be multiples of the number of classes in a dataset.

**Table 14: Hyper-parameter search space of DC-GNN for datasets in Tab. 1.**
|  | Penn94 | Cornell5 | Amherst41 | Genius | US-election | Wisconsin | Cora | Citeseer | Pubmed |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| lr | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.01 | 0.005 | 0.001, 0.002 | 0.005 |
| $\lambda$ | 2 | 2, 5 | 2 | 2 | 2,5 | 2 | 2 | 2 | 2 |
| $T^{\Omega}$ | 10,5,3 | 10,5,3 | 10,5,3 | 10,5,3 | 10,5,3 | 10,5,3 | 10,5,3 | 10,5,3 | 10,5,3 |
| $T^{\Gamma}$ | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 | 5, 3, 1 |
| $|\Gamma_{i}|$ | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| $|\Omega|$ | 2,4, 8, 16, 30 | 2,4, 8, 16, 30 | 2,4, 8, 16, 30 | 2,4,8 | 2,4, 8, 16, 30 | 5, 10, 20 | 6,12,24,48 | 7, 14, 28 | 6, 12 |
| $\alpha$ | 0. 0.2, 0.5, 0.8, 1 | 0. 0.2, 0.5, 0.8, 1 | 0. 0.2, 0.5, 0.8, 1 | 0. 0.2, 0.5, 0.8, 1 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 |
| $\beta$ | 0, 0.2, 0.5, 0.8 | 0, 0.2, 0.5, 0.8 | 0, 0.2, 0.5, 0.8 | 0, 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 |
| # layers in MLP | 1,2 | 1,2 | 1,2 | 1,2 | 1,2 | 1,2 | 1,2 | 1,2 | 1,2 |
| $L:$ # layers | 2, 5 | 2, 5 | 2, 5 | 2, 4, 8,16 | 2, 5 | 2, 5 | 2,4,8,10 | 2,4,8 | 2,4,8 |
| $\omega_{1}$ | 0.001, 0.01 | 0.001, 0.01 | 0.001, 0.01 | 0.001,0.01 | 0.001, 0.01 | 0.001, 0.01 | 0, 0.001, 0.01, 0.05 | 0.001, 0.01 | 0.001, 0.01 |
| $\omega_{2}$ | 0.005, 0.05 | 0.005, 0.05 | 0.005, 0.05 | 0,0.005,0.05 | 0.005, 0.05 | 0.005, 0.05 | 0.005, 0.05, 0.08, 0.1 | 0.005, 0.05 | 0.005, 0.05 |
| epochs | 30 | 30 | 30 | 3000 | 500 | 200 | 200 | 50 | 500 |
| weight_decay | 5e-4 | 5e-4 | 5e-4 | 5e-4 | 5e-4 | 1e-3 | 5e-4 | 5e-4 | 5e-4 |
| aggregation | mean, sum | mean, sum | mean, sum | mean, sum | mean, sum | mean, sum | mean, sum | mean, sum | mean, sum |
| dropout | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 |
| normalization | None | None | None | LN | None | None | None | None | None |
| hidden_channels | 16, 32, 64, 128 | 16, 32, 64, 128 | 16, 32, 64, 128 | 16, 32 | 16, 32, 64, 128 | 16, 32, 64, 128 | 16, 32, 64, 128 | 16, 32, 64, 128 | 16, 32, 64, 128 |

**Table 15: Hyper-parameter search space of DC-GNN for datasets in Tab. 2.**
|  | Roman-Empire | Amazon-Ratings | Minesweeper | Tolokers | Questions |
| --- | --- | --- | --- | --- | --- |
| lr | 0.005 | 0.005 | 0.005 | 0.005 | 0.001 |
| $\lambda$ | 2 | 2 | 2 | 2 | 2 |
| $T^{\Omega}$ | 5 | 5 | 5 | 5 | 5 |
| $T^{\Gamma}$ | 3 | 3 | 3 | 3 | 3 |
| $|\Gamma_{i}|$ | 2 | 2 | 2 | 2 | 2 |
| $|\Omega|$ | 18 | 5, 10 | 2,4,8 | 8, 16 | 2,4,8 |
| $\alpha$ | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 |
| $\beta$ | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 | 0.2, 0.5, 0.8 |
| # layers in MLP | 1,2 | 1,2 | 1,2 | 1,2 | 1,2 |
| $L:$ # layers | 2, 4, 8,16,20 | 2, 4, 8 | 2, 4, 8 | 2, 4, 8 | 2, 4, 8, 16 |
| $\omega_{1}$ | 0.001, 0.01 | 0.001, 0.01 | 0.001, 0.01 | 0.02 | 0.001, 0.01 |
| $\omega_{2}$ | 0.005, 0.05 | 0.005, 0.05 | 0.005, 0.05 | 0.001 | 0.001 |
| epochs | 3000 | 3000 | 3000 | 3000 | 3000 |
| weight_decay | 5e-4 | 5e-4 | 5e-4 | 5e-4 | 5e-4 |
| aggregation | mean, sum | mean, sum | mean, sum | mean, sum | mean, sum |
| dropout | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 |
| normalization | None, LN | None, LN | None, LN | None, LN | None, LN |
| hidden_channels | 16, 32, 64 | 16, 32, 64 | 16, 32, 64 | 16, 32, 64 | 16, 32, 64 |

##### Hyperparameters for baselines

We used the code provided by GloGNN [] to reproduce the baseline results for dataset Cornell5, Amherst41, and US-election. The grid search space for hyper-parameters are listed below. Note that some hyper-parameters only apply to a subset of baselines. All other baselines results are obtained from Li et al. [] and Platonov et al. [].

- •
MLP: hidden dimension $\in\{16,32,64\}$, number of layers $\in\{2,3\}$. Activation function is ReLU.
- •
GCN: lr $\in\{.01,.001\}$, hidden dimension $\in\{4,8,16,32,64\}$. Activation function is ReLU.
- •
GAT: lr $\in\{.01,.001\}$. hidden channels $\in\{4,8,12,32\}$ and gat heads $\in\{2,4,8\}$. number of layers $\in\{2\}$. We use the ELU as activation.
- •
MixHop: hidden dimension $\in\{8,16,32\}$, number of layers $\in\{2\}$.
- •
GCNII: number of layers $\in\{2,8,16,32,64\}$, strength of initial residual connection $\alpha\in\{0.1,0.2,0.5\}$, hyperparameter for strength of the identity mapping $\theta\in\{0.5,1.0,1.5\}$.
- •
H2GCN: hidden dimension $\in\{16,32\}$, dropout $\in\{0,.5\}$, number of layers $\in\{1,2\}$. Model architecture follows Section 3.2 of Zhu et al. [].
- •
WRGAT: lr $\in\{.01\}$, hidden dimension $\in\{32\}$.
- •
GPR-GNN: lr $\in\{.01,.05,.002\}$, hidden dimension $\in\{16,32\}$.
- •
GGCN: lr $\in\{.01\}$, hidden channels $\in\{16,32,64\}$, number of layers $\in\{1,2,3\}$, weight decay $\in\{1e^{-7},1e^{-2}\}$, decay rate $\in\{0,1.5\}$, dropout rate $\in\{0,.7\}$,
- •
ACM-GCN: lr $\in\{.01\}$, weight decay $\in\{5e^{-5},5e^{-4},5e^{-3}\}$, dropout $\in\{0.1,0.3,0.5,0.7,0.9\}$, hidden channels $\in\{64\}$, number of layers $\in\{2\}$, display step $\in\{1\}$.
- •
LINKX: hidden dimension $\in\{16,32,64\}$, number of layers $\in\{1,2\}$. Rest of the hyper-parameter settings follow Lim et al. [].
- •
GloGNN++: lr $\in\{.001,.005,.01\}$ , weight decay $\in\{0,.01,.1\}$, dropout $\in\{0,.5,.8\}$, hidden channels $\in\{128,256\}$, number of layers $\in\{1,2\}$,
$\alpha\in\{0,1\}$,
$\beta\in\{0.1,1\}$,
$\gamma\in\{0.2,0.5,0.9\}$,
$\delta\in\{0.2,0.5\}$,
number of normalization layers $\in\{1,2\}$,
orders $\in\{1,2,3\}$.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 6
presentation: 4
questions: 1.  How to choose model hyper-parameters in practice?
2. I expect the authors to provide the training time and loss curves for the model.
soundness: 2
strengths: 1. The paper is clearly written and easy to read.
2. The proposed end-to-end differentiable model is convincing and the proposed model is trying to solve the important problems encountered in graph representation learning models such as over-squashing and heterophily.
3. Experimental results show that the model works well and experimental details are provided in the appendix.
summary: This paper presents a differentiable end-to-end clustering-based graph neural network for node classification tasks. The proposed model attempts to address over-squashing and heterogeneity. The model is carefully designed to make it differentiable and the authors provide theoretical guarantees on convergence and time complexity.
weaknesses: 1. The model has multiple hyperparameters, which is very confusing for the potential user of the proposed model to select the optimal hyperparameters.
2. The proposed model seems difficult to train and converge despite Theorem 3.3 can provide some guarantee for convergence. I'm not sure if the model can converge only in a very narrow range of hyperparameters, and the code is not open-sourced.
3. $|\Omega|$ can not be removed from asymptotic time complexity $O(T|\mathcal{V}||\Omega|)$ simply because of $|\Omega| \ll |\mathcal{V}|$. ). I suggest that the authors keep the original complexity with these hyperparameters, and then provide a simplified complexity when these hyperparameters are considered as constants. I suggest that the authors report in their experiments the training time and running time of DCGNN in comparison with classical models such as GCN.
4. For this reason and concerns about convergence speed, I am worried about the actual training-to-convergence time might be longer than expect. I note that run times are provided in Appendix D.3, but I would expect the authors to report training times and loss curves.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: see weaknesses
soundness: 2
strengths: 1. Converting the graph structure (adjacency or graph Laplacian) into a bipartite graph as GNN input is interesting.
summary: This paper proposes a GNN that uses cluster-based messaging passing to address the over-squashing and heterophily problems in graph representation learning.The author claims that their cluster-based method (DC-GNN) can handle both long-range dependencies and heterophilous aggregation, by projecting graphs and node neighbors into a bipartite graph of global clusters and local clusters. DC-GNN is different than the current methods by taking the bipartite graph as input instead of using the adjacency or graph Laplacian in the network.
weaknesses: 1. No model access, no code access. (major)

2. Most baseline results reported in Table1 do not align with the results reported in original papers. Where did you get the baseline results? If you run the baselines, please provide codes/loggers or any proof. If not, please cite the sources that you used. (major)

3. Overall contribution is not much, not so exciting. (minor)

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 1
questions: 1.	The authors mention “cluster patterns” many times, but it is unclear what the definition of cluster patterns is?
2.	How does the construction of the bipartite graph scale for extremely large datasets?
3.	The paper uses soft cluster assignments but it lacks the sensitivity analysis of the entropic regularization \lambda and the number of clusters?
4.	From the ablation studies, the contributions of local and global clustering vary. It would be better to provide further intuition or guidelines on when one should prioritize local or global clusters for new datasets.
5.	As the bipartite graph grows with the addition of cluster-nodes, what are the memory implications?
soundness: 2
strengths: 1.	The paper addresses two critical challenges—over-squashing and heterophilous neighborhood aggregation with a unified framework.
2.	The iterative optimization with soft cluster assignment makes it possible to learn both node and cluster-node embeddings efficiently.
3.	Experimental results are sufficient to demonstrate the effectiveness of the proposed model.
summary: The paper proposes a differentiable cluster GNN to deal with the long-rage information propagation and local heterophilous neighborhood aggregation problem in graph domain. Specifically, the framework considers a clustering inductive bias into the message propagation by using additional cluster-nodes. Besides, authors adopt an iterative process to optimize the cluster assignments and node/cluster-node embeddings. Extensive experimental results show the effectiveness of the proposed framework.
weaknesses: 1.	The complexity analysis is very rough, the preprocessing of constructing bipartite graph may cost time, so it would be better to supplement the complexity of preprocessing data.
2.	The introduction of so-called “cluster nodes” is basically assigning a pseudo label to each node, which is quite common among methods on graph heterophily. 
3.	The novelty of directly integrate clustering into the message-passing mechanism is somehow weak.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: See Weakness.
soundness: 3
strengths: 1. This work integrates clustering into message passing, allowing for both global and local pattern capture, which is effective for diverse graphs.
2.  This work demonstrates effectiveness on multiple datasets.
summary: This paper incorporates a clustering inductive bias, combining nodes with "cluster nodes" in a bipartite structure and leveraging optimal transport for differentiable clustering. The framework iteratively optimizes node and cluster-node embeddings and integrates clustering directly within message passing, effectively capturing both global and local patterns.
weaknesses: 1. The author claims that facilitating long-range interactions across distant nodes. Thus it should be tested on long-range graph datasets in [1].
2. The proposed method does not achieve the best performance in all cases in Table 2. More analysis on why performing poorly should be added.
3. [2] and [3] are related to this work. I think quoting them is needed.
4. The OT-based clustering and iterative optimization may add implementation complexity.

[1] Dwivedi, Vijay Prakash, et al. "Long range graph benchmark." Advances in Neural Information Processing Systems 35 (2022): 22326-22340.
[2] Kosmala, Arthur, et al. "Ewald-based long-range message passing for molecular graphs." International Conference on Machine Learning. PMLR, 2023.
[3] Chen, Dexiong, Till Hendrik Schulz, and Karsten Borgwardt. "Learning Long Range Dependencies on Graphs via Random Walks." arXiv preprint arXiv:2406.03386 (2024).

Meta-review: This paper proposes a differentiable cluster-based GNN (DC-GNN) to address over-squashing and heterophily issues in graph neural networks. It integrates clustering into the message-passing mechanism, optimizing both node and cluster-node embeddings iteratively, and provides theoretical guarantees on convergence. The authors present experimental results demonstrating the model’s effectiveness across multiple datasets.

While the proposed DC-GNN model shows promising results, the reviewers have identified several weaknesses that need to be addressed:

1. The model's requirement for additional nodes and updates regarding clustering and assignment results in inefficiencies in both time and memory usage. This significantly limits its ability to handle large-scale datasets, making the trade-off between efficiency and performance unfavorable.
2. The paper lacks a comprehensive analysis of the convergence speed. The model may require more training rounds than initially expected. Additionally, no comparison of convergence behavior has been made between the proposed model and the baseline GNNs.
3. The experimental results are not convincing, as some baseline results are lower than commonly expected benchmarks. The model's improvement, particularly on heterophilic graphs, appears marginal. Additionally, the necessary artifact of the work is not submitted promptly for checking and reproduction.

Based on these weaknesses, we recommend rejecting this paper. We hope this feedback helps the authors improve their paper.

Decision: Reject