# Higher-Order GNNs Meet Efficiency: Sparse Sobolev Graph Neural Networks

Jhony Giraldo, Aref Einizade, Andjela Todorovic, Jhon Castro-Correa, Mohsen Badiey, Thierry Bouwmans, Fragkiskos D. Malliaros

► To cite this version:

Jhony Giraldo, Aref Einizade, Andjela Todorovic, Jhon Castro-Correa, Mohsen Badiey, et al. Higher-Order GNNs Meet Efficiency: Sparse Sobolev Graph Neural Networks. IEEE Transactions on Signal and Information Processing over Networks, 2024, 11, pp.11-22. 10.1109/TSIPN.2024.3503416. hal-04801077

HAL Id: hal-04801077

https://hal.science/hal-04801077v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Higher-Order GNNs Meet Efficiency: Sparse Sobolev Graph Neural Networks

Jhony H. Giraldo, Aref Einizade, Andjela Todorovic, Jhon A. Castro-Correa, Mohsen Badiey, Thierry Bouwmans, and Fragkiskos D. Malliaros

**Abstract**—Graph Neural Networks (GNNs) have shown great promise in modeling relationships between nodes in a graph, but capturing higher-order relationships remains a challenge for large-scale networks. Previous studies have primarily attempted to utilize the information from higher-order neighbors in the graph, involving the incorporation of powers of the shift operator, such as the graph Laplacian or adjacency matrix. This approach comes with a trade-off in terms of increased computational and memory demands. Relying on graph spectral theory, we make a fundamental observation: the regular and the Hadamard power of the Laplacian matrix behave similarly in the spectrum. This observation has significant implications for capturing higher-order information in GNNs for various tasks such as node classification and semi-supervised learning. Consequently, we propose a novel graph convolutional operator based on the sparse Sobolev norm of graph signals. Our approach, known as Sparse Sobolev GNN (S2-GNN), employs Hadamard products between matrices to maintain the sparsity level in graph representations. S2-GNN utilizes a cascade of filters with increasing Hadamard powers to generate a diverse set of functions. We theoretically analyze the stability of S2-GNN to show the robustness of the model against possible graph perturbations. We also conduct a comprehensive evaluation of S2-GNN across various graph mining, semi-supervised node classification, and computer vision tasks. In particular use cases, our algorithm demonstrates competitive performance compared to state-of-the-art GNNs in terms of performance and running time.

**Index Terms**—Graph neural networks, sparse graph convolutions, higher-order convolutions, graph spectrum, Sobolev norm

## I. INTRODUCTION

GRAPH representation learning and its applications have garnered significant attention in recent years [1]–[4]. Notably, Graph Neural Networks (GNNs) have emerged as a powerful extension of Convolutional Neural Networks (CNNs) for modeling non-Euclidean data as graphs. GNNs have

Jhony H. Giraldo, Aref Einizade, and Andjela Todorovic are with the LTCI, Télécom Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France e-mail: jhony.giraldo@telecom-paris.fr, aref.einizade@telecom-paris.fr, andela.todorovic@ip-paris.fr.

Jhon A. Castro-Correa and Mohsen Badiey are with the Department of Electrical and Computer Engineering, University of Delaware, Newark, DE, USA. E-mail: jcastro@udel.edu, badiey@udel.edu.

Thierry Bouwmans is with the laboratoire MIA, Mathématiques, Image et Applications, La Rochelle Université, 17000 La Rochelle, France e-mail: tbouwman@univ-lr.fr.

Fragkiskos D. Malliaros is with Université Paris-Saclay, CentraleSupélec, Inria, Centre for Visual Computing (CVN), Gif-Sur-Yvette, France. E-mail: fragkiskos.malliaros@centralesupelec.fr.

This work was supported in part by the ANR (French National Research Agency) under the JCJC project GraphIA (ANR-20-CE23-0009-01), by the Office of Naval Research ONR (Grant No. N00014-21-1-2760), and by the center Hi! PARIS.

been successfully applied in various domains, including semi-supervised learning [3], clustering [5], point cloud semantic segmentation [6], misinformation detection [7], and molecular modeling [8], among others.

Most GNNs update their node embeddings by performing specific operations within the neighborhood of each node through message passing [9]. However, this updating rule has limitations when it comes to capturing higher-order<fnref n="1" /> relationships between nodes since it only leverages information from the immediate neighbors (1-hop) of each vertex. To address this limitation, previous methods in GNNs have attempted to capture higher-order connections by incorporating powers of the sparse shift operator of the graph, such as the adjacency or Laplacian matrix [2], [10]. Nevertheless, these methods encounter computational and memory issues, attributed either to the densification of the shift operator [10] or a bottleneck during the implementation of recursive graph diffusion operations [2].

This work explores a somewhat radical departure from prior approaches in higher-order graph convolutions for GNNs. We observe that the eigenvalues of the higher-order Laplacian matrix in a weighted graph exhibit similar behavior when subjected to regular and Hadamard power operations (the Hadamard power is also known as element-wise power). To establish the precise relationship between the spectra of the regular and sparse higher-order Laplacian matrix, we rely on tools from spectral graph theory [11] and the Schur product theorem [12]. This fundamental observation enables us to design accurate and efficient sparse graph convolutions for GNNs. More specifically, we propose a novel approach called Sparse Sobolev GNN (S2-GNN), which employs a cascade of sparse higher-order filtering operations. This allows for the computation of a more diverse set of functions at each layer. To achieve this goal, we introduce a new sparse Sobolev norm, drawing inspiration from the literature of Graph Signal Processing (GSP) [13], [14].

In S2-GNN, the information of the higher-order filtering operations is aggregated using either a linear combination or a Multi-Layer Perceptron (MLP) fusion layer. We conduct thorough evaluations of S2-GNN in various semi-supervised learning tasks across several domains, including tissue phenotyping in colon cancer histology images [15], text classification of news [16], activity recognition with sensors [17], and recognition of spoken letters [18]. Furthermore, we test

<footnote n="1">¹We use the term “higher-order” to refer to higher-order hop information (more than one hop).</footnote>
our algorithm in node classification tasks using benchmark networks commonly employed in the literature [3], [19]. S2-GNN demonstrates superior or competitive performance compared to a broad range of GNN methods that employ first-order or higher-order graph convolutions [2], [3], [10], various types of multi-head attentions [4], [20], and graph transformer architectures [21].

In this work, we build upon and improve our preliminary study [22]. Specifically, we provide a more comprehensive theoretical explanation, along with an extensive experimental evaluation and in-depth analysis of novel findings. Additionally, we add a rigorous theoretical study about the stability properties of S2-GNN against possible graph perturbations, which is an important and practical issue in real-world scenarios [23], [24]. Furthermore, we empirically observe that utilizing smooth-learned graphs [25], [26] improves the performance of both S2-GNN and baseline GNN methods compared to the commonly used *k*-Nearest Neighbors (*k*-NN) technique for graph construction. Smoother graphs exhibit better results as they promote homophily within the graph.

The main contributions of this work can be summarized as follows:

1) We show that the spectrum of higher-order weighted graphs, using the Laplacian matrix, exhibits similar behavior to their sparse counterparts. This finding enables the design of GNN architectures that deliver superior or competitive performance while maintaining efficiency.

2) We propose S2-GNN, a novel GNN architecture that utilizes a cascade of higher-order sparse filters. S2-GNN effectively leverages sparse higher-order operations without excessively increasing complexity.

3) We rigorously analyze the stability properties of S2-GNN, shedding light on its robustness against possible graph perturbations.

4) We conduct extensive experimental evaluations on multiple publicly available benchmark datasets and compare S2-GNN against nine state-of-the-art GNNs. Our algorithm outperforms or achieves competitive results compared to prior methods.

The remainder of the paper is organized as follows: Section II presents an overview of the related work, while Section III provides an introduction to the preliminary concepts and foundations relevant to this study. Section IV provides a detailed explanation of the proposed S2-GNN model and its theoretical stability analysis. In Section V, we present the experimental framework, results, and discussion. Finally, Section VI offers the concluding remarks of the paper.

## II. RELATED WORK

The study of graphs is a well-established field in machine learning and mathematics [11], [27]. With the emergence of GNNs in 2014, motivated by the success of CNNs in regular-structured data like images, a wide range of GNN models have been proposed to learn representations of graph-structured data. Bruna *et al.* proposed the first modern GNN by extending the convolutional operator of CNNs to graphs [1]. Defferrard *et al.* incorporated concepts from GSP to propose localized

spectral filtering [2], while Kipf and Welling approximated the spectral filtering operation to enable efficient Graph Convolution Networks (GCNs) [3]. These works demonstrated the potential of GNNs and have served as inspiration for the development of various GNN models.

Veličković *et al.* [4] proposed Graph Attention Networks (GATs), which employ an attention mechanism to learn different weights for different nodes in the graph. This enables GATs to effectively model complex relationships between nodes, albeit with increased computational complexity. Subsequent works have further explored attention mechanisms in GNNs, such as [20], [28].

The Simple Graph Convolution (SGC) model was introduced by Wu *et al.* [29], which simplifies the GCN by removing the non-linear activation functions, resulting in improved efficiency but reduced expressivity. Later, Chiang *et al.* [30] proposed ClusterGCN, a scalable GNN model that leverages graph clustering to enhance efficiency on large-scale graphs. By partitioning the graph into clusters and applying GNNs to each cluster in parallel, ClusterGCN can learn representations of large-scale graphs. However, ClusterGCN sacrifices global graph information and the exchange of messages between clusters. Frasca *et al.* proposed Scalable Inception Graph Neural Networks (SIGN) [10], which leverage higher-order relationships in graphs to improve expressivity at the expense of scalability. SIGN computes the powers of the adjacency matrix, making it unsuitable for large-scale problems due to its complexity.

More recently, graph transformers have emerged as an alternative to GNNs for learning graph-structured data. Some graph transformers in the literature, similar to GATs, employ attention mechanisms computed in the neighborhood of each node, as shown in [21]. Another approach in graph transformers is to compute attention mechanisms (with or without positional embeddings) by leveraging fully connected graphs [31]. However, the latter approach becomes computationally prohibitive for large-scale applications, as the complexity grows quadratically with the number of nodes in the graph. For a comprehensive review of GNNs, we recommend referring to the survey by Zhang *et al.* [9].

Previous studies have focused on improving scalability at the expense of expressivity [29], [30] or vice versa [2], [10]. Although GAT, its derivatives, and Transformers have been explored, they are not optimal solutions to this problem due to the high computational burden of the attention mechanism and the multi-head function. As a result, outside the GNN community, the GCN model [3] is the preferred choice due to its simplicity and high efficiency. In this paper, we propose an accurate and efficient GNN architecture based on higher-order sparse convolutions. Our algorithm improves accuracy while maintaining low computational and memory footprints.

## III. PRELIMINARIES

### A. Mathematical Notation

In this work, calligraphic letters such as $V$ denote sets, and $|V|$ represents the cardinality of the set. Uppercase boldface letters such as $A$ represent matrices, while lowercase boldface letters such as $x$ denote vectors. $I$ denotes the identity
matrix, and **1** represents a vector of ones with appropriate dimensions. The superscripts (·)ᵀ and (·)ᵀ correspond to transposition and Hermitian transpose, respectively. diag(**x**) denotes a diagonal matrix with entries given by the vector **x** = [x₁, x₂, ..., xₙ]ᵀ ∈ ℝⁿ. tr(·) is the trace of a matrix. The ℓ₂-norm of a vector is denoted by ||·||₂, the Frobenius norm is denoted by ||·||ₚ, and the entry-wise norm-1 of a matrix is denoted by ||·||₁,₁. The symbols ◦ and ⊗ represent the Hadamard and Kronecker products between matrices, respectively.

## B. Graph Signals

A graph is a mathematical entity represented as $G = (\mathcal{V}, \mathcal{E})$, where $\mathcal{V} = \{1, \dots, N\}$ is the set of $N$ nodes and $\mathcal{E} \subseteq \{(i, j) | i, j \in \mathcal{V} \text{ and } i \neq j\}$ is the set of edges between nodes $i$ and $j$. The weighted adjacency matrix of the graph is denoted as $\mathbf{A} \in \mathbb{R}^{N \times N}$, where $\mathbf{A}(i, j) = a_{i,j} \in \mathbb{R}_+$ represents the weight of the edge $(i, j)$, and $\mathbf{A}(i, j) = 0 \quad \forall (i, j) \notin \mathcal{E}$. For undirected graphs, $\mathbf{A}$ is symmetric. A graph signal is a function $x: \mathcal{V} \to \mathbb{R}$ and is represented as $x \in \mathbb{R}^N$, where $x(i)$ is the graph signal evaluated on the $i$th node. The degree matrix of $G$ is a diagonal matrix given by $\mathbf{D} = \text{diag}(\mathbf{A}1)$. Different definitions exist for the Laplacian matrix, including the combinatorial Laplacian $\mathbf{L} = \mathbf{D} - \mathbf{A}$, and the symmetric normalized Laplacian $\mathbf{\Delta} = \mathbf{I} - \mathbf{D}^{-\frac{1}{2}}\mathbf{A}\mathbf{D}^{-\frac{1}{2}}$ [32]. The Laplacian matrix is a positive semi-definite matrix for undirected graphs, with eigenvalues<fnref n="2" /> $0 = \lambda_1 \le \lambda_2 \le \dots \le \lambda_N$ and corresponding eigenvectors $\{\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_N\}$. The adjacency and Laplacian matrices, as well as their normalized forms, are examples of the shift operator $\mathbf{S} \in \mathbb{R}^{N \times N}$ in GSP.

