Efficient Learning of Linear Graph Neural Networks via Node Subsampling
Authors: Seiyun Shin
Ilan Shomorony
Han Zhao
Venue: NeurIPS 2023 poster (2023)
Keywords: Graph neural networks
Random sampling
Regression

Abstract: Graph Neural Networks (GNNs) are a powerful class of machine learning models with applications in recommender systems, drug discovery, social network analysis, and computer vision. One challenge with their implementation is that GNNs often take large-scale graphs as inputs, which imposes significant computational/storage costs in the training and testing phases. In particular, the message passing operations of a GNN require multiplication of the graph adjacency matrix $A \in \mathbb{R}^{n \times n}$ and the data matrix $X \in \mathbb{R}^{n \times d}$, and the $O(n^2 d)$ time complexity can be prohibitive for large $n$. Thus, a natural question is whether it is possible to perform the GNN operations in (quasi-)linear time by avoiding the full computation of $A X$. To study this question, we consider the setting of a regression task on a two-layer Linear Graph Convolutional Network (GCN). We develop an efficient training algorithm based on (1) performing node subsampling, (2) estimating the leverage scores of $A X$ based on the subsampled graph, and (3) performing leverage score sampling on $A X$. We show that our proposed scheme learns the regression model observing only $O(nd\epsilon^{-2}\log n)$ entries of $A$ in time $O(nd^2 \epsilon^{-2}\log n)$, with the guarantee that the learned weights deviate by at most $\epsilon$ under the $\ell_2$ norm from the model learned using the entire adjacency matrix $A$. We present empirical results for regression problems on real-world graphs and show that our algorithm significantly outperforms other baseline sampling strategies that exploit the same number of observations.

