Leveraging Personalized PageRank and Higher-Order Topological Structures for Heterophily Mitigation in Graph Neural Networks

1 Introduction

Graph Representation Learning (GRL) tackles the challenge of capturing complex relational information in non-Euclidean structured data, enabling advanced machine learning applications on graphs [Hamilton et al., 2017; Hamilton, 2020]. Graph Neural Networks (GNNs), a cornerstone of GRL, were originally developed under the homophily assumption, where connected nodes share similar properties. Early GNN models focused on homophilic graphs [Veliˇckovi´c et al., 2017;

Xu et al., 2018], yet many real-world networks exhibit heterophily, with connections spanning dissimilar nodes. The performance of GNNs often depends on the graph's homophily level [Pei et al., 2020], driving recent efforts to address heterophily [Jin et al., 2021; Zhu et al., 2022; Wu et al., 2023b; Qiu et al., 2024; Wo et al., 2024].

Despite these advances, most prior work centers on pairwise node relationships [Majhi et al., 2022], which fail to capture the multi-scale, complex interactions prevalent in realworld scenarios while higher-order structures better reflect real-world complex systems [Benson et al., 2016]. To overcome this limitation, higher-order network methods, such as hypergraphs and simplicial complexes (SCs), have emerged. Hypergraphs generalize graphs by allowing edges to connect multiple nodes but often neglect intra-hyperedge relationships, limiting their robustness to noise compared to SCs [Antelmi et al., 2023]. Grounded in algebraic topology [Spanier and Spanier, 1989], SCs model multidimensional relationships, offering richer node context [Hajij et al., 2021] and proving effective in domains like propagation [Iacopini et al., 2019], sensor networks [Schaub et al., 2021], and social systems [Fu et al., 2024; Benson et al., 2016].

A central challenge in higher-order neural networks lies in constructing node relationships using SCs. Recent approaches leverage the Hodge Laplacian from SC theory [Ebli et al., 2020; Roddenberry et al., 2021; Yang et al., 2022], but its computation can reach O(n3 ) complexity in dense graphs, posing scalability issues. Meanwhile, Personalized PageRank (PPR), a random walk variant with teleportation [Page, 1999], efficiently captures local node information. Recent studies integrate PPR with GNNs to enhance representations while preserving efficiency [Gasteiger et al., 2018; Chien et al., 2020; Jack et al., 2023], though they often overlook multi-level, higher-order interactions.

In this paper, we propose HPGNN, a novel framework that combines higher-order topological structures with Personalized PageRank to advance graph representation learning, particularly for heterophilic graphs. Unlike methods limited to node-to-node interactions, HPGNN incorporates relationships across multiple simplicial complexes. Our approach comprises two key components: 1) Higher-order Personalized PageRank (HiPPR), which efficiently computes longrange interactions while reducing noise, and 2) Higher-order Adaptive Spectral Convolution (HiASC), which extends spectral convolution to capture higher-dimensional relationships via SCs. Our contributions are threefold:

• HPGNN pioneers the integration of higher-order information into PPR approximation, bridging simplicial complex theory with personalized random walks.
• We introduce an adaptive PPR matrix operator that encodes interactions among higher-order graph structures.
• Extensive experiments on real-world datasets demonstrate that HPGNN achieves better performance than five out of seven state-of-the-art models.

2 RelatedWork

In this section, we provide a brief discussion of related work on learning from graphs with heterophily and personalized pageRank.

2.1 GNNs with Heterophily

Recent GNNs for heterophilous graphs break traditional homophily assumptions, APPNP [Gasteiger et al., 2018] decouples feature transformation and propagation, while GPRGNN [Chien et al., 2020] adaptively combines graph filters for task-specific aggregation. They enable long-range information propagation and reduce over-smoothing, preserving node-specific features. BernNet expands by leveraging Bernstein polynomial filters to effectively balance low and high frequency signals to handle label inconsistency in heterophily [Chen et al., 2022]. Most recently, POLYGCL [Chen et al., 2024] improves graph contrastive learning via structure-aware sampling and polynomial propagation, preserving structural roles and handling heterophily. Nonetheless, these methods generally lack a high-order perspective to stably suppress boundary noise caused by conflicting class in- formation across nodes.