The graph Fourier basis of $\mathbf{L}$ is given by the spectral decomposition $\mathbf{L} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^\mathrm{T}$ [32], where $\mathbf{U} = [\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_N]$ and $\mathbf{\Lambda} = \text{diag}([\lambda_1, \lambda_2, \dots, \lambda_N]^\mathrm{T})$. The Graph Fourier Transform (GFT) $\hat{\mathbf{x}}$ of the graph signal $\mathbf{x}$ is defined as $\hat{\mathbf{x}} = \mathbf{U}^\mathrm{T}\mathbf{x}$, while the inverse GFT is given by $\mathbf{x} = \mathbf{U}\hat{\mathbf{x}}$ [32]. Although we use the spectral definitions of graphs to establish the relationship between the spectrum of the higher-order Laplacian and its sparse counterpart in this work, the eigenvalue decomposition is not required to implement S2-GNN.

## C. Graph Convolutional Filters

The concept of graph convolutional filters was introduced in the literature of GSP [33].

**Definition 1.** Given a set of parameters $\mathbf{q} = [q_0, q_1, \dots, q_K]^T$ and a graph signal $\mathbf{x} \in \mathbb{R}^N$, a graph convolutional filter of order $K$ is a linear mapping $Q: \mathbb{R}^N \to \mathbb{R}^N$ comprising a linear combination of $K$ shifted signals:

$$Q(\mathbf{x}) = \sum_{k=0}^{K} q_k \mathbf{S}^k \mathbf{x} = \mathbf{Q}(\mathbf{S})\mathbf{x}, \quad (1)$$

where $\mathbf{Q}(\mathbf{S}) = \sum_{k=0}^{K} q_k \mathbf{S}^k$ is the $N \times N$ polynomial filtering matrix.

The output signal $Q(\mathbf{x})$ at node $i$ is a linear combination of signal values located within up to $K$-hops away from node $i$.

given by $q_0\mathbf{x}(i)+q_1[\mathbf{S}\mathbf{x}](i)+\dots+q_K[\mathbf{S}^K\mathbf{x}](i)$, where $[\mathbf{S}^K\mathbf{x}](i)$ is the ith value of the vector $\mathbf{S}^K\mathbf{x}$.

## IV. SPARSE SOBOLEV GRAPH NEURAL NETWORK

Figure 1 provides an overview of our S2-GNN algorithm. S2-GNN utilizes either a k-Nearest Neighbors (k-NN) method or the procedure described in [25] to establish a graph representation for a specific problem when the structure is unavailable. In cases where the structure is already known, step a) in Fig. 1 is not required. Consequently, we employ the proposed GNN architecture to address the primary problem at hand<fnref n="3" />. Within S2-GNN, a modified sparse Sobolev term serves as the shift operator in the propagation function of the GNN.

### A. Sobolev Norm

The Sobolev norm in GSP has been used as a regularization term to solve problems in 1) video processing [34], 2) modeling of infectious diseases [35], and 3) interpolation of graph signals [13], [14].

**Definition 2** (Pesenson [13]). For fixed parameters $\epsilon \ge 0$ and $\rho \in \mathbb{R}$, the Sobolev norm is defined as $||\mathbf{x}||_{\rho, \epsilon} \triangleq ||(\mathbf{L} + \epsilon\mathbf{I})^{\rho/2}\mathbf{x}||$. When $\mathbf{L}$ is symmetric, we can rewrite $||\mathbf{x}||_{\rho, \epsilon}^2$ as follows:

$$||\mathbf{x}||_{\rho, \epsilon}^2 = \mathbf{x}^\mathrm{T}(\mathbf{L} + \epsilon\mathbf{I})^\rho\mathbf{x}. \quad (2)$$

We divide the analysis of (2) into two steps: 1) when $\epsilon = 0$, and 2) when $\rho = 1$. For $\epsilon = 0$ in (2) and considering that $\mathbf{U}$ is orthonormal, we obtain:

$$||\mathbf{x}||_{\rho, 0}^2 = \mathbf{x}^\mathrm{T}\mathbf{L}^\rho\mathbf{x} = \mathbf{x}^\mathrm{T}\mathbf{U}\mathbf{\Lambda}^\rho\mathbf{U}^\mathrm{T}\mathbf{x} = \hat{\mathbf{x}}^\mathrm{T}\mathbf{\Lambda}^\rho\hat{\mathbf{x}} = \sum_{i=1}^{N} \hat{\mathbf{x}}^2(i)\lambda_i^\rho. \quad (3)$$

Notice that the spectral components $\hat{\mathbf{x}}(i)$ are penalized with powers of the eigenvalues $\lambda_i^\rho$ of $\mathbf{L}$. Since the eigenvalues are ordered in increasing order, the higher frequencies of $\hat{\mathbf{x}}$ are penalized more than the lower frequencies when $\rho = 1$, leading to a smooth function in $G$. For $\rho > 1$, the GFT $\hat{\mathbf{x}}$ is penalized with a more diverse set of eigenvalues. Similarly, we can analyze the adjacency matrix $\mathbf{A}$ using the eigenvalue decomposition $\mathbf{A}^\rho = (\mathbf{V}\mathbf{\Sigma}\mathbf{V}^\mathrm{T})^\rho = \mathbf{V}\mathbf{\Sigma}^\rho\mathbf{V}^\mathrm{T}$, where $\mathbf{V}$ is the matrix of eigenvectors, and $\mathbf{\Sigma}$ is the matrix of eigenvalues of $\mathbf{A}$. For $\mathbf{A}$, the GFT can be defined as $\hat{\mathbf{x}} = \mathbf{V}^\mathrm{T}\mathbf{x}$.

For the second analysis of (2), when $\rho = 1$, we have:

$$||\mathbf{x}||_{1, \epsilon}^2 = \mathbf{x}^\mathrm{T}(\mathbf{L} + \epsilon\mathbf{I})\mathbf{x}. \quad (4)$$

The term $(\mathbf{L} + \epsilon\mathbf{I})$ in (4) is associated with a better condition number<fnref n="4" /> than using $\mathbf{L}$ alone. Better condition numbers are associated with faster convergence rates in gradient descent methods, as demonstrated in [14]. For the Laplacian matrix $\mathbf{L}$, we know that $\kappa(\mathbf{L}) = \frac{|\lambda_{\max}(\mathbf{L})|}{|\lambda_{\min}(\mathbf{L})|} \approx \frac{\lambda_{\max}(\mathbf{L})}{0} \to \infty$, where $\kappa(\mathbf{L})$ is the condition number of $\mathbf{L}$, $\lambda_{\max}(\mathbf{L})$ and $\lambda_{\min}(\mathbf{L})$ are the maximum and minimum eigenvalues of $\mathbf{L}$, respectively. Since $\kappa(\mathbf{L}) \to \infty$, relying solely on the Laplacian matrix can result in a poorly conditioned problem, particularly when solving an

<footnote n="2">²$\lambda_N \le 2$ in the case of the symmetric normalized Laplacian $\mathbf{\Delta}$.</footnote>

<footnote n="3">³In this paper, we focus on the node and some graph classification tasks. However, our framework can be easily extended to other problem domains, such as link prediction.</footnote>

<footnote n="4">⁴The condition number $\kappa(\mathbf{L})$ associated with the square matrix $\mathbf{L}$ is a measure of how well- or ill-conditioned is the inversion of $\mathbf{L}$.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: The pipeline of our S2-GNN algorithm. S2-GNN can be used in a broad range of data such as images, text, and videos, among others. However, the step of mapping the original dataset to the data matrix $X \in \mathbb{R}^{N \times M}$ could be different in each case. Our framework is composed of (a) inference of the graph topology and (b) the S2-GNN architecture.</figcaption></figure>

optimization problem [14]. On the other hand, for the Sobolev term, we observe that:

$$L + \epsilon I = U\Lambda U^T + \epsilon I = U(\Lambda + \epsilon I)U^T. \quad (5)$$

Therefore, $\lambda_{\min}(L + \epsilon I) = \epsilon$, i.e., $L + \epsilon I$ is positive definite ($L + \epsilon I > 0$) for $\epsilon > 0$, and:

$$\kappa(L + \epsilon I) = \frac{|\lambda_{\max}(L + \epsilon I)|}{|\lambda_{\min}(L + \epsilon I)|} = \frac{\lambda_{\max}(L) + \epsilon}{\epsilon} < \kappa(L); \quad \forall \epsilon > 0. \quad (6)$$

Namely, $L + \epsilon I$ has a better condition number than $L$. The significance of a better condition number in the context of GNNs might not be immediately apparent since the inverses of the Laplacian or adjacency matrices are not required for performing the propagation rules. However, several studies have highlighted the adverse effects of poorly conditioned matrices. For instance, Kipf and Welling [3] employed a renormalization trick ($A + I$) in their filtering operation to avoid exploding/vanishing gradients. Similarly, Wu et al. [29] demonstrated that adding the identity matrix to $A$ shrinks the graph spectral domain, resulting in a low-pass-type filter.

The previous theoretical analysis shows the benefits of the Sobolev norm in two aspects: 1) the computation of diverse frequencies in (3), and 2) the better condition number in (6).

## B. Sparse Sobolev Norm

The use of $L$, $A$, or their normalized forms in GNNs is computationally efficient due to their sparsity. This allows us to perform a small number of sparse matrix operations, which is advantageous for computational efficiency. However, when considering the Sobolev norm, the term $(L + \epsilon I)^\rho$ can rapidly become a dense matrix for large values of $\rho$, leading to computational and memory issues as illustrated in Fig. 2. To address this problem, we employ a sparse Sobolev norm that preserves the sparsity level, ensuring computational feasibility and efficient memory usage.

**Definition 3.** Let $L \in \mathbb{R}^{N \times N}$ be the Laplacian matrix of a graph $G$. For fixed parameters $\epsilon \ge 0$ and $\rho \in \mathbb{N}$, the sparse Sobolev term for GNNs is introduced as the $\rho$ Hadamard

<figure><img src="image_3.png" /><figcaption>Fig. 2: Process of using higher-order relationships in GNNs vs. higher-order sparse convolutions.</figcaption></figure>

multiplications of $(L + \epsilon I)$ (also known as the Hadamard powers) such that:

$$(L + \epsilon I)^{(\rho)} = \underbrace{(L + \epsilon I) \circ (L + \epsilon I) \circ \dots \circ (L + \epsilon I)}_{\rho \text{ times}}. \quad (7)$$

For example, $(L + \epsilon I)^{(2)} = (L + \epsilon I) \circ (L + \epsilon I)$. Thus, the sparse Sobolev norm is given by:

$$\|x\|_{(\rho), \epsilon} \triangleq \|(L + \epsilon I)^{(\rho/2)} x\|. \quad (8)$$

Let $\langle x, y \rangle_{(\rho), \epsilon} = x^T(L + \epsilon I)^{(\rho)}y$ be the inner product between two graph signals $x$ and $y$ that induces the associated sparse Sobolev norm.

**Theorem 1.** The sparse Sobolev norm $\|x\|_{(\rho), \epsilon} \triangleq \|(L + \epsilon I)^{(\rho/2)}x\|$ satisfies the basic properties of vector norms for $\epsilon > 0$ (for $\epsilon = 0$, we obtain a semi-norm).

*Proof: See Appendix A.*

The sparse Sobolev term in (7) possesses the property of maintaining the same sparsity level as $L + \epsilon I$ for any value of $\rho$. It is worth noting that $(L + \epsilon I)^\rho$ is equivalent to the sparse Sobolev term under two conditions: 1) restricting $\rho$ to be in $\mathbb{N}$, and 2) replacing the conventional matrix multiplication with
<figure><img src="image_4.png" /><figcaption>Fig. 3: Eigenvalues penalization for the non-sparse and sparse matrix multiplications of the combinatorial Laplacian matrix.</figcaption></figure>

the Hadamard product. While the theoretical properties of the Sobolev norm in (3) and (6) do not directly extend to its sparse counterpart, we can gain some theoretical insights by leveraging the concepts of Kronecker products and the Schur product theorem [12].

**Theorem 2.** Let $\mathbf{L}$ be any Laplacian matrix of a graph with eigenvalue decomposition $\mathbf{L} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^T$. We observe that:

$$ \mathbf{L}^{(\rho)} = \mathbf{P}_{(\rho)}^{\mathrm{T}} (\mathbf{U} \otimes \mathbf{U}_{(\rho-1)})(\mathbf{\Lambda} \otimes \mathbf{\Lambda}_{(\rho-1)})(\mathbf{U}^{\mathrm{T}} \otimes \mathbf{U}_{(\rho-1)}^{\mathrm{T}})\mathbf{P}_{(\rho)}, \quad (9) $$

