YOSO: You-Only-Sample-Once via Compressed Sensing for Graph Neural Network Training

## Abstract

Abstract Graph neural networks (GNNs) have become essential tools for analyzing non-Euclidean data across various domains. During training stage, sampling plays an important role in reducing latency by limiting the number of nodes processed, particularly in large-scale applications. However, as the demand for better prediction performance grows, existing sampling algorithms become increasingly complex, leading to significant overhead. To mitigate this, we propose YOSO (You-Only-Sample-Once), an algorithm designed to achieve efficient training while preserving prediction accuracy. YOSO introduces a compressed sensing (CS)-based sampling and reconstruction framework, where nodes are sampled once at input layer, followed by a lossless reconstruction at the output layer per epoch. By integrating the reconstruction process with the loss function of specific learning tasks, YOSO not only avoids costly computations in traditional compressed sensing (CS) methods, such as orthonormal basis calculations, but also ensures high-probability accuracy retention which equivalent to full node participation. Experimental results on node classification and link prediction demonstrate the effectiveness and efficiency of YOSO, reducing GNN training by an average of 75% compared to state-of-the-art methods, while maintaining accuracy on par with top-performing baselines.

## 1 Introduction

Graph Neural Networks (GNNs) (Kipf & Welling, ; Hamilton et al., ; Veličković et al., ; Chen et al., ; Chiang et al., ; Zou et al., ) have become pivotal in analyzing graph data across various domains, such as social network (Guo & Wang, ), protein interactions (Réau et al., ), and transportation systems (Liu et al., ). As graphs rapidly grow, long training time becomes a crucial factor impeding the wide utilization of GNNs in real world. To mitigate this issue, various sampling strategies such as node-wise (Hamilton et al., ; Chen et al., ), layer-wise (Chen et al., ; Zou et al., ; Huang et al., ), and subgraph-based methods (Chiang et al., ; Zeng et al., ) have been developed. These sampling strategies reduce the amount of data that required to sustain training and potentially shorten the training time. However, with the increasing complexity of sampling algorithms, GNNs have struggled to maintain training efficiency in large-scale applications, such as IGB dataset (Khatua et al., ).

Theoretically, the model accuracy loss caused by sampling algorithms stems from the bias and variance introduced by estimating the overall data based on the samples (Huang et al., ). Unlike the unbiased and variance-free GCN (Kipf & Welling, ; Huang et al., ) that utilize all training nodes, low time complexity sampling methods struggle to accurately estimate both graph structure and embeddings (Jin et al., ), potentially degrading the outcomes. As a result, recent sampling algorithms focus solely on improving accuracy while overlooking high computational cost, have become increasingly complex. This highlights a significant gap as shown in Figure 1(a): finding a method that achieves both high accuracy and efficiency.

To reveal the large overhead introduced by sampling in GNN training, we conduct empirical evaluations for state-of-the-art (SOTA) sampling schemes with Reddit dataset (Hamilton et al., ). As shown in Figure 1(b), we break down the total training time into three non-overlapping components: (1) Sampling, (2) Mem2GPU: refers to transferring data to GPU memory, and (3) Computation: all processes on GPU. Our results indicate that sampling stage occupies 35.7% to 64% of the total training time across various sampling algorithms, making it a significant overhead when considering both training efficiency and model accuracy. For instance, as a representative of layer-wise, AS-GCN (Huang et al., ) spends 55.6% of the total training time on sampling but only achieves suboptimal accuracy. Subgraph-based sampling methods, although achieving the highest model accuracy, incur the most significant overhead with sampling stage accounting for up to 64% of the total training time. Node-wise sampling falls between these two paradigms in terms of overall performance. For example, VR-GCN spends 685.72 seconds on sampling and achieves a Micro-F1 score of 0.962.

To fill identified research gap, we propose YOSO (You-Only-Sample-Once), a novel approach applies compressed sensing (CS) (Candes & Tao, ) technique to GNN sampling. YOSO reimagines the feature matrix as multi-channel signal and utilize adapted CS to reduce the amount of computation involved in the training by transferring the feature matrix to another domain with high sparsity. YOSO enables training with only $M$ nodes sampled from the graph with $N$ nodes where $M\ll N$, followed by a nearly lossless reconstruction which guarantees the model accuracy closely aligns with zero bias and variance as if all $N$ nodes were used for training. Moreover, sampling in YOSO is designed to occur only once at the beginning of the training. This involves determining the sampling matrix, represented as $\mathbf{\Phi}$, based on the characteristics of the graph dataset. Subsequently, the reconstruction process takes place after each forward propagation, achieved by integrating
the reconstruction process with the loss function of specific learning tasks, i.e., cross-entropy loss in node classification, to guide the backward propagation. Thus, YOSO streamlines the entire training process by eliminating the need for continuous resampling through entire training and ensures every step of learning is informed by an optimally reconstructed data state, significantly enhancing both the efficiency and efficacy of the model training. We summarize our contributions below.

- •
We propose a novel sampling method called YOSO, which significantly reduces GNN training time by performing one-time sampling for the entire training while maintaining strong prediction accuracy through a nearly lossless reconstruction of the embedding matrix.
- •
YOSO eliminates the need for expensive computations typically associated with combining CS with GNN sampling, thereby making the sampling process highly efficient.
- •
Experimental results demonstrate the effectiveness of YOSO on both node classification and link prediction tasks. Specifically, YOSO significantly reduces overall training time by an average of around 75% while preserving model accuracy. Ablation studies further reveal that YOSO achieves near-zero bias and variance, effectively reconstructing the embedding matrix with minimal error.

Figure: (a) Total training time v.s. model accuracy
Refer to caption: extracted/5987727/motivation2.png

## 2 Preliminaries

Graph Neural Networks (GNNs) operate on a graph, represented as $G=\{V,E,\hat{\mathbf{A}},\mathbf{X}\}$, where $V=\{1,2,\dots,N\}$ represents the set of nodes, $E=\{(i,j)\mid i,j\in V\}$ defines the edges, and $\hat{\mathbf{A}}\in\mathbb{R}^{N\times N}$ is a matrix that encoding the connection properties between nodes, i.e., adjacency matrix or normalized Laplacian matrix. $\mathbf{X}\in\mathbb{R}^{N\times d}$ is feature matrix, where $d$ is the feature dimension. GNNs aim to learn node embeddings $\mathbf{H}^{(l)}$ through the layer-specific transformations governed by parameters $\theta^{(l)}$, expressed as $\mathbf{H}^{(l)}=f_{\theta^{(l)}}(\mathbf{H}^{(l-1)},\hat{\mathbf{A}}),l=1,2,.
..,L$, where $L$ represents the number of layers, with the initial embedding $\mathbf{H}^{(0)}=\mathbf{X}$.

Sampling is employed in GNNs to manage computational complexity, where a subset $V^{{}^{\prime}}\subset V$ of nodes is selected based on certain sampling rules $\mathcal{P}$, such as importance sampling and Monte Carlo estimation (Chen et al., ). Thus, The embeddings are estimated as $f_{\theta^{(l)}}(\mathbf{H}_{[V^{{}^{\prime}}]}^{(l-1)},\hat{\mathbf{A}})$, where $[V^{{}^{\prime}}]$ denotes the indices corresponding to $V^{{}^{\prime}}$, reducing the data need to be processed. However, existing sampling algorithms have increasingly complicated the computation of $\mathcal{P}$ to achieve more accurate approximations, leading to a growing overhead in sampling time.

Compressed Sensing (CS) is originally developed for sampling and reconstructing signal (Candes & Tao, ). This technique can be applied to $\mathbf{H}^{(l)}$ if it exhibits sparsity in some specific transform domain. If $\mathbf{H}^{(l)}=\mathbf{U}^{(l)}\hat{\mathbf{H}}^{(l)}$, where $\hat{\mathbf{H}}^{(l)}$ is sparse, i.e., $\hat{\mathbf{H}}$ contains at most $k$ non-zero rows, noted as $\|\hat{\mathbf{H}}^{(l)}\|_{0,row}\leq k$. The indices corresponding to the dense rows in $\hat{\mathbf{H}}^{(l)}$ is called support, which holds the most significant information and can be used to effectively reconstruct $\mathbf{H}^{(l)}$. The orthonormal basis $\mathbf{U}^{(l)}\in\mathbb{R}^{N\times N}$ is used to transform $\mathbf{H}^{(l)}$ into the sparse domain: $\hat{\mathbf{H}}^{(l)}=\mathbf{U}^{T}\mathbf{H}^{(l)}$. Its existence is a necessary condition of CS. Fortunately, the orthonormal basis $\mathbf{U}^{(l)}$ that satisfying $\mathbf{H}^{(l)}=\mathbf{U}^{(l)}\hat{\mathbf{H}}^{(l)}$ where $\|\hat{\mathbf{H}}^{(l)}\|_{0,row}\leq k$, always exists (Isufi et al., ; Bo et al., ) and can be derived from the graph’s structural properties (Tsitsvero et al., ; Puy et al., ; Chen et al., ), i.e., normalized Laplacian matrix. Let $\mathbf{T}^{(l)}\in\mathbb{R}^{M\times d}$ where $M\ll N$, be the measurement matrix, computed as:

$$ $\mathbf{T}^{(l)}=\mathbf{\Phi}^{(l)}\mathbf{U}^{(l)}\hat{\mathbf{H}}^{(l)}$ (1) $$

Here, $\mathbf{\Phi}^{(l)}\in\mathbb{R}^{M\times N}$ is known as the sampling matrix. The measurement matrix represents the specific numerical values that can be directly observed during the computational process. To reconstruct the original sparse $\hat{\mathbf{H}}^{(l)}$, the following optimization problem need to be solved:

$$ $\tilde{\mathbf{H}}^{(l)}=\text{min}_{\hat{\mathbf{H}}^{(l)}}\|\hat{\mathbf{H}} ^{(l)}\|_{2,1}\quad\text{subject to}\quad\mathbf{T}^{(l)}=\mathbf{\Phi}^{(l)} \mathbf{U}^{(l)}\hat{\mathbf{H}}^{(l)}$ (2) $$

where $\|\cdot\|_{2,1}$ is $l_{2,1}$ norm (Liu et al., ). Accurate reconstruction requires that the matrix $\mathbf{\Phi}^{(l)}\mathbf{U}^{(l)}$ satisfies the Restricted Isometry Property (RIP) (Candes & Tao, ), formulated as

