SAGMAN: Stability Analysis of Graph Neural Networks on the Manifolds

## Abstract

Abstract Modern graph neural networks (GNNs) can be sensitive to changes in the input graph structure and node features, potentially resulting in unpredictable behavior and degraded performance. In this work, we introduce a spectral framework known as SAGMAN for examining the stability of GNNs. This framework assesses the distance distortions that arise from the nonlinear mappings of GNNs between the input and output manifolds: when two nearby nodes on the input manifold are mapped (through a GNN model) to two distant ones on the output manifold, it implies a large distance distortion and thus a poor GNN stability. We propose a distance-preserving graph dimension reduction (GDR) approach that utilizes spectral graph embedding and probabilistic graphical models (PGMs) to create low-dimensional input/output graph-based manifolds for meaningful stability analysis. Our empirical evaluations show that SAGMAN effectively assesses the stability of each node when subjected to various edge or feature perturbations, offering a scalable approach for evaluating the stability of GNNs, extending to applications within recommendation systems. Furthermore, we illustrate its utility in downstream tasks, notably in enhancing GNN stability and facilitating adversarial targeted attacks.

## 1 Introduction

The advent of Graph Neural Networks (GNNs) has sparked a significant shift in machine learning (ML), particularly in the realm of graph-structured data (Keisler, ; Hu et al., ; Kipf & Welling, ; Veličković et al., ; Zhou et al., ). By seamlessly integrating graph structure and node features, GNNs yield low-dimensional embedding vectors that maximally preserve the graph structural information (Grover & Leskovec, ). Such networks have been successfully deployed in a broad spectrum of real-world applications, including but not limited to recommendation systems
(Fan et al., ), traffic flow prediction (Yu et al., ), chip placement (Mirhoseini et al., ), and social network analysis (Ying et al., ). However, the enduring challenge in the deployment of GNNs pertains to their stability, especially when subjected to perturbations in the graph structure
(Sun et al., ; Jin et al., ; Xu et al., ). Recent studies suggest that even minor alterations to the graph structure (encompassing the addition, removal, or rearrangement of edges)
can have a pronounced impact on the performance of GNNs (Zügner et al., ; Xu et al., ). This phenomenon is particularly prominent in tasks such as node classification (Yao et al., ; Veličković et al., ; Bojchevski & Günnemann, ). The concept of stability here transcends mere resistance to adversarial attacks, encompassing the network’s ability to maintain consistent performance despite inevitable variations in the input data (graph structure and node features).

In the literature, while there are studies primarily focused on developing more stable GNN architectures (Wu et al., ; Zhao et al., ; Song et al., ; Gravina et al., ), a few attempts to analyze GNN stability. Specifically, (Keriven et al., ) first studied the stability of graph convolutional networks (GCN) on random graphs under small deformation. Later, (Gama et al., ) and (Kenlay et al., ) explored the robustness of various graph filters, which are then used to measure the stabilities of the corresponding (spectral-based) GNNs. However, these prior methods are limited to either synthetic graphs or specific GNN models.

In this work, we present SAGMAN, a novel framework devised to quantify the stability of GNNs through individual nodes. This is accomplished by assessing the resistance-distance distortions incurred by the nonlinear map (GNN model) between low-dimensional input and output graph-based manifolds. It is crucial to note that this study aims to offer significant insights for understanding and enhancing the stability of GNNs. The key technical contributions of this work are outlined below:

$\bullet$ This study introduces a spectral framework (SAGMAN) for measuring the stability of GNN models at the node level. This is achieved by effectively assessing the distance distortions caused by the maps between the input and output smooth manifolds.

$\bullet$ To construct the input smooth manifold for stability analysis in SAGMAN, we introduce a nonlinear graph dimension reduction (GDR) framework to transform the original input graph (with node features) into a low-dimensional graph-based manifold that can well preserve the original graph’s spectral properties, such as
effective resistance distances between nodes.

$\bullet$ SAGMAN has been empirically evaluated and shown to be effective in assessing the stability of individual nodes across various GNN models in realistic graph datasets. Moreover, SAGMAN allows for more powerful adversarial targeting attacks and greatly improving the stability (robustness) of the GNNs.

$\bullet$ SAGMAN has a near-linear time complexity and its data-centric nature allows it to operate across various GNN variants, independent of label information, network architecture, and learned parameters, demonstrating its wide applicability.

## 2 Background

### 2.1 Stability Analysis of ML Models on the Manifolds

The stability of a machine learning (ML) model refers to its ability to produce consistent outputs despite small variations or noise in the input data (Szegedy et al., ). To assess this stability, we utilize the *Distance Mapping Distortion* (DMD) metric (Cheng et al., ). For two input data samples $p$ and $q$, the DMD metric $\delta^{M}(p,q)$ is defined as the ratio of their distance on the output manifold to the one on the input manifold:

$$ $\delta^{M}(p,q)\overset{\mathrm{def}}{=}\frac{d_{Y}(p,q)}{d_{X}(p,q)}.$ (1) $$

By evaluating $\delta^{M}(p,q)$ for each pair of data samples, we can assess the stability of the ML model. Specifically, if two nearby data samples on the input manifold are mapped to distant points on the output manifold, this indicates a large $\delta^{M}(p,q)$ or equivalently a large local Lipschitz constant, and thus poor stability of the model near these samples; On the other hand, a small $\delta^{M}(p,q)$ implies that the model is stable in that region.

### 2.2 Previous Investigations on the Stability of GNNs

The stability of a GNN refers to its output stability in the presence of edge/node perturbations (Sun et al., ). This includes maintaining the fidelity of predictions and outcomes when subjected to changes such as edge alterations or feature attacks. A desired GNN model is expected to exhibit good stability, wherein every predicted output or the graph embeddings do not change drastically in response to the aforementioned minor perturbations (Jin et al., ; Zhu et al., ). Several recent studies have underlined the importance of analyzing the stability of GNNs. For instance,
Sharma et al. () examined the task and model-agnostic vulnerabilities of GNNs, demonstrating that these networks remain susceptible to adversarial perturbations regardless of the specific downstream tasks. Furthermore, Huang et al. () investigated robust graph representation learning via predictive coding, proposing a method that enhances GNN stability by reconstructing input data to mitigate the effects of adversarial attacks.

Moreover, while recent studies have investigated the stability issues of Graph Neural Networks (GNNs) on synthetic graphs or specific models (Keriven et al., ; Kenlay et al., ), they have not provided a unified framework for evaluating GNN stability. The latest state-of-the-art method (Li et al., ) develops such a unified framework but does not consider feature perturbations. This omission leaves a significant gap in understanding how GNNs respond to changes in node features, which is crucial for applications where feature data may be noisy or subject to perturbations.

## 3 The SAGMAN Framework for Stability Analysis of GNNs

#### Challenges in Applying DMD Metrics to GNN Stability Analysis.

When adopting the DMD metric for the stability analysis of GNN models, a natural approach is to use graph-based manifolds (details in Appendix A.1) for calculating DMDs. However, directly using the input graph structure as the input graph-based manifold may not produce satisfactory results as shown in our empirical results presented in Table 2 and Appendix A.8. This inadequacy stems from the fact that the original input graph data (including node features) may not reside near a low-dimensional manifold (Bruna et al., ), while meaningful DMD-based stability analysis requires both the input and output data samples to lie near low-dimensional manifolds (Cheng et al., ). Therefore, a naive application of DMD metrics on the original graph structure is insufficient for assessing the stability of GNNs.

### 3.1 Overview of SAGMAN

Figure: Figure 1: The proposed SAGMAN framework for stability analysis of GNNs on the manifolds.
Refer to caption: x1.png

To extend the applicability of the DMD metric to GNN settings, we introduce SAGMAN, a spectral framework for stability analysis of GNNs. A key component of SAGMAN is a novel distance-preserving Graph Dimensionality Reduction (GDR) algorithm that transforms the original input graph data—including node features and graph topology, which may reside in high-dimensional space—into a low-dimensional graph-based manifold.

As illustrated in Figure 1, the SAGMAN framework comprises three main phases:

- •
Phase 1: Creation of the input graph embedding matrix based on both node features and spectral properties of the graph. This embedding matrix is essential for the subsequent GDR step.
- •
Phase 2: Construction of low-dimensional input and output graph-based manifolds using a Probabilistic Graphical Model (PGM) approach.
- •
Phase 3: Node stability evaluation using the DMD metric, leveraging a spectral graph embedding scheme that utilizes generalized Laplacian eigenvalues and eigenvectors.

Detailed descriptions of each phase are provided in the following sections. The complete algorithmic flow of SAGMAN is shown in Section 3.5.

### 3.2 Phase 1: Embedding Matrix Construction for GDR

The calculation of the DMD metric inherently relies on pairwise distances between nodes. For graph-based manifolds, two widely used metrics for evaluating pairwise node distances are: (1) the shortest-path distance and (2) the effective resistance distance. As discussed in Appendix A.1, the effective resistance distance is more closely related to the graph’s structural (spectral) properties, providing a more meaningful measure of connectivity between nodes. Moreover, previous studies have demonstrated a significant correlation between effective resistance distances and the stability of machine learning models (Cheng et al., ). Therefore, our focus is on performing dimensionality reduction on the input graph while preserving its original effective resistance distances.

#### Graph Dimensionality Reduction via Laplacian Eigenmaps.

To achieve dimensionality reduction of the input graph, we utilize the widely recognized nonlinear dimensionality reduction algorithm, *Laplacian Eigenmaps* (Belkin & Niyogi, ). The Laplacian Eigenmaps algorithm begins by constructing an undirected graph where each node represents a high-dimensional data sample, and edges encode the similarities between data samples. It then computes the eigenvectors corresponding to the smallest eigenvalues of the graph Laplacian matrix, which are used to map each node into a low-dimensional space while preserving the local relationships between the data samples.

#### Spectral Embedding with Eigengaps.

For a graph with $N$ nodes, a straightforward approach to apply Laplacian Eigenmaps is to compute an $N\times N$ spectral embedding matrix using the complete set of graph Laplacian eigenvectors and eigenvalues (Ng et al., ), representing each node with an $N$-dimensional vector. However, computing the full set of eigenvalues and eigenvectors is computationally prohibitive for large graphs.