2.2 Higher-order Simplicial Neural Networks

Simplicial Neural Networks focus on efficiently modeling higher-order structures and embedding them into neural networks. BScNets use block operators to capture dependencies across multi-order structures. [Chen et al., 2022]. Another approach extends the block matrix concept, defining an inter-order adjacency matrix for transformations between simplices [Zhou et al., 2024]. SCN introduces random walks on simplices, enabling comprehensive simplex learning [Wu et al., 2023a]. HiGCN extends SCN with two-step random walks, capturing distant higher-order relationships [Huang et al., 2024]. Despite these advancements, traditional random walk methods face issues like rapid information diffusion and lack of personalized control, hindering effective node importance differentiation.

3 Preliminary

This subsection provides the foundational concepts and problem setting for our work.

3.1 Background

Our work builds on foundational concepts of graphs, graph heterophily, simplicial complexes, and clique complex lifting.

Graph. Consider an undirected, unweighted graph G = (V,E), where V = {v1, v2, . . . , vN} is a finite set of N = |V | nodes, and E ⊆ V × V is the set of M = |E| edges. Let X ∈ RN×F denote the node feature matrix, where F is the feature dimension, and A ∈ {0, 1}N×N represent the adjacency matrix, with Aij = 1 if (vi, vj) ∈ E and 0 otherwise.

Graph Heterophily. Graph homophily reflects the tendency of connected nodes to share similar attributes, while heterophily indicates connections between dissimilar nodes [McPherson et al., 2001]. We quantify homophily using the node homophily ratio:

where yv is the label of node v, and dv is its degree. The ratio Hnode ∈ [0, 1] indicates the level of homophily, with lower values signifying stronger heterophily. This metric is critical for understanding network properties in heterophilic settings [Pei et al., 2020].

Simplicial Complexes. A simplicial complex K is a finite collection of simplices satisfying: 1) a k-simplex σk = {v0, v1, . . . , vk} is the convex hull of k + 1 vertices; 2) if σk ∈ K, all its subsets (faces) are in K; and 3) K is finite [Battiston et al., 2020]. Simplicial complexes model higher-order interactions, where a 0-simplex is a vertex, a 1-simplex an edge, a 2-simplex a triangle, and so forth, enabling the study of topological relationships.

3.2 Problem Setting

In semi-supervised node representation learning, the goal is to learn an encoder f(·) : RN×F × RN×N → RN×D that maps the node feature matrix X and adjacency matrix A to a representation matrix Z = f(X,A) = {z1, z2, . . . , zN}, where zi ∈ RD is the D-dimensional representation of node vi. These representations are optimized for downstream tasks, such as node classification or clustering, particularly in graphs with varying levels of homophily or heterophily.

4 Methodology

Overview. We propose the Higher-order Personalized PageRank Graph Neural Network (HPGNN) to address heterophily and model complex node interactions in graph representation learning. As shown in Figure 1, HPGNN integrates three components: (1) Clique complex lifting [Schaub et al., 2020] to transform the original pairwise graphs into simplicial complexes, providing more comprehensive views of the graph structure; (2) Higher-order Personalized PageRank (HiPPR) to model long-range dependencies while mitigating noise in heterophilic graphs; and (3) Higher-order Adaptive Spectral Convolution (HiASC) to model higher-dimensional relationships via spectral methods. By combining HiPPR and HiASC, HPGNN produces robust node representations that excel in heterophilic settings and remain competitive in homophilic ones, enabling superior performance in downstream tasks.

4.1 Constructing Higher-order Network

To capture higher-order interactions, we transform a graph G = (V,E) into simplicial complexes(SCs) via clique complex lifting [Schaub et al., 2020]. Let C(G) denote the set of all cliques in G, where a clique Ci ⊆ V is a complete subgraph. Each p-clique in C(G) is mapped to a (p-1)-simplex, where 2-simplices are the 3-cliques in the graph, forming the SC K = ⋃Ci∈C(G) σi. This process encodes multi-node interactions beyond pairwise edges.

4.2 Higher-order Personalized PageRank

PPR is a variant based on random walks, it modified the original PageRank algorithm to incorporate a "teleport"mechanism that allows us to personalize the random walk starting from a specific root node [Page, 1999]. This adaptation ensures the influence of the root node is reflected, preserving the local structure around it even in the long run.