where $\mathbf{U}_{(\rho-1)}$ and $\mathbf{\Lambda}_{(\rho-1)}$ are respectively the matrices of eigenvectors and eigenvalues of $\mathbf{L}^{(\rho-1)}$, $\mathbf{P}_{(\rho)} \in \{0, 1\}^{N^2 \times N}$ is a partial permutation matrix, and:

$$ \mathbf{L}^{(2)} = \mathbf{P}_{N}^{\mathrm{T}}(\mathbf{U} \otimes \mathbf{U})(\mathbf{\Lambda} \otimes \mathbf{\Lambda})(\mathbf{U}^{\mathrm{T}} \otimes \mathbf{U}^{\mathrm{T}})\mathbf{P}_{N}. \quad (10) $$

*Proof: See Appendix B.*

Theorem 2 provides the closed-form solution for the spectrum of the Hadamard power of the Laplacian matrix. It reveals that the spectrum of the Hadamard product is a compressed representation of the Kronecker product of its spectral components. In the case of the sparse Sobolev term utilized in S2-GNN, denoted by $(\mathbf{L} + \epsilon\mathbf{I})^{(\rho)}$, the spectral components of the graph undergo changes for each value of $\rho$, as shown in (9).

By applying the Schur product theorem, which states that the Hadamard product between two positive definite matrices is also positive definite, we can analyze the condition number of the Hadamard powers [12]. We know that $(\mathbf{L} + \epsilon\mathbf{I})^{(\rho)} \succ 0; \forall \epsilon > 0$ since $(\mathbf{L} + \epsilon\mathbf{I}) \succ 0; \forall \epsilon > 0$, and therefore $\kappa((\mathbf{L} + \epsilon\mathbf{I})^{(\rho)}) < \infty$. For the adjacency matrix $\mathbf{A}$, its eigenvalues lie within the interval $[-d, d]$, where $d$ is the maximum degree of the vertices in $G$ (see Theorem 8.5 in [36]). Therefore, we can bound the eigenvalues of $\mathbf{A}$ into the interval $[-1, 1]$ by normalizing the adjacency matrix such that $\mathbf{A}_N = \mathbf{D}^{-\frac{1}{2}}\mathbf{A}\mathbf{D}^{-\frac{1}{2}}$. Thus, we can establish that $\mathbf{A}_N + \epsilon\mathbf{I} \succ 0; \forall \epsilon > 1$, and subsequently, $(\mathbf{A}_N + \epsilon\mathbf{I})^{(\rho)} \succ 0; \forall \epsilon > 1$. We can assert that the theoretical developments presented in Section IV-A for the sparse Sobolev norm hold to some extent, including a more diverse set of frequencies and a better condition number. Figure 3 depicts five normalized eigenvalue penalizations for $\mathbf{L}^{\rho}$ (non-sparse) and $\mathbf{L}^{(\rho)}$ (sparse), revealing that the normalized spectrum of $\mathbf{L}^{\rho}$ and $\mathbf{L}^{(\rho)}$ exhibit a high degree of similarity.

### C. Graph Sparse Sobolev Layer

In the S2-GNN architecture, we introduce the propagation rule using the graph convolution filter in Definition 1. Let $\mathbf{S} \in \mathbb{R}^{N \times N}$

be the shift operator of a graph, and let $\mathbf{H}^{(l)} \in \mathbb{R}^{N \times H_l}$ be the output of the previous layer of the GNN (matrix of node embeddings), such that $\mathbf{H}^{(1)} = \mathbf{X} \in \mathbb{R}^{N \times M}$ is an $M$-dimensional graph signal (matrix of features). The propagation rule of the $\rho$th branch of S2-GNN is defined as:

$$ \mathbf{B}_{\rho}^{(l+1)} = \sigma \left( \mathbf{S}_{\rho} \mathbf{H}^{(l)} \mathbf{W}_{\rho}^{(l)} \right), \quad (11) $$

where $\sigma(\cdot)$ is some activation function, $\mathbf{S}_{\rho} = \mathbf{S}^{(\rho)}$, and $\mathbf{W}_{\rho}^{(l)} \in \mathbb{R}^{H_l \times H_{l+1}}$ is the matrix of trainable weights in the $\rho$th branch of the layer $l$. The matrix $\mathbf{W}_{\rho}^{(l)}$ corresponds to a Multi-Layer Perceptron (MLP) with $H_{l+1}$ hidden units. In practice, we typically include a bias term in the MLP; however, we have omitted it here for simplicity. Notice that (11) is an approximation of the graph convolutional filter in Definition 1 when substituting the regular matrix power by the Hadamard power, $q_i = 1$ for $i = \rho$, $q_i = 0; \forall i \neq \rho$, the input signal is a $H_l$-dimensional graph signal, and adding an MLP at the end of the filtering process.

The S2-GNN architecture consists of multiple layers, where each layer involves computing a cascade of $\alpha + 1$ propagation rules (11) with increasing Hadamard powers. Finally, the outputs from all branches are combined using some fusion mechanism, such as a linear combination layer or an MLP fusion technique.

**Linear combination layer:** Let $\mathbf{B}_{\rho}^{(l+1)} \in \mathbb{R}^{N \times H_{l+1}}$ be the output of the $\rho$th branch in the layer $l$ of S2-GNN. The output of the $l$th layer of S2-GNN using the linear combination layer is given by:

$$ \mathbf{H}^{(l+1)} = \sum_{i=0}^{\alpha} \mu_{l,i} \mathbf{B}_{i}^{(l+1)}, \quad (12) $$

where $\mu_{l,i}$ is the learnable weight associated with $\mathbf{B}_{i}^{(l+1)}$, and we consider $\mathbf{B}_{0}^{(l+1)} = \sigma \left( \mathbf{H}^{(l)} \mathbf{W}_{0}^{(l)} \right) = \sigma \left( \mathbf{H}^{(l)} \mathbf{W}_{0}^{(l)} \right)$.

**MLP fusion:** Let $\mathbf{B}_{\rho}^{(l+1)} \in \mathbb{R}^{N \times H_{l+1}}$ be the output of the $\rho$th branch in the layer $l$ of S2-GNN. The output of the $l$th layer of S2-GNN using the MLP fusion is given by:

$$ \mathbf{H}^{(l+1)} = \left[ \mathbf{B}_{0}^{(l+1)}, \mathbf{B}_{1}^{(l+1)}, \dots, \mathbf{B}_{\alpha}^{(l+1)} \right] \mathbf{W}_{\mathrm{MLP}}^{(l)}, \quad (13) $$

where $[\mathbf{B}_{0}^{(l+1)}, \mathbf{B}_{1}^{(l+1)}, \dots, \mathbf{B}_{\alpha}^{(l+1)}]$ is a concatenation of the matrices $\mathbf{B}_{i}^{(l+1)} \forall 0 \le i \le \alpha$, $\mathbf{B}_{0}^{(l+1)} = \sigma \left( \mathbf{H}^{(l)} \mathbf{W}_{0}^{(l)} \right)$, and $\mathbf{W}_{\mathrm{MLP}}^{(l)} \in \mathbb{R}^{(\alpha+1)H_{l+1} \times H_{l+1}}$ is the matrix of learnable weights to combine the $\alpha + 1$ branches.

It is worth noting that a layer of S2-GNN, as described in (12) or (13), can be viewed as a mapping $H: \mathbb{R}^{N \times H_l} \to \mathbb{R}^{N \times H_{l+1}}$. In practice, the choice of the fusion mechanism is determined by a hyperparameter optimization algorithm that selects the most suitable mechanism for a given dataset. Other techniques, like attention-based fusion methods [37] could also be applied into our model. However, such techniques could increase the overall complexity of S2-GNN due to the computation of the attention mechanism.

1) Re-normalization Trick and Sobolev Term: In practice, we utilize a normalized Sobolev version of the adjacency matrix as the shift operator for computing our propagation rule. Specifically, we set $\mathbf{S}_{(\rho)} = \bar{\mathbf{D}}_{\rho}^{-\frac{1}{2}}\bar{\mathbf{A}}_{\rho}\bar{\mathbf{D}}_{\rho}^{-\frac{1}{2}}$, where $\bar{\mathbf{A}}_{\rho} =$
<figure><img src="image_5.png" /><figcaption>Fig. 4: Basic configuration of our S2-GNN architecture with $n$ layers and $\alpha$ branches per layer.</figcaption></figure>

($\mathbf{A} + \epsilon\mathbf{I}$)$^{(\rho)}$ represents the $\rho$th sparse Sobolev term of $\mathbf{A}$ and $\bar{\mathbf{D}}_{\rho}$ denotes the degree matrix of $\bar{\mathbf{A}}_{\rho}$. This operation can be interpreted from a message-passing perspective, where we aggregate the embeddings from the 1-hop neighborhood with weights proportional to the edge attributes, node degrees, and the Hadamard power involved in that operation. Unlike the GCN operation, our approach incorporates the Hadamard power term. Our operation does not compute attention coefficients like GAT or graph Transformers, but it does learn the importance of each branch through the fusion layer, either by using a linear combination or an MLP operation.

Figure 4 illustrates the basic configuration of S2-GNN, serving as the core for constructing more complex architectures based on our filtering operation. Additional components such as batch normalization and residual connections can be added. Notice that our graph convolution is efficiently computed since the term $\bar{\mathbf{D}}_{\rho}^{-\frac{1}{2}}\bar{\mathbf{A}}_{\rho}\bar{\mathbf{D}}_{\rho}^{-\frac{1}{2}}$; $\forall \rho \in \{1, 2, \dots, \alpha\}$ remains the same across all layers (allowing for offline computation), and these terms remain sparse for any value of $\rho$ (given the sparsity of $\mathbf{A}$). S2-GNN can employ conventional activation functions like ReLU or softmax, as well as typical loss functions $\mathcal{L}$ suited to the specific task. The basic architecture of S2-GNN is defined by the number of filters $\alpha$ in each layer, the parameter $\epsilon$, the number of hidden units in each $\mathbf{W}_{\rho}^{(l)}$, and the number of layers $n$. When constructing weighted graphs with Gaussian kernels, the edge weights fall within the interval $[0, 1]$. As a result, large values of $\rho$ could lead to $\bar{\mathbf{A}}_{\rho} = 0$, and the diagonal elements of $\bar{\mathbf{D}}_{\rho}^{-\frac{1}{2}}$ may become $\infty$. Similarly, when $\alpha$ is very large, the resulting architectures become wide with a significant number of parameters, making it desirable to maintain a reasonable value for $\alpha$. In our experiments, we evaluate S2-GNN architectures up to $\alpha = 6$.

## D. Relationship with Graph Convolutional Network and Computational Complexity

Kipf and Welling [3] proposed one of the most successful yet simple GNNs, called Graph Convolutional Networks (GCNs), using the following propagation rule:

$$ \bar{\mathbf{H}}^{(l+1)} = \sigma(\tilde{\mathbf{D}}^{-\frac{1}{2}}\tilde{\mathbf{A}}\tilde{\mathbf{D}}^{-\frac{1}{2}}\bar{\mathbf{H}}^{(l)}\mathbf{W}^{(l)}). \quad (14) $$

Here $\tilde{\mathbf{A}} = \mathbf{A} + \mathbf{I}$, $\tilde{\mathbf{D}}$ represents the degree matrix of $\tilde{\mathbf{A}}$, $\bar{\mathbf{H}}^{(l)} \in \mathbb{R}^{N \times H_l}$ denotes the output of the previous layer $l$, where $\bar{\mathbf{H}}^{(1)} = \mathbf{X} \in \mathbb{R}^{N \times M}$ is matrix of input features, and

$\mathbf{W}^{(l)} \in \mathbb{R}^{H_l \times H_{l+1}}$ is the matrix of trainable weights in layer $l$. The motivation behind the propagation rule in (14) is based on the first-order approximation of localized spectral filters on graphs [2]. We can observe that GCN is an approximation of the graph convolutional filter in Definition 1, where $q_1 = 1$, $q_i = 0$; $\forall i \neq 1$, the input signal is an $H_l$-dimensional graph signal, and an MLP is added at the end of the filtering process. Notice that our S2-GNN layer is a generalization of the GCN layer when only the branch $\rho = 1$ is used in (11) without the fusion layer, and the shift operator is the normalized adjacency matrix.

The computational complexity of a GCN with $n$ layers is $O(n|\mathcal{E}|)$. The GCN model is widely used in practice due to its simplicity, good performance, and low computational complexity. For S2-GNN, the computational complexity of each branch remains $O(|\mathcal{E}|)$ since our sparse Sobolev term maintains the same sparsity level in the graph representation for any value of $\rho$. When considering $\alpha$ branches, the complexity becomes $O(\alpha|\mathcal{E}| + 2\alpha)$ per layer, taking into account the linear combination layer and the branch for $\rho = 0$. With $n$ layers, the computational complexity of our architecture is $O(n(\alpha|\mathcal{E}| + 2\alpha)) = O(n\alpha|\mathcal{E}| + 2n\alpha)$. Therefore, S2-GNN achieves a computational complexity similar to GCN while delivering better performance on several benchmark datasets. It is worth noting that the exact complexity of both methods, S2-GNN and GCN, also depends on factors such as the feature dimension, the hidden units, and the number of nodes in the graph, which have been omitted for simplicity and clarity.