To address this issue, we leverage a theoretical result from spectral graph clustering (Peng et al., ), which shows that the existence of a significant gap between consecutive eigenvalues—known as an *eigengap*—implies that the graph can be well represented in a lower-dimensional space. Specifically, the eigengap is defined as $\Upsilon(k)=\frac{\lambda_{k+1}}{\rho(k)}$, where $\rho(k)$ denotes the $k$-way expansion constant, and $\lambda_{k+1}$ is the $(k+1)$-th smallest eigenvalue of the normalized Laplacian matrix. A significant eigengap indicates the existence of a $k$-way partition where each cluster has low conductance, meaning the graph is well-clustered. Based on this, we can approximate the spectral embedding using only the first $k$ eigenvalues and eigenvectors. We define the weighted spectral embedding matrix as follows:

###### Definition 3.1 .

For a connected graph $G=(V,E,w)$ with its $k$ smallest nonzero Laplacian eigenvalues denoted by $0<\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{k}$ and corresponding eigenvectors $u_{1},u_{2},\dots,u_{k}$, the weighted spectral embedding matrix is defined as
$U_{k}=\left[\frac{u_{1}}{\sqrt{\lambda_{1}}},\dots,\frac{u_{k}}{\sqrt{\lambda_
{k}}}\right]\in\mathbb{R}^{|V|\times k}$.

For a graph with a significant eigengap $\Upsilon(k)$, this embedding matrix allows us to represent each node with a $k$-dimensional vector such that the effective resistance distance between any pair of nodes can be well approximated by
$d^{\mathrm{eff}}(p,q)\approx\|U_{k}^{\top}e_{p,q}\|_{2}^{2}$,
where $e_{p}\in\mathbb{R}^{|V|}$ is the standard basis vector with a 1 at position $p$ and zeros elsewhere, and $e_{p,q}=e_{p}-e_{q}$.

###### Definition 3.1 .

#### Using Eigengaps for Graph Dimension Estimation.

Determining the precise graph dimension required to embed a graph into Euclidean space while preserving certain properties (e.g., unit edge lengths) is an NP-hard problem (Erdös et al., ; Schaefer, ). However, the presence of a significant eigengap $\Upsilon(k)$ suggests that the graph can be well represented in a $k$-dimensional space (Peng et al., ), making $k$ an approximate measure of the graph’s intrinsic dimensionality. While computing the exact value of $\Upsilon(k)$ may be challenging in practice, we can use the identified eigengap as an indicator of the suitability of SAGMAN for a given graph: graphs with significant eigengaps are more suitable to our framework since they can be effectively represented in low-dimensional spaces. Empirically, for datasets with $c$ classes, we can approximate $k$ as $k\approx 10c$ to effectively capture significant eigengaps (Deng et al., ).

### 3.3 Phase 2: Manifold Construction via PGMs

While the original Laplacian Eigenmaps algorithm suggests constructing graph-based manifolds using $k$-nearest-neighbor graphs, we find that these manifolds do not adequately preserve the original effective resistance distances. Consequently, they are not suitable for GNN stability analysis, as empirically demonstrated in Appendix A.6.

#### Probabilistic Graphical Models (PGMs) for Graph-based Manifold Learning.

PGMs, also known as Markov Random Fields (MRFs), are powerful tools in machine learning and statistical physics for representing complex systems with intricate dependency structures (Roy et al., ). PGMs encode the conditional dependencies between random variables through an undirected graph structure (see Appendix A.3 for more details). Recent studies have shown that the graph structure learned through PGMs can have resistance distances that encode the Euclidean distances between their corresponding data samples (Feng, ).

In our context, each column vector in the embedding matrix $U_{k}$ (as defined in Definition 3.1) corresponds to a data sample used for graph topology learning. By constructing the low-dimensional graph-based manifold via PGMs, we can effectively preserve the resistance distances from the original graph. Empirical evidence, detailed in Appendix A.6, supports that preserving these distances is essential for distinguishing between stable and unstable nodes. Therefore, our methodology leverages PGMs to maintain accurate effective resistance distances.

However, existing methods for learning PGMs may require numerous iterations to achieve convergence (Feng, ), limiting their applicability to large graphs.

#### Scalable PGM via Spectral Sparsification.

In the proposed SAGMAN framework, we employ PGMs to create low-dimensional input graph-based manifold $G_{X}=(V,E_{X})$ using the embedding matrix $U_{k}$ from Definition 3.1, and output manifold $G_{Y}=(V,E_{Y})$ using the GNN’s post-softmax vectors, as illustrated in Figure 1. Below, we detail the construction of the input manifold; the output manifold can be constructed similarly. Given the input embedding matrix $X=U_{k}\in\mathbb{R}^{|V|\times k}$, the maximum likelihood estimation (MLE) of the precision matrix $\Theta$ (PGM) can be obtained by solving the following convex optimization problem (see Appendix A.3 for more details) (Dong et al., ):

$$ $\max_{\Theta}\quad F(\Theta)=\log\det(\Theta)-\frac{1}{k}\operatorname{Tr}(X^{ \top}\Theta X),$ (2) $$

where $\Theta=\mathcal{L}+\frac{1}{\sigma^{2}}I$, $\operatorname{Tr}(\cdot)$ denotes the trace of a matrix, $\mathcal{L}$ is a valid Laplacian matrix, $I$ is the identity matrix, and $\sigma^{2}>0$ is a prior feature variance. To solve this, we give the following theorem:

###### Theorem 3.2 .

Maximizing the objective function in Equation 2 can be achieved in nearly-linear time via the following edge pruning strategy equivalent to spectral sparsification of the initial dense nearest-neighbor graph. Specifically, edges with small distance ratios

$$ $\rho_{p,q}=\frac{d^{\mathrm{eff}}(p,q)}{d^{\mathrm{dat}}(p,q)}=w_{p,q}\,d^{ \mathrm{eff}}(p,q)$ $$

are pruned, where $d^{\mathrm{eff}}(p,q)$ is the effective resistance distance between nodes $p$ and $q$, $d^{\mathrm{dat}}(p,q)=\|X_{p}-X_{q}\|_{2}^{2}$ is the data distance between the embeddings of nodes $p$ and $q$, and $w_{p,q}=\dfrac{1}{d^{\mathrm{dat}}(p,q)}$ is the weight of edge $(p,q)$.

The proof for Theorem 3.2 is available in Appendix A.4.

###### Theorem 3.2 .

#### Spectral Sparsification via Graph Decomposition.

Computing the edge sampling probability $\rho_{p,q}$ for each edge $(p,q)$ requires solving Laplacian matrices multiple times (Spielman & Srivastava, ), making the original sparsification method computationally expensive for large graphs. An alternative approach employs a short-cycle graph decomposition scheme (Chu et al., ), which partitions an unweighted graph $G$ into multiple disjoint cycles by removing a fixed number of edges while ensuring a bound on the length of each cycle. However, such methods are limited to unweighted graphs.

###### Lemma 3.3 .

Spectral sparsification of an undirected graph $G$, with Laplacian $L_{G}$, can be achieved by leveraging a short-cycle decomposition algorithm that returns a sparsified graph $H$, with Laplacian $L_{H}$, such that for all real vectors $x$, $x^{\top}L_{G}x\approx x^{\top}L_{H}x$ (Chu et al., ).

To extend these methods to weighted graphs, we introduce an improved spectral sparsification algorithm, illustrated in Figure 2. Our approach utilizes a *low-resistance-diameter* (LRD) decomposition scheme to limit the length of each cycle as measured by the effective resistance metric. This method is particularly effective for sparsifying weighted graphs.

Figure: Figure 2: The proposed spectral sparsification algorithm. (a) The initial graph. (b) LRD decomposition for graph clustering. (c) LSSTs for pruning non-critical edges within clusters. (d) The final graph-based manifold with two inter-cluster edges.
Refer to caption: x2.png

The key idea is to efficiently compute the effective resistance of each edge (see Appendix A.9 for details) and employ a multilevel framework to decompose the graph into several disjoint clusters bounded by an effective resistance threshold. Importantly, the inter-cluster edges identified during this process can be inserted back into the original graph to significantly enhance the stability of GNN models, as demonstrated in Section 4.5.

###### Lemma 3.3 .

### 3.4 Phase 3: Stability Analysis on the Manifolds

In this phase, we perform stability analysis of GNNs by quantifying the distortions between the input and output graph-based manifolds using the Distance Mapping Distortion (DMD) metric.

#### DMD with Effective Resistance Distance Metric.

Let $M$ be the mapping function of a machine learning model that transforms input data $X$ into output data $Y$, i.e., $Y=M(X)$. To assess the stability of $M$, we employ the DMD metric, which measures how distances between data samples are distorted by $M$. On the input and output graph-based manifolds, we use the *effective resistance distance* as the metric between nodes. The effective resistance distance between nodes $p$ and $q$ on the input manifold $G_{X}$ is computed as $d^{\mathrm{eff}}_{X}(p,q)=e_{p,q}^{\top}L_{X}^{+}e_{p,q},$ and similarly on the output manifold $G_{Y}$ as $d^{\mathrm{eff}}_{Y}(p,q)=e_{p,q}^{\top}L_{Y}^{+}e_{p,q},$ where $L_{X}^{+}$ and $L_{Y}^{+}$ denote the Moore–Penrose pseudoinverses of the Laplacian matrices $L_{X}$ and $L_{Y}$ of the input and output manifolds, respectively, and $e_{p,q}=e_{p}-e_{q}$ with $e_{p}$ being the standard basis vector corresponding to node $p$. The DMD between nodes $p$ and $q$ is defined as the ratio of the output distance to the input distance:

$$ $\delta^{M}(p,q)=\frac{d^{\mathrm{eff}}_{Y}(p,q)}{d^{\mathrm{eff}}_{X}(p,q)}= \frac{e_{p,q}^{\top}L_{Y}^{+}e_{p,q}}{e_{p,q}^{\top}L_{X}^{+}e_{p,q}}.$ (3) $$