On pairwise graphs, The root node s is defined via the teleport vector is, which is a one-hot indicator vector. The PPR of node s can be computed as: πppr(is) = αis + (1 - α) ˜ Aπppr(is), with the teleport (or restart) probability α ∈ (0, 1]. By solving the equation and substituting the indicator vector is with the unit matrix Is, we obtain the PPR matrix:

where ˜ A is the normalized adjacency matrix, defined as ˜ A = AD-1 , A is the adjacency matrix of the graph and D is the degree matrix. In is the identity matrix of size n × n, and n is the number of nodes in the graph. The matrix Πppr thus represents the personalized influence of node s on all other nodes in the graph.

Higher-order Information Integration in PPR

After incorporating the higher-order information, we design to treat the node as a simplice to obtain higher-order PPR. The HiPPR vector Πppr for a given simplix σ can be derived from the concept of a random walk with restarts. Specifically, for the p-order simplicial complex, the corresponding Higher-order Personalized PageRank (HiPPR) Matrix can be calculated as:

where I p is the identity matrix corresponding to the number of simplices of order p and ˜ Ap is the normalized adjacency matrix incorporating the higher-order simplicial structure, representing the adjacency matrix between 0simplices(nodes) and p-simplices (where p > 0) . This innovative formula extends the basic idea of PPR to account for higher-order relationships within the simplicial complex.

In practical computations, approximate algorithms balance speed and accuracy. Recent key approximate approaches [Yang et al., 2024] include Monte Carlo, Power Iteration, and Forward Push. We chose to improve the PwrPushSOR method [Chen et al., 2023], which combines Power Iteration and Forward Push, achieving better error and runtime on undirected graphs. It computes PPR vectors for each node v with O(n · d(v)avg) time and O(n) space complexity. The innovative algorithm HiPwrPushSOR used to approximate the higher-order PPR vector is presented in Algorithm 1. The input is a single root node s, and the output is the PPR vector we aim to compute. For SCs with order p, we compute p PPR vctor for the root node.

The input includes the SC K, the source node s and parameters α, ϵ, λ and ω. Where α controls the transition probability, which affects the range of information propagation in the PPR algorithm; ϵ determines the node activation criteria for queue insertion; λ regulates the stopping condition to guarantee computational precision; and ω ∈ (0, 2) is the Successive Over-Relaxation(SOR) parameter to control the intensity of the update at each iteration. In this study, we adopt the value recommended for undirected graphs [Chen et al., 2023].

In the new algorithm, we iterate all nodes (0 to n-1), treating each as the source node s to calculate its PPR vector, which reflects its relationships with all other nodes. By com- Algorithm 1 HiPwrPushSOR: Higher-order Personalized PageRank Approximation Input: Simplicial complex K, source node s ∈ K, number of

nodes n, number of simplices m Parameter: Precision ϵ, teleport probability α ∈ (0, 1), error threshold λ > 0, relaxation parameter ω > 0 Output: Estimated PPR vector ˆπs ∈ Rn 1: Initialize: epochNum← 8, scanThreshold ← ⌊n/4⌋ 2: Set ˆπ(s, σ) ← 0, r(s, σ) ← 0 for all σ ∈ K; r(s, s) ← 1 3: Initialize FIFO queue Q← ∅; Q.append(s) 4: Set rmax ← λ/m; rsum = 1 5: while Q ̸= ∅ and |Q| ≤ scanThreshold and rsum > λ do 6: σ ←Q.pop() 7: ˆπ(s, σ) ← ˆπ(s, σ) +ωαr(s, σ) 8: rsum ← rsum -ωαr(s, σ) 9: for each τ ∈ Nout(σ) do 10: r(s, τ ) ← r(s, τ ) + ω(1-α)r(s,σ) dσ 11: if τ /∈ Q then 12: Q.append(τ ) 13: if rsum > λ then 14: /* Switch to sequential scan */ 15: for i = 1 to epochNum do 16: r′max ← λi/epochNum/m 17: while rsum >m · r′max do 18: for each σ ∈ K do 19: if |r(s, σ)| ≥ ϵdσ then 20: ˆπ(s, σ) ← ˆπ(s, σ) + ωαr(s, σ) 21: for each τ ∈ Nout(σ) do 22: r(s, τ ) ← r(s, τ ) + ω(1-α)r(s,σ) dσ 23: return ˆπs = [ˆπ(s, σ)]σ∈K