$$ $(1-\delta_{k})\|\hat{\mathbf{H}}^{(l)}\|^{2}_{F}\leq\|{\mathbf{\Phi}^{(l)} \mathbf{U}^{(l)}\hat{\mathbf{H}}^{(l)}}\|^{2}_{F}\leq(1+\delta_{k})\|\hat{ \mathbf{H}}^{(l)}\|^{2}_{F}$ (3) $$

where $0<\delta_{k}<1$, and $\|\cdot\|_{F}$ is the Frobenius norm. After obtaining $\tilde{\mathbf{H}}^{(l)}$ through Equation (2), the original $\mathbf{H}^{(l)}$ can be reconstructed as:

$$ $\mathbf{H}^{(l)}=[\mathbf{U}^{(l)}]^{T}\tilde{\mathbf{H}}^{(l)}$ (4) $$

## 3 Compressed Sensing as Sampling in GNNs

As discussed in Section 2, CS reduces the amount of data required for computation by transforming $\mathbf{H}^{(l)}\in\mathbb{R}^{N\times d}$ into a much smaller matrix $\mathbf{T}\in\mathbb{R}^{M\times d}$ since $M\ll N$. This reduction depends on converting $\mathbf{H}^{(l)}$ into a sparse domain where its basis is $\mathbf{U}$, resulting in $\hat{\mathbf{H}}^{(l)}$, such that $\|\hat{\mathbf{H}}^{(l)}\|_{0,row}\leq k$. Efficient reconstruction is possible if such orthonormal basis $\mathbf{U}$ and sampling matrix $\mathbf{\Phi}$ exist and satisfy the RIP (Equation (3)).

When applied to GNNs, CS offers two main advantages over other schemes: (1) $\mathbf{H}\in\mathbb{R}^{N\times d}$ can be sampled into a much smaller $\mathbf{T}\in\mathbb{R}^{M\times d}$, significantly reducing computation time while retaining essential information; (2) CS enables lossless reconstruction at the output layer, allowing $\mathbf{T}$ to be accurately expanded back to $\mathbf{H}$ with high probability, as if all nodes were involved in the computation. Thus, a smaller sampled set can emulate the full training set, achieving high accuracy and reduced sampling time. This lossless property ensures that the model retains all information, thereby enhancing accuracy. Specifically:

$$ $\mathbf{H}^{(l)}=f_{\theta^{(l)}}\left(Rec\left\{\mathbf{T}^{(l-1)}\right\}, \hat{\mathbf{A}}\right)$ (5) $$

where $Rec\{\cdot\}$ represents the processing of reconstruction (Equation (2) and Equation (4)). However, the iterative processes in Equation (5) is highly inefficient and has the following challenges:

- •
Expensive Computations of $\mathbf{U}^{(l)}$ and $\mathbf{\Phi}^{(l)}$. Determining appropriate orthonormal bases $\mathbf{U}^{(l)}$ and sampling matrices $\mathbf{\Phi}^{(l)}$ for $l=1,...,L$, is time-consuming. While Section 2 theoretically confirms the existence of $\mathbf{U}^{(l)}$, practical computation is costly since it requires matrix decompositions with an average time complexity of $O(n^{3})$. More seriously, $\mathbf{H}^{(l)}$ changes across GNN layers, therefore a fixed $\mathbf{U}^{(l)}$ is unlikely to meet the sparsity requirements for all layers, which necessitating $(L+1)$ separate decompositions. Similarly, $\mathbf{\Phi}^{(l)}$ must adapt to changes in $\mathbf{U}^{(l)}$ to maintain RIP, requiring an additional $(L+1)$ adjustments. In summary, determining $\mathbf{U}^{(l)}$ and $\mathbf{\Phi}^{(l)}$ involves $2(L+1)O(n^{3})$ costly computations during training.
- •
Accurate but Time-inefficient Reconstruction. As in Equation (5), to minimize error propagation, we reconstruct $\mathbf{H}$ at every layer before proceeding to the next layer. However, this incurs significant computational overhead. The fastest known reconstruction algorithm has an average time complexity of $O(nm)$ (Maleki, ), where $n$ is the signal dimension and $m$ is the measurement length. For GNNs, this translates to an average reconstruction time complexity of $O(dM)$ per layer, resulting in a total cost of $O(dML)$ for an $L$-layer GNN. Such overhead greatly reduces training efficiency.

Consequently, directly applying CS to GNN sampling introduces significant time complexities. To effectively integrate CS into GNNs and ensure its efficiency, we must overcome the two obstacles:

- I.
Working with Unknown $\mathbf{U}^{(l)}$ and Universal $\mathbf{\Phi}$. Given the high computational cost of determining $\mathbf{U}^{(l)}$, we need to satisfy or approximate CS’s necessary and sufficient condition without explicitly knowing $\mathbf{U}^{(l)}$. Without $\mathbf{U}^{(l)}$, identifying the support and determining essential nodes for reconstruction becomes challenging, complicating the construction of $\mathbf{\Phi}^{(l)}$. Since $\mathbf{\Phi}^{(l)}$ is layer-specific, calculating it for each layer is impractical. Thus, we require a method that works with an unknown $\mathbf{U}$ using a universal sampling matrix $\mathbf{\Phi}$, ensuring $\mathbf{\Phi}$ remains adaptable to any $\mathbf{U}$ while satisfying compressed sensing conditions.
- II.
Balancing computational efficiency with the need for accurate reconstruction. If we sample once at the input layer and use these results throughout the GNN computation, followed by reconstruction only at the output layer, this approach requires just one sampling and reconstruction step for the entire training process. Although it may introduce some accuracy loss due to reduced intermediate layer information, it remains efficient if this loss is controllable with a known upper bound, allowing a balance between computational efficiency and model accuracy.

## 4 Methodology

The overall YOSO algorithm is presented in Section 4.1, where we also address the challenge of working with the unknown $\mathbf{U}$. Followed by design of universal sampling matrix $\mathbf{\Phi}$ in Section 4.2.

### 4.1 YOSO Design

Figure: Algorithm 1 Forward and Backward Propagation of YOSO

YOSO proposes a CS-based sampling and reconstruction framework, where nodes are sampled once at the input layer, followed by a lossless reconstruction at the output layer during each epoch. As shown in Algorithm 1, the entire training process of YOSO consists of forward propagation, loss computation, and backward propagation, similar to the conventional GNN training. Unlike the standard process, YOSO operates within a specific sparse domain instead of the original data domain. Initially, the data is transformed into the sparse domain (Line 3), where the one-time sampling is also performed by using the sampling matrix $\mathbf{\Phi}$. The subsequent steps–forward propagation (Lines 4-7), loss computation (Lines 8-10), and backward propagation (Lines 11-17)–are all executed within this sparse domain. The detailed description is as follows:
One time sampling (Line 3). Given a graph $G=\{V,E,\hat{\mathbf{A}},\mathbf{X}\}$, where specific $\hat{\mathbf{A}}$ is the normalized Laplacian matrix. We perform the sampling stage only once using the sampling matrix $\mathbf{\Phi}$ on the sparsity domain $\hat{\mathbf{X}}$ as $\mathbf{\Phi}\mathbf{U}\hat{\mathbf{X}}$, resulting in $\mathbf{T}^{(0)}\in\mathbb{R}^{M\times d}$, where $M\ll|V|=N$. This process involves the construction of the sampling matrix $\mathbf{\Phi}$, for details, please refer to Section 4.2.
Forward propagation (Lines 4-7). The forward propagation of YOSO can be expressed as:

$$ $\left\{\begin{array}[]{lr}\mathbf{T}^{(l)}=\sigma\left(\mathbf{\Phi}\hat{ \mathbf{A}}\mathbf{W}^{(l)}\mathbf{T}^{(l-1)}\right)&1\leq l\leq L-1\\ \mathbf{U},\hat{\mathbf{H}}^{(L)}=Rec\left\{\sigma\left(\mathbf{\Phi}\hat{ \mathbf{A}}\mathbf{W}^{(L)}\mathbf{T}^{(L-1)}\right)\right\}&l=L\end{array}\right.$ (6) $$

where $\sigma(\cdot)$ is the activation function, $\mathbf{W}^{(l)},l=1,...,L$ is the $l-$th layer’s trainable parameters, $\mathbf{U}$ is the unknown orthonormal basis and the method for addressing this (working with unknown $\mathbf{U}$) will be discussed in the following.
Loss function and working with unknown $\mathbf{U}$ (Lines 8-10). First, we discuss the construction of YOSO’s loss function in (1), and then in (2), we explain why this construction effectively addresses the challenge of working with the unknown $\mathbf{U}$.
(1) Loss function. The $Rec\{\cdot\}$ in Equation (6) is equal to solve the following optimization problem:

$$ $\underset{\hat{\mathbf{H}}^{(L)},\mathbf{U}}{\min}\ \frac{1}{2}\left\|\mathbf{ Z}-\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}^{(L)}\right\|_{F}^{2}+\lambda\left \|\hat{\mathbf{H}}^{(L)}\right\|_{2,1}\quad\text{s.t. }\mathbf{U}\mathbf{U}^{T }=\mathbf{U}^{T}\mathbf{U}=\mathbf{I}$ (7) $$

where $\mathbf{Z}=\sigma(\mathbf{\Phi}\hat{\mathbf{A}}\mathbf{W}^{(L)}\mathbf{T}^{(L-
1)})$ represents the sampled measurement matrix at output layer, and $\lambda$ is a hyperparameter controlling the balance between data fidelity and sparsity. Equation (7) is a non-trivial optimization problem involving both $\hat{\mathbf{H}}^{(L)}$ and $\mathbf{U}$ due to non-convexity introduced by orthogonality constraint ($\mathbf{U}\mathbf{U}^{T}=\mathbf{U}^{T}\mathbf{U}=\mathbf{I}$) and the interaction between variables. To overcome it, we perform joint optimization of Equation (7) with the GNN’s specific loss function (e.g., cross-entropy). Let the GNN’s loss function be $\mathcal{L}_{GNN}^{\Theta}$, where $\Theta=\{\mathbf{W}^{(1)},...,\mathbf{W}^{(L)}\}$ represents the set of all trainable parameters. The joint optimization objective function is defined as:

$$ $\underset{\hat{\mathbf{H}}^{(L)},\mathbf{U},\Theta}{\min}\left\{\alpha\left( \frac{1}{2}\left\|\mathbf{Z}-\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}^{(L)} \right\|_{F}^{2}+\lambda\left\|\hat{\mathbf{H}}^{(L)}\right\|_{2,1}\right)+ \beta\mathcal{L}_{GNN}^{\Theta}(\mathbf{Z})\right\}\quad\text{s.t. }\mathbf{U} \mathbf{U}^{T}=\mathbf{U}^{T}\mathbf{U}=\mathbf{I}$ (8) $$