## E. Theoretical Stability Analysis

Here, we aim to study the effect of graph perturbations on the prediction performance of the proposed network. Similarly, we analyze the stability [23] of S2-GNN against the existence of erroneous edges. Particularly, we consider the deviation of the adjacency matrix $\mathbf{A}$ with the error (or perturbation) matrix $\mathbf{E}$ as follows:

$$ \hat{\mathbf{A}} = \mathbf{A} + \mathbf{E}; \mathbf{E} \in \mathcal{P}, \quad (15) $$

$\mathcal{P} = \{\mathbf{E} \in \mathbb{R}^{N \times N} | \mathbf{E} = \mathbf{E}^{\top}, |e_{i,j}| \le a_{i,j}; \forall i, j, e_{i,i} = 0; \forall i\}.$

Note that the symmetry of the error matrix $\mathbf{E}$ stems from the underlying graph being undirected. The conditions on $\mathcal{P}$ imply that: 1) the perturbations cannot cause negative edge weights, and 2) we do not consider edge perturbations causing self-loops.
For the sake of simplicity and without loss of generality, we consider the true and perturbed output of a simplified [29] version of the S2-GNN (SS2-GNN), which can be stated as:

$$ 
\begin{aligned}
\mathbf{Y} &= \sigma(\mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W}), \\
\hat{\mathbf{Y}} &= \sigma(\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}}),
\end{aligned}
\qquad (16) $$

where $\hat{\mathbf{L}}$ is the perturbed Laplacian matrix. The next theorem obtains an upper bound on the distance between the true and perturbed outputs (i.e., $d(\mathbf{Y}, \hat{\mathbf{Y}}) = \|\mathbf{Y} - \hat{\mathbf{Y}}\|$) for the SS2-GNN in the case of having true and perturbed versions of the graph and graph filter coefficients. First, we need to mention the definition of being Lipschitz for a general function $\sigma(\cdot)$ in $\mathbb{R}$.

**Definition 4.** We call a function $\sigma(\cdot)$ as Lipschitz if there exists a positive constant $\phi$ such that:

$$ \forall x_1, x_2 \in \mathbb{R}: |\sigma(x_1) - \sigma(x_2)| \le \phi|x_1 - x_2|. \qquad (17) $$

**Theorem 3.** With the assumptions of the non-linear function $\sigma(\cdot)$ being Lipschitz, the node feature matrix $\mathbf{X}$ is column-wisely normalized as $\|\mathbf{X}\|_F = \sqrt{v}, \|\mathbf{W} - \hat{\mathbf{W}}\| \le \delta_W, \|\mathbf{E}\| \le \xi$, and $\|\mathbf{D}_E\| \le \xi_D$, where $\mathbf{D}_E$ is the degree matrix of $\mathbf{E}$, the difference between the true and perturbed outputs of an SS2-GNN is bounded by:

$$ 
\begin{aligned}
d(\mathbf{Y}, \hat{\mathbf{Y}}) &= \|\sigma(\mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W}) - \sigma(\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}})\| \\
&\le \sqrt{v} \left[ (\hat{d}^{\max} + \epsilon)^\rho + \|\hat{\mathbf{A}}^{(\rho)}\| \right] \delta_W + \mathcal{O}(\xi^2) + \mathcal{O}(\xi_D^2) \\
&\quad + \rho\sqrt{v} \left[ \eta \|\mathbf{E}\| + (d^{\max} + \epsilon)^{\rho-1} d_e^{\max} \right] \|\mathbf{W}\|,
\end{aligned}
\qquad (18) $$

where $\eta \triangleq \min(r_1(\mathbf{A}^{(\rho-1)}), c_1(\mathbf{A}^{(\rho-1)}))$, $r_1(\cdot)$ and $c_1(\cdot)$ denote the row and column of the involved matrix with highest norm-2 among the other rows and columns, respectively. Proof: See Appendix C.

The obtained bound in (18) reveals interesting aspects of the S2-GNN stability against possible perturbations. Precisely, this bound relies on the Hadamard order ($\rho$), which implies that a high degree of Hadamard multiplications could make the network more susceptible to perturbations. Besides, the absolute dependency of this bound on $\|\mathbf{E}\|, \xi$, and $\xi_D$, which are respectively representations of the perturbation (noise) power, depicts the effect of perturbation amplitudes on the network stability. The existence of $\|\mathbf{W}\|$ and $\delta_W$ emphasizes the choice of the network size and selection of optimization approach (e.g., stochastic gradient descent) for reaching the suboptimal solution for $\|\mathbf{W}\|$, which has a direct impact on the distance (i.e., $\delta_W$) from the true values of $\mathbf{W}$. Also, as can be seen in (18), an unnecessarily high value of $\epsilon$ could have considerably negative effects on the network robustness, especially in the case of high Hadamard order ($\rho$). Finally, due to the existence of $\eta$ and more importantly, $d^{\max}$, in the upper bound (18), the smaller values for $d^{\max}$ (i.e., sparse graphs) help in stability and robustness of the network. Therefore, the proposed network can effectively benefit from the sparse underlying structures, which are widely found in real-world

applications. We empirically validate the stability analysis of S2-GNN in Appendix D.

## V. EXPERIMENTAL EVALUATION

S2-GNN is compared to nine state-of-the-art methods in GNNs. These methods include Chebyshev graph convolutions (Cheby) [2], Graph Convolutional Networks (GCN) [3], Graph Attention Networks (GAT) [4], Simple Graph Convolution (SGC) [29], ClusterGCN [30], Scalable Inception Graph Neural network (SIGN) [10], SuperGAT [28], a graph Transformer [21], and GAT version 2 (GATv2) [20].

### A. Datasets

We evaluate the performance of S2-GNN on multiple datasets, which include tissue phenotyping in colon cancer histology images [15], text classification of news [16], activity recognition using sensors [17], recognition of spoken letters [18], as well as the commonly used node classification benchmarks Cora [38], Citeseer [39], Pubmed [40], and OGBN-proteins [19]. These diverse datasets provide a comprehensive evaluation of S2-GNN's performance across different domains. The **colorectal cancer histology** dataset comes from ten anonymized tissue slides from the University Medical Center Mannheim [15]. The dataset consists of eight classes: tumor epithelium, simple stroma, complex stroma, immune cells, debris, normal mucosal glands, adipose tissue, and background (no tissue). This work does not explore the feature representation, and only uses Local Binary Patterns (LBP) to represent each image. We explore two partitions of this dataset: 1) Cancer Binary (Cancer-B) considers only tumor epithelium and simple stroma, and 2) Cancer Multi-class (Cancer-M) tests all classes.

The **20 Newsgroups** dataset (20News) contains approximately 20,000 news articles categorized into 20 classes [16]. In this work, we use a subset of ten categories: computer graphics, windows operative system, IBM hardware, MAC hardware, autos, motorcycles, cryptography, electronics, medicine, and space. Each news document is represented using the Term Frequency - Inverse Document Frequency (TFIDF) as described in [41].

The **Human Activity Recognition** dataset (HAR) comprises recordings of 30 individuals engaged in various activities of daily living [17]. The participants wore waist-mounted smartphones with embedded sensors, and the dataset consists of six activity classes: standing, sitting, laying down, walking, walking downstairs, and walking upstairs. The recordings are represented as in [17].

The **Isolated Letter recognition** dataset (Isolet) includes sound recordings from the 26 letters of the English alphabet spoken by 30 speakers. Isolet contains 6,238 samples represented as 617-dimensional feature vectors. The feature representation method is explained in [18].

The **Cora**, **Citeseer**, and **Pubmed** datasets are well-known citation networks. The datasets consist of sparse bag-of-words feature vectors for each document (node) and a list of citation links between documents. In this work, we treat the citation links as undirected edges as in [3]. To fulfill the requirement of
<table><caption>TABLE I: Statistics of the datasets in the experimental framework of this work.</caption><thead><tr><th></th><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th><th>Cora</th><th>Citeseer</th><th>Pubmed</th><th>OGBN-proteins</th></tr></thead><tbody><tr><th>Nodes</th><td>1, 250</td><td>5,000</td><td>3,000</td><td>3,000</td><td>3,000</td><td>2,485</td><td>2,120</td><td>19,717</td><td>30,000</td></tr><tr><th>Edges*</th><td>25,734</td><td>101,439</td><td>62,399</td><td>87,032</td><td>43,330</td><td>205,928</td><td>183,572</td><td>1,723,360</td><td>13,251,798</td></tr><tr><th>Edges†</th><td>19,163</td><td>81,462</td><td>47,496</td><td>61,400</td><td>31,130</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><th>Features</th><td>38</td><td>38</td><td>2,000</td><td>561</td><td>617</td><td>1,433</td><td>3,703</td><td>500</td><td>8</td></tr><tr><th>Classes</th><td>2</td><td>8</td><td>10</td><td>6</td><td>26</td><td>7</td><td>6</td><td>3</td><td>112</td></tr><tr><th>Undirected</th><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

weighted graphs, as mentioned in Section III, we preprocess the *unweighted* citation networks using the Graph Diffusion Convolution (GDC) method [42].

**The Open Graph Benchmark Proteins** dataset (OGBN-Proteins) is a comprehensive graph-based dataset designed for protein-protein interaction networks [19]. The dataset features more than 132,000 nodes (proteins) and 39 million edges, reflecting a diverse range of biological relationships. Each node in the graph is associated with 8-dimensional feature vectors, capturing essential biological characteristics of the proteins. The primary task of the dataset is to predict protein functions, with each protein potentially linked to one or more of the 112 function categories. Due to computational constraints, the scope of the data used in the experiments is reduced to approximately 25% of the full dataset, which equates to around 30,000 nodes. This subset selection ensures manageable computational demands while still providing a representative sample for analysis.

## B. Experiments

For the Cancer-B, Cancer-M, 20News, HAR, and Isolet datasets, we generate graphs using two methods: *k*-Nearest Neighbors (*k*-NN) and the algorithm described in [25]. In our results, we refer to the graphs generated by the latter method as *learned graphs*. Regarding Cora, Citeseer, and Pubmed datasets, we consider the largest connected component of the graph following common practices [42]. OGBN-proteins is a dataset that comprises a single weighted graph, where edges represent the confidence level of interactions between proteins, eliminating the necessity for preprocessing steps like GDC. This simplifies the data preparation phase and allows for a more direct application of the model. Table I presents a comprehensive summary of the statistics for all the graphs used in the study. To optimize the hyperparameters of citation networks, we employ a random search procedure that maximizes the average accuracy on the validation sets. For OGBN-proteins, we empirically obtain the hyperparameters, focusing on the critical importance of learning rate and hidden channels over other hyperparameters. We split the datasets into training, validation, and testing sets. Initially, we divide the data into a development set and a test set, ensuring that the test set is not used during the hyperparameter optimization process. Following the experimental framework of [42], we optimize the hyperparameters for each dataset-preprocessing combination separately using a random search over five data splits, with multiple repetitions for enhanced robustness. For

the Cora, Citeseer, and Pubmed datasets, the development set consists of 1,500 nodes, while the remaining nodes are used for testing. Similarly, the train set includes 20 nodes from each class, and the remaining nodes are utilized for validation. As for the other datasets, we employ a 10/45/45 split, meaning that 10% of the nodes are assigned for training, 45% for validation, and 45% for testing. Finally, we report average accuracies on the test set accompanied by 95% confidence intervals calculated via bootstrapping with 1,000 samples, with each model being tested using 50 different seeds to ensure robustness and reliability of the results.

## C. Implementation Details

S2-GNN and the state-of-the-art methods are implemented using PyTorch Geometric (PyG) [43]. The learned graph methods are implemented using the GSP toolbox [44]. The value of *k* in the *k*-NN method is set to 30. Rectified Linear Unit (ReLU) and log-softmax are used as activation functions in all methods. The hyperparameters such as the number of GNN layers, hidden units, learning rate, weight decay, and dropout rate are optimized for each method using the validation set in each experiment. For Cheby we optimize the Chebyshev filter size *K*. For SIGN we optimize the number of powers. For attention-based techniques like GAT and Transformers, we optimize the number of heads. For S2-GNN we optimize $\alpha$, $\epsilon$, and whether the fusion mechanism is a linear combination layer or an MLP fusion as introduced in (12) or (13). The hyperparameter search space for each method is defined as follows: 1) learning rate $lr \in [0.005, 0.02]$, 2) weight decay $wd \in [0.0001, 0.001]$, 3) hidden units of each graph convolutional layer $hu \in \{16, 32, 64\}$, 4) dropout $d \in [0.3, 0.7]$, 5) the number of layers $L \in \{2, 3, 4, 5\}$, 6) Chebyshev filter size $K \in \{1, 2, 3\}$, 7) heads of attention mechanism $hd \in \{1, 2, \dots, 6\}$, 8) powers for SIGN $pw \in \{1, 2, 3\}$, 9) $\alpha \in \{1, 2, \dots, 6\}$, 10) $\epsilon \in [0.5, 2]$, 11) fusion layer $fl \in \{\text{Linear}, \text{MLP}\}$. For the GDC method applied to Cora, Citeseer, and Pubmed datasets, personalized PageRank is used with a teleport probability of 0.05, and the top 128 entries with the highest mass per column are selected to sparsify the output of GDC. The code and constructed graphs will be publicly available upon publication under the MIT license<fnref n="5" />.