puting PPR vectors for every node, we obtain the PPR matrix. After computing the PPR matrix for each order p, ΠpHiPPR is obtained. Moreover, excessive propagation of higher-order information in node classification can cause redundancy and increased computational cost. To address this, we introduce a constraint in HiPPR, modeling nodes propagating to higherorder neighbors and back, thus efficiently controlling higherorder information flow.

4.3 Higher-order Adaptive Spectral Convolution

Spectral graph convolution

Served as the foundation for many spectral GNNs [Bianchi et al., 2021], spectral graph convolution is typically unified under the expression Z = g(L) · X, where Z represents the output feature matrix and g(L) is the graph filter applied to the graph Laplacian matrix L. The core of spectral graph convolution lies in the design of the graph filter g, which governs how graph signals are manipulated in the spectral domain.

For a given graph Laplacian matrix L = UΛU ⊤, the eigenvalue matrix Λ = diag(λ1, . . . , λm) captures the frequency components of the graph signal, while U = [u1, u2, . . . , um] represents the eigenvectors corresponding to these frequencies. The graph Fourier transform for a graph signal x is given by ˆx = U ⊤x. The goal of spectral graph convolution is to manipulate the frequency response of the graph signal by designing an appropriate graph filter g(L). The graph convolution in the spectral domain can then be expressed as:

However, computing the spectral decomposition of largescale networks is computationally prohibitive. To address this challenge, a practical approach is to approximate the graph filter using polynomial expansions, which avoids direct spectral decomposition. Specifically, the graph filter g(λ) is often approximated as a truncated polynomial: g(λ) = ∑Kk=0 αkλk , where αk are learnable coefficients and K denotes the highest polynomial order considered. This polynomial approximation allows for efficient spectral convolution, bypassing expensive eigen-decompositions.

Incorporating Higher-order Structures

Expanding upon the standard spectral graph convolution, we introduce higher-order spectral convolutions by integrating SCs. In this context, the graph laplacian is generalized to higher-order laplacians, which model interactions within simplices-structures that generalize edges to higherdimensional relationships like triangles and tetrahedra. For each order of the SCs, we compute a PPR laplacian matrix, which captures the connectivity relationships at different orders. The convolution operation for each Higher-order Laplacian is defined as:

where gp(Lp) represents the graph filter applied to the p-th order Laplacian, p corresponds to the order of the SCs. βp,k are learnable coefficients that capture the influence of neighbors at different hops in the higher-order SCs, and k determines the polynomial degree for each PPR laplacian matrix.

The convolution across higher-order graphs allows us to encode richer information from neighbors at varying distances within the graph. For a signal X ∈ Rn×d with d characteristics per node, we first extracts hidden state features for each node and then uses HiPPR to propagate them:

where Z p represents the output combining contributions from the k-th hop, βk are k learnable coefficients, ˜ΠpHiPPR captures neighborhood information at order p, ˜ΠpHiPPR is the normalised version of ΠpHiPPR, formally defined as ˜ΠpHiPPR = D-1/2ΠpHiPPRD-1/2 , where D is the diagonal degree matrix. X is the input feature matrix, ρ is the simplified linear functions suggested by[Wang and Zhang, 2022], and ∥∥ denotes concatenation.

Higher-order Spectral Convolution Formulation

To implement the process, we refer to the work of [Wang and Zhang, 2022], who simplified the functions ρ into linear forms. Building on this, we integrate our novel design of the PPR operator into Higher-order spectral convolutions, leading to the following formulation for the model output:

where βp,k, Θp and W are learnable parameters, ∥∥ denotes concatenation, and P refers to the highest order of the simplicial complex under consideration.

By incorporating higher-order interactions through our proposed PPR-based filter, this framework extends the capabilities of spectral graph convolution, enabling the model to capture complex dependencies and multiscale information from the graph structure. This adaptive spectral convolution method not only enhances the representation of graph signals but also improves the scalability and efficiency of the model.

5 Experiments