To quantify the maximum distortion introduced by $M$, we consider the maximum DMD over all pairs of distinct nodes: $\delta_{\mathrm{max}}^{M}$.
According to Lemma A.2 in (Cheng et al., ), the optimal Lipschitz constant $K^{*}$ of the mapping $M$ is bounded by the largest eigenvalue of $L_{Y}^{+}L_{X}$ and the maximum DMD:

$$ $\delta_{\mathrm{max}}^{M}\leq K^{*}\leq\lambda_{\max}(L_{Y}^{+}L_{X}).$ (4) $$

This relationship allows us to assess the stability of $M$ using spectral properties of the Laplacian matrices.

#### Node Stability Score via Spectral Embedding.

Building on the theoretical insights from (Cheng et al., ), we utilize the largest generalized eigenvalues and their corresponding eigenvectors of $L_{Y}^{+}L_{X}$ to evaluate the stability of individual nodes in GNNs. We compute the weighted eigensubspace matrix $V_{s}\in\mathbb{R}^{|V|\times s}$ for spectral embedding of the input manifold $G_{X}=(V,E_{X})$, where $|V|$ is the number of nodes. The matrix $V_{s}$ is defined as: $V_{s}=\left[v_{1}\sqrt{\zeta_{1}},\,v_{2}\sqrt{\zeta_{2}},\,\dots,\,v_{s}\sqrt
{\zeta_{s}}\right]$,
where $\zeta_{1}\geq\zeta_{2}\geq\dots\geq\zeta_{s}$ are the largest $s$ eigenvalues of $L_{Y}^{+}L_{X}$, and $v_{1},v_{2},\dots,v_{s}$ are the corresponding eigenvectors. Using $V_{s}$, we embed the nodes of $G_{X}$ into an $s$-dimensional space by representing each node $p$ with the $p$-th row of $V_{s}$. The stability of an edge $(p,q)\in E_{X}$ can then be estimated by computing the spectral embedding distance between nodes $p$ and $q$: $|V_{s}^{\top}e_{p,q}\|_{2}^{2}$. To assess the stability of individual nodes, we define the stability score of node $p$ as the average embedding distance to its neighbors in the input manifold:

$$ $\text{score}(p)=\frac{1}{|\mathbb{N}_{X}(p)|}\sum_{q\in\mathbb{N}_{X}(p)}\|V_{ s}^{\top}e_{p,q}\|_{2}^{2},$ (5) $$

where $\mathbb{N}_{X}(p)$ denotes the set of neighbors of node $p$ in $G_{X}$. Since $\|V_{s}^{\top}e_{p,q}\|_{2}^{2}$ is proportional to $\left(\delta^{M}(p,q)\right)^{3}$, the node stability score effectively serves as a surrogate for the local Lipschitz constant, analogous to $\|\nabla_{X}M(p)\|$ under the manifold setting (Cheng et al., ). For a more detailed derivation and theoretical justification, please refer to Appendix A.11.

### 3.5 Algorithm Flow of SAGMAN

The Algorithm 1 shows the key steps in SAGMAN.

Figure: Algorithm 1 SAGMAN Algorithm

### 3.6 Time Complexity of SAGMAN

The proposed SAGMAN framework is designed to be efficient and scalable for large graphs. Below, we analyze the time complexity of its key components.
We utilize fast multilevel eigensolvers to compute the first $c$ Laplacian eigenvectors. These eigensolvers operate in nearly linear time, $O(c|V|)$, without loss of accuracy (Zhao et al., ), where $|V|$ denotes the number of nodes in the graph.
To construct the initial graph for manifold learning, we employ the $k$-nearest neighbor algorithm, which has a nearly linear computational complexity of $O(|V|\log|V|)$ (Malkov & Yashunin, ).
The spectral sparsification step, performed via Low-Resistance-Diameter (LRD) decomposition, has a time complexity of $O(|V|dm)$, where $d$ is the average degree of the graph, and $m$ is the order of the Krylov subspace used in the computation.
By leveraging fast generalized eigensolvers (Koutis et al., ; Cucuringu et al., ), we can compute all DMD values in $O(|E|)$ time, where $|E|$ denotes the number of edges in the graph.

Overall, the dominant terms in the time complexity are nearly linear with respect to the size of the graph. This near-linear scalability allows SAGMAN to handle large graph datasets efficiently. Results for runtime scalability are available in Appendix A.13

## 4 Experiments

We validate our proposed SAGMAN framework through comprehensive experiments. First, we compare exact resistance distances from the original graph with their approximations derived from the constructed manifold to demonstrate the manifold’s fidelity. Next, we present numerical experiments showcasing the effectiveness of our metric in quantifying GNN stability. We then highlight the efficacy of our SAGMAN-guided approach in executing graph adversarial attacks. Finally, we demonstrate how leveraging the low-dimensional input manifold created by SAGMAN significantly enhances GNN stability. Details of our experimental setup are provided in Appendix A.5.

### 4.1 Evaluation of Graph Dimension Reduction (GDR)

A natural concern arises regarding how much the constructed input graph-based manifold might change the original graph’s structure. To this end, we compare the exact resistance distances calculated using the complete set of eigenpairs with approximate resistance distances estimated using the embedding matrix $U_{k}$ (as defined in Definition 3.1), considering various values of $k$ as shown in Table 1. Our empirical results demonstrate that using a small number of eigenpairs can effectively approximate the original effective-resistance distances. Furthermore, Table 2 shows that the SAGMAN-guided stability analysis with GDR consistently distinguishes between stable and unstable nodes, whereas the analysis without GDR fails to do so. Additional related results for various datasets and GNN architectures are provided in Appendix A.8.

**Table 1: Resistance-distance preservation for the Cora graph, evaluating $100$ randomly selected node pairs. Larger correlation coefficients (CC) indicate more accurate estimations.**
| $k$ | 20 | 30 | 50 | 100 | 200 | 400 | 500 |
| --- | --- | --- | --- | --- | --- | --- | --- |
| CC | 0.69 | 0.78 | 0.82 | 0.87 | 0.93 | 0.97 | 0.99 |

**Table 2: Cosine similarities between original and perturbed node embeddings for stable/unstable nodes under Nettack adversarial attacks on the Cora dataset using GCN. Higher cosine similarities for stable nodes and lower for unstable nodes indicate better distinction. Better results are in bold.**
| Nettack Level | 1 | 2 | 3 |
| --- | --- | --- | --- |
| w/o GDR | 0.90/0.96 | 0.84/0.93 | 0.81/0.91 |
| w/ GDR | 0.99/0.90 | 0.98/0.84 | 0.97/0.81 |

**Table 3: Comparison of Nettack/FGA error rates for 40 nodes selected using Nettack’s recommendation, confidence ranking, and SAGMAN. All nodes chosen by SAGMAN-guided methods are correctly classified before perturbation. Better results are highlighted in bold.**
| Selection | Nettack’s default | Confidence Ranking | SAGMAN |
| --- | --- | --- | --- |
| Cora | $0.725/0.850$ | $0.925/0.775$ | 0.975/0.975 |
| Cora-ml | $0.750/0.850$ | $0.800/0.700$ | 0.950/0.950 |
| Citeseer | $0.800/0.875$ | $0.925/0.950$ | 1.000/0.975 |
| Pubmed | $0.750/0.875$ | $0.750/0.925$ | 0.825/0.950 |

**Table 4: Error rates for the Cora and Citeseer datasets before/after Nettack and IG-attack evasion attacks, with and without SAGMAN enhancement.‘Size’ indicates the fraction of the nodes selected. Better results are highlighted in bold.**
| Nettack Attack Budget Size | 1% | 5% | 10% | 15% | 20% |  |
| --- | --- | --- | --- | --- | --- | --- |
| Cora | w/o SAGMAN | 0.00/0.79 | 0.01/0.86 | 0.01/0.90 | 0.02/0.93 | 0.02/0.92 |
| w/ SAGMAN | 0.00/0.00 | 0.00/0.00 | 0.01/0.01 | 0.06/0.06 | 0.08/0.08 |  |
| Citeseer | w/o SAGMAN | 0.00/0.76 | 0.03/0.82 | 0.04/0.85 | 0.04/0.86 | 0.03/0.85 |
| w/ SAGMAN | 0.09/0.42 | 0.10/0.37 | 0.12/0.37 | 0.11/0.35 | 0.08/0.28 |  |
| IG-attack Budget Size | 1% | 5% | 10% | 15% | 20% |  |
| Citeseer | w/o SAGMAN | 0.00/0.76 | 0.02/0.74 | 0.02/0.74 | 0.02/0.72 | 0.02/0.75 |
| w/ SAGMAN | 0.00/0.05 | 0.02/0.07 | 0.02/0.07 | 0.02/0.04 | 0.02/0.09 |  |
| Cora | w/o SAGMAN | 0.00/0.96 | 0.02/0.83 | 0.02/0.81 | 0.02/0.82 | 0.02/0.81 |
| w/ SAGMAN | 0.00/0.00 | 0.02/0.02 | 0.02/0.03 | 0.02/0.02 | 0.02/0.04 |  |

**Table 5: Baseline and Perturbed error rates for the Citeseer dataset using GOOD-AT, SAGMAN, and both. Our experiments’ attack method and backbone architecture follow those from Li et al. (). The best results are marked in bold, and the second-best results are underlined.**
| Edge Attack Budget | 25 | 109 | 219 | 302 | 410 | 550 |
| --- | --- | --- | --- | --- | --- | --- |
| GCN(baseline) | 0.2968 | 0.3211 | 0.3590 | 0.3780 | 0.4064 | 0.4301 |
| GOOD-AT | 0.2808 | 0.2820 | 0.2850 | 0.2838 | 0.2855 | 0.2915 |
| SAGMAN | 0.2808 | 0.2802 | 0.2808 | 0.2820 | 0.2814 | 0.2838 |
| GOOD-AT+SAGMAN | 0.2684 | 0.2695 | 0.2684 | 0.2690 | 0.2690 | 0.2701 |

### 4.2 Metrics for GNN Stability Evaluation

Figure: Figure 3: The horizontal axes, denoted by $X$, represent the perturbation applied. ‘Random EdgesPT’ refers to the DICE attack. ‘Random FeaturePT’ indicates the application of Gaussian noise perturbation, expressed as $X\eta$ perturbation, where $\eta$ represents Gaussian noise. SAGMAN Stable/Unstable denotes the samples classified as stable or unstable by SAGMAN, respectively.
Refer to caption: x3.png