## D. Results and Discussion

Table II provides a summary of the comparison between S2-GNN and state-of-the-art methods in terms of Accuracy for

<footnote>* This number of edges corresponds to the *k*-NN method for the constructed graphs, and the outputs of the GDC method for the citation networks.</footnote>

<footnote>† This number of edges corresponds to the method in [25] for the constructed graphs.</footnote>

<footnote n="5">⁵https://github.com/jhonygiraldo/S2-GNN</footnote>
<table><caption>TABLE II: Accuracy (in %) for the state-of-the-art methods and our S2-GNN algorithm in five datasets, inferring the graphs with *k*-NN or the method for learning smooth graphs in [25].</caption><thead><tr><th rowspan="2">Model</th><th colspan="5"><i>k</i>-NN Graphs</th><th colspan="5">Learned Graphs</th></tr><tr><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th></tr></thead><tbody><tr><td>Cheby [2]</td><td>82.55±4.95</td><td>59.52±1.42</td><td>69.53±1.47</td><td>85.64±2.92</td><td>71.52±1.42</td><td>87.64±3.69</td><td>57.50±2.12</td><td><b>72.49</b>±0.82</td><td>87.07±3.57</td><td>59.14±2.50</td></tr><tr><td>GCN [3]</td><td>71.31±5.72</td><td>50.87±2.39</td><td>63.08±2.36</td><td>64.28±4.90</td><td>63.61±1.76</td><td>80.50±4.77</td><td>52.45±3.00</td><td>69.83±1.25</td><td>82.85±3.00</td><td>71.27±2.01</td></tr><tr><td>GAT [4]</td><td>68.67±5.24</td><td>48.87±1.94</td><td>48.90±2.17</td><td>70.13±3.65</td><td>58.93±2.03</td><td>77.69±4.52</td><td>51.23±2.04</td><td>54.30±2.13</td><td>75.45±3.97</td><td>58.78±1.75</td></tr><tr><td>SGC [29]</td><td>72.00±4.56</td><td>42.22±1.88</td><td>54.24±2.40</td><td>39.39±2.72</td><td>39.90±2.05</td><td>70.82±4.87</td><td>43.97±1.97</td><td>72.13±1.31</td><td>38.22±2.71</td><td>43.21±1.85</td></tr><tr><td>ClusterGCN [30]</td><td>83.03±4.43</td><td>56.85±1.25</td><td>67.95±1.55</td><td>72.67±4.70</td><td>59.05±3.20</td><td>85.63±4.26</td><td>58.58±1.76</td><td>62.12±1.59</td><td>65.98±4.45</td><td>69.54±2.37</td></tr><tr><td>SIGN [10]</td><td><b>89.67</b>±0.39</td><td><b>62.45</b>±0.17</td><td><b>72.34</b>±0.25</td><td><b>91.58</b>±0.24</td><td><b>84.11</b>±0.24</td><td><b>90.30</b>±0.38</td><td><b>64.99</b>±0.20</td><td><b>74.10</b>±0.30</td><td><b>92.88</b>±0.23</td><td><b>84.17</b>±0.27</td></tr><tr><td>SuperGAT [28]</td><td>69.54±5.09</td><td>40.04±1.84</td><td>54.40±1.95</td><td>72.66±3.76</td><td>60.57±2.17</td><td>78.15±4.42</td><td>52.73±2.07</td><td>59.53±1.96</td><td>74.93±4.00</td><td>61.49±1.68</td></tr><tr><td>Transformer [21]</td><td>74.23±5.65</td><td>54.20±1.70</td><td>59.06±2.75</td><td>79.66±3.34</td><td>65.68±2.42</td><td>79.23±4.83</td><td>52.07±2.21</td><td>56.19±2.83</td><td>78.30±2.81</td><td>65.36±2.18</td></tr><tr><td>GATv2 [20]</td><td>62.89±5.21</td><td>47.15±2.07</td><td>46.56±2.04</td><td>73.28±4.02</td><td>64.59±1.40</td><td>79.28±4.01</td><td>45.80±1.78</td><td>51.17±2.14</td><td>75.09±4.15</td><td>66.46±1.87</td></tr><tr><td>S2-GNN (ours)</td><td><b>93.39</b>±0.41</td><td><b>68.48</b>±0.28</td><td><b>72.09</b>±0.35</td><td><b>93.92</b>±0.22</td><td><b>86.57</b>±0.27</td><td><b>93.50</b>±0.39</td><td><b>68.88</b>±0.25</td><td>72.19±0.41</td><td><b>94.50</b>±0.19</td><td><b>86.61</b>±0.29</td></tr></tbody></table>

The best and second-best performing methods on each dataset are shown in red and blue, respectively.

<table><caption>TABLE III: Accuracy (in %) comparison in the citation networks using the graph diffusion convolution and OGBN-proteins.</caption><thead><tr><th>Model</th><th>Cora</th><th>Citeseer</th><th>Pubmed</th><th>OGBN-proteins</th></tr></thead><tbody><tr><td>Cheby [2]</td><td><b>80.79</b>±1.43</td><td><b>69.41</b>±1.09</td><td><b>78.25</b>±1.57</td><td><b>63.55</b>±2.89</td></tr><tr><td>GCN [3]</td><td>70.22±2.91</td><td>58.85±2.88</td><td>56.64±5.81</td><td>58.86±7.63</td></tr><tr><td>GAT [4]</td><td>60.29±4.25</td><td>56.18±2.76</td><td>56.57±5.11</td><td>41.58±4.42</td></tr><tr><td>SGC [29]</td><td>64.23±3.66</td><td>56.74±2.82</td><td>70.35±2.89</td><td>56.87±6.42</td></tr><tr><td>ClusterGCN [30]</td><td>75.95±2.53</td><td>63.93±1.96</td><td>59.24±6.31</td><td>50.01±0.02</td></tr><tr><td>SuperGAT [28]</td><td>65.71±3.69</td><td>55.80±3.27</td><td><b>OOM</b></td><td>58.59±0.79</td></tr><tr><td>Transformer [21]</td><td>66.86±3.47</td><td>62.85±2.19</td><td>61.83±5.46</td><td>39.75±4.76</td></tr><tr><td>GATv2 [20]</td><td>62.16±3.29</td><td>53.24±2.25</td><td>56.39±5.53</td><td>50.06±0.13</td></tr><tr><td>S2-GNN (ours)</td><td><b>82.47</b>±0.42</td><td><b>68.95</b>±0.55</td><td><b>77.64</b>±1.05</td><td><b>70.26</b>±4.08</td></tr></tbody></table>

OOM: Out Of Memory in a GPU Nvidia A100 40GB.

<table><caption>TABLE IV: Sparsity (%) for the non-sparse Sobolev term $(L + \epsilon I)^\rho$ for several values of $\rho$.</caption><thead><tr><th>Power</th><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th></tr></thead><tbody><tr><td>ρ = 1</td><td>96.63</td><td>99.17</td><td>98.58</td><td>98.03</td><td>99.01</td></tr><tr><td>ρ = 2</td><td>80.33</td><td>94.05</td><td>65.65</td><td>82.29</td><td>92.73</td></tr><tr><td>ρ = 3</td><td>53.82</td><td>82.48</td><td>0.615</td><td>59.74</td><td>75.92</td></tr><tr><td>ρ = 4</td><td>29.32</td><td>67.67</td><td>0</td><td>46.05</td><td>46.19</td></tr><tr><td>ρ = 5</td><td>12.46</td><td>52.55</td><td>0</td><td>35.07</td><td>20.41</td></tr></tbody></table>

Zero sparsity percentage means that the matrix is completely dense.

Cancer-B, Cancer-M, 20News, HAR, and Isolet. In most cases, S2-GNN demonstrates the best performance, except for the 20News dataset. Both Cheby and SIGN methods, which rely on higher-order graph convolutions, also achieve competitive results. However, the disadvantage of these non-sparse higher-order methods is that the graph matrix representations (either the adjacency or Laplacian matrices) become denser with higher powers of the shift operator. S2-GNN achieves better performances while maintaining sparsity in the graph matrix representations and providing a good confidence interval in the results. Furthermore, Table II highlights that learned graphs on the right side generally outperform the *k*-NN graphs on the left for most methods (Appendix E presents a study in the homophily level of *k*-NN and learned graphs). The performance of the learned graphs is consistently better, with only a few specific cases where the *k*-NN graphs perform better. However, these instances are rare exceptions.

Table III provides the results of S2-GNN and state-of-the-art methods on the citation network datasets, after being pre-processed with the GDC method, and the OGBN-proteins

dataset. S2-GNN achieves the best performance in Cora and OGBN-proteins datasets and very competitive performance in Citeseer and Pubmed datasets. Tables II and III collectively demonstrate that S2-GNN exhibits strong performance in cases where we construct graphs for the downstream task (using either *k*-NN or smooth graphs) as well as when the graph is naturally provided by the task, as seen in the citation networks. We provide additional experiments on graph classification tasks in the Appendix F.

Table IV shows the sparsity percentage of the non-sparse Sobolev term $(L + \epsilon I)^\rho$ for different values of $\rho$. We observe that the non-sparse filter loses sparsity quickly for higher values of $\rho$, resulting in increased memory and computational complexity since more edges are considered for message passing.

**Memory consumption and inference time.** We thoroughly analyze the memory consumption and inference time of S2-GNN with synthetic graphs. To this end, we first build underlying graphs using Erdős–Rényi (ER) models with varying number of nodes $N \in \{500, 1000, 5000, 10000\}$ and edge probabilities $p_{ER} \in \{0.03, 0.04, 0.05, 0.06\}$. Therefore, by feeding an identical initial feature matrix $X \in \mathbb{R}^{N \times F}$, where $F = 16$, we plotted the averaged memory usage (in terms of GB) and inference times (in seconds) over 30 random seeds in Figure 5. We observe that S2-GNN is among the most efficient models in terms of memory and time complexity when the number of nodes increases and the sparsity of the graphs decreases.

### E. Ablation Study

We conduct an ablation study to investigate whether the performance improvement in S2-GNN is attributed to the additional parameters introduced by having multiple branches or the Hadamard operation in the graph filter. To explore this, we create a similar model to S2-GNN but exclude the Hadamard powers in each branch. Instead, we use the identity operation in the first branch and the term $\tilde{D}^{-\frac{1}{2}} \tilde{A} \tilde{D}^{-\frac{1}{2}}$ in the remaining branches. By performing hyperparameter optimization as outlined in Section V-B, we obtain the results summarized in Table V. The findings suggest that, for nearly all datasets, incorporating the Hadamard powers using either *k*-NN or learned graphs leads to better performance compared to simply adding extra parameters to the GNN architecture.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Comparison of inference times (in seconds) and memory consumption (in terms of GB) on the underlying ER graphs with varying number of nodes $N \in \{500, 1000, 5000, 10000\}$ and edge probabilities $p_{ER} \in \{0.03, 0.04, 0.05, 0.06\}$.</figcaption></figure>

<table><caption>TABLE V: Ablation study in the selection of Hadamard power for the filtering operation.</caption><thead><tr><th rowspan="2">Hadamard Power</th><th colspan="5">k-NN Graphs</th><th colspan="5">Learned Graphs</th></tr><tr><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th></tr></thead><tbody><tr><td>✗</td><td>91.20±0.50</td><td>66.29±0.35</td><td>72.08±0.37</td><td>93.08±0.34</td><td>85.66±0.31</td><td>88.49±3.67</td><td>67.97±0.24</td><td><strong>72.63±0.31</strong></td><td>93.82±0.25</td><td>86.44±0.32</td></tr><tr><td>✓</td><td><strong>93.39</strong>±0.41</td><td><strong>68.48</strong>±0.28</td><td><strong>72.09</strong>±0.35</td><td><strong>93.92</strong>±0.22</td><td><strong>86.57</strong>±0.27</td><td><strong>93.50</strong>±0.39</td><td><strong>68.88</strong>±0.25</td><td>72.19±0.41</td><td><strong>94.50</strong>±0.19</td><td><strong>86.61</strong>±0.29</td></tr></tbody></table>

The best results on each dataset are shown in **bold**.