We first give the experimental setup, and then compare the performance of our method with baselines on node classification. After that, we analyse the performance of different components of HPGNN and give the parameter analysis. Finally, we provide a visualization task with t-SNE.

5.1 Experimental Setup

Datasets

We adopt seven real world datasets to evaluate the performance of our proposed HPGNN, including two categories:

1. Homophilic Graphs: Citation graphs Cora and Citeseer [Sen et al., 2008]; Social network Photo [McAuley et al., 2015].
2. Heterophilic Graphs: Cornell, Texas, and Wisconsin webpage graphs [Pei et al., 2020]; Actor collaboration network Actor [Pei et al., 2020].

We tabulated the basic statistics of the datasets in Table 1. The number of triangles, represents the count of 2-simplex in the graph, where three nodes are mutually connected.

Baselines

We compare HPGNN with three categories of baselines:

1. Homophily-aware Models: Classical graph convolutional networks based on homophily assumption: ChebNet [Defferrard et al., 2016], GCN [Xu et al., 2018], GAT [Veliˇckovi´c et al., 2017].
2. Heterophily-aware Models: APPNP [Gasteiger et al., 2018] mitigates the limitations of local aggregation by employing a fixed personalized pagerank-based propagation scheme, GPRGNN [Chien et al., 2020] generalizes by introducing learnable propagation weights, enabling the model to flexibly capture multi-hop dependencies through adaptive PageRank coefficients, BernNet [He et al., 2021] adjusts node influences using a probabilistic model to handle heterophily, HiGCN [Huang et al., 2024] uses flower-petals Laplacian operators to mitigate heterophily challenges, and PolyGCL [Chen et al., 2024]integrates polynomial filters into GCL to effectively capture both local and global structural patterns, leading to improved performance on graphs with varying degrees of homophily.

(3) Higher-order Models: S2V [Billings et al., 2019] and SNN [Ebli et al., 2020] embeds simplicial complexes to capture higher-order relationships, SGAT and SGAT-EF [Lee et al., 2022] adapts aggregation weights using attention mechanisms for heterophilic graphs.

Implementation Details

We implement our proposed HPGNN with Pytorch framework on four NVIDIA RTX 3090 GPUs. Following [He et al., 2021], For all datasets, we randomly split nodes into 60%, 20%, and 20% for training, validation, and testing, all methods share the same 10 random splits. We train the proposed HiGCN model with the learning rate lr ∈ {0.01, 0.05, 0.1, 0.2, 0.3} and the weight decay wd ∈ {0.0, 0.0001, 0.001, 0.005, 0.1}. For all methods, we run 100 times with the same partition and report the average and standard deviation of accuracy values.

5.2 Main Experiment Results

Tables 2 shows the comparison of our proposed HPGNN with SOTAs for node classification task on seven real-world network datasets. It can be concluded that: (1) On heterophilic graphs Cornell, Actor, Texas and Wisconsin, HPGNN shows substantial improvement, particularly in datasets with more complex relationships. This suggests its effectiveness in capturing multiscale node dependencies, which is crucial for heterophilic graph structures. (2) HPGNN maintaining competitive performance on homophilic graphs Cora, Citeseer and Photo, showcasing its versatility and robustness across different types of graph structures, a point further supported by the average standard deviation in the table. (3) Compared to higher-order simplicial complex neural networks, our proposed algorithm shows an overwhelming advantage, mainly because it combines multi-scale features and adaptive weights to avoid over-smoothing. (4) Though HPGNN performs competitively on Actor and Texas, exceeding the third-best method by at least 1.49 and 1.65 percentage points, respectively. (5) Notably, the top two results across dataset, aside from ours, come from four different algorithms, demonstrating the stability and consistency of our method, highlighting the impact of higher order techniques in achieving this stability.

5.3 Ablation Study