We empirically demonstrate SAGMAN’s ability to distinguish between stable and unstable samples under DICE attack (Waniek et al., ) and Gaussian perturbation, as shown in Figure 3. Additional results on various datasets—including large-scale datasets—and different GNN architectures under Nettack attacks (Zügner et al., ) are provided in Appendix A.10, further illustrating the effectiveness of our metric.

### 4.3 Stability of GNN-based Recommendation Systems

We evaluate the stability of GNN-based recommendation systems using the PinSage framework (Ying et al., ) on the MovieLens 1M dataset (Harper & Konstan, ). To construct the input graph, we first homogenize the various node and edge types into a unified format. We then calculate the weighted spectral embedding matrix $U_{k}$ (as defined in Definition 3.1) and extract user-type samples to build a low-dimensional input (user) graph-based manifold. For the output graph-based manifold, we utilize PinSage’s final output, which includes the top-10 recommended items for each user. We construct this output (user) graph-based manifold based on Jaccard similarity measures of these recommendations. Table 6 demonstrates the effectiveness of SAGMAN in distinguishing between stable and unstable users.

**Table 6: Comparison of the mean Jaccard similarity (MJS) between SAGMAN-selected stable/unstable users at a perturbation level ($l$) ranging from $1$ to $5$ where each selected user is connected to $l$ randomly chosen new items. The MJS is computed over $20$ iterations for each perturbation.**
| Perturbation Level | 1 | 2 | 3 | 4 | 5 |
| --- | --- | --- | --- | --- | --- |
| Stable Users | 0.8513 | 0.8837 | 0.8295 | 0.8480 | 0.8473 |
| Unstable Users | 0.7885 | 0.7955 | 0.8167 | 0.8278 | 0.7794 |

### 4.4 SAGMAN-guided Adversarial Targeted Attack

We adapt our GNN training methodology from previous work (Jin et al., ). Using GCN as our base model for the Citeseer, Cora, Cora-ML, and Pubmed datasets, we employ Nettack and FGA (Chen et al., ) as benchmark attack methods. For target node selection, we compare Nettack’s recommendation (Zügner et al., ), SAGMAN-guided strategies, and a heuristic confidence ranking (Chang et al., ).

Table 3 presents the error rates after applying Nettack and FGA attacks. The results demonstrate that SAGMAN-guided attacks outperform both Nettack’s recommendation and the confidence ranking, leading to more effective adversarial attacks.

### 4.5 SAGMAN-guided GNN Stability Enhancement

A naive approach to improving the stability of a GNN model is to replace the entire input graph with the low-dimensional graph-based manifold for GNN predictions. However, due to the significantly increased densities in the graph-based manifold, the GNN prediction accuracy may be adversely affected. To achieve a flexible trade-off between model stability and prediction accuracy, we select only the inter-cluster edges from the input graph-based manifold, as shown in Figure 2(d), and insert them into the original graph. Since the inter-cluster (bridge-like) edges are typically spectrally critical—with high sampling probabilities as defined in Theorem 3.2—adding them to the original graph significantly alters its structural (spectral) properties.

Whereas previous work (Deng et al., ) focuses on improving the robustness of GNNs against poisoning attacks, our work centers on enhancing robustness against evasion attacks. In Table 4, we present the error rates for both the original and the enhanced graphs, where Nettack was employed to perturb SAGMAN-selected most unstable samples within each graph. Table 5 shows the error rates for the state-of-the-art GOOD-AT (Li et al., ) method and our proposed SAGMAN approach. Since SAGMAN is a versatile plug-in method that can be combined with other robustness techniques, we also report results for the combined application of GOOD-AT and SAGMAN. Notably, reintegrating selected edges into the original graph significantly reduces the error rate when subjected to adversarial evasion attacks.

## 5 Conclusions

In this work, we introduced SAGMAN, a novel framework for analyzing the stability of GNNs at the individual node level by assessing the resistance-distance distortions between low-dimensional input and output graph-based manifolds. A key component of SAGMAN is the proposed Graph Dimensionality Reduction (GDR) approach for constructing resistance-preserving manifolds, which enables effective stability analysis.

Our experimental results demonstrate that SAGMAN effectively quantifies GNN stability, leading to significantly enhanced targeted adversarial attacks and improved GNN robustness. The current SAGMAN framework is particularly effective for graphs that can be well represented in low-dimensional spaces and exhibit large eigengaps.

## Appendix A Appendix

### A.1 Graph-based Manifolds and Distance Metrics

#### Graph-based Manifolds.

A manifold is a topological space that locally resembles Euclidean space near each point (Lee, ). In our work, we utilize *graph-based manifolds* to represent complex data structures. Specifically, we represent each manifold as an undirected (connected) graph $G=(V,E)$, where $V$ is the set of vertices corresponding to data points, and $E$ is the set of edges encoding relationships (e.g., conditional dependencies) between these points (Tenenbaum et al., ). This representation is particularly effective when the underlying data structure can be approximated by a network of discrete points, as is common in spectral clustering and manifold learning (Belkin & Niyogi, ).

#### Mappings Between Manifolds.

Understanding the *mappings between manifolds* is essential for transferring and comparing information across different data representations. In our context, these mappings correspond to transformations between the input and output graph-based manifolds of an ML model. Formally, a mapping $\varphi:M_{X}\rightarrow M_{Y}$ from the input manifold $M_{X}$ to the output manifold $M_{Y}$ allows us to analyze how the model transforms data points and their relationships (Pan & Yang, ). This is crucial for assessing the model’s stability with respect to small perturbations in the input.

#### Distance Metrics on Manifolds.

Measuring dissimilarities between points on manifolds requires appropriate distance metrics. While *geodesic distances*, representing the shortest paths between points on a manifold, are commonly used, they can be computationally intensive for large datasets and may not capture global structural information due to their inherently local nature (Tenenbaum et al., ).

#### Effective Resistance Distances on Graph-based Manifolds.

To overcome these limitations, we employ the *effective resistance distance* from electrical network theory. By modeling the undirected graph as a resistive network, the effective resistance between two nodes captures both local relationships and global structural properties of the graph (see Appendix A.1 for more details). This metric effectively combines the advantages of geodesic and global distances, providing a more comprehensive measure of similarity in graph-based manifolds (Tenenbaum et al., ). Additionally, it is computationally efficient for large-scale graphs (Klein & Zhu, ), making it practical for our analysis.

To illustrate the difference between geodesic distance and effective resistance distance, consider the example in Figure 4. In all three graphs, nodes $A$ and $B$ have the same geodesic distance (i.e., the shortest path length is identical). However, their effective resistance distances $\Omega_{AB}$ vary significantly due to differences in the global structure of each graph. This example demonstrates that while geodesic distance only accounts for the shortest path, effective resistance distance incorporates the overall connectivity and network topology, making it a more informative metric for assessing similarities in graph-based manifolds.

Figure: Figure 4: Examples of three different graph structures with nodes $A$ and $B$. Despite having the same geodesic distance (shortest path length), the effective resistance distances $\Omega_{AB}$ vary due to the different global structures of the graphs.
Refer to caption: x4.png

By incorporating effective resistance distances into our stability analysis, we can more accurately assess the distance distortions introduced by the mappings between the input and output manifolds (as discussed in Section 2.1). This provides a solid foundation for evaluating the stability of ML models on graph-based manifolds, enhancing both the theoretical rigor and practical applicability of our approach.

### A.2 Spectral Graph Theory

Spectral graph theory is a branch of mathematics that studies the properties of graphs through the eigenvalues and eigenvectors of matrices associated with the graph (Chung, ).
Let $G=(V,E,w)$ denote an undirected graph $G$, $V$ denote a set of nodes (vertices), $E$ denote a set of edges and $w$ denote the corresponding edge weights. The adjacency matrix can be defined as:

$$ ${A}(i,j)=\begin{cases}w(i,j)&\text{ if }(i,j)\in E\\ 0&\text{otherwise }\end{cases}\vspace{-0pt}$ (6) $$

The Laplacian matrix of $G$ can be constructed by ${L={D}-{A}}$, where $D$ denotes the degree matrix.

###### Lemma A.1 .

(Courant-Fischer Minimax Theorem) The $k$-th largest eigenvalue of the Laplacian matrix $L\in\mathbb{R}^{|V|\times|V|}$ can be computed as follows:

$$ $\lambda_{k}(L)=\min_{dim(U)=k}\,{\max_{\begin{subarray}{c}u_{k}\in U\\ u_{k}\neq 0\end{subarray}}{\frac{u_{k}^{\top}Lu_{k}}{u_{k}^{\top}u_{k}}}}$ (7) $$

Lemma A.1 is the Courant-Fischer Minimax Theorem (Golub & Van Loan, ) for solving the eigenvalue problem: $Lu_{k}=\lambda_{k}u_{k}$.
The generalized Courant-Fischer Minimax Theorem for solving generalized eigenvalue problem $L_{X}v_{k}=\lambda_{k}L_{Y}v_{k}$ can
be expressed as follows:

###### Lemma A.2 .

(The Generalized Courant-Fischer Minimax Theorem) Given two Laplacian matrices $L_{X},L_{Y}\in\mathbb{R}^{|V|\times|V|}$ such that $null(L_{Y})\subseteq null(L_{X})$, $L_{Y}^{+}$ denotes the Moore–Penrose pseudoinverse of $L_{Y}$, the $k$-th largest eigenvalue of $L_{Y}^{+}L_{X}$ can be computed under the condition of $1\leq k\leq rank(L_{Y})$ by:

$$ $\lambda_{k}(L_{Y}^{+}L_{X})=\min_{\begin{subarray}{c}dim(U)=k\\ U\bot null(L_{Y})\end{subarray}}\,{\max_{v_{k}\in U}{\frac{v_{k}^{\top}L_{X}v_ {k}}{v_{k}^{\top}L_{Y}v_{k}}}}.$ (8) $$

###### Lemma A.1 .

###### Lemma A.2 .

### A.3 Graph-based Manifold Learning via PGMs