where $\alpha$ and $\beta$ is the hyperparameters to balance the reconstruction loss and GNN loss.
(2) Working with unknown $\mathbf{U}$. To address the challenge of unknown $\mathbf{U}$, we treat $\mathbf{U}$ as an optimization target. Using Equation (8), we obtain a total loss, which is then used to generate gradients for updating $\mathbf{U}$ through all training process (Detailed calculation of the gradient of the loss in Equation (8) with respect to $\mathbf{U}$ can be found in Appendix B.1).
Backward Propagation (Lines 11-17). The backward propagation process uses the loss generated by Equation (8) to update three parameters, which are $\mathbf{U}$, $\hat{\mathbf{H}}^{(L)}$, and $\Theta=\{\mathbf{W}^{(1)},\dots,\mathbf{W}^{(L)}\}$ through gradient descent. This process results in three gradients, namely $\nabla_{\mathbf{U}}\mathcal{L}$, $\nabla_{\hat{\mathbf{H}}^{(L)}}\mathcal{L}$, and $\nabla_{\Theta}\mathcal{L}$, each corresponding to three learning rates $\eta_{\mathbf{U}}$, $\eta_{\hat{\mathbf{H}}^{(L)}}$, and $\eta_{\Theta}$, respectively. For the detailed setting of hyperparameters used here, i.e., $\alpha$ and $\eta_{\mathbf{U}}$, please refer to Appendix A.4 and the detailed gradient computation list in Appendix B.1.

Through Algorithm 1, we obtain both $\mathbf{U}$ and $\hat{\mathbf{H}}^{(L)}$. With $\mathbf{U}$ now determined, we can apply Equation (4) to reconstruct $\mathbf{H}^{(L)}$, which can then be utilized for downstream tasks, such as link prediction. Compared to Equation (5), the process described in Algorithm 1 trades some accuracy for improved efficiency, and importantly, this accuracy loss is bounded. For detailed statements and proofs, please refer to Appendix B.4.

### 4.2 Construction of Sampling Matrix 𝚽 𝚽 \mathbf{\Phi} bold_Φ

When the orthonormal basis $\mathbf{U}$ remains unspecified before training, we encounter the challenge of computing $\mathbf{T}^{(0)}=\mathbf{\Phi}\mathbf{U}\mathbf{\hat{X}}$ in Equation (6) since the absence of knowledge about $\mathbf{U}$ complicates the design of $\mathbf{\Phi}$. In traditional CS, $\mathbf{U}$ maps data into a sparse domain where the support (i.e., the indices of non-zero rows) is clearly identifiable, and these non-zero rows contain the crucial information. This clarity allows $\mathbf{\Phi}$ to be designed in a targeted manner based on the support. Without knowing $\mathbf{U}$, it becomes challenging to design a $\mathbf{\Phi}$ that effectively captures the essential information. Therefore, the main difficulty lies in designing an effective and universal sampling matrix $\mathbf{\Phi}$ that not only accurately captures the essential characteristics of the graph data but also works with any $\mathbf{U}$ without violating the RIP.

To address challenge, we propose an approach that integrates the design of a matrix $\hat{\mathbf{S}}\in\mathbb{R}^{M\times N}$, derived from the graph structure, with the construction of the sampling matrix $\mathbf{\Phi}$, i.e., $\mathbf{\Phi}=\hat{\mathbf{S}}\otimes\mathbf{\Sigma}$ where $\mathbf{\Sigma}\in\mathbb{R}^{M\times N}$ is a random matrix and $\otimes$ is element-wise production.

$\hat{\mathbf{S}}$ remains unchanged through the entire training and is determined only once during pre-processing phase. Design of $\hat{\mathbf{S}}$ is graph-structure-based for two reasons: first, the graph structure is invariant, and second, it reflects the importance of certain nodes, which is crucial for the GNN message-passing process. For the sampling matrix $\mathbf{\Phi}$, it is essential to be row full rank. Intuitively, $\mathbf{\Phi}$ serves to linearly combine the features or embeddings of nodes according to weights corresponding to the indices of the support (non-zero rows). If $\mathbf{\Phi}$ is row over-ranked, it results in redundant information, whereas a row under-ranked $\mathbf{\Phi}$ leads to information loss. Thus, ensuring a row full-rank sampling matrix is crucial for effectively capturing the necessary information.
Construction of $\hat{\mathbf{S}}$. Considering the normalized Laplacian matrix $\hat{\mathbf{A}}=\mathbf{I}-\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}$ where $\mathbf{D}$ and $\mathbf{A}$ are the degree matrix and adjacency matrix, respectively. The $N$ nodes correspond to $N$ eigenvalues from $\hat{\mathbf{A}}$’s spectral decomposition, denoted as $\{\lambda_{1},\ldots,\lambda_{N}\}$ and $\lambda_{i}\geq 0$ for any $i$ holds. These eigenvalues often reflect the important structural properties of the graph. For example, larger eigenvalues correspond to more influential nodes within the graph. To construct the sampling probability distribution, we define $P(i)=\frac{\lambda_{i}}{\sum_{j=1}^{N}\lambda_{j}}$, where node $i$ has a sampling probability proportional to its eigenvalue relative to the total eigenvalue sum. Using this probability distribution, we sample $M$ times to form the $M$ rows of $\hat{\mathbf{S}}$. Suppose node $i$ is sampled; the corresponding row in $\hat{\mathbf{S}}$ will include node $i$’s 1-hop neighbors. Assume node $i$ has $N(i)$ neighbors, each neighbor is randomly sampled with a probability of $\frac{1}{N(i)}$. This construction ensures that $\hat{\mathbf{S}}$ will not contain any all-zero rows, thanks to the self-loop added by the normalized Laplacian. Consequently, the matrix $\mathbf{\Phi}=\hat{\mathbf{S}}\mathbf{\Sigma}$ will be row full rank (detailed proof in Appendix B.2), avoiding any issues with row rank deficiency.
Construction of $\mathbf{\Sigma}$. Some studies have highlighted the importance of randomness in achieving the RIP (Baraniuk et al., ). Therefore, we define $\Sigma$ as a random matrix. Intuitively, since we do not have precise knowledge of the support, we randomly sample $M$ nodes based on eigenvalue weights to estimate the support. The matrix $\mathbf{\Sigma}$ should reflect the contribution level of each node $i$ to the non-zero rows (i.e., the support). For instance, if node $k$ is shared by both nodes $i$ and $j$, we need to determine how much node $k$ contributes to node $i$ and to node $j$. This is crucial for ensuring accurate reconstruction and satisfying the RIP. For any column $j$ in $\hat{\mathbf{S}}$, assume it contains $g(j)$ non-zero elements. We assign the corresponding elements in $\Sigma$ as random values drawn from a Gaussian distribution $N(0,\frac{1}{g(j)})$. This design helps capture the contribution levels effectively, which is important for achieving the Restricted Isometry Property (detailed proof in Appendix B.3).

## 5 Related Work

A widely accepted criterion (Liu et al., ) divides current different sampling methods into three categories: node-wise sampling, layer-wise sampling, and subgraph-based sampling, depending on the granularity of the sampling operation during mini-batch generation.
Node-wise Sampling: This fundamental approach, pioneered by works such as GraphSage (Hamilton et al., ) and others (Ying et al., ; Chen et al., ; Dai et al., ), involves sampling at the individual node level. Each node’s neighbors are selected according to specific probabilities, often using a uniform distribution. For example, GraphSage samples $k-$hop neighbors at varying depths, with the sampling size, or fanout, for each depth tailored to optimize model performance. This approach, while simple and effective, has been criticized for its exponential increase in sampling time complexity as the number of GNN layers grows.
Layer-wise Sampling: Developed to address the exponential growth in computational complexity as network depth increases in node-wise sampling, this method samples multiple nodes simultaneously in one step. Techniques like FastGCN (Chen et al., ) reframe GNN loss functions as integral transformations and utilize importance sampling and Monte-Carlo approximation to manage variance. Further developments, such as AS-GCN (Huang et al., ) and LADIES (Zou et al., ), focus on maintaining sparse connections between sampled nodes to aid convergence. However, these methods tend to introduce additional complexity and computational cost.
Subgraph-based Sampling: These methods form mini-batch training subgraphs using graph partitioning algorithms. Cluster-GCN (Chiang et al., ) partitions the full graph into clusters, sampling these clusters to create subgraphs for training batches. GraphSAINT (Zeng et al., ) dynamically estimates sampling probabilities for nodes and edges to form subgraphs over which the full GNN model is trained. While these techniques typically improve model accuracy, they also lead to longer training time.

## 6 Experiments

In Section 6.2, we evaluate the training time along with model accuracy across two learning tasks: node classification and link prediction. Also, to investigate convergence performance, we assess the convergence of both the baselines and YOSO in Section 6.3. Finally, we conduct an ablation study on the proposed compensations in Section 6.4. Details on the dataset, baselines, experimental hardware and software configuration can be found in Section 6.1 and Appendix A.1.

### 6.1 Experimental Settings

Datasets. For the node classification task, we selected Reddit (Hamilton et al., ), ogbn-arxiv and ogbn-products (Hu et al., ). For the link prediction task, we used ogbl-ppa, and ogbl-citation2 (Hu et al., ). For detailed dataset statistics, data splits and metrics, please refer to Appendix A.2.
Baselines and Implementation. The baselines used in this paper include node-wise sampling methods (GraphSage (Hamilton et al., ) and VR-GCN (Chen et al., )), layer-wise sampling methods (FastGCN (Chen et al., ), AS-GCN (Huang et al., ) and LADIES (Zou et al., )) and subgraph-based sampling methods (Cluster-GCN (Chiang et al., ) and GraphSAINT (Zeng et al., )). Notably, several baseline models lacked implementations for link prediction, prompting us to modify them accordingly. Detailed information on the source code for these baselines, the YOSO implementation, and other related materials can be found in Appendix A.3.
Hyperparameter Setting. All experiments are conducted using a two-layer GNN. Detailed hyperparameter settings are described in Appendix A.4.

### 6.2 Overall Comparison

In this section, we evaluate baselines and YOSO with two key metrics: model accuracy (varies with different datasets and tasks) and total training time. The training time is broken down into three non-overlapped parts: Sampling, Mem2GPU, and Computation.