As shown in Figure 2, ablation studies were performed on four datasets: two homophilic graphs (Cora and Citeseer) and two heterophilic graphs (Actor and Wisconsin). To validate the effectiveness of each component in our proposed model, we conducteded several ablation experiments comparing our model with different variants: 1) HPGNN(w/o High): This variant excludes the higher-order structure, using only pairwise PPR weights and convolution for signal filtering between nodes. For this variant, we set the model hyperparameter order to 1 in different datasets. It shows that removing higher-order structures degrades performance, especially on Wisconsin and Actor, highlighting their role in complex heterophilic graph interactions. 2) HPGNN(w/o High and PPR): This variant removes the PPR mechanism from the pairwise graph, replacing the PPR matrix Π with ˜D-1/2 ˜ A ˜D-1/2 as suggested in [Wu et al., 2019]. It can be analyzed that excluding PPR causes a notable drop in performance, particularly on Wisconsin and Actor, underlining its importance in noise filtering and long-range dependency.

In summary, HPGNN outperforms both variants, demonstrating the effectiveness of the proposed modules. In heterophilic graphs, performance degradation is more significant when higher-order and PPR modules are removed, emphasizing their importance in such environments. Combining both higher-order structures and PPR yields the best results, demonstrating their complementary roles in improving performance.

5.4 Parameter Analysis

We estimate the sensitivity of three important hyperparameters, that is, the order p, the teleportation parameter α and the convolution layer k.

Analysis of Order p. The order p represents the order in higher-order networks SCs. Table 3) shows the results of order p on node classification across Cora, Citeseer, Actor and Cornell four datasets. As observed in the table and also supported by findings in the previous literature [Huang et al., 2024], second order generally yields better performance. Hence, we consistently adopt second-order information in all our experiments to ensure fair and consistent evaluation.

Analysis of Teleportation Parameter α. The teleportation parameter α represents the restart probability in HiPPR. All lines show relatively stable accuracy across the tested range, with minor fluctuations. The results on Wisconsin and Cora exhibit slight peaks near α = 0.15, while Citeseer remains consistently flat. In reference to common practice in most papers, we selected α = 0.15 as the optimal value.

Analysis of Layer k. The layer k represents the layer in higher-order adaptive spectral graph convolutions. We investigated the impact of k in HPGNN, evaluating its influence on model performance across Photo and Actor datasets. As shows in Table 3(b), both lines show a general upward trend as k increases, with minor fluctuations. we selected k = 10 as the optimal value, as it yields the highest accuracy while maintaining stability across the tested range.

5.5 Computationally Efficient

Table 4 illustrates that our model significantly outperforms Hodge-based SNN in terms of speed while maintaining performance comparable to pairwise methods. The runtime and memory usage curves for the Actor dataset, shown in Figure 4 highlight distinct trends. Specifically, it demonstrates that as the graph size increases, both the training time per epoch and peak memory usage exhibit exponential growth. In contrast, with an increase in simplex density, the growth in training time per epoch and peak memory usage follows a more linear trend.

5.6 Visualization

To evaluate the clustering behavior and effectiveness of our proposed model, HPGNN, we present t-SNE visualizations for the Cora and Wisconsin datasets in Figure 5. These plots show how well the models separate classes in a lowdimensional space, with edges omitted for clarity.

As shows Figure 5(c), for the Cora dataset, HPGNN clusters nodes more clearly than GPRGNN, indicating that its integration of high-order information and PPR mechanisms improves the separation of classes in homophilic graphs. Figure 5(e) shows that in the Wisconsin dataset, HPGNN also outperforms other models in clustering heterophilic graphs, with better class separation, showing its robustness in handling complex structures.

6 Conclusion

In this paper, we propose HPGNN, a novel graph neural network model that combines higher-order node relationships with the personalized PageRank (PPR) mechanism to improve node classification, especially in complex graph structures. Unlike previous methods focusing on pairwise relationships, our approach captures multilevel interactions, crucial for heterophilic environments.

Experiments on seven real-world datasets show that HPGNN outperforms state-of-the-art methods in 5 out of 7 cases, especially on heterophilic graphs. Its ability to handle noisy and complex relationships sets it apart, offering a more robust solution for graph representation learning. Further optimization is needed for simpler or sparse higher-order graphs to improve performance across all scenarios.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No.62472305 and No.62272338), the Key R&D Project in Hainan Province, China (No.ZDYF2024SHFZ051), and the Hainan Tropical Ocean Institute Yazhou Bay Innovation Research Institute Major Science and Technology Plan Project (No.2023CXYZD001).

Contribution Statement

Yumeng Wang and Zengyi Wo contributed equally to this work. Minglai Shao is the corresponding author.