Given $M$ samples of $N$-dimensional vectors stored in a data matrix $X\in{\mathbb{R}^{N\times M}}$, the recent graph topology learning methods (Kalofolias & Perraudin, ; Dong et al., ) estimate graph Laplacians from $X$ for achieving the following desired characteristics:

Smoothness of Graph Signals. The graph signals corresponding to the real-world data should be sufficiently smooth on the learned graph structure: the signal values will only change gradually across connected neighboring nodes.
The smoothness of a signal $x$ over an undirected graph $G=(V,E,w)$ can be measured with the following Laplacian quadratic form:
${x^{\top}}Lx=\sum\limits_{\left({p,q}\right)\in E}{{w_{p,q}}{{\left({x\left(p
\right)-x\left(q\right)}\right)}^{2}}}$,
where $w_{p,q}$ denotes the weight of edge ($p,q$), $L=D-W$ denotes the Laplacian, $D$ denotes the diagonal (degree) matrix, and $W$ denotes the adjacency matrix of $G$, respectively.
The smaller quadratic form implies the smoother signals across the edges in the graph.
The smoothness ($Q$) of a set of signals $X$ over graph $G$ is computed using the following matrix trace (Dong et al., ):
$Q(X,L)=Tr({X^{\top}LX})$,
where $Tr(\bullet)$ denotes the matrix trace.

Sparsity of the Estimated Graph. Graph sparsity is another critical consideration in graph learning. One of the most important motivations of learning a graph is to use it for downstream computing tasks. Therefore, more desired graph topology learning algorithms should allow better capturing and understanding the global structure (manifold) of the data set, while producing sufficiently sparse graphs that can be easily stored and efficiently manipulated in the downstream algorithms, such as circuit simulations, network partitioning, dimensionality reduction, data visualization, etc.

Problem Formulation. Consider a random vector $x\sim N(0,\Sigma)$ with probability density function:

$$ $f(x)=\frac{\exp{\left(-\frac{1}{2}x^{\top}\Sigma^{-1}x\right)}}{(2\pi)^{N/2} \det(\Sigma)^{(1/2)}}\propto\det(\Theta)^{1/2}\exp{\left(-\frac{1}{2}x^{\top} \Theta x\right)},$ (9) $$

where $\Sigma=\mathbb{E}[xx^{\top}]\succ 0$ denotes the covariance matrix, and $\Theta=\Sigma^{-1}$ denotes the precision matrix (inverse covariance matrix). Prior graph topology learning methods aim at estimating sparse precision matrix $\Theta$ from potentially high-dimensional input data, which fall into the following two categories:

(A) The graphical Lasso method aims at estimating a sparse precision matrix $\Theta$ using convex optimization to maximize the log-likelihood of $f(x)$ (Friedman et al., ):

$$ $\max_{\Theta}:\log\det(\Theta)-Tr(\Theta S)-{\beta}{{\|\Theta\|}}_{1},$ (10) $$

where $\Theta$ denotes a non-negative definite precision matrix, $S$ denotes a sample covariance matrix, and $\beta$ denotes a regularization parameter.
The first two terms together can be interpreted as the log-likelihood under a Gaussian Markov Random Field (GMRF). $\|\bullet\|$ denotes the entry-wise $\ell_{1}$ norm, so ${\beta}{{\|\Theta\|}}_{1}$ becomes the sparsity promoting regularization term. This model learns the graph structure by maximizing the penalized log-likelihood.
When the sample covariance matrix $S$ is obtained from $M$ i.i.d (independent and identically distributed) samples $X=[x_{1},...,x_{M}]$ where $X\sim N(0,S)$ has an $N$-dimensional Gaussian distribution with zero mean, each element in the precision matrix $\Theta_{i,j}$ encodes the conditional dependence between variables $X_{i}$ and $X_{j}$. For example, $\Theta_{i,j}=0$ implies that variables $X_{i}$ and $X_{j}$
are conditionally independent, given the rest.

(B) The Laplacian estimation methods have been recently introduced for more efficiently solving the following convex problem (Dong et al., ; Lake & Tenenbaum, ):

$$ ${\max_{\Theta}}:F(\Theta)=\log\det(\Theta)-\frac{1}{M}Tr({X^{\top}\Theta X})- \beta{{\|\Theta\|}}_{1},$ (11) $$

where $\Theta={L}+\frac{1}{\sigma^{2}}I$, ${L}$ denotes the set of valid graph Laplacian matrices, $I$ denotes the identity matrix, and $\sigma^{2}>0$ denotes prior feature variance. It can be shown that the three terms in (11) are corresponding to $\log\det(\Theta)$, $Tr(\Theta S)$ and ${\beta}{{\|\Theta\|}}_{1}$ in (10), respectively. Note that the second term also promotes graph sparsity, so the ${\beta}{{\|\Theta\|}}_{1}$ can be dropped without impacting the final solution. Since $\Theta={L}+\frac{1}{\sigma^{2}}I$ correspond to symmetric and positive definite (PSD) matrices (or M matrices) with non-positive off-diagonal entries, this formulation will lead to the estimation of attractive GMRFs (Dong et al., ; Slawski & Hein, ).
In case $X$ is non-Gaussian, formulation (11) can be understood as
Laplacian estimation based on minimizing the Bregman divergence between positive definite matrices induced by the function $\Theta\mapsto-\log\det(\Theta)$ (Slawski & Hein, ).

### A.4 Proof for Theorem 3.2

In the SAGMAN framework, we aim to construct low-dimensional manifolds for both the input and output of a GNN. For the input manifold, we use the embedding matrix $X=U_{k}\in\mathbb{R}^{|V|\times k}$ (as defined in Definition 3.1), where $|V|$ is the number of nodes in the graph and $k$ is the dimensionality of the embedding space. The goal is to learn a precision matrix $\Theta$ that captures the underlying graph structure reflected in the embeddings.

#### Maximum Likelihood Estimation (MLE) of the PGM (Precision Matrix).

We start by formulating the MLE of the precision matrix $\Theta$ as a convex optimization problem (Dong et al., ):

$$ $\max_{\Theta}\quad F(\Theta)=\log\det(\Theta)-\frac{1}{k}\operatorname{Tr}(X^{ \top}\Theta X),$ (12) $$

where:

- •
$\Theta=\mathcal{L}+\frac{1}{\sigma^{2}}I$,
- •
$\mathcal{L}$ is the graph Laplacian matrix,
- •
$I$ is the identity matrix,
- •
$\sigma^{2}>0$ is a prior variance term,
- •
$\operatorname{Tr}(\cdot)$ denotes the trace of a matrix.

#### Expanding the Laplacian Matrix.

The graph Laplacian matrix $\mathcal{L}$ can be decomposed as:

$$ $\mathcal{L}=\sum_{(p,q)\in E}w_{p,q}e_{p,q}e_{p,q}^{\top},$ (13) $$

where:

- •
$E$ is the set of edges in the graph,
- •
$w_{p,q}$ is the weight of edge $(p,q)$,
- •
$e_{p,q}=e_{p}-e_{q}$, with $e_{p}$ being the standard basis vector corresponding to node $p$.

The edge weights are defined as:

$$ $w_{p,q}=\frac{1}{\|X^{\top}e_{p,q}\|_{2}^{2}}=\frac{1}{\|X_{p}-X_{q}\|_{2}^{2}},$ (14) $$

where $X_{p}$ and $X_{q}$ are the embeddings of nodes $p$ and $q$, respectively.

#### Expressing the Objective Function.

Substituting $\Theta$ and $\mathcal{L}$ into the objective function, we can split $F(\Theta)$ into two parts:

$$ $F=F_{1}-\frac{1}{k}F_{2},$ (15) $$

where:

- 1.
First Term ($F_{1}$):
$F_{1}=\log\det(\Theta)=\sum_{i=1}^{|V|}\log\left(\lambda_{i}+\frac{1}{\sigma^{
2}}\right),$
(16)
with $\lambda_{i}$ being the $i$-th eigenvalue of the Laplacian $\mathcal{L}$.
- 2.
Second Term ($F_{2}$):
$F_{2}=\operatorname{Tr}(X^{\top}\Theta X)=\frac{\operatorname{Tr}(X^{\top}X)}{
\sigma^{2}}+\sum_{(p,q)\in E}w_{p,q}\|X^{\top}e_{p,q}\|_{2}^{2}.$
(17)
The term $\frac{\operatorname{Tr}(X^{\top}X)}{\sigma^{2}}$ is constant with respect to $w_{p,q}$ and can be ignored for optimization over $w_{p,q}$.

#### Computing Partial Derivatives.

To optimize $F$ with respect to the edge weights $w_{p,q}$, we compute the partial derivatives of $F_{1}$ and $F_{2}$:

#### Derivative of F 1 subscript 𝐹 1 F_{1} italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT with respect to edge weight:

$$ $\frac{\partial F_{1}}{\partial w_{p,q}}=\sum_{i=1}^{|V|}\frac{1}{\lambda_{i}+ \frac{1}{\sigma^{2}}}\frac{\partial\lambda_{i}}{\partial w_{p,q}}.$ (18) $$

Since $\frac{\partial\lambda_{i}}{\partial w_{p,q}}=v_{i}^{\top}\frac{\partial
\mathcal{L}}{\partial w_{p,q}}v_{i}$, where $v_{i}$ is the eigenvector corresponding to $\lambda_{i}$, and $\frac{\partial\mathcal{L}}{\partial w_{p,q}}=e_{p,q}e_{p,q}^{\top}$, we have:

$$ $\frac{\partial F_{1}}{\partial w_{p,q}}=\sum_{i=1}^{|V|}\frac{(v_{i}^{\top}e_{ p,q})^{2}}{\lambda_{i}+\frac{1}{\sigma^{2}}}.$ (19) $$

When $\sigma$ approaches infinity, this expression is known as the effective resistance distance between nodes $p$ and $q$:

$$ $d^{\mathrm{eff}}(p,q)=e_{p,q}^{\top}\Theta^{-1}e_{p,q}.$ (20) $$

So,

$$ $\frac{\partial F_{1}}{\partial w_{p,q}}=d^{\mathrm{eff}}(p,q).$ (21) $$