Node Classification Task:
First, YOSO achieves the shortest total training time with an average of 75.3% reduction across all datasets compared to all baselines as shown in Figure 2. For example, YOSO reduces around 95% total training time from 233.22 seconds (ogbn-arxiv/AS-GCN) and 12,387.2 seconds (ogbn-products/AS-GCN) to 199.02 and 8,013.23 seconds, respectively. The main reason is that YOSO significantly reduces the sampling time while introducing a little re-construction overhead. As shown in Figure 2(a)-(c), the most substantial sampling time reduction occurs on the Reddit dataset, where YOSO achieved a 99% decrease, cutting the sampling time from 1149.02 seconds for GraphSAINT-EDGE and 1107.54 seconds for Random Walk to just 15.13 seconds. On average, YOSO reduced sampling time by approximately 95.7% compared to all other baselines.

For model accuracy shown in Table 1, YOSO consistently matches or closely approaches the top performers. For example, YOSO obtains an accuracy of 0.71 on ogbn-arxiv, just 0.01 below GraphSage. On Reddit, it achieves the highest score of 0.967, matching GraphSAINT-Random Walk, and on ogbn-products, it reaches 0.787, slightly trailing GraphSAINT-EDGE’s 0.792.

Figure: (a) ogbn-arxiv
Refer to caption: x1.png

**Table 1: Model accuracy results for different sampling schemes on node classification and link prediction tasks. For specific evaluation metrics on each dataset, please refer to Table 2.**
| Different<br>Sampling<br>Schemes | Dataset |  |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Node Classification | Link Prediction |  |  |  |  |
| ogbn-arxiv | Reddit | ogbn-products | ogbl-ppa | ogbl-citation2 |  |
| GraphSage | 0.72 | 0.949 | 0.772 | 0.1704 | 0.8054 |
| VR-GCN | 0.697 | 0.962 | 0.699 | 0.1704 | 0.7967 |
| FastGCN | 0.438 | 0.927 | 0.404 | 0.1088 | 0.6555 |
| AS-GCN | 0.687 | 0.964 | 0.51 | 0.1245 | 0.6593 |
| LADIES | 0.649 | 0.927 | 0.501 | 0.1131 | 0.6693 |
| Cluster-GCN | 0.653 | 0.966 | 0.769 | 0.2053 | 0.7904 |
| GraphSAINT-EG | 0.702 | 0.967 | 0.792 | 0.2143 | 0.8039 |
| GraphSAINT-RW | 0.701 | 0.967 | 0.783 | 0.2263 | 0.8054 |
| YOSO | 0.72 | 0.967 | 0.787 | 0.2238 | 0.8025 |

Link Prediction Task:
For total training time, similar to the node classification task, YOSO achieves the best training time with a 72.13% average training time decrease across all datasets for the link prediction. For example, YOSO decreases the training time for the ogbl-ppa dataset from 44.53 seconds with AG-GCN to 21.42 seconds, and for the ogbl-citation2 dataset from 8423.06 seconds with AG-GCN to 455.35 seconds.This improvement is consistent with the node classification task, where YOSO achieves considerable reductions in sampling time while introducing minimal reconstruction overhead. As depicted in Figure 2(d)-(e), YOSO achieves an average sampling time reduction of about 80.5% across all datasets. As for model accuracy, outlined in Table 1, YOSO maintained results with only a very small gap–0.0025 on ogbn-arxiv and 0.0029 on ogbl-citation2–compared to the best results achieved by GraphSAINT-Random Walk and GraphSage, respectively.

In summary, for both tasks of node classification and link prediction, by combining high accuracy with substantial reductions in sampling and total training time, YOSO demonstrates its efficiency in GNN training and significantly improves both sampling and total training times across all datasets while maintaining competitive accuracy, highlighting its effectiveness compared to the baselines on the node classification task.

### 6.3 Convergence Comparison

We investigate YOSO’s convergence performance compared to other baselines. Specifically, we select ogbn-arxiv and ogbl-ppa as representatives for node classification and link prediction, respectively. The training loss-epoch curves are shown in Figure 3.

Figure: (a) ogbn-arxiv
Refer to caption: extracted/5987727/curve_1_1.png

In both experiments, YOSO consistently outperformed the baselines in terms of convergence speed and stability. On the ogbn-arxiv dataset, YOSO reached a lower training loss more rapidly than GraphSAGE, GraphSAINT-EDGE, and FastGCN, with significantly fewer oscillations, indicating a more stable and efficient training process. Similarly, on the ogbl-ppa dataset, YOSO demonstrated faster convergence and maintained a smoother training loss curve, while the baselines, especially FastGCN, exhibited more fluctuations. These results suggest that YOSO not only accelerates the convergence process but also ensures a more stable training path compared to existing sampling methods, highlighting its effectiveness in GNN training

### 6.4 Ablation Study

In this subsection, we explore how YOSO’s total training time and model accuracy vary with different sampling sizes $M$ and evaluate reconstruction effectiveness by comparing the $\mathbf{H}^{(L)}$ matrix generated without sampling to the $\tilde{\mathbf{H}}^{(L)}$ matrix produced by YOSO’s sampling-reconstruction process, with the differences visualized with heatmaps.

Varying sampling size $M$:
We examine how total training time (including breakdown) and model accuracy vary with $M$ values, specifically $M=\{64,128,256,1024,2048\}$, as shown in Figure 4. The results indicate that YOSO’s sampling time remains stable across different $M$, ranging from 107.94 to 111.53 seconds on ogbn-products and 143.56 to 149.65 seconds on ogbl-citation2, showing minimal impact from $M$. In contrast, as $M$ decreases, computation time increases, reflecting more iterations needed for convergence (e.g., rising from 275.98s at $M=2048$ to 301.94s at $M=64$ on ogbn-products, with a similar trend on ogbl-citation2). Model accuracy improves with larger $M$, eventually stabilizing; it rises from 0.597 to 0.7873 on ogbn-products and from 0.312 to 0.8025 on ogbl-citation2. These findings highlight YOSO’s efficient sampling and improved accuracy and convergence with larger $M$.

Figure: Figure 4: Total training time (including its breakdown) and model accuracy for YOSO with different sampling sizes: (a) for the node classification learning task on the ogbn-products dataset, and (b) for the link prediction learning task on the ogbl-citation2 dataset.
Refer to caption: extracted/5987727/product_M_varies.png

Figure: Figure 5: Reconstruction effectiveness visualized via heatmap. Using the ogbn-products dataset, 10 nodes are randomly selected from the training set, and for each node, 10 embedding dimensions are randomly picked. The heatmap shows the absolute differences between original and reconstructed embeddings for these elements. $M$ is the size of the sampling set.
Refer to caption: x6.png

Reconstruction effectiveness: The heatmap in Figure 5 shows the reconstruction effectiveness for different sampling sizes $M$. Each $10\times 10$ block represents the absolute difference between reconstructed embeddings from our two-layer GNN sampling and those computed with all neighbors (without sampling). As $M$ increases, reconstruction accuracy improves, enhancing overall model accuracy. However, beyond a certain point, such as $M=512$ in Figure 5, further increases in $M$ offer diminishing returns in both reconstruction quality and model accuracy. This suggests there is an optimal $M$ that balances reconstruction quality and computational efficiency.

## 7 conclusion

In this paper, we introduce YOSO (You Only Sample Once), a novel algorithm aimed at significantly enhancing the efficiency of GNN training without sacrificing prediction accuracy. By leveraging a compressed sensing-based sampling and reconstruction framework, YOSO performs node sampling only once at the input layer, followed by a lossless reconstruction at the output layer during each training epoch. Our experimental results demonstrate that YOSO can achieve up to 75% reduction of existing state-of-the-art methods while achieving accuracy comparable to top-performing baselines.

Ethics Statement: In this paper, we present a technique grounded in compressed sensing that addresses the growing computational demands of GNN sampling schemes. Our approach significantly reduces sampling time and overall GNN training duration without compromising model accuracy, thereby enhancing the efficiency of graph neural network training. This improvement holds potential for a wide range of applications, such as recommendation systems and social network analysis, and bioinformatics. We believe that our method contributes positively to the advancement of machine learning research by promoting computational efficiency. Although we do not anticipate any immediate negative ethical implications or societal concerns from our approach, it’s important to acknowledge that machine learning technologies, including graph-based methods, have broader impacts. Therefore, responsible implementation is crucial to ensure that such technologies are applied in a manner that promotes fairness and beneficial societal outcomes.

## Appendix A Details about experiments

### A.1 Hardware and Software Configuration

We evaluate all baselines and our design on a Linux Desktop running Ubuntu 18.04.6 LTS, equipped with an NVIDIA GTX 1060Ti (6GB memory) using CUDA version 11.8 and PyTorch version 2.0.0. The system features a AMD Ryzen 5 5500 CPU with 64 GB DDR4 RAM, and the Python version used is 3.9.0.

### A.2 Datasets

Data splitting: We adopt strategies consistent with previous works (Hamilton et al., ; Hu et al., ). Specifically, for the Reddit dataset, we follow the data splitting used in GraphSage (Hamilton et al., ), and for the OGB series (ogbn and ogbl), we maintain the splitting described in (Hu et al., ).

The basic summary information of the datasets we use is provided in Table 2, and detailed descriptions are as follows:
ogbn-arxiv: This dataset is a directed citation network of Computer Science (CS) arXiv papers from the Microsoft Academic Graph (MAG) (Wang et al., ). Each node represents a paper, with directed edges indicating citations. The task is to classify unlabeled papers into primary categories using labeled papers and node features, which are derived by averaging word2vec embeddings (Mikolov et al., ) of paper titles and abstracts.
Reddit: Originally from GraphSage (Hamilton et al., ), this Reddit dataset is a post-to-post graph where each node represents a post, and edges indicate shared user comments. The task is to classify posts into communities using GloVe word vectors (Pennington et al., ) from post titles and comments, along with features such as post scores and comment counts.
ogbn-products: This undirected, unweighted graph represents an Amazon product co-purchasing network, where nodes are products and edges indicate frequent co-purchases. Node features are derived from bag-of-words features of product descriptions, reduced to 100 dimensions via Principal Component Analysis (Dunteman, ).
ogbl-ppa: This undirected, unweighted graph has nodes representing proteins from 58 species, with edges indicating biologically meaningful associations. Each node features a 58-dimensional one-hot vector for the protein’s species. The task is to predict new association edges, evaluated by ranking positive test edges over negative ones.
ogbl-citation2: This dataset is a directed graph representing a citation network among a subset of papers from Microsoft Academic Graph (MAG), similar to ogbn-arxiv. For each source paper, two references are randomly removed, and the task is to rank these missing references above 1,000 randomly selected negative references, which are sampled from all papers not cited by the source paper.

