On the Importance of Sampling in Training GCNs: Convergence Analysis and Variance Reduction
Authors: Weilin Cong
Morteza Ramezani
Mehrdad Mahdavi
Venue:  (2021)
Keywords: Graph neural network
large-scale machine learning
convergence analysis

Abstract: Graph Convolutional Networks (GCNs) have achieved impressive empirical advancement across a wide variety of graph-related applications. Despite their great success, training GCNs on large graphs suffers from computational and memory issues. A potential path to circumvent these obstacles is sampling-based methods, where at each layer a subset of nodes is sampled. Although recent studies have empirically demonstrated the effectiveness of sampling-based methods, these works lack theoretical convergence guarantees under realistic settings and cannot fully leverage the information of evolving parameters during optimization. In this paper, we describe and analyze a general \textbf{\textit{doubly variance reduction}} schema that can accelerate any sampling method under the memory budget. The motivating impetus for the proposed schema is a careful analysis for the variance of sampling methods where it is shown that the induced variance can be decomposed into node embedding approximation variance (\emph{zeroth-order variance}) during forward propagation and layerwise-gradient variance (\emph{first-order variance}) during backward propagation. We theoretically analyze the convergence of the proposed schema and show that it enjoys an $\mathcal{O}(1/T)$ convergence rate.  We complement our theoretical results by integrating the proposed schema in different sampling methods and applying them to different large real-world graphs.

Under review as a conference paper at ICLR 2021
ON THE IMPORTANCE OF SAMPLING IN TRAINING GCNS: CONVERGENCE ANALYSIS AND VARIANCE REDUCTION
Anonymous authors Paper under double-blind review
ABSTRACT
Graph Convolutional Networks (GCNs) have achieved impressive empirical ad- vancement across wide variety of graph-related applications. Despite their great success, training GCNs on large graphs suffers from computational and memory issues. A potential path to circumvent these obstacles is sampling-basee methods, where at each layer a subset of nodes is sampled. Although recent studies have em- pirically demonstrated the effectiveness of sampling- based methods, these works lack theoretical convergence guarantees under realistic settings and cannot fully leverage the information of evolving parameters during optimization. In this pa- per, we describe and analyze a general doubly variance reduction schema that can accelerate any sampling method under the memory budget. The motivating impetus for the proposed schema is careful analysis for the variance of sampling methods where is shown that the induced variance can be decomposed into node embedding approximation variance (zeroth-order variance) during forward prop- agation and layerwise- gradienn variance (firstorder variance) during backward propagation. We theoretically analyze the convergence of the proposed schema and show that it enjoys an OL/TT) convergence rate. We complement our theo- retical results by integrating the proposed schema in different sampling methods and applying them to different large real-world graphs
INTRODUCTION
In the past few years, graph convolutional networks (GCNs) have achieved great success in many graph-related applications, such as semi-supervised node classification (Kipf & Welling, 2016), su- pervised graph classification (Xu et al., 2018), protein interface prediction (Fout et al., 2017), and knowledge graph (Schlichtkrull et al., 2018; Wang et al., 20.7). However, most works on GCNs fo- on relatively small graphs, and scaling GCNs for large-scale graphs is not straight forward. Due to the dependency of the nodes in the graph, we need to consider a large receptive-fiele to calculate the representation of each node in the mini-batch, while the receptive field grows exponentially with respect to the number of layers. To alleviate this issue, sampling-based methods, such as node-wise sampling (Hamilton al., 2017: Ying al., 2018; Chen et al., 2017), layer-wise sampling (Chen et al., 2018; Zou et al., 2019), and subgraph sampling (Chiang et al., 2019; Zeng et al., 2019) are proposed for mini-batch GCN training.
Although empirical results show that sampling- -based methods can scale GCN training to large graphs, these methods suffer from few key issues. First, the theoretical understanding of sampling- based methods is still lacking. Second, the aforementioned sampling strategies are only based on the structure of the graph. Although most recent works (Huang et al., 2018; Cong et al., 2020) propose to utilize adaptive importance sampling strategies to constantly re-evaluate the relative importance of nodes during training (c.g., current gradient or representation of nodes), finding the optimal adap- tive sampling distribution is computationally inadmissible, it requires to calculate the full gradient or node representations in each iteration. This necessitates developing alternative solutions that can efficiently be computed and that come with theoretical guarantees.
In this paper, we develop a novel variance reduction schema that can be applied to any sampling strategy to significantly reduce the induced variance. The key idea is to the historical node

---

Under review as a conference paper at ICLR 2021
LADIES
LADIES (Fint-onder
LADIES (Zeniih-odder VR)
DDoully VR)
Trinnng Loss
Viliaaion Loas
of Crndient
100
50
150
1.50
0.0
150
Epoch
Epach
Epach
Figure 1: The effect of doubly variance reduction on training loss, validation loss, and mean-square error (MSE) of gradient on Flickr dataset using LADIES proposed in Zou et al. (2019).
embeddings and the historical layerwise gradient of each graph convolution layer as control vari- ants. The main motivation behind the proposed schema stems from our theoretical analysis of the sampling methods' variance in training GCNs. Specifically, we show that due to the composite struc- ture of training objective, any sampling strategy introduces two types of variance in estimating the stochastic gradients: node embedding approximation variance (zeroth-order variance) which results from embeddings approximation during forward propagation, and layerwise-gradient variance (first- order variance) which results from gradient estimation during backward propagation. In Figure 1, we exhibit the performance of proposed schema when utilized in the sampling strategy introduced in (Zou et al., 2019). The plots show that applying our proposal can lead to significant reduction variance; hence faster convergence rate and better test accuracy. We can also see that both zeroth- order and first-order methods are equally important and demonstrate significant improvement when applied jointly (i.e, doubly variance reduction).
Contributions. We summarize the contributions of this paper as follows:
We provide the theoretical analysis for sampling-baset GCN training (SGCN) with a non- asymptotic convergence rate. We show that due to the node embedding approximation variance, SGCNs suffer from residual error that hinders their convergence.
We mathematically show that the aforementioned residual error can be resolved by em- ploying zeroth-order variance reduction to node embedding approximation (dubbed as SGCN+), which explains why VRGCN (Chen al., 2017) enjoys better convergence than GraphSAGE (Hamilton et al., 2017), even with less sampled neighbors.
We extend the algorithm from node embedding approximation to stochastic gradient approximation, and propose generic and efficient doubly variance reduction schema (SGCN++). SGCN++ can be integrated with different sampling-bas baaedd methods to signif- icantly reduce both zeroth- and first-order variance, and resulting in faster convergence rate and better generalization.
. We theoretically analyze the convergence of SGCN++ and obtain an O(1/T) rate, which significantly improves the best known bound O(1/VT). We empirically verify + through various experiments on several real-world datasets and different sampling methods, where it demonstrates significant improvements over the original sampling methods.
RELATED WORKS
Training GCNs via sampling. The full-batch training of a typical GCN is employed in Kipf & Welling (2016) which necessities keeping the whole graph data and intermediate nodes' represen- tations in the memory. This is the key bottleneck that hinders the scalability of full-batch GCN training. To overcome this issues, sampling-bai
GCN training methods (Hamilton et al., 2017: Chen et al., 2017; Chiang et al., 2019; Chen et al., 2018; Huang et al., 2018) are proposed to train GCNs based on mini-batch of nodes, and only aggregate the embeddings of a sampled subset of neighbors of nodes in the mini- batch For example, GraphSAGE (Hamilton et al., 2017) restricts the computation complexity by uniformly sampling a fixed number of neighbors from the previous layer nodes. However, significant computational overhead is introduced when GCN goes deep. VRGCN (Chen et al., 2017) further reduces the neighborhood size and uses history activation of the previous layer to reduce variance. However, they require to perform a full-batch graph convolu- tional operation on history activation during each forward propagation, which is computationally expensive. Another direction applies layerwise importance sampling to reduce variance. For exam- plc, FastGCN (Chen et al., 2018) independently sample a constant number of nodes in all layers using importance sampling. However, the sampled nodes too sparse to achieve high accuracy.
2

---