#### Derivative of F 2 subscript 𝐹 2 F_{2} italic_F start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT with respect to edge weight:

$$ $\frac{\partial F_{2}}{\partial w_{p,q}}=\|X^{\top}e_{p,q}\|_{2}^{2}=\|X_{p}-X_ {q}\|_{2}^{2}.$ (22) $$

This is the data distance between nodes $p$ and $q$:

$$ $d^{\mathrm{dat}}(p,q)=\|X_{p}-X_{q}\|_{2}^{2}=\frac{1}{w_{p,q}}.$ (23) $$

Since $w_{p,q}=\frac{1}{d^{\mathrm{dat}}(p,q)}$, we have:

$$ $\frac{\partial F_{2}}{\partial w_{p,q}}=\frac{1}{w_{p,q}}.$ (24) $$

#### Optimization Strategy.

The gradient of $F$ with respect to $w_{p,q}$ is:

$$ $\frac{\partial F}{\partial w_{p,q}}=\frac{\partial F_{1}}{\partial w_{p,q}}- \frac{1}{k}\frac{\partial F_{2}}{\partial w_{p,q}}=d^{\mathrm{eff}}(p,q)-\frac {1}{k}\frac{1}{w_{p,q}}.$ (25) $$

To maximize $F$, we can:

- •
Increase $w_{p,q}$ when $d^{\mathrm{eff}}(p,q)>\dfrac{1}{k}\dfrac{1}{w_{p,q}}$.
- •
Decrease $w_{p,q}$ when $d^{\mathrm{eff}}(p,q)<\dfrac{1}{k}\dfrac{1}{w_{p,q}}$.

However, since $w_{p,q}\geq 0$, decreasing $w_{p,q}$ effectively means pruning the edge $(p,q)$.

#### Edge Pruning Strategy.

We aim to prune edges where:

- •
Effective Resistance Distance is Small ($d^{\mathrm{eff}}(p,q)$ is small): The nodes are well-connected in terms of the graph structure.
- •
Data Distance is Large ($d^{\mathrm{dat}}(p,q)$ is large): The embeddings of the nodes are far apart.

This leads us to consider the distance ratio $\rho_{p,q}$:

$$ $\rho_{p,q}=\frac{d^{\mathrm{eff}}(p,q)}{d^{\mathrm{dat}}(p,q)}=w_{p,q}\,d^{ \mathrm{eff}}(p,q).$ (26) $$

Edges with small $\rho_{p,q}$ are candidates for pruning because they contribute less to maximizing $F$. By pruning such edges, we focus on retaining edges that have a significant impact on the graph’s spectral properties.

#### Connection to Spectral Graph Sparsification.

The ratio $\rho_{p,q}$ corresponds to the edge sampling probability used in spectral graph sparsification (Spielman & Teng, ). Spectral sparsification aims to approximate the original graph with a sparser graph that preserves its spectral (Laplacian) properties.

In spectral sparsification:

- •
Edges are sampled with probability proportional to $w_{p,q}d^{\mathrm{eff}}(p,q)$.
- •
Edges with higher $\rho_{p,q}$ are more likely to be included in the sparsified graph.

Therefore, our edge pruning strategy—performing spectral sparsification on the initial graph—is equivalent to maximize the objective function in Equation 2. This ensures that the essential structural properties of the graph are maintained while reducing complexity.

### A.5 Experimental Setup

See Appendix A.10 for detailed descriptions of all graph datasets used in this work. We employ the most popular backbone GNN models including GCN (Kipf & Welling, ), GPRGNN (Chien et al., ), GAT (Veličković et al., ), APPNP (Gasteiger et al., ), ChebNet (Defferrard et al., ), and Polynormer (Deng et al., ). The recommendation system is based on PinSage (Ying et al., ). Perturbations include Gaussian noise evasion attacks and adversarial attacks (DICE (Waniek et al., ), Nettack (Zügner et al., ), and FGA (Chen et al., )). The input graph-based manifolds are constructed using graph adjacency and node features, while the output graph-based manifold is created using post-softmax vectors. To showcase SAGMAN’s effectiveness in differentiating stable from unstable nodes, we apply SAGMAN to select $1\%$ of the entire dataset as stable nodes, and another $1\%$ as unstable nodes.
This decision stems from that only a portion of the dataset significantly impacts model stability (Cheng et al., ; Hua et al., ; Chang et al., ). Additional evaluation results on the entire dataset can be found in Appendix A.12. We quantify output perturbations using cosine similarity and Kullback-Leibler divergence (KLD). Additional insights on cosine similarity, KLD, and accuracy can be found in Appendix A.7. For a large-scale dataset “ogbn-arxiv", due to its higher output dimensionality, we focus exclusively on accuracy comparisons. This decision is informed by the KLD estimator’s $n^{-\frac{1}{d}}$ convergence rate (Roldán & Parrondo, ), where $n$ is the number of samples and $d$ is the dimension. In this paper, our spectral embedding method consistently utilizes the smallest $50$ eigenpairs for all experiments, unless explicitly stated otherwise.

### A.6 Measure DMD without PGMs

**Table 7: The stable (unstable) nodes identified based on traditional graph-based manifolds (Belkin & Niyogi, )**
| Attack | Embedded Matrix Cosine Similarity (Stable / Unstable) |
| --- | --- |
| DICE level 1 | 0.89 / 0.89 |
| Gaussian noise 1.0 | 0.96 / 0.96 |

As shown in Table 7, we cannot obtain meaningful DMDs with traditional graph-based manifolds (kNN graphs) adopted in the Laplacian Eigenmaps framework (Belkin & Niyogi, ).

### A.7 Metrics for Assessing GNN Stability

In the context of single-label classification in graph nodes, consider an output vector $\mathbf{y}=[y_{1},y_{2},...,y_{k}]$ corresponding to an input $\mathbf{x}$, where $k$ represents the total number of classes. The model’s predicted class is denoted as $\hat{y}=\text{argmax}_{i}(y_{i})$. Now, let’s assume that the output vector transforms to $\mathbf{y^{\prime}}=[y^{\prime}_{1},y^{\prime}_{2},...,y^{\prime}_{k}]$, while preserving the ordinality of the elements, i.e., if $y_{i}>y_{j}$, then $y^{\prime}_{i}>y^{\prime}_{j}$. This condition ensures that $\hat{y^{\prime}}=\text{argmax}_{i}(y^{\prime}_{i})=\hat{y}$.

Relying solely on model accuracy can be deceptive, as it is contingent upon the preservation of the ordinality of the output vector elements, even when the vector itself undergoes significant transformations. This implies that the model’s accuracy remains ostensibly unaffected as long as the ranking of the elements within the output vector is conserved. However, this perspective neglects potential alterations in the model’s prediction confidence levels.

In contrast, the cosine similarity provides a more holistic measure as it quantifies the angle between two vectors, thereby indicating the extent of modification in the output direction. This method offers a more granular insight into the impact of adversarial attacks on the model’s predictions.

Moreover, it is crucial to consider the nature of the output space, $Y$. In situations where $Y$ forms a probability distribution, a common occurrence in classification problems, the application of a distribution distance measure such as the Kullback-Leibler (KL) divergence is typically more suitable. Unlike the oversimplified perspective of accuracy, these measures can provide a nuanced understanding of the degree of perturbation introduced in the predicted probability distribution by an adversarial attack. This additional granularity can expose subtle modifications in the model’s output that might be missed when solely relying on accuracy as a performance metric.

### A.8 GNN Stability Analysis without GDR

In this study, we present the outcomes of stability quantification using original input and output graphs, as depicted in Figure 5. Our experimental findings underscore a key observation: SAGMAN without GDR does not allow for meaningful estimations of the GNN stability.

Figure: Figure 5: The horizontal axes, denoted by $X$, represent the magnitude of perturbation applied. ‘Random EdgesPT’ refers to the DICE adversarial attack scenario, in which pairs of nodes with different labels are connected and pairs with the same label are disconnected, with the number of pairs being equal to $X$. ‘Random FeaturePT’ indicates the application of Gaussian noise perturbation, expressed as $FM+X\eta$, where $FM$ denotes the feature matrix and $\eta$ represents Gaussian noise. The upper and lower subfigures illustrate the cosine similarity and the Kullback–Leibler Divergence (KLD). ‘Stable/Unstable’ denotes the samples that are classified as stable or unstable without GDR, respectively.
Refer to caption: x5.png

### A.9 Fast Effective-Resistance Estimation for LRD-based Graph Decomposition

The effective-resistance between nodes $(p,q)\in|V|$ can be computed using the following equation:

$$ $d^{eff}({p,q})=\sum\limits_{i=2}^{N}\frac{(u_{i}^{\top}e_{p,q})^{2}}{u_{i}^{ \top}L_{G}u_{i}},$ (27) $$

where $u_{i}$ represents the eigenvector corresponding to $\sigma_{i}$ eigenvalue of $L_{G}$ and $e_{p,q}=e_{p}-e_{q}$. To avoid the computational complexity associated with computing eigenvalues/eigenvectors, we leverage a scalable algorithm that approximates the eigenvectors by exploiting the Krylov subspace. In this context, given a nonsingular matrix $A_{N\times N}$ and a vector $c\neq 0\in\mathbb{R}^{N}$, the order-$(m)$ Krylov subspace generated by $A$ from $c$ is defined as:

$$ $\kappa_{m}(A,c):=span(c,Ac,A^{2}c,...,A^{m-1}c),$ (28) $$

where $c$ denotes a random vector, and $A$ denotes the adjacency matrix of graph $G$. We compute a new set of vectors denoted as $x^{(1)}$, $x^{(2)}$, …, $x^{(m)}$ by ensuring that the Krylov subspace vectors are mutually orthogonal with unit length. We estimate the effective-resistance between node $p$ and $q$ using Equation 27 by exploring the eigenspace of $L_{G}$ and selecting the vectors that capture various spectral properties of $G$:

$$ $d^{eff}({p,q})\approx\sum\limits_{i=1}^{m}\frac{(x^{(i)^{\top}}e_{p,q})^{2}}{x ^{(i)\top}L_{G}x^{(i)}},$ (29) $$