**Table 2: Statistics and metrics of the dataset**
| Dataset | #Node | #Edge | #Dim. | Metric |  |
| --- | --- | --- | --- | --- | --- |
| Node Property<br>Prediction | ogbn-arxiv | 169,343 | 1,166,243 | 128 | Accuracy |
| Reddit | 232,965 | 11,606,919 | 602 | Mirco-F1 |  |
| ogbn-products | 2,449,029 | 61,859,140 | 100 | Accuracy |  |
| Link Property<br>Prediction | ogbl-ppa | 576,289 | 30,326,273 | 128 | Hits@100 |
| ogbl-citation2 | 2,927,963 | 30,561,187 | 128 | MRR |  |

### A.3 Baselines and Implementation

**Table 3: Baselines and their public available source code link**
| Method | Available Link |
| --- | --- |
| GraphSage | https://github.com/williamleif/graphsage-simple |
| VR-GCN | https://github.com/THUDM/cogdl/tree/master/examples/VRGCN |
| FastGCN | https://github.com/gmancino/fastgcn-pytorch |
| AS-GCN | https://github.com/Gkunnan97/FastGCN_pytorch |
| LADIES | https://github.com/acbull/LADIES |
| Cluster-GCN | https://github.com/benedekrozemberczki/ClusterGCN |
| GraphSAINT | https://github.com/GraphSAINT/GraphSAINT |