Under review as a conference paper at ICLR 2021
LADIES (Zou et al., 2019) further restrict the candidate nodes in the union of the neighborhoods of the sampled nodes in the upper layer. However, significant overhead may be incurred due to the ex- pensive sampling algorithm. In addition, subgraph sampling methods such GraphSAINT (Zeng et al., 2019) constructs mini-batches by importance sampling, and apply normalization techniques to eliminate bias and reduce variance. However, the sampled subgraphs are usually sparse and requires a large sampling size to guarantee the performance.
Theoretical analysis. Despite many algorithmic progresses over the years, the theoretical under- standing of the convergence for SGCNs training method is still limited. VRGCN provides a conver- gence analysis under strong assumption that the stochastic gradient due to sampling is unbiased and achieved convergence rate of (L/VVT). However, the convergence analysis is limited to VRGCN, and the assumption is not true due to the composite structure of training objective as will be elaborated. Chen & Luss (2018) provides another convergence analysis for FastGCN under a strong assumption that the stochastic gradient of GCN converges to the consistent gradient exponen- tially fast with respect to the sample size, and results the same convergence rate as unbiased ones, i.e. ((/vT)). Most recently, Sato et al. (2020) provides PAC learning-style bounds on the node embedding and gradient estimation for SGCNs training. Another direction of theoretical research focuses on analyzing the expressive power of GCN (Garg et al., 2020; Chen et al., 2019; Zhang et al., 2020), which is not the focus of this paper and omitted for brevity.
Connection to composite optimization. The proposed doubly variance reduction algorithm shares the same spirit with the variance reduced composite optimization problem considered in Zhang & Xiao (2019a); Hu et al. (2020); Tran-Dinh et al. (2020); Zhang & Xiao (2019c;b), but with two main differences. Firstly, the objective function is different. In composite optimization, objective function that only the first composite layer has the trainable parameters is considered, where the output of the lower-level function acting as the parameter for the higher-level function. However, in GCN model, the output of one graph convolutional layer is the input node embedding matrix for the next layer. As result, when analyzing the convergence of variance reduced neural network model, we have to explicitly handle both the evolving node embedding matrices and the trainable parameters at all layers. Secondly, the data points in these works are sampled independently. but the data points (nodes) in SGCN are sampled node- or layer-dependent according to the graph structure. In our analysis, we provide a sampled -gaaph structure dependent convergence rate by bridging the connection between the convergence rate of GCN to the graph Laplacian matrices.
3 SGCN: A TIGHIT ANALYSIS OF SGD FOR GCN TRAINING
Full-batch GCN training. We begin by introducing the basic mathematical formulation of training GCNs. In this paper, we consider training GCNs in semi-supervised multi-class classification set- ting. Given an undirected graph G = (V. E) with N | and |&| edges and the adjacency matrix A E {0, 1 we assume that each node is associated with feature vector æ; E I and label yi- We use X [x...... xN] E R Nxd and y = ly uN E IRM to denote the node feature matrix and label vector, respectively. The Laplacian matrix calculated as L = 112 AD A or 1/2
D---- 'A where D E NxN is the degree matrix. We use 0 (W(1) ...., the stacked weight parameters of a L-layer GCN. The training of full-batch GCN (Ful (illlGM) as an empirical risk minimization problem aims at minimizing the loss L(0) over all training data
,..., w(1) denote
= 2 EN Losshh (L) .)),, .8),,, H = =o( (L..... (Lo (LXW() w(O) ..()))
((0)=
Z
where h (c) is the ith row of node embedding matrix H(O = o oZ(0), Z = L that corresponds to embedding of ith node at 9th layer (hop), Loss(-,)) -) is the loss function (e.g.. cross- entropy loss) to measure the discrepancy between the prediction of the GCN and its ground truth label, and (-) is the activation function (e.g., ReLU function).
campling-based GCN training. When the graph is large, the computational complexity of forward and backward propagation could be very high. One practical solution to alleviate this issue is to sample a subset of nodes and construct a sparser normalized Laplacian matrix ((0) for each layer with (!)) < supp(L), and perform forward and backward propagation only based on the sampled Laplacian matrices. The sparse Laplacian matrix construction algorithms can be roughly
3

---

Under review as a conference paper at ICLR 2021
classified nodewise sampling, layerwise sampling, and subgraph sampling. A detailed discussion of different sampling strategies can be found in Appendix D. To apply Stochastic Gradient Descent (SGD) (Bottou et al., 2018) to train GCN (SGCN), all nodes with size B | and construct the set of sparser Laplacian matrices {L ) J based on nodes sampled each layer and compute the stochastic gradient update parameters
sample a mini-batch of nodes VB S V from
vū(0) B Liev, (Loss(h; vs). H = a (I .... (L(3) L w(e) ...W))
20
Key challenges. Compared to vanilla SGD, the key challenge of theoretical understanding for SGCN training is the biasedness of stochastic gradient due sampling of nodes at inner layers. Let ac(0) denote FullGCN's full-batch gradient as VL(() = {G = aw!) V/=1 and SGCN's stochastic
gradient as Vā(() {G(0 =
Swan Ve. By the chain rule, we can compute the full-batch
gradient (()) w.r.t the éth layer weight matrix (( as
G = |LH 11" (D ( o o )z))))), D = L' (D)
(D (+))) c oz w(O, D
2210. 0H12)
(1)
and compute stochastic gradient G utilized in SGCN for the fth layer w.r.t W( as - " "D o o
)).") D [I (1,T((;;*1) oo W/) B
aã1e. oH
(2)
For any layer & E [L]. the stochastic gradient G.!² is a biased estimator of full-batch gradient ((
as it is computed from H.()) and z which are not available in SGCN since Hi and Z are used as an approximation during training. Recently, Chen et al. (2017) established a convergence rate under the strong assumption that the stochastic gradient of SGCN is unbiased and Chen & Luss (2018) provided another analysis under the strong assumption that the stochastic gradient converges to the consistent gradient exponentially fast as the number of sampled nodes increases. While both studies establish the same convergence rate of O(1/VT) however, these assumptions do not hold in reality due to the composite structure of the training objectives and sampling of nodes at inner layers. Motivated by this, we aim at providing tight analysis without the aforementioned strong assumptions on the stochastic gradient. Our analysis is inspired by the bias and variance decompo- sition of the mean--quare error of stochastic gradient, which has been previously used in Cong et al. (2020) to analysis the stochastic gradient in GCN. Formally, we can decompose mean-square error of stochastic gradient as
E||l E VC) B| - -
(3)
Bius [
Variance ||||D||]]
where the bias terms TE b is mainly due to the node embedding approximation variance (zeroth- order variance) during forward propagation and the variance term ||[||| is mainly due to the layerwise gradient variance (first-order variance) during backward propagation. Before proceeding analysis, we make the following standard assumptions on the Lipschitz-continue and smoothness of the loss function Loss(,,) :) and activation function (.)).
Assumption 1. The loss function Loss(;))) .) is the input node embedding vector, i.e., ()))ss(H(()) .y) Loss(h" ( .v)ll2 < Closs (1))
continuous and Lions tmos-smottsesss w.r.f. h (T) ||2
(l
and || Losss h (L), (,)) - LLoss(h" (L) .3))l2 S Lious L - h (1)) 1|2-
Assumption 2. The activation function o()) is Co-Lipschitz continuous and Lo to-mmothhnsss, i.e. [|a()(()) - -z (2) < C= - 2x) [|2 and (|a'(=(()) - -=(=((()) S Loll = 2000
We also make the following assumptions on the norm of weight matrices, Laplacian matrices, and node feature matrix, which are used in the generalization analysis of GNNs Garg al. (2020).
Assumption 3. For any l E [L], the norm of weight matrices, Laplacian matrices, input node feature matrix are bounded: ( 1pp S Bw, [ I BLA, [L]| I BLA and |X|| S Bh.
We use tilde symhol Ã for their stochastic form
4

---

Under review as a conference paper at ICLR 2021
Proposition 1. For any E [L], there exist constants Bи and Bp such that the norm of node embedding matrices and the gradient with respect to the input node embedding matrices satisfy
|H S Bи. |H S Bu, Bo 1WO))
1 BD. and (o H 1w(()) 11 BD-
OH
8H
Before presenting the convergence of SGCN, we introduce the notation of propagation matrices which are defined as the column-wise expectation of the sparser Laplacian matrices. Note that this notation is only for presenting the theoretical results, and are not used the practical training algorithms. By doing so, we can decompose the difference between ( and L as the summation of column-wise difference [[]( p()) l and row-wise difference
Iē=1: (P(O)
(P()) -
In the following theorem, we show that the upper bound of the bias and variance of stochastic gradient is closely related to the expectation of column-wise difference E[ p and row-
wise difference B(][P(() || which can significantly impact the convergence of SGCN.
Theorem (Convergence of SGCN), Suppose Assumptions I, 2, 3 hold and apply SGCN with learn- ing rate chosen as 1 = min Lp, 1/vT) where LE is the smoothness constant. Let An and Ab. denote the upper bound on the variance and bias of stochastic gradients as:
1. Σ
A
OE
+ O 1(0
LIED,
Ab = -
(4)
Then, the output of SGCN satisfies
2(((0)) - 2(0)))
Lp VT
min BIIVCC))1 EETT
+ Ab.
(5)
VT
The exact value of key parameters LF. An. and An are computed in Lemma 1, Lemma 2, and Lemma 3 respectively and can be found in Appendix G. Theorem implies that after T iterations the gradient norm of SGCN is most ((An/VT) + Åb. which suffers from constant residual error Ab that is not decreasing as the number of iterations' I' increases. Without the bias" we recover the convergence of vanilla SGD. Of course, this type of convergence is only useful if Ab and A, are small enough. We note that existing SGCN algorithms propose to reduce Ab by increasing the number of neighbors sampled at each layer (e.g., GraphSAGE), or applying importance sampling (e.g., FastGCN, LADIES and GraphSAINT).
4 S ZEROTH-ORDED VARIANCE REDUCTION
An important question to answer is: can we eliminate the residual ervor without using all neigh- bors during forwurd-prepagatton
A remarkable attempt to answer this question has been recently made in VRGCN (Chen et al., 2017) where they propose to use historical node embeddings as an approximation estimate the true node embeddings. More specifically, the graph convolution in VRGCN is defined H. a 1-1 of historical node embeddings, VRGCN requires less sampled neighbors and results in significant less computation overhead during gradient computation. Although VRGCN achieves significant speed up and better performance compared other SGCNs, it involves using the full Laplacian matrix at each iteration, which can be computationally prohibitive. Moreover, since both SGCNs and VRGCN are approximating the exact node embeddings calculated using all neighbors, it is still not clear why VRGCN achieves a better convergence result than SGCNs using historical node embeddings.
i i (( -H "(""))((()) Taking advantage
To fill in these gaps, introduce zeroth-order variance reduced sampling based GCN training method dubbed as SGCN+. As shown in Algorithm 1, SGCN+ has two types of forward propagation: the forward propagation at the snapshot steps and the forward propagation the regular steps. At the snapshot step (t. mod K = 0), full Laplacian matrix is utilized:
Z = LH (8-1) (!
! = oz 2, 4 Z
(6)
During the regular steps (l mod K + 0), the sampled Laplacian matrix is utilized:
7
2 - i 8 1 w( 1 "( (-11
H (Z
We have Ab 0 if all neighhor are used calculate the exact node embeddings, i.e., p()) L, V8 [7].
5

---

Under review as a conference paper at ICLR 2021
Algorithm S Zeroth-order variance reduction (Detailed version in Algorithm 4)
1:
Input: Learning 1) 0. snapshot gap K
for t= I...., do
if mod K then
4:
Calculate node embeddings using Eq. 6
5:
Calculate full-batch gradient V ((())) as Eq. and update 0.11 0: -VV((())
6:
else
Calculate node embeddings using Eq. 7
8:
Calculate stochastic gradient VC(0.) Eq. and update 0.+1 0. -vVV(()))
end if
10:
end for
11:
Output: Model with parameter 0y+1
Algorithm ++ Doubly variance reduction (Detailed version Algorithm 5)
I: Input: Learning rate > 0, snapshot gap K >> 0
2: for t ., do
3: if mod K then
Calculate node embeddings using Eq. 6
5:
Calculate full-batch gradient V (0,) usng Eq. and update as 0,-1 = 0, -7)((())))
Save the per layerwise gradient (6) +- ( ( + D VEE [4]
6:
7:
else
8:
Calculate node embeddings using Eq. 7
9:
Calculate stochastic gradient VZ(),) usng Eq. 10 and update as 0,41 = 0. n)V2(0))
10:
end if
11:
end for
12:
Output: Model with parameter 0r+1
Comparing with VRGCN, the proposed SGCN+ only requires one full Laplacian graph convolution operation every K iterations, where K is an additional parameter to be tuned.
In the following theorem, we introduce the convergence result of SGCN+. Recall that the node em- bedding approximation variance (zeroth-order variance) determines the bias of stochastic gradient B|||b|| []] Applying SGCN+ can significantly reduce the bias of stochastic gradient, such that its value is small enough that will not deteriorate the convergence.
Theorem 2 (Convergence of SGCN+), Suppose Assumptions 1, 2, 3 hold and apply SGCN+ with learning rate chosen as mi 1/VT) where Lp is the smoothness constant. Let An and A denote the upper bound for the variance and bias of stochastic gradient as:
An O - +(0
P -LIE)),
(8)
A
A where A = o(K Σ
B O I
Then, the oulput of SGCN+ satisfies
min B(I 2(C(0)) - C(())) NETT]
LpAn
A T
+
VT
The exact value of key parameters Lp, An. and A5 are computed in Lemma 1, Lemma 2, and Lemma respectively, and can be found in Appendix G, H. Theorem implies that after iterations the gradient norm of SGCN+ is at most O(A,/vT) + ((At/T)) When using all neighbors for calculating the exact node embeddings, we have p(t) = L such that A O, which leads to con- vergence rate of SGD. Comparing with vallina SCCN, the bias of SGCN- is scaled by learning rate reduce the negative effect of bias by choose learning rate
OI)vTT). This also explains why SGCN+ achieves a significantly better convergence rate compared to SGCN.
77. Therefore,

---

Under review as a conference paper at ICLR 2021
5 SGCN++: DOUBLY VARIANCE REDUCTION
Algorithm applies zeroth-order variance reduction on node embedding matrices and results in a faster convergence. However, both SGCN and SGCN+ suffer from the same stochastic gradient vari- ance An. which can be only reduced either by increasing the mini-batch size of SGCN or applying variance reduction on stochastic gradient. An interesting question that arises is: can we further accelerate the convergence by simultaneously employing zeroth-order variance reduction on node embeddings and first-oder variance reduction on layerwise grudient? To answer this question, we propose doubly variance reduction algorithm SGCN++, that extends the variance reduction algo- rithm from node embedding approximation to layerwise gradient estimation.
As shown in Algorithm 2, the main idea of SGCN++ is to use the historical gradient as control variants for current layerwise gradient estimation. More specifically, similar to SGCN+ + that has two types of forward propagation steps, SGCN+ also has two types of backward propagation: at the snapshot steps and at the regular steps. The snapshot steps (f mod K = 0) backward propagation are full-batch gradient computation is defined in Eq. 1, and the computed full-batch gradient are saved as control variants for the following regular steps. The backward propagation (t. mod K * 0) at the regular steps defined as
G = 6 HH ( (DietD) oZ)(2)) - [L [D o oo o((2i-1) 5 = (0 L0O1T (D.ttu) o o' [w!"] [I (D o' o'(Zz-)) 1
(10)
Next, in the following theorem, we establish the convergence rate of SGCN++. Recall that the mean-square error of the stochastic gradient can be decomposed into bias E[[| that is due to node embedding approximation and variance [[|[]|]]]]] that is due to layerwise gradient estimation. Applying doubly variance reduction on node embedding and layerwise gradient simultaneously can significantly reduce mean- meen-square error of stochastic gradient and speed up convergence.
Theorem 3 (Convergence of SGCN++). Suppose Assumptions 1, 2, 3 hold, and denote LE as the smoothness constant and A+t. as the upper-bound of mean-square error of stochastic gradient
= n+b = Ei
[B[]]
1( - |L|S|
(11)
Apply + + in Algorithm 2 with learning rale as
Then it holds that
Let 11 +40
+ Σ EINCCONNI (Lv L3 + 44 (C(O)) (C(01)-C(0') C(8"))
(12)
The exact value of key parameter LP and A n4b are computed in Lemma and Lemma 12 re- spectively, and can be found in Appendix G, Theorem 2 implies that applying doubly variance reduction can scale the mean-square error O11² K) times smaller. As result, after T' iterations the norm of gradient of solution obtained by SGCN+ is at most ((A+t//)) which enjoys the same rate vanilla variance reduced SGD (Reddi al., 2016; Fang et al., 2018).
Scalability of SGCN++. One might doubt whether the computation at the snapshot step with full- batch gradient will hinder the scalability of SGCN+ for extremely large graphs? Heuristically, we can approximate the full-batch gradient with the gradient calculated on a large-batch using all neighbors. The intuition stems from matrix Bernstein inequality (Gross, 2011), 1), where the probabil- ity of the approximation eтror violating the desired accuracy decreases exponentially as the number of samples increase. Please refer Algorithm 6 for the full-batch free SGCN++ and explanation on why large-batch approximation is feasible using tools from matrix concentration. Moreover, we provide the empirical evaluation on the large-batch size instead of full-batch in Appendix We remark that large-batch batch approximation can also be utilized in SGCN+ to further reduce the memory requirement for historical node embeddings.
Connection to composite optimization. Although we formulate sampling-based GCNs spe- cial case of the composite optimization problem. is worth noting that compared to the classical composite optimization, there are few key differences that make the utilization of variance reduc- tion methods for composite optimization non-trive triviali (a) different objective function that makes the GCN analysis challenging: (b) different gradient computation, analysis, and algorithm which make
7

---

Under review as a conference paper at ICLR 2021
Table 1: Comparison of the accuracy (F1-score) of SGCN, SGCN+, and 5GCN+
METHOD / VR
PPI
PPI-LARGE
FLICKR
REDDIT
YELP
SGCN
77.61
78.10
52.30
95.07
59.99
EXACT
SGCN+ (DOUBLY)
82.18
88.98
52.88
95.17
62.09
SGCN+ (ZEROTH)
77.65
77.82
52.57
95.17
61.29
VRGCN
SGCN+ (DOUBLY)
82.50
88.65
52.53
95.17
62.64
SGCN
71.40
69.51
51.13
94.73
59.58
GRAPHSAGE
SGCN+ (ZEROTH)
72.16
69.67
51.13
94.89
58.62
SGCN++ (DOUBLY)
79.63
85.41
52.65
95.18
61.75
SGCN
63.51
59.60
50.74
87.36
55.75
FASTGCN
SGCN+ (ZEROTH)
72.32
72.30
51.07
94.54
56.86
SGCN++ (DOUBLY)
81.92
85.04
52.57
94.99
60.63
SGCN
62.46
60.74
50.29
94.11
59.84
LADIES
SGCN+ (ZEROTH)
72.16
74.03
51.83
94.39
57.06
SGCN++ (DOUBLY)
81.58
84.18
52.09
95.05
60.96
SGCN
61.51
38.68
50.10
93.68
54.65
GRAPHSAINT
SGCN+ (ZEROTH)
66.94
41.40
50.66
84.61
55.42
SGCN++ (DOUBLY)
79.63
79.71
50.94
94.18
57.00
FULLGCN
N/A
82.14
52.99
95.15
62.77
directly applying multi-level variance reduction methods such as SPIDER (Zhang & Xiao, 2019b)
nontrivial; (c) different theoretical results and novel intuition for sampling- baased GCN training. Due to the space limit, we defer the detail discussion to the Appendix A.
EXPERIMENTS
Experimental setup. We evaluate our proposed methods in semi-supervised learning setting on various classification datasets, summarized in Appendíx B. In addition different sampling mech- anisms, we introduce Exact sampling that takes all neighbors for node embedding computation during mini-batch training (no zeroth-order variance), which can used to explicitly test the im- portance of first-order variance reduction. We add SGCN+ (Zeroth) and SGCN+ (Doubly) on top of each sampling method illustrate how zeroth-order and doubly variance reduction affect GCN training. All implementation details are deferred to Appendix B. By default, we train 2-layer GCNs with hidden dimension of 256 and snapshot gap K 10. We use all nodes for snapshot large-graph, computation on Flickr, PPI, PPl-large datasets. To scale variance reduced algorithm we employ snapshot large-batch approximation 50% of nodes for Reddit and 15% of nodes for Yelp dataset. We update the model with mini-bater batch size of B 512 and Adam optimizer with a learning rate of 1 = 0.01. We conduct training 3 times for 200 epochs and report the average results. We choose the model with the lowest validation error as the convergence point. A summary of experiment configurations and data statistic can be found in Table 2 and Table in Appendix B. Due to the space limit, more experiments can be found in Appendix C.
both SGCN+ and SGCN+ by randomly selecting
Experiment results. In Table 1 and Figure of SGCN, SGCN+, and SGCN++. We remark that multi-class classification tasks prefer a more sta- ble node embedding and gradient than single-class classification tasks. Therefore, even the vanilla Exact, GraphSAGE and VRGCN already outperforms other baseline methods on PPI, PPI-large, addition, we ob-
show the accuracy and convergence comparison
and Yelp. Applying variance reductions can further improve its performance. serve that the effect of variance reduction depends on its base sampling algorithms. Even though the performance base sampling algorithm various significantly, the doubly variance reduction can bring their performance to similar level. Moreover, SGCNs suffers an residual error as discussed in Theorem 1, and the residual error is proportional to node embedding approximation variance (zeroth-order variance), where VRGCN has less vari- ance than GraphSAGE because of its zeroth-order variance reduction, and GraphSAGE has less variance than LADTES because more nodes sampled for node embedding approximation.
observe from the loss curves that
8

---

Under review as a conference paper at ICLR 2021
GraphBAdE
GragABAGE (Daubly)
VBOCH (Zeroch)
VROCK (Doubly)
LADIES
LADIES (Doubly)
FallGEN
PP
PPI-1 ango
Flickr
Reddit
Valklation Loss
116
d:55
0.50
0.33
0.4
0,10
0
500
1000
1500
0
500
1000
1500
0
500
1000
1500
300
600
9000
500
1000
1500
lterations
Iterations
Iserations
Iterations
Iterations
Figure 2: Comparing the validation loss of SGCN and SGCN++ on real world datasets.
7 CONCLUSION
In this work, we develop a theoretical framework for analyzing the convergence of sampling based mini-batch GCNs training. We show that the node embedding approximation variance and layerwise gradient variance are two key factors that slow down the convergence of these methods. Further- more, we propose doubly variance reduction schema and theoretically analyzed its convergence. Experimental results on benchmark datasets demonstrate the effectiveness of proposed schema to significantly reduce the variance of different sampling strategies to achieve better generalization.
REFERENCES
Léon Bottou, Frank E Curtis, and Jorge Nocedal. Optimization methods for large-scale machine learning. Siam Review, 60(2):223-311, 2018.
Tianle Cai, Shengjie Luo, Keyulu Xu, Di He, Tie-yan Liu, and Liwei Wang. Graphnorm: A prin- cipled approach to accelerating graph neural network training. arXiv preprint arXiv:2009.03294, 2020.
Jianfei Chen, Jun Zhu, and Le Song. Stochastic training of graph convolutional networks with variance reduction. arXiv preprint 11X1v:1710.10668, 2017.
Jie Chen and Ronny Luss. Stochastic gradient descent with biased but consistent gradient estimators. arXiv preprint arXiv: 1807. / 1880, 2018
Jie Chen, Tengfei Ma, and Cao Xiao. Fastgen: fast learning with graph convolutional networks via importance sampling. arXiv preprint arXiv: 10247, 2018.
Minshuo Chen, Xingguo Li, and Tuo Zhao. On generalization bounds a family of recurrent neural networks. arXiv preprint arXiv: arXiv:/910. /2947, 2019.
Tianyi Chen, Yuejiao Sun, and Wotao Yin. Solving stochastic compositional optimization is nearly as easy as solving stochastic optimization. arXiv preprint arXiv:2008.108// 2020.
Wei-Lin Chiang, Xuanqing Liu, Si Si, Yang Li, Samy Bengio, and Cho-Jui Hsieh. Cluster-gen: An efficient algorithm for training deep and large graph convolutional networks. arXiv preprint arXiv:1905.07953 2019.
Weilin Cong, Rana Forsati, Mahmut Kandemir, and Mehrdad Mahdavi. Minimal variance sampling with provable guarantees for fast training of graph neural networks. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 1393- 1403, 2020.
Vijay Prakash Dwivedi, Chaitanya K Joshi, Thomas Laurent, Yoshua Bengio, and Xavier Bresson. Benchmarking graph neural networks. arXiv preprint arXiv:2003.00988. 2020.
Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. Spider: Near-optimal non-convex op- timization via stochastic path-integrated differential estimator. In Advances in Neural Information Processing Systems, pp. 689-699, 2018.
Alex Fout, Jonathon Byrd, Basir Shariat, and Asa Ben-Hur, Protein interface prediction using graph convolutional networks. In NIPS, 2017.
Vikas K Garg, Stefanie Jegelka, and Tommi Jaakkola. Generalization and representational limits of graph neural networks. arXiv preprint arXiv:2002.06//7 2020.

---

Under review as a conference paper at ICLR 2021
David Gross. Recovering low-rank matrices from few coefficients in any basis. IEEE Transactions on Information Theory, 57(3):1548 1566, 2011.
Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. In Advances Neural Information Processing Systems, pp. 1024-1034, 2017.
Yifan Hu, Siqi Zhang, Xin Chen, and Niao He. Biased stochastic gradient descent for conditional stochastic optimization. arXiv preprint arXiv:2002.10790 2020.
Wenbing Huang, Tong Zhang, Yu Rong, and Junzhou Huang. Adaptive sampling towards fast graph representation learning. In Advances in Neural Information Processing Systems, pp. 4558-4567, 2018.
Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional net- works. arXiv preprint arXiv:1609.02907 2016.
Sashank Reddi, Ahmed Hefny, Suvrit Sra, Barnabás Póczos, and Alex Smola. Stochastic variance reduction for nonconvex optimization. In International conference on muchine learning, pp. 314- 323, 2016.
Ryoma Sato, Makoto Yamada, and Hisashi Kashima. Constant time graph neural networks, 2020. URL
net nwt//frunmtt id=rkgkW64PPP
Michael Schlichtkrull, Thomas N Kipf, Peter Bloem, Rianne Van Den Berg, Ivan Titov, and Max Welling. Modeling relational data with graph convolutional networks. In European Semantic Web Conference, pp. 593-607. Springer, 2018.
Quoe Tran-Dinh, Deyi Liu, and Lam M Nguyen. Hybrid variance-reduced sgd algorithms for monconvex-concavy minimax problems. arXiv preprint arXiv:2006. 2020.
Quan Wang, Zhendong Mao, Bin Wang, and Guo. Knowledge graph embedding: A survey of approaches and applications. IEEE Transactions on Knowledge and Data Engineering, 29(12): 2724-2743, 2017.
Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful are graph neural networks? arXiv preprint 18X1v:1810.00826. 2018.
Shuoguang Yang, Mengdi Wang, and Ethan X Fang. Multilevel stochastic gradient methods for nested composition optimization. SIAM Journal on Optimization, 29(1)611-6659 2012.
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 Knowledge Discovery & Data Mining, pp. 974- 983. ACM, 2018.
Hanqing Zeng, Hongkuan Zhou, Ajitesh Srivastava, Rajgopal Kannan, and Viktor Prasanna. Graph- saint: Graph sampling based inductive learning method. arXiv preprint arXiv: 907.04931, 2019.
Junyu Zhang and Lin Xiao. A composite randomized incremental gradient method. In International Conference Machine Learning, pp. 7454-7462, 2019a.
Junyu Zhang and Lin Xiao. Multi-level composite stochastic optimization via nested variance duction. arXiv preprint arXiv: / /468, 2019b.
Junyu Zhang and Lin Xiao. A stochastic composite gradient method with incremental variance 2019c.
reduction. In Advances in Neural Information Processing Systems, pp.
Yuyu Zhang, Xinshi Chen, Yuan Yang, Arun Ramamurthy, Bo Li, Yuan Qi, and Le Song. Efficient probabilistic logic reasoning with graph neural networks. arXiv preprint a0Xiv:2001.118500 2020.
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. In Advances in Neural Information Processing Systems, pp. 1 11249- 11259, 2019.
10

---

Under review as a conference paper at ICLR 2021
Appendix
Table of Contents
A Connection to composite optimization
12
B Experiment configurations
13
Additional experiments
15
Different sampling strategies
E Detailed algorithms
21
E.1
Description
21
E,2
SGCN
21
E,3
SGCN+
21
E,4
SGCN++
22
SGCN++ without full-batch
23
F
Notations, important propositions and lemmas
24
F.1
Notations for gradient computation
24
F2
Upper-bounded the node embedding matrices and layerwise gradients
25
F3
Lipsehitz continuity and smoothness property of graph convolution layers
26
F4
Lipschitz continouity of the gradient graph convolutional network.
27
G Proof Theorem
G.1 Supporting lemmas
31
G.2
Remaining steps toward Theorem
34
H Proof of Theorem 2
36
H.1 Supporting lemmas
36
H.2
Remaining steps toward Theorem
40
Proof Theorem
41
I.1
Supporting lemmas
41
I.2
Remaining steps toward Theorem
52
J Reproducing experiment results
53

---

Under review as a conference paper at ICLR 2021
A CONNECTION TO COMPOSITE OPTIMIZATION
In this section, we formally compare the optimization problem in training GCNs to the standard composite optimization and highlight the key differences that necessitates developing a completely different variance reduction schema and convergence analysis compared to the composite optimiza- tion counterparts ge.g., see Fang et al. (2018)).
Different objective function. In composite optimization, the output of the lower-level function is treated as the parameter of the outer-level function. However GCN, the output of the lower-level function is used as the input of the outer-level function, and the parameter of the outer-level function is independent of the output of the inner-layer result.
More specifically, a two-level composite optimization problem can be formulated
N
M
F(O)
=
N
J
M
9 gy(w)) j=1
0 = {w),
(13)
1=1
where f:(.) is the outer-level function computed on the ith data point, gj(-) is the inner-level function computed on the jth data point, and us is the parameter. We denote (.)))) and Vgy(+) as the gradient. Then, the gradient for Eq. 13 is computed as
VF(O)
1 N
1 M
Σ j=1
95 (
(
Σ j=1
(g)(
0 [
(14)
where the dependency between inner- and outer-level sampllin are not considered. One can indepen- dently sample inner layer data estimate g R
Ej=1 Oj (w) and Vğ Dド= Vg (w). sam-
M
M
ple outer layer data to estimate V
V fi()) then estimate V VF by using [V] Vg.
By casting the optimizaion problem GCN as composite optimization problem in Eq. 13, we have
((0)= Σ
Loss(h (L) 9:). 0 (w())
iEVB
H = a XW(()) W(4 o (i [L--2 x x... . [ (I"xwa))... ...).
()),
(15)
W
which is different from the vanilla GCN model. To see this, we note that in vanilla GCNs, since the sampled nodes at the fth layer are dependent from the nodes sampled at the (e + 1)th layer, have []] ( * L. However in Eq. 15, since the sampled nodes have no dependency on the weight matrices or nodes sampled at other layers, we can easily obtain [[ = L. These key differences makes the analysis more involved and are reflected in all three theorems, that give us different results.
Different gradient computation and algorithm. The stochastic gradients update the parameters in Eq. 15 are computed as
aC(0) AC(0)
OW
πI
(16)
Ow(3) OW())
j=(+1
However in GCN, there are two types of gradient at each layer (i.e., ( and G(O) that are fused with each other (i.e., (( is part of (e-1) and ((0) is part of ((e))) but with different functionality. ((0 is passing gradient between different layers, (()) is passing gradient weight matrices.
These two types of gradient and their coupled relation make both algorithm and analysis different from Zhang & Xiao (2019b). For example in Zhang & Xiao (2019b), the zeroth-order variance
reduction is applied to ((()) in Eq. 13 (please refer Algorithm 3 in Zhang & Xiao (2019b)).
12

---

Under review as a conference paper at ICLR 2021
where W W 1-1 is used a control variant reduce the variance of W i.e.
W
+ x - o aL
(17)
However in +- the zeroth-order variance reduction is applied H Because the node sampled at the tth and (t- iteration are unlikely the same, we cannot directly use H to reduce
the variance of H. Instead, the control variant in SGCN++ is computed by applying historical weight w(() t-1 on the historical node embedding from previous layer H (t-1) i...
H = - H (H (!") o (L, H
1-))
(18)
These changes are not simply heuristic modifications, but all reflected in the analysis and the result.
Different theoretical results and intuition. The aforementioned differences further result in novel analysis of Theorem 1, where we show that the vanilla sampling-bass GCNs suffer a residual enror Ab that is not decreasing as the number of iterations T increases, and this residual error is strongly connected to the difference between sampled and full Laplacian matrices. This is one of our novel observations for GCNs, when compared to (1) multi-level composite optimization with layerwise changing learning rate Yang et al. (2019); Chen et al. (2020), (2) variance reduction based methods Zhang& & Xiao (2019b), and (3) the previous analysis et al. (2018); Chen Luss (2018). Our observation can be used as a theoretical motivation on using first-order and doubly variance reduction, and can mathematically explain why VRGCN outperform GraphSAGE, even with fewer nodes during training. Furthermore, as the algorithm and gradient computation are different, the theoretical results in Theorems 2 and are also different.
the convergence of GCNs Chen
B EXPERIMENT CONFIGURATIONS
Hardware specification and environment. We run our experiments on single machine with Intel i5-7500, NVIDIA GTX 1080 GPU (8GB memory) and, 32GB RAM memory. The code is written in Python 3.7 and we use PyTorch 1.4 on CUDA 10.1 train the model on GPU. During each epoch, we randomly construct 10 mini-batches in parallel.
Implementation details. To demonstrate the effectiveness of doubly variance reduction, we mod- ified the PyTorch implementation of GCN (Kipf & Welling, 2016)3 to add LADIES (Zou et al., 2019), Fast (Chen et al., 2018), GraphSAGE (Hamilton et al., 2017), GraphSAINT (Zeng et al., 2019), VRGCN (Chen et al., 2017), and Exact sampling mechanism. Then, we implement SGCN+ and + on the top of each sampling method to illustrate how zeroth-order variance reduction and doubly variance reduction help for GCN training.
By default, we train 2-layer GCNs with hidden state dimension of 256, element-wise ELU as the activation function and symmetric normalized Laplacian matrix = D-1/2 AD-12 We use mean-aggregation for single-class classification task and classification. The default mini-batch batch size and sampled node size are summarized in Table 2. We update the model using Adam optimizer with learning rate of 0.01. For SGCN++, historical node embeddings are first calculated on GPUs and transfer to CPU memory using PyTorch com- mand Tensor. to (device)
for multi-class
Therefore, no extra GPU memory is required when training with SGCN++. To balance the staleness of snapshot model and the computational efficiency, as default choose snapshot gap K 10 and early stop inner-loop if the Euclidean distance between current step gradient to snapshot gradient is larger than 0.002 times the norm of snapshot gradient.
Table 2: Configuration of different sampling algorithms during training
GraphSAGE
VRGCN
Exact
FastGON
LADIES
GraphSAINT
Mini-batch size Sampled neighbors
512
512
512
512
512
2048
All
Samples in layerwise
4096
512
2048
'https://g httt:::///gttbbb..o///t/////////
13

---

Under review as a conference paper at ICLR 2021
PPI: Training loss Epochs
P7: Validation Epechs
0.7
(SG:
2.10
GrapnSage
GrapnSage 0.38
SruphSage
Graph
0.6
0.50
Draph Fl 0.72
0.5
0.55
0.50
0.4
0.40
0.3
0.35
500
1000
1500
2000
500
1010
1500
2000
Spocts
Epochs
PPI: Training loss Epochs
P7:: Validation loss Epochs
0.0
Exact (Deubty) SOD:
0.6
(Diullyy
0.05
0.50
SSD;
Fl 0.76
0.42
0.5
0.55
0.4
0.40
0.40
0.35
0.2
530
1000
1500
2000
2500
3000
502
1000
1500
2000
2500
3000
Spocts
Esochs
Figure 3: Comparison of doubly variance reduction and vanilla sampling- GCN training on PPI dataset with SGD (learning rate 0.1) and Adam optimizer (learning rate 0.01). All other config- urations are as default.
During training, for each epoch we construct 10 wultiprocessing and perform training on the sampled 10 mini-batches. To achieves a fair comparison of different sampling strategies in terms of sampling complexity, we implement all sampling algorithms using numpy. random and
in parallel using Python package
package.
We have to emphasize that, in order to better observe the impact of sampling on convergence, we have not use any augmentation methods (e.g., "layer normalization". "skip-comnection" and "atten- tion""), which have been proven to impact the GCN performance in Cai et al. (2020); Dwivedi et al. (2020). Notice that we are not criticizing the usage of these augmentations. Instead, we use the most primitive network structure better explore the impact of sampling and variance reduction on convergence.
Comparison of SGD and Adam. is worth noting that Adam optimizer is used as the default op- timizer during training. We choose Adam optimizer over SGD optimizer for the following reasons:
(a) Baseline methods training with SGD cannot converge when using a constant learning rate due to the bias and variance in stochastic gradient (Adam has some implicit variance re- duction effect, which can alleviate the issue). The empirical result of SGD trained baseline models has a huge performance gap to the one trained with Adam, which makes the com- parison meaningless. For example in Figure 3, we compare Adam and SGD optimizer on PPI dataset. For Adam optimizer we use PyTorch's 's default learning rate 0.01, and for SGD optimizer we choose learning rate as 0.1, which is selected as the most stable learning rate from range [0.01, 1] for this dataset. Although the SGD is using a learning rate 10 times larger than Adam, it requires 100 times more iterations than Adam to reach the early stop point (valid loss do not decrease for 200 iterations), and suffers giant performance gap when comparing to Adam optimizer.
(b) Most public implementation of GCNs, including all implementations in Py Torch Geometric and DGL packages, use Adam optimizer instead of SGD optimizer.
(c) In this paper, we mainly focus on how to estimate a stabilized stochastic gradient, instead of how to take the existing gradient for weight update. We employ Adam optimizer for all algorithms during experiment, which lead to a fair comparison.
14

---

Under review as a conference paper at ICLR 2021
Dataset statisties. We summarize the dataset statistics in Table 3.
Table 3: Summary of dataset statistics. m stands for multi-class classification, and stands for single-class.
Dataset
Nodes
Edges
Degree
Feature
Classes
Train/Val/Test
PPI
14, 755
225, 270
15
50
121(m)
06%/12%/222%
PPI-Large
56, 944
818, 716
14
50
121(m)
19% 10%
Flickr
89, 250
899, 756
10
500
7 (s)
50%/25% /25%
Reddit
232, 965
106 919
50
602
41(s)
66%/10%/220
Yelp 716, 847
6,977, 410
10
300
100(m)
75%/10%/155
C ADDITIONAL EXPERIMENTS
Effective of variance reduction. In Figure 4, we empirically evaluate the effectiveness of variance reduction by comparing the mean-square error of the stochastic gradient and training loss curve of different sampling strategies on Reddit t dataset.
Because Exact sampling using all neighbors for the node embedding approximation, it only affected by layerwise gradient variance (first-order variance). Therefore, employing first-order vari- ance reduction on Exact sampling can significantly reduce the mean-square error of stochastic gradient to full- full-gradient, and speed up the convergence.
Different from Exact sampling that the exact node embeddings are available during training, layer- wise sampling algorithm LADIES and nodewise sampling GraphSAGE approximating the true node embeddings using a subset of nodes (neighbors). Therefore, these methods both suffer from node embedding approximation variance (zeroth-order variance) and layerwise gradient variance (first-order variance). As a result, applying zeroth-order variance reduction and first-order variance reduction simultaneously necessary to reduce the rquare error of the stochastic gradient and speed up the convergence.
LoR deral 2006
Eigir Capk sisect eren
1.0
4.50
Comparing the mean -quaee error of stochastic gradient to full gradient and training loss of SCCN, SGCN+, +GCN+ in the first 200 iterations of training process on Reddit it dataset.
Figure 4:
GPU memory usage. In Figure 5, We compare the GPU memory usage of SGCN and + We calculate the allocated memory by torch. cuda. www.mmmmr..... aeemoryyalllcatte.
which is the current GPU memory occupied by tensors in bytes for given device. We calculate the maximum allocated memory by torch. cuda max memory all allocated, which is the maximum GPU memory oc- cupied by tensors in bytes for a given device.
15

---

Under review as a conference paper at ICLR 2021
From Figure 5, we observe that neither running full-batch GCN nor saving historical node embed- dings and gradients will significantly increase the computation overhead during training. Besides, since all historical activations stored outside GPU, that +GCN++ only requires several megabytes to transfer data between GPU memory to the host, which can be ignored compared to the memory usage of calculation itself.
monery elecatoe Spoch
2135
2150
2112
UVd
LVD -iead eitro
Libil pred exed
0.
5.2
16.5
15.6
0.0
25
12.3
P11 VMaoryy alocatedd Epech
-aa mamary alocated Epech
LODVE
Libees (Derth yeer.
via 17A
Spack
Figure Comparison of GPU memory usage of SGCN and SGCN++ on Flickr and PPI dataset.
Evaluation of total time. In Table 4 and Table 5, we report the average time of doubly variance duced LADIES+ and vanilla LAD IES. We classify the wall clock time during the training process into five categories:
Snapshot step sampling time: The time used to construct the snapshot full-batch or the snapshot large- batch. practice, we directly use full-batch training for the smaller datasets (e.g., PPI, Par-large, and Flickr) and use sampled snapshot bate-bbtch for large datasets (e.g., Reddit and Yelp). When constructing snapshot large-batcle the Exact sampler has to go through all neighbors of each node using for-loops based the graph structure, such that it time-consuming.
Snapshot step transfer time: The time required to transfer the sampled snapshot batch nodes and Laplacian matrices to the GPUs.
Regular step sampling time: The time used construct the mini-batches using layerwise
LAC IES sampler.
The time required transfer the sampled mini-batch nodes and Laplacian matrices to GPUs, and the time to transfer the historical node embeddings and the stochastic gradient between GPUs and CPUs.
Regular step transfer time:
Computation time: The time used for forward- and backward wwwaard-roppagatin...
Notice that we are reporting the total time per iteration because the vanilla sampling-bor based method cannot reach the same accuracy as the doubly variance reduced algorithm (due to the residual error as shown in Theorem 1).
From Table 4 and Table 5, we can observe that the most time-consuming process in sampling- based GCN training data sampling and data transfer. The extra computation time introduces by employing the snapshot step is negligible when comparing to the mini-batch sampling time during each regular step. Therefore, promising future direction for large-scale graph training developing provable sampling algorithm with low sampling complexity.
Evaluation of snapshot gap for SGCN+ and SGCN+ Doubly variance reduced SGCN++ requires performing full-batch (large-batch) training periodically to calculate the snapshot node embeddings and gradients. A larger snapshot gap K can make training faster, but also might make the snapshot node embeddings and gradients too stale for variance reduction. In this experiment,
evaluate the effect of snapshot grap on training by choosing mini-bater size as B = 512 and change the inner- loop intervals from K mini-batches K 20 mini-batches. In Figure 6 and Figure 7, we show
16

---

Under review as a conference paper at ICLR 2021
Table 4: Comparison of average time (1 snapshot step and 10 regular steps) of doubly variance reduced + with regular step batch size as 512. Full-batch is used for snapshot step on PPI. PPI-Large, and Flickr. 50% training set nodes are sampled for the snapshot step on Reddit, and 15% training set nodes are sampled for the snapshot step on Yelp.
Time (second)
PPI
PPI-Large
Flickr
Reddit
Yelp
Snapshot step sampling
0.182
0.355
0.221
18.446
21.909
Snapshot step transfer
0.035
0.070
0.036
0.427
0.176
Regular step sampling
1.128
1.322
0.899
9.499
9.102
Regular step transfer
0.393
0.459
0.250
0.550
0.372
Computation
0.215
0.196
0.136
0.399
0.139
Total time
1.954
2.377
1.442
29.321
31.697
Table 5: Comparison of average time (10 regular steps) of LADIES with regular step batch size as 512
Time (second)
PPI
PPI-Large
Flickr
Reddit
Yelp
Regular step sampling
1.042
1.077
0.977
9.856
9.155
Regular step transfer
0.036
0.047
0.016
0.496
0.041
Computation
0.077
0.068
0.034
0.029
0.082
Total time
1.156
1.192
1.028
10.381
9.278
the comparison of training loss and validation loss with different number of inner-loop intervals for SGCN++ and SCCN+ on Reddit dataset, respectively. We can observe that the model with a smaller snapshot gap requires less number of iterations reach the same training and validation loss, and gives us better generalization performance (FI (FI-score).
4.0
0.8
LESPLES Ext5T) unend Mi0e
3.14
Cloon Pretiy Plur whi (svereaand
6.5
6.9
475
0.2
1.74
415
4.6
8.14
Peratioss
475
Shenip sar pellaae
2911
0.31
011
0.35
429
4.20
420
0.71
422
010
Soarders
Figure 6: Comparison of training loss, validation loss, and F1-score of SGCN++ with different snapshot gap on Reddit dataset.
Evaluation of large-batch size for SGCN+ and SGCN+ The full-batch gradient calculation at each snapshot step is computationally expensive. Heuristically, we can approximate the full-batch gradient by using the gradient computed on large-batch of nodes. Besides, it is worth noting that large-batch approximation can be also used for the node embedding approximation in zeroth-order variance reduction. In SGCN+, saving the historical node embeddings for all nodes in an extreme large graph can be computationally prohibitive. An alternative strategy is sampling a large-batch during the snapshot step, computing the node embeddings for all nodes in the large-batch, and saving the freshly computed node embeddings sampled from the large-batch during the regular steps. Let denote B' as the snapshot step large- batch size and B denote the regular step mini-batch size. By default, we choose snapshot gap as
the storage. After that, mini-batci nodes
17

---

Under review as a conference paper at ICLR 2021
Fuddi: The ring bett berntione
Feeds: Traiiig Yeer Fiernores
Fedd:: Thainng loen lervitiee
omarllapy ligh danch Roda coeaa e
Fastack daalie 8--22 Lasare
0.3
3.4
4:40
0:35
300
6.5
2003
1270
205
0.00
302
350
50%
220
Vilidation lous Barallor
Sdation Sanu/ Seratlora
Mdation haae/ Beruilora
Cranidact Sarth And5 Forttnn
LAUES SUree
Fartich Quntt
4.43
3.10
0.40
3
0.2-
6.00
2002
12.0
10.0
2207
Figure 7: Comparison of training loss, validation loss, and Fi-score of SCCN+ with different snap-
shot gap on Reddit dataset.
K = 10, fix the regular step batch size as B = 512, and change the snapshot step batch size B' from 20, 000 (20K) 80, 000 (80K). In Figure 8 and Figure 9, we show the comparison of training loss and validation loss with different snapshot step large-batch size B' for SGCN++ and SGCN+, respectively.
545
kodelt Trai Reretiore
2.15
Feeds Thareng Fivertors
Tedil Resstees
Gesetint Avecivs treel EE scan Uos der seroee 20
(Me)) Jou del trodii El soenc peit dor den
(eick Snalle Loarns (All rwn le
3.2
4.30
4.3
3
6:30
4:2
250
250
190
4.24
18575
0:2
Figure 8: Comparison of training loss, validation loss, and F1-score of SGCN++ with different snapshot large-batch size on Reddit dataset.
The effect of mini-batch size. In Figure 10, show the comparsion of training loss and validation loss with different regular step mini-batch size. By default, we choose the snapshot gap as K = 10 fix the snapshot step batch size B' = 80, 000, and change the regular step mini-batch size B from 256 to 2, 048. Besides, we note that subgraph sampling algorithm GraphSAINT requires an extreme large mini-batcl size every iterations. In Figure 11, we explicitly compare the effectiveness of mini-batch size on doubly variance reduced SraphSAINT++ and vanilla GraphSAING and show that a smaller mini-batch is required by GraphSAINT+
Evaluation of increasing snapshot gap. Snapshot gap size K serves as a budge hyper- hyper- parameter that balances between training speed and the quality of variance reduction. During training, as the number of iterations increases, the GCN models convergences to a saddle point. Therefore, it is interesting to explore whether increasing the snapshot gap K during the training process can obtain
18

---

Under review as a conference paper at ICLR 2021
Foski: Thamieg ROG/ Tenal 675
LANEES 10404
01.15
LASET Dhvel Budtict 20413
0.350
0.43
6.39
0.33
0,000
037
0.4
6:13
140
0.00
25
125
375
TacilC: V. W.A.Ann houu Beselinno
Poodd L. VoBiiation Iaso (Taeraaare
Fudddic Via RAAUUr Saual Resollano
Fasttio
0.4
0.44
0.24
180
122
125
Figure 9: Comparison of training loss, validation loss, and Fi-score of SCCN+ with different snap- shot large-batch size on Reddit t dataset.
LOINNIOON CO OA cubinn ders Coundius
EALARI 8-01271 sarad
0.4
0.0
0.5
Perger
Lyte
2420
4.6
Be1
0.4
4.4
474
0.2
174
Figure 10: Comparison of training loss, validation loss, and F1-score of SGCN++ with different mini-batch size on Reddit dataset.
a speed boost. Figure 12, we show the comparison of validation loss of fixed snapshot gap K = 10 and gradually increasing snapshot gap K = 10 + 0.1 x 8,8 1, where is the number of snapshot steps has been computed. Recall that the key bottleneck for SGCN++ is memory budget and sampling complexity, rather than snapshot computing. Dynamically increasing snapshot gap can reduce the number of snapshot steps, but cannot significantly reduce the training time but might lead to performance drop.
D DIFFERENT SAMPLING STRATEGIES
In this section, highlight the difference between node-wise sampling, layer-wise sampling, and subgraph sampling algorithms.
Node-wide sampling. The main idea of node-wise sampling is first sample all the nodes needed for the computation using neighbor sampling (NS), then train the GCN based on the sampled nodes. For each node in the éth GCN layer, NS randomly samples of its neighbors at the (f 1)th GCN layer and formulate [(0) by
[NMI x Li,j.
if E (1) otherwise
(19)
0.
19

---

Under review as a conference paper at ICLR 2021
Craphesim: Grap traphealtte--ee -ioubby with Baich
Grepitaa 522 Teat
Grepntient fool Plrd scr0
Gregptia Ter
110
135
110
Figure 11: Comparing the validation loss and Fl-score of GraphSAINT and GraphSAINT with different mini- batch size on Reddit dataset
Effectiveness of gradually increasing snapshot gap K during training on wallclock time (second) and accuracy on PPI dataset. We choose snapshot gap K = 10 for fxed-K. For increasing-Ae K, we choose snapshot gap K 10 + 0.1 x s, 1,2...., where is the number of snapshot steps.
Figure 12:
where (i) is the full set of ith node neighbor, (i) is the sampled neighbors of node for éth GCN layer. GraphSAGE Hamilton et al. (2017) follows the spirit of node-wise sampling where it performs uniform node sampling the previous layer neighbors for fixed number of nodes to bound the mini-batch computation complexity.
Liyer-wise sampling. To avoid the neighbor explosion issue, layer-wise sampling is introduced to controls the size of sampled neighborhoods in each layer. For the êth GCN layer, layer-wise sampling methods sample set of nodes ((() C V size s under the distribution p to approximate the Laplacian by
I
exps 0,
Li.js
if E B(O)
(20)
otherwise
Existing work FastGCN Chen et al. (2018) and LADIES Zou et al. (2019) follows the spirit of layer-wise sampling. FastGCN performs independently node sampling for each layer and applies important sampling to reduce variance and results in a constant sample size in all layers. However, mini-batches potentially become too sparse to achieve high accuracy. LADIES improves F FastGCN by layer-dependent sampling. Based on the sampled nodes in the upper layer, it selects their neigh- borhood nodes, constructs bipartite subgraph, and computes the importance probability accord- ingly. Then, it samples a fixed number of nodes based on the calculated probability, and recursively conducts such a procedure per layer to construct the whole computation graph.
Subgraph sampling. Subgraph sampling is similar to layer-wise sampling by restricting the sam- pled Laplacian matrices at each layer are identical
Lags if B
ZO = =...= I(()
sxps 0.
(21)
otherwise
For example, GraphSAINT Zeng et al. (2019) can be viewed as special case of layer-wise sam- pling algorithm FastGCN by restricting the nodes sampled the 1-st to (L-1)th layer the same as the nodes sampled at the Lth layer. However, GraphSAINT requires significant large mini-bated size compared to other layer-wise sampling methods. We leave this as potential future direction to explore.
20

---

Under review as a conference paper at ICLR 2021
E DETAILED ALGORITHMS
E.1 DESCRIPTION
In order to help readers better compare the difference of different algorithms, we summarize the vanilla sampling-base basen GCN training algorithm SGCN in Algorithm 3, zeroth-order variance re- duced algorithm SGCN+ in Algorithm 4, and doubly variance reduced algorithm SGCN++ in Al- gorithm 5. In addition, we illustrate Figure 13 the relationship hetween node embedding ap- proximation variance and layerwise gradient variance to the forward- and wookkward-roopggaiiii We remark that zeroth-order variance reduction SGCN+ is only applied during forward-propagation and doubly (zeroth- and first-order) variance reduction SGCN+ applied during forward- and backward- propagation, simultaneously.
Forward- Porwand-propagatiio (Node emhedding approximation variance)
m xw
( a (()
i oi
i1()) [Rg R...
1 layer
nd layer
1 layer
Backward-propagation (Layerwise gradient variance)
107g
Figure 13: Relationship between the two types of variance with the training process, where em- bedding approximation variance (zeroth-order variance) happens during forward- prward-proppgation and layerwise gradient variance (first-order variance) happens during packward-propagaatio.
E.2 SGCN
Algorithm SGCN: Vanilla sampling-bass based GCN training method
1: Input: Learning
2: for t = 1,...,T's do
3: Sample mini-batch V88 V
4: Calculate node embeddings using
H = H (1-1) where H X.
(22)
5: Calculate loss as [(0,) = B Eiev, Loss (L)
:/i)
6: Calculate stochastic gradient VC(0.) G (8) | as
G [
T
D (++1) vo(2.)))
(+1) vo(2,)) w!o) D((+1) Oī(0.)
(23)
5 = [[] (ī(
0H())
7: Update parameters 0.+1 0, -
8: end for
9: Output: Model with parameter 0r+1
E.3 SGCN+
Algorithm S Zeroth-order variance reduction (Detailed version of Algorithm 4)
1: Input: Learning rate > 0, snapshot gap K
2: for t 1,....7 T do
3: if mod K then
4: Snapshot steps
21

---

Under review as a conference paper at ICLR 2021
5: Calculate node embeddings and update historical node embeddings using
z -LH ( H. = )) Z z
(24)
6: Calculate loss as ((0.) NV Ei (hs (L) .36))
7:
Calculate full-batch gradient (C())) {G (G
=t=1
G = LH 1 T D ((+1) (o(z
(25)
D ( : L T (D. ( vo(z, w! ((++1)_ OC(O.)
(H(4)
8:
Update parameters 0,+1 0, nVL(()))
9:
else
Regulaa steps
11:
Sample mini-batch VB C CV Calculate node embeddings using
2,0 = L H ( w!O LH w 1-1 H. =(
(26)
Calculate loss Σ(0.) 1/8 Eieny Loss (hs (L) .//
13:
14:
Calculate the stochastic gradient VL():) = (
[==1
G == [ ((
D ((+1) vo(2.)))
(27)
D ( = I 1 ((t+t)) Vo(Z.") w! D (1+1) OC(0.)
(H())
Update parameters as 0.11 0: (
16: end if
17:
end for
18:
Output: Model with parameter 0r+1
E.4 +
Algorithm SCCN++: Doubly variance reduction (Detailed version of Algorithm 5)
1: Input: Learning rate > 0, snapshot gap K
2: for = 1,... T do
3: if mod K then
4:
Snapshot steps 5: Calculate node embeddings and update historical node embeddings using
z LH! D H! = oz) 2 + -
(28)
6: Calculate loss as C(O.) 1 Ei Loss(h (1) ,16)
7:
Calculate the full-batch gradient VL(O,) = {G 1
18=1 as
G = LH (1 I/T ( vo(z(())
(29)
D (0) L (D. (+11 vo(z(()) (z) w! D.d (L+1) OC(0.)
0H
Save the per layerwise gradient 7 + (( D D for all & [4] Update parameters as 6.41 = 0, 1) v((,)
8:
9:
10:
else
11:
& Reqular steps
12:
Sample mini-batch V/8
22

---

Under review as a conference paper at ICLR 2021
13:
Calculate node embeddings using
7 = 2 + + (2 w! -L H. =0(
(30)
Calculate loss as [(0.) B Eiery Loss (L) :1/i)
14:
15:
Calculate the stochastic gradient V((0,) = (G 1
1==1 as
( G 1-1 G
D (++1) (o
6-
D (++1)
Vo Vo
D D (0
(0)
[
D (e+1)
va(Z.)) (w
(31)
Vo 1) |w D (1++1) Oā(0.)
1 (Di
OH!
16: Update parameters 0,+1 0: (
17:
end if
18:
end for
19:
Output: Model with parameter 0r+1
E.5 S WITHOUT FULL-BATCHI
Furthermore, in Algorithm 6, provide alternative version of SGCN++ that does not require full- batch forward- and backward- backwarr--ropaagttnn propagation at the snapshot step. The basic idea to approximate the full-batch gradient by sampling large mini-batch V/s of size B' V using Exact sampling. then compute the node embedding matrices and stochastic gradients on the sampled large-batch Va-
Algorithm 6 SGCN++ (without full-batchi Doubly variance reduction
1: Input: Learning rate 0, snapshot gap K >>
2: for = ,,...,,77 do
3: if' mod K = then
4: Snapshot steps
5: Sample large-batch V/8 of size B' and construct the Laplacian matrices (() for each layer
using all neighbors, i.e.,
(
[4.5
i E
(32)
1. lo
otherwise
6: Calculate node embeddings and update historical node embeddings using
z i = H w!( H!!) = =z( 2 < (4)
(33)
7: Calculate loss as L(,.) 1.- Eveyy (hss(h) .32))
8:
Calculate the approximated snapshot gradient VL(0,) (G(1) 34m
G ( = L HH TIT (D!tt)oVOoZ vo(Z
(34)
((+1) vo(z!() w D. (L+1) aC(O.)
D ( = L [ (D(ttt)
2H!L)
Save the per layerwise gradient ( + G D D vl e (L) Update parameters 0 = 0, V((,))
9:
10:
11:
else
12:
Regular steps Sample mini-batch V88 C V/8 Caleulate node embeddings using
13:
14:
Z = I
H. (Z
(35)
23

---

Under review as a conference paper at ICLR 2021
15:
Calculate loss as Ž(0,) B Eiey, Loss (h [hs
yi)
Calculate the stochastic gradient VL( ( ] as
12=1
(0) G 1-1 +
H T [ ((+1)
vo(2.)
vo(Z.-1)
D 0 D + 1 [ (++)) vo(Z.) !
(36)
D D 111
vo(Z)-1)) w
D (2+1) Oā(0.)
aH!")) (T)
17:
Update parameters 0,+1 0; v(
18:
end if
19:
end for
20:
Output: Model with parameter 0r+1
The intuition of snapshot step large-batch approximation stems from matrix Bernstein inequal- ity Gross (2011). More specifically, suppose given G, Rdxd be the stochastic gradient computed by using the ith node with Exact sampling (all neighbors are used to calculate the exact node em- beddings). Suppose the different between G; and full-gradient [[G,)] is uniformly bounded and the variance is bounded:
|| - E s u, EI| B 31 s²
(37)
Let G' as the snapshot step gradient computed on the sampled large batch
G'= 1
Σ G
(38)
B
SEV3
By matrix Bernstein inequality, we know the probability of l larger than some constant decreases exponentially the size of the sampled large-batch size 3' increase, i.c..
Pr[ - E(G, 2 c) S 2dexp (-n min
2 E
(39)
4 2μ
Therefore, by choosing a large enough snapshot step batch size B', we can obtain good approxi- mation of full- gradient.
F NOTATIONS, IMPORTANT PROPOSITIONS AND LEMMAS
F.1 NOTATIONS FOR GRADIENT COMPUTATION
We introduce the following notations to simplify the representation and make it easier for readers to understand. Let formulate each GCN layer in FullGCN as a function
H (0 (H(--), w(0) = ((-I E Nxde
(40)
2
and its gradient w.r.t. the input node embedding matrix D E Nx is computed
O Vn n (0) (++))) ((--)), w()) = [L] (D(++)) od d' (LH (LH [w())]]
(41)
and its gradient w.r.t. the weight matrix ( E RR 1 xde is computed as
G = (0 ((111) w(0) = (H D ((e+t) oo (Lo (L (0 (l- [((o))
(42)
24

---

Under review as a conference paper at ICLR 2021
Similarly, formulate the calculation of node embedding matrix H layer in SGCN as
E E RNxx at each GCN
H = (F0 w o HL
(43)
Z
and its gradient w.r. the input node embedding matrix ( E R Nxdx is computed as ( V 7 (D(et)) (e)-)), w(f) == 1 ((++1) I w(()) IWI)TT
(44)
and its gradient w.r.t. the weight matrix G E [R xd is computed
G ( = V J D (++1) H w()) = I ( 1H IT (ittt) oo
(45)
Let us denote the gradient of loss w.r.t. the final node embedding matrix as
D (1++1) Loss (L),
(H()) E Nxdz, D (
1
Losssh (L)
Ui)
ERd2
N
oh (L)
(46)
((++1) Öoss (1) v)
E Nxd2
D
1 (ievo)
(h (L)
y/s
E Rd
UH
B
oh! (L)
Notice that B(L+1) is a N x dr matrix with the row number i VB are non-zero vectors. Then we write the gradient for the 9th weight matrix in FullGCN and SGCN as
G = (V " (D(+)) H (1!))) w( (.) w H ((-1), w())
H O w (e-1), w())
G Vw Fl0 (V ((tt)) (H..VH) D (L+1) H w("))... ..., ...
(47)
F.2 UPPER-BOUEE BOUNDED ON THE NODE EMBEDDING MATRICES AND LAYERWISE GRADIENTS
Based on the Assumption 3, we first derive the upper-bound on the node embedding matrices and the gradient passing from /th layer node embedding matrix the (E 1)th layer node embedding matrix.
Proposition 2 (Detailed version of Proposition 1). For any [T]) the Frobenius norm of node embedding matrices, gradient passing from the Êth layer node embeddings to the (E 1)th are bounded
|H S BH, do (o(L || BD, ||
do (o ( ||F S Bn
1 S Bu,
@H(-1)
0H(E-1)
(48)
where
T--t Bи (CoBLABW) Bh, Bp (BLACoBW)
Class
(49)
Proof.
H (0 = l (-
CCABLA (H
(50)
(C,,Bra B 1IXI1r (C.Bra Bur
H H ((
s C.BBLA (
(51)
<(CoBIABW) ) "1XIIr s (CoBraBur) Bu
25

---

Under review as a conference paper at ICLR 2021
DD | D (D od
( (
S BLACO Bw [D(tt))
(52)
(BLACoBW Bw
(D S
(BEAC Bw )
Closs
[ || = || (D (Dittt) o0 "O (
BLACo BLACCBW B
(53)
D ( ( S (BLACoBW)
Live Closs
F.3 LIPSCHITZ CONTINUITY AND SMOOTHNESS PROPERTY OF GRAPH CONVOLUTION LAYERS
Then, we derive the Lipschitz continuity of (..)) and its gradient. Notice that the following result hold for deterministic function (..)) as well.
Proposition 3. (.,-) is CH-Lipschitz continuous w.r.t the input node embedding matrix where Cи = BLABA
Proof.
1 w(0) 7 (H) w
[o[ w(() o(L H (())
(54)
C. | (1-1) - H |r||
S BLABIB H (t-1) H F
Proposition 4. ((.)) is Ow-Lipschitz continuous w.t. the weight matrix where Cw C₂BLABr BH
Proof.
(H(C)) - 7 (H w || -I( - (( w
(55)
BL.BLABH ( w/
Proposition 5. Vn (..,,. LH -ipschitzz continuous where
Lи = max max (BLACoBw Bw Bi, BrA + BL LoBB
(56)
Proof.
7o (++1) ((-1) - - (D. (e+1) 1 w(")
[[[ D ((-))))))))))) WOOJT
(57)
-[D ( (1-1)
(1)) [w(0)
< BrACoBW | ((+1) D ((+1)
||
26

---

Under review as a conference paper at ICLR 2021
H -
((tt)) H w
s <I[[<0]
TO
(58)
11 D
H W()) |
B2, BD By Lo
1) H 1
Vu
w
Vn 100) ((++1) H w
[[[
(11))) w 1
H
(59)
- (
(BLABDC + B BD Bw L.g Bh BH) W -
L
Proposition 6. Vw (io is Lw-Lipsccitz Tnw-Lipxchiit contimuuus where
Lw = max BLA + Bi, BoLo
(60)
Proof.
Vw (D ((11),
w Vw (D( ((-1)), w(0)|
s ((1))
D (e+1)
(61)
-
D
(((tt))
H1 ((O))) |F
< BLABH Co | ((+1) D Di ((+1)
Nw (1+1) H (((-)) w()) Vw Fi (D(+))) - w())lF
< II[ 1 D (e+1)) o (["
(62)
-H (
w())
S (B) BBC. BI BB
H
H. - ||
Nw (1+1) H( w
V
DD (e+1) (e w )r
((++)) 11 D
(63)
- D ((+)))
((" |F
< Bi ( w |F
F.4 LIPSCHITZ CONTINOUITY OF THE GRADIENT OF GRAPH CONVOLUTIONAL NETWORK
Let first recall the parameters and gradients of a L-layer GCN is defined
0 = (W()) W('), VC(0) = ( G(0)
(64)
where ( is defined as the gradient w.r.t. the 9th layer weight matrix. Let us slight abuse of notation and define the distance between two set of parameters 0,, O2 and its gradient as
10: 8211F
|w(
(
|p. [V((()) VL
Σ
I G |F
(65)
Then, we derive the Lipschitz continuous constant of the gradient of a L-layer graph convolutonal network. Notice that the above result also hold for sampling-based GCN training.
27

---

Under review as a conference paper at ICLR 2021
Lemma 1. The gradient of an L-layer GCN is Le-Lipschitz continuous with Lp L(LU max m c + U max L) i.e.
[VL(())) V & Lr|101 0211
(66)
where
Unax = max Cin,w,
je je(0,..,b-1)
(67)
Tw Tw Tit x max {ı, Tuaaa
U U =
min LjE(0...b---11)
Proof. We first consider the gradient w.r.t. the 9th graph convolutional layer weight matrix
|G] _ MF
|V (VH
(+1)
H 11-11) w(
H
w
H 1 w(
-w
(1-1) w (0)
w H (-11
(++1)
w! w(()
11 1
(
H (((+1)
Vw Vw
H
)[F
+IV 10) A
Vw
(1-1) w
H W
H
w |||
+ VI f10 (V
H
(1-1) (())
w
Vw
VHf
H
((1))
H
w
}||
+ +V (0)
(++))
H (t-1) w(()
Vw
H (-") w(
VV ( (D (++1)
H ( w Vw (+1) H (-11) w
(68)
By the Lipschitz continuity of Vw ((0 (-) and V 4 (.)
| 1 <
L I--
Loss H (1) H (1) |F
+ Lw 4 8- (--1) - H 1--1)
1 || + | (4) (4) (1)) +...
(69)
+ Lw LH (H
H
W
(++1) w! (1e
+ Lw L ((-1) H - w
w! lp))
Let define Unax as then we can rewrite the above equation as
U Umax
min LE
Tw T2g1 max 11oss)
(70)
| G/0 ||F I Umax
Σ H H11111)) Umax &
|| -w(')r
(71)
Then, let consider the upper bound of [H((0 -
H H ||r (1 (p
(1) (x, w(").
w(--1), w(T)
-110)
(x. w
w((-")) !
$ [1(00))((-)) (X
(x. w(").
(( "y, ((9)
(72)
- (x, w!l))
w((-)) w
+1 (x, w
w((-1), w(9)
- (x, w(l)) w((-1)) w/
28

---

Under review as a conference paper at ICLR 2021
By the Lipschitz continuity of (.,)) have H
H
1| CH Cw |w -W2" || + + Cw /w/ W
C/,Ow] [ (E |w()) w l1r)
(73)
max
Let define Unax is defined as
U Umax
max jE (0....I.--1)
CinCW
(74)
Plugging it back we have
1
| G || S (LUmax mmmxx Umax L.)
|w/)
(75)
Summing both size from == to = L we have
[VC(())) VC((2)
||
Σ t=1
| - 1l
(76)
< LL max L LUmax C + Umax 1)
Σ (w()) wl('llF)
1=1
< LLLUmaxx I L.max + Umax L) |0 - 02 1F
29

---

Under review as a conference paper at ICLR 2021
G PROOF OF THEOREM
By bias-variance decomposition, decompose the mean-square error of stochastic gradient as
L Σ BI|@ - G || E 10 G |1 1=1
- E|
(77)
E[|
Therefore, have to explicitly define the computation of L which requires computing D( = E[D ( [J and G E|
Let defined general form of the sampled Laplacian matrix L E NRN as
if E (( and j E B(t-1)
(78)
otherwise
where oij is the weighted constant depends the sampling algorithms.
The expectation of (.(() iij is computed as
Eg EjeBre-D
[24] B(0)]
(79)
In order to compute the expectation of SGCN's node embedding matrices, let define the propagation matrix p(t) E x as
p! Bie
[ l B(o]
(80)
where the expectation is taken over row indices i. The above equation implies that under the condi- tion that knowing the ith node B((), we have p!d Ls.s. Vj = {l...., N.. Let consider the mean-aggregation for the ith node
N
Σ
(1-1))
(81)
j=1
Then, under the condition ith node is in B(t), we can replace 7,(!) by P!D which gives us
N
p
(1-1)
(82)
3
As result, we can write the expectation of æ (e)
with respect to the indices i as
BiEBtolzz (8)
l BL0] Eiesto
7(0 1
x( 1 BB())
p(t) 20 (1-1) E B(0]
Bieson
(83)
N
p(t) x (e-1))
jml
Then define ( E Nxde as the node embedding of using full-batch but subset of neighbors for neighbor aggregation, i.e.,
O = -
(84)
where all rows in O non-zero vectors.
30

---

Under review as a conference paper at ICLR 2021
Using the notations defined above, can compute D(1+1) D E N as
R Nxdd, G Rdl-ixd, and
Losss(H
1 Losss(h (1 (7,))) 31)
D(L+)) =E
R Nxde d.
Rd
(85)
(H())
N oh!
and
= D R0 (D (Et))) ( [ ( (D(4+)) oo (P() (
(86)
and
G Vw D 1++1) H (1-1) w()) LH ()
( (87)
As result, we represent G E | as
GO = Vw (n V (++1) H(--1,, ) H w(+h)), (e-1) w()))
(88)
G.1 SUPPORTING LEMMAS
We derive the upper-bound of the bias and variance of the stochastic gradient in the following lem- mas.
Lemma (Upper-bound on variance). We can upper-bound the variance of stochastic gradient in SGCN as
EI EIGH S O p(0 + O([P((( - -(
(89)
Proof. By definition, write the variance in SGCN
B - ||
= f( (V 71641) Vn H(-D)()) H w(t+ H (olm (o
(D (L-1)) w(4) H w(tt)), ((1))) fll
- -ww (V ((tt)) VH
+1+ + V 710
O w(++))) 1) (6-1
(())
w(+))) H  w
-WwF (. HA
(D
((-1)
0 w((+1)) ( w())
- FO ( HA
HO w(+))) (e-1) w
+1+ C+E|)E|IV 70 ((041) DD ((22) H ((+)))
((-1) (0l)
- -V (V
( ((++))), H (1-1)
+ 1[.+ 70)
(11) w Vw f (++1) 1) w(()
< 11
( ||
+ (L+ + 1) 2(
D (++1) (((-1) - v f4 (D(+1)) H -(-1)
1
+ (L+ + 1) B11l V 70
H w(++))
v / F (++1) D (2+)) H() w(+
+ (L+ + 1) V 710 ((+1) (1-1,, w()) F D ((+1) (e-1) w +
(90)
From the previous equation, we know that there are three key factors that will affect the variance:
The difference of gradient with respect to the last layer node representations
([I[D((+1) (4+1)
(91)
31

---

Under review as a conference paper at ICLR 2021
The difference of gradient with respect to the input node embedding matrix at each graph convolutional layer
E f(t1)) (D 2 H
w
Vn (D ((+2)) H w ||]]
(92)
The difference gradient with respect to the weight matrix at each graph convolutional layer
V fle) (D(C+)) ((-1),, w()) V f(0 (D ((-11), w
(93)
First, Let consider the upper-bound of Eq. 91.
[ E[| Loss H (1) y) Losss (H
E||
(
[
(H))
(H
Liose
(L)
(94)
(()) 1
< Lioee Co B2 1 H (1-1) p -
S Lioss C B2By B B|| - |(
Then, let consider the upper-bound of Eq. 92.
w H
Vn (0) ((+1)) H (
D
W (-())
H w(())
[L]
D ((tt))
( [ [
2E|||
D ((+1)
/ ( [H
+2 DD
L0 [L]
) (
< $2B4 BI, By
2 -
< B3
2 23 - p())
-P()) [
1 LI
< OBBII = pd (2) 1 O - L1
(95)
Finally, let consider the upper-bound of Eq. 93.
V w 7()) (D (4+1) H w())
(++1) H)- w(0))
- ( ((Piir((-))) [
$
2E
+ +2
)() LL
+
+2
(1-)))
C L
< 2 B3 D -P(0
+ +2B -
< _ - LI
-
(96)
32

---

Under review as a conference paper at ICLR 2021
Combining the result from Eq. 91, 92, 93 we have |I|
> O ]|
[(]) + +... + - 1())
(97)
O -
LI(i) + O(B[[I(L)) L L
Lemma (Upper-bound on bias). We can upper-bound the bias of stochastic gradient in SGCN as
B| G
O(I - -
(98)
Proof. By definition, we can write the bias of stochastic gradient in SGCN as
E [[| 1
GO
( H(--1), w()) H w(++))) (e!),, ((0)
V Fro (V
A.W. -Vw
(... (.
((-1) w( H w(tht)) ( ( l(l))) |
AAIIBBCT -
Vn 1il
w(1) HO) W ( (( w(o)
((-11 w H w((+))), H ( M (
V 10 MIBAII 1711
, HO, w )() ( w())
- 44 (. SAweAAA
1 w(o)
Ho wichty H ( ((0)) +
4111 Millll V
- Sew/or
0 W ((tt  -.
w())
(--1 | w())
+ 11 1) V
1
V V
(D (++1) H(( w()))
D f
< (L+ (1L+1)4 11 2
+ +I.+ 1) 1.211-
(1+1) H (111)) w(l) Vn f ( (D(+)))
H(L-1), w
+ +L+ E/I ((+))) (D ((+2) H w
VH ( VHz
(e+2) H w(++))
(L + V
V F (E+1) H ((-1), w()) Vw D ( (e+1) H!- w(
(99)
From the previous equation, know that there are three key factors that will affect the bias:
The difference of gradient with respect to the last layer node representations
([I[D(+++) (
(100)
The difference of gradient with respect to the input node embedding matrix at each graph convolutional layer
f([V) HOO w((+1))- - ( (C+22) H()) w(++1) (101)
The difference of gradient with respect to the weight matrix at each graph convolutional layer
f (D (C++1) ((-)), w(1)
Vw (
(1) (((+1)) ((-)), w(
(102)
Firstly, let consider the upper-bound of Eq. 100.
E D D ((++))2) = Ell (L(s) .y) @Loss(H())..) .v)
OH()) [|2]
(L)))
(103)
< L3 B3
(
H (1)
33

---

Under review as a conference paper at ICLR 2021
The upper-bound for |
H 1
[
H (0
- (LH -
= - w(0)
s C. H((-)
LH (6-1) LH
LH 1
(104)
- 200 B p(0) LI 22 Bi,, H -
- 1
- ) + + - -
Therefore, have
||]1 $ -
((+1))
- + -((
(105)
Then, let consider the upper- bound of Eq. 101.
B f0 (D (e+1),
H(-1),
w Vz (D (e+1) H (-1)) w
[[[[[]] D ((++1)
(P ( (()) 1
D (c+1) o o (LH (LH(( ) [W(O]]]1
$ BLA Bi L p (1-1) LH LH LH
LH ( |
2B B3 B 12 B p(0 - [ + BB B 1
H (t-1) H (
=(IIp() - + -
(106)
Finally, let consider the upper-bound of Eq. 102.
B w (()) (D (e+1) ((-1)) w(f) - Vw ( ((-1) H w())
=[[] T D (et))) (P (P())()-)(()))))) -
" ( (( [
H. 1 T (D(+1) oo
2E|| (11
D ((+1)
HT Llie ((mppp---- |
+ +
LH (111)) D
LH Lli-sH
BB
+2B B
LH (H
LH (1
2 B B3 B3 R H H (6
(()) < O( - L)) + +..+++++++ O[ - L
+ BBB p()) L
(107)
Combining the result from Eq. 100, 101, 102 we have
B| 1 > OEE| L +... + O[EIII(( - -L
(108)
G.2 REMAINING STEPS TOWARD THEOREM
By the smoothness of C(O)), we have
T. |1 0,11 2
C(O,+1) < ((0,) + (VL(0.) 0:11 -,))
(109)
-((,) 7u V((1) + n 2
|V
Let = V (1 Note that the weight parameters 0, is a function of history of the generated random process and hence is random. Taking expectation both sides condition
34

---

Under review as a conference paper at ICLR 2021
on Fi and using 7 < 1/Ly have
B[C(0=+)]FF.)
$ (())) XIV + 7LL1
- (UITI-cc0---70000000011111111110
2
C(O)) (ZCC Fil) +
2
C(())) 7 2.V V. + + ]L Ebbb + 2
بااي B |7
" -(()) + 2/4 + B|
Lا 2
(110)
Denote Ab as the upper bound bias of stochasite gradient as shown in Lemma and Denote An the upper bound of bias of stochastic gradient shown in Lemma 2. Plugging in the upper bound of bias and variance, taking expectation over Ft, and rearranging the term we have
2 -E[CCO)) + nLlAn + Ab
V -
(1II)
Summing up from 1 to rearranging we have
2
Σ
[E
(IC(O++)))) + + Ab
7 1-1
(112)
2
(C(0)) - (()") + nL, +Ab
n/
where the inequality (a) due to L()')= S E(((OT+1)l-
By selecting learning rate as = 1/VT. we have
E [ 2(((0,) L()'))
L An Ab
(113)
VT
VT
35

---

Under review as a conference paper at ICLR 2021
H PROOF OF THEOREM 2
H.1 SUPPORTING LEMMAS
In the following lemma, we derive the upper-bound on the node embedding approximation error of each GCN layer in SGCN+ This upper-bound plays important role in the analysis of the upper- bound of the bias term for the stochastic gradient. Suppose the imput node embedding matrix for the é GCN layer as H. (t-1) the forward propagation for the 7th layer in SGCN+ is defined as
(H w()) [ p( H ( w(())
(114)
and the forward propagation for the 4th layer in FullGcN is defined as
f (H,
w(f) L (
(115)
In the following, we derive the upper-bound of (H w(O) - (H, (e)
II EIIH
(116)
Lemma 4. Let denote E E 10 - 1] as the current epoch, let t be the current step. Therefore, for any i {EK + I.... EK + KJ. we have
E LH
w!! 2z 11:1
(117)
$ n x 0 [ |||
1|1
ILII
E E
I
Proof.
LLH [--1) w! (1) (000 1b
=I[LH, (1) LA 1-1
LH 6- LH
-ILH,
LH
|H - (0) 1 + 1| 2,2,1
(118)
+ 2LH w LH (H(-!)) (
-2(LH,
0 LH w/ 1-11 z 2
-2(LH ((-1) W (0) 2 7 Z (
Recall that by the update rule, we have
z Z (6)
H f- (w(( ((O
(119)
and
E (0) 7000 [F = LH.
w!O
(0)
(120)
Taking expectation on both side condition on Fr, we have
E LA
w(() Z (
LH w/ 2 + B11Z 2.2
(121)
LA w LA (--) 2
Then take the expectation over Fr, we have
LH ( 1) w( (
|| < LA
2,2
E1
2
1
(122)
LH.
- LH
36

---

Under review as a conference paper at ICLR 2021
Since know E (EK + 1, ................................ EK + K), denote t = EK + k. & I K such that
B|]LHH - - B|||LH
(1-1)
(1)
EK+ ((6) +
EK
LH
1)
W
EK EK 1F
(0)
(2)
EK
(A)
EK+K Σ
+
BIIZZ
(1) - 2,2 - LIH
-L w 1 11
t= EK+1
(123)
Knowing that we are using all neighbors at the snapshot step (t mod K) = 0, we have (A) 0. As result, we have
E BI
(2-1)
w
EK+ 2 -
EK-A
EK+K
E E
(0) - 2,00 1
ILA w(O
LH
(124)
(=EK+1
(B)
Let take a closer look at term (B).
E -
())
7
LH,
LH.
E[]
LH
w W-11 1F
H
((( (1)
(125)
s
[E
[+]
|
(
H
(C)
Let take closer look at term (C).
H
H - 1)
w
2
E| H
H
w
(
H
(126)
< E
!
(1) --1
|
+
2E
-H
By induction, we have
H
10)
!
2B E W 222,11 - 2
(1-) - W +...
(127)
2 B
(
w
By the update rule of weight matrices, we know
[[] w(O [121
w = B
(128)
Therefore, we have
LH (f-1) w
(1) xx=++
EK+ EK-A
EK+K
Σ
2O
())())(())
(0 [2:]
LI G
1 E - || x E
11-1
t=
(129)
By the definition of ( we have that
B
(130)
37

---

Under review as a conference paper at ICLR 2021
Plugging it back, we have
LH! w! (2) 2, 11:]
E
= El
(6-1)
(1-1)
BK+ EK+*
7
(131)
K x
B P ||
|
E|
Based on the upper-bound of node embedding approximation error of each graph convolutional layer, we derived the upper-bound on the bias of stochastic gradient SGCN. Lemma (Upper-bound on bias). We can upper-bound the bias of stochastic gradient in SGCN+ as
B|| 1 S T
K
(
B
- IILIP)
(132)
t=1
Proof. From the decomposition of bias as shown in previously in Eq. 108, we have
| G
-
(L++)) 2 ((()))) - D(L+1))1 [
+ +L++ 1 1 2 2 B 1 (D (L+1)) 1+1) H (L-11) w()) - ( H(l-1), w(t)
+
+ (L+ 1 E|Il (e+1)) (2+2) H w(+t))) VH ((tt))) (e+2) HO, w((+))
+ (L + V F (D H(1-1), w(0) Vw ( ((+t)) H w
V
(133)
From the previous equation, we know that there are three key factors that will affect the bias:
The difference of gradient with respect to the last layer node representations
B||[((++))
(1
(134)
The difference of gradient with respect to the input node embedding matrix at each graph convolutional layer
E (D (6+2) H w((+)))
Vn (D (6+2) H
w [|))
(135)
The difference of gradient with respect the weight matrix at each graph convolutional layer
W ( (D (1+1) H)- w()) Vw (9) (0+)) H (
(136)
Firstly, let consider the upper-bound of Eq. 134.
(IID((+11) -D(++)112 = E[l ÖLoss(His (L)) y)
(Loss(H(1) y)
(H)))
0H(L)
L2 Lise E. (l] (4) [2]
(137)
< L2 (E ((11
We can decompose E [|[ [()) (Z()) [2] as
E (1) Z (z)))) (1))
= L] LH (1((-1))(())1
< 2|[] ((4) ((4)
LH w (L)
+ 2E 2E LH -
((
(138)
< 2E[ (L)) LH w
|
+ +B2 AB
B - (1-
Bi
< Σ OEI|2 0 LH 4- w(()))
38

---

Under review as a conference paper at ICLR 2021
Using result from Lemma 4, we have
B|ID ((++
p S
Σ A=1
K x O
[B|]]]]]]] 1 - -I
(139)
Then, let consider the upper-bound of Eq. 135.
f H w N (()) (D (e+1) ()) w(
= =B[I[[]T (D (D(++)) o o' ((
((
W - [L] D (f+1) oo (LH - ( (I
BBL By 12
LH w()) [
B3
LH
+ 2B B 13 El|
LH w
(A)
(o41)
where - 2 H ( w! P( w(
1-1-
Let take closer look at term (A), we have
E ((() (H( | & BI (c-)) LH w((-1)
I B C
LH w
+BB By C.B
LH
((()1)
(141) [
Therefore, by induction we have
E|| h (0)) (D (e+1) (e-1) (0) VH
(++1) (++1)
H (t-1) w )
S ( 7 8 LH
[ + OE] LH W
(t-1) (E)
(142)
+ +[E ( - LH (2) ( + ((II LX
Using result from Lemma 4, we have
E||l E H. (00 (D (++1) H ((-1), w()) - H (0) (0 (e+1) H (--1) w
S < x 0 [ f
(143)
[i - || +... + K x O | | || - (
Finally, let consider the upper-bound of Eq. 136.
E V w ( (it91
H ((-11), wr(0) V ( (D (++1) ((111) w())
= 1 HH Llit (( o (Z (2110) ( L D (((+1)
od
[i]
llae <
(E+ o (2(3)
[L 1 D ((+1) oZ [|
23
D (e+1))
o o(
[LH
o o' (L ( H (1
w())) -
IP-1
< 2B3 B
E
H
|
+4 B3 B3
LH - |
+ +AB B3 L ( (
w1 (0)
(144)
By definition, we can write the term (B) as
E|| HH H (()) < C (0) LH 11
< 20 (0) - LH 2) (2>))
(145) (
+ 2C LH w
LH
22)
39

---

Under review as a conference paper at ICLR 2021
Plugging term (B) back and using Eq. 141 and Lemma we have
||| W f (D (
(t-1),, w()) Vw f(() (D(+t)) ((1-1), |
(146)
S s K x O
P - ||با + + x || -
Combining the result from Eq. 134, 135, 136 we have
B||E[G - G (8)
|| S 1 K
(BIP
(147)
H.2 REMAINING STEPS TOWARD THEOREM 2
Now we are ready to prove Theorem By the smoothness [(O)), we have
LF
s + (VL(0))) 0 - -0,) + 10 - 0ell
(148)
= ((0,) - 7 (VC())) v20)) + 17 LF |v
2
1........ (B 1 J. Note that the weight parameters O. is function of history of the generated random process and hence is random. Taking expectation on both sides condition on F, and using 7) < L/Lp we have
Let F; (1B(0))
LF 2
$ ((0)) (VV +
-
-
+ LE
C(.)) (L.
[=]
FF]11
2
LF E n 1| Fd
C(0.) 7( - 2(VC VC bb, (F))
ア| LF E
& ((0.) + 3/2 - || E + +
n
3|
(149)
Plugging in the upper bound of bias and variance, taking expectation over F, and rearranging the we have
E [ < 2 -BLCC)))-EV(((0H)) + MLFAn + 7
(150)
Summing up from = 1 to T, rearranging we have
Σ [E(C(())] - B[C(O++1))) ALFA, + + T/ t=1
(151)
(C())) - C((")) + + +
(0)
np
where the inequality (a) due to L()) E(
By selecting learning rate as = 1/VT. we have
E [ 2(((0)) -(0)))
Lp n
A#
(152)
TT
T
T
40

---

Under review as a conference paper at ICLR 2021
1 PROOF OF THEOREM 3
1.1
SUPPORTING LEMMAS
In the following lemma, we decompose the mean-square error of stochastic gradient the 4th layer E||[@0
GO [2]
as the summation of
The difference between the gradient with respect to the last layer node embedding matrix
BIID((++) 11
(153)
The difference of gradient passing from the (e+ 1) layer node embedding to the 4th layer node embedding
E| V (c+2) ( w V f f((+))
(D((22 (D+
H w (154)
The difference of gradient passing from the êth layer node embedding to the éth layer weight matrix
B w (D (++1) H w - ( (D(+))) (
(155)
Lemma 6. The mean-squar error of stochastic gradient at the tth layer can be decomposed as
BlIG@ -
O(E [(E1][((4))) - (4+1) [+])
+ O( (1)) D (1+1) H ((()) / (D w(2)
O(B -( (D (e+2) H w((+))) (n H w(++
+ / ((-1), w(9) F (D4))) ((1),, w()))
(156)
Proof. By definition, can write down the mean-square error of stochastic gradient
El| - 1
w (V (.. اu ( ((1-)) HO, w((+)) H 1, w
- P (((t))) /HA V (D (L-1)) w(4)) H w(th1), H((-))) ; " 160 [
() + (1+11BB1T IJEIIIV 1 700
( w(++))) ( [- M (0o HO, w H((-)) !
1111 MAMIIIaA 1... - Pw////
 w w -
IA
(1-1))())))
o " ((t ا ا w
- Pw//y
Ho w(+))), )) do (
(D (c+ (e
Ho w (-)) (
(L Ilaat V
M V
w((+)), H (-1)
(())
+ (L 11E V (++) H (-1), w(0) Vw ( (( ()) w
(lllaa1 - (1+1) [1])
+ + 7/ ( ((()) ( H(L-) (( +...
O V
0 VV ((t H(0) w((+))
+ O (++1) 1)1), w(i) F(0 ( H( w())
(157)
Recall the definition stochastic gradient for all model parameters Vū(0.) = (G |/m where G the gradient for the εth weight matrix, i.e.,
(4)
(
(4) t-1 (0
(158)
41

---

Under review as a conference paper at ICLR 2021
In the following lemma, derive the upper-bound on the difference of the gradient passing from the 7th to (e 1)th layer given the same inputs D ((++1) for the éth layer in SGCN++ is defined as
H ((--1), where the backward propagation
V 70 [D] (e+1) H ( w(
(159)
(0) + [ T (++1) ( o' (2!")) w/
L (D (6+1) oo (Z ()
and the backward propagation for the 9th layer in FullGGN is defined
(D
(e+1) H. (L-1) w! L (D! o' (Z, w!
(160)
Lemma 7. Let suppose t E {EK 1..... EK + K). The upper-bound the difference of the gradient with respect to the input node embedding matrix at the fth graph convolutional layer given the same input D (e+1) and H
,....
(1) is defined
B H
o (D
(e+1)
H (1-1)
w(O) VHI
(++1) H (D)
)1
EK+K
(161)
Σ
④ B (
L x =
EK+1
Proof. To simplify the presentation, let us denote D = vn 70 D (D. [+1) H(--1) w Then, by definition have
D = 1 [
(
(D)
((+1)
o (0)
((+1)
o
(162)
Therefore, we know that
L DD (+1)
w (
D 1
= IIL
LT (++1)
o / (Z
w
- - D ( D (1)
-L LL
00 (Z
(163)
(Ai)
(e+1) (e+1) T (D t-1 oodd00oo0 o w)
D t-1
| +
D
D
Ag
Aa
2(A,,, A2) - A3) - 2(A2,
conditioon on F, on both side, and using the fact that Taking expectation condition
EDD ( 1
EE[L (D) oo w! [D
(D) (Z
(164)
L L ((1))
(D((t)) od d(( ))) (
L LD (++1) of (Z
the following incquality holds
E T (D.
|L (D. (e+1) o W!() - -D - (165)
| (D) (111) (Z) w(( -L LTD oZ
Then, taking expectation over Fú, have
B||l (D) ((+1) o' Z 11)
1 D 1 [|]]
E L (D (E+1)
t-1 o ( - (0)
+ D - 12
(166)
- [ D (+1) ( w!O
(D(+t)) o o (W
42

---

Under review as a conference paper at ICLR 2021
Let suppose t {EK + +,, I...., EK K). Then can denote t EK + k for some à K such that
B (D, (8-1) o (Z, (8)
D (2)
(D
(t-1) (7
EK+A -D
EK+* BK+k
EK+K
<BILL (D
((+1) FK
od
(8)
W
(1)
- D
Σ
-
EK
EK
EK
FK
t=EK-
-IILT (+-))
(Z W L LD
(++1)
(Z
:-1
[B[ (D (2+1)
o'
[EEK] EK
W( D
E
EK+K
E
(D, (E+1) o oZ w TD
(e+1) 0o (2 ))
1=EK+1
B T (D) (++1)
( oo (Z, w/ L (D oo 00 w(
(167)
Knowing that we are taking full-batch gradient descent when (1 mod K) = 0, we have
B
(
(E+1)
(Z
w/
1
ERAK
Σ
[BIL
1 [1] I
EDI + o' - (Z.
EEK+
(168)
Let take closer look at term (B).
B[I[DD (8-1) oZ w)(()) - (D (++1)
1-1 o o ( (1) w
S 3C2BB |D (E+1) (6+1) 3 BB 1 (
+30
W
t-
--1
(C2)
(169)
For term (Ci) by definition know
[ (+1) (++1) 12
= [.
(4,
w((())
[ (D
(D ((+2) (2 (1+1) w((+))
||
3|
[L (D)
(L
(E+1)
(170)
3|
((t11)
(D
w (+1)
+ 3|
((+2)
((+1)
(LT ((+2) (1-1) w i-1
(++1))
s SOIDD D
1 | (1 ((+1)
w (1
By induction, we have
(D((t))
-
D "
[ (1+1) (L+1)) || 0
(++) 2 (1 + 0(12 - -
(D:)
(Dx)
((++)) (2) +... !! - w!
(171)
43

---

Under review as a conference paper at ICLR 2021
For term (D1) we have
[ (T(
(1++1)
|1
)L(0.)
ac(0z-1)
2-
ow!!
aw
(172)
loss C2 1
2,
For term (D2) we have
|| (++1)
2 ((+1)
| s LH, w (++1) - LH W
(41)
1-1
S C. H (((1)) -H W (++1)
H (0) ((+1)
-
H (
- BC
1 BI
w(+)))
(+1))
H
w
1-1
s 2C, BL Bi [ - + 20 B
BB W
(++1)
w
(++1)
(173)
By induction we have
IZ (z((+1) - z (2+1) < (6+1 w -/ O -
(174)
For term (C2) by definition we have
E Z (6) 7 | E
W W -
<B Bl
C (t-1) + 2 W
w
(175)
By induction we have
E1|Ž!0 - < O(BIl
+... - (176)
Plugging (Di), (D2) back (Ci) and (C2), (Ca) back to (B), we have
(D ((+1) o oZ(( w!
d - (e+1)
Ii wh 2 + O (2) w(1
-D
EK+K
(O(ElI (1) w
(177)
t= EK+1
BK*K
Σ
20
B||| [|]]Vž((
BI
Then plugging term (B) back to Eq. 169 conclude the proof.
Using the previous lemma, provide the upper-bound of Eq. 154, which is of the three key factors that affect the mean-square error of stochastic gradient at the Eth layer.
Lemma 8. Let suppose I. E {EK + 1,..., F.K + K). The upper-bound the difference of the gradient with respect to the input node embedding matrix at the lth graph convolutional layer given the same inpu D (2+1) but different input H. (e-)) H is defined
B
F0 DD.
H
VHJ D
(++1) H
(f-1) W.
(
EK+K
(178)
Σ
20
x =
EK
44

---

Under review as a conference paper at ICLR 2021
Proof. For the gradient w.r.t. the node embedding matrices, we have
B H F (D)
(E+1)
H w
VH
(D.
(++1)
H 1)
Elll D
(++)
w! T (++1) o (2
(I T (D. (+1))
|ا
2Bll [D 1 T (D) (++1) o (2, w!
(--1) .d o d(
(1)
LL (D (D( oo o' (Z
(4)
+ 2 2lll (T( ((+1) oo (Z)") (Z w!) (D (1+1) o oz () [P:
(179)
Let first take closer look at term (A). Let suppose t E {EK + 1,..,, EK + K), where E = t mod K is the current epoch number and K is the inner-loop size. By the previous lemma, term (A) can be bounded by
EK+K E
(A) S
Fo
(
| x ElJ V(
(180)
I=EK+1
Then take closer look at term (B).
[I (D (E+1)
(z
w!
IL. T.T(D) (D
(++1)
o o )) [1
S BI 122 E 7 Z 0
(181)
(C)
The term (C') can be decomposed as
[|| Z (0)
- (4)
[
|1 LH
(182)
s 2 (0)
LA.
| + 2RIILL
w
LH
w
s 2E 7
LH.
|| + 2B
By
C
Z (6-1)
By induction, we have
El Z
z
+ 12 <
O
E
LH.
12)
(0..+0
-LX 1PI) (183)
The upper-bound for term (D) is similar to one we have in the proof of Lemma 4.
|LH ( (6-1)
w( (
2
=I[LH, w!
LH
1-1 1 + LA
--1 -
- - - -
-IILH w(O
LH 1-1
LA
1-1 2 + |7 7 (184)
+ +L
(- w
()) 2(L
w ( -00
-2 -0, LH, 1 (2)
Taking expectation condition on F, and using
E (6)
7 1 LH.
()) w
LH!
(185)
45

---

Under review as a conference paper at ICLR 2021
have
E LH
w 2
[
ILH
(0)
+
ILH
w!
LH.
(-1 IP (186)
Take expectation over Fz have
E L
W ス (
|
(0) -|LHH w 1-1 - 1122 -
E112
-
ILH
w!
LH
1-
0
l (187)
t-1
Let suppose & E {EK + ..... EK K}. Then we can denote t EK + k for some k K such that
E LI
(
= =E||LH
(1)
70 +
EK +
EK+*
EK+K
|| w
((-1)
(2)
(0)
||
Σ
E
LLH 6 w!
LH
EK
EK
EK+K Σ
E 1
L
w
222
t=EK+1
(188)
Plugging (D) to (C) and (C) (B), have
EK K
(B) < Σ Σ
- x B( - w21111
(189)
F= E
Combing with term (A) we have
Bill V R 700 (D (++1) ((-1)
! v ( ( (++1) H (L-1) )
1 ÆK+K
( | x E|| w 2 w(1ll)) 11
(190)
1=1 t= EK+1
EK+K Σ
[21 | x =(I
t=E
In the following lemma, we derive the upper-bound on the difference of the gradient with respect to the weight matrix each graph convolutional layer. Suppose the input node embedding matrix for the Åth GCN layer is defined as H (e-1) the gradient calculated for the εth weight matrix in + defined as
Vw Fw (D
(191)
6 LL "I (D( (et)) of (Z,)))
((+1) o (T (Z()))
and the backward propagation for the 9th layer in Ful1GCN is defined as
Vw ( (E+1) (((e1)) w = ( (e+)) oz (Z
(192)
Lemma 9. Let suppose t E {EK + +,..., EK + K). The upper-boune on the difference of the gradient with respect the lth graph convolutional layer given the input D (2+1) and H (((-1))
46

---

Under review as a conference paper at ICLR 2021
is defined as
Vw (6+1) H (()1)) w(9)
(D.
V "
w
EK+K
(193)
< Σ 70 [ |2
||ا
=EK+
Proof. To simplify the presentation, let us denote G.O)) by definition, we have
(w
(++1) , w(()). Then,
G = G I H
(DD
(1) oo -
""""
(+)) oo o' (2(,)) (194)
Therefore, we know that [LA D1
o' (2 1
[
(Z.)) - LH LH -
D
(Z
+ | با (D. ( o d(
o0
(7
6
_
(D
(41)
od
(195)
(A1)
LH "T (D(+))
(++1) o (Z 2 |
G
9
G
+ 2((A), A2) 2(A), Aa) - A(A2, Aa)
Taking expectation condition on Fz on both side, and using the fact that
B|G G [1) = ( (TTTTTT((S) (D. ((+1) oz - " (++1)
(z oZ
1-1
|L (D. (++1) d' (Z - (++1) o (Z(2))
(196)
we have the following inequality holds
B[][ (2)T (D(( (e+1)
(Z) ()) 2 [
LH ((-1) (D (6+1) o ( 2 Eाा.. - -
(197)
| tt- "J (D)
(+11) o - (<| o
o 2Z
Then, taking expectation over Fr, have
Elll [H (D. (1-1) d' ( (1 - ( [|] 2
< B|I[LA (t-
1-1 T DD (f+1) od 7 (0) w - 20, BII
E - -(
(198)
- [|| (D. (+1)) c o (7 | DD |- o od oZ
(+1)
47

---

Under review as a conference paper at ICLR 2021
Let suppose t {EK + 1, .,..., EK + K). Then we can denote t EK + k for some A K such
that
|
1++1)
(Z - [
E|| EK*
((+1)
(k)+ GO 1
EK*
EK+K
< [[
EK
EK
(++))
oo
(Z
(2)
EK
)) G
Σ
BII (2)
- G2
- -II[LA((1] (D.
(+1)
(Z( LH
(D
( (++1) oo (2
s B (D
((+1)
oo (2)
10)
D (4)
2
EK
EK
EK
EK
EK+K
Σ
1 (D) ((+1) 2Z L ( H ( (DD (t+1) o (2
Bill L
[D] (++1) o o (Z!"))
L (
(++) o o' (2
(199)
Knowing that we are taking full-batch gradient descent when (1 mod K)= = 0, we have
B
( (E+1)
(Z!)))
EK+K
Σ
[BIL 1 [1] I x
B (f-))
(D) 1-1 o oZ(2
(++1)
EK
(B)
(200)
Let take closer look at term (B).
B T (D. ; (++) (D,
o (Z(
- [ 11 (D (+11)
1-1 0
(2
< 3 + B B3 13 |
(t-1)
- - 333
C
El
(+1)
D
(4 [|]
(201)
(Ci)
(C2)
For term (Cz) by definition we know
|ID] (++1) D " ||
[L T (e-
(++1)
w ((111
(D
(4+2)
(++1)
((( |
(-1
3|
((+1)
(E+2)
(1-1) |/
LL.
w!
3| LL
w(
(LT
(--1)) w( |1
s O ((+)
(++2)
(2) 0 (+1)
(+1)
(6+1)
(202)
By induction, we have
[D] (++1)
D D
< O
[
(1+1) ((+)))
D D
| 0
(f+1) 2
O(1Z (4)
Z (1)
(D.)
(D²)
([ w
[ + ... O[|
(
w
(203)
For term (D,) we have
D (L+1) - D (1+1) ||
aC(O.) Oc
aw(L) (L)
ow t-1
(204)
1 02 12
(1) -
48

---

Under review as a conference paper at ICLR 2021
For term (D2) we have
[| Z 11+1)
| - C [H W (+1))
LH (
W ((+1)
Ca H ((+1) H
W. (6+1)
H
t-1
- 20 -H
H
s 20g BL By |z!
-
11 20 B B3
W ((+1)
w
(1+1)
[-1
(205)
By induction we have
(++1) (2 - [-1 | $ [(1
- 1-1 1. w
(206)
For term (C2) by definition we have
E 7 (
| = BlJi
W
;- w 1
< BW (f-1) -
BI w!
(207)
By induction we have
E[Iī - < O( w/
w( +
+ +B - (208)
Plugging (D)), (D2) to (C2) and C₁ (C2) to (B), we have
EK+K
E |H ( (D. (D( ((+1) od ( [2] S
Σ
E=
(209)
Using the previous lemma, we provide the upper-bound of Eq. 155, which is one of the three key factors that affect the mean-square error of stochastic gradient at the éth layer.
Lemma 10.
Let suppose E {EK + L,..,, + K}. The upper-bound on the difference of the
gradient with respect to the weight of the εth graph convolutional layer given the same input D (f+1) but different input H (f-1)
..! (t-1)
is defined as
Blll 70
D (e+1)
H
F( (D. (e-1) H (! w! |
EK+K
(210)
Σ
700
(
=
E
49

---

Under review as a conference paper at ICLR 2021
Proof. For the gradient w.r.t the weight matrices, have
V f (D.
G H
(E+1)
w! (
(
(
H. ( 1) w/
=Ell
DD. (++1 o' (Z.0)) - H (D) D 00 (2.21)
[L (+1))
|
s 2 laal 2 [ H 1 T (D. (D.ttt ((+t) oo oo (Z())) - [D (D od d
(LH (D. ((tt) o' ( (Z1 |
(4)
2 Bill (LU (D(t))) o (2,
(DD(t") oo oz( ( [1]]
(211)
Let first take closer look at term (A). Let suppose t E (EK + 1,.... EK + K), where E = t mod K is the current epoch number and K is the inner-loop size. By the previous lemma, term (A) can be bounded by
EK+K Σ
(A) S
(oulw)
(ollw"-WWii + ......O0 + + - -
w
(212)
(EK+1
Then take closer look at term (B).
(D. (+1)
(Z
[LH (1- 1/
(D. (++1) oz(z!")) 1
BB B
Z (2)
(213)
(C)
Let suppose t {EK + 1...., EK + K}.
Then we can denote
EK + l for some k: S K.
From Eq. 220, we know that
(C) = EBILLA
2 ||
EK*K
(214)
Σ
E
I
w(.112)
(=E+1
Plugging (C) back to (B), combing with (A) we have
B[llV w f' (D. (++1)
H
((-1)
w(O)) (w
DD.
(--1) H w(
EK+K
(215)
Σ
120 (70 E
|
E=
In the following lemma, provide the upper-bound of Eq. 153, which is one of the three key factors that affect the mean- mean-square error stochastic gradient at the éth layer.
Lemma 11. Let suppose t E (EK + +1, 1,...., FK + K}. Then the upper bound of
EK+K
E D (1+) -D.
D (1+1) <
Σ
n O [
[2
LI
E|IVVC((-)11111) (216)
50

---

Under review as a conference paper at ICLR 2021
Proof. By definition we have
[ (D(+t)
- (1--1)
- aā((.)
0C(0,)
0 FF
x Loss [ (4)
H
(217)
2 oss C2
(Z(40) z 1
(A)
Let take closer look at term (A):
1 " = =I (1) , w(
(1-1)) -
- (H w
+ 2 21 w(() 400 (H. (1-1) w())
=2 w(()) (H, (1 w
(218)
+ -
||
$ <211 - (-1))
{ w
20 By ((1-11 z 11
By induction, we have
17, - 2(
O (1) HH (1-1)
-(0) (H. (H())))((0))
- f(0
w())
(219)
t=1
Let suppose {EK + L,....,
EK + Kj. Then we can denote
EK + for some S K.
From Eq. 220, we know that
(
w(
(H)
w
EK+K
(220)
Σ
|
E
|
x B||
w 11))
E=EK+1
Therefore, we know
[D((+t)) D (1+1))
2
Liss 10
1
L E. FK+K K
ILIEI x W! 0 w(ilf)
(221)
1=EK+1
BK+K
Σ
n² O (
| ExxEEI [
EK+1
which conclude the proof.
Combing the upper-bound of Eq. 153, 154, 155, we provide the upper-bound of mean-suqare error of stochastic gradient in SGCN++.
Lemma 12. Let suppose t {EK+ + I,..., EK + K}. Then we can denote k < K such that
-EK + k for some
EK+K Σ
Ell Vī -C
120 E 1
|
E V((0--)))
(222)
E
51

---

Under review as a conference paper at ICLR 2021
Proof. From Lemma 6, we have
B G
-O ((++)) I?:))
+ O 7 (D 11+11 ( w()) V (1+1)) H H w(()))
+ +0 ((-1) (D (222) H V H w((+
+ VO 1 (D (++1) ((-1) w()) Vw (D ((+1)) ((-1), (()
(
(223)
Plugging the result from support lemmas, i.e., Lemma 8, 10, 11 and using the definition of stochastic gradient for all model parameters V (O,)) (
Jf=1
have
EK+K Σ
E|| Vā())) - v=(0,) <
O
P
=|I (224)
t=EK +1
1.2 REMAINING STEPS TOWARD THEOREM 3
By the smoothness of L(:)). we have
Lf 10++1 0:1
[(0r+1) ((),) (VL(0.), - 0₂)
L2 | ((
[(0)) (VL,)))) Vā(O.))
C(0,) - 7 (c.)) + " (V vī(0,)1 - 13 Lm?
(a) 2
2
(225)
where equality (a) is due to the fact 2(x, v) = || |1 || - u1 for any y.
Take expectation on both sides, we have
B[((Or+)) s - + 1 L
E V 2(0
2
2
(226)
By summing over t I,..., T where T' is the inner-loop size, we have
r
[(
2 [ + (BBII - (1 (1- [L_L,m)[[|V((0.)]1
<
2 -I(((r)))) + [[ - - (1-Lm) E v v2
17
t=1
(227)
where L(8*) is the global optimal solution.
52

---

Under review as a conference paper at ICLR 2021
Let us consider each inner-loop with E T/K and t {EK 1...., EK + K}. Using
Lemma 12 we have
T/K-1 EK K
7/K-1 EK+K
E F=
Σ
(())
(1 (1- (1-
Σ E
tm ER+1
F=
t= EK+1
T/K-1 KK t-EK
E
|L|
Σ
1=1
E=U B=E y=2
T/K-1
EK+K
- (1 Lim)
Σ
Σ
(228)
E=0
-=
T/K-1 EK-K
(n) K 0
E
| - (11-
Σ
Σ
E|I B
1=1
E=0
EK-1
T/K-1 EK*K
In A
bin - (1 - Lyn)
Σ
Σ
E=
EK+1
Notice that 1j
is root of equation 17 r (1 Iyn) 0. Therefore
Lyt 13 btn
have
El CCO 2 (BIC(0))
(229)
which implies
Σ
_ T V
(B(C(O)) B[C(0"))
(230)
T
J REPRODUCING EXPERIMENT RESULTS
To reproduce the results reported in the paper, we provide link for dataset download, the bash script to reproduce the experiment results, and jupyter notebook file for quick visualization and GPU utilization calculation. It is worth noting that due to the existence of randomness, the obtained results (e.g., loss curve) may be slightly different. However, it is not difficult find that the overall trend of loss curves and conclusions will remains the same.
This implementation is based *Py Torch using Python 3. We notice that Python 2 might results in a wrong gradient update, even for vanilla SGCNs.
Install dependencies:
create virtual environment
$ virtualenv env
$ source env/bin/activate
install dependencies
$ pip install -r requi
Experiments are produced on PPI, PPI-Large, Flickr, Reddit, and Yelp datasets. The utilized datasets can be downloaded from 'Google drive.
create folders that save experiment results and datasets $ mkdir /rresulss $ mkdir //atta please download the dataset and put them inside this folder
https::/drrie....ggle.....///// htt.://w...gog..eccccee///e//e// 15Ph0 htw.:///ooarrrt//oooooooooooonnn h00 usp-sharing
53

---

Under review as a conference paper at ICLR 2021
To reproduce the results, please run the following commands:
$ python train.py sample_method
'ladies' --dataset 'reddit'
$ python train.py S ample_method
'fastgen' --dataset 'reddit'
$ Python train.py sample_methhd
'graphsage' --dataset 'reddit'
$ python train.py sample_method
'vrgen' --dataset 'reddit'
$ python train.py sample_method
'graphsaint! --dataset 'reddit'
$ python train.py sample_method
'exact' --dataset 'reddit'
$ python train.py sample_method
'ladies' --dataset 'ppi' 'fastqen' --dataset 'ppi'
$ python train.py ample_method
$ python train.py sample_method
'graphsage' --dataset 'ppi'
$ python train.py s sample_method
'vrgen' --dataset 'ppi' 'graphsaint' --dataset 'ppi' 'exact' --dataset 'ppi'
$ python train.py
$ python train.py ample_method
$ python train.py S sample_method
'ladies' dataset 'flickr' 'fastgen' --dataset 'flickr'
$ python train.py sample_method
$ python train.py sample_method
'graphsage' --dataset 'flickr'
$ python train.py sample_method
'vrgen' --dataset 'flickr' 'graphsaint --dataset 'flickr'
$ python train.py s ample_methoo
$ python train.py
amppeemethod
'exact' --dataset 'flickr'
$ python train.py sample_method
'ladies' --dataset 'ppi-large'
$ python train.py sample_method
'fastgcn' --dataset 'ppi-large'
$ python train.py sample_nethod
'graphsage' --dataset 'ppi-large'
$ python train.py sample_method
'vrgen' --dataset 'ppi-large'
$ python train.py
hoo $ python train.py amppe_meehoe
'graphsaint --dataset 'ppi-large'
'exact' --dataset 'ppi-large'
$ python train.py sample_method
'ladies' --dataset 'yelp' 'fastgen' --dataset 'yelp' $ python train.py sample_method 'graphsage' --dataset 'yelp'
$ python train.py sample_metho
$ python train.py sample_method 'vrgen' --dataset 'yelp'
$ python train.py sample_methoo 'graphsaint' --dataset 'yelp'
$ python train.py -- --ample_method 'exact' --dataset 'yelp'
54


Reviews:
[{"review_id": "NakK4-wffNP", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "##########################################################################\n\nSummary: \nThis paper presents a novel variance reduction method which can adapt to any sampling-based GCN methods (inductive GCNs). The paper draws the idea from VRGCN that integrates the historical latent representations of nodes computed with full Laplacian to approximate the that computed with sampled sparse Laplacian. The variance reduction is implemented on both node embedding approximation, as well as layer-wise gradient computation in back-propagation. The resulting algorithms lead to faster convergence rate.\n\n##########################################################################\n\nReasons for score:\nOverall, I vote for accepting. The proposed variance reduction techniques can successfully accelerate convergence of any sampling method, according to the experiments, while also enjoys theoretical guarantee. Yet the novelty of the proposed SGCN+/SGCN++ algorithms themselves is a little limited to some extent. \n\n##########################################################################\n\nPros:\n1. The authors introduced a doubly variance reduction which can effectively reduce the node approximation variance the layer-wise gradient variance of the existing sampling based GCN methods and accelerate convergence.\n2. This paper also provides thorough theoretical analysis and convergence guarantee of the proposed algorithms.\n3. The  authors have conducted comprehensive experiments using a variety of sampling-based GCN methods as building blocks. The quantitative results clearly demonstrate the effectiveness of the proposed algorithms. The authors also provide detailed empirical analysis on the training time / GPU memory usage of the proposed method.\n\n##########################################################################\n\nCons:\n1. To better illustrate the idea the variance reduction, the authors could compare the proposed algorithms with a vanilla full-batch GCN (instead of the mini-batch training version Exact used in this paper, and it could be evaluated on smaller datasets like Cora). This baseline may serve as a theoretical upper bound of SGCN+/SGCN++.\n2. The proposed SGCN+ and SGCN++ requires a full-batch forward/backward computation every k step. As the authors suggest, this might hinder the scalability of SGCN++ on extremely large graphs. The authors hence propose a variant of SGCN++ which applies a large-batch approximation. The authors could also provide  an alternative version of SGCN+ without full-batch that only reduces the zeroth-order variance, and evaluate how SGCN+ without full-batch snapshot computation would impact on the zeroth-order and first-order variance.\n3. Since the snapshot gap K serve as a budge hyper parameter balances between training speed and quality of variance reduction. As the model converges w.r.t. increasing number of epochs, I would like to know whether we can dynamically increase K during the training process to obtain some speed boost.\n4. The paper is well-written in general. However, there is still some typos. For example, in Eq. (1) and Eq. (2), in the computation of gradient G_t^(l), the superscript in D_t should be (l+1) instead of (l), since we require the gradient of the loss w.r.t. the upper layer.\n\n#########################################################################\n\nQuestions during rebuttal period: \nPlease address and clarify the cons above\n\n#########################################################################\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "Mtqp-HbQYf-", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper provides a convergence analysis for GCN training and proposes two variance reduction algorithms to ensure convergence. \n\nThe authors first draw a connection between sampling-based GCN training and compositional optimization, and they divide the error of gradient into a bias term and a variance term. This is actually first established in [Cong et al. Minimal Variance Sampling with Provable Guarantees for Fast Training of Graph Neural Networks]. Although there is some inaccuracy in their formulation, the authors should clarify that this observation is first made by previous work. \n\nBased on the variance and bias decoupling, the authors then establish a convergence result for biased SGD (Theorem 1). The analysis itself is trivial and is not new. The authors should consider adding some references. \n\nThe authors then propose an algorithm to reduce the bias term (Algorithm 1). Again, this is not new. It is exactly the nested SPIDER method for solving multi-level stochastic compositional optimization [Zhang et al. Multi-Level Composite Stochastic Optimization via Nested Variance Reduction]. Although the paper is the first to apply the algorithm to GCN training (as far as I know), it is important to have correct references to clarify the contribution. Also, it is worth noting that the algorithm is based on the assumption that the full neighbor aggregation result is available and the storage of Z and H is possible. I doubt the practicality of the algorithm for training on large graphs. \n\nThe experimental results validate that the algorithm works for GCN training. However, I notice that some values in Table 1 for the baseline (SGCN) are not the best results reported in other papers. For example, GraphSAINT achieves 96% on Reddit graph, but the authors report 93.68% as the baseline. This makes me suspect that the good results reported in the paper may be due to some parameter tuning instead of the algorithm itself. \n\n\nA minor point: \nThe paper says “the data points in these works are sampled independently, but the data points (nodes) in SGCN are sampled node- or layer-dependent according to the graph structure.”\nThe data samples in compositional optimization do not need to be independent. The only requirement is that the stochastic function with the data samples is unbiased. \n", "extra_scores": {}, "extra_text": {}}, {"review_id": "Hdq09Ra3jhW", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper studies the convergence of stochastic training methods for graph neural networks. Here, this paper views GNN as a compositional optimization problem. Then, to reduce the variance incurred by the neighbor sampling, this paper uses SPIDER to reduce the variance to accelerate the convergence speed. It provides theoretical convergence analysis for SPIDER used on GNNs, showing that the proposed method has a better convergence rate compared with the traditional gradient descent method. At last, this paper conducts experiments to verify the proposed algorithm. \n\n1. Overall, this is a new application of SPIDER on GNNs. It shows how to bound the variance under the setting of GNNs, which is then used for the convergence proof. However, some assumptions are too strong, such as assumption 3. With these strong assumptions, the theoretical analysis is simplified too much.\n\n2. Some parts are not very clear. In eq.(3), this paper claims that the bias term is mainly caused by node embedding approximation in the forward pass while the variance term is caused by gradient varinace in the backward pass. It is very confusing. The node embedding approximation affects both the forward and backward pass. Why does it only affect the forward pass?\n\n3. For figure 3, the exact sampling method uses all neighbors. Therefore, it is actually the full gradient desenct method. Why does there exist variance? \n\n4. The most weak point of this paper is that the theory is NOT consistent with experiments. In detail, the theorems study the convergence rate of (variance-reduced) SGD. However, in experiments, this paper uses Adam. Thus, the experimental results cannot support the claim of theories. \n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "6gteQiaJymJ", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Node sampling is a crucial point in making GCNs efficient. While several sampling methods have been proposed previously, the theorectical convergence analysis is still lacking. This paper finds that the convergence speed is related to not only the function approximation error but also the layer-gradient error. Based on this finding, the authors suggest to take historical hidden features and historical gradients to do doubly variance reduction. Experiments are done on 5 datasets for 7 baseline sampling-based GCNs.\n\nPros:\n\n1. The core contribution of this paper lies in Theorem 1, which reveals the relationship between convergence speed and the function approximation error and the layer-gradient error. \n\n2. The idea of doubly variance reduction is reasonable.\n\nCons:\n\nThe biggest weakness of the proposed method is that it requires to compute snapshot features and gradients over all nodes (Line 5, Alg. 1 & Line 5 Alg. 2) before the sampling process. As this paper aims at enhancing sampling based GCNs, we should assumes no computation access/memory of performing full GCN. The authors have provided the related analyses in the appendix. It will be better if the experiments without full-batch shapshot are added in Table 1, as such we can check how it influences the final performance and if certain approximation will work well. \n\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper provides variance reduction techniques for GCN training. When training a GCN it is common to sample nodes as in SGD, but also subsample the nodes’ neighbors, due to computational reasons. The entire mechanism introduces both bias and variance to the gradient estimation. The authors decompose the gradient estimate into its variance and bias error, allowing them to apply more targeted variance (and bias) reduction techniques.

The results and improvement over existing GCN methods seem to be solid. The main weakness of the paper is its novelty. As pointed out in the reviews the techniques seem to be quite close to papers [5],[11] (referring to the authors posted list).
It therefore boils down to the question of whether the authors simply applied existing techniques, achieving a better implementation than previous art, or did they develop a truly new algorithm that will encourage further research and deepen the understanding of GCNs. Given the decisive opinions of reviewers 1 and 4, that remained after taking the response into account, I tend to believe that the improvement provided here is either too incremental or not stated in a crisp enough manner in order to be published in its current form

Decision: Reject