Efficient Learning of Linear Graph Neural Networks via Node Subsampling
Seiyun Shin', Ilan Shomorony's and Han Zhao? Department of ECE and "Department of CS University of Illinois at Urbana- Chaamaa-aampiinn IL {seiyuns2, ilans, hanzhao] @i @eeellnniiiiiie edu
Abstract
Graph Neural Networks (GNNs) are powerful class of machine learning models with applications in recommender systems, drug discovery, social network analysis, and computer vision. One challenge with their implementation is that GNNs often take large-scale graphs as inputs, which imposes significant computational/storagg costs in the training and testing phases. particular, the message passing opera- tions of GNN require multiplication of the graph adjacency matrix A E R"*= and the data matrix X E Raad and the (nndd) time complexity can be prohibitive for large n. Thus, natural question is whether it is possible to perform the GNN operations in (quasi-) linner time by avoiding the full computation of AX. To study consider the setting of regression task on two-layer Linear Graph Convolutional Network (GCN). We develop an efficient training algorithm based on (1) performing node subsampling, (2) estimating the leverage scores of AX based on the subsampled graph, and (3) performing leverage score sampling AX. We show that our proposed scheme learns the regression model observing only (nds 2 log n) entries of in time O log n). with the guarantce that the learned weights deviate by at most under the lg norm from the model learned using the entire adjacency matrix A. We present empirical results for regression problems real-world graphs and show that algorithm significantly outperforms other baseline sampling strategies that exploit the same number observations.
this question,
Introduction
Graph Neural Networks (GNNs) have gained popularity powerful machine learning method for graph-structured data. By learning rich representations of graph data, GNNs can solve variety of prediction tasks on graphs (177]. GNNs have delivered impressive results across many different areas, including social network analysis 8 bioinformatics M--M. and recommendation systems [12].
Given their remarkable performance, being able to train GNNs efficiently is important task. However, training GNNs can be quite challenging in the context of large-scale graphs, which impose significant computational costs. In particular, the message passing scheme in GNNs requires, for each node, summing up the feature vectors of all neighboring nodes into one feature vector. For multi-layer Graph Convolutional Network (GCN) 3 the layer-wise propagation rule is
H (f+1)) = ar ((H)
where (t) denotes the feature representation at the 9th layer, with O X (i.c., the data matrix consisting of nodes' features of an input graph G); w-( denotes the weight matrix; a(-) non-linear activation function like the ReLU; and A denotes the adjacency matrix of G.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

If the graph has n nodes and the feature dimension is d, computing the matrix multiplication AX requires O(n2d) time and can be prohibitive in big data settings. A natural question whether is possible to train GCN to avoid the quadratic complexity scaling with n.
Another motivation for avoiding the full computation of AX is applications where the adjacency matrix A is not fully known a priori, and it must be learned via node/edge queries. For example, in large social networks (the current number of social media users is over 4.89 billion (I3), one may need to access the adjacency matrix A by querying the adjacency list of specific nodes. As another example, certain biological networks must be learned through the physical probing of pairwise interactions, which may make obtaining the entire adjacency matrix A prohibitively expensive. For example. the mapping of the neural network of living organisms (such as the connectome of C. elegans (14)) requires physical probing of the connectivity between neurons.
A natural approach avoid the quadratic complexity of computing AX is via the subsampling of the graph G 15 6 16 [7]. If the resulting subgraph is sparse enough, A can be computed very efficiently. But how many entries of A need to be observed in order to guarantee that the message passing step AX can be computed accurately enough in the context of GNNs? What kinds of graph subsampling strategies are amenable to theoretical performance guarantees on the training of GNNs?
We consider the setting of a regression task to be learned via a GCN. Let A E R be the weighted adjacency matrix of graph G = (V, A) with [V| n. Each node 2/2 E V has an associated feature vector of dimension d (< n)) and a label, denoted by (x,, U1) E RRd x IR, which can also be represented as a n x d data matrix X and label vector y. Training the GCN corresponds to minimizing the loss
C(W, := 2 lly SOCN (W,, A)
on the training data {(x;, 1))) where W denotes the GCN network weights, and faCN AW X) denotes a feed-forward computation of a GCN's output. As a first step to studying whether @ can be solved accurately a sparse subgraph of G, we focus no non-linearity o (.). Specifically, the feed-forward output is given by JGCN (w, A, == AXw, where we use w E Rd (instead of W E R xxd to indicate that the learnable parameters are in the form of vector. Hence, our goal is to solve
simple linear GCN, where there is
min
fly A
Note that one can view this optimization problem as a graph-weighted linear regression problem.
The setting of linear regression provides natural suggestions for a subsampling strategy. it known that leverage score sampling allows one to solve the linear regression problem with (1 + E) accuracy
using only O(de 2 log n) subsampled rows of the data matrix X and in time O(nd + dd
22
Nevertheless, to apply this strategy to (3). one would need to first compute AX, requiring O d)
time, so that the row leverage scores of AX can then be computed. This motivates us to propose a
two-step approach that (1) performs node subsampling to obtain an approximate observation of AX
and estimate the leverage scores of the rows of the augmented matrix [AX] - y] E R=*((+1) and (2)
performs leverage score sampling on the rows of [AX] - y] using the estimated scores.
This two-step approach is illustrated in Figure T In order to obtain an estimator of AX, our sampling scheme in the first step builds Olog n) rank-1 matrices by sampling the same amount of columns of A and the corresponding rows of X. The key idea behind this approach is taking a random subset of nodes and propagating their feature vectors to their neighbors. We provide spectral approximation guarantee for the resulting estimate of AX and an approximation guarantee for the leverage scores of the augmented matrix [AX] y] computed on this estimate. In the second step, we adopt the standard leverage score sampling for sampling rows of A, using the leverage score estimates of [AX] - y obtained from the first step. With Odde log n) sampled rows of Á and T. the training algorithm then computes AX and uses it to solve the regression problem. We show that our proposed scheme learns the regression model with learned weights deviated by most from those with full information, by observing only O(nde-2 2og n) entries of A, and in time Ond?g log n..
On real-world benchmark datasets, we demonstrate the performance improvements of the two proposed schemes over other baseline sampling schemes via numerical simulations.
MMore precisely, the time complexity is Ó((nn (AX) )++ + d²) log(s where Õ(+) hides factors polyloga-
rithmic in (n, d).
2

---

Step Estimate leverage scores
Step 2: leverage score estimate for training
AX
([4X
y]
X
X
architecture
linear GNN
t
sample Ole "lognj dy Xy'ss w.p.
sample Oldce " "logi) rows leverage sampling
perform regression task
Figure I: Two-step Algorithm: (1) perform node subsampling to obtain an estimate AX, from which we can compute ([AX) -y for = 1,... n, and (2) use the leverage score estimates to perform leverage score sampling on and y and perform regression using AX and y instead of AX and
Notation: Throughout the paper, we use A= (as;) Rax" to denote the graph adjacency matrix, where each entry value nonnegative and bounded (i.e., 0 S 0ij S M for some constant >> O). We write A; and A.j to indicate the ith row and jth column of the matrix A respectively. In addition, X 1 denotes the data (design) matrix, where the ith row in X, x,, corresponds to the node vector at the ith node V the graph G. We also assume that the absolute value of each entry value of X is bounded by M (i.e., XXy] S. M). We denote vectors by lowercase bold letter (e.g.. ", and by default, all the vectors will be in column format, hence X = (xT
)T, where (+) denotes the transpose. We indicate the ith entry of x with either [x]; or Ii, and we use [x]s denote the concatenation of x,S for all E S. We also let [m]] := {
.... (1.2,.... n). We denote by y) := -, T,YI E IR the inner product between and y. Unless otherwise mentioned, we use |/v the f2-norm for vector V. We let ||-||| max wi| |Av| denote the operator norm of matrix,
and ||A|| == V Eig ai denote the Frobenius norm. We denote the standard basis vector by e,. We write Bern (p) for Bernoulli distribution with parameter p. Lastly, we denote by nnz(A) the number of non-zero entries in matrix A.
Related Work: We provide detailed discussion of related works in Appendix A
2
Motivation:
Leverage Score Sampling for Linear Regression
First, we briefly describe key results on leverage score sampling that will be relevant to our discussion. We refer to [23 for a detailed discussion on the subject, and we provide key proofs in Appendix
Ix; where (-)f denotes Definition 1. The leverage score of the ith row of X is ((X)=== x3 (X XXX))xx
the Moore-Penrose pseudoinverse.
Intuitively, row's leverage score measures how important it is in composing the row space of X. If row has a component orthogonal to all the other rows, its leverage score is 1. Removing it would decrease the rank of X, completely changing its row space. The leverage scores are the diagonal entries the projection matrix X(X XX-XXx which can be used to show that E (X) rank(X). The following alternative characterizations of (X)) will be useful.
Proposition 1. Let X T T be the singular value decomposition (SVD) of X, where U E Rxd dadd and V ddd Then (X) |U where Ug. is the ith row of U.
Proposition 2. The leverage score can be alternatively computed as (x() maxy
[X v)
Through the use of leverage scores, it is possible to approximately solve linear regression task using a subset of the data points [18-22]. In particular, leverage score sampling allows the spectral approximation of matrix.
Definition 2. A matrix X E Baxd is e-spectral approximation of X E guxd if, for all v E Rdd
(1 (1-6) 8) |Xv|2a < || | < (1++) ||Xv||2-
Suppose we are given data matrix X :x. and vector = (ui,... Un)) of leverage score overestimates of X (i.c., (xX) < us, Vi [n.). We create matrix X by including the ith row
3

---

of X in X with probability P: min muçe 2lg ni] for some constant scaled by 1/vPs. The
following is slight modification of Lemma in [1] (which we prove Appendix D):
Lemma 1. With probability at least - n (20) for any
have
(1 - |xv|| - [Xv] - (1 c)) ([vv|
Notice that the expected number of rows in X is c||u|| log n. Hence, as long as ||u|] O(d), we only need to sample (de log n) rows of X to obtain a good spectral approximation.
Now suppose we want to solve the OLS problem minw 2 lly Xwilž where y E R" is the response
vector corresponding to X. The optimal solution is given by w (x xx- Exactly computing w* requires O(nd + d3) time.
An alternative approach using leverage score sampling would be follows. First, we note that the OLS objective can be rewritten using lly - Xw| = |[[]] y| [Y]ll- If we now have leverage score overestimates Uj...., ,, Un for the augmented matrix [X] y] E (d*(d+)), we can construct the subsampled matrix [Xsl - yss (dx(1+1) (with rows rescaled by 1/ /Pv) Then, by Lemma] solution to
arg min
[xs
ysl
(5)
1
would satisfy lyss - XsW/: 2| S lyys - ll S E)ly Xw
The exact leverage scores of X can be computed in time O(nd²) (see Algorithm in Appendix and, since |S] = (dls log n.), the reduced linear system (5) can be solved time O : log n The overall complexity for solving the OLS problem in this way is given by (nd + d3 -- log n). Notice that there is no real accelcration obtained by this procedure. But also notice that here we are using the exact leverage scores for sampling while, from Lemma [ know that this approach works as long as we have good leverage score overestimates.
3 Efficient Training of Linear GCN via Leverage Scores
In this section, we introduce an efficient training algorithm that approximately solves the graph- weighted regression task (3) by observing only (nde- n) entries A, and time O(nd2 2 log n). The necessity for introducing the efficient training algorithm comes from the
Onnd)) time complexity computing AX via matrix multiplication.
Motivated by the discussion in Section for standard linear regression, a natural approach is to attempt to compute leverage score overestimates for AX efficiently without the full computation AX. Our proposed algorithm uses the following three steps to approximately solve the graph-weighted linear regression problem.
EBTTMATTBLEVERR.GSSCOEEEAAA X) (See Algorithm I and Algorithm 2
2. LEVERKACESSCORISAAPPLIII/AA (E ([AX[ y)), Vi (m]}) (See Algorithm in Ap- pendix B
3. REGRRESSIONSOOVVERXXX (See Algorithm in Appendix
B
In ESSIMATTELIVBRAGGSSCRRSS X), a node subsampling technique is used to build an estimate AX of AX with an approximation guarantee. The estimate AX can then be used to produce provably good overestimates of the leverage scores of the augmented matrix [AX] - y]- Next, LEVERAGESCOORSS AMPLING(AA. AA (&([AX]] - yl). Vi [m]}} uses these overestimates to produce matrix Ā and scaled labels consisting of reduced number of rows of A and respectively. Using Lemma I we then show that ÅX provides spectral approximation for AX, which allows approximately solving the regression problem 3 using REGGRESSIONSOLVERR AX,ý), We note that none of the algorithmic procedures requires fully computing AX and the end-to-end algorithm has an (n log n) time dependence

---

Algorithm 1 ESSTMAATELEVERRGEESSE
A, X) via Uniform Sampling
Input: Adjacency matrix A E IRux,, data matrix X E IRaxd budget B, threshold e >> 0 Output: Leverage score overestimates for AX
1: Draw 1 ~ Bern(p), Vj E [n] independently, where p p= min n
2: Ag 1, Xj Vj
[n]
3: AX t Ej= A.
D (Bd))
y) 4: + LEVERAGESCORE((AN (1++)*
4 See Alg. in App.
O(nz [AX] - y]) + d
5: ([.([AX]- - y)), vi E [n]] + (1-a)²
6: return [.([AX]- - yl), vi [n]
3.1 First Stage: Uniform Sampling of A.
In order to efficiently estimate the leverage scores of AX, a key observation we make is that AX can be decomposed as a sum of rank-1 matrices as AX -== Az where A and Xj: denote jth column of the adjacency matrix A and jth row of the data matrix X. Notice that one can view Xy; as the effect of node j's feature vector propagated to other nodes after the one-step message passing operation. Hence, natural approximation strategy is to take a subset of these rank-1 1 matrices. First, we consider a uniform sampling strategy that samples node indices j uniformly at random and computes Xy: See Algorithm for the detailed procedure.
Given budget B (the total number of nodes that should be observed), the algorithm draws sampling indicator variables I ~ Bern(p) independently with p min [B] [B/n, 1]. For j = 1,..., n, we build
1 Ar p
A,
and then set AX = Ej AzXj: Notice that AX is an unbiased estimate of AX, because for all i,j.
1 TA = aurky = [AX]y
E
= Σ k=1
(6)
P
In addition, one can readily see that rank-1 matrices are sampled in expectation, and, with high probability, the actual number of sampled rank-1 matrices is less than (1 + a for >> 0.
Next, present spectral approximation guarantee for AX. In order establish that result, make the following mild assumption on the density of the input graph. Assumption 1. We assume that the input graph G satisfies
||a|| 2 1/222/1/2
Strictly speaking, ? is not just an assumption about the density of G, as it takes into consideration the interaction between A and X. Nevertheless, this assumption would be expected hold in cases where the graph G is dense. In particular, since ||X|| 2 d=1//|||AX|]= F. Assumption holds if ||AX]]| || 2 3/22d. Assuming that A has ön² non-zero entries and that entries of X drawn i.i.d. from some distribution with positive second moment 8², independently of A, we have ||AX]] F z /8n = 3/22 Treating d as a constant and rescaling X by a constant, we have that the assumption holds. More generally, if d a constant and
{(i.)) [n] [d] : e, AXe; 2 am)) 2 c2n
for constants Ci and (2, then the assumption holds (after a constant rescaling of the data matrix). The following stronger assumption will be needed to establish guarantees for the leverage score estimates: Assumption 2. For all (i.j) E [n] [d] e AXe; 2 cin for some constant Cj.
Notice that if the graph G is instead sparse, the problem we focus on due to the computation of AX is less severe since one can efficiently compute AX in time that depends on nnz(A).
The proof of the following theorem is in Appendix
5

---

Theorem 1. Under Assumption | if we run Algorithm with budget B = Ce log n for large
enough C 0, the estimate AX satisfies
(1 - ---)--1-e [14v]|| & ||ll||||XAv s ( VvER".
(8)
with probability at least n 2(1)
Given the spectral approximation guarantee in Theorem 1 we expect the leverage scores computed on the augmented matrix [AX] y] to be close to the leverage scores of [AX Since we want an overestimate of the leverage scores, we the the multiplicative constant
obtain
1+ 2 (12)
E(AAX] y))
4((AA]] Y)),
Vi E [ni]
The following lemma, whose proof is in Appendix asserts that [[A[ -y]) is indeed an overestimate of 4[[AX]] yy)).
Lemma 2. Under Assumption 2 Algorithm with budget Ce log for large enough outputs leverage score overestimates satisfying
1 1+E 1-e)
4(AAX| -1) s [[[AX] y)) 5
4[[AX] y)), vi [n],
(10)
with probability at least
n -911)
3.2 Second Stage: Leverage Score Sampling and Approximate Regression
The second step of the proposed algorithm exploits the leverage score estimates {[AX -y)) Vi [n] to perform the standard leverage score sampling. With Lemma [ one can readily obtain (1+ approximatii guarantee for the graph-weighted linear regression task.
Theorem 2. Under Assumption suppose we run Algorithm with budget B = C'= log for large enough C' > O, to obtain (MAAI - y)), Vi E [n]]. Subsequently, if run AWVERAGGSSOORESAMPPICC
(E<[[AX] - y)). Vi E [n]]), we obtain À gO(de x)gn)xn
and y E O(ds logn) such that REGRESSIONSOLVER AAx, provides w satisfying
y Axw 2 S (1+)) s) min lly AXwI2.
Theorem captures our main contribution. By only observing log n) rows A, we can obtain a reduced adjacency matrix Á (and corresponding label vector y. Running REGRESSIONSOLVER (ÁXX,) ý) then yields solution w with an approximation guarantee. We refer to detailed descriptions of LEVERACESSORRESAMPLNN REGreSSIoNSSOLvr (AXX, f) in Appendix B
(&([AX) - -y)), Vi E [m]}) and
Sample complexity and run-time: We note that the sample complexity O(nde log log n). First, X) exploits O log n)) columns of A and the same number of rows of X, thus accessing O(ne -2 +logn) entries of A and O n) of X. In addition, SEREAAGESCOORESAMMPLN
([AX[ - y]). Vi e [n]}) exploits O(de 2 log n) rows and thus uses Onde log n) entries. Hence, the total number of entry observations O(nde 2log n).
For the time complexity, the run time ESTMATTELaVERACESSOORRSS that of AATEERRGESOORSS A AMPLING (A, {M([AX] - y)), vi E [m]]) is O(nd2g-2 log n). Finally, we note that the run time of SOLVER (AX, ý) is O(rd + d Hence, the total run time is O =2 logn) given the setting that consider has d n.
X) is O(nde log n):
We refer
24 28 that introduce variety of acceleration algorithms for regression.

---

Extension to the more-than-two-laye case: Note that without non -on-linear,,ies, considering multi-layer network is equivalent to performing additional matrix multiplication by A to the output of the previous layer. Since our proposed scheme yields spectral approximation guarantee for matrix multiplication, an analysis similar to that of Algorithm Ianihmmmssh should be possible for this extension to obtain a (1 + e) ( apaarooximttionn guarantee, where I. denotes the number of layers. Since the approximation guarantee comes with high probability, a union bound could be used to obtain a high probability guarantee for the multi-layer case (since the order error for each approximate matrix multiplication is the same).
4 Variance Minimization via Data-dependent Random Sampling
As stated in the previous section, Algorithm D produces an unbiased estimate of AX. Yet the uniform sampling strategy used in Algorithm [ does not take the variance of the estimator into consideration. This motivates us to develop an estimator of AX that further reduces the estimator variance while maintaining unbiasedness. Due to the unbiasedness of AX, the total variance of the estimator coincides with the mean square error (MSE): Var
IEJAAX I MSB(AX).
In order to reduce the estimator variance Var(AX), we consider data-dependent sampling probabil- ity Pj that depends on the jth column of A and the jth row of X. As we view Ay as capturing the impact of node j's features after a one-step message passing operation, one way to set the sampling probability for node j is to take the norm of Ax into account. Specifically, we set
||I,I B E, E,
Pj min
(11)
Algorithm 2 describes the detailed procedure, where the only distinction relative to Algorithm is the choice of pj. We claim that this sampling probability choice minimizes the total variance (as well as the MSE) the estimator. First notice that the total variance can be expressed as
Var (AX)=
Var
[AX]] ij
Fkj
C2 (C2 = ndC1)
Pk
(1t)
1|
xx:113
Σ -=
Pk
C2
pk
(18))
||
-C2,
where (i) follows from the fact that
Var (AX](s)
Var
a
Var (1s)
12
1 Pk
a
2
Pk
Pk
Pk
2 Ci
Σ
C1 (for some constant C1)
Pk
k=1 pk
and (ii) follows from using the fact that ΣyP; 1; and (ii) follows from Cauchy-Schwarz inequality. Notice that, in order for Cauchy-Schwar to hold with equality, we must have ps proportional to | x::ll This implies that (1) is an optimal choice of the sampling probabilities.
As it turns out, the same theoretical guarantees obtained for Algorithm I Thold in this case. See Ap- pendix for the detailed proof. While the theoretical performance guarantee is the same, in settings where the budget B is very small, the data-dependent sampling probabilities p, lead to a better estimate AX,, which improves the overall performance of the algorithm.
One important observation is that is not obvious how to use the data-dependent sampling probabili- ties pj while avoiding the an Onn²) computation complexity. In particular, naive computation of ||A E [n] would require observing all n² entries of A. However, standard sampling techniques could be used in order to also estimate || ll. E [m] fairly accurate, and use those estimates to compute the sampling probabilities Pj. [n].
7