Table 3 presents the baselines used in this paper along with their publicly available source code links. Since some baselines were not originally implemented in PyTorch, we standardized the framework for fair comparison. If a PyTorch version involved the original authors, we selected that source code (e.g., FastGCN (Chen et al., )). Otherwise, we chose the most popular implementation based on the number of stars. Notably, the repository linked for AS-GCN (Huang et al., ) in the table includes implementations of both FastGCN and AS-GCN, but we only used the AS-GCN version, while the FastGCN implementation was taken from the source listed in the table.
YOSO’s Implementation: The base code of YOSO(^1^11https://anonymous.4open.science/r/YOSO-B49B) is built on GCN (Kipf & Welling, ), with the link available at https://github.com/tkipf/pygcn. The sampling stage in YOSO occurs on the CPU and main memory since it involves calculations related to the entire feature matrix and the regularized Laplacian matrix. After sampling, the relevant data is migrated to GPU memory for computation. Throughout the training process, multiple data exchanges occur between main memory and GPU memory, such as in link prediction tasks where node embeddings need to be updated.
Modification: All baselines support updating node embeddings and performing node classification tasks. For node classification, if a baseline did not originally use the cross-entropy loss function, we adjusted it to adopt this loss function. For the link prediction task, the following loss function is applied:

$$ $\mathcal{L}=\frac{1}{N^{+}}\sum_{(i,j)\in E^{+}}\left(1-\frac{\mathbf{h}_{i}^{ (L)}\cdot\mathbf{h}_{j}^{(L)}}{\|\mathbf{h}_{i}^{(L)}\|\|\mathbf{h}_{j}^{(L)} \|}\right)+\frac{1}{N^{-}}\sum_{(i,j)\in E^{-}}\text{max}\left(0,\gamma-\left( 1-\frac{\mathbf{h}_{i}^{(L)}\cdot\mathbf{h}_{j}^{(L)}}{\|\mathbf{h}_{i}^{(L)} \|\|\mathbf{h}_{j}^{(L)}\|}\right)\right)$ $$

where $N^{+}$ and $N^{-}$ represent the number of positive and negative samples, respectively, and $E^{+}$ and $E^{-}$ denote the sets of positive and negative edges. The parameter $\gamma$ is a hyperparameter, set to 0.5 in this study. As the ogbl-ppa and ogbl-citation2 datasets provide corresponding negative edges by default, we used these pre-defined negative edges for our calculations.

### A.4 Hyper-parameter Setting

The hyperparameter settings for both YOSO and the baselines are provided in Table 4 and Table 5 for node classification and link prediction datasets, respectively. All experiments were conducted using a two-layer GCN with official configurations. When certain parameters were not clearly specified in some papers, we fine-tuned them for optimal accuracy. The recorded hyperparameters include the sampling size (per node/layer/subgraph), the optimizer, and the learning rate. For YOSO, the sampling size is denoted as $M$; for example, on the ogbl-ppa dataset (Table 5), $M=128$.

**Table 4: Node classification hyperparamter setting for baselines and YOSO on different datasets.**
|  | ogbn-arxiv | Reddit | ogbn-products |
| --- | --- | --- | --- |
| GraphSage | 25&10 / Adam / 0.7 | 25&10 / Adam / 0.01 | 50&20 / Adam / 0.01 |
| VR-GCN | 8 / Adam / 0.01 | 16 / Adam / 0.01 | 32 / Adam / 0.01 |
| FastGCN | 64 / Adam / 0.01 | 128 / Adam / 0.001 | 256 / Adam / 0.001 |
| AS-GCN | 128 / Adam / 0.001 | 512 / Adam / 0.01 | 1000 / Adam / 0.01 |
| LADIES | 64 / Adam / 0.001 | 128 / Adam / 0.001 | 256 / Adam / 0.001 |
| Cluster-GCN | - / Adam / 0.01 | - / Adam / 0.005 | - / Adam / 0.005 |
| GraphSAINT-EG | 300 / Adam / 0.01 | 600 / Adam / 0.01 | 4000 / Adam / 0.01 |
| GraphSAINT-RW | 4000 / Adam / 0.01 | 8000 / Adam / 0.01 | 10000 / Adam / 0.01 |
| YOSO | 128 / Adam / 0.01 | 256 / Adam / 0.01 | 512 / Adam / 0.01 |

**Table 5: Link prediction hyperparamter setting for baselines and YOSO on different datasets.**
|  | ogbl-ppa | ogbl-citation2 |
| --- | --- | --- |
| GraphSage | 25&10 / Adam / 0.7 | 50&20 / Adam / 0.01 |
| VR-GCN | 8 / Adam / 0.01 | 32 / Adam / 0.01 |
| FastGCN | 64 / Adam / 0.01 | 256 / Adam / 0.001 |
| AS-GCN | 128 / Adam / 0.001 | 1000 / Adam / 0.01 |
| LADIES | 64 / Adam / 0.001 | 256 / Adam / 0.001 |
| Cluster-GCN | - / Adam / 0.01 | - / Adam / 0.005 |
| GraphSAINT-EG | 300 / Adam / 0.01 | 4000 / Adam / 0.01 |
| GraphSAINT-RW | 4000 / Adam / 0.01 | 10000 / Adam / 0.01 |
| YOSO | 128 / Adam / 0.01 | 512 / Adam / 0.01 |

## Appendix B Computation and Proof

### B.1 Gradient Computation

#### B.1.1 Computation of ∇ Θ ℒ subscript ∇ Θ ℒ \nabla_{\Theta}\mathcal{L} ∇ start_POSTSUBSCRIPT roman_Θ end_POSTSUBSCRIPT caligraphic_L :

$\nabla_{\Theta}\mathcal{L}=\alpha\nabla_{\Theta}\mathcal{L}_{recon}+\beta
\nabla_{\Theta}\mathcal{L}_{GNN}^{\Theta}(\mathbf{Z})=\frac{\partial\mathcal{L
}_{recon}}{\partial\mathbf{Z}}\cdot\frac{\partial\mathbf{Z}}{\partial\Theta}+
\frac{\partial\mathcal{L}_{GNN}^{\Theta}(\mathbf{Z})}{\partial\mathbf{Z}}\cdot
\frac{\partial\mathbf{Z}}{\partial\Theta}$

- •
$\frac{\partial\mathcal{L}_{recon}}{\partial\mathbf{Z}}=(\mathbf{Z}-\mathbf{
\Phi}\mathbf{U}\hat{\mathbf{H}}^{(L)})$
- •
Consider the $g^{(L)}$ which is the gradient at the output layer, and we have $g^{(L)}=\frac{\partial\mathcal{L}_{recon}}{\partial\mathbf{Z}}\odot\sigma^{{}^
{\prime}}(\mathbf{S}^{(L)})$ where $\odot$ denotes element-wise multiplication, $\sigma^{{}^{\prime}}(\mathbf{S}^{(L)})$ is the derivation of the activation function at layer $L$ and $\mathbf{S}^{(L)}$ is the pre-activation input at layer $L$. Therefore, for $l=L,L-1,...,1$, we have $g^{(l-1)}=\nabla_{\mathbf{W}^{(L)}}\mathcal{L}_{recon}\odot\sigma^{{}^{\prime}
}(\mathbf{S}^{(l-1)})=(\mathbf{\Phi}\hat{\mathbf{A}}\mathbf{W}^{(l)})^{T}g^{(l
)}\odot\sigma^{{}^{\prime}}(\mathbf{S}^{(l-1)})$. By iteratively executing this process, we can obtain $\frac{\partial\mathbf{Z}}{\partial\Theta}$
- •
$\frac{\partial\mathcal{L}_{GNN}^{\Theta}}{\partial\mathbf{Z}}$ depends on the specific loss function used.

#### B.1.2 Computation of ∇ 𝐔 ℒ subscript ∇ 𝐔 ℒ \nabla_{\mathbf{U}}\mathcal{L} ∇ start_POSTSUBSCRIPT bold_U end_POSTSUBSCRIPT caligraphic_L

$\nabla_{\mathbf{U}}\mathcal{L}=\alpha\nabla_{\mathbf{U}}\mathcal{L}_{recon}+
\beta\nabla_{\mathbf{U}}\mathcal{L}_{GNN}^{\Theta}(\mathbf{Z})=\alpha\nabla_{
\mathbf{U}}\mathcal{L}_{recon}+\beta(\frac{\partial\mathcal{L}_{GNN}^{\Theta}(
\mathbf{Z})}{\partial\mathbf{Z}}\cdot\frac{\partial\mathbf{Z}}{\partial\mathbf
{U}})$

- •
$\nabla_{\mathbf{U}}\mathcal{L}_{recon}=-\mathbf{\Phi}^{T}(\mathbf{Z}-\mathbf{
\Phi}\mathbf{U}\hat{\mathbf{H}}^{(L)})(\hat{\mathbf{H}}^{(L)})^{T}$
- •
As in Section B.1.1, $\frac{\partial\mathcal{L}_{GNN}^{\Theta}}{\partial\mathbf{Z}}$ depends on specific loss function and easy to compute.
- •
For $\frac{\partial\mathbf{Z}}{\partial\mathbf{U}}$, it need to be computed recursively. Since $\mathbf{T}^{(0)}=\mathbf{\Phi}\mathbf{U}\hat{\mathbf{X}}$, $\frac{\partial\mathbf{T}^{(0)}}{\partial\mathbf{U}}=\mathbf{\Phi}\hat{\mathbf{
X}}$. The gradient propagates from $\mathbf{Z}$ back to $\mathbf{U}$: $\nabla_{\mathbf{U}}\mathcal{L}_{GNN}^{\Theta}(\mathbf{Z})=(\frac{\partial
\mathcal{L}_{GNN}^{\Theta}(\mathbf{Z})}{\partial\mathbf{Z}}\cdot\frac{\partial
\mathbf{Z}}{\partial\mathbf{T}^{(L-1)}}\cdots\frac{\partial\mathbf{T}^{(1)}}{
\partial\mathbf{T}^{(0)}}\cdot\frac{\partial\mathbf{T}^{(0)}}{\partial\mathbf{
U}})$. As we know that $\mathbf{T}^{(l)}=\sigma(\mathbf{S}^{(l)})$ and $\mathbf{S}^{(l)}=\mathbf{\Phi}\hat{\mathbf{A}}\mathbf{W}^{(l)}\mathbf{T}^{(l-1)}$, therefore $\frac{\partial\mathbf{T}^{(l-1)}}{\partial\mathbf{T}^{(l-1)}}=(\mathbf{\Phi}
\hat{\mathbf{A}}\mathbf{W}^{(l)})^{T}\text{diag}(\sigma^{{}^{\prime}}(\mathbf{
S}^{(l)}))$

#### B.1.3 Computation of ∇ 𝐇 ^ ( L ) ℒ subscript ∇ superscript ^ 𝐇 𝐿 ℒ \nabla_{\hat{\mathbf{H}}^{(L)}}\mathcal{L} ∇ start_POSTSUBSCRIPT over^ start_ARG bold_H end_ARG start_POSTSUPERSCRIPT ( italic_L ) end_POSTSUPERSCRIPT end_POSTSUBSCRIPT caligraphic_L

$\nabla_{\hat{\mathbf{H}}^{(L)}}\mathcal{L}=\alpha\nabla_{\hat{\mathbf{H}}^{(L)
}}\mathcal{L}_{recon}=-\mathbf{U}^{T}\mathbf{\Phi}^{T}(\mathbf{Z}-\mathbf{\Phi
}\mathbf{U}\hat{\mathbf{H}}^{(L)})+\lambda\partial\|\hat{\mathbf{H}}^{(L)}\|_{
2,1}$ where $\partial\|\hat{\mathbf{H}}^{(L)}\|_{2,1}$ is the subgradient of the $l_{2,1}$ norm and computed as $(\partial\|\hat{\mathbf{H}}^{(L)}\|_{2,1})_{i}=\frac{\hat{\mathbf{H}}^{(L)}_{i
,:}}{\|\hat{\mathbf{H}}^{(L)}_{i,:}\|_{2}}$ if and only if $\hat{\mathbf{H}}^{(L)}_{i,:}\neq 0$, otherwise, $(\partial\|\hat{\mathbf{H}}^{(L)}\|_{2,1})_{i}=0$

### B.2 Full Rank of 𝚽 𝚽 \mathbf{\Phi} bold_Φ

Theorem 1: Let $\hat{\mathbf{S}}\in\mathbb{R}^{M\times N}$ be a binary sampling matrix derived from the graph’s structure, where each entry $\hat{\mathbf{S}}_{i,j}\in\{0,1\}$ and each row has at least one non-zero entry. Let $\mathbf{\Sigma}\in\mathbb{R}^{M\times N}$ be a random matrix with entries drawn independently from a continuous probability distribution. Define $\mathbf{\Phi}=\hat{\mathbf{S}}\otimes\mathbf{\Sigma}$, where $\otimes$ denotes element-wise multiplication. Then, with probability $1$, the matrix $\mathbf{\Phi}$ has full row rank $M$.
Proof: First, we know that the structure of $\mathbf{\Phi}$ satisfies the following conditions:

- •
Each entry of $\mathbf{\Phi}$ is given by $\mathbf{\Phi}_{i,j}=\hat{\mathbf{S}}_{i,j}\cdot\mathbf{\Sigma}_{i,j}$.
- •
The $i$-th row of $\mathbf{\Phi}$ is $\mathbf{\Phi}_{i,:}=\hat{\mathbf{S}}_{i,:}\otimes\mathbf{\Sigma}_{i,:}$.
- •
Non-zero entries in $\mathbf{\Phi}_{i,:}$ correspond to positions where $\hat{\mathbf{S}}_{i,j}=1$.

Assume there exist scalars $c_{1},c_{2},\dots,c_{M}$, not all zero, such that $\sum_{i=1}^{M}c_{i}\mathbf{\Phi}_{i,:}=\mathbf{0}$. This implies that for each $j=1,\dots,N$, we have $\sum_{i=1}^{M}c_{i}\hat{\mathbf{S}}_{i,j}\mathbf{\Sigma}_{i,j}=0$. Let $I_{j}=\{i\mid\hat{\mathbf{S}}_{i,j}=1\}$; then $\sum_{i\in I_{j}}c_{i}\mathbf{\Sigma}_{i,j}=0$.

Since the $\mathbf{\Sigma}_{i,j}$ values are independently drawn from continuous distributions, the probability that this equation holds for any non-zero set of $\{c_{i}\}$ is zero unless all $c_{i}$ in $I_{j}$ are zero. Therefore, for the equation to be valid, $c_{i}=0$ for all $i$ where $\hat{\mathbf{S}}_{i,j}=1$.

As each row $i$ contains at least one entry with $\hat{\mathbf{S}}_{i,j}=1$, it follows that $c_{i}=0$ for all $i$. This contradicts the assumption that not all $c_{i}$ are zero. Hence, the only solution is $c_{i}=0$ for all $i$, indicating that the rows of $\mathbf{\Phi}$ are linearly independent with probability 1. Thus, $\text{rank}(\mathbf{\Phi})=M$ with probability 1.

### B.3 Sampling Matrix 𝚽 𝚽 \mathbf{\Phi} bold_Φ and RIP

Theorem 2: Let $\hat{\mathbf{S}}\in\mathbb{R}^{M\times N}$ be a selection matrix derived from the graph’s structure, where each entry $\hat{\mathbf{S}}_{i,j}\in\{0,1\}$ indicates whether node $j$ is included in the $i$-th measurement. Let $\mathbf{\Sigma}\in\mathbb{R}^{M\times N}$ be a matrix whose entries $\mathbf{\Sigma}_{i,j}$ are independent sub-Gaussian random variables with mean zero and variance $\frac{1}{g(j)}$, where $g(j)>0$. Define the sampling matrix $\mathbf{\Phi}=\hat{\mathbf{S}}\otimes\mathbf{\Sigma}$, where $\otimes$ denotes element-wise multiplication. Then, for any $0<\delta_{k}<1$, there exists a constant $c>0$ such that if $M\geq c\cdot k\log\left(\frac{N}{k}\right)$, then with probability at least $1-e^{-cM}$, the matrix $\mathbf{\Phi}\mathbf{U}$ satisfies the Restricted Isometry Property (RIP) of order $k$ with constant $\delta_{k}$; that is, for all $\hat{\mathbf{H}}\in\mathbb{R}^{N\times d}$ with $\|\hat{\mathbf{H}}\|_{0,\text{row}}\leq k$,

$$ $(1-\delta_{k})\|\hat{\mathbf{H}}\|_{F}^{2}\leq\|\mathbf{\Phi}\mathbf{U}\hat{ \mathbf{H}}\|_{F}^{2}\leq(1+\delta_{k})\|\hat{\mathbf{H}}\|_{F}^{2}.$ $$

Proof: To demonstrate that $\mathbf{\Phi}\mathbf{U}$ satisfies the Restricted Isometry Property (RIP) of order $k$ with high probability, we consider $\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}=(\hat{\mathbf{S}}\otimes\mathbf{\Sigma
})\mathbf{U}\hat{\mathbf{H}}$. For each row $i$ and column $r$, the entry $(\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}})_{i,r}$ can be expressed as $\sum_{j=1}^{N}\hat{\mathbf{S}}_{i,j}\mathbf{\Sigma}_{i,j}(\mathbf{U}\hat{
\mathbf{H}})_{j,r}$. This sum only involves terms where $\hat{\mathbf{S}}_{i,j}=1$. Therefore, $(\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}})_{i,r}=\sum_{j\in\mathcal{S}_{i}}
\mathbf{\Sigma}_{i,j}(\mathbf{U}\hat{\mathbf{H}})_{j,r}$, where $\mathcal{S}_{i}=\{j\mid\hat{\mathbf{S}}_{i,j}=1\}$.

The variables $\mathbf{\Sigma}_{i,j}$ are independent sub-Gaussian random variables with mean zero and variance $\frac{1}{g(j)}$. Therefore, the expectation of $\|\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}$ can be computed as follows:

$$ $\mathbb{E}\left[\|\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}\right]= \sum_{i=1}^{M}\sum_{r=1}^{d}\mathbb{E}\left[\left(\sum_{j\in\mathcal{S}_{i}} \mathbf{\Sigma}_{i,j}(\mathbf{U}\hat{\mathbf{H}})_{j,r}\right)^{2}\right]$ $$

Expanding this and leveraging the independence of $\mathbf{\Sigma}_{i,j}$, we have:

$$ $\mathbb{E}\left[\left(\sum_{j\in\mathcal{S}_{i}}\mathbf{\Sigma}_{i,j}(\mathbf{ U}\hat{\mathbf{H}})_{j,r}\right)^{2}\right]=\sum_{j\in\mathcal{S}_{i}}\mathbb{ E}\left[\mathbf{\Sigma}_{i,j}^{2}\right]\left((\mathbf{U}\hat{\mathbf{H}})_{j, r}\right)^{2}$ $$

Since $\mathbb{E}[\mathbf{\Sigma}_{i,j}^{2}]=\frac{1}{g(j)}$, the expectation simplifies to:

$$ $\mathbb{E}\left[\|\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}\right]= \sum_{i=1}^{M}\sum_{j\in\mathcal{S}_{i}}\frac{1}{g(j)}\sum_{r=1}^{d}\left(( \mathbf{U}\hat{\mathbf{H}})_{j,r}\right)^{2}$ $$

If we assume $p(j)=\frac{g(j)}{G}$, where $G=\sum_{j=1}^{N}g(j)$ serves as a normalization factor, the expected measurement count for each node $j$ is $Mp(j)=M\frac{g(j)}{G}$. Thus:

$$ $\mathbb{E}\left[\|\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}\right]= \sum_{j=1}^{N}M\frac{g(j)}{G}\cdot\frac{1}{g(j)}\|(\mathbf{U}\hat{\mathbf{H}}) _{j,:}\|_{2}^{2}=\frac{M}{G}\|\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}$ $$

By setting $G=M$, we have:

$$ $\mathbb{E}\left[\|\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}\right]=\| \hat{\mathbf{H}}\|_{F}^{2}$ $$

Now, define $Z_{i,r}=\sum_{j\in\mathcal{S}_{i}}\mathbf{\Sigma}_{i,j}(\mathbf{U}\hat{\mathbf
{H}})_{j,r}$, which are sub-Gaussian random variables. Applying Bernstein’s inequality, we obtain:

$$ $\mathbb{P}\left(\left|\|\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}\|_{F}^{2}-\| \hat{\mathbf{H}}\|_{F}^{2}\right|\geq\delta_{k}\|\hat{\mathbf{H}}\|_{F}^{2} \right)\leq 2\exp\left(-c\cdot\frac{\delta_{k}^{2}\|\hat{\mathbf{H}}\|_{F}^{4} }{\sum_{i,r}\sigma_{i,r}^{2}}\right)$ $$

where $\sigma_{i,r}^{2}=\sum_{j\in\mathcal{S}_{i}}\frac{1}{g(j)}\left((\mathbf{U}\hat
{\mathbf{H}})_{j,r}\right)^{2}$. By bounding the total variance, we conclude that the probability of RIP failing is very low. This confirms that $\mathbf{\Phi}\mathbf{U}$ satisfies the RIP for all sparse $\hat{\mathbf{H}}$ with $\|\hat{\mathbf{H}}\|_{0,\text{row}}\leq k$ with high probability.

### B.4 Error Bound

Theorem 3: Let $\mathbf{H}^{(L)}$ be the output embeddings obtained by the standard GNN computation with full reconstruction at each layer as per Equation (5). Let $\tilde{\mathbf{H}}^{(L)}$ be the output embeddings obtained by Algorithm 1, which performs sampling once at the input layer and reconstructs only at the output layer. Assume that the activation function $\sigma$ is Lipschitz continuous with Lipschitz constant $L_{\sigma}$, and the sampling matrix $\mathbf{\Phi}\mathbf{U}$ satisfies the Restricted Isometry Property (RIP) of order $k$ with constant $\delta_{k}$ (i.e., $0<\delta_{k}<1$). Then, the error between $\tilde{\mathbf{H}}^{(L)}$ and $\mathbf{H}^{(L)}$ can be bounded as:

$$ $\left\|\tilde{\mathbf{H}}^{(L)}-\mathbf{H}^{(L)}\right\|_{F}\leq\left(\frac{L_ {\sigma}}{1-\delta_{k}}\right)^{L}\left\|\mathbf{E}\right\|_{F},$ $$

where $\mathbf{E}=\mathbf{Z}-\mathbf{\Phi}\mathbf{U}\hat{\mathbf{H}}^{(L)}$ is the reconstruction error at the output layer, and $L$ is the number of layers in the GNN.
Proof: We aim to bound the error $\left\|\tilde{\mathbf{H}}^{(L)}-\mathbf{H}^{(L)}\right\|_{F}$ between the output embeddings of the standard GNN computation and those obtained by Algorithm 1.

Assume the activation function $\sigma$ is Lipschitz continuous with a constant $L_{\sigma}$, such that

$$ $\left\|\sigma(\mathbf{X})-\sigma(\mathbf{Y})\right\|_{F}\leq L_{\sigma}\left\| \mathbf{X}-\mathbf{Y}\right\|_{F}\quad\forall\mathbf{X},\mathbf{Y}.$ $$

Further, let the sampling matrix $\mathbf{\Phi}\mathbf{U}$ satisfy the RIP of order $k$ with constant $\delta_{k}$, meaning

$$ $(1-\delta_{k})\left\|\hat{\mathbf{H}}\right\|_{F}^{2}\leq\left\|\mathbf{\Phi} \mathbf{U}\hat{\mathbf{H}}\right\|_{F}^{2}\leq(1+\delta_{k})\left\|\hat{ \mathbf{H}}\right\|_{F}^{2},$ $$

for all $\hat{\mathbf{H}}$ with $\left\|\hat{\mathbf{H}}\right\|_{0,\text{row}}\leq k$. We also have $\mathbf{H}^{(l)}=\mathbf{U}\hat{\mathbf{H}}^{(l)}$, where $\hat{\mathbf{H}}^{(l)}$ has at most $k$ non-zero rows.

We will prove by induction on $l=1,2,\dots,L$ that

$$ $\left\|\tilde{\mathbf{H}}^{(l)}-\mathbf{H}^{(l)}\right\|_{F}\leq\left(\frac{L_ {\sigma}}{1-\delta_{k}}\right)^{l}\left\|\tilde{\mathbf{H}}^{(0)}-\mathbf{H}^{ (0)}\right\|_{F}.$ $$

For the base case $l=0$, at the input layer, we have $\tilde{\mathbf{H}}^{(0)}=\mathbf{U}\hat{\mathbf{X}}$ and $\mathbf{H}^{(0)}=\mathbf{X}$. The initial error $\left\|\tilde{\mathbf{H}}^{(0)}-\mathbf{H}^{(0)}\right\|_{F}$ is assumed.

Assume that for some $l\geq 0$,

$$ $\left\|\tilde{\mathbf{H}}^{(l)}-\mathbf{H}^{(l)}\right\|_{F}\leq\left(\frac{L_ {\sigma}}{1-\delta_{k}}\right)^{l}\left\|\tilde{\mathbf{H}}^{(0)}-\mathbf{H}^{ (0)}\right\|_{F}.$ $$

We aim to show that

$$ $\left\|\tilde{\mathbf{H}}^{(l+1)}-\mathbf{H}^{(l+1)}\right\|_{F}\leq\left( \frac{L_{\sigma}}{1-\delta_{k}}\right)^{l+1}\left\|\tilde{\mathbf{H}}^{(0)}- \mathbf{H}^{(0)}\right\|_{F}.$ $$

For Algorithm 1, $\tilde{\mathbf{T}}^{(l)}=\sigma\left(\mathbf{\Phi}\hat{\mathbf{A}}\mathbf{W}^{
(l+1)}\tilde{\mathbf{T}}^{(l-1)}\right)$. At the output layer $l=L$, we perform reconstruction:

$$ $\tilde{\mathbf{H}}^{(L)}=\mathbf{U}\hat{\mathbf{H}}^{(L)},$ $$

where $\hat{\mathbf{H}}^{(L)}$ is obtained by solving

$$ $\min_{\hat{\mathbf{H}}^{(L)}}\ \frac{1}{2}\left\|\mathbf{Z}-\mathbf{\Phi} \mathbf{U}\hat{\mathbf{H}}^{(L)}\right\|_{F}^{2}+\lambda\left\|\hat{\mathbf{H} }^{(L)}\right\|_{2,1},$ $$

with $\mathbf{Z}=\tilde{\mathbf{T}}^{(L)}$. Due to the optimization and the RIP condition, we have

$$ $\left\|\hat{\mathbf{H}}^{(L)}-\hat{\mathbf{H}}^{(L)}_{\text{true}}\right\|_{F} \leq C_{\text{rec}}\left\|\mathbf{E}\right\|_{F},$ $$

where $\hat{\mathbf{H}}^{(L)}_{\text{true}}$ is the true sparse representation of $\mathbf{H}^{(L)}$, and $C_{\text{rec}}=\frac{2\delta_{k}}{1-\delta_{k}}$. Since $\mathbf{U}$ is orthonormal,

$$ $\left\|\tilde{\mathbf{H}}^{(L)}-\mathbf{H}^{(L)}\right\|_{F}=\left\|\hat{ \mathbf{H}}^{(L)}-\hat{\mathbf{H}}^{(L)}_{\text{true}}\right\|_{F},$ $$

implying

$$ $\left\|\tilde{\mathbf{H}}^{(L)}-\mathbf{H}^{(L)}\right\|_{F}\leq\frac{2\delta_ {k}}{1-\delta_{k}}\left\|\mathbf{E}\right\|_{F}.$ $$

Given the Lipschitz continuity of $\sigma$, the error accumulates multiplicatively through $L$ layers:

$$ $\left\|\tilde{\mathbf{H}}^{(L)}-\mathbf{H}^{(L)}\right\|_{F}\leq\left(\frac{L_ {\sigma}}{1-\delta_{k}}\right)^{L}\left\|\tilde{\mathbf{H}}^{(0)}-\mathbf{H}^{ (0)}\right\|_{F}.$ $$

If the initial error $\left\|\tilde{\mathbf{H}}^{(0)}-\mathbf{H}^{(0)}\right\|_{F}=0$, the primary source of error is from the reconstruction at the output layer, yielding

$$ $\left\|\tilde{\mathbf{H}}^{(L)}-\mathbf{H}^{(L)}\right\|_{F}\leq\left(\frac{L_ {\sigma}}{1-\delta_{k}}\right)^{L}\left\|\mathbf{E}\right\|_{F}.$ $$

## Peer Reviews

confidence: 5
contribution: 3
initial_rating: 3
presentation: 3
questions: + Does the (sparse reconstruction) optimization problem have to be solved at every layer, for every new graph? Or is the assumption that the graph is always fixed (between training and inference)?

+ If the parameters of the model significantly change, isn't the re-sampling necessary? Why not sample with some periodicity (as model drifts), instead of sample only once at-start?

+ Is it at-all possible to think about the GNN in the inductive setting in your framework?
soundness: 3
strengths: * GNNs are indeed popular architectures, solving problems across many domains. Even though there are many methods are proposes to speed-up training of GNNs, making things even faster, should realize some advantages.

* Their sampling method takes less time than competition, yet reaches SOTA metrics.