<table><caption>TABLE VI: Ablation study regarding regular versus sparse Sobolev norm.</caption><thead><tr><th rowspan="2">Sparse Norm</th><th colspan="5">k-NN Graphs</th><th colspan="5">Learned Graphs</th></tr><tr><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th><th>Cancer-B</th><th>Cancer-M</th><th>20News</th><th>HAR</th><th>Isolet</th></tr></thead><tbody><tr><td>✗</td><td>87.05±4.25</td><td>64.75±0.34</td><td><strong>72.29</strong>±0.34</td><td>92.80±0.30</td><td>85.81±0.35</td><td>92.92±0.48</td><td>66.41±0.37</td><td><strong>73.02</strong>±0.41</td><td>93.32±0.35</td><td>86.25±0.30</td></tr><tr><td>✓</td><td><strong>93.39</strong>±0.41</td><td><strong>68.48</strong>±0.28</td><td>72.09±0.35</td><td><strong>93.92</strong>±0.22</td><td><strong>86.57</strong>±0.27</td><td><strong>93.50</strong>±0.39</td><td><strong>68.88</strong>±0.25</td><td>72.19±0.41</td><td><strong>94.50</strong>±0.19</td><td><strong>86.61</strong>±0.29</td></tr></tbody></table>

This implies that the Hadamard operation plays a crucial role in enhancing the model's effectiveness.

Another ablation experiment studies the usage of the regular and sparse Sobolev norms in S2-GNN. Table VI shows the comparison results, where we observe the sparse norm provides better results on all datasets except on 20News.

**Sensibility analysis of $\epsilon$.** Figure 6 shows a sensibility analysis of the hyperparameter $\epsilon \in [0, 2]$ (with standard deviations) for different datasets using k-NN graphs. We observe that $\epsilon$ has a relatively important impact on the Cancer-M dataset. For the rest of the datasets, S2-GNN is relatively robust to the choice of $\epsilon$. We note that the sparsity level of the Cancer-M dataset is the highest among the compared datasets. Therefore, there is a higher chance for the sparsest graph of Cancer-M to have an ill-condition Laplacian and be improved by the Sobolev coefficient $\epsilon$, as discussed in (6).

## F. Limitations

One of the main limitations of S2-GNN is its reliance on weighted graphs as input when using the adjacency matrix. This requirement can pose a limitation in certain scenarios where only unweighted graphs are available or when the weighted graph construction process is challenging or unreliable. Additionally, compared to traditional GCN, S2-GNN introduces three additional hyperparameters: $\alpha$, $\epsilon$, and the fusion layer. The inclusion of these hyperparameters increases the complexity of the model and makes the hyperparameter optimization process more challenging compared to GCN. Finding an optimal set of hyperparameters for S2-GNN through random search becomes more difficult due to the expanded search space. These limitations emphasize the importance of careful consideration and understanding of the graph structure and hyperparameter selection when applying
<figure><img src="image_7.png" /><figcaption>Fig. 6: Mean and standard deviation for the accuracy in different datasets for the sensibility of the hyperparameter $\epsilon$.</figcaption></figure>

S2-GNN in practical scenarios.

## VI. CONCLUSION

In this work, we introduced a new higher-order sparse graph filter for GNNs. We proposed a novel sparse Sobolev norm inspired by the GSP literature and used Hadamard powers to compute efficient graph filters. Sections IV-A, IV-B, and IV-E provide theoretical insights into our filtering operation. To capture higher-order dependencies in GNNs, we introduced the S2-GNN model. This model incorporates a cascade of increasing sparse filtering operations on each layer, along with a fusion layer that selects the most effective operations. By leveraging these mechanisms, S2-GNN effectively captures and utilizes higher-order information in graph-structured data. In addition, we explored the inference of smooth graphs as an alternative when the underlying graph topology is unavailable for a given problem. Our observations indicate that smooth inferred graphs improve the performance of baseline GNN architectures by promoting homophily. We thoroughly evaluated S2-GNN in several datasets in a semi-supervised regimen. S2-GNN achieves competitive performance compared to state-of-the-art GNNs across all tasks while maintaining considerably lower complexity than the best-performing methods.

This work opens several research directions. The first direction is the extension of S2-GNN to the setting of ClusterGCN [30]. This extension would explore the possibility of removing the explicit need for inter-cluster communication, potentially simplifying the model architecture. Another promising avenue is the study of shift operators that do not require explicit weighted graphs. This direction would address the limitation of S2-GNN, which requires weighted graphs as input. By developing alternative shift operators, we can broaden the applicability of the model and overcome this constraint. Furthermore, neural architecture search can be studied in the context of S2-GNN to approach other problems such as graph classification and link prediction, while addressing the challenge of hyperparameter search in S2-GNN.

## REFERENCES

[1] J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun, “Spectral networks and locally connected networks on graphs,” in *ICLR*, 2014.

[2] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in *NeurIPS*, 2016.

[3] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in *ICLR*, 2017.

[4] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio, “Graph attention networks,” in *ICLR*, 2018.

[5] A. Duval and F. Malliaros, “Higher-order clustering and pooling for graph neural networks,” in *ACM CIKM*, 2022.

[6] G. Li, M. Muller, A. Thabet, and B. Ghanem, “DeepGCNs: Can GCNs go as deep as CNNs?” in *IEEE/CVF ICCV*, 2019.

[7] A. Benamira, B. Devillers, E. Lesot, A. K. Ray, M. Saadi, and F. D. Malliaros, “Semi-supervised learning and graph neural networks for fake news detection,” in *IEEE/ACM ASONAM*, 2019.

[8] A. Duval, V. Schmidt, A. H. Garcia, S. Miret, F. D. Malliaros, Y. Bengio, and D. Rolnick, “FAENet: frame averaging equivariant GNN for materials modeling,” in *ICML*, 2023.

[9] Z. Zhang, P. Cui, and W. Zhu, “Deep learning on graphs: A survey,” *IEEE T-KDE*, vol. 34, no. 1, pp. 249–270, 2022.

[10] F. Frasca, E. Rossi, D. Eynard, B. Chamberlain, M. M. Bronstein, and F. Monti, “SIGN: Scalable inception graph neural networks,” in *ICML-W*, 2020.

[11] F. R. K. Chung, *Spectral graph theory*. American Mathematical Society, 1997, no. 92.

[12] R. A. Horn and C. R. Johnson, *Matrix analysis*. Cambridge university press, 2012.

[13] I. Pesenson, “Variational splines and Paley–Wiener spaces on combinatorial graphs,” *Constructive Approximation*, vol. 29, no. 1, pp. 1–21, 2009.

[14] J. H. Giraldo, A. Mahmood, B. Garcia-Garcia, D. Thanou, and T. Bouwmans, “Reconstruction of time-varying graph signals via Sobolev smoothness,” *IEEE T-SIPN*, vol. 8, pp. 201–214, 2022.

[15] J. N. Kather, C.-A. Weis, F. Bianconi, S. M. Melchers, L. R. Schad, T. Gaiser, A. Marx, and F. G. Zöllner, “Multi-class texture analysis in colorectal cancer histology,” *Scientific Reports*, vol. 6, p. 27988, 2016.

[16] K. Lang, “NewsWeeder: Learning to filter netnews,” in *JMLR*, 1995.

[17] D. Anguita et al., “A public domain dataset for human activity recognition using smartphones,” in *ESANN*, 2013.

[18] M. Fanty and R. Cole, “Spoken letter recognition,” in *NeurIPS*, 1991.

[19] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec, “Open graph benchmark: Datasets for machine learning on graphs,” in *NeurIPS*, 2020.

[20] S. Brody, U. Alon, and E. Yahav, “How attentive are graph attention networks?” in *ICLR*, 2022.

[21] Y. Shi, Z. Huang, S. Feng, H. Zhong, W. Wang, and Y. Sun, “Masked label prediction: Unified message passing model for semi-supervised classification,” in *IJCAI*, 2021.

[22] J. H. Giraldo, S. Javed, A. Mahmood, F. D. Malliaros, and T. Bouwmans, “Higher-order sparse convolutions in graph neural networks,” in *IEEE ICASSP*, 2023.

[23] F. Gama, J. Bruna, and A. Ribeiro, “Stability properties of graph neural networks,” *IEEE T-SP*, vol. 68, pp. 5680–5695, 2020.

[24] A. Parada-Mayorga, Z. Wang, F. Gama, and A. Ribeiro, “Stability of aggregation graph neural networks,” *IEEE T-SIPN*, 2023.

[25] V. Kalofolias and N. Perraudin, “Large scale graph learning from smooth signals,” in *ICLR*, 2019.

[26] X. Pu, S. L. Chau, X. Dong, and D. Sejdinovic, “Kernel-based graph learning from smooth signals: A functional viewpoint,” *IEEE T-SIPN*, vol. 7, pp. 192–207, 2021.

[27] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, “Geometric deep learning: going beyond euclidean data,” *IEEE SPM*, vol. 34, no. 4, pp. 18–42, 2017.

[28] D. Kim and A. Oh, “How to find your friendly neighborhood: Graph attention design with self-supervision,” in *ICLR*, 2021.

[29] F. Wu, A. Souza, T. Zhang, C. Fifty, T. Yu, and K. Weinberger, “Simplifying graph convolutional networks,” in *ICML*, 2019.

[30] W.-L. Chiang, X. Liu, S. Si, Y. Li, S. Bengio, and C.-J. Hsieh, “Cluster-GCN: An efficient algorithm for training deep and large graph convolutional networks,” in *ACM SIGKDD*, 2019.

[31] D. Kreuzer, D. Beaini, W. Hamilton, V. Létourneau, and P. Tossou, “Rethinking graph transformers with spectral attention,” in *NeurIPS*, 2021.
[32] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, "Graph signal processing: Overview, challenges, and applications," *Proceedings of the IEEE*, vol. 106, no. 5, pp. 808–828, 2018.

[33] A. Sandryhaila and J. M. Moura, "Discrete signal processing on graphs," *IEEE T-SP*, vol. 61, no. 7, pp. 1644–1656, 2013.

[34] J. H. Giraldo, S. Javed, and T. Bouwmans, "Graph moving object segmentation," *IEEE T-PAMI*, vol. 44, no. 5, pp. 2485–2503, 2022.

[35] J. H. Giraldo and T. Bouwmans, "On the minimization of Sobolev norms of time-varying graph signals: Estimation of new Coronavirus disease 2019 cases," in *IEEE MLSP*, 2020.

[36] B. Nica, *A brief introduction to spectral graph theory*. European Mathematical Society Publishing House, 2018.

[37] Z. Peng, H. Liu, Y. Jia, and J. Hou, "Attention-driven graph clustering network," in *ACM MM*, 2021.

[38] A. K. McCallum, K. Nigam, J. Rennie, and K. Seymore, "Automating the construction of internet portals with machine learning," *Information Retrieval*, vol. 3, no. 2, pp. 127–163, 2000.

[39] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad, "Collective classification in network data," *AI Magazine*, vol. 29, no. 3, pp. 93–93, 2008.

[40] G. Namata, B. London, L. Getoor, B. Huang, and U. Edu, "Query-driven active surveying for collective classification," in *International Workshop on Mining and Learning with Graphs*, 2012.

[41] A. Gadde, A. Anis, and A. Ortega, "Active semi-supervised learning using sampling theory for graph signals," in *ACM KDD*, 2014.

[42] J. Klicpera, S. Weißenberger, and S. Günnemann, "Diffusion improves graph learning," in *NeurIPS*, 2019.

[43] M. Fey and J. E. Lenssen, "Fast graph representation learning with PyTorch Geometric," in *ICLR-W*, 2019.

[44] N. Perraudin, J. Paratte, D. Shuman, L. Martin, V. Kalofolias, P. Vandergheynst, and D. K. Hammond, "GSPBOX: A toolbox for signal processing on graphs," *arXiv preprint arXiv:1408.5781*, 2014.

[45] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains," *IEEE SPM*, vol. 30, no. 3, pp. 83–98, 2013.

[46] G. Visick, "A quantitative version of the observation that the Hadamard product is a principal submatrix of the Kronecker product," *LAA*, vol. 304, no. 1-3, pp. 45–68, 2000.

[47] X. Zhan, "Inequalities for the singular values of Hadamard products," *SIAM Journal on Matrix Analysis and Applications*, vol. 18, no. 4, pp. 1093–1095, 1997.

[48] H. Pei, B. Wei, K. C.-C. Chang, Y. Lei, and B. Yang, "Geom-GCN: Geometric graph convolutional networks," in *ICLR*, 2020.

[49] C. Morris, N. M. Kriege, F. Bause, K. Kersting, P. Mutzel, and M. Neumann, "TUDataset: A collection of benchmark datasets for learning with graphs," in *ICML-W*, 2020.
## APPENDIX A

### PROOF OF THE THEOREM 1

*Proof.* 1) For the triangle inequality we have that:

$$ \|\mathbf{x} + \mathbf{y}\|_{(\rho), \epsilon}^2 = \langle \mathbf{x} + \mathbf{y}, \mathbf{y} + \mathbf{x} \rangle_{(\rho), \epsilon} = \|(\mathbf{L} + \epsilon\mathbf{I})^{(\rho/2)}(\mathbf{x} + \mathbf{y})\|^2 \\ \rightarrow \|\mathbf{x} + \mathbf{y}\|_{(\rho), \epsilon}^2 = \|\mathbf{x}\|_{(\rho), \epsilon}^2 + \|\mathbf{y}\|_{(\rho), \epsilon}^2 + 2\langle\mathbf{x}, \mathbf{y}\rangle_{(\rho), \epsilon}. \quad (19) $$