We control the diameter of each cycle by propagating effective resistances across multiple levels. Let $G=(V,E)$ represent the graph at the $\delta$-th level, and let the edge $(p,q)\in E$ be a contracted edge that creates a supernode $\vartheta\in V^{(\delta+1)}$ at level $\delta+1$. We denote the vector of node weights as $\eta^{(\delta)}\in\mathbb{R}_{\geq 0}^{V^{(\delta)}}$, which is initially set to all zeros for the original graph. The update of $\eta$ at level $\delta+1$ is defined as follows:

$$ $\eta_{\vartheta}:=\eta(p^{(\delta)})+\eta(q^{(\delta)})+d_{eff}^{(\delta)}(p,q).$ (30) $$

Consequently, the effective-resistance diameter of each cycle is influenced not only by the computed effective-resistance ($d_{eff}^{(\delta)}$) at the current level but also by the clustering information acquired from previous levels.

Figure: Figure 6: Graph decomposition results with respect to effective-resistance (ER) diameter
Refer to caption: extracted/5894680/figures/ER_fig.png

The graph decomposition results with respect to effective-resistance (ER) diameter are illustrated in Figure 6. The figure demonstrates that selecting a larger ER diameter leads to the decomposition of the graph into a smaller number of partitions, with more nodes included in each cluster. On the left side of the figure, the graph is decomposed into seven partitions: $P_{1},...,P_{7}$, by choosing a smaller ER diameter. Conversely, increasing the ER diameter on the right side of the figure results in the graph being partitioned into three clusters: $P_{1},P_{2},$ and $P_{3}$.

### A.10 Additional Results for GNN Stability Evaluation and Statistics of Datasets

We present the additional results in Figure 8, Figure 9, Figure 10, Figure 7, and Table 9. Table 8 summarizes the datasets utilized.

Figure: Figure 7: SAGMAN-selected stable/unstable samples for the "ogbn-arxiv" dataset. We report the GCN prediction accuracies under different levels of Gaussian Noise Perturbation $FM+X\eta$, where $FM$ denotes the feature matrix, $\eta$ represents Gaussian noise, and $X$ is the noise perturbation level.
Refer to caption: x6.png

**Table 8: Nettack adversarial attack targeting selected Cora samples in GCN**
| Nettack Level | Cosine Similarities:<br>Stable/Unstable |
| --- | --- |
| 1 | 0.99/0.90 |
| 2 | 0.98/0.84 |
| 3 | 0.97/0.81 |

Figure: Figure 8: Figures represent cosine similarities and KL divergence. "Random EdgesPT" corresponds to the DICE edge evasion attack. "Random FeaturePT" refers to Gaussian noise evasion perturbation $X+\xi\eta$, where $X$ is feature matrix, $\eta$ is Gaussian noise, $\xi$ is noise level controls
Refer to caption: x7.png

Figure: Figure 9: The horizontal axes, denoted by $X$, represent the magnitude of perturbation applied. ‘Random EdgesPT’ refers to the DICE adversarial attack scenario, in which pairs of nodes with different labels are connected and pairs with the same label are disconnected, with the number of pairs being equal to $X$. ‘Random FeaturePT’ indicates the application of Gaussian noise perturbation, expressed as $FM+X\eta$, where $FM$ denotes the feature matrix and $\eta$ represents Gaussian noise. The upper and lower subfigures illustrate the cosine similarity and the Kullback–Leibler Divergence (KLD). ‘SAGMAN Stable/Unstable’ denotes the samples that are classified as stable or unstable by SAGMAN, respectively.
Refer to caption: x10.png

Figure: Figure 10: The horizontal axes, denoted by $X$, represent the magnitude of perturbation applied. ‘Random EdgesPT’ refers to the DICE adversarial attack scenario, in which pairs of nodes with different labels are connected and pairs with the same label are disconnected, with the number of pairs being equal to $X$. ‘Random FeaturePT’ indicates the application of Gaussian noise perturbation, expressed as $FM+X\eta$, where $FM$ denotes the feature matrix and $\eta$ represents Gaussian noise. The upper and lower subfigures illustrate the cosine similarity and the Kullback–Leibler Divergence (KLD). ‘SAGMAN Stable/Unstable’ denotes the samples that are classified as stable or unstable by SAGMAN, respectively.
Refer to caption: x12.png

**Table 9: Summary of datasets used in our experiments**
| Dataset | Type | Nodes | Edges | Classes | Features |
| --- | --- | --- | --- | --- | --- |
| Cora | Homophily | 2,485 | 5,069 | 7 | 1,433 |
| Cora-ML | Homophily | 2,810 | 7,981 | 7 | 2,879 |
| Pubmed | Homophily | 19,717 | 44,324 | 3 | 500 |
| Citeseer | Homophily | 2,110 | 3,668 | 6 | 3,703 |
| Chameleon | Heterophily | 2,277 | 62,792 | 5 | 2,325 |
| Squirrel | Heterophily | 5,201 | 396,846 | 5 | 2,089 |
| ogbn-arxiv | Homophily | 169,343 | 1,166,243 | 40 | 128 |

### A.11 Why Generalized Eigenpairs Associate with DMD

(Cheng et al., ) propose a method to estimate the maximum distance mapping distortion (DMD), denoted as $\delta^{M}_{max}$, by solving the following combinatorial optimization problem:

$$ $\max{\delta^{M}}=\mathop{\max_{\forall p,q\in V}}_{p\neq q}\frac{e_{p,q}^{\top }L_{Y}^{+}e_{p,q}}{e_{p,q}^{\top}L_{X}^{+}e_{p,q}}$ (31) $$

When computing $\delta^{M}_{max}$ via effective-resistance distance, the stability score is an upper bound of $\delta^{M}_{max}$ (Cheng et al., ).

A function $Y=M(X)$ is called Lipschitz continuous if there exists a real constant $K\geq 0$ such that for all $x_{i}$, $x_{j}\in X$:

$$ $dist_{Y}(M(x_{i}),M(x_{j}))\leq Kdist_{X}(x_{i},x_{j}),$ (32) $$

where $K$ is the Lipschitz constant for the function $M$. The smallest Lipschitz constant, denoted by $K^{*}$, is called the best Lipschitz constant. Let the resistance distance be the distance metric, then (Cheng et al., ):

$$ $\lambda_{max}(L_{Y}^{+}L_{X})\geq K^{*}\geq\delta^{M}_{max}.$ (33) $$

Equation 33 indicates that the $\lambda_{max}(L_{Y}^{+}L_{X})$ is also an upper bound of the best Lipschitz constant $K^{*}$ under the low dimensional manifold setting. A greater $\lambda_{max}(L_{Y}^{+}L_{X})$ of a function (model) implies worse stability since the output will be more sensitive to small input perturbations. A node pair $(p,q)$ is deemed non-robust if it exhibits a large DMD, i.e., $\delta^{M}(p,q)\approx\delta^{M}_{\max}$. This suggests that a non-robust node pair consists of nodes that are adjacent in the $G_{X}$ but distant in the $G_{Y}$. To effectively identify such non-robust node pairs, the Cut Mapping Distortion (CMD) metric was introduced. For two graphs $G_{X}$ and $G_{Y}$ sharing the same node set $V$, let $S\subset V$ denote a node subset and $\bar{S}$ denote its complement. Also, let $cut_{G}(S,\bar{S})$ denote the number of edges crossing $S$ and $\bar{S}$ in graph $G$. The CMD $\zeta(S)$ of node subset $S$ is defined as (Cheng et al., ):

$$ $\zeta(S)\overset{\mathrm{def}}{=}\frac{cut_{G_{Y}}(S,\bar{S})}{cut_{G_{X}}(S, \bar{S})}.$ (34) $$

A small CMD score indicates that node pairs crossing the boundary of $S$ are likely to have small distances in $G_{X}$ but large distances in $G_{Y}$.

Given the Laplacian matrices $L_{X}$ and $L_{Y}$ of input and output graphs, respectively, the minimum CMD $\zeta_{\min}$ satisfies the following inequality:

$$ $\zeta_{\min}=\mathop{\min_{\forall S\subset V}}\zeta(S)\geq\frac{1}{\sigma_{ \max}(L_{Y}^{+}L_{X})}$ (35) $$

Equation 35 establishes a connection between the maximum generalized eigenvalue $\sigma_{\max}(L_{Y}^{+}L_{X})$ and $\zeta_{\min}$, indicating the ability to exploit the largest generalized eigenvalues and their corresponding eigenvectors to measure the stability of node pairs. Embedding $G_{X}$ with generalized eigenpairs. We first compute the weighted eigensubspace matrix $V_{s}\in{\mathbb{R}^{N\times s}}$ for spectral embedding on $G_{X}$ with $N$ nodes:

$$ $V_{s}\overset{\mathrm{def}}{=}\left[{v_{1}}{\sqrt{\sigma_{1}}},...,{v_{s}}{ \sqrt{\sigma_{s}}}\right],$ (36) $$

where $\sigma_{1},\sigma_{2},...,\sigma_{s}$ represent the first $s$ largest eigenvalues of $L_{Y}^{+}L_{X}$ and $v_{1},v_{2},...,v_{s}$ are the corresponding eigenvectors. Consequently, the input graph $G_{X}$ can be embedded using $V_{s}$, so each node is associated with an $s$-dimensional embedding vector. We can then quantify the stability of an edge $(p,q)\in E_{X}$ by measuring the spectral embedding distance of its two end nodes $p$ and $q$. Formally, we have the edge stability score defined for any edge $(p,q)\in E_{X}$ as $stability^{M}(p,q)\overset{\mathrm{def}}{=}\|V_{s}^{\top}e_{p,q}\|_{2}^{2}$
Let $u_{1},u_{2},...,u_{s}$ denote the first $s$ dominant generalized eigenvectors of $L_{X}L_{Y}^{+}$. If an edge $(p,q)$ is dominantly aligned with one dominant eigenvector $u_{k}$, where $1\leq k\leq r$, the following holds:

$$ $(u^{\top}_{i}e_{p,q})^{2}\approx\begin{cases}\alpha_{k}^{2}\gg 0&\text{ if }(i =k)\\ 0&\text{ if }(i\neq k).\end{cases}$ (37) $$