* Their sampling is computed only once, whereas competing methods draw subgraph samples for every input example every epoch. Since sampling time is significant (often the dominant), then sampling only once significantly saves resources.
summary: Paper is interested in sampling graphs, for the purpose of training graph neural networks, in cases where the input graph is large. Rather than taking gradient-steps using entire graphs, sampling-based training of GNNs can take gradient steps using subgraphs (more compute efficient).

Paper shows while indeed existing graph-sampling methods *can* scale learning onto larger graphs, *however*, they significant time in sampling itself (i.e., data prep) rather than on gradient calculation (i.e., training).
weaknesses: # Related Work

* Paper does not address non-sampling-based methods for speeding-up GNNs. At least, they should be mentioned in related work (IMO). One missed family is historical embeddings (e.g., https://arxiv.org/abs/1710.10568, https://arxiv.org/abs/2106.05609, https://arxiv.org/abs/2305.12322), another family is some "linearization" of models (a.k.a "decoupled" GNNs), e.g., https://arxiv.org/pdf/1902.07153, https://arxiv.org/abs/2004.11198, https://arxiv.org/abs/2111.06312)

* Other methods also sample-once. E.g., ClusterGCN.

# Sampling is inherently parallelizable
 
Paper is founded upon a statement that sampling takes more time than the actual training step. I believe the statement is correct (because the graph may be large and does not fit in memory). However, sampling is trivially distributable -- in fact, many papers, such as, PinSAGE, TFGNN, ..., propose and implement distributed or multi-threaded sampling. I think this could also be mentioned in the related work.


# Possible Writing Improvements

* First line of Intro: "analyzing" should probably become "modeling"
* The abstract says "lossless" compression but line 90 says "nearly lossless". It is better to be consistent.
* The term "embedding matrix" is used a few times in the first 2 pages, without clear definition of what it is -- is it the node input features? is it the hidden representations (between layers)? Is it the output of the GNN? I feel that it is all of those [but only after reading a couple of more pages].


# Missing analysis

The paper talks about bias and variance. This is usually done **either** in the forward pass -- let $\widetilde{z}$ be the latent values obtain via the graph sample and $\mathbf{z}$ be the latent values when full graph is used, the zero-biased analysis should show $\mathbb{E}[\widetilde{z}] = \mathbf{z}$; or in the backward pass e.g. $\mathbb{E}[\frac{\partial loss(\widetilde{z})}{\partial \theta}] = \frac{\partial loss(\mathbf{z})}{\partial \theta}$. This paper does neither of those, nor it shows bounds on the variance of these quantities. I suggest the authors remove bias&variance arguments or else show some analysis to back the arguments.

# Math inaccuracies
* Adjacency matrix is square (N x N). How is it multiplied on the right by the model parameter matrix? Are GNN parameters a function of the input graph? This happens in Eq6 and Lines 5&7 in Alg1.

* Line 1 in Alg1 does not explain the initialization process. Crucially, is $U$ initialized to be a(n orthonormal) basis?

* line 161 "exists" must be a function of k and/or the rank of $\widehat{H}$ -- for example, what if I choose $k=0$, then there is no orthonormal basis $U$ that can recover $H$. In fact, $rank(\widehat{H})$ must be at least equal the rank of $H$ for recovery to be possible.

* The sampling matrix (described below Eq1) can be better described. Can you explain its structure? or its motive? Even a google search on "compressed sensing sampling matrix" does not pull-up the answer immediately. What is its rank?

* Most-crucially, Eq2, which the rest of the work is founded upon, reads wrong/incomplete. Most-likely, *there are trivial typos, however, I will not do deductions and I expect all authors to ensure correctness of their **main** equation*. $\widetilde{H}$ reads as a scalar (the minimum norm).

* Line 339: eigenvalues do not correspond to nodes. They correspond to eigenvectors. For instance, one can do a low-rank representation of the adjacency matrix (effectively using few eigenvalues).

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Please see the weakness part.
soundness: 2
strengths: 1.	The paper is well-organized and easy to follow. It provides a thorough background on GNNs and compressed sensing, equipping readers with the necessary foundational knowledge. Additionally, the authors offer a detailed discussion on the use of compressed sensing for sampling and the proposed YOSO algorithm.
2.	The compressed sensing-based sampling method is both novel and effective in reducing computational overhead associated with sampling.
3.	The authors present experiments demonstrating YOSO's efficiency on large-scale graph datasets across multiple tasks, supporting the practical impact of their approach.
summary: This paper introduces YOSO, a novel compressed sensing-based sampling approach for Graph Neural Networks (GNNs) that samples nodes only once per training. By performing a one-time sampling at the input layer and a lossless reconstruction at the output layer, YOSO aims to minimize the computational burden associated with repeated sampling in GNN training. The proposed method is evaluated on node classification and link prediction tasks, demonstrating significant training time reduction (up to 75%) while maintaining comparable accuracy to several baselines.
weaknesses: 1.	The motivation for using compressed sensing-based sampling is unclear. Graph distillation or condensation methods could also be feasible for generating smaller graph datasets while preserving data distribution. What specific advantages does compressed sensing sampling offer over graph condensation?
2.	Compressed sensing relies on the Restricted Isometry Property (RIP) for accurate reconstruction. In the context of graph neural networks, it is unclear if RIP holds for the input feature matrix, representation matrix, and output matrix. Preliminary experiments would be beneficial to validate RIP applicability within GNNs.
3.	This paper proposes using an unknown $U^{(l)}$ and an universal $\Phi$ to enhance the efficiency. However, a key question remains: does the universal $\Phi$ really exist? More justification is needed. Additionally, what’s the disadvantage compared with layer-wise $\Phi$? Is the accuracy loss significant? 
4.	For Line 257, I am confused on the forward propagation. Suppose the shape of $\Phi$ is M*N, the normalized Laplacian matrix is N*N, trainable parameters W is d*d, $T^{(l-1)}$ is M*d. It seems that the forward equation is incorrect with an unmatched shape. Is it like $\Phi A \Phi^\top$ to transform the adjacency matrix?
5.	The limitation discussion of YOSO is missing. (1) is YOSO robust over the initialization of sampling matrix $\Phi$? What’s the performance on a random initialized measurement matrix? Why do the authors design a handcraft sampling matrix? (2) YOSO requires more memory during GNN training? What’s the memory consumption of YOSO compared with other methods? 

I am open to revising my score if the authors can address my concerns.

---------------------------------------Post Rebuttal----------------------------------

After reading the response, I increase the score to 6.

---

confidence: 3
contribution: 4
initial_rating: 6
presentation: 2
questions: 1. According to the article, YOSO operates on sparse data rather than the original dataset. Could this approach lead to potential information loss?
2. If there is data bias, how will it handle the same?
3. As the framework is tested for just two tasks, will it be generic for all type of graph related downstream tasks?
soundness: 3
strengths: 1. The paper introduces a fairly novel method to train the GNNs, which is a significant as per the details provided.
2. The mathematical representation provided is fair enough to make the readers understand.
3. The results and the studies are significant to show the efectiveness of the framework. As per the problem, large-scale datasets are considered for the study, which is also appreciable.
summary: This paper presents YOSO (You-Only-Sample-Once), an efficient algorithm designed to streamline graph neural network (GNN) training while maintaining predictive accuracy in downstream tasks. YOSO addresses the high computational overhead in existing sampling methods by introducing a novel compressed sensing-based framework. In this framework, nodes are sampled once at the input layer and then losslessly reconstructed at the output layer each epoch, eliminating costly operations like orthonormal basis calculations and ensuring high-probability accuracy retention comparable to full-node sampling. Experimental evaluations demonstrate YOSO’s effectiveness in reducing GNN training time by approximately 75% on tasks such as node classification and link prediction, achieving accuracy levels similar to leading baselines. This approach positions YOSO as a resource-efficient alternative with potential for significant performance benefits in large-scale GNN training.
weaknesses: 1. The work is evaluated for only two tasks, i.e., Node classification and link prediction, in which the link prediction results are also just comparable to the baselines.
2. The performance metrics considered is accuracy and loss only. There should me more like F1-score as there may be data bias in case of such datasets.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: Please refer to the points I mentioned in the weakness part.
soundness: 2
strengths: 1. The concept of "sampling only once" presented in this paper is really interesting and novel, offering a fresh perspective on GNN sampling.
2. The author provides theoretical support to guarantee the effectiveness of YOSO, which is very convincing.
3. The author provides code to ensure the reproducibility of the paper's results.
summary: This paper introduces YOSO, a novel algorithm for efficient GNN training using a compressed sensing approach. GNNs are crucial for analyzing structured data, but existing sampling methods often introduce significant computational overhead. YOSO addresses this by sampling nodes only once at the input layer and using a nearly lossless reconstruction at the output layer for each epoch. This method reduces training time by about 75% compared to state-of-the-art techniques while maintaining high accuracy.
weaknesses: 1. According to Formula 1, the generation of $\mathbf{T}$ is closely related to the feature matrix $\mathbf{H}$. However, in YOSO, $\mathbf{U}$ is randomly initialized, and the generation of $\mathbf{\Phi}$ is only related to the graph structure. Does this approach overlook too much feature information?
2. The methodology section of the paper lacks a description of critical steps. The resulting representations $\mathbf{Z}$  do not have the scale as the input matrix $\mathbf{X}$. So, how should the nodes in $\mathbf{X}$ correspond to the representations in $\mathbf{Z}$? How can we obtain the representations for all nodes in the entire graph? What nodes' loss does $L_{GNN}(\mathbf{Z})$ calculate? Does this mean that it is necessary to sample nodes from the training set?
3. Some important related works are missing. The paper only discusses related work on graph sampling, while relevant works on Compressed Sensing are not included. What are the objectives and common practices of this technique? Why can it be applied to the GNN sampling problem? These questions need further discussion.
4. The mathematical notation in the paper is confusing, making it difficult to understand. For example, in line 155, $\hat{\mathbf{H}}$ should be $\hat{\mathbf{H}}^{(l)}$; in line 159, $\mathbf{U}^{\top}$ should be ${\mathbf{U}^{(l)}}^\top$; and in Equation 2, $\min$ should be $\arg \min$.

Meta-review: This paper proposes a new sampling algorithm for efficient GNN training using compressed sensing. Although there has been extensive work in this area, the proposed sampling method is different from previous ones and has the benefit of sampling only once. Although this is an important area and the proposed method seems promising, there’s plenty of clarification issues in the paper, and as reviewers pointed out in the final discussion, even the revision submitted by the authors contain many errors in key formulas and symbols. Therefore, we think the current paper does not meet the ICLR standard.

Decision: Reject