Using the Cauchy–Schwarz inequality we obtain:

$$ \|\mathbf{x} + \mathbf{y}\|_{(\rho), \epsilon}^2 \le \|\mathbf{x}\|_{(\rho), \epsilon}^2 + 2\|\mathbf{x}\|_{(\rho), \epsilon}\|\mathbf{y}\|_{(\rho), \epsilon} + \|\mathbf{y}\|_{(\rho), \epsilon}^2 \\ \rightarrow \|\mathbf{x} + \mathbf{y}\|_{(\rho), \epsilon}^2 \le (\|\mathbf{x}\|_{(\rho), \epsilon} + \|\mathbf{y}\|_{(\rho), \epsilon})^2. \quad (20) $$

2) Let $s$ be a scalar, then for the absolute homogeneity property we have:

$$ \|s\mathbf{x}\|_{(\rho), \epsilon} = \|(\mathbf{L} + \epsilon\mathbf{I})^{(\rho/2)}s\mathbf{x}\| = |s|\|(\mathbf{L} + \epsilon\mathbf{I})^{(\rho/2)}\mathbf{x}\| = |s|\|\mathbf{x}\|_{(\rho), \epsilon}. \quad (21) $$

3) For the positive definiteness property, we rely on the Schur product theorem, which states that the Hadamard product between two positive definite matrices is also positive definite [12]. When $\epsilon = 0$ and $\rho = 1$, the sparse Sobolev norm reduces to the graph Laplacian quadratic form [45]. In this case, $\|\mathbf{x}\|_{(1),0}^2 = \mathbf{x}^\mathrm{T}\mathbf{L}\mathbf{x} = 0 \iff \mathbf{x} = \tau\mathbf{1}$ (with $\tau$ some constant value). As a result, $\|\mathbf{x}\|_{(\rho),\epsilon}$ becomes a semi-norm because it satisfies only the first two conditions of vector norms, i.e., the triangle inequality and absolute homogeneity. However, for $\epsilon > 0$, we have $\mathbf{L} + \epsilon\mathbf{I} > 0$ according to (6), and then $(\mathbf{L} + \epsilon\mathbf{I})^{(\rho)} > 0$ by using the Schur product theorem. As a result:

$$ \|\mathbf{x}\|_{(\rho), \epsilon}^2 = \mathbf{x}^\mathrm{T}(\mathbf{L} + \epsilon\mathbf{I})^{(\rho)}\mathbf{x} = 0 \iff \mathbf{x} = 0; \forall \epsilon > 0. \quad (22) $$

As a consequence, the sparse Sobolev norm satisfies the properties of vector norms for all $\epsilon > 0$ according to (20), (21), and (22). Similarly, the sparse Sobolev norm definition satisfies the properties of semi-norms for $\epsilon = 0$ according to (20) and (21). □

## APPENDIX B

### PROOF OF THEOREM 2

*Proof.* We can compute the spectral decomposition of the Kronecker product between matrices using concepts from product graphs. The spectrum of the Kronecker product of multiple graphs is obtained by taking the Kronecker products of their respective spectral components. For instance, the spectral decomposition of $\mathbf{L} \otimes \mathbf{L}$ can be expressed as:

$$ \mathbf{L} \otimes \mathbf{L} = (\mathbf{U} \otimes \mathbf{U})(\mathbf{\Lambda} \otimes \mathbf{\Lambda})(\mathbf{U}^{\mathrm{T}} \otimes \mathbf{U}^{\mathrm{T}}). \quad (23) $$

Similarly, there is a relationship between the Kronecker and the Hadamard products as follows:

**Theorem 4.** (Visick [46]) For $\mathbf{S}, \mathbf{T} \in \mathbb{R}^{n \times m}$, we have:

$$ \mathbf{S} \circ \mathbf{T} = \mathbf{P}_n^\mathrm{T}(\mathbf{S} \otimes \mathbf{T})\mathbf{P}_m, \quad (24) $$

where $\mathbf{P}_n \in \{0, 1\}^{n^2 \times n}$ and $\mathbf{P}_m \in \{0, 1\}^{m^2 \times m}$ are partial permutation matrices. If $\mathbf{S}, \mathbf{T} \in \mathbb{R}^{n \times n}$ are square matrices, we have that $\mathbf{S} \circ \mathbf{T} = \mathbf{P}_n^\mathrm{T}(\mathbf{S} \otimes \mathbf{T})\mathbf{P}_m$.

*Proof:* see [46].

Therefore, by using (23) and Theorem 4, we can derive a general form of the spectrum of the Hadamard product (specifically, for the Hadamard power of order 2) as follows:

$$ \mathbf{L} \circ \mathbf{L} = \mathbf{L}^{(2)} = \mathbf{P}_N^\mathrm{T}(\mathbf{U} \otimes \mathbf{U})(\mathbf{\Lambda} \otimes \mathbf{\Lambda})(\mathbf{U}^\mathrm{T} \otimes \mathbf{U}^\mathrm{T})\mathbf{P}_N. \quad (25) $$

We can then conclude that the spectrum of the Hadamard power of order $\rho$ is such that:

$$ \mathbf{L}^{(\rho)} = \mathbf{P}_{(\rho)}^{\mathrm{T}} (\mathbf{U} \otimes \mathbf{U}_{(\rho-1)})(\mathbf{\Lambda} \otimes \mathbf{\Lambda}_{(\rho-1)})(\mathbf{U}^{\mathrm{T}} \otimes \mathbf{U}_{(\rho-1)}^{\mathrm{T}})\mathbf{P}_{(\rho)}, \quad (26) $$

where $\mathbf{U}_{(\rho-1)}$ and $\mathbf{\Lambda}_{(\rho-1)}$ are respectively the matrices of eigenvectors and eigenvalues of $\mathbf{L}^{(\rho-1)}$, and $\mathbf{P}_{(\rho)} \in \{0, 1\}^{N^2 \times N}$ is a partial permutation matrix for $\mathbf{L}^{(\rho)}$. □

## APPENDIX C

### PROOF OF THEOREM 3

*Proof.* Firstly, we express the perturbed Laplacian matrix $\hat{\mathbf{L}}$ based on the perturbation model in (15) as:

$$ \hat{\mathbf{L}} = \hat{\mathbf{D}} - \hat{\mathbf{A}} = (\mathbf{D} + \mathbf{D}_E) - (\mathbf{A} + \mathbf{E}), \quad (27) $$

where

$$ \mathbf{D}_E = \mathrm{diag}(\mathbf{E}\mathbf{1}) = \mathrm{diag}(d_{e_1}, \dots, d_{e_N}). \quad (28) $$

With the notation of

$$ \hat{\mathbf{L}}^{(\rho)} = \underbrace{(\hat{\mathbf{L}} + \epsilon\mathbf{I}) \circ \dots \circ (\hat{\mathbf{L}} + \epsilon\mathbf{I})}_{\rho \text{ times}}, \quad (29) $$

the next lemma expresses the elements of $\mathbf{L}^{(\rho)}$ and $\hat{\mathbf{L}}^{(\rho)}$, i.e., $\{(\mathbf{L}^{(\rho)})_{ij}\}_{i,j=1}^N$ and $\{(\hat{\mathbf{L}}^{(\rho)})_{ij}\}_{i,j=1}^N$, respectively, based on the edge weights and node degrees as follows:

**Lemma 1.** The elements of $\mathbf{L}^{(\rho)}$ and $\hat{\mathbf{L}}^{(\rho)}$, i.e., $\{(\mathbf{L}^{(\rho)})_{ij}\}_{i,j=1}^N$ and $\{(\hat{\mathbf{L}}^{(\rho)})_{ij}\}_{i,j=1}^N$, can be written as:

$$ \forall i \neq j : (\mathbf{L}^{(\rho)})_{ij} = (-1)^{\rho} a_{ij}^{\rho}; \quad (\hat{\mathbf{L}}^{(\rho)})_{ij} = (-1)^{\rho}(a_{ij} + e_{ij})^{\rho} \\ (\mathbf{L}^{(\rho)})_{ii} = (d_i + \epsilon)^{\rho}; \quad (\hat{\mathbf{L}}^{(\rho)})_{ii} = (d_i + d_{e_i} + \epsilon)^{\rho}. \quad (30) $$

*Proof.* Using the convention (29) and zero-diagonality of $\mathbf{A}$, one can write:

$$ \begin{aligned} \mathbf{L}^{(\rho)} &= \underbrace{(\mathbf{L} + \epsilon\mathbf{I}) \circ \dots \circ (\mathbf{L} + \epsilon\mathbf{I})}_{\rho \text{ times}} \\ &= \underbrace{(\mathbf{D} - \mathbf{A} + \epsilon\mathbf{I}) \circ \dots \circ (\mathbf{D} - \mathbf{A} + \epsilon\mathbf{I})}_{\rho \text{ times}} \\ &= \underbrace{(\mathbf{D} + \epsilon\mathbf{I})^{\rho} + (-1)^{\rho} \mathbf{A} \circ \dots \circ \mathbf{A}}_{\rho \text{ times}} \\ &\rightarrow (\mathbf{L}^{(\rho)})_{ii} = (d_i + \epsilon)^{\rho}; \quad (\mathbf{L}^{(\rho)})_{ij} = (-1)^{\rho} a_{ij}^{\rho}. \end{aligned} \quad (31) $$

Similarly for $\hat{\mathbf{L}}^{(\rho)}$ and by considering the zero-diagonality of $\mathbf{A}$ and $\mathbf{E}$, one can write:

$$ \begin{aligned} \hat{\mathbf{L}}^{(\rho)} &= \underbrace{(\hat{\mathbf{L}} + \epsilon\mathbf{I}) \circ \dots \circ (\hat{\mathbf{L}} + \epsilon\mathbf{I})}_{\rho \text{ times}} \\ &= \underbrace{(\mathbf{D} + \mathbf{D}_E + \epsilon\mathbf{I} - (\mathbf{A} + \mathbf{E})) \circ \dots \circ (\mathbf{D} + \mathbf{D}_E + \epsilon\mathbf{I} - (\mathbf{A} + \mathbf{E}))}_{\rho \text{ times}} \\ &= \underbrace{(\mathbf{D} + \mathbf{D}_E + \epsilon\mathbf{I})^{\rho} + (-1)^{\rho} (\mathbf{A} + \mathbf{E}) \circ \dots \circ (\mathbf{A} + \mathbf{E})}_{\rho \text{ times}} \\ &\rightarrow (\hat{\mathbf{L}}^{(\rho)})_{ii} = (d_i + d_{e_i} + \epsilon)^{\rho}; \quad (\hat{\mathbf{L}}^{(\rho)})_{ij} = (-1)^{\rho}(a_{ij} + e_{ij})^{\rho}. \end{aligned} \quad (32) $$
<figure><img src="image_8.png" /></figure>

Besides, the next lemma gives upper bounds for $\|\hat{\mathbf{L}}\|$ and Laplacian perturbation distance $\|\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)}\|$. Note that $r_1(M)$ and $c_1(M)$ denote the maximum Euclidean row and column norms of the matrix $M$.

**Lemma 2.** The Laplacian perturbation distance $\|\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)}\|$ and $\|\hat{\mathbf{L}}\|$ are upper bounded by:

$$ \|\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)}\| \le \sum_{m=1}^{\rho} \binom{\rho}{m} \left( \eta_m \|\mathbf{E}^{(m)}\| + (d^{max} + \epsilon)^{\rho-m} (d_e^{max})^m \right), \quad (33) $$

and

$$ \|\hat{\mathbf{L}}^{(\rho)}\| \le (\hat{d}^{max} + \epsilon)^{\rho} + \|\hat{\mathbf{A}}^{(\rho)}\|, \quad (34) $$

where $\hat{d}^{max} = d^{max} + d_e^{max}$, $d^{max} = \max_i d_i$, and $d_e^{max} = \max_i d_{e_i}$.

*Proof.* For the sake of simplicity, we start with $\|\hat{\mathbf{L}}^{(\rho)}\|$. Using (32) and by exploiting the triangular inequality principle of the norms, one can write:

$$ \begin{aligned} \|\hat{\mathbf{L}}^{(\rho)}\| &\le \|(\mathbf{D} + \mathbf{D}_E + \epsilon \mathbf{I})^{\rho} + (-1)^{\rho} \underbrace{(\mathbf{A} + \mathbf{E}) \circ \dots \circ (\mathbf{A} + \mathbf{E})}_{\rho \text{ times}}\| \\ &\le \|(\mathbf{D} + \mathbf{D}_E + \epsilon \mathbf{I})^{\rho}\| + \underbrace{\|(\mathbf{A} + \mathbf{E}) \circ \dots \circ (\mathbf{A} + \mathbf{E})\|}_{\rho \text{ times: } \hat{\mathbf{A}}^{(\rho)}} \\ &\le \underbrace{(d^{max} + d_e^{max})}_{\hat{d}^{max}} + \epsilon^{\rho} + \|\hat{\mathbf{A}}^{(\rho)}\|. \end{aligned} \quad (35) $$