Then its edge stability score has the following connection with its DMD computed using effective-resistance distances (Cheng et al., ):

$$ $\|V_{s}^{\top}e_{p,q}\|_{2}^{2}\propto\left(\delta^{M}(p,q)\right)^{3}.$ (38) $$

The stability score of an edge $(p,q)\in E_{X}$ can be regarded as a surrogate for the directional derivative $\|\nabla_{v}M(x)\|$ under the manifold setting, where $v=\pm(x_{p}-x_{q})$. An edge with a larger stability score is considered more non-robust and can be more vulnerable to attacks along the directions formed by its end nodes.

Last, the node stability score can be calculated for any node (data sample) $p\in V$ as follows:

$$ $score(p)=\frac{1}{|\mathbb{N}_{X}(p)|}\sum_{q_{i}\in\mathbb{N}_{X}(p)}\left(\| V_{s}^{\top}e_{p,q}\|^{2}_{2})\right)\propto\frac{1}{|\mathbb{N}_{X}(p)|}\sum_ {q_{i}\in\mathbb{N}_{X}(p)}\left(\delta^{M}(p,q_{i})\right)^{3}$ (39) $$

where $q_{i}\in\mathbb{N}_{X}(p)$ denotes the $i$-th neighbor of node $p$ in graph $G_{X}$, and $\mathbb{N}_{X}(p)\in V$ denotes the node set including all the neighbors of $p$. The DMD score of a node (data sample) $p$ can be regarded as a surrogate for the function gradient $\|\nabla_{x}M(p)\|$ where $x$ is near $p$ under the manifold setting. A node with a larger stability score implies it is likely more vulnerable to adversarial attacks.

### A.12 Various Sampling Schemes for SAGMAN-guided Perturbations

Previous works (Cheng et al., ; Hua et al., ; Chang et al., ) highlighted that only part of the dataset plays a crucial role in model stability, so we want to focus on the difference between the most "stable" and "unstable" parts. However, it is certainly feasible to evaluate the entire graph. Table 10 shows the result regarding the Pubmed dataset in GPRGNN under Gaussian noise perturbation. Samples were segmented based on SAGMAN ranking, with the bottom 20% being the most "stable", the middle 60% as intermediate, and the top 20% representing the most "unstable". As anticipated, the "stable" category (representing the bottom 20%) should exhibit the lowest average KL divergences. This is followed by the intermediate category (covering the mid 60%), and finally, the "unstable" category (comprising the top 20%) should display the highest divergences.

**Table 10: KLD across varying Gaussian noise perturbations, expressed as $FM+X\eta$, where $FM$ denotes the feature matrix, $\eta$ represents Gaussian noise, and $X$ denotes the perturbation level. The dataset is divided into three segments based on the stability ranking of nodes as determined by SAGMAN.**
| Perturbation Level | KL divergence<br>(bottom 20%) | KL divergence<br>(mid 60%) | KL divergence<br>(top 20%) |
| --- | --- | --- | --- |
| 0.4 | 0.01 | 0.03 | 0.03 |
| 0.8 | 0.09 | 0.16 | 0.19 |
| 1.2 | 0.43 | 0.56 | 0.59 |

### A.13 SAGMAN’s Runtime Across Datasets

Figure 11 demonstrates the high efficiency of SAGMAN in processing large graphs, attributing its performance to the near-linear time complexity of its components.

Figure: Figure 11: SAGMAN’s runtime across datasets
Refer to caption: x13.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. In Figure 1, do the three graphs represent internal views of the same input graph within the framework, or are they independent examples with no connection to each other?

2. How about the efficiency of the proposed algorithm?



-----
After rebuttal:
Thanks for your effort for the reply, and I will keep my score. My main concern remains the insufficient discussion and comparsion of related works in graph robustness area [1]. Thanks for your understanding.

[1] A comprehensive survey on trustworthy graph neural networks: Privacy, robustness, fairness, and explainability.
soundness: 2
strengths: It is interesting to find another metric to measure the GNN stability, and cooperating with manifold via PGMs is natural.
summary: This paper proposes a complicated framework called SAGMAN to measure the stability of trained GNN model. 
It has three phases. 1. Construct the graph embedding based on node features and spectral features.
2. Use probabilistic graphical models (PGMs) to build low-dimension input graphs in the manifold. 
	At the same time, the node embeddings of GNN outputs are also used to build low-dimension graphs in another manifold via PGM.
3. Use the  distance mapping distortion (DMD) calculations to measure the GNN stability.
weaknesses: 1. Graph robustness is a widely discussed problems.
However, a comprehensive discussion of related works is needed in this literature.
Moreover, more baselines and benchmarks are needed to be provided in the experimental part.
Survey paper [1] can help you gain more understanding about more recent works.
[1] A comprehensive survey on trustworthy graph neural networks: Privacy, robustness, fairness, and explainability.

2. For this metric of GNN stability, is there any conclusion from the experiments? 
Does various graph robustness algorithms have been applied with this new metric?
Any comparison between this metric to other graph stability metric? For example, the prediction drop.

3. From my opinion, applying DMD on graphs is a problem as the manifold on graph is hard to define.
Although the PGMs can provide a solution, how good this PGMs can form the graph manifold is still a considerable questions.
Will the change of graph topology influence the mapping between graph structures and labels? This is also need to be analyzed.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: SAGMAN’s effectiveness depends on the input graph being well-represented in a low-dimensional space and having a significant eigengap. Is there any fix for this?
soundness: 3
strengths: The paper is well-structured and well-written.

SAGMAN can be used to quantify stability at the node level, which I think is quite useful.

The paper has a significant theoretical component, backed up by some experimental data as well.
summary: The authors propose a new framework based on spectral theory, called SAGMAN, for stability analysis of GNNs. They analyze the distance distortions between the input graph manifold and output graph manifold as a measure of stability. To this end, the authors propose a distance preserving graph dimensionality reduction method to obtain a low dimensional input/output graph manifolds and use the DMD metric to characterize the stability. The authors show that this approach is scalable by showing that the time complexity is near-linear in the graph size. In the experimental section, the authors show that SAGMAN can be used in recommendation systems. Further, the authors show that this framework can be used to facilitate adversarial targeted attacks.
weaknesses: The paper evaluates SAGMAN mainly on node classification and recommendation tasks. Testing on other GNN tasks, like link prediction or community detection, could strengthen the generalizability of the findings​

SAGMAN needs to be compared with other robustness techniques to get a better idea of what specific advantages it offers versus what are the limitations

I think that ignoring d (which in theory can be O(|V|) for a complete graph) from the complexity analysis and claiming that the method is near-linear in time complexity in the size of the graph is slightly misleading.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: In addition to the above weaknesses:
How the framework performs against more sophisticated graph structure and attribute attacks (Metattack, PGD). Could you perform analysis on adaptive attacks [R2].
How the model performs on OGBN datasets.
Empirical comparison with Lipschitz-based stability methods
soundness: 3
strengths: The paper studies an important problem of estimating the stability of GNNs
The framework is easy to understand and the computational complexity is linear w.r.t the graph dimensions.
Experimental evaluations validate the efficacy
summary: The paper concerns the stability analysis of GNNs. The idea is inspired from the previously proposed Distance Mapping Distortion (DMD) metric. To match the requirement of having a near low-dimensional manifold for both input and output the authors have proposed to perform graph dimensionality reduction step for estimating the probabilistic graphical model. Additionally for scalability requirements the authors have also proposed spectral sparsification.
weaknesses: The proposed framework although effective but largely uses a combination of many existing approaches.
The authors have clearly stated the distinction from the previous methods for analyzing GNN stability. However, the theoretical methods have guarantees. The current framework does not provide any strong theoretical guarantees. 
Theoretical results related to DMD could be extended however the authors need to show theoretical guarantees on their approach of estimating the low-dimensional manifolds.
The authors have mentioned that previous literature lack analysis on stability towards feature perturbations. However, the results in  [R1] could be seen as stability bounds of GNN w.r.t. Feature perturbations. The authors are suggested to include it in the discussion.
The framework is restricted only to evasion attacks and does not considers a practical perturbation setup where the input graph itself is poisoned

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: At the current stage, I have the following questions:

1. More content on how PGM acts on the graph data may be needed here. Also, if PGMs create a manifold for graphs, what is the metric of the manifold? Furthermore, does this (discrete) metric change over the GNN propagation? 

2. When presenting the notation of PGM and the reason for preserving the resistance distance is essential, the introduction of stable or unstable nodes needs to be enriched in Appendix A.6; I couldn't find any formal definition or threshold of identifying what exactly is stable/unstable nodes.

3. I think the stability induced by the algorithm is limited to the GNNs that smooth the node features. What will be the conclusion/influence for those GNNs that typically fit the heterophilic graphs or can induce a sharpening dynamic to the node features? As in this case, each row of X will tend to be dissimilar to each other. 

4. Perhaps a hyperparameter study for the algorithm with several components is necessary here.

5. Only to my view, as the proposed algorithm contains many components, they seems not be presented as a whole pipeline in the paper. Therefore, the author is suggested to refine the structure of the paper, e.g., adding some inter-contextual contents between each step. 

Some minor changes are suggested as follows: 
1. On row 234, the input embedding matrix $X = U_k$, I think $X$ should be the initial feature?

2. On row 759, the Appendix A.1 is self-refered.
soundness: 2
strengths: This paper well-organized the computational flow of inducing new DMD distance for evaluating GNN stability. 

The newly proposed algorithm has nearly linear time complexity and empirical study shows the GNNs under the guidance of the proposed model own greater stability compared to their original counterparts.
summary: This paper introduces a framework for evaluating the GNNs' stability. Specifically, the metrics for defining the so-called DMD distance are selected as effective resistance distance. The initial graph is first embedded via its spectral characteristics, and a secondary graph is constructed by the KNN algorithm followed by the graph pruning (sparsification) approach. Hence, for any pair of nodes, e.g., p and q, their distance quotient on the embedded manifold is leveraged as the evaluation metric for measuring GNN (node level) stability. Experimental studies show promising results via graph data under adversarial attacks.
weaknesses: Some of the contents are unclear, making the paper difficult to follow. 

Some technical parts of the paper may need further justification.