---

Algorithm 2 EESTIMATELLVVER
A, X) via Data-dependent Sampling
Input: Adjacency matrix A E IRux,, data matrix X E R budget B, threshold e >> 0 Output: Leverage score estimates [[A[ y)). Vi E [n]
1: Draw 1 ~ Bern Vj E [a] independently, where Pj = min
X B Σ, x
2: AgX3: -- 4 Pj
Ag xj: Vj E [n]
3: AX t Ej=1 Ag
D (Bd)
4: 5 LEVERAGESCORE((AA yyl (1+z)²
> Ónnz [AA -yl) + dd.)
5: ((.([AX]] - y)), vi E [n]] + (I-2)= S
6: return [[[AX]- - y)), Vie [n]
Empirical Results
We support our theoretical findings via numerical experiment to validate the performance of Algo-
rithm and Algorithm 2 in a non-asymptotic setting.
Dataset and evaluation methods: We consider benchmark datasets from the Open Graph Bench- mark (OGB) (29), Stanford Network Analysis Project (SNAP) (30), and House dataset [31] We
defer detailed explanations of the datasets Appendix H For evaluation methods, we compute the mean squared eπror (MSE), wall-clock run-time, and peak memory usage of our two proposed schemes and five baselines: (1) REGRESSIONSOLVERR EEEEEER(AAA (AX, yy with fully known and X; (2)
REGRESSIONSOLVER (AX,y) y) with parts of rows of A, say Á, obtained from sampling rows of A uniformly at random; (3) REGRESSIONSOLLERR (ÄX, y) with the exact leverage score sampling for sampling rows of A to obtain partial A: (4) REGRESSIONSOLVER (AX, y) with Algorithm | (5) REGRESSIONSOLVER (AX,3) with Algorithm 2 (6) REGRESSIONSOLVER (ĀX, with a sam- pling scheme inspired by GraphSage 6 that selects the neighborhood of feature aggregation for each node; and (7) REGRESSIONSOLVER (AX,)) with sampling scheme motivated by Graph- Saint (32) that selects subgraph of the original graph. We note that (6) and (7) are not the tools themselves 6 32 but graph subsampling strategies inspired by these tools.
Results and analysis:
Figure plots the mean squared error with respect to the budget allocation. First, we observe that Algorithm | and Algorithm 2 have a similar performance to the one with exact leverage scores sampling, and the performance of the two proposed schemes tends to be similar to that of REGreSSIoNSoLVER with full (AX,y), as the budget increases. Particularly seen from Figure [2] the error reduction on Algorithm 22 from Algorithm iis 84% when the budget is at 3%. Even for the sparse datasets, as seen from Figure 2b and Figure 2c the corresponding improvements of 86% and 67% exist when the budget is at 3% and 7% respectively. We also observe that the error decrease for Algorithm over Algorithm [] for instance Figure 26 demonstrates the existence of 91% improvement at 5% observation budget. On the other hand, we observe that other subsampling approaches-including those employed in uniform sampling, GraphSage, and GraphSaint do not perform the graph-weighter regression task as well, especially in the low-budget regime.
Mean squared error:
ovee
Extension to nonlinear GCN: While our theoretical results are for linear GCNs, we also assess the efficacy of our proposed algorithm on nonlinear GCNs. Using the ReLU activation function and one hidden layer, Figure 2d plots the mean squared error as a function of the observation budget, shown as a percentage the number of observed nodes the graph. Interestingly, we observe that running the non-linear GCN with uniformly sampled (AX, y) fails to perform the graph-weighted regression task, especially in the low-budget regime. On the other hand, we observe that running the non-linear GCN with our proposed sampling technique has a similar performance to the one with exact leverage scores sampling, and the performance of the two proposed schemes tends to be similar to that of Regression with exact (AX, y), as the budget increases. Particularly, the error reduction on Algorithm from Regression with uniformly sampled (AX. y) is 60% when the budget is at 5%
Our code publicly available online suthttps:
htt.e///wssiuuunsssgnnm.m....mee subsampling
8

---

unifomily samplka
MDE sumpled
leveraee
samppiing
Aggrntim
HISE gartthm
104
142
10
Budger
(a) Linear GCN on ogbl-ddi dataset 29 and X (b) Lincar GCN ago-Facebook dataset 30 and
generated according the Cauchy distribution.
X generated by the Cauchy distribution.
101
NISE
Mapritim
KID ganithmn
MSG Alg Aigum
MISE via Awiitthnn
VISI vie ShaptSaage rapires sempingg
via Graphiaa et inopired samsling
300
20
Buspet
(c) Linear GCN on house dataset 33
(d) Jonelinear GCN on house dataset
33
Figure 2: MSE w.r.t. observation budget (%): (1) Regression with full AX; (2) one with partially observed AX obtained by sampling rows of A uniformly at random; (3) one with partial AX based on the exact leverage score sampling; (4) and (5) one with Algorithm| and Algorithm 2 and (6) and
(7) ones with GraphSage [6] inspired sampling and GraphSaint 32 inspired sampling algorithms.
Table 1: Wall-clock time comparison end-to-end processes
Dataset
# Nodes
Edges
# Features
Wall-clock time (sec)
use of full AX
use of Algorithm and leverage score sampling
ogbl-ddi (29
4.3K
1.3M
100
1.49
1.39
ogbn-arxiv (29
169.3K
1.2M
128
299.48
7.40
Synthetic data (Gaussian)
50.0K
625.0N
500
27.28
5.77
Synthetic data (Gaussian)
100.0K
2,5B
500
107.10
8.97
Synthetic data (Gaussian)
150.0K
5.6B
500
247.70
9.96
Wall-clock time: In Table 1 we provide run-time comparison for the end-to-end training process for a regression task. We compare the wall-clock time of performing a regression task with full AX computation and that with Algorithm| [ that uses partial observations of A and X. The results show that our proposed scheme requires orders magnitude less wall-clock time for large-scale graphs than the regression with exact AX computation. In particular, for the ogbn-arxiv dataset [29), our algorithm runs about 40x faster than the regression with the exact computation of AX.
Peak memory usage: We also demonstrate the efficacy on the memory usage of our proposed algorithms. See Table Appendix Hooor details.
6 Concluding Remarks
Motivated by the prohibitive computationaalsstoage costs of running GNNs with large-scale graphs, considered the problem of subsampling nodes' information for performing regression task.

---

Our main contribution lies in providing an efficient sampling algorithm for learning linear GNNs with large number of nodes with a theoretical guarantee. Specifically, for two-layer linear GNN with adjacency matrix A E Raxa and data matrix X E IRaxd, we show that it is possible to learn the model accurately while observing only (nd log n) entries of A. This in turn yields a run-time of O(nd log n), avoiding the computation time O(n²d) of using the full matrix multiplication AX. While we view our main contribution to be of a theoretical flavor, on real-world benchmark datasets, we also demonstrate the run- tun-iime/memory improvements of our proposed schemes over other baselines via wall-clock time and peak-memory usage comparisons and MSE comparisons with respect to the observation budget.
7 Discussion and Future Work
Our result on learnability with the run-time complexity reduction is particularly significant when the number of nodes n is large and the feature dimension is much smaller compared to n (i.e., d < n). In addition, we note that the result holds even when the graph adjacency matrix and the data matrix are dense (e.g., weighted graphs where the dij entries of the adjacency matrix are non-negative real values, and not restricted to {0, 1)). Examples include pairwise similarity matrices and Erdos-Renyi graphs with constant edge probability p. which have a number of edges that scales as n², In such cases, sparse matrix multiplication would not help much, and our approach could be very useful.
Our proposed scheme is not only useful for speeding up the training of linear GCNs. It is also useful in partial observation settings where the algorithm only has access to partial view of the graph, and we provide theoretical guarantees that this partial view is still sufficient accurately perform the regression task. On the contrary, conventional speed-up methods such as sparse matrix multiplication are purely a technique for computational purposes and still require complete information on A and X. For instance, consider social network setting where, in order to know friendships between users, a learner needs to use an API query the friendships. Note that one needs fewer queries for using our method (since it requires partial A) than for using the sparse matrix multiplication (which needs full knowledge of A). the graph is large, the query complexity gain becomes larger.
Several extensions of our results are possible. First, establishing theoretical guarantees for the non- linear GCN framework, and for other target tasks, such as graph link prediction and node classification, would be of significant interest. Based on the empirical observation that our training scheme works well for a two-layer non-linear GCN, extending our theoretical results the non-linear GCN under some assumptions (as in (4 on non-linear activation functions seems to be promising. Note that our proposed sampling scheme (based approximating the leverage score) and the theoretical guarantee we made are tailored for the node-level regression task in the training phase. The reason we focused on this specific setting was so that we could take advantage of the theoretical results on leverage score sampling and linear regression in order to obtain theoretical results for GNN training using a subsampled graph. In particular, we are helped by the fact that in the linear GCN setting, the optimal solution for the linear regression is a matrix-vector product which can be expressed as (AX)) T y, where A, X, and y denote the adjacency matrix, data matrix, and the labels respectively.
For graph link prediction and classification tasks, on the other hand, since the optimal solution cannot be expressed as a matrix-vector product (due to non-linearities and different cost functions), our proposed algorithm does not extend the classification problems similarly, our theoretical guarantee may not hold anymore. For future work, we want to study whether our same subsampling process can have any guarantee for node classification task. For more general tasks, we believe that we need different sampling strategies other than the leverage score sampling and hence we may need other tools for developing a theoretical guarantee for sampling algorithms.
straightforward way, and
Another direction for future work is to establish generalization guarantees beyond the approximation errors (35]. Lastly, motivated by the fast acceleration of gradient descent algorithms (24)228 and iterative schemes that sample a subset of the entries of 21 36), developing algorithms that adaptively select subset (A, X) observe based on the information each gradient descent step would be an interesting direction to explore.
10

---

Acknowledgments
We thank Ruosong Wang, Yun Yang for the helpful discussion and the anonymous reviewers for their constructive feedback. Han Zhao and Seiyun Shin were partly supported by the Defense Advanced Research Projects Agency (DARPA) under Cooperative Agreement Number: HR001 12320012, an IBM-IL Discovery Accelerator Institute research award, and Amazon AWS Cloud Credit. The work of Ilan Shomorony was supported part by the National Science Foundation (NSF) under grant CCF-2046991.
References
[1] Marco Gori, Gabriele Monfardini, and Franco Scarselli. A new model for learning in graph domains. In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume 2, pages 729-734. IEEE, 2005.
[2] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The graph neural network model. IEEE transactions on neural networks, 20(1):61-80, 2008.
[3] Max Welling and Thomas N Kipf. Semi-supervised classification with graph convolutional networks. In J. International Conference on Learning Representations (ICLR 2017), 2016.
[4] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. Advances in neural information processing systems, 29, 2016.
[5] Yujia Li, Daniel Tarlow, Mare Brockschmidt, and Richard Zemel. Gated graph sequence neural networks. arXiv preprint 15 2015.
[6] Will Hamilton, Zhitao Ying. and Jure Leskovec. Inductive representation learning on large graphs. Advances in neural information processing systems, 30, 2017.
[7] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful are graph neural networks? arXiv preprint arXiv:1810.0088 2018.
[8] Carl Yang. Jieyu Zhang, Haonan Wang, Sha Li, Myungwan Kim, Matt Walker, Yiou Xiao, and Jiawei Han. Relation learning social networks with multi-modal graph edge variational autoencoders. In Proceedings of the 13th International Conference on Web Search and Data Mining, pages 699-707, 2020.
[9] Marinka Zitnk and Jure Leskovee. Predicting multicellular function through multi-layer tissue networks. Bioinformatics, 33(14):110-i118, 2017.
[10] Xiang Yue, Zhen Wang, Jingong Huang. Srinivasan Parthasarathy, Soheil Moosavinasab, Yungui Huang, Simon M Lin, Wen Zhang, Ping Zhang, and Huan Sun. Graph embedding on biomedical networks:
methods, applications and evaluations. Bioinformatics, 36(4): (2044)):21-11555 2020.
[11] Zhenqin Wu, Alexandro E Trevino, Eric Wu, Kyle Swanson, Honesty Kim, H Blaize D' Angio, Ryan Preska, Gregory W Charville, Piero D Dalerba, Ann Marie Egloff, et al. Graph deep learning for the characterization of tumour microenvironmentt from spatial protein profiles in tissue specimens. Nature Biomedical Engineering, pages 1-14, 2022.
[12] Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec. Graph convolutional neural networks for web-scale recommender systems. In Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, pages 974 983, 2018.
[13] Statista. Number of social media users worldwide from 2018 to 2022, with forecasts from 2023 to 2027, 2022. URL httpp:://ww.ssatista.u ofmberr-of-orlddiddd social maatallnetwor---ssee http://www..aaaaaacaaaaaaaaaaaaa
[14] Steven J Cook, Travis A Jarrell, Christopher A Brittin, Yi Wang, Adam E Bloniarz, Maksim A Yakovlev, Ken CQ Nguyen, Leo T-H Tang, Emily A Bayer, Janet S Duerr, et al. Whole-animal connectomes of both cacnorhabditis elegans sexes. Nature, 571 2019.
11

---

[15] Wenbing Huang, Tong Zhang, Yu Rong, and Junzhou Huang. Adaptive sampling towards fast graph representation learning. Advances neural information processing systems, 31, 2018.
[16] Jie Chen, Tengfei Ma, and Cao Xiao. Fastgen: fast learning with graph convolutional networks via importance sampling. arXiv preprint arXiv: /801 10247, 2018.
[17] Ziqi Liu, Zhengwei Wu, Zhiqiang Zhang, Jun Zhou, Shuang Yang, Le Song, and Yuan Qi. Bandit samplers for training graph neural networks. Advances in Neural Information Processing Systems, 33:6878-6888 2020.
[18] Kenneth L. Clarkson and David P Woodruff. Low-rank approximation and regression in input sparsity time. Journal of the ACM (JACM), 63(6): 2017.
[19] Mu Li, Gary L Miller, and Richard Peng. Iterative row sampling. In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, pages 127-136. IEEE, 2013.
[20] Jelani Nelson and Huy L Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings. In 2013 ieee 54th annual symposium on foundations of computer science, pages 117-126. IEEE, 2013.
[21] Michael B Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, and Aaron Sidford. Uniform sampling for matrix approximation. Proceedings of the 2015 Conference on Innovations in Theoretical Computer Science, pages 181 181-190, 2015.
[22] Michael Cohen. Nearly tight oblivious subspace embeddings by trace inequalities. In Proceedings of the twenty-seventh annual ACM-SIAM symposium on Discrete algorithms, pages 278-287. SIAM, 2016.
[23] David P Woodruff et al. Sketching as a tool for numerical linear algebra. Foundations and Trends@ in Theoretical Computer Science, 10(1-2):1-157, 2014.
[24] Zeyuan Allen-Zhu, Zheng Qu, Peter Richtárik, and Yang Yuan. Even faster accelerated coordinate descent using non-uniform sampling. In International Conference Machine Learning, pages 1110-1119. PMLR, 2016.
[25] Yurii Nesterov and Sebastian U Stich. Efficiency of the accelerated coordinate descent method on structured optimization problems. SIAM Journal on Optimization, 27(1): 10-123, 2017.
[26] Roy Frostig, Rong Ge. Sham Kakade, and Aaron Sidford. Un-regularizing: approximate proximal point and faster stochastic algorithms for empirical risk minimization. International Conference on Machine Learning, pages 2540-2548. PMLR, 2015.
[27] Hongzhou Lin, Julien Mairal, and Zaid Harchaoui. A universal catalyst for first-order optimiza- tion. Advances in neural information processing systems, 28, 2015.
[28) Naman Agarwal, Sham Kakade, Rahul Kidambi, Yin Tat Lee, Praneeth Netrapalli, and Aaron Sidford. Leverage score sampling for faster accelerated regression and erm. arXiv preprint arXiv:1711.08426. 2017.
[29] Weihua Hu, Matthias Fey, Marinka Zimik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec. Open graph benchmark: Datasets for machine learning on graphs, 2020. URL wttthttss:/// http:///bbbb..t.. stanf...ssaoorr.oo.oo.ooooooooo. fed.uuuurrrddud
[30] Jure Leskovec and Julian Mcauley. Learning to discover social circles in ego networks. Advances in neural information processing systems, 25, 2012. URL htttps:/snapp stanford. edu/
index. html
[31] R Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3):291-297, 1997.
[32] Hanqing Zeng. Hongkuan Zhou, Ajitesh Srivastava, Rajgopal Kannan, and Viktor Prasanna. Graphsaint: Graph sampling based inductive learning method. arXiv preprint arXiv: 190,0419,, 2019.
12

---

[33] Sergei Ivanov and Liudmila Prokhorenkova. Boost then convolve: Gradient boosting meets graph neural networks. arXiv preprint arXiv:2101.08443 2021.
[34] Aarshvi Gajjar, Christopher Musco, and Chinmay Hegde. Active learning for single neuron models with lipschitz non-linearities. In International Conference Artificial Intelligence and Statistics, pages
PMLR, 2023.
[35] Jason D Lee, Ruoqi Shen, Zhao Song, Mengdi Wang, et al. Generalized leverage score sampling for neural networks. Advances in Neural Information Processing Systems, 33:10775-10787, 2020.
[36] Seiyun Shin, Han Zhao, and Ilan Shomorony. Adaptive power method: Eigenvector estimation from sampled data. In International Conference on Algorithmic Learning Theory, pages 1387-1410. PMLR, 2023.
[37] David C Hoaglin and Roy E. Welsch. The hat matrix in regression and anova. The American Statistician, 32(1): 1978.
[38] Samprit Chatterjee and Ali S Hadi. Influential observations, high leverage points, and outliers in linear regression. Statistical science, pages 379-393, 1986.
[39] Petros Drineas, Michael W Mahoney, and Shan Muthukrishnan. Sampling algorithms for regression and applications. In Proceedings of the seventeenth annual ACM-SIAM symposium on Discrele algorithm, pages 127-1113, 2006.
[40] Daniel Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. In Proceedings of the fortieth annual ACM symposium on Theory of computing, pages 563-568, 2008.
[41] Christos Boutsidis, Michael W Mahoney, and Petros Drineas. An improved approximation algorithm for the column subset selection problem. In Proceedings of the twentieth annual ACM-SIAM symposium on Discrete algorithms, pages 968-977. SIAM, 2009.
[42] Xin Liu, Mingyu Yan, Lei Deng, Guoqi Li, Xiaochun Ye, and Dongrui Fan. Sampling meth- ods for efficient training of graph convolutional networks: A survey. IEEE/CAA Journal of Automatica Sinica, 9(2):005-23 234, 2021.
[43] Difan Zou, Ziniu Hu, Yewen Wang, Song Jiang. Yizhou Sun, and Quanquan Gu. Layer- dependent importance sampling for training deep and large graph convolutional networks. Advances in neural information processing systems, 32, 2019.
[44] Jianfei Chen, Jun Zhu, and Le Song. Stochastic training of graph convolutional networks with variance reduction. arXiv preprint arXiv: 10. 10568, 2017.
[45] Wei-Lin Chiang, Xuanqing Liu, Si Si, Yang Li, Samy Bengio, and Cho-Jui Hsich. Cluster-gen: An efficient algorithm for training deep and large graph convolutional networks. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pages 257-266, 2019.
[46] Koen Jochmans and Martin Weidner. Fixed-effect regressions on network data. Econometrica. 1745:15533 2019.
[47] Arun Venkitaraman, Hermina Petric Maretic, Saikat Chatterjee, and Pascal Frossard. Supervised linear regression for graph learning from graph signals. arXiv preprint arXiv:7811.015555 2018.
[48] Can M Le and Tianxi Li. Linear regression and its inference on noisy network-linked data. Journal of the Royal Statistical Society Series B: Statistical Methodology, 8465):851-11855, 2022.
[49] Joel A Tropp et al. An introduction to matrix concentration incqualities. Foundations and Trendy@ in Machine Learning, 8(1-2):1-230, 2015.
13

---

A Related Work
Leverage score sampling for linear regression: The history of leverage scores traces back to linear regression problems in statistics [37 [37
38). The statistical leverages of a matrix indicate the importance of rows in solving a least-squares problem. It has been well-know that leverage score sampling [39 40 uses leverage scores to sample a subset of the rows of the matrix and solve the regression task with (1 + g) accuracy using only Ode= log d) subsampled rows of the data matrix X and in time (nd + ) and there have been extensive works on reducing the time complexity well the sample complexity [41 [8 (2). However, linear regression, which involves a matrix multiplication of A and X, taking ((n²d) time. To the best of our knowledge, our work is the first one that addresses the setting of graph-weighted linear regression and proposes an algorithm with subquadratic time complexity in n.
note that our setting is a graph-weighted
Sampling on Graph Convolutional Networks (GCNs): Most of the related works focus on reducing the receptive field size via subsampling neighbors of nodes in the input graph. The motivation is that the neighborhood size can grow exponentially with respect to the number of layers due to the message passing operation of each layer. In contrast, we focus the effect of partial observation one hidden layer and show that node subsampling can reduce the run time as well as the sample complexity while ensuring the desired performance guarantee. We also highlight that our work provides a theoretical guarantee, whereas many of the cited works below provide empirical results without theoretical guarantee. For readers interested in works related sampling GCNs, we refer to the comprehensive survey in [42]. Here we provide a brief summary.
To mitigate the "neighbor explosion" problem, most of the related works 6
take layer- wise sampling approaches sample neighbors the nodes in the previous layer in a top-down manner. Specifically, layer-wise wise samplers [16 15 [3] adopt importance sampling sample neighbors in the lower layer given each node in the upper layer, assuming that the sample size for each layer is independent of each other.
16
15
43
Second, node wise samplers randomly sample neighbors the lower layer given each node in the upper layer. In addition. 7V proposed bandit-based approach to perform node-wist wiie sampling in order to optimize the variance in the vein of layer sampling approaches. See also [4), which proposed a variance-reduced estimator to reduce the receptive field of single node.
Third, there have been several works [45 32
that adopt graph-wise sampling approaches. In particular, the proposed schemes first partition the node set of a graph 45 or sample subgraphs (22, and then train models on those partitions subgraphs in batch mode [3). They show that the training time of each epoch may be much faster compared with layer-wise sampling approaches.
Graph-weighted linear regression: We note that there have been a few works that address linear regression network-linked data 46.48). Yet neither sample complexity analysis nor time complexity analysis is provided.
B
Basic Algorithms
Algorithm 3 LEVERAGESCOREE(X
Input: Data matrix X E Raxd
Output:
Leverage score (X Vi [n]
1: Compute K == X
b O(d2)
2: +
O
3: U X(VE-1)
(nd
4: return 6(x)= | Vie [n]
b ((dd)
'The time complexity is (nd + dd) in total, and more precisely, it is O((nnz(AX) + dd) log(s>1));
|8 22 for details.
14

---

Algorithm 4 LLEVERAGESSORRESSMPPINNG (( - yD,, vi E [n]})
Input: Adjacency matrix A E R leverage estimate of AX (E]([AX] - yl), VE E [n]]
Output: Ā E Ro(d= logm)) scaled labels y E IRo log
I: Draw 1 Bern(py) Vj [n] independently, where Py ((AX -yl)e-2logm
2: - Σ j=1 VPa 0j4j
n 1
3: y - اj=1 /05 yje
4: refurn (A,
Algorithm 5 REGRESSIONSOLVER(AXX
Input: Matrix AX E ORd logg)xdd scaled labels y E RRdee log Output: Weight vector w E TRd
1: w + ((A)) (A
2: return w
C Proof of Proposition 2
Define max.
Hence, it follows that
Xvi
Vv, T (xxx,) v S t. T (XT x)
Equivalently, this means that t is the smallest value such that
x;x; L6-7 X X.
Let w Xv. Since X has full column rank, know that X1w Without loss of generality, assume w 2 1, then we can equivalently transform the above constraint as
Vw, ||w||2 = 1, w ((X)) x.x X < l. X is rank-one matrix, hence suffices to choose such that / 2 tr ((xt) x TexxxX' xt) = tr (x xxx]xxxxx)= x.) = tr (x7 (X (Xtx)x.)
However, note that
= tr { (X (xx xx)x xx xx (X).
This completes the proof.
D Proof of Lemma 1
In this section, we prove Lemma D which yields good approximation by sampling with probability proportional to the leverage score overestimates:
Raad and vector (w1, un) of leverage score overestimates of X (i.e., (X) S u, Vi E [n],, suppose we create matrix X by including the ith row of in X with probability D3 = min[1, cu,e log n] for some constant e, scaled by 1/ /Ps. With probability at least n -2(1) for any v, we have
Lemma 3. Given data matrix X (x];"..;xx)"
(1 =) ||xv|| < |[v]|| S (1 (1+ (|vv].
To prove the lemma, we rely the following matrix Chernoff bound, which is variant of the matrix Chernoff bound in (21 Lemma 11] (also, see )49).
Suppose Wi....y Wn, are independent, random positive semidefinite matrices with size d' x d. Define Z Σド= W. If W, = R E[Z], then
Pr 2 (1 + EZ) < d exp
(12)
3R
62
Pr 3 (1 - E] < d exp
(13)
2R
15

---

where A z B is equivalent to saying v TovvvTAA- (A- Byv 2 0 for matrices A and B and for all v.
To apply this inequality, define W3 == 1 x,,x), which is sampled with probability pj- We first know
that Z is an unbiased estimate of X X since
E
Pj Pj
Xjxj
- E = =x
First, consider the case where Pj < 1. This implies that
((x) 1 Ug S E E log
In order to invoke the matrix Chernoff bound, we only need to verify that W, (2 0) is spectrally upper bounded:
1,(x)-Xxx
X
W cuje- log x,X 3
clj(X)e- log xjx
(14)
logn
CE log
where the last inequality holds due to the proof of Proposition
2 Now, the statement follows from
the matrix Chernoff bound with Wi (x (x-'x)=11 1 X) /2 and R cs llogn' to conclude
that
(1-6) X
W (1+0) X j=1
with probability at least 1 c//3 d Note that
n 0(1), for some c and the regime where union-bound technique derive the desired high probability guarantee.
c/3 + n -/2
Now consider the case where P, 1, where W, := x,x is sampled with probability In this case, [14 generally does not hold. However, since sampling W, with probability 1 is equivalent to
selecting W()). ..... wrle n) each equal to logg x,x with probability 1, it follows that
X
( y
i=
logn.
22-21og1 logn
Here we assume that cs log n is integer for simplicity. Replacing Y, with E log yr(8) yields the desired result since it does not change B[Z]. This completes the proof.
E Proof of Theorem 1
Theorem 6.1.1]: Suppose We will use the following statement of Matrix Bernstein's 's inequality [49 W..... W,, are independent, random matrices with common dimension d, x d2 satisfying
E[W5] = o, and ||y $ L, Vje [m].
Define Z = -E W,. Then for all / >>> 0,
1
|||
Σ 1 (2((2) log(d1 d2) + log(d1 + dz); 3
(15)
1/2
Pr
Σ W3
21 5 (dl + dz) exp (~x(2) 1/2 Tt
(16)
where
}.
max Σ
E [ wall
u(z) max { | 2z 1|
||
[z
W WT]]
Σ
16

---

To apply this inequality, first recall that Algorithm La Xx with probability
samples Agx to get rank-1 matrix A ==
B n
P. min
Defining W == Az Az
j:. we note that B[W,, | = 0 as shown in (6. In addition, assuming
that entries of A E RRxx and X E Baadd satisfy AAy
_ .M and X,yl I M, we then have
1 L - A
1
[wyl =
1
A X
=
Pj
Pj
1
- Pj
Xx:11
max
n
1
||A
||,
||g
|Xy:ll
(1) 3/2
B
where (i) follows from the fact that are interested in the regime where E = =({); and (ii) follows
from the fact that ||A S vnM and || S VdM. We can also bound | T Z]|l as
B 21|1
Σ
W W
Σ
E
W
w,ll
j=1
Σ
E
1/
-1) 2 1) )
XTAJAAX
Pj
1
||x ||
j=1 Pj
n2 dM
B
and B[ZZ [ can be similarly bounded, implying that v(2) <
n aaa B
Setting & = s||AX|| yields
1A
Pr
Σ W
2 s-11AX1]
s 2(n + d) exp
nddM4 1
8//2
|AXI|
B
= 2(n + d) exp
7ddM4
|AXX
|AX
(Mi)
2(n + d) exp
M + MMPE
(in)
4n exp
M + 1 M
where (zi2) follows from Assumption 1 that ||X]] 2 n //2 and (iv) follows from the setting that
we focus on d < n. Hence, as long as we set B C= log n for some constant (> 3 > 0). the exponent of the Matrix Bernstein inequality can be bounded as -C" log n for some C" > 0, and thus can bound the deviation probability to be n (2(1)
Now using the reverse triangle inequality,
|XX| - || S S|AX
AXII S EXIAAXII,
17

---

we can conclude that the following spectral approximation holds
(1 [Ix & ||AX| - (1 + +)]l ||,
with probability at least 1 11) This completes the proof.
Proof of Lemma
We will use the following statement of Bernstein's 's inequality to get the desired result: If W1,-... Wi..... Wn, are independent zero-mean random variables such that [W;] s M with probability 1. Then for all
1 0,
1 [ Mt
Pr
Σ
W&
2t
S
exp
(17)
To apply this inequality, we first provide the following proposition which will serve as key ingredients for proving the lemma.
Proposition 3. Fix (i, [n] x [d]. Let Wi = de pk
A Then for all is [n].
the following two hold:
1. W S nM² B
2. E W
nM B
Note that these two hold for all (i.5) E [m] [d].
Proof. We first bound [W as:
Ia
Is
W
[A Xn:les
[A
Xa:ls
Xades
Pk
Pk
Xx:lisl
Pk
Ti
1
-22
max
B
nM2
B
where (i) follows from the fact that we are interested in the regime where =(1). We also bound the second moment of W& as:
(Le
-1)
E
W
LA.xXA::1
Ps
XR:
Pk
n.M*
B
where we used the facts that E
Pj- This completes the proof.
Now setting = = e, AXe, and applying the Bernstein's 's inequality now yields:
Pr (i.j)
Σ We
2t
= Pr (1.3) E [m] x [d] | AXe, - -e exAx,, 20 AXe;)
< 2 Pr AXe; o AXe3 e. AXe,
18

---

2 (e ,Xe,) )²
2n2 exp
Tt BM
+ 3
nA
(e AXe,
B
(i)
B
2n exp
M + 3 M
where (i) follows from a union-bound technique; and (ii) follows from Assumption 2 that o, AXe, = 22()), v(i,j)
E [n] x [d]. Setting B loog n for sufficiently large C > M > 0.
3
one can observe that the last bound is upper bounded by n n-0(1), Since the other direction holds similarly, we conclude that with probability 1 - n -((1), with budget of B go n,
(1 e) eg AXe, AXe, (1 +) e, AXe5, v(i,)) [n] [d].
For the augmented matrix [AX] y]. the above then yields
(1 a) e, [AX] yie, f c [AX| yle, s (1 + E) o |AX| yyee, v(i,)) [n] x [d].
Pick v E Rd+1, where it can be represented as
d+1
for some 01,,... O.411 From the above, one can readily see that with probability at least - n 9(1,
(e, [AX]
y/v)
1 +
(e, [AX ylv)²
1 +
(e [AX] y]v)²
<
Vv.
||AX| v
[[[AX
y)v
|[AXX|
vlv12
(e] -yvv* and w arg max
(e) [AX y/w)² from the above, ||
Choosing v² arg maxy
1 >/w
we have
1 + 1 -
[AX y)vt)2
I (1-3)
[[AAx-y)vvi ||AX| - ylv
1 )
(e. y)w*)2
||[AX|
y/v+13
|[AX| y/w"
where the last formula is ( by the choice of w". In similar manner, we obtain
(o [AX
y/w*)
1+
(c. [AX[AX] y/w")²
(c. (c. CAAAX - y]v*)2
|[AX| y/w 12
1-
|[[AX] y]w*1B
|[X - y]v=
Referring to Proposition
and the leverage score estimate (( satisfying
1
1 1+8
(e) AX y]v)²
{[AAX] -y)) =
4((AX) - -y)
max
1-E)
1-6
|[AX
y/vll3
we get the desired result of
1+E 1-
4(AAX|- - yl) < [[[AX| - yy)) S
4[[AX] -y),,
This completes the proof.
G Theoretical Guarantees for Algorithm 2
In this section, we provide theoretical guarantees for Algorithm As mentioned Section 4 we obtain the same series of theoretical guarantees as those for Algorihhm For completeness, we leave the detailed statements below.
Theorem 3. Under Assumption if we run Algorithm 2 with the budget B = C² log for large enough C > 0, the estimate AX satisfies
(1 - - ||Ax||| $ ||Xv|| (1+ (1+t)-11AXv]], Vv E Rd
(18)
with probability at least - n- R(1)
19

---

Lemma 4. Under Assumption |
with the budget B = C logn n for large enough C > 0 outputs the leverage score overestimates with Lte multiplicative ratio; i.e.
Algorithm 2
4 1 ++ 1-
4[AAX] y)) s ([[A yl) =
4([AX] -y)),
Vi [n],
(19)
with probability at least n n-1(1)
Theorem 4. Under Assumption suppose we run Algorithm 2 with the budget B C'e - log for large enough C > 0, to obtain {4[AAX] - y)). Vi E [n]]] Subsequently, if we run LEVERAGESCCRESAMMLINE and y E R l) such that REGRESSIONSOLVER A provides
A, (É([AX] yl), Vi E [n]]), we obtain A RO(de 1ogn)xn
min
Ax S (1 + + min lly AXw
We note that it suffices to provide the proof of Theorem the rest follows in similar manner to those for Algorithm
Proof of Theorem 3 As in Appendix E we will use the Matrix Bernstein's 's incquality described
in (16). Recall that Algorithm
2 samples Ay to get rank-1 matrix A, Pi
12 Agx with
probability as in
|Am ΣS || |||l
Pj min
B
such that Σ,p; - B.
Defining Wj = A.j. A X,: and we note that E[W;] = 0 as shown in (6). addition, assuming
that entries of A E TRaxn and X E Rnadd satisfy |A, S M and X,sl I M, we then have
1 X, AイXx1 Pj
12 X
[[
Pj
1 15:1l
1
Σ, | X||l
|
1 || ||x;:1
max
B
T
max
Ø=1
(
n 1/221/2 n 1/
<
max
B
n 3/2 M2
B
where (i) follows from the fact that ||A S VrM and || S vdM; and (ii) follows from the fact that we are interested the regime where B = o(n). We can also bound E[Z Z as
LA -AX,)-AAXX))
1/ Liyx - (LA
E I
||||
E
Pj
Dj
(A pj 20A2(,A)))))) ( (AgX,) + (A),)))))))
p3
20

---

1) (1l X3x
Pj
-
(11 xx)x3:)
B
||
||I|
Σ
A
A ||
B
X
2
Σ
|
|
B
n dM
B
Also, one can readily see that E| can be similarly bounded, implying that v22) S Setting & ||| yields:
ddst
B
Pr
Σ
||AX|||
s 2((n+ d) exp
2
"ddd1 B
13/221//2
M
E
=|AAX||
$ 2(n + d) exp
ndd |AX
13/2d-12
+
|XXX
(sit) 2(n d) exp
_
M4 + M
(i0)
2 1
4n cxp
M + M
where (iái) follows from Assumption that ||A]] > 222 and (iv) follows from the setting that
we focus d < n. Hence, as long as we set B Ce-2 log n for some constant (> 3M > 0). the exponent of the Matrix Bernstein incquality thus can bound the deviation probability to be n -Q(1)_
be bounded -C" log for some C" >> 0, and
Now using the reverse triangle inequality,
|
||AX]||
IIAX AXII S Ee||AX||,
can conclude that the following spectral approximation holds
[1= < ||A|| s (1 + =)] AXX
with probability at least n Q(1), This completes the proof.
H Further Implementation Details
In this section, provide additional details the explanations of the datasets, wall-clock time with the adoption of sparse matrix multiplication, and peak memory usage comparisons.
Datasets:
We consider benchmark datasets from the Open Graph Benchmark (OGB) 29) Stanford Network Analysis Project (SNAP) [30], and House dataset [31]. For the MSE comparisons, we consider (1) ogbl-ddi dataset from OGB, (2) ego-Facebook dataset from SNAP, and (3) House dataset. We note that while the third dataset has its data matrix (i.e., concatenation of feature vectors) and labels, the first two datasets do not have feature matrices and target variables. Accordingly, conduct controlled experiments, synthetically generate data matrix X, weight vector w, and noisy labels
21

---

Table 2: Wall-clock time comparison on end-to-end processes with sparse matrix multiplication
Dataset
# Nodes
# Edges
# Features
Wall-clock time (sec)
of full AX
use of Algorithm | and leverage score sampling
ogbl-ddi (29)
4.3K
1.3M
100
1.02
0.84
ogbn-arxiv (29
169.3
1.2M
128
6.69
0,86
Synthetic data (Gaussian)
50.0K
625.0M
500
2452.14
58,30
Table 3: Peak memory usage comparison on sampling and regression steps
Dataset
# Nodes
Edges
# Features
Peak memory (MiB)
of full AX
of Algorithm
and leverage score sampling
ogbl-ddi 29
4.3K
1.3M
100
1.45
0.17
ogbn-arxiv 29
169.3K
1.2M
128
330.75
160.45
Synthetic data (Gaussian)
50.0K
625.0M
500
240.43
0.17
Synthetic data (Gaussian)
100.0K
2.5B
500
431.82
0.52
Synthetic data (Gaussian)
150.0K
5.6B
500
865.25
19.39
y by setting the linear relationship between labels and features: y := Xw + n. Here n denotes an additive Gaussian noise with parameters (M, o)) (1, 10). Based upon the size of the two datasets described below, we consider the data matrices X E Raxd (with n E {4267, 4039) and d = 100), where each row (i.e., each node's feature vector) is drawn according to the Cauchy distribution with parameters (æo, 7) (10, 100), in an i.i.d. manner. Here zo denotes the location parameter, indicating the location of the peak of the distribution, and denotes the scale parameter, specifying the half-width at half-maximum (HWHM). Also, consider IRd generated by the Cauchy distribution with (zo, 7) = (10, 100).
ogbl-ddi dataset [29 : This dataset represents a undirected and unweighted graph which models the drug-drup interactions. Each node corresponds to an FDA-approved or experimental drug. An edge indicates interactions between drugs that can be interpreted as the difference between the joint effect of taking the two drugs together and the expected effect when drugs act independently of each other. As there are neither feature vectors nor labels that are available, we synthetically generate those as mentioned above.
ego-Facebook dataset [30]: This dataset models social network that represents friendship of "circles" between users from Facebook. Examples of such relationships between users include students of common universities, sports teams, relatives, etc. Notice that the graph adjacency matrix more sparse compared ogbi-ddi dataset. As above, since there exists feature vectors and labels that are available, we synthetically generate those to perform regression task.
House dataset [31]: This dataset represents graph in which nodes correspond the properties; the existence of edges means that two properties are located close to each other; and the target label is the property's 's price. In addition, there exist node features representing MedIne, HouseAge, AveRooms, AveBedrms, Population, and AveOccup. In particular, we use the refined dataset from 33
Furthermore, for the wall-clock time and peak memory usage comparisons, additionally use ogbl- arxiv (29 dataset, representing the citation network between all Computer Science (CS) arXiv papers (although the dataset is tailored for classification of 40 subject areas) and synthetic dataset where all X, y) are generated according to Gaussian distribution with size specified Table Table and Table
Additional results and analysis:
Wall-clock time with the adoption of sparse matrix multiplication: Notice that as compared to results in Table [] Table demonstrates that the computation time of both our approach and using exact AX become faster due to the adoption of sparse matrix multiplication and the reduction time was larger for using exact AX (especially for the ogbn-arxiv dataset which is sparse). However, we observe that our proposed scheme is faster than the regression by using the exact multiplication of AX. Furthermore, for the synthetic dataset where A and X are generated according to a Gaussian
22

---

distribution and hence A is dense, we observe that the sparse matrix multiplication does not help speed up the computation time and the matrix multiplication with seipy.sparse package took even longer than the matrix multiplication with numpy package.
Peak memory usage: In addition to run-time, expect proposed approach to offer significant improvements in terms of memory usage. As simple experiment in terms of memory usage, we provide the peak memory usage comparison for the same datasets during the sampling and regression steps of the algorithm. As shown in Table g in the best case, the proposed algorithm requires 1414x less memory than the algorithm using the exact computation of AX.
Computing architecture for wall-clock time and peak-memory usage comparisons: For fair comparison on Table | Table 2 and Table we use the same regression solver and use the same specification of 48 cores of an x86 64 processor with 503.74 34GB memory.
23


Reviews:
[{"review_id": "rv4wL9friq", "reviewer": "Reviewer_S2HF", "summary": "The paper presents an efficient GNN training solution through node subsampling and leverage score sampling, which is proven to be efficient in learning a regression model with bounded entry access and running time.", "questions": "N/A", "limitations": "N/A", "rating": 3, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "* The proposed technique leads to a proven efficient approach for GNN training, with potentially many real applications. \n\n * Experimental study on real graph data validates the proposed technique on real world graphs, showing significant improvement over baseline sampling strategies.", "weaknesses": "* It is interesting to see the proposed approach for GNN training, as efficiency of GNN on big graphs is a headache for many real world problems. However, although being different approaches, we can find the idea of subsampling used in various GNN models, but it is not well compared in this paper. For example, GraphSage performs neighbor subsampling and GraphSaint subsamples subgraphs from a big graph. The authors should compare the proposed approach with those different sampling approaches to show the pros and cons. \n\n* The experimental results show the significant efficiency of the proposed method. However, since it is for GNN training, the only evaluation on MSE would draw some concerns, if MRS is the only metrics GNN should be concerned about. If not, the authors should expand the experimental study for a comprehensive understanding/evaluation. \n\n* Compared to MSE, it is critical to see how much typical GNN machine learning tasks would be impacted by the subsampling. For example, graph link prediction, node classification, label prediction, etc. This paper does not include evaluation form those perspectives. \n\n* Although the big-O notation of complexity analysis is quite helpful, the experimental study of efficiency should be provided, if the training time is also significantly reduced with comparable performance. \n\n* It is a bit concerned that, for studying GNN efficiency, the authors selected two small graphs from OGB, although much bigger graphs exist.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "W1. Comparison with other baselines (e.g., GraphSage and GraphSaint) and the pros and cons:\n\nWe thank the reviewer for the suggestion. Please see the comments in the common response section and the attached pdf for additional experiments.\n\nW2. MSE was the only evaluation criterion:\nSince we focused on a regression task, we believe that the MSE is a natural performance metric. Furthermore, our theoretical guarantee is in terms of the L2 error in the regression task, and the main goal of our experiments was to validate the theoretical guarantees. In future work, once we consider other learning tasks, such as node classification, different performance metrics may be needed.\n\nW3. Applicability of our proposed algorithm to graph link prediction and node classification:\nWe thank the reviewer for suggesting interesting directions for future work. Our proposed sampling scheme (based on approximating the leverage score) and the theoretical guarantee we made are tailored for the node-level regression task in the training phase. The reason we focused on this specific setting was so that we could take advantage of the theoretical results on leverage score sampling and linear regression in order to obtain theoretical results for GNN training using a subsampled graph. In particular, we are helped by the fact that in the linear GCN setting, the optimal solution for the linear regression is a matrix-vector product which can be expressed as $((AX)^{\\top}(AX))^{\\top}\\cdot (AX)^{\\top} \\mathbf{y}$, where $A$, $X$, and $\\mathbf{y}$ denote the adjacency matrix, data matrix, and the labels respectively.\n\nFor graph link prediction and classification tasks, however, since the optimal solution cannot be expressed as a matrix-vector product (due to non-linearities), our proposed algorithm does not extend to the classification problems in a straightforward way and similarly our theoretical guarantee may not hold anymore. That being said, we agree that both node and edge classification problems are interesting. For future work, we want to study whether our same subsampling process can have any guarantee for a node classification task. For more general tasks,  we believe that we need different sampling strategies other than the leverage score sampling and hence we may need other tools for developing a theoretical guarantee for sampling algorithms. \n\nW4. Lack of run-time comparison:\n\nWe thank the reviewer for pointing out this issue. We leave comments in the common response section and provide additional experimental results in the attached pdf.\n\nW5. Experiments on two small graphs:\n\nWe would like to first highlight our main contribution lies in providing an efficient sampling algorithm for learning linear GNNs with a large number of nodes with a theoretical guarantee (please see the common response for details). We agree that we include experiments on rather smaller graphs, so we have run further experiments to provide wall-clock runtimes for the experiments and also empirical results on nonlinear GNNs for various sizes of graphs. We refer to the common response above and the attached pdf file for detailed experiment results. We will include additional experiments on large graphs in the main draft if the paper is accepted."}, {"review_id": "duyXUyqlCE", "reviewer": "Reviewer_J15Y", "summary": "The authors proposed a sampling method to train graph neural networks efficiently. However, the current implementation of graph neural networks is based on the sparse matrix multiplication and thus the complexity is not O(n^2d). The authors only analyze the complexity theoretically without experimental supports. Both of major issues make the contribution weak.", "questions": "1. What is the computational time of this approach compared with sparse matrix multiplication?\n2. Besides MSE, how about applying this to the real semi-supervised classification tasks? Will this lead to the same prediction accuracy?\n3. Is it possible to use this for link prediction tasks when dropping edges for validation? Is there still any theory guarantee in this case?\n", "limitations": "N/A", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. Theoratically, the proposed method is efficient compared with the dense implementation of graph neural networks.\n2. The experiments are conducted on the large scale datasets with different variant of the methods.", "weaknesses": "1. The statement about the complexity of the graph neural network is not correct. The current implementation is based on the sparse multiplication, so the complexity is not O(n^2d).\n2. The computational time comparision is missing, which is essential for this work.\n3. There is lack of baselines, it only compares with its own variants.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "W1. The statement about the complexity of the graph neural network is not correct. The current implementation is based on sparse multiplication, so the complexity is not O(n^2d).\n\nWe thank the reviewer for pointing this out. The main motivation for our graph subsampling procedure is settings where the adjacency matrix is dense (with $\\Theta(n^2)$ non-zero entries), and sparse multiplication would still require time $O(n^2d)$. The reviewer is correct that, whenever $A$ is sparse, sparse matrix multiplication could avoid the $O(n^2d)$ complexity. We will correct the paper to be clear on this point.\n\nW2. Lack of run-time comparison:\n\nWe thank the reviewer for pointing out this issue. Please see the common response and the attached pdf file for further experiments.\n\nW3. Comparison with other baselines:\n\nThanks for the suggestion. Please see the common response and the attached pdf file for details.\n\nQ1. What is the computational time of this approach compared with sparse matrix multiplication?\n\nWe believe that a fair comparison would not be comparing our proposed algorithm with sparse matrix multiplication since our algorithm includes a multiplication step after reducing the size of the graph adjacency matrix via leverage score sampling. Hence our algorithm could potentially be combined with sparse matrix multiplication. For this reason, we would expect a similar tendency as demonstrated in Table 1 in the attached pdf that even with adopting sparse matrix multiplication; that is, the computation time of our approach will be faster than the regression by using the exact multiplication of $AX$.\n\nQ2. Besides MSE, how about applying this to the real semi-supervised classification tasks? Will this lead to the same prediction accuracy?\n\nWe thank the reviewer for suggesting interesting directions for future work. Our proposed sampling scheme (based on approximating the leverage score) and the theoretical guarantee we made are tailored for the node-level regression task in the training phase. The reason we focused on this specific setting was so that we could take advantage of the theoretical results on leverage score sampling and linear regression in order to obtain theoretical results for GNN training using a subsampled graph. In particular, we are helped by the fact that in the linear GCN setting, the optimal solution for the linear regression is a matrix-vector product which can be expressed as $((AX)^{\\top}(AX))^{\\top}\\cdot (AX)^{\\top} \\mathbf{y}$, where $A$, $X$, and $\\mathbf{y}$ denote the adjacency matrix, data matrix, and the labels respectively.\n\nFor graph link prediction and classification tasks (under possibly semi-supervised settings), however, since the optimal solution cannot be expressed as a matrix-vector product (due to non-linearities), our proposed algorithm does not extend to the classification problems in a straightforward way and similarly our theoretical guarantee may not hold anymore. That being said, we agree that both node and edge classification problems are interesting. For future work, we want to study whether our same subsampling process can have any guarantee for a node classification task. For more general tasks,  we believe that we need different sampling strategies other than the leverage score sampling and hence we may need other tools for developing a theoretical guarantee for sampling algorithms. \n\nQ3. Is it possible to use this for link prediction tasks when dropping edges for validation? Is there still any theory guarantee in this case?\n\nThat is an interesting question. As mentioned above, our proposed sampling scheme is a node subsampling scheme that selects some representative nodes based on leverage scores and uses all their edges and is tailored for predicting labels (continuous value) of each node. It is not straightforward to apply our algorithm to predicting links between nodes nor to extend the theoretical results to that setting.  However, standard GNN approaches to link prediction still rely on message passing on the graph, requiring the computation of $AX$, and could benefit from a subsampling technique similar to the one we proposed in order to reduce the computational costs. But since link prediction is a fundamentally different inference task, we will need to find an alternative to the leverage score sampling-based approach that we considered. This is an interesting direction for future work.\n"}, {"review_id": "lUVcS07ooJ", "reviewer": "Reviewer_F1DR", "summary": "This paper proposes a method to overcome the computational difficulties often encountered when using graph neural networks (GNNs) for large-scale datasets, by subsampling in the adjacency and feature matrices. By assuming a two-layer linear GNN for the regression problem and performing subsampling based on the leverage score, the authors show that the MSE of the reduced problem generated by their proposed method is bounded from above by the true MSE times $(1+\\epsilon)$.\n", "questions": "- Concerning to the experimental section, how about the computational time in practice? The proposed method comes from the computational issue and hence it is better to show the computational time/cost as the result. \n\n- Figures 2b and 2d are the magnified views of 2a and 2c, respectively. Is this right? If so, it is better to mention this point.\n", "limitations": "The authors address their research limitations well. I think there is no concern about the potential societal impact. \n", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Originality & Significance: The proposed method basically follows the common idea in subsampling algorithms for linear regression in a standard setting; the same is true for the derivation of the upper bound of the reduced problem MSE. However, in the current problem, the regression coefficients are multiplied with the product of the adjacency matrix and the feature matrix, and the computation of the product itself is expensive and should be avoided. A new twist is introduced at this point: the proposed algorithm breaks the problem into two steps; in the first step the leverage score is approximately estimated by subsampling the rows of the adjacency matrix and the corresponding columns of the feature matrix; in the second step, the reduced adjacency matrix is generated according to the estimated leverage score. In this way, the direct computation of the product is avoided. I think this idea is interesting. Subsampling algorithms with theoretical guarantees are currently rare in this field, and this makes the proposed method valuable. The originality and significance of the proposed method are considered to be thus high.\n\n\nQuality: The explanation of the algorithm and the derivation of the theoretical guarantee are very clear. The numerical experiment is also convincing. The quality of the paper is thus high.\n\n\nClarity: The presentation is nice and clear, except that the references are missing in the main manuscript (they are in the supplementary material). I think the explanations of the numerical experiment result should be more detailed, though the important parts are understandable. Some typos are also found. Although there are small flaws like these, the overall clarity is high.\n\nList of typos found:\n\nLine 95: discussio -> discussion\n\nLine 186: unnecessary space seems to exist before ``Algorithm''. \n\nLine 225: unnecessary space seems to exist before ``Algorithm''. \n\nLine 226: Appendix Appendix G -> Appendix G \n\nLine 245: 4039 -> 4039 nodes\n\nLine 248: 100),in -> 100), in \n", "weaknesses": "- The linear GNN is not employed in practice, and the extension to nonlinear GNN is desired. I think it is at least possible to apply the proposed algorithm to the nonlinear case, although there is no theoretical guarantee in that case. Some more discussion or experimental results in such a nonlinear case are welcome.\n\n- The extension to the more-than-two layer case is also nontrivial. Some more discussions are welcome.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "W1. Extension to nonlinear GCN:\n\nWe thank the reviewer for suggesting interesting directions for future work. We leave comments in the common response section and provide additional experimental results for nonlinear GCN in the attached pdf.\n\nW2. Extension to the more-than-two-layer case:\n\nWe thank the reviewer for pointing out that we did not discuss this properly. Note that without non-linearities, considering a multi-layer network is equivalent to performing additional matrix multiplication by $A$ and the output of the previous layer. Since our proposed scheme yields a spectral approximation guarantee for matrix multiplication, we believe that an analysis similar to that of Algorithm 1 should be possible for this extension to obtain a $(1+\\epsilon)^{L}$-approximation guarantee, where $L$ denotes the number of layers. Since the approximation guarantee comes with high probability, a union bound could be used to obtain a high probability guarantee for the multi-layer case (since the order or error for each approximate matrix multiplication is the same). If the paper is accepted, we will include a discussion on this extension.\n\nQ1. Concerning to the experimental section, how about the computational time in practice? The proposed method comes from the computational issue and hence it is better to show the computational time/cost as the result.\n\nPlease see the common response and additional run-time experiments described in the attached pdf.\n\nQ2. Figures 2b and 2d are the magnified views of 2a and 2c, respectively. Is this right? If so, it is better to mention this point.\n\nYes, that is correct. We will add a sentence to mention that point.\n\n\n\n"}, {"review_id": "VEhAOLmtD0", "reviewer": "Reviewer_XHGm", "summary": "The work tackles the problem of scaling GNNs used for regression tasks to large graphs by subsampling given nodes. The technique consists of first performing node sub-sampling to estimate the leverage scores of $AX$ and then performs leverage score sampling on $AX$. The authors show that this technique is a good spectral approximation to the full $AX$ computation, but avoids the $O(n^2d)$ matrix multiplication cost. ", "questions": "While the work mentions linear GNNs, it seems to be the case that the work only considers one-hop linear GNNs. While a multi-layer linear GNN is still a linear one, do these results apply to multi-hop linear GNNs as well? It would be useful to have more details on the number of layers and exact model used in section 5 as well.\n\nIs there a reason in figure 2 (b) and (d) that MSE with uniformly sampled $AX$ is not included in the experiments? \n\n", "limitations": "The authors acknowledge a main limitation that is the fact that this work only studies linear GNNs. A further limitation is that while the authors acknowledge the improved theoretical run-times of their technique, they do not give practical wall-clock runtimes for the experiments. While the technique may be theoretically faster, in practice many operations such as matrix operation are extremely optimized. ", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper tackles the interesting and practically useful problem of scaling GNNs to large graphs. The proposed technique has useful theoretical guarantees and as far as I am aware it is a novel approach. The approximation to the regression problem is interesting as it has useful spectral guarantees and could be of general interest from an algorithmic perspective outside of the GNN community.   ", "weaknesses": "There are a few weaknesses of the work, some of which the authors point out. The most glaring one is the fact that this work targets only linear GNNs. While analysing non-linear models is clearly more challenging, it would be interesting to have some comments or results for non-linear models with more easy to study non-linearities (such as ReLU). \n\nI found following the different run-times a bit challenging. It may be useful to have a summary table that compares the various computational/memory complexities for the different techniques used in the experiments in Section 5. \n\nAs this problem is a practical one, I believe Section 5 could be stronger. While the authors argue that the matrix multiplication is a bottleneck, modern GPUs are extremely optimized for exactly these sort of operations and as a consequence it may be the case that the full matrix multiplication is still faster than sampling. As such, I would be very interested to see some results on the runtimes of the various techniques. Further, it would be interesting to test regimes in which the graph is too large to fit in a GPU and sampling is absolutely necessary. \n\nIt may also be useful to run experiments with non-linear GNNs using the proposed sampling technique.\n\nI also spotted some typos and formatting issues in the manuscript (line 60, line 95). ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "W1 & W2. Lack of run-time comparison:\n\nThank you for pointing out this issue. Please see the common response and the additional experiments in the attached pdf above.\n\nW3. Extension to nonlinear GCN:\n\nPlease see the common response.\n\nW4. Some typos (line 60, line 95):\n\nThanks for pointing out these typos. We fixed those.\n\nQ1. Extension to the more-than-two-layer case:\n\nWe thank the reviewer for pointing out that we did not discuss this properly. Note that without non-linearities, considering a multi-layer network is equivalent to performing additional matrix multiplication by $A$ and the output of the previous layer. Since our proposed scheme yields a spectral approximation guarantee for matrix multiplication, we believe that an analysis similar to that of Algorithm 1 should be possible for this extension to obtain a $(1+\\epsilon)^{L}$-approximation guarantee, where $L$ denotes the number of layers. Since the approximation guarantee comes with high probability, a union bound could be used to obtain a high probability guarantee for the multi-layer case (since the order or error for each approximate matrix multiplication is the same). If the paper is accepted, we will include a discussion on this extension.\n\nQ2. Is there a reason in Figures 2 (b) and (d) that MSE with uniformly sampled AX is not included in the experiments?\n\nThank you for pointing out the need for clarification. Figures 2b and 2d are the magnified views of 2a and 2c and hence the uniformly sampled curve does not appear within the axes limits, respectively.\n\nL1. The authors acknowledge a main limitation which is the fact that this work only studies linear GNNs. A further limitation is that while the authors acknowledge the improved theoretical run-times of their technique, they do not give practical wall-clock runtimes for the experiments. While the technique may be theoretically faster, in practice many operations such as matrix operation are extremely optimized.\n\nWe agree that the focus on linear GNNs is a limitation of this work, but we view our work as studying a specific regime where existing theoretical tools allow the derivation of theoretical guarantees. We also agree that practical wall-clock run-times for the experiments are important here. We have run additional experiments to provide wall-clock runtimes for the experiments and also empirical results on nonlinear GNNs. Please see the common response above and the attached pdf file for detailed experiment results.\n\n\n"}]

Author rebuttal:
We thank the reviewers for their valuable comments. For weaknesses/questions raised by multiple reviewers, we provide common responses. For the remaining comments, we provide point-by-point responses below.

We would like to first highlight what we see as the main contribution of our paper. Our paper provides an efficient sampling algorithm for learning linear GNNs with a large number of nodes with a theoretical guarantee. Specifically, for a two-layer linear GNN with adjacency matrix $A \in \mathbb{R}^{n \times n}$ and data matrix $X \in \mathbb{R}^{n \times d}$, we show that it is possible to learn the model accurately while observing only $O(nd \log(n))$ entries of $A$. This in turn yields a run-time of $O(nd^2\log n)$, avoiding the computation time $O(n^2d)$ of using the full matrix multiplication $AX$. This reduction in run-time complexity is particularly significant when the number of nodes $n$ is large and the feature dimension is much smaller compared to $n$ (i.e., $d \ll n$).
We would also like to highlight that this result holds even when the graph adjacency matrix and the data matrix are dense in which case sparse matrix multiplication would not help much. Hence the proposed algorithm bears practical importance as well.
While we view our main contribution to be of a theoretical flavor, we agree with several of the reviewers' comments requesting additional experimental validation, which we discuss next.

W1. Lack of run-time comparison:

We note that the experiments in the original submission focused on the tradeoff between the number of observed entries and regression accuracy. In the attached pdf, we provide the computation time comparison for the end-to-end training process for a regression task. In particular, we compared the wall-clock time of performing a regression task with full $AX$ computation and that with Algorithm 1 (our proposed algorithm) that uses partial observations of $A$ and $X$. For a fair comparison, we use the same regression solver and the use the same specification of 48 cores of a x86 64 processor with 503.74GB memory. The results  (see Table 1 in attached pdf) show that our proposed scheme requires orders of magnitude less wall-clock time when the graph is large than the regression with exact $AX$ computation. In particular, for the “ogbn-arxiv dataset” (from Open Graph Benchmark) having 17K nodes and 1M edges,  our algorithm runs about 40x faster than the regression with the exact computation of $AX$.

In addition to run-time, we expect our proposed approach to offer significant improvements in terms of memory usage. As a simple experiment in terms of memory usage, in Table 2 in the attached pdf,  we provide the peak memory usage comparison for the same datasets during the sampling and regression steps of the algorithm. In the best case, the proposed algorithm requires 1414× less memory than the algorithm using the exact computation of $AX$ in the best case.

W2. Comparison with other baselines (e.g., GraphSage and GraphSaint) and the pros and cons:

Thanks for the suggestion. In the original submission, we focused on comparisons with two baselines (full computation of $AX$ and reduced computation via uniform sampling) to validate and highlight the conceptual contributions of our approach. But we agree that comparing our approach to popular methods such as GraphSage and GraphSaint is helpful. The pros and cons between our algorithm and other baselines (GraphSage and GraphSaint) are summarized below:
 
1. Both our algorithm and the baselines use some form of subsampling of graphs motivated by the large computational/storage complexity of running GNNs with large-scale graphs. However, our algorithm comes with a theoretical guarantee of getting $(1+\epsilon)$-approximation guarantee to the optimal MSE (for the specific setting of a two-layer linear GCN), while only observing $O(nd \log n)$ entries of $A$. GraphSage and GraphSaint, to the best of our knowledge, do not provide similar theoretical guarantees.

2. To obtain the performance guarantee, our algorithm employs a more complicated sampling strategy based on leverage score sampling, which leads to a non-uniform subsampling of the graph. GraphSage uses uniform sampling for choosing the neighborhood of a node when performing a feature aggregation.

Furthermore, we demonstrate comparisons between our sampling approach and the sampling strategies utilized in (1) GraphSage and (2) GraphSaint (not the tools themselves), in addition to applying regression solvers. Figure 2 in the attached file demonstrates the MSE of our proposed algorithm, the baselines we have considered, and GraphSage and GraphSaint. Figures 2(a) and 2(b) demonstrate that the sampling ideas from GraphSage and GraphSaint do not seem to outperform our proposed sampling technique (tailored for the graph-weighted linear regression).

W3. Extension to nonlinear GCN:

We thank reviewers for suggesting interesting directions for future work. As illustrated in Figure 2(d) in the attached pdf, we have run additional experiments to demonstrate our proposed algorithm's efficacy on nonlinear GCNs. Using ReLU activation function and one hidden layer, Figure 2(d) plots the mean squared error as a function of the observation budget, shown as a percentage of the number of observed nodes in the graph. Interestingly, we observe that running the non-linear GCN with our proposed sampling techniques have a similar performance to the one with exact $(AX, \mathbf{y})$, as the budget increases. Also, the error reduction on Algorithm 1 from Regression with uniformly sampled $(AX, \mathbf{y})$ is 60% when the budget is at 5%.
Besides, as compared to Figure 2(c), we observe that although our proposed algorithm is designed to optimize the MSE on linear GCNs, the performance works even better for the non-linear case. It would be interesting to see whether a slight modification of our node subsampling technique will lead us to prove a theoretical guarantee.

Meta-review: This paper proposes an efficient sampling method for linear Graph Neural Networks. Overall, this is a well-written theoretical paper that proposes a new algorithm of node subsampling for GNNs and evaluates the reasonable upper bound of the approximation error. The analyzed model and task are limited to a linear model and a regression for enabling theoretical analysis.

This paper has both positive and negative reviews, and thus this is a borderline paper. There is some limitation (it is a linear model), however, the idea itself is interesting and theoretically well-supported. Moreover, the reviewer who agrees to accept is excited about the paper. One of the weaknesses of this paper is that it does not compare with other nonlinear GNN models, as pointed out by a reviewer. If this paper proposes a nonlinear GNN method, this could be a reason for rejection. However, this paper considers a linear setting; it sounds OK not to compare with nonlinear models.

Based on the above reasons, I vote for acceptance.
Decision: Accept (poster)