Regarding $\|\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)}\|$, it is easier to work on the Laplacian elements. We have:

$$ \begin{aligned} \forall i \neq j: & (\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)})_{ij} = (-1)^{\rho+1}((a_{ij} + e_{ij})^{\rho} - a_{ij}^{\rho}) \\ &= (-1)^{\rho+1} \left( \sum_{m=1}^{\rho} \binom{\rho}{m} a_{ij}^{\rho-m} e_{ij}^{m} \right), \end{aligned} \quad (36) $$

$$ \begin{aligned} (\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)})_{ii} &= (d_i + \epsilon)^{\rho} - (d_i + \epsilon + d_{e_i})^{\rho} \\ &= -\sum_{m=1}^{\rho} \binom{\rho}{m} (d_i + \epsilon)^{\rho-m} d_{e_i}^{m}. \end{aligned} \quad (37) $$

Therefore, using (36) and (37), we can obtain a matrix-form to express $\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)}$ as:

$$ \mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)} = \sum_{m=1}^{\rho} \binom{\rho}{m} \left( (-1)^{\rho+1} \mathbf{A}^{(\rho-m)} \circ \mathbf{E}^{(m)} + (\mathbf{D} + \epsilon \mathbf{I})^{(\rho-m)} \circ \mathbf{D}_E^{(m)} \right) \quad (38) $$

Next, by applying the triangular principle of the norm, one can write:

$$ \|\mathbf{L}^{(\rho)} - \hat{\mathbf{L}}^{(\rho)}\| \le \sum_{m=1}^{\rho} \binom{\rho}{m} \left( \|\mathbf{A}^{(\rho-m)} \circ \mathbf{E}^{(m)}\| + \|(\mathbf{D} + \epsilon \mathbf{I})^{(\rho-m)} \circ \mathbf{D}_E^{(m)}\| \right). \quad (39) $$

The next Lemma makes it possible to take further steps for simplifying the current bound in (39).

**Lemma 3.** (Lemma 2c, [47]). For any symmetric matrices $\mathbf{A}_1$ and $\mathbf{A}_2$, there is an upper bound on $\|\mathbf{A}_1 \circ \mathbf{A}_2\|$ as:

$$ \|\mathbf{A}_1 \circ \mathbf{A}_2\| \le \min\{c_1(\mathbf{A}_1), r_1(\mathbf{A}_1)\} \|\mathbf{A}_2\| \quad (40) $$

Using Lemma 3, we can better simplify the bound in (39). We first apply Lemma 3 on the first term of the bound in (39):

$$ \|(\mathbf{A}^{(\rho-m)} \circ \mathbf{E}^{(m)})\| \le \underbrace{\min\{r_1(\mathbf{A}^{(\rho-m)}), c_1(\mathbf{A}^{(\rho-m)})\}}_{\eta_m} \|\mathbf{E}^{(m)}\|. \quad (41) $$

Similarly, the second term can be written as:

$$ \begin{aligned} & \|(\mathbf{D} + \epsilon \mathbf{I})^{(\rho-m)} \circ \mathbf{D}_E^{(m)}\| \\ &\le \underbrace{\min\{r_1((\mathbf{D} + \epsilon \mathbf{I})^{(\rho-m)}), c_1((\mathbf{D} + \epsilon \mathbf{I})^{(\rho-m)})\}}_{(d_e^{max}+\epsilon)^{\rho-m}} \underbrace{\|\mathbf{D}_E^{(m)}\|}_{\|\hat{\mathbf{A}}^{(\rho)}\|}. \end{aligned} \quad (42) $$

Finally, by plugging (41) and (42) into (39), the desired result in (33) is obtained, and it completes the proof of Lemma 2. □

Now, using the triangular principle of the norms, the property of being Lipschitz for the nonlinearity function $\sigma(\cdot)$, and the results in Lemma 2, we come back to find an upper bound for $d(\mathbf{Y}, \hat{\mathbf{Y}})$ in (18) as follows:

$$ \begin{aligned} & \|\sigma(\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}}) - \sigma(\mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W})\| \le \|\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}} - \mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W}\| \\ &= \|\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}} - \hat{\mathbf{L}}^{(\rho)}\mathbf{X}\mathbf{W} + \hat{\mathbf{L}}^{(\rho)}\mathbf{X}\mathbf{W} - \mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W}\| \\ &\le \|\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}} - \hat{\mathbf{L}}^{(\rho)}\mathbf{X}\mathbf{W}\| + \|\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\mathbf{W} - \mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W}\| \\ &\le \|\hat{\mathbf{L}}^{(\rho)}\| \|\mathbf{X}\| \|\mathbf{W} - \hat{\mathbf{W}}\| + \|\hat{\mathbf{L}}^{(\rho)} - \mathbf{L}^{(\rho)}\| \|\mathbf{X}\| \|\mathbf{W}\| \\ &\le \sqrt{v} \left( \left[ (\hat{d}^{max} + \epsilon)^{\rho} + \|\hat{\mathbf{A}}^{(\rho)}\| \right] \delta_{\mathbf{W}} + \sum_{m=1}^{\rho} \binom{\rho}{m} \left( \eta_m \|\mathbf{E}^{(m)}\| + (d^{max} + \epsilon)^{\rho-m} (d_e^{max})^m \right) \|\mathbf{W}\| \right). \end{aligned} \quad (43) $$

Then, using $\|\mathbf{E}\| \le \xi$ and $\|\mathbf{D}_E\| \le \xi_D$, where $\xi$ and $\xi_D$ are sufficiently small, one can use the first-order Taylor approximation (by considering only $m=1$) to further simplify the upper bound as:

$$ \begin{aligned} & \|\sigma(\hat{\mathbf{L}}^{(\rho)}\mathbf{X}\hat{\mathbf{W}}) - \sigma(\mathbf{L}^{(\rho)}\mathbf{X}\mathbf{W})\| \\ &\le \sqrt{v} \left( \left[ (\hat{d}^{max} + \epsilon)^{\rho} + \|\hat{\mathbf{A}}^{(\rho)}\| \right] \delta_{\mathbf{W}} \right. \\ &\qquad \left. + \rho \left[ \eta \|\mathbf{E}\| + (d^{max} + \epsilon)^{\rho-1} d_e^{max} \right] \|\mathbf{W}\| \right) + \mathcal{O}(\xi^2) + \mathcal{O}(\xi_D^2). \end{aligned} \quad (44) $$

and this completes the proof of Theorem 3. Note that, in the last line, we used $\|\mathbf{X}\| \le \|\mathbf{X}\|_F = \sqrt{v}$. □

# APPENDIX D

## EXPERIMENTAL STABILITY ANALYSIS

We validate Theorem 3 through experiments by analyzing how the theoretical bounds on the Right Hand Side (RHS) and Left Hand Side (LHS) depend on: 1) the Signal-to-Noise-Ratio (SNR) in both adjacency and weight matrices, 2) the Hadamard order $\rho$, 3) the graph sparsity (relevant to $p_{ER}$),
<figure><img src="image_9.png" /><figcaption>Fig. 7: Evaluation of Theorem 3 across different configurations for the Sobolev norm when certain perturbations are made to the adjacency matrix.</figcaption></figure>

and 4) the Sobolev coefficient $\epsilon$. More precisely, we build the underlying graph using the ER model with $N$ nodes and edge probability $p_{ER}$. We then generate true and perturbed elements using (16) with varying SNRs $\in \{5, 10, 20, 30, 40\}$ to perturb the underlying adjacency matrix and weight matrix $\mathbf{W}$. Both the initial feature matrix $\mathbf{X} \in \mathbb{R}^{N \times F_0}$ and weight matrix $\mathbf{W} \in \mathbb{R}^{F_0 \times F_1}$ are drawn from normal distributions. In all scenarios, we set $F_0 = 16$ and $F_1 = 2$, and the results are averaged over 100 random seeds. We analyze the following scenarios:

* Varying $\rho$: we test $\rho$ in $\{2, 3\}$ by letting $N = 10, p_{ER} = 0.3$, and $\epsilon = 0.5$.

* Varying $p_{ER}$: we test $p_{ER}$ in $\{0.1, 0.3, 0.5\}$ by letting $N = 10, \rho = 2$, and $\epsilon = 0.5$.

* Varying $\epsilon$: we test $\epsilon$ in $\{0.5, 5, 10\}$ by letting $N = 10, \rho = 2$, and $p_{ER} = 0.3$.

Figure 7 shows the averaged values of RHS and LHS across different SNRs and varying $\rho$, $p_{ER}$, and $\epsilon$. We find that higher values of $\rho$ make the network more sensitive to perturbations, especially with low SNRs. Additionally, the network benefits from graph sparsity, as it is more robust to perturbations with lower values of $p_{ER}$. We also observe that high values of $\epsilon$ can negatively affect the robustness of S2-GNN. Finally, as the SNR increases, the upper bound becomes tighter, which is related to smaller values of $\delta_W$, $\|\mathbf{E}\|$, $d_e^{max}$, and also $\zeta$ and $\zeta_D$.

## APPENDIX E
### HOMOPHILY INDEX

The homophily index for the datasets used to construct the graphs is presented in Table VII. $H(G) \to 1$ when graphs exhibit strong homophily, while $H(G) \to 0$ for graphs with strong heterophily [48]. Thus, low homophily corresponds to high heterophily, and vice versa. The learned graphs consistently exhibit stronger homophily, which directly correlates with the results in Table II, where methods applied to the learned graphs generally perform better.

## APPENDIX F
### GRAPH CLASSIFICATION EXPERIMENTS

We perform additional comparisons between the state-of-the-art methods and S2-GNN for graph classification tasks.

<table><caption>TABLE VII: Homophily index $H(G)$ for the constructed graphs with different values of $k$.</caption><thead><tr><th rowspan="2">Dataset</th><th colspan="4">k-NN Graphs</th><th colspan="4">Learned Graphs</th></tr><tr><th>k = 10</th><th>k = 20</th><th>k = 30</th><th>k = 40</th><th>k = 10</th><th>k = 20</th><th>k = 30</th><th>k = 40</th></tr></thead><tbody><tr><td>Cancer-B</td><td>0.8625</td><td>0.8447</td><td>0.8334</td><td>0.8230</td><td>0.8838</td><td>0.8670</td><td>0.8568</td><td>0.8463</td></tr><tr><td>Cancer-M</td><td>0.6267</td><td>0.5954</td><td>0.5765</td><td>0.5629</td><td>0.6534</td><td>0.6242</td><td>0.6054</td><td>0.5921</td></tr><tr><td>20News</td><td>0.5557</td><td>0.4968</td><td>0.4615</td><td>0.4354</td><td>0.6103</td><td>0.5348</td><td>0.4919</td><td>0.4628</td></tr><tr><td>HAR</td><td>0.8849</td><td>0.8478</td><td>0.8211</td><td>0.8007</td><td>0.9087</td><td>0.8851</td><td>0.8660</td><td>0.8500</td></tr><tr><td>Isolet</td><td>0.7346</td><td>0.6851</td><td>0.6454</td><td>0.6105</td><td>0.7619</td><td>0.7255</td><td>0.6943</td><td>0.6664</td></tr></tbody></table>

<table><caption>TABLE VIII: Accuracy (in %) comparison for graph classification.</caption><thead><tr><th>Model</th><th>ENZYMES</th><th>MUTAG</th><th>PROTEINS</th></tr></thead><tbody><tr><td>Cheby [2]</td><td><strong>29.39</strong><sup>±1.21</sup></td><td>76.11±2.54</td><td>67.24±0.59</td></tr><tr><td>GCN [3]</td><td>22.69±1.18</td><td>73.14±3.86</td><td><strong>68.15</strong><sup>±0.80</sup></td></tr><tr><td>GAT [4]</td><td>22.12±1.04</td><td>74.98±5.53</td><td>64.89±1.23</td></tr><tr><td>SGC [29]</td><td>19.71±0.84</td><td>66.73±3.86</td><td>61.29±0.73</td></tr><tr><td>ClusterGCN [30]</td><td>22.56±0.81</td><td><strong>81.14</strong><sup>±1.11</sup></td><td><u>67.54</u><sup>±0.57</sup></td></tr><tr><td>SuperGAT [28]</td><td>21.97±0.97</td><td>71.70±6.24</td><td>65.61±1.11</td></tr><tr><td>Transformer [21]</td><td><u>24.74</u><sup>±1.07</sup></td><td><u>80.78</u><sup>±2.33</sup></td><td>64.71±1.18</td></tr><tr><td>GATv2 [20]</td><td>20.78±1.24</td><td>71.43±6.00</td><td>64.47±1.13</td></tr><tr><td>S2-GNN (ours)</td><td>23.66±0.90</td><td>73.60±1.46</td><td>67.20±0.78</td></tr></tbody></table>

Table VIII shows the results in three datasets for graph classification: ENZYMES, MUTAG, and PROTEINS [49]. We observe that S2-GNN presents competitive results regarding previous state-of-the-art models. However, we believe that further theoretical and experimental studies should be considered in future work for graph classification.