summaries: 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.  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. This 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.
abstractText: 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 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 (zeroth-order variance) during forward propagation and layerwise-gradient variance (first-order variance) during backward propagation. We theoretically analyze the convergence of the proposed schema and show that it enjoys an 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.
references.authors: Léon Bottou Frank E Curtis Jorge Nocedal
references.title: Optimization methods for large-scale machine learning
references.venue: Siam Review,
references.authors: Tianle Cai Shengjie Luo Keyulu Xu Di He Tie-yan Liu Liwei Wang
references.title: Graphnorm: A principled approach to accelerating graph neural network training
references.venue: arXiv preprint arXiv:2009.03294,
references.authors: Jianfei Chen Jun Zhu Le Song
references.title: Stochastic training of graph convolutional networks with variance reduction
references.venue: arXiv preprint arXiv:1710.10568,
references.authors: Jie Chen Ronny Luss
references.title: Stochastic gradient descent with biased but consistent gradient estimators
references.venue: arXiv preprint arXiv:1807.11880,
references.authors: Jie Chen Tengfei Ma Cao Xiao
references.title: Fastgcn: fast learning with graph convolutional networks via importance sampling
references.venue: arXiv preprint arXiv:1801.10247,
references.authors: Minshuo Chen Xingguo Li Tuo Zhao
references.title: On generalization bounds of a family of recurrent neural networks
references.venue: arXiv preprint arXiv:1910.12947,
references.authors: Tianyi Chen Yuejiao Sun Wotao Yin
references.title: Solving stochastic compositional optimization is nearly as easy as solving stochastic optimization
references.venue: arXiv preprint arXiv:2008.10847,
references.authors: Wei-Lin Chiang Xuanqing Liu Si Si Yang Li Samy Bengio Cho-Jui Hsieh
references.title: Cluster-gcn: An efficient algorithm for training deep and large graph convolutional networks
references.authors: Weilin Cong Rana Forsati Mahmut Kandemir Mehrdad Mahdavi
references.title: Minimal variance sampling with provable guarantees for fast training of graph neural networks
references.venue: In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining,
references.authors: Vijay Prakash Dwivedi Chaitanya K Joshi Thomas Laurent Yoshua Bengio Xavier Bresson
references.title: Benchmarking graph neural networks
references.venue: arXiv preprint arXiv:2003.00982,
references.authors: Cong Fang Chris Junchi Li Zhouchen Lin Tong Zhang
references.title: Spider: Near-optimal non-convex optimization via stochastic path-integrated differential estimator
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Alex Fout Jonathon Byrd Basir Shariat Asa Ben-Hur
references.title: Protein interface prediction using graph convolutional networks
references.venue: In NIPS,
references.authors: Vikas K Garg Stefanie Jegelka Tommi Jaakkola
references.title: Generalization and representational limits of graph neural networks
references.venue: arXiv preprint arXiv:2002.06157,
references.authors: David Gross
references.title: Recovering low-rank matrices from few coefficients in any basis
references.venue: IEEE Transactions on Information Theory,
references.authors: Will Hamilton Zhitao Ying Jure Leskovec
references.title: Inductive representation learning on large graphs
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Yifan Hu Siqi Zhang Xin Chen Niao He
references.title: Biased stochastic gradient descent for conditional stochastic optimization
references.venue: arXiv preprint arXiv:2002.10790,
references.authors: Wenbing Huang Tong Zhang Yu Rong Junzhou Huang
references.title: Adaptive sampling towards fast graph representation learning
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Thomas N Kipf Max Welling
references.title: Semi-supervised classification with graph convolutional networks
references.venue: arXiv preprint arXiv:1609.02907,
references.authors: Sashank J Reddi Ahmed Hefny Suvrit Sra Barnabás Póczos Alex Smola
references.title: Stochastic variance reduction for nonconvex optimization
references.venue: In International conference on machine learning,
references.authors: Ryoma Sato Makoto Yamada Hisashi Kashima
references.title: Constant time graph neural networks, 2020
references.venue: URL https://openreview.net/forum?id=rkgKW64FPH
references.authors: Michael Schlichtkrull Thomas N Kipf Peter Bloem Rianne Van Den Berg Ivan Titov Max Welling
references.title: Modeling relational data with graph convolutional networks
references.venue: In European Semantic Web Conference,
references.authors: Quoc Tran-Dinh Deyi Liu Lam M Nguyen
references.title: Hybrid variance-reduced sgd algorithms for nonconvex-concave minimax problems
references.venue: arXiv preprint arXiv:2006.15266,
references.authors: Quan Wang Zhendong Mao Bin Wang Li Guo
references.title: Knowledge graph embedding: A survey of approaches and applications
references.venue: IEEE Transactions on Knowledge and Data Engineering,
references.authors: Keyulu Xu Weihua Hu Jure Leskovec Stefanie Jegelka
references.title: How powerful are graph neural networks
references.venue: arXiv preprint arXiv:1810.00826,
references.authors: Shuoguang Yang Mengdi Wang Ethan X Fang
references.title: Multilevel stochastic gradient methods for nested composition optimization
references.venue: SIAM Journal on Optimization,
references.authors: Rex Ying Ruining He Kaifeng Chen Pong Eksombatchai William L Hamilton Jure Leskovec
references.title: Graph convolutional neural networks for web-scale recommender systems
references.venue: In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining,
references.authors: Hanqing Zeng Hongkuan Zhou Ajitesh Srivastava Rajgopal Kannan Viktor Prasanna
references.title: Graphsaint: Graph sampling based inductive learning method
references.authors: Junyu Zhang Lin Xiao
references.title: A composite randomized incremental gradient method
references.venue: In International Conference on Machine Learning,
references.authors: Junyu Zhang Lin Xiao
references.title: Multi-level composite stochastic optimization via nested variance reduction
references.venue: arXiv preprint arXiv:1908.11468,
references.authors: Junyu Zhang Lin Xiao
references.title: A stochastic composite gradient method with incremental variance reduction
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Yuyu Zhang Xinshi Chen Yuan Yang Arun Ramamurthy Bo Li Yuan Qi Le Song
references.title: Efficient probabilistic logic reasoning with graph neural networks
references.venue: arXiv preprint arXiv:2001.11850,
references.authors: Difan Zou Ziniu Hu Yewen Wang Song Jiang Yizhou Sun Quanquan Gu
references.title: Layer-dependent importance sampling for training deep and large graph convolutional networks
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Chen
references.title: 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
references.authors: Cai
references.title: tion”), which have been proven to impact the GCN performance
references.authors: FastGCN Chen
references.title: 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
references.authors: GraphSAINT Zeng
references.title: 2019) can be viewed as a special case of layer-wise sampling algorithm FastGCN by restricting the nodes sampled at the 1-st to (L−1)th layer the same as the nodes sampled at the Lth layer. However, GraphSAINT requires a significant large mini-batch size compared to other layer-wise sampling methods
references.authors: The intuition of snapshot step large-batch approximation stems from matrix Bernstein inequality Gross
references.title: More specifically, suppose given G̃i ∈ Rd×d be the stochastic gradient computed by using the ith node with Exact sampling (all neighbors are used to calculate the exact node embeddings)
references.venue: Suppose the different between G̃i and full-gradient
sections.text: 1 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), supervised graph classification (Xu et al., 2018), protein interface prediction (Fout et al., 2017), and knowledge graph (Schlichtkrull et al., 2018; Wang et al., 2017). However, most works on GCNs focus 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-field 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 et al., 2017; Ying et 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 a few key issues. First, the theoretical understanding of samplingbased 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 (e.g., current gradient or representation of nodes), finding the optimal adaptive sampling distribution is computationally inadmissible, as 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 use the historical node
  embeddings and the historical layerwise gradient of each graph convolution layer as control variants. 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 structure 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 (firstorder 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 a significant reduction in variance; hence faster convergence rate and better test accuracy. We can also see that both zerothorder 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-based GCN training (SGCN) with a nonasymptotic 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 employing zeroth-order variance reduction to node embedding approximation (dubbed as SGCN+), which explains why VRGCN (Chen et al., 2017) enjoys a 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 a generic and efficient doubly variance reduction schema (SGCN++). SGCN++ can be integrated with different sampling-based methods to significantly reduce both zeroth- and first-order variance, and resulting in a 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/ √ T ). We empirically verify SGCN++
  through various experiments on several real-world datasets and different sampling methods, where it demonstrates significant improvements over the original sampling methods.
  2 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’ representations in the memory. This is the key bottleneck that hinders the scalability of full-batch GCN training. To overcome this issues, sampling-based 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, a 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 convolutional operation on history activation during each forward propagation, which is computationally expensive. Another direction applies layerwise importance sampling to reduce variance. For example, FastGCN (Chen et al., 2018) independently sample a constant number of nodes in all layers using importance sampling. However, the sampled nodes are too sparse to achieve high accuracy.
  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 expensive sampling algorithm. In addition, subgraph sampling methods such as 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 understanding of the convergence for SGCNs training method is still limited. VRGCN provides a convergence analysis under a strong assumption that the stochastic gradient due to sampling is unbiased and achieved a convergence rate of O(1/ √ T ). 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 exponentially fast with respect to the sample size, and results in the same convergence rate as unbiased ones, i.e., O(1/ √ T ). 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, an objective function that only the first composite layer has the trainable parameters is considered, where the output of the lower-level function is 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 a 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-graph 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 setting. Given an undirected graph G = (V, E) with N = |V| and |E| edges and the adjacency matrix A ∈ {0, 1}N×N , we assume that each node is associated with a feature vector xi ∈ Rd and label yi. We use X = [x1, . . . ,xN ] ∈ RN×d and y = [y1, . . . , yN ] ∈ RN to denote the node feature matrix and label vector, respectively. The Laplacian matrix is calculated as L = D−1/2AD−1/2 or L = D−1A where D ∈ RN×N is the degree matrix. We use θ = {W(1), . . . ,W(L)} to denote the stacked weight parameters of a L-layer GCN. The training of full-batch GCN (FullGCN) as an empirical risk minimization problem aims at minimizing the loss L(θ) over all training data
  L(θ) = 1N ∑N i=1 Loss(h (L) i , yi), H (L) = σ ( L . . . σ ( Lσ ( LXW(1)︸ ︷︷ ︸
  Z(1)
  ) W(2) ) . . .W(L) ) where h(`)i is the ith row of node embedding matrix H
  (`) = σ(Z(`)),Z(`) = LH(`−1)W(`) that corresponds to embedding of ith node at `th layer (hop), Loss(·, ·) is the loss function (e.g., crossentropy 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). Sampling-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 L̃(`) for each layer with supp(L̃(`)) supp(L), and perform forward and backward propagation only based on the sampled Laplacian matrices. The sparse Laplacian matrix construction algorithms can be roughly
  classified as 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), we sample a mini-batch of nodes VB ⊆ V from all nodes with size B = |VB|, and construct the set of sparser Laplacian matrices {L̃(`)}L`=1 based on nodes sampled at each layer and compute the stochastic gradient to update parameters as 1
  ∇L̃(θ) = 1B ∑ i∈VB ∇Loss(h̃ (L) i , yi), H̃ (L) = σ ( L̃(L) . . . σ ( L̃(2)σ ( L̃(1)XW(1)︸ ︷︷ ︸
  Z̃(1)
  ) W(2) ) . . .W(L) ) Key challenges. Compared to vanilla SGD, the key challenge of theoretical understanding for SGCN training is the biasedness of stochastic gradient due to sampling of nodes at inner layers. Let denote FullGCN’s full-batch gradient as ∇L(θ) = {G(`) = ∂L(θ)
  ∂W(`) }L`=1 and SGCN’s stochastic
  gradient as ∇L̃(θ) = {G̃(`) = ∂L̃(θ) ∂W(`) }L`=1. By the chain rule, we can compute the full-batch gradient G(`)t w.r.t. the `th layer weight matrix W (`) as
  G (`) t = [LH (`−1) t ] >(D(`+1)t ◦ σ′(Z(`)t )), D(`)t = L>(D(`+1)t ◦ σ′(Z(`)t ))W(`)t , D(L+1)t = ∂L(θt)∂H(L) (1) and compute stochastic gradient G̃(`)t utilized in SGCN for the `th layer w.r.t. W (`) as
  G̃ (`) t = [L̃ (`)H̃ (`−1) t ] >(D̃(`+1)t ◦ σ′(Z̃(`)t )), D̃(`)t = [L̃(`)]>(D̃(`+1)t ◦ σ′(Z̃(`)t ))W(`)t , D̃(L+1)t = ∂L̃(θt)∂H̃(L) (2) For any layer ` ∈ [L], the stochastic gradient G̃(`)t is a biased estimator of full-batch gradient G (`) t , as it is computed from H(L)t and Z (`) t , which are not available in SGCN since H̃ (L) t and Z̃ (`) t 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/ √ T ), 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 a tight analysis without the aforementioned strong assumptions on the stochastic gradient. Our analysis is inspired by the bias and variance decomposition of the mean-square 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̃(θ)−∇L(θ)‖2F] = E[‖E[∇L̃(θ)]−∇L(θ)‖ 2 F]︸ ︷︷ ︸
  Bias E[‖b‖2F]
  +E[‖∇L̃(θ)− E[∇L̃(θ)]‖2F]︸ ︷︷ ︸ Variance E[‖n‖2F]
  (3)
  where the bias terms E[‖b‖2F] is mainly due to the node embedding approximation variance (zerothorder variance) during forward propagation and the variance term E[‖n‖2F] is mainly due to the layerwise gradient variance (first-order variance) during backward propagation. Before proceeding to analysis, we make the following standard assumptions on the Lipschitz-continuity and smoothness of the loss function Loss(·, ·) and activation function σ(·).
  Assumption 1. The loss function Loss(·, ·) is Closs-Lipschitz continuous and Lloss-smoothness w.r.t. to the input node embedding vector, i.e., ‖Loss(h(L), y)− Loss(h′(L), y)‖2 ≤ Closs‖h(L)−h′(L)‖2 and ‖∇Loss(h(L), y)−∇Loss(h′(L), y)‖2 ≤ Lloss‖h(L) − h′(L)‖2.
  Assumption 2. The activation function σ(·) is Cσ-Lipschitz continuous and Lσ-smoothness, i.e., ‖σ(z(`))− σ(z′(`))‖2 ≤ Cσ‖z(`) − z′(`)‖2 and ‖σ′(z(`))− σ′(z′(`))‖2 ≤ Lσ‖z(`) − z′(`)‖2.
  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 et al. (2020).
  Assumption 3. For any ` ∈ [L], the norm of weight matrices, Laplacian matrices, input node feature matrix are bounded: ‖W(`)‖
  F ≤ BW , ‖L̃(`)‖F ≤ BLA, ‖L‖F ≤ BLA, and ‖X‖F ≤ BH .
  1We use a tilde symbol ̃ for their stochastic form
  Proposition 1. For any ` ∈ [L], there exist constants BH and BD such that the norm of node embedding matrices and the gradient with respect to the input node embedding matrices satisfy ‖H(`)‖F ≤ BH , ‖H̃(`)‖F ≤ BH , ‖ ∂σ(LH(`−1)W(`)) ∂H(`−1) ‖F ≤ BD, and ‖ ∂σ(L̃(`)H̃(`−1)W(`)) ∂H̃(`−1) ‖ F ≤ BD.
  Before presenting the convergence of SGCN, we introduce the notation of propagation matrices {P(`)}L`=1, 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 in the practical training algorithms. By doing so, we can decompose the difference between L̃(`) and L as the summation of column-wise difference ‖L̃(`) −P(`)‖2F and row-wise difference ‖P(`) − L‖2F. 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[‖L̃(`) −P(`)‖2F] and rowwise difference E[‖P(`) − L‖2F] which can significantly impact the convergence of SGCN. Theorem 1 (Convergence of SGCN). Suppose Assumptions 1, 2, 3 hold and apply SGCN with learning rate chosen as η = min{1/LF, 1/ √ T} where LF is the smoothness constant. Let ∆n and ∆b denote the upper bound on the variance and bias of stochastic gradients as:
  ∆n = L∑ `=1 O(E[‖L̃(`) −P(`)‖2 F ]) +O(E[‖P(`) − L‖2F]), ∆b = L∑ `=1 O(E[‖P(`) − L‖2 F ]) (4)
  Then, the output of SGCN satisfies
  min t∈[T ]
  E[‖∇L(θt)‖2F] ≤ 2(L(θ1)− L(θ?))√ T + LF∆n√ T + ∆b. (5)
  The exact value of key parameters LF, ∆n, and ∆b are computed in Lemma 1, Lemma 2, and Lemma 3 respectively and can be found in Appendix G. Theorem 1 implies that after T iterations the gradient norm of SGCN is at most O(∆n/ √ T ) + ∆b, which suffers from a constant residual error ∆b that is not decreasing as the number of iterations T increases. Without the bias2 we recover the convergence of vanilla SGD. Of course, this type of convergence is only useful if ∆b and ∆n are small enough. We note that existing SGCN algorithms propose to reduce ∆b by increasing the number of neighbors sampled at each layer (e.g., GraphSAGE), or applying importance sampling (e.g., FastGCN, LADIES and GraphSAINT).
  4 SGCN+: ZEROTH-ORDER VARIANCE REDUCTION
  An important question to answer is: can we eliminate the residual error without using all neighbors during forward-propagation? 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 to estimate the true node embeddings. More specifically, the graph convolution in VRGCN is defined as H̃(`)t = σ ( LH̃ (`−1) t−1 W (`) + L̃(`)(H̃ (`−1) t − H̃ (`−1) t−1 )W (`) ) . Taking advantage 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 to 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.
  To fill in these gaps, we 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 at the regular steps. At the snapshot step (t mod K = 0), a full Laplacian matrix is utilized:
  Z (`) t = LH (`−1) t W (`) t , H (`) t = σ(Z (`) t ), Z̃ (`) t ← Z (`) t (6)
  During the regular steps (t mod K 6= 0), the sampled Laplacian matrix is utilized:
  Z̃ (`) t = Z̃ (`) t−1 + L̃ (`)H̃ (`−1) t W (`) t − L̃(`)H̃ (`−1) t−1 W (`) t−1, H̃ (`) t = σ(Z̃ (`)) (7)
  2We have ∆b = 0 if all neighbor are used to calculate the exact node embeddings, i.e., P(`) = L, ∀` ∈ [L].
  Algorithm 1 SGCN+: Zeroth-order variance reduction (Detailed version in Algorithm 4) 1: Input: Learning rate η > 0, snapshot gap K > 0 2: for t = 1, . . . , T do 3: if t mod K = 0 then 4: Calculate node embeddings using Eq. 6 5: Calculate full-batch gradient∇L(θt) as Eq. 1 and update as θt+1 = θt − η∇L(θt) 6: else 7: Calculate node embeddings using Eq. 7 8: Calculate stochastic gradient∇L̃(θt) as Eq. 2 and update as θt+1 = θt − η∇L̃(θt) 9: end if 10: end for 11: Output: Model with parameter θT+1
  Algorithm 2 SGCN++: Doubly variance reduction (Detailed version in Algorithm 5) 1: Input: Learning rate η > 0, snapshot gap K > 0 2: for t = 1, . . . , T do 3: if t mod K = 0 then 4: Calculate node embeddings using Eq. 6 5: Calculate full-batch gradient∇L(θt) usng Eq. 1 and update as θt+1 = θt − η∇L(θt) 6: Save the per layerwise gradient G̃(`)t ← G (`) t , D̃ (`) t ← D (`) t , ∀` ∈ [L]
  7: else 8: Calculate node embeddings using Eq. 7 9: Calculate stochastic gradient∇L̃(θt) usng Eq. 10 and update as θt+1 = θt − η∇L̃(θt)
  10: end if 11: end for 12: Output: Model with parameter θT+1
  Comparing with VRGCN, the proposed SGCN+ only requires one full Laplacian graph convolution operation every K iterations, where K > 0 is an additional parameter to be tuned.
  In the following theorem, we introduce the convergence result of SGCN+. Recall that the node embedding approximation variance (zeroth-order variance) determines the bias of stochastic gradient E[‖b‖2F]. 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 η = min{1/LF, 1/ √ T} where LF is the smoothness constant. Let ∆n and ∆+b denote the upper bound for the variance and bias of stochastic gradient as:
  ∆n = L∑ `=1 O(E[‖L̃(`) −P(`)‖2 F ]) +O(E[‖P(`) − L‖2F]), ∆+b = η 2∆+′b where ∆ +′ b = O ( K
  L∑ `=1 |E[‖P(`)‖2F]− ‖L‖2F| ) (8)
  Then, the output of SGCN+ satisfies
  min t∈[T ]
  E[‖∇L(θt)‖2F] ≤ 2(L(θ1)− L(θ?))√ T + LF∆n√ T + ∆+′b T . (9)
  The exact value of key parameters LF, ∆n, and ∆+b are computed in Lemma 1, Lemma 2, and Lemma 5 respectively, and can be found in Appendix G, H. Theorem 2 implies that after T iterations the gradient norm of SGCN+ is at most O(∆n/ √ T ) + O(∆+′b /T ). When using all neighbors for calculating the exact node embeddings, we have P(`) = L such that ∆+′b = 0, which leads to convergence rate of SGD. Comparing with vallina SGCN, the bias of SGCN+ is scaled by learning rate η. Therefore, we can reduce the negative effect of bias by choose learning rate as η = O(1/ √ T ). This also explains why SGCN+ achieves a significantly better convergence rate compared to SGCN.
  5 SGCN++: DOUBLY VARIANCE REDUCTION
  Algorithm 1 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 variance ∆n, 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-order variance reduction on layerwise gradient? To answer this question, we propose doubly variance reduction algorithm SGCN++, that extends the variance reduction algorithm 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 (t mod K = 0) backward propagation are full-batch gradient computation as 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 6= 0) at the regular steps are defined as
  G̃ (`) t = G̃ (`) t−1 + [L̃ (`)H̃ (`−1) t ]
  > ( D̃
  (`+1) t ◦ σ′(Z̃t) ) − [L̃(`)H̃(`−1)t−1 ]> ( D̃ (`+1) t−1 ◦ σ′(Z̃t−1) ) D̃
  (`) t = D̃ (`) t−1 + [L̃
  (`)]> ( D̃
  (`+1) t ◦ σ′(Z̃t)
  ) [W
  (`) t ]− [L̃(`)]>
  ( D̃
  (`+1) t−1 ◦ σ′(Z̃t−1)
  ) [W
  (`) t−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[‖b‖2F] that is due to node embedding approximation and variance E[‖n‖2F] that is due to layerwise gradient estimation. Applying doubly variance reduction on node embedding and layerwise gradient simultaneously can significantly reduce mean-square error of stochastic gradient and speed up convergence. Theorem 3 (Convergence of SGCN++). Suppose Assumptions 1, 2, 3 hold, and denote LF as the smoothness constant and ∆++n+b as the upper-bound of mean-square error of stochastic gradient
  ∆++n+b = η 2∆++′n+b = η 2O ( K ∑L `=1 |E[‖L̃(`)‖2F]− ‖L‖2F| ) (11)
  Apply SGCN++ in Algorithm 2 with learning rate as η = 2 LF+ √ L2F+4∆ ++′ n+b . Then it holds that
  1 T T∑ t=1 E[‖∇L(θt)‖2] ≤ 1 T ( LF + √ L2F + 4∆ ++′ n+b )( L(θ1)− L(θ?) ) . (12)
  The exact value of key parameter LF and ∆++′n+b are computed in Lemma 1 and Lemma 12 respectively, and can be found in Appendix G, I. Theorem 2 implies that applying doubly variance reduction can scale the mean-square error O(η2K) times smaller. As a result, after T iterations the norm of gradient of solution obtained by SGCN++ is at most O(∆++′n+b/T ), which enjoys the same rate as vanilla variance reduced SGD (Reddi et al., 2016; Fang et al., 2018).
  Scalability of SGCN++. One might doubt whether the computation at the snapshot step with fullbatch 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), where the probability of the approximation error violating the desired accuracy decreases exponentially as the number of samples increase. Please refer to 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 C. We remark that large-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 as a special case of the composite optimization problem, it is worth noting that compared to the classical composite optimization, there are a few key differences that make the utilization of variance reduction methods for composite optimization non-trivial: (a) different objective function that makes the GCN analysis challenging; (b) different gradient computation, analysis, and algorithm which make
  directly applying multi-level variance reduction methods such as SPIDER (Zhang & Xiao, 2019b) nontrivial; (c) different theoretical results and novel intuition for sampling-based GCN training. Due to the space limit, we defer the detail discussion to the Appendix A.
  6 EXPERIMENTS
  Experimental setup. We evaluate our proposed methods in semi-supervised learning setting on various classification datasets, summarized in Appendix B. In addition to different sampling mechanisms, we introduce Exact sampling that takes all neighbors for node embedding computation during mini-batch training (no zeroth-order variance), which can be used to explicitly test the importance of first-order variance reduction. We add SGCN+(Zeroth) and SGCN++(Doubly) on top of each sampling method to 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 computation on Flickr, PPI, PPI-large datasets. To scale variance reduced algorithm on large-graph, we employ snapshot large-batch approximation for both SGCN+ and SGCN++ by randomly selecting 50% of nodes for Reddit and 15% of nodes for Yelp dataset. We update the model with a mini-batch size of B = 512 and Adam optimizer with a learning rate of η = 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 3 in Appendix B. Due to the space limit, more experiments can be found in Appendix C.
  Experiment results. In Table 1 and Figure 2, we show the accuracy and convergence comparison of SGCN, SGCN+, and SGCN++. We remark that multi-class classification tasks prefer a more stable 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, and Yelp. Applying variance reductions can further improve its performance. In addition, we observe that the effect of variance reduction depends on its base sampling algorithms. Even though the performance of base sampling algorithm various significantly, the doubly variance reduction can bring their performance to a similar level. Moreover, we can observe from the loss curves that 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 variance than GraphSAGE because of its zeroth-order variance reduction, and GraphSAGE has less variance than LADIES because more nodes are sampled for node embedding approximation.
  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. Furthermore, 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 Lé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 principled 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 arXiv:1710.10568, 2017.
  Jie Chen and Ronny Luss. Stochastic gradient descent with biased but consistent gradient estimators. arXiv preprint arXiv:1807.11880, 2018.
  Jie Chen, Tengfei Ma, and Cao Xiao. Fastgcn: fast learning with graph convolutional networks via importance sampling. arXiv preprint arXiv:1801.10247, 2018.
  Minshuo Chen, Xingguo Li, and Tuo Zhao. On generalization bounds of a family of recurrent neural networks. arXiv preprint arXiv:1910.12947, 2019.
  Tianyi Chen, Yuejiao Sun, and Wotao Yin. Solving stochastic compositional optimization is nearly as easy as solving stochastic optimization. arXiv preprint arXiv:2008.10847, 2020.
  Wei-Lin Chiang, Xuanqing Liu, Si Si, Yang Li, Samy Bengio, and Cho-Jui Hsieh. Cluster-gcn: 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.00982, 2020.
  Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. Spider: Near-optimal non-convex optimization 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.06157, 2020.
  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 in 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 networks. arXiv preprint arXiv:1609.02907, 2016.
  Sashank J Reddi, Ahmed Hefny, Suvrit Sra, Barnabás Póczos, and Alex Smola. Stochastic variance reduction for nonconvex optimization. In International conference on machine learning, pp. 314– 323, 2016.
  Ryoma Sato, Makoto Yamada, and Hisashi Kashima. Constant time graph neural networks, 2020. URL https://openreview.net/forum?id=rkgKW64FPH.
  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.
  Quoc Tran-Dinh, Deyi Liu, and Lam M Nguyen. Hybrid variance-reduced sgd algorithms for nonconvex-concave minimax problems. arXiv preprint arXiv:2006.15266, 2020.
  Quan Wang, Zhendong Mao, Bin Wang, and Li 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 arXiv: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):616–659, 2019.
  Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec. Graph convolutional neural networks for web-scale recommender systems. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 974– 983. ACM, 2018.
  Hanqing Zeng, Hongkuan Zhou, Ajitesh Srivastava, Rajgopal Kannan, and Viktor Prasanna. Graphsaint: Graph sampling based inductive learning method. arXiv preprint arXiv:1907.04931, 2019.
  Junyu Zhang and Lin Xiao. A composite randomized incremental gradient method. In International Conference on Machine Learning, pp. 7454–7462, 2019a.
  Junyu Zhang and Lin Xiao. Multi-level composite stochastic optimization via nested variance reduction. arXiv preprint arXiv:1908.11468, 2019b.
  Junyu Zhang and Lin Xiao. A stochastic composite gradient method with incremental variance reduction. In Advances in Neural Information Processing Systems, pp. 9078–9088, 2019c.
  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 arXiv:2001.11850, 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. 11249–11259, 2019.
  Appendix
  Table of Contents A Connection to composite optimization 12
  B Experiment configurations 13
  C Additional experiments 15
  D Different sampling strategies 19
  E Detailed algorithms 21 E.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 E.2 SGCN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 E.3 SGCN+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 E.4 SGCN++ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 E.5 SGCN++ without full-batch . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
  F Notations, important propositions and lemmas 24 F.1 Notations for gradient computation . . . . . . . . . . . . . . . . . . . . . . . . 24 F.2 Upper-bounded on the node embedding matrices and layerwise gradients . . . . 25 F.3 Lipschitz continuity and smoothness property of graph convolution layers . . . . 26 F.4 Lipschitz continouity of the gradient of graph convolutional network . . . . . . . 27
  G Proof of Theorem 1 30 G.1 Supporting lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 G.2 Remaining steps toward Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . 34
  H Proof of Theorem 2 36 H.1 Supporting lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 H.2 Remaining steps toward Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . 40
  I Proof of Theorem 3 41 I.1 Supporting lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 I.2 Remaining steps toward Theorem 3 . . . . . . . . . . . . . . . . . . . . . . . . 52
  J Reproducing experiment results 53
  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 optimization counterparts (e.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 in 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 as
  F (θ) = 1
  N N∑ i=1 fi ( 1 M M∑ j=1 gj(w) ) , θ = {w}, (13)
  where fi(·) is the outer-level function computed on the ith data point, gj(·) is the inner-level function computed on the jth data point, and w is the parameter. We denote ∇fi(·) and ∇gj(·) as the gradient. Then, the gradient for Eq. 13 is computed as
  ∇F (θ) = [ 1 N N∑ i=1 ∇fi ( 1 M M∑ j=1 gj(w) )]( 1 M M∑ j=1 ∇gj(w) ) , θ = {w}, (14)
  where the dependency between inner- and outer-level sampling are not considered. One can independently sample inner layer data to estimate g̃ ≈ 1M ∑M j=1 gj(w) and ∇g̃ ≈ 1 M ∑M j=1∇gj(w), sam-
  ple outer layer data to estimate∇f̃ ≈ 1N ∑N i=1∇fi(g̃), then estimate∇F (θ) by using [∇f̃ ]>∇g̃.
  By casting the optimizaion problem in GCN as composite optimization problem in Eq. 13, we have
  L(θ) = 1 B ∑ i∈VB Loss(h(L)i , yi), θ = {W (1)}
  H(L) = σ(L̃(L)XW̃(L)),W̃(L) = σ ( L̃(L−1)Xσ ( L̃(L−2)X . . . σ ( L̃(1)XW(1) ) ︸ ︷︷ ︸
  W̃(2)
  . . . ) ,
  (15)
  which is different from the vanilla GCN model. To see this, we note that in vanilla GCNs, since the sampled nodes at the `th layer are dependent from the nodes sampled at the (` + 1)th layer, we have E[L̃(`)] = P(`) 6= 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 E[L̃(`)] = 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 to update the parameters in Eq. 15 are computed as
  ∂L(θ) ∂W̃(`) = ∂L(θ) ∂W̃(L) ( L∏ j=`+1 ∂W̃(j) ∂W̃(j−1) ) . (16)
  However in GCN, there are two types of gradient at each layer (i.e., D̃(`) and G̃(`)) that are fused with each other (i.e., D̃(`) is a part of G̃(`−1) and D̃(`) is a part of D̃(`−1)) but with different functionality. D̃(`) is passing gradient between different layers, G̃(`) is passing gradient to 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 W̃(`)t in Eq. 13 (please refer to Algorithm 3 in Zhang & Xiao (2019b)),
  where W̃(`)t−1 is used as a control variant to reduce the variance of W̃ (`) t , i.e.,
  W̃ (`+1) t = W̃ (`+1) t−1 + σ(L̃ (`) t XW̃ (`) t )− σ(L̃ (`) t XW̃ (`) t−1). (17)
  However in SGCN++, the zeroth-order variance reduction is applied to H̃(`)t . Because the node sampled at the tth and (t−1)th iteration are unlikely the same, we cannot directly use H(`)t−1 to reduce the variance of H(`)t . Instead, the control variant in SGCN++ is computed by applying historical weight W(`)t−1 on the historical node embedding from previous layer H (`−1) t−1 , i.e.,
  H̃ (`) t = H̃ (`) t−1 + σ(L̃ (`) t H (`−1) t W (`) t )− σ(L̃ (`) t H (`−1) t−1 W (`) t−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 a novel analysis of Theorem 1, where we show that the vanilla sampling-based GCNs suffer a residual error ∆b 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 on the convergence of GCNs Chen 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 3 are also different.
  B EXPERIMENT CONFIGURATIONS
  Hardware specification and environment. We run our experiments on a single machine with Intel i5-7500, NVIDIA GTX1080 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 to 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 modified the PyTorch implementation of GCN (Kipf & Welling, 2016)3 to add LADIES (Zou et al., 2019), FastGCN (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 SGCN++ 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 L = D−1/2AD−1/2. We use mean-aggregation for single-class classification task and concatenate-aggregation for multi-class classification. The default mini-batch batch size and sampled node size are summarized in Table 2. We update the model using Adam optimizer with a learning rate of 0.01. For SGCN++, historical node embeddings are first calculated on GPUs and transfer to CPU memory using PyTorch command Tensor.to(device). Therefore, no extra GPU memory is required when training with SGCN++. To balance the staleness of snapshot model and the computational efficiency, as default we choose snapshot gapK = 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.
  During training, for each epoch we construct 10 mini-batches in parallel using Python package multiprocessing 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 scipy.sparse 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-connection”, and “attention”), 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 to better explore the impact of sampling and variance reduction on convergence.
  Comparison of SGD and Adam. It is worth noting that Adam optimizer is used as the default optimizer 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 reduction 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 comparison meaningless. For example in Figure 3, we compare Adam and SGD optimizer on PPI dataset. For Adam optimizer we use PyTorch’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 a giant performance gap when comparing to Adam optimizer. (b) Most public implementation of GCNs, including all implementations in PyTorch 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.
  Dataset statistics. We summarize the dataset statistics in Table 3.
  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 dataset.
  Because Exact sampling is using all neighbors for the node embedding approximation, it is only affected by layerwise gradient variance (first-order variance). Therefore, employing first-order variance reduction on Exact sampling can significantly reduce the mean-square error of stochastic gradient to full-gradient, and speed up the convergence.
  Different from Exact sampling that the exact node embeddings are available during training, layerwise sampling algorithm LADIES, and nodewise sampling GraphSAGE are 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 is necessary to reduce the mean-square error of the stochastic gradient and speed up the convergence.
  GPU memory usage. In Figure 5, We compare the GPU memory usage of SGCN and SGCN++. We calculate the allocated memory by torch.cuda.memory allocated, which is the current GPU memory occupied by tensors in bytes for a given device. We calculate the maximum allocated memory by torch.cuda.max memory allocated, which is the maximum GPU memory occupied by tensors in bytes for a given device.
  From Figure 5, we observe that neither running full-batch GCN nor saving historical node embeddings and gradients will significantly increase the computation overhead during training. Besides, since all historical activations are stored outside GPU, we see that SGCN++ 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.
  Evaluation of total time. In Table 4 and Table 5, we report the average time of doubly variance reduced LADIES++ and vanilla LADIES. 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. In practice, we directly use full-batch training for the smaller datasets (e.g., PPI, PPI-large, and Flickr) and use sampled snapshot large-batch for large datasets (e.g., Reddit and Yelp). When constructing snapshot large-batch, the Exact sampler has to go through all neighbors of each node using for-loops based on the graph structure, such that it is 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 to construct the mini-batches using layerwise LADIES sampler. • Regular step transfer time: The time required to 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. • Computation time: The time used for forward- and backward-propagation.
  Notice that we are reporting the total time per iteration because the vanilla sampling-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 samplingbased GCN training is 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, a promising future direction for large-scale graph training is developing a 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, we evaluate the effect of snapshot grap on training by choosing mini-batch size as B = 512 and change the innerloop intervals fromK = 5 mini-batches toK = 20 mini-batches. In Figure 6 and Figure 7, we show
  the comparison of training loss and validation loss with different number of inner-loop intervals for SGCN++ and SGCN+ on Reddit dataset, respectively. We can observe that the model with a smaller snapshot gap requires less number of iterations to reach the same training and validation loss, and gives us a better generalization performance (F1-score).
  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 a 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 on the storage. After that, mini-batch nodes are sampled from the large-batch during the regular steps. Let denote B′ as the snapshot step largebatch size and B denote the regular step mini-batch size. By default, we choose snapshot gap as
  K = 10, fix the regular step batch size as B = 512, and change the snapshot step batch size B′ from 20, 000 (20K) to 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.
  The effect of mini-batch size. In Figure 10, we 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 as 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-batch size every iterations. In Figure 11, we explicitly compare the effectiveness of mini-batch size on doubly variance reduced GraphSAINT++ and vanilla GraphSAINT, 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-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
  a speed boost. In Figure 12, we show the comparison of validation loss of fixed snapshot gapK = 10 and gradually increasing snapshot gap K = 10 + 0.1 × s, s = 1, 2, . . ., where s 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 a performance drop.
  D DIFFERENT SAMPLING STRATEGIES
  In this section, we 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 to 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 s of its neighbors at the (` − 1)th GCN layer and formulate L̃(`) by
  L̃ (`) i,j =
  { |N (i)| s × Li,j , if j ∈ Ñ (`)(i)
  0, otherwise (19)
  where N (i) is the full set of ith node neighbor, Ñ (`)(i) is the sampled neighbors of node i for `th GCN layer. GraphSAGE Hamilton et al. (2017) follows the spirit of node-wise sampling where it performs uniform node sampling on the previous layer neighbors for a fixed number of nodes to bound the mini-batch computation complexity.
  Layer-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 a set of nodes B(`) ⊆ V of size s under the distribution p to approximate the Laplacian by
  L̃ (`) i,j =
  { 1
  s×pj × Li,j , if j ∈ B (`) 0, otherwise (20)
  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 FastGCN by layer-dependent sampling. Based on the sampled nodes in the upper layer, it selects their neighborhood nodes, constructs a bipartite subgraph, and computes the importance probability accordingly. 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 sampled Laplacian matrices at each layer are identical
  L̃ (1) i,j = . . . = L̃ (L) i,j =
  { 1
  s×pj × Li,j , if j ∈ B (`) 0, otherwise (21)
  For example, GraphSAINT Zeng et al. (2019) can be viewed as a special case of layer-wise sampling algorithm FastGCN by restricting the nodes sampled at the 1-st to (L−1)th layer the same as the nodes sampled at the Lth layer. However, GraphSAINT requires a significant large mini-batch size compared to other layer-wise sampling methods. We leave this as a potential future direction to explore.
  E DETAILED ALGORITHMS
  E.1 DESCRIPTION
  In order to help readers better compare the difference of different algorithms, we summarize the vanilla sampling-based GCN training algorithm SGCN in Algorithm 3, zeroth-order variance reduced algorithm SGCN+ in Algorithm 4, and doubly variance reduced algorithm SGCN++ in Algorithm 5. In addition, we illustrate in Figure 13 the relationship between node embedding approximation variance and layerwise gradient variance to the forward- and backward-propagation. We remark that zeroth-order variance reduction SGCN+ is only applied during forward-propagation, and doubly (zeroth- and first-order) variance reduction SGCN++ is applied during forward- and backward-propagation, simultaneously.
  E.2 SGCN
  Algorithm 3 SGCN: Vanilla sampling-based GCN training method 1: Input: Learning rate η > 0 2: for t = 1, . . . , T do 3: Sample mini-batch VB ⊂ V 4: Calculate node embeddings using
  H̃(`) = σ(L̃(`)H̃(`−1)W(`)), where H̃(0) = X, (22)
  5: Calculate loss as L̃(θt) = 1B ∑ i∈VB Loss(h̃ (L) i , yi) 6: Calculate stochastic gradient∇L̃(θt) = {G̃(`)}L`=1 as
  G̃ (`) t := [L̃ (`)H̃ (`−1) t ] > ( D̃ (`+1) t ◦ ∇σ(Z̃ (`) t ) ) ,
  D̃ (`) t := [L̃
  (`)]> ( D̃
  (`+1) t ◦ ∇σ(Z̃ (`) t ) ) W (`) t , D̃ (L+1) t =
  ∂L̃(θt) ∂H̃(L)
  (23)
  7: Update parameters as θt+1 = θt − η∇L̃(θt) 8: end for 9: Output: Model with parameter θT+1
  E.3 SGCN+
  Algorithm 4 SGCN+: Zeroth-order variance reduction (Detailed version of Algorithm 4) 1: Input: Learning rate η > 0, snapshot gap K > 0 2: for t = 1, . . . , T do 3: if t mod K = 0 then 4: % Snapshot steps
  5: Calculate node embeddings and update historical node embeddings using
  Z (`) t = LH (`−1) t W (`) t , H (`) t = σ(Z (`) t ), Z̃ (`) t ← Z (`) t (24)
  6: Calculate loss as L(θt) = 1N ∑N i=1 Loss(h (L) i , yi) 7: Calculate full-batch gradient∇L(θt) = {G(`)}L`=1 as
  G (`) t := [LH (`−1) t ] > ( D (`+1) t ◦ ∇σ(Z (`) t ) ) ,
  D (`) t := L
  > ( D
  (`+1) t ◦ ∇σ(Z (`) t ) ) W (`) t , D (L+1) t =
  ∂L(θt) ∂H(L)
  (25)
  8: Update parameters as θt+1 = θt − η∇L(θt) 9: else
  10: % Regular steps 11: Sample mini-batch VB ⊂ V 12: Calculate node embeddings using
  Z̃ (`) t = Z̃ (`) t−1 + L̃ (`)H̃ (`−1) t W (`) t − L̃(`)H̃ (`−1) t−1 W (`) t−1, H̃ (`) t = σ(Z̃ (`)) (26)
  13: Calculate loss as L̃(θt) = 1B ∑ i∈VB Loss(h̃ (L) i , yi) 14: Calculate the stochastic gradient∇L̃(θt) = {G̃(`)}L`=1 as
  G̃ (`) t := [L̃ (`)H̃ (`−1) t ] > ( D̃ (`+1) t ◦ ∇σ(Z̃ (`) t ) ) ,
  D̃ (`) t := [L̃
  (`)]> ( D̃
  (`+1) t ◦ ∇σ(Z̃ (`) t ) ) W (`) t , D̃ (L+1) t =
  ∂L̃(θt) ∂H̃(L)
  (27)
  15: Update parameters as θt+1 = θt − η∇L̃(θt) 16: end if 17: end for 18: Output: Model with parameter θT+1
  E.4 SGCN++
  Algorithm 5 SGCN++: Doubly variance reduction (Detailed version of Algorithm 5) 1: Input: Learning rate η > 0, snapshot gap K > 0 2: for t = 1, . . . , T do 3: if t mod K = 0 then 4: % Snapshot steps 5: Calculate node embeddings and update historical node embeddings using
  Z (`) t = LH (`−1) t W (`) t , H (`) t = σ(Z (`) t ), Z̃ (`) t ← Z (`) t (28)
  6: Calculate loss as L(θt) = 1N ∑N i=1 Loss(h (L) i , yi) 7: Calculate the full-batch gradient∇L(θt) = {G(`)}L`=1 as
  G (`) t := [LH (`−1) t ] > ( D (`) t ◦ ∇σ(Z (`) t ) ) ,
  D (`) t := L
  > ( D
  (`+1) t ◦ ∇σ(Z (`) t ) ) W (`) t , D (L+1) t =
  ∂L(θt) ∂H(L)
  (29)
  8: Save the per layerwise gradient G̃(`)t ← G (`) t , D̃ (`) t ← D (`) t for all ` ∈ [L]
  9: Update parameters as θt+1 = θt − η∇L(θt) 10: else 11: % Regular steps 12: Sample mini-batch VB ⊂ V
  13: Calculate node embeddings using
  Z̃ (`) t = Z̃ (`) t−1 + L̃ (`)H̃ (`−1) t W (`) t − L̃(`)H̃ (`−1) t−1 W (`) t−1, H̃ (`) t = σ(Z̃ (`)) (30)
  14: Calculate loss as L̃(θt) = 1B ∑ i∈VB Loss(h̃ (L) i , yi) 15: Calculate the stochastic gradient∇L̃(θt) = {G̃(`)}L`=1 as
  G̃ (`) t = G̃ (`) t−1 + [L̃ (`)H̃ (`−1) t ]
  > ( D̃
  (`+1) t ◦ ∇σ(Z̃t) ) − [L̃(`)H̃(`−1)t−1 ]> ( D̃ (`+1) t−1 ◦ ∇σ(Z̃t−1)
  ) D̃
  (`) t = D̃ (`) t−1 + [L̃
  (`)]> ( D̃
  (`+1) t ◦ ∇σ(Z̃t)
  ) [W
  (`) t ]
  − [L̃(`)]> ( D̃
  (`+1) t−1 ◦ ∇σ(Z̃t−1)
  ) [W
  (`) t−1], D̃ (L+1) t = ∂L̃(θt) ∂H̃
  (L) t
  (31)
  16: Update parameters as θt+1 = θt − η∇L̃(θt) 17: end if 18: end for 19: Output: Model with parameter θT+1
  E.5 SGCN++ WITHOUT FULL-BATCH
  Furthermore, in Algorithm 6, we provide an alternative version of SGCN++ that does not require fullbatch forward- and backward-propagation at the snapshot step. The basic idea is to approximate the full-batch gradient by sampling a large mini-batch V ′B of size B′ = |V ′B| using Exact sampling, then compute the node embedding matrices and stochastic gradients on the sampled large-batch V ′B.
  Algorithm 6 SGCN++ (without full-batch): Doubly variance reduction 1: Input: Learning rate η > 0, snapshot gap K > 0 2: for t = 1, . . . , T do 3: if t mod K = 0 then 4: % Snapshot steps 5: Sample a large-batch V ′B of sizeB′ and construct the Laplacian matrices L(`) for each layer
  using all neighbors, i.e.,
  L (`) i,j = { Li,j , if j ∈ N (`)(i) 0, otherwise
  (32)
  6: Calculate node embeddings and update historical node embeddings using
  Z (`) t = L (`)H (`−1) t W (`) t , H (`) t = σ(Z (`) t ), Z̃ (`) t ← Z (`) t (33)
  7: Calculate loss as L(θt) = 1B′ ∑ i∈V′B Loss(h(L)i , yi) 8: Calculate the approximated snapshot gradient∇L(θt) = {G(`)}L`=1 as
  G (`) t := [L (`)H (`−1) t ] > ( D (`+1) t ◦ ∇σ(Z (`) t ) ) ,
  D (`) t := [L
  (`)]> ( D
  (`+1) t ◦ ∇σ(Z (`) t ) ) W (`) t , D (L+1) t =
  ∂L(θt) ∂H(L)
  (34)
  9: Save the per layerwise gradient G̃(`)t ← G (`) t , D̃ (`) t ← D (`) t , ∀` ∈ [L]
  10: Update parameters as θt+1 = θt − η∇L(θt) 11: else 12: % Regular steps 13: Sample mini-batch VB ⊂ V ′B 14: Calculate node embeddings using
  Z̃ (`) t = Z̃ (`) t−1 + L̃ (`)H̃ (`−1) t W (`) t − L̃(`)H̃ (`−1) t−1 W (`) t−1, H̃ (`) t = σ(Z̃ (`)) (35)
  15: Calculate loss as L̃(θt) = 1B ∑ i∈VB Loss(h̃ (L) i , yi) 16: Calculate the stochastic gradient∇L̃(θt) = {G̃(`)}L`=1 as
  G̃ (`) t = G̃ (`) t−1 + [L̃ (`)H̃ (`−1) t ]
  > ( D̃
  (`+1) t ◦ ∇σ(Z̃t) ) − [L̃(`)H̃(`−1)t−1 ]> ( D̃ (`+1) t−1 ◦ ∇σ(Z̃t−1)
  ) D̃
  (`) t = D̃ (`) t−1 + [L̃
  (`)]> ( D̃
  (`+1) t ◦ ∇σ(Z̃t)
  ) [W
  (`) t ] − [L̃(`)]> ( D̃
  (`+1) t−1 ◦ ∇σ(Z̃t−1)
  ) [W
  (`) t−1],
  D̃ (L+1) t = ∂L̃(θt) ∂H̃
  (L) t
  (36)
  17: Update parameters as θt+1 = θt − η∇L̃(θt) 18: end if 19: end for 20: Output: Model with parameter θT+1
  The intuition of snapshot step large-batch approximation stems from matrix Bernstein inequality Gross (2011). More specifically, suppose given G̃i ∈ Rd×d be the stochastic gradient computed by using the ith node with Exact sampling (all neighbors are used to calculate the exact node embeddings). Suppose the different between G̃i and full-gradient E[G̃i] is uniformly bounded and the variance is bounded:
  ‖G̃i − E[G̃i]‖F ≤ µ, E[‖G̃i − E[G̃i]‖2F ] ≤ σ2 (37)
  Let G̃′ as the snapshot step gradient computed on the sampled large batch
  G̃′ = 1
  B′ ∑ i∈V′B G̃i (38)
  By matrix Bernstein inequality, we know the probability of ‖G̃′−E[G̃i]‖F larger than some constant decreases exponentially as the size of the sampled large-batch size B′ increase, i.e.,
  Pr(‖G̃′ − E[G̃i]‖F ≥ ) ≤ 2d exp ( −n ·min { 2
  4σ2 , 2µ
  }) (39)
  Therefore, by choosing a large enough snapshot step batch size B′, we can obtain a good approximation 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(`) = [f (`)(H(`−1),W(`)) := σ(LH(`−1)W(`)︸ ︷︷ ︸ Z(`) )] ∈ RN×d` (40)
  and its gradient w.r.t. the input node embedding matrix D(`) ∈ RN×d`−1 is computed as D(`) = [ ∇Hf (`)(D(`+1),H(`−1),W(`)) := [L]> ( D(`+1) ◦ σ′(LH(`−1)W(`)) ) [W(`)]> ] (41)
  and its gradient w.r.t. the weight matrix G(`) ∈ Rd`−1×d` is computed as G(`) = [ ∇W f (`)(D(`+1),H(`−1),W(`)) := [LH(`−1)]> ( D(`+1) ◦σ′(L(`)H(`−1)W(`)) )] (42)
  Similarly, we can formulate the calculation of node embedding matrix H̃(`) ∈ RN×d` at each GCN layer in SGCN as
  H̃(`) = [f̃ (`)(H̃(`−1),W(`)) := σ(L̃(`)H̃(`−1)W(`)︸ ︷︷ ︸ Z̃(`) )] (43)
  and its gradient w.r.t. the input node embedding matrix D̃(`) ∈ RN×d`−1 is computed as D̃(`) = [ ∇H f̃ (`)(D̃(`+1), H̃(`−1),W(`)) := [L̃(`)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)) ) [W(`)]> ] (44) and its gradient w.r.t. the weight matrix G̃(`) ∈ Rd`−1×d` is computed as
  G̃(`) = [ ∇W f̃ (`)(D̃(`+1), H̃(`−1),W(`)) := [L̃(`)H̃(`−1)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)) )] (45)
  Let us denote the gradient of loss w.r.t. the final node embedding matrix as
  D(L+1) = ∂Loss(H(L),y)
  ∂H(L) ∈ RN×dL , [D(L+1)]i =
  1
  N
  ∂Loss(h(L)i , yi)
  ∂h (L) i
  ∈ RdL
  D̃(L+1) = ∂Loss(H̃(L),y)
  ∂H̃(L) ∈ RN×dL , [D̃(L+1)]i =
  1 B 1{i∈VB}
  ∂Loss(h̃(L)i , yi)
  ∂h̃ (L) i
  ∈ RdL (46)
  Notice that D̃(L+1) is a N × dL matrix with the row number i ∈ VB are non-zero vectors. Then we can write the gradient for the `th weight matrix in FullGCN and SGCN as
  G(`) = ∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1),H(L−1),W(L)) . . . ,H(`),W(`+1)),H(`−1),W(`))
  G̃(`) = ∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̃(L+1), H̃(L−1),W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`)) (47)
  F.2 UPPER-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 to the (` − 1)th layer node embedding matrix.
  Proposition 2 (Detailed version of Proposition 1). For any ` ∈ [L], the Frobenius norm of node embedding matrices, gradient passing from the `th layer node embeddings to the (` − 1)th are bounded
  ‖H(`)‖F ≤ BH , ‖H̃(`)‖F ≤ BH , ‖ ∂σ(LH(`−1)W(`))
  ∂H(`−1) ‖F ≤ BD, ‖
  ∂σ(L̃(`)H̃(`−1)W(`))
  ∂H̃(`−1) ‖F ≤ BD
  (48) where
  BH = ( CσBLABW )` BH , BD = ( BLACσBW )L−` Closs (49)
  Proof.
  ‖H(`)‖F = ‖σ(LH(`−1)W(`))‖F ≤ CσBLABW ‖H(`−1)‖F
  ≤ ( CσBLABW )` ‖X‖F ≤ ( CσBLABW )` BH
  (50)
  ‖H̃(`)‖F = ‖σ(L̃(`)H̃(`−1)W(`))‖F ≤ CσBLABW ‖H̃(`−1)‖F
  ≤ ( CσBLABW )` ‖X‖F ≤ ( CσBLABW )` BH
  (51)
  ‖D(`)‖F = ∥∥∥[L]>(D(`+1) ◦ σ′(LH(`−1)W(`)))[W(`)]>∥∥∥
  F
  ≤ BLACσBW ‖D(`+1)‖F ≤ ( BLACσBW )L−` ‖D(L+1)‖F ≤ ( BLACσBW )L−` Closs
  (52)
  ‖D̃(`)‖F = ∥∥∥[L̃(`)]>(D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)))[W(`)]>∥∥∥
  F
  ≤ BLACσBW ‖D̃(`+1)‖F ≤ ( BLACσBW )L−` ‖D̃(L+1)‖F ≤ ( BLACσBW )L−` Closs
  (53)
  F.3 LIPSCHITZ CONTINUITY AND SMOOTHNESS PROPERTY OF GRAPH CONVOLUTION LAYERS
  Then, we derive the Lipschitz continuity of f̃ (`)(·, ·) and its gradient. Notice that the following result hold for deterministic function f (`)(·, ·) as well.
  Proposition 3. f̃ (`)(·, ·) is CH -Lipschitz continuous w.r.t. the input node embedding matrix where CH = CσBLABW
  Proof.
  ‖f̃ (`)(H(`−1)1 ,W(`))− f̃ (`)(H (`−1) 2 ,W (`))‖F = ‖σ(L̃(`)H(`−1)1 W(`))− σ(L̃(`)H (`−1) 2 W
  (`))‖F ≤ Cσ‖L̃(`)‖F‖H(`−1)1 −H (`−1) 2 ‖F‖W(`)‖F ≤ CσBLABW ‖H(`−1)1 −H (`−1) 2 ‖F
  (54)
  Proposition 4. f̃ (`)(·, ·) is CW -Lipschitz continuous w.r.t. the weight matrix where CW = CσBLABH
  Proof.
  ‖f̃ (`)(H(`−1),W(`)1 )− f̃ (`)(H(`−1),W (`) 2 )‖F
  = ‖σ(L̃(`)H(`−1)W(`)1 )− σ(L̃(`)H(`−1)W (`) 2 )‖F
  ≤ Cσ‖L̃(`)‖F‖H(`−1)‖F‖W(`)1 −W (`) 2 ‖F
  ≤ CσBLABH‖W(`)1 −W (`) 2 ‖F
  (55)
  Proposition 5. ∇H f̃ (`)(·, ·, ·) is LH -Lipschitz continuous where LH = max { BLACσBW , B 2 LABDB 2 WLσ, BLABDCσ +B 2 LABDBWLσBH } (56)
  Proof.
  ‖∇H f̃ (`)(D̃(`+1)1 , H̃(`−1),W(`))−∇H f̃ (`)(D̃ (`+1) 2 , H̃ (`−1),W(`))‖F ≤ ‖[L̃(`)]> ( D̃
  (`+1) 1 ◦ σ′(L̃(`)H̃(`−1)W(`))
  ) [W(`)]>
  − [L̃(`)]> ( D̃
  (`+1) 2 ◦ σ′(L̃(`)H̃(`−1)W(`)) ) [W(`)]>‖F
  ≤ BLACσBW ‖D̃(`+1)1 − D̃ (`+1) 2 ‖F
  (57)
  ‖∇H f̃ (`)(D̃(`+1), H̃(`−1)1 ,W(`))−∇H f̃ (`)(D̃(`+1), H̃ (`−1) 2 ,W (`))‖F ≤ ‖[L̃(`)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)1 W(`)) ) [W(`)]>
  − [L̃(`)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)2 W(`)) ) [W(`)]>‖F
  ≤ B2LABDB2WLσ‖H̃ (`−1) 1 − H̃ (`−1) 2 ‖F
  (58)
  ‖∇H f̃ (`)(D̃(`+1), H̃(`−1),W(`)1 )−∇H f̃ (`)(D̃(`+1), H̃(`−1),W (`) 2 )‖F ≤ ‖[L̃(`)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)1 ) ) [W (`) 1 ] >
  − [L̃(`)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)2 ) ) [W (`) 2 ] >‖F
  ≤ (BLABDCσ +B2LABDBWLσBH)‖W (`) 1 −W (`) 2 ‖F
  (59)
  Proposition 6. ∇W f̃ (`)(·, ·, ·) is LW -Lipschitz continuous where LW = max { BLABHCσ, B 2 LAB 2 HHDLσ, BLABDCσ +B 2 LAB 2 HBDLσ } (60)
  Proof. ‖∇W f̃ (`)(D̃(`+1)1 , H̃(`−1),W(`))−∇W f̃ (`)(D̃ (`+1) 2 , H̃ (`−1),W(`))‖F
  ≤ ‖[L̃(`)H̃(`−1)]> ( D̃
  (`+1) 1 ◦ σ′(L̃(`)H̃(`−1)W(`)) ) − [L̃(`)H̃(`−1)]> ( D̃
  (`+1) 2 ◦ σ′(L̃(`)H̃(`−1)W(`)) ) ‖F
  ≤ BLABHCσ‖D̃(`+1)1 − D̃ (`+1) 2 ‖F
  (61)
  ‖∇W f̃ (`)(D̃(`+1), H̃(`−1)1 ,W(`))−∇W f̃ (`)(D̃(`+1), H̃ (`−1) 2 ,W (`))‖F ≤ ‖[L̃(`)H̃(`−1)1 ]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)1 W(`)) ) − [L̃(`)H̃(`−1)2 ]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)2 W(`)) ) ‖F
  ≤ (B2LAB2DCσ +B2LABHBDLσBW )‖H̃ (`−1) 1 − H̃ (`−1) 2 ‖F
  (62)
  ‖∇W f̃ (`)(D̃(`+1), H̃(`−1),W(`)1 )−∇W f̃ (`)(D̃(`+1), H̃(`−1),W (`) 2 )‖F ≤ ‖[L̃(`)H̃(`−1)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)1 ) ) − [L̃(`)H̃(`−1)]> ( D̃(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)2 ) ) ‖F
  ≤ B2LAB2HHDLσ‖W (`) 1 −W (`) 2 ‖F
  (63)
  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 as
  θ = {W(1), . . . ,W(L)}, ∇L(θ) = {G(1), . . . ,G(L)} (64)
  where G(`) is defined as the gradient w.r.t. the `th layer weight matrix. Let us slight abuse of notation and define the distance between two set of parameters θ1,θ2 and its gradient as
  ‖θ1 − θ2‖F = L∑ `=1 ‖W(`)1 −W (`) 2 ‖F, ‖∇L(θ1)−∇L(θ2)‖F = L∑ `=1 ‖G(`)1 −G (`) 2 ‖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.
  Lemma 1. The gradient of an L-layer GCN is LF-Lipschitz continuous with LF = L(LU2maxLU 2 maxC + U 2 maxL), i.e.,
  ‖∇L(θ1)−∇L(θ2)‖2F ≤ LF‖θ1 − θ2‖2F (66) where
  UmaxC = max j∈{0,...,L−1}
  CjHCW ,
  UmaxL = [
  min j∈{0,...,L−`−1}
  LWL j H
  ] ×max { 1, Lloss } (67) Proof. We first consider the gradient w.r.t. the `th graph convolutional layer weight matrix
  ‖G(`)1 −G (`) 2 ‖F
  = ‖∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)1 ,H (L−1) 1 ,W (L) 1 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 )
  −∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)2 ,H (L−1) 2 ,W (L) 2 ) . . . ,H (`) 2 W (`+1) 2 ),H (`−1) 2 ,W (`) 2 )‖F
  ≤ ‖∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)1 ,H (L−1) 1 ,W (L) 1 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 )
  −∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)2 ,H (L−1) 1 ,W (L) 1 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 )‖F + ‖∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)2 ,H (L−1) 1 ,W (L) 1 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 ) −∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)2 ,H (L−1) 2 ,W (L) 1 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 )‖F + ‖∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)2 ,H (L−1) 2 ,W (L) 1 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 ) −∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1)2 ,H (L−1) 2 ,W (L) 2 ) . . . ,H (`) 1 ,W (`+1) 1 ),H (`−1) 1 ,W (`) 1 )‖F + . . . + ‖∇W f (`)(D(`+1)2 ,H (`−1) 1 ,W (`) 1 )−∇W f (`)(D (`+1) 2 ,H (`−1) 2 ,W (`) 1 )‖F + ‖∇W f (`)(D(`+1)2 ,H (`−1) 2 ,W (`) 1 )−∇W f (`)(D (`+1) 2 ,H (`−1) 2 ,W (`) 2 )‖F
  (68)
  By the Lipschitz continuity of∇W f (`)(·) and ∇Hf (`)(·)
  ‖G(`)1 −G (`) 2 ‖F ≤ LWL L−`−1 H Lloss‖H (L) 1 −H (L) 2 ‖F
  + LWL L−`−1 H (‖H (L−1) 1 −H (L−1) 2 ‖F + ‖W (L) 1 −W (L) 2 ‖F) + . . .
  + LWLH(‖H(`)1 −H (`) 2 ‖F + ‖W (`+1) 1 −W (`+1) 2 ‖F) + LW (‖H(`−1)1 −H (`−1) 2 ‖F + ‖W (`) 1 −W (`) 2 ‖F)
  (69)
  Let define UmaxL as UmaxL = [
  min j∈{0,...,L−`−1}
  LWL j H
  ] ×max { 1, Lloss } (70)
  then we can rewrite the above equation as
  ‖G(`)1 −G (`) 2 ‖F ≤ UmaxL ( L∑ j=1 ‖H(j)1 −H (j) 2 ‖F ) + UmaxL ( L∑ j=1 ‖W(j)1 −W (j) 2 ‖F ) (71)
  Then, let consider the upper bound of ‖H(`)1 −H (`) 2 ‖F
  ‖H(`)1 −H (`) 2 ‖F = ‖f (`)(f (`−1)(. . . f (1)(X,W (1) 1 ) . . . ,W (`−1) 1 ),W (`) 1 )
  − f (`)(f (`−1)(. . . f (1)(X,W(1)2 ) . . . ,W (`−1) 2 ),W (`) 2 )‖F
  ≤ ‖f (`)(f (`−1)(. . . f (1)(X,W(1)1 ) . . . ,W (`−1) 1 ),W (`) 1 )
  − f (`)(f (`−1)(. . . f (1)(X,W(1)2 ) . . . ,W (`−1) 1 ),W (`) 1 )‖F + . . . + ‖f (`)(f (`−1)(. . . f (1)(X,W(1)2 ) . . . ,W (`−1) 2 ),W (`) 1 ) − f (`)(f (`−1)(. . . f (1)(X,W(1)2 ) . . . ,W (`−1) 2 ),W (`) 1 )‖F
  (72)
  By the Lipschitz continuity of f (`)(·, ·) we have
  ‖H(`)1 −H (`) 2 ‖F ≤ C L−1 H CW ‖W (1) 1 −W (1) 2 ‖F + . . .+ CW ‖W (`) 1 −W (`) 2 ‖F
  ≤ [
  max j∈{0,...,L−1}
  CjHCW ]( L∑ j=1 ‖W(j)1 −W (j) 2 ‖F ) (73) Let define UmaxC is defined as
  UmaxC = max j∈{0,...,L−1}
  CjHCW (74)
  Plugging it back we have
  ‖G(`)1 −G (`) 2 ‖F ≤ (LUmaxLUmaxC + UmaxL) ( L∑ j=1 ‖W(j)1 −W (j) 2 ‖F ) (75)
  Summing both size from ` = 1 to ` = L we have
  ‖∇L(θ1)−∇L(θ2)‖F = L∑ `=1 ‖G(`)1 −G (`) 2 ‖F
  ≤ L(LUmaxLUmaxC + UmaxL) ( L∑ j=1 ‖W(j)1 −W (j) 2 ‖F ) ≤ L(LUmaxLUmaxC + UmaxL)‖θ1 − θ2‖F
  (76)
  G PROOF OF THEOREM 1
  By bias-variance decomposition, we can decompose the mean-square error of stochastic gradient as L∑ `=1 E[‖G̃(`) −G(`)‖2F] = L∑ `=1 [ E[‖E[G̃(`)]−G(`)‖2F]︸ ︷︷ ︸
  bias E[‖b‖2F]
  +E[‖G̃(`) − E[G̃(`)]‖2F]︸ ︷︷ ︸ variance E[‖n‖2F]
  ] (77)
  Therefore, we have to explicitly define the computation of E[G̃(`)], which requires computing D̄(L+1) = E[D̃(L+1)], D̄(`) = E[D̃(`)], and Ḡ(`) = E[G̃(`)].
  Let defined a general form of the sampled Laplacian matrix L̃(`) ∈ RN×N as
  L̃ (`) i,j =
  { Li,j αi,j if i ∈ B(`) and j ∈ B(`−1)
  0 otherwise (78)
  where αi,j is the weighted constant depends on the sampling algorithms.
  The expectation of L̃(`)i,j is computed as E[L̃(`)i,j ] = Ei∈B(`) [ Ej∈B(`−1) [L̃ (`) i,j | i ∈ B (`)] ]
  (79)
  In order to compute the expectation of SGCN’s node embedding matrices, let define the propagation matrix P(`) ∈ RN×N as
  P (`) i,j = Ei∈B(`)
  [ L̃
  (`) i,j | i ∈ B
  (`) ]
  (80)
  where the expectation is taken over row indices i. The above equation implies that under the condition that knowing the ith node is in B(`), we have P (`)i,j = L̃i,j , ∀j = {1, . . . , N}. Let consider the mean-aggregation for the ith node as
  x (`) i = σ ( N∑ j=1 L̃ (`) i,jx (`−1) j ) (81)
  Then, under the condition ith node is in B(`), we can replace L̃(`)i,j by P (`) i,j , which gives us
  x (`) i = σ ( N∑ j=1 P (`) i,j x (`−1) j ) (82)
  As a result, we can write the expectation of x(`)i with respect to the indices i as
  Ei∈B(`) [x (`) i | i ∈ B (`)] = Ei∈B(`) [ σ ( N∑ j=1 L̃ (`) i,jx (`−1) j ) | i ∈ B(`) ]
  = Ei∈B(`) [ σ ( N∑ j=1 P (`) i,j x (`−1) j ) | i ∈ B(`) ]
  = σ ( N∑ j=1 P (`) i,j x (`−1) j
  ) (83)
  Then define H̄(`) ∈ RN×d` as the node embedding of using full-batch but a subset of neighbors for neighbor aggregation, i.e.,
  H̄(`) = σ(P(`)H̄(`−1)W(`)) (84)
  where all rows in H̄(`) are non-zero vectors.
  Using the notations defined above, we can compute D̄(L+1) ∈ RN×dL , Ḡ(`) ∈ Rd`−1×d` , and D̄(`) ∈ RN×d`−1 as
  D̄(L+1) = E [∂Loss(H̄(L))
  ∂H̄(L)
  ] ∈ RN×dL , d̄i = 1
  N
  ∂Loss(h̄(L)i , yi)
  ∂h̄ (L) i
  ∈ RdL (85)
  and D̄(`) = [ ∇H f̄ (`)(D̄(`+1), H̄(`−1),W(`)) := [L]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]> ] (86) and
  Ḡ(`) = [ ∇W f̄ (`)(D̄(`+1), H̄(`−1),W(`)) := [LH̄(`−1)]> ( D̄(`+1)◦σ′(P(`)H̄(`−1)W(`)) )] (87)
  As a result, we can represent Ḡ(`) = E[G̃(`)] as
  Ḡ(`) = ∇W f̄ (`)(∇H f̄ (`+1)(. . .∇H f̄ (L)(D̄(L+1), H̄(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`))) (88)
  G.1 SUPPORTING LEMMAS
  We derive the upper-bound of the bias and variance of the stochastic gradient in the following lemmas.
  Lemma 2 (Upper-bound on variance). We can upper-bound the variance of stochastic gradient in SGCN as
  L∑ `=1 E[‖G̃(`) − E[G̃(`)]‖2F] ≤ L∑ `=1 O(E[‖L̃(`) −P(`)‖2F]) +O(‖P(`) − L‖2F) (89)
  Proof. By definition, we can write the variance in SGCN as
  E[‖G̃(`) − E[G̃(`)]‖2F]
  = E[‖∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̃(L+1), H̃(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`)) −∇W f̄ (`)(∇H f̄ (`+1)(. . .∇H f̄ (L)(D̄(L+1), H̄(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1)),W(`)‖2F]
  ≤ (L+ 1)E[‖∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̃(L+1), H̃(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`))
  −∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̄(L+1), H̃(L−1)W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`))‖2F]
  + (L+ 1)E[‖∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̄(L+1), H̃(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`))
  −∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̄ (L)(D̄(L+1), H̄(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`))‖2F] + . . .
  + (L+ 1)E[‖∇W f̃ (`)(∇H f̃ (`+1)(D̄(`+2), H̃(`),W(`+1)), H̃(`−1),W(`))
  −∇W f̃ (`)(∇H f̄ (`+1)(D̄(`+2), H̄(`),W(`+1)), H̃(`−1),W(`))‖2F]
  + (L+ 1)E[‖∇W f̃ (`)(D̄(`+1), H̃(`−1),W(`))−∇W f̄ (`)(D̄(`+1), H̄(`−1),W(`))‖2F]
  ≤ (L+ 1)L2WL 2(L−`−1) H E[‖D̃ (L+1) − D̄(L+1)‖2F]
  + (L+ 1)L2WL 2(L−`−2) H E[‖∇H f̃ (L)(D̄(L+1), H̃(L−1)W(L))−∇H f̄ (L)(D̄(L+1), H̄(L−1)W(L))‖2F] + . . .
  + (L+ 1)L2WE[‖∇H f̃ (`+1)(D̄(`+2), H̃(`),W(`+1))−∇H f̄ (`+1)(D̄(`+2), H̄(`),W(`+1))‖2F]
  + (L+ 1)E[‖∇W f̃ (`)(D̄(`+1), H̃(`−1),W(`))−∇W f̄ (`)(D̄(`+1), H̄(`−1),W(`))‖2F] (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
  E[‖D̃(L+1) − D̄(L+1)‖2F] (91)
  • The difference of gradient with respect to the input node embedding matrix at each graph convolutional layer
  E[‖∇H f̃ (`+1)(D̄(`+2), H̃(`),W(`+1))−∇H f̄ (`+1)(D̄(`+2), H̄(`),W(`+1))‖2F] (92) • The difference of gradient with respect to the weight matrix at each graph convolutional
  layer
  E[‖∇W f̃ (`)(D̄(`+1), H̃(`−1),W(`))−∇W f̄ (`)(D̄(`+1), H̄(`−1),W(`))‖2F] (93) First, Let consider the upper-bound of Eq. 91.
  E[‖D̃(L+1) − D̄(L+1)‖2F] = E[‖ ∂Loss(H̃(L),y) ∂H̃(L) − ∂Loss(H̄ (L),y) ∂H̄(L) ‖2F]
  ≤ L2lossE[‖H̃(L) − H̄(L)‖2F]
  ≤ L2lossE[‖σ(L̃(L)H̃(L−1)W(L))− σ(P(`)H̄(L−1)W(L))‖2F]
  ≤ L2lossC2σB2WE[‖L̃(L)H̃(L−1) −P(L)H̄(L−1)‖2F]
  ≤ L2lossC2σB2WB2HE[‖L̃(L) −P(L)‖2F]
  (94)
  Then, let consider the upper-bound of Eq. 92.
  E[‖∇H f̃ (`)(D̄(`+1), H̃(`−1),W(`))−∇H f̄ (`)(D̄(`+1), H̄(`−1),W(`))]‖2F] = E[‖[L̃(`)]> ( D̄(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)) ) [W(`)]> − [L]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]>‖2F]
  ≤ 2E[‖[L̃(`)]> ( D̄(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)) ) [W(`)]> − [L̃(`)]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]>‖2F]
  + 2E[‖[L̃(`)]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]> − [L]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]>‖2F]
  ≤ 2B2LAB2DB4WL2σE[‖L̃(`)H̃(`−1) −P(`)H̄(`−1)‖2F] + 2B2DC2σB2WE[‖L̃(`) − L‖2F]
  ≤ 2B2LAB2DB2HB4WL2σE[‖L̃(`) −P(`)‖2F] + 2B2DC2σB2WE[‖L̃(`) −P(`) + P(`) − L‖2F] ≤ 2 ( B2LAB 2 DB 2 HB 4 WL 2 σ + 2B 2 DC 2 σB 2 W ) E[‖L̃(`) −P(`)‖2F] + 4B2DC2σB2WE[‖P(`) − L‖2F]
  ≤ O(E[‖L̃(`) −P(`)‖2F]) +O(E[‖P(`) − L‖2F]) (95)
  Finally, let consider the upper-bound of Eq. 93.
  E[‖∇W f̃ (`)(D̄(`+1), H̃(`−1),W(`))−∇W f̄ (`)(D̄(`+1), H̄(`−1),W(`))‖2F] ≤ E[‖[L̃(`)H̃(`−1)]> ( D̄(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)) ) − [LH̄(`−1)]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) ‖2F]
  ≤ 2E[‖[L̃(`)H̃(`−1)]> ( D̄(`+1) ◦ σ′(L̃(`)H̃(`−1)W(`)) ) − [L̃(`)H̃(`−1)]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) ‖2F]
  + 2E[‖[L̃(`)H̃(`−1)]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) − [LH̄(`−1)]> ( D̄(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) ‖2F]
  ≤ 2B2LAB2HB2DB2WL2σE[‖L̃(`)H̃(`−1) −P(`)H̄(`−1)‖2F] + 2B2DC2σE[‖L̃(`)H̃(`−1) − LH̄(`−1)‖2F] ≤ 2 ( B2LAB 2 HB 2 DB 2 WL 2 σ +B 2 DC 2 σ ) E[‖L̃(`)H̃(`−1) −P(`)H̄(`−1)‖2F]
  + 2B2DB 2 HC 2 σE[‖P(`) − L‖2F] ≤ 2 ( B2LAB 4 HB 2 DB 2 WL 2 σ +B 2 HB 2 DC 2 σ ) E[‖L̃(`) −P(`)‖2F]
  + 2B2DB 2 HC 2 σE[‖P(`) − L‖2F]
  ≤ O(E[‖L̃(`) −P(`)‖2F]) +O(E[‖P(`) − L‖2F]) (96)
  Combining the result from Eq. 91, 92, 93 we have
  E[‖G̃(`) − E[G̃(`)]‖2F] ≤ O(E[‖L̃(`) −P(`)‖2F]) + . . .+O(E[‖L̃(L) −P(L)‖2F]) +O(E[‖P(`) − L‖2F]) + . . .+O(E[‖P(L) − L‖2F])
  (97)
  Lemma 3 (Upper-bound on bias). We can upper-bound the bias of stochastic gradient in SGCN as L∑ `=1 E[‖E[G̃(`)]−G(`)‖2F] ≤ L∑ `=1 O(‖P(`) − L‖2F) (98)
  Proof. By definition, we can write the bias of stochastic gradient in SGCN as
  E[‖E[G̃(`)]−G(`)‖2F] ≤ E[‖∇W f̄ (`)(∇H f̄ (`+1)(. . .∇H f̄ (L)(D̄(L+1), H̄(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`)) −∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1),H(L−1),W(L)) . . . ,H(`),W(`+1)),H(`−1),W(`))‖2F]
  ≤ (L+ 1)E[‖∇W f̄ (`)(∇H f̄ (`+1)(. . .∇H f̄ (L)(D̄(L+1), H̄(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`)) −∇W f̄ (`)(∇H f̄ (`+1)(. . .∇H f̄ (L)(D(L+1), H̄(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`))‖2F]
  + (L+ 1)E[‖∇W f̄ (`)(∇H f̄ (`+1)(. . .∇H f̄ (L)(D(L+1), H̄(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`)) −∇W f̄ (`)(∇H f̄ (`+1)(. . .∇Hf (L)(D(L+1),H(L−1),W(L)) . . . , H̄(`),W(`+1)), H̄(`−1),W(`))‖2F] + . . . + (L+ 1)E[‖∇W f̄ (`)(∇H f̄ (`+1)(D(`+2), H̄(`),W(`+1)), H̄(`−1),W(`)) −∇W f̄ (`)(∇Hf (`+1)(D(`+2),H(`),W(`+1)), H̄(`−1),W(`))‖2F] + (L+ 1)E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]
  ≤ (L+ 1)L2WL 2(L−`−1) H E[‖D̄ (L+1) −D(L+1)‖2F]
  + (L+ 1)L2WL 2(L−`−2) H E[‖∇H f̄ (L)(D(L+1), H̄(L−1),W(L))−∇Hf (L)(D(L+1),H(L−1),W(L))‖2F] + . . . + (L+ 1)L2WE[‖∇H f̄ (`+1)(D(`+2), H̄(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F] + (L+ 1)E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]
  (99)
  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
  E[‖D̄(L+1) −D(L+1)‖2F] (100) • The difference of gradient with respect to the input node embedding matrix at each graph
  convolutional layer
  E[‖∇H f̄ (`+1)(D(`+2), H̄(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F] (101) • The difference of gradient with respect to the weight matrix at each graph convolutional
  layer
  E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F] (102) Firstly, let consider the upper-bound of Eq. 100.
  E[‖D̄(L+1) −D(L+1)‖2F] = E[‖ ∂Loss(H̄(L),y) ∂H̄(L) − ∂Loss(H (L),y) ∂H(L) ‖2F]
  ≤ L2lossE[‖H̄(L) −H(L)‖2F] (103)
  The upper-bound for E[‖H̄(`) −H(`)‖2F] as
  ‖H̄(`) −H(`)‖2F = ‖σ(P(`)H̄(`−1)W(`))− σ(LH(`−1)W(`))‖2F ≤ C2σB2W ‖P(`)H̄(`−1) − LH̄(`−1) + LH̄(`−1) − LH(`−1)‖2F ≤ 2C2σB2WB2H‖P(`) − L‖2F + 2C2σB2WB2LA‖H̄(`−1) −H(`−1)‖2F ≤ O(‖P(1) − L‖2F) + . . .+O(‖P(`) − L‖2F)
  (104)
  Therefore, we have
  E[‖D̄(L+1) −D(L+1)‖2F] ≤ O(‖P(1) − L‖2F) + . . .+O(‖P(L) − L‖2F) (105)
  Then, let consider the upper-bound of Eq. 101.
  E[‖∇H f̄ (`)(D(`+1), H̄(`−1),W(`))−∇Hf (`)(D(`+1),H(`−1),W(`))‖2F] = E[‖[L]> ( D(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]> − [L]> ( D(`+1) ◦ σ′(LH(`−1)W(`)) ) [W(`)]>‖2F]
  ≤ B2LAB2DB4WL2σE[‖P(`)H̄(`−1) − LH̄(`−1) + LH̄(`−1) − LH(`−1)‖2F] ≤ 2B2LAB2DB4WL2σB2HE[‖P(`) − L‖2F] + 2B4LAB2DB4WL2σE[‖H̄(`−1) −H(`−1)‖2F] ≤ O(‖P(1) − L‖2F) + . . .+O(‖P(`) − L‖2F)
  (106)
  Finally, let consider the upper-bound of Eq. 102.
  E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F] = E[‖[LH̄(`−1)]> ( D(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) − [LH(`−1)]> ( D(`+1) ◦ σ′(L(`)H(`−1)W(`)) ) ‖2F]
  ≤ 2E[‖[LH̄(`−1)]> ( D(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) − [LH(`−1)]> ( D(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) ‖2F]
  + 2E[‖[LH(`−1)]> ( D(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) − [LH(`−1)]> ( D(`+1) ◦ σ′(L(`)H(`−1)W(`)) ) ‖2F]
  ≤ 2B2DC2σB2LAE[‖H̄(`−1) −H(`−1)‖2F] + 2B2LAB 2 HB 2 DL 2 σB 2 WE[‖P(`)H̄(`−1) − LH̄(`−1) + LH̄(`−1) − LH(`−1)‖2F]
  ≤ 2 ( B2DC 2 σB 2 LA +B 4 LAB 2 HB 2 DL 2 σB 2 W ) E[‖H̄(`−1) −H(`−1)‖2F]
  + 2B2LAB 4 HB 2 DL 2 σB 2 WE[‖P(`) − L‖2F]
  ≤ O(‖P(1) − L‖2F) + . . .+O(‖P(`) − L‖2F) (107)
  Combining the result from Eq. 100, 101, 102 we have
  E[‖E[G̃(`)]−G(`)‖2F] ≤ O(E[‖P(1) − L‖2F]) + . . .+O(E[‖P(L) − L‖2F]) (108)
  G.2 REMAINING STEPS TOWARD THEOREM 1
  By the smoothness of L(θt), we have
  L(θt+1) ≤ L(θt) + 〈∇L(θt),θt+1 − θt〉+ Lf 2 ‖θt+1 − θt‖2F
  = L(θt)− η〈∇L(θt),∇L̃(θt)〉+ η2Lf
  2 ‖∇L̃(θt)‖2F
  (109)
  Let Ft = {{B(`)1 }L`=1, . . . , {B (`) t−1}L`=1}. Note that the weight parameters θt is a function of history of the generated random process and hence is random. Taking expectation on both sides condition
  on Ft and using η < 1/Lf we have E[L(θt+1)|Ft]
  ≤ L(θt)− η〈∇L(θt),E[∇L̃(θt)|Ft]〉+ η2Lf
  2
  ( E[‖∇L̃(θt)− E[∇L̃(θt)|Ft]‖2F|Ft] + E[‖E[∇L̃(θt)|Ft]‖2F|Ft] ) = L(θt)− η〈∇L(θt),∇L(θt) + E[bt|Ft]〉+
  η2Lf 2
  ( E[‖nt‖2F|Ft] + ‖∇L(θt) + E[bt|Ft]‖2F ) ≤ L(θt) + η
  2
  ( − 2〈∇L(θt),∇L(θt) + E[bt|Ft]〉+ ‖∇L(θt) + E[bt|Ft]‖2F ) + η2Lf
  2 E[‖nt‖2F|Ft]
  = L(θt) + η
  2
  ( − ‖∇L(θt)‖2F + E[‖bt‖2F|Ft] ) + η2Lf
  2 E[‖nt‖2F|Ft]
  (110)
  Denote ∆b as the upper bound of bias of stochasitc gradient as shown in Lemma 3 and Denote ∆n as the upper bound of bias of stochastic gradient as shown in Lemma 2. Plugging in the upper bound of bias and variance, taking expectation over Ft, and rearranging the term we have
  E[‖∇L(θt)‖2F] ≤ 2
  η
  ( E[L(θt)]− E[L(θt+1)] ) + ηLf∆n + ∆b (111)
  Summing up from t = 1 to T , rearranging we have
  1
  T T∑ t=1 E[‖∇L(θt)‖2F] ≤ 2 ηT T∑ t=1 (E[L(θt)]− E[L(θt+1)]) + ηLf∆n + ∆b
  ≤ (a)
  2
  ηT (L(θ1)− L(θ?)) + ηLf∆n + ∆b
  (112)
  where the inequality (a) is due to L(θ?) ≤ E[L(θT+1)]. By selecting learning rate as η = 1/ √ T , we have
  1
  T T∑ t=1 E[‖∇L(θt)‖2F] ≤ 2(L(θ1)− L(θ?))√ T + Lf∆n√ T + ∆b (113)
  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 an important role in the analysis of the upperbound of the bias term for the stochastic gradient. Suppose the input node embedding matrix for the ` GCN layer as H̄(`−1)t , the forward propagation for the `th layer in SGCN+ is defined as
  f̃ (`)(H̄ (`−1) t ,W (`)) = [ Z̄ (`) t := P (`)H̄ (`) t W (`) ]
  (114)
  and the forward propagation for the `th layer in FullGCN is defined as
  f (`)(H̄ (`−1) t ,W (`)) = LH̄ (`−1) t W (`) (115)
  In the following, we derive the upper-bound of
  E[‖f̃ (`)(H̄(`−1)t ,W(`))− f (`)(H̄ (`−1) t ,W (`))‖2F] = E[‖LH̄ (`−1) t W (`) t − Z̄ (`) t ‖2F] (116)
  Lemma 4. Let denote E ∈ [0, T/K − 1] as the current epoch, let t be the current step. Therefore, for any t ∈ {EK + 1, . . . , EK +K}, we have
  E[‖LH̄(`−1)t W (`) t − Z̄ (`) t ‖2F] ≤ η2K ×O ( |E[‖P(`)‖2F]− ‖L‖2F|+ . . .+ |E[‖P(1)‖2F]− ‖L‖2F| )) (117)
  Proof.
  ‖LH̄(`−1)t W (`) t − Z̄ (`) t ‖2F
  = ‖[LH̄(`−1)t W (`) t − LH̄ (`−1) t−1 W (`) t−1] + [LH̄ (`−1) t−1 W (`) t−1 − Z̄ (`) t−1]− [Z̄ (`) t − Z̄ (`) t−1]‖2F = ‖LH̄(`−1)t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F + ‖LH̄ (`−1) t−1 W (`) t−1 − Z̄ (`) t−1‖2F + ‖Z̄ (`) t − Z̄ (`) t−1‖2F
  + 2〈LH̄(`−1)t W (`) t − LH̄ (`−1) t−1 W (`) t−1,LH̄ (`−1) t−1 W (`) t−1 − Z̄ (`) t−1〉 − 2〈LH̄(`−1)t W (`) t − LH̄ (`−1) t−1 W (`) t−1, Z̄ (`) t − Z̄ (`) t−1〉 − 2〈LH̄(`−1)t−1 W (`) t−1 − Z̄ (`) t−1, Z̄ (`) t − Z̄ (`) t−1〉
  (118)
  Recall that by the update rule, we have
  Z̄ (`) t − Z̄ (`) t−1 = P (`)H̄ (`−1) t W (`) t −P(`)H̄ (`−1) t−1 W (`) t−1 (119)
  and
  E[Z̄(`)t − Z̄ (`) t−1|Ft] = LH̄ (`−1) t W (`) t − LH̄ (`−1) t−1 W (`) t−1 (120)
  Taking expectation on both side condition on Ft, we have
  E[‖LH̄(`−1)t W (`) t − Z̄ (`) t ‖2F|Ft] ≤ ‖LH̄ (`−1) t−1 W (`) t−1 − Z̄ (`) t−1‖2F + E[‖Z̄ (`) t − Z̄ (`) t−1‖2F|Ft]
  − ‖LH̄(`−1)t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F
  (121)
  Then take the expectation over Ft, we have
  E[‖LH̄(`−1)t W (`) t − Z̄ (`) t ‖2F] ≤ ‖LH̄ (`−1) t−1 W (`) t−1 − Z̄ (`) t−1‖2F + E[‖Z̄ (`) t − Z̄ (`) t−1‖2F]
  − ‖LH̄(`−1)t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F
  (122)
  Since we know t ∈ {EK + 1, . . . , EK +K}, we can denote t = EK + k, k ≤ K such that
  E[‖LH̄(`−1)t W (`) t − Z̄ (`) t ‖2F] = E[‖LH̄ (`−1) EK+kW (`) EK+k − Z̄ (`) EK+k‖ 2 F]
  = E[‖LH̄(`−1)EK W (`) EK − Z̄ (`) EK‖ 2 F]︸ ︷︷ ︸
  (A)
  + EK+K∑ t=EK+1 ( E[‖Z̄(`)t − Z̄ (`) t−1‖2F]− ‖LH̄ (`−1) t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F ) (123)
  Knowing that we are using all neighbors at the snapshot step (t mod K) = 0, we have (A) = 0. As a result, we have
  E[‖LH̄(`−1)EK+kW (`) EK+k − Z̄ (`) EK+k‖ 2 F]
  ≤ EK+K∑ t=EK+1 ( E[‖Z̄(`)t − Z̄ (`) t−1‖2F]− ‖LH̄ (`−1) t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F︸ ︷︷ ︸
  (B) ) (124) Let take a closer look at term (B).
  E[‖Z̄(`)t − Z̄ (`) t−1‖2F]− ‖LH̄ (`−1) t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F
  = E[‖P(`)H̄(`−1)t W (`) t −P(`)H̄ (`−1) t−1 W (`) t−1‖2F]− ‖LH̄ (`−1) t W (`) t − LH̄ (`−1) t−1 W (`) t−1‖2F ≤ E[‖P(`)(H̄(`−1)t W (`) t − H̄ (`−1) t−1 W (`) t−1)‖2F]− ‖L(H̄ (`−1) t W (`) t − H̄ (`−1) t−1 W (`) t−1)‖2F
  ≤ ( E[‖P(`)‖2F]− ‖L‖2F ) E[‖H̄(`−1)t W (`) t − H̄ (`−1) t−1 W
  (`) t−1‖2F]︸ ︷︷ ︸
  (C)
  (125)
  Let take a closer look at term (C).
  E[‖H̄(`−1)t W (`) t − H̄ (`−1) t−1 W (`) t−1‖2F]
  = E[‖H̄(`−1)t W (`) t − H̄ (`−1) t W (`) t−1 + H̄ (`−1) t W (`) t−1 − H̄ (`−1) t−1 W (`) t−1‖2F]
  ≤ 2B2HE[‖W (`) t −W (`) t−1‖2F] + 2E[‖H̄ (`−1) t W (`) t−1 − H̄ (`−1) t−1 W (`) t−1‖2F]
  (126)
  By induction, we have
  E[‖H̄(`−1)t W (`) t − H̄ (`−1) t−1 W (`) t−1‖2F]
  ≤ 2B2HE[‖W (`) t −W (`) t−1‖2F] + 22B4HE[‖W (`−1) t −W (`−1) t−1 ‖2F] + . . .
  + 2`B2`HE[‖W (1) t −W (1) t−1‖2F]
  (127)
  By the update rule of weight matrices, we know
  E[‖W(`)t −W (`) t−1‖2F] = η2E[‖Ḡ (`) t−1‖2F] (128)
  Therefore, we have
  E[‖LH̄(`−1)EK+kW (`) EK+k − Z̄ (`) EK+k‖ 2 F]
  ≤ EK+K∑ t=EK+1 η2O ( |E[‖P(`)‖2F]− ‖L‖2F| × E[‖Ḡ (`) t−1‖2F] + . . .+ |E[‖P(1)‖2F]− ‖L‖2F| × E[‖Ḡ (1) t−1‖2F] )) (129)
  By the definition of Ḡ(`)t−1, we have that
  E[‖Ḡ(`)t−1‖2F] ≤ B2LAB2HB2DC2σ (130)
  Plugging it back, we have
  E[‖LH̄(`−1)t W (`) t − Z̄ (`) t ‖2F]
  := E[‖LH̄(`−1)EK+kW (`) EK+k − Z̄ (`) EK+k‖ 2 F] ≤ η2K ×O ( |E[‖P(`)‖2F]− ‖L‖2F|+ . . .+ |E[‖P(1)‖2F]− ‖L‖2F| )) (131)
  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 in SGCN.
  Lemma 5 (Upper-bound on bias). We can upper-bound the bias of stochastic gradient in SGCN+ as L∑ `=1 E[‖E[G̃(`)]−G(`)‖2F] ≤ η2K L∑ `=1 O ( |E[‖P(`)‖2F]− ‖L‖2F| ) (132)
  Proof. From the decomposition of bias as shown in previously in Eq. 108, we have
  E[‖E[G̃(`)]−G(`)‖2F]
  ≤ (L+ 1)L2WL 2(L−`−1) H E[‖ ¯D(L+1) −D (L+1)‖2F]
  + (L+ 1)L2WL 2(L−`−2) H E[‖∇H f̄ (L)(D(L+1), H̄(L−1),W(L))−∇Hf (L)(D(L+1),H(L−1),W(L))‖2F] + . . . + (L+ 1)L2WE[‖∇H f̄ (`+1)(D(`+2), H̄(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F] + (L+ 1)E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]
  (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
  E[‖D̄(L+1) −D(L+1)‖2F] (134) • The difference of gradient with respect to the input node embedding matrix at each graph
  convolutional layer
  E[‖∇H f̄ (`+1)(D(`+2), H̄(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F] (135) • The difference of gradient with respect to the weight matrix at each graph convolutional
  layer
  E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F] (136) Firstly, let consider the upper-bound of Eq. 134.
  E[‖D̄(L+1) −D(L+1)‖2F] = E[‖ ∂Loss(H̄(L),y) ∂H̄(L) − ∂Loss(H (L),y) ∂H(L) ‖2F]
  ≤ L2lossE[‖H̄(L) −H(L)‖2F] ≤ L2lossC2σE[‖Z̄(L) − Z(L)‖2F]
  (137)
  We can decompose E[‖Z̄(L) − Z(L)‖2F] as
  E[‖Z̄(L) − Z(L)‖2F] = E[‖Z̄(L) − LH(L−1)W(L)‖2F] ≤ 2E[‖Z̄(L) − LH̄(L−1)W(L)‖2F] + 2E[‖LH̄(L−1)W(L) − LH(L−1)W(L)‖2F] ≤ 2E[‖Z̄(L) − LH̄(L−1)W(L)‖2F] + 2B2LAB2WC2σE[‖Z̄(L−1) − Z(L−1)‖2F]
  ≤ L∑ `=1 O(E[‖Z̄(`) − LH̄(`−1)W(`)‖2F])
  (138)
  Using result from Lemma 4, we have
  E[‖D̄(L+1) −D(L+1)‖2F] ≤ L∑ `=1 η2K ×O ( |E[‖P(`)‖2F]− ‖L‖2F| ) (139)
  Then, let consider the upper-bound of Eq. 135.
  E[‖∇H f̄ (`)(D(`+1), H̄(`−1),W(`))−∇Hf (`)(D(`+1),H(`−1),W(`))‖2F] = E[‖[L]> ( D(`+1) ◦ σ′(P(`)H̄(`−1)W(`)) ) [W(`)]> − [L]> ( D(`+1) ◦ σ′(LH(`−1)W(`)) ) [W(`)]>‖2F]
  ≤ B2LAB2DB2WL2σE[‖Z̄ (`) t − LH(`−1)W(`)‖2F]
  ≤ 2B2LAB2DB2WL2σE[‖Z̄ (`) t − LH̄(`−1)W(`)‖2F]
  + 2B2LAB 2 DB 2 WL 2 σ E[‖LH̄(`−1)W(`) − LH(`−1)W(`)‖2F]︸ ︷︷ ︸
  (A)
  (140) where Z̄(`)t = Z̄ (`) t−1 + P (`)H̄ (`−1) t W (`) t −P(`)H̄ (`−1) t−1 W (`) t−1.
  Let take a closer look at term (A), we have
  E[‖LH̄(`−1)W(`) − LH(`−1)W(`)‖2F] ≤ B2LAB2WC2σE[‖Z̄(`−2) − LH(`−2)W(`−1)‖2F] ≤ 2B2LAB2WC2σE[‖Z̄(`−2) − LH̄(`−2)W(`−1)‖2F]
  + 2B2LAB 2 WC 2 σE[‖LH̄(`−2)W(`−1) − LH(`−2)W(`−1)‖2F]
  (141)
  Therefore, by induction we have
  E[‖∇H f̄ (`)(D(`+1), H̄(`−1),W(`))−∇Hf (`)(D(`+1),H(`−1),W(`))‖2F] ≤ O(E[‖Z̄(`) − LH̄(`−1)W(`)‖2F]) +O(E[‖Z̄(`−1) − LH̄(`−2)W(`−1)‖2F]) + . . .
  +O(E[‖Z̄(2) − LH̄(1)W(2)‖2F]) +O(E[‖Z̄(1) − LXW(1)‖2F])
  (142)
  Using result from Lemma 4, we have
  E[‖∇H f̄ (`)(D(`+1), H̄(`−1),W(`))−∇Hf (`)(D(`+1),H(`−1),W(`))‖2F] ≤ η2K ×O ( |E[‖P(`)‖2F]− ‖L‖2F| ) + . . .+ η2K ×O ( |E[‖P(1)‖2F]− ‖L‖2F| ) (143) Finally, let consider the upper-bound of Eq. 136.
  E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F] = E[‖[LH̄(`−1)]> ( D(`+1) ◦ σ′(Z̄(`)) ) − [LH(`−1)]> ( D(`+1) ◦ σ′(L(`)H(`−1)W(`)) ) ‖2F]
  ≤ 2E[‖[LH̄(`−1)]> ( D(`+1) ◦ σ′(Z̄(`)) ) − [LH(`−1)]> ( D(`+1) ◦ σ′(Z̄(`)) ) ‖2F]
  + 2E[‖[LH(`−1)]> ( D(`+1) ◦ σ′(Z̄(`)) ) − [LH(`−1)]> ( D(`+1) ◦ σ′(L(`)H(`−1)W(`)) ) ‖2F]
  ≤ 2B2DC2σB2LA E[‖H̄(`−1) −H(`−1)‖2F]︸ ︷︷ ︸ (B) +4B2LAB 2 HB 2 DL 2 σE[‖Z̄(`) − LH̄(`−1)W(`)‖2F]
  + 4B2LAB 2 HB 2 DL 2 σE[‖LH̄(`−1)W(`) − LH(`−1)W(`)‖2F]
  (144) By definition, we can write the term (B) as
  E[‖H̄(`−1) −H(`−1)‖2F] ≤ C2σE[‖Z̄(`) − LH(`−2)W(`−1)‖2F] ≤ 2C2σE[‖Z̄(`) − LH̄(`−2)W(`−1)‖2F]
  + 2C2σE[‖LH̄(`−2)W(`) − LH(`−2)W(`−1)‖2F]
  (145)
  Plugging term (B) back and using Eq. 141 and Lemma 4, we have
  E[‖∇W f̄ (`)(D(`+1), H̄(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F] ≤ η2K ×O ( |E[‖P(`)‖2F]− ‖L‖2F| ) + . . .+ η2K ×O ( |E[‖P(1)‖2F]− ‖L‖2F| ) (146) Combining the result from Eq. 134, 135, 136 we have
  E[‖E[G̃(`)]−G(`)‖2F] ≤ η2K L∑ `=1 O ( |E[‖P(`)‖2F]− ‖L‖2F| ) (147)
  H.2 REMAINING STEPS TOWARD THEOREM 2
  Now we are ready to prove Theorem 2. By the smoothness of L(θt), we have
  L(θt+1) ≤ L(θt) + 〈∇L(θt),θt+1 − θt〉+ LF 2 ‖θt+1 − θt‖2
  = L(θt)− η〈∇L(θt),∇L̃(θt)〉+ η2LF
  2 ‖∇L̃(θt)‖2
  (148)
  Let Ft = {{B(`)1 }L`=1, . . . , {B (`) t−1}L`=1}. Note that the weight parameters θt is a function of history of the generated random process and hence is random. Taking expectation on both sides condition on Ft and using η < 1/LF we have E[∇L(θt+1)|Ft]
  ≤ L(θt)− η〈∇L(θt),E[∇L̃(θt)|Ft]〉+ η2LF
  2
  ( E[‖∇L̃(θt)− E[∇L̃(θt)|Ft]‖2|Ft] + E[‖E[g|Ft]‖2|Ft] ) = L(θt)− η〈∇L(θt),∇L(θt) + E[bt|Ft]〉+
  η2LF 2
  ( E[‖nt‖2|Ft] + ‖∇L(θt) + E[bt|Ft]‖2 ) ≤ L(θt) + η
  2
  ( − 2〈∇L(θt),∇L(θt) + E[bt|Ft]〉+ ‖∇L(θt) + E[bt|Ft]‖2 ) + η2LF
  2 E[‖nt‖2|Ft]
  ≤ L(θt) + η
  2
  ( − ‖∇L(θt)‖2 + E[‖bt‖2|Ft] ) + η2LF
  2 E[‖nt‖2|Ft]
  (149)
  Plugging in the upper bound of bias and variance, taking expectation over Ft, and rearranging the term we have
  E[‖∇L(θt)‖2] ≤ 2
  η
  ( E[L(θt)]− E[∇L(θt+1)] ) + ηLF∆n + η 2∆+′b (150)
  Summing up from t = 1 to T , rearranging we have
  1
  T T∑ t=1 E[‖∇L(θt)‖2] ≤ 2 ηT T∑ t=1 (E[L(θt)]− E[L(θt+1)]) + ηLF∆n + η2∆+′b
  ≤ (a)
  2
  ηT (L(θ1)− L(θ?)) + ηLF∆n + η2∆+′b
  (151)
  where the inequality (a) is due to L(θ?) ≤ E[L(θT+1)]. By selecting learning rate as η = 1/ √ T , we have
  1
  T T∑ t=1 E[‖∇L(θt)‖2] ≤ 2(L(θ1)− L(θ?))√ T + LF∆n√ T + ∆+′b T
  (152)
  I PROOF OF THEOREM 3
  I.1 SUPPORTING LEMMAS
  In the following lemma, we decompose the mean-square error of stochastic gradient at the `th layer E[‖G̃(`) −G(`)‖F 2] as the summation of
  • The difference between the gradient with respect to the last layer node embedding matrix
  E[‖D̃(L+1) −D(L+1)‖2F] (153) • The difference of gradient passing from the (`+1)th layer node embedding to the `th layer
  node embedding
  E[‖∇H f̃ (`+1)(D(`+2), H̃(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F] (154) • The difference of gradient passing from the `th layer node embedding to the `th layer
  weight matrix
  E[‖∇W f̃ (`)(D(`+1), H̃(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F] (155) Lemma 6. The mean-square error of stochastic gradient at the `th layer can be decomposed as
  E[‖G̃(`) −G(`)‖F 2]
  ≤ O(E[‖D̃(L+1) −D(L+1)‖2F])
  +O(E[‖∇H f̃ (L)(D(L+1), H̃(L−1)W(L))−∇Hf (L)(D(L+1),H(L−1)W(L))‖2F]) + . . .
  +O(E[‖∇H f̃ (`+1)(D(`+2), H̃(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F])
  +O(E[‖∇W f̃ (`)(D(`+1), H̃(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]) (156)
  Proof. By definition, we can write down the mean-square error of stochastic gradient as
  E[‖G̃(`) −G(`)‖2F]
  = E[‖∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̃(L+1), H̃(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`)) −∇W f (`)(∇Hf (`+1)(. . .∇Hf (L)(D(L+1),H(L−1),W(L)) . . . ,H(`),W(`+1)),H(`−1)),W(`)‖2F]
  ≤ (L+ 1)E[‖∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D̃(L+1), H̃(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`))
  −∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D(L+1), H̃(L−1)W(L)) . . . ,H(`),W(`+1)),H(`−1),W(`))‖2F]
  + (L+ 1)E[‖∇W f̃ (`)(∇H f̃ (`+1)(. . .∇H f̃ (L)(D(L+1), H̃(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`))
  −∇W f̃ (`)(∇H f̃ (`+1)(. . .∇Hf (L)(D(L+1),H(L−1)W(L)) . . . , H̃(`),W(`+1)), H̃(`−1),W(`))‖2F] + . . .
  + (L+ 1)E[‖∇W f̃ (`)(∇H f̃ (`+1)(D(`+2), H̃(`),W(`+1)), H̃(`−1),W(`))
  −∇W f̃ (`)(∇Hf (`+1)(D(`+2),H(`),W(`+1)), H̃(`−1),W(`))‖2F]
  + (L+ 1)E[‖∇W f̃ (`)(D(`+1), H̃(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]
  ≤ O(E[‖D̃(L+1) −D(L+1)‖2F])
  +O(E[‖∇H f̃ (L)(D(L+1), H̃(L−1)W(L))−∇Hf (L)(D(L+1),H(L−1)W(L))‖2F]) + . . .
  +O(E[‖∇H f̃ (`+1)(D(`+2), H̃(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F])
  +O(E[‖∇W f̃ (`)(D(`+1), H̃(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]) (157)
  Recall the definition of stochastic gradient for all model parameters ∇L̃(θt) = {G̃(`)t }L`=1 where G̃
  (`) t is the gradient for the `th weight matrix, i.e.,
  W (`) t = W (`) t−1 − ηG̃ (`) t−1 (158)
  In the following lemma, we derive the upper-bound on the difference of the gradient passing from the `th to (` − 1)th layer given the same inputs D(`+1)t , H̃ (`−1) t , where the backward propagation for the `th layer in SGCN++ is defined as
  ∇H f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )
  = D̃ (`) t−1 + [L̃ (`)]>(D (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t − [L̃(`)]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1]
  (159)
  and the backward propagation for the `th layer in FullGCN is defined as
  ∇Hf (`)(D(`+1)t , H̃ (`−1) t ,W (`) t ) = L >(D (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t (160)
  Lemma 7. Let suppose t ∈ {EK + 1, . . . , EK + K}. The upper-bound on the difference of the gradient with respect to the input node embedding matrix at the `th graph convolutional layer given the same input D(`+1)t and H̃ (`−1) t is defined as
  E[‖∇H f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇Hf (`)(D (`+1) t , H̃ (`−1) t ,W (`) t )‖2F]
  ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (161)
  Proof. To simplify the presentation, let us denote D̃(`)t = ∇H f̃ (`)(D (`+1) t , H̃ (`−1) t ,W (`) t ). Then, by definition we have
  D̃ (`) t := D̃ (`) t−1 + [L̃ (`)]>(D (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t − [L̃(`)]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1 (162)
  Therefore, we know that
  ‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − D̃ (`) t ‖2F
  = ‖[L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1]
  + [L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1 − D̃ (`) t−1]− [D̃ (`) t − D̃ (`) t−1]‖2F
  ≤ ‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1︸ ︷︷ ︸
  (A1)
  ‖2F
  + ‖L>(D(`+1)t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1 − D̃ (`) t−1︸ ︷︷ ︸
  A2
  ‖2F + ‖ D̃ (`) t − D̃ (`) t−1︸ ︷︷ ︸
  A3
  ‖2F
  + 2〈A1, A2〉 − 2〈A1, A3〉 − 2〈A2, A3〉
  (163)
  Taking expectation condition on Ft on both side, and using the fact that
  E[D̃(`)t − D̃ (`) t−1|Ft]
  = E[[L̃(`)]>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − [L̃(`)]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1|Ft]
  = L>(D (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1
  (164)
  the following inequality holds
  E[‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − D̃ (`) t ‖2F|Ft]
  ≤ ‖L>(D(`+1)t−1 ◦ σ′(Z̃ (`) t−1))W (`) t − D̃ (`) t−1‖2F + E[‖D̃ (`) t − D̃ (`) t−1‖2F|Ft]
  − ‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F
  (165)
  Then, taking expectation over Ft, we have
  E[‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − D̃ (`) t ‖2F]
  ≤ E[‖L>(D(`+1)t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1 − D̃ (`) t−1‖2F] + E[‖D̃ (`) t − D̃ (`) t−1‖2F]
  − [‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F]
  (166)
  Let suppose t ∈ {EK + 1, . . . , EK +K}. Then we can denote t = EK + k for some k ≤ K such that
  E[‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − D̃ (`) t ‖2F]
  E[‖L>(D(`+1)EK+k ◦ σ ′(Z̃ (`) EK+k))W (`) EK+k − D̃ (`) EK+k‖ 2 F]
  ≤ E[‖L>(D(`+1)EK ◦ σ ′(Z̃ (`) EK))W (`) EK − D̃ (`) EK‖ 2 F] + EK+K∑ t=EK+1 ( E[‖D̃(`)t − D̃ (`) t−1‖2F]
  − [‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F] ) ≤ E[‖L>(D(`+1)EK ◦ σ ′(Z̃ (`) EK))W (`) EK − D̃ (`) EK‖ 2 F]
  + EK+K∑ t=EK+1 ( E[‖[L̃(`)]>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − [L̃(`)]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F]
  − E[‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − L>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F] ) (167)
  Knowing that we are taking full-batch gradient descent when (t mod K) = 0, we have
  E[‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − D̃ (`) t ‖2F]
  ≤ EK+K∑ t=EK+1 ( |E[‖L̃(`)‖2F]− ‖L‖2F| ) × E[‖(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − (D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F]︸ ︷︷ ︸
  (B)
  (168)
  Let take closer look at term (B).
  E[‖(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − (D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F]
  ≤ 3C2σB2W ‖D (`+1) t −D (`+1) t−1 ‖2F︸ ︷︷ ︸
  (C1)
  +3B2DB 2 WL 2 σ E[‖Z̃ (`) t − Z̃ (`) t−1‖2F]︸ ︷︷ ︸
  (C2)
  +3C2σB 2 DE[‖W (`) t −W (`) t−1‖2F]
  (169)
  For term (C1) by definition we know
  ‖D(`+1)t −D (`+1) t−1 ‖2F = ‖ ( L>(D
  (`+2) t ◦ σ′(Z (`+1) t ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t−1 ) ‖2F
  ≤ 3‖ ( L>(D
  (`+2) t ◦ σ′(Z (`+1) t ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t ))W (`+1) t ) ‖2F
  + 3‖ ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t ) ‖2F
  + 3‖ ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t−1 ) ‖2F
  ≤ O(‖D(`+2)t −D (`+2) t−1 ‖2F) +O(‖Z (`+1) t − Z (`+1) t−1 ‖2F) +O(‖W (`+1) t −W (`+1) t−1 ‖2F)
  (170)
  By induction, we have
  ‖D(`+1)t −D (`+1) t−1 ‖2F ≤ O(‖D (L+1) t −D (L+1) t−1 ‖2F)︸ ︷︷ ︸
  (D1)
  +O(‖Z(`+1)t − Z (`+1) t−1 ‖2F)︸ ︷︷ ︸
  (D2)
  + . . .+O(‖Z(L)t − Z (L) t−1‖2F)
  +O(‖W(`+1)t −W (`+1) t−1 ‖2F) + . . .+O(‖W (L) t −W (L) t−1‖2F)
  (171)
  For term (D1) we have
  ‖D(L+1)t −D (L+1) t−1 ‖2F = ‖ ∂L(θt) ∂W
  (L) t
  − ∂L(θt−1) ∂W
  (L) t−1
  ‖2F
  ≤ L2lossC2σ‖Z (L) t − Z (L) t−1‖2F
  (172)
  For term (D2) we have
  ‖Z(`+1)t − Z (`+1) t−1 ‖2F ≤ C2σ‖LH (`) t W (`+1) t − LH (`) t−1W (`+1) t−1 ‖2F
  ≤ C2σB2LA‖H (`) t W (`+1) t −H (`) t−1W (`+1) t + H (`) t−1W (`+1) t −H (`) t−1W (`+1) t−1 ‖2F
  ≤ 2C2σB2LA‖H (`) t W (`+1) t −H (`) t−1W (`+1) t ‖2F + 2C2σB2LA‖H (`) t−1W (`+1) t −H (`) t−1W (`+1) t−1 ‖2F
  ≤ 2C4σB2LAB2W ‖Z (`) t − Z (`) t−1‖2F + 2C2σB2LAB2H‖W (`+1) t −W (`+1) t−1 ‖2F
  (173) By induction we have
  ‖Z(`+1)t − Z (`+1) t−1 ‖2F ≤ O(‖W (`+1) t −W (`+1) t−1 ‖2F) + . . .+O(‖W (1) t −W (1) t−1‖2F) (174)
  For term (C2) by definition we have
  E[‖Z̃(`)t − Z̃ (`) t−1‖2F] = E[‖L̃(`)H̃ (`−1) t W (`) t − L̃(`)H̃ (`−1) t−1 W (`) t−1‖2F]
  ≤ 2B2LAB2WC2σE[‖Z̃ (`−1) t − Z̃ (`−1) t−1 ‖2F] + 2B2LAB2HE[‖W (`) t −W (`) t−1‖2F]
  (175) By induction we have
  E[‖Z̃(`)t − Z̃ (`) t−1‖2F] ≤ O(E[‖W (`) t −W (`) t−1‖2F]) + . . .+O(E[‖W (1) t −W (1) t−1‖2F]) (176)
  Plugging (D1), (D2) back to (C1) and (C1), (C2), (C3) back to (B), we have
  E[‖(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − (D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F]
  ≤ EK+K∑ t=EK+1 ( O(E[‖W(1)t −W (1) t−1‖2F]) + . . .+O(E[‖W (L) t −W (L) t−1‖2F]) )
  = EK+K∑ t=EK+1 η2O ( E[‖∇L̃(θt−1)‖2F] ) (177)
  Then plugging term (B) back to Eq. 169 we conclude the proof.
  Using the previous lemma, we provide the upper-bound of Eq. 154, which is one of the three key factors that affect the mean-square error of stochastic gradient at the `th layer.
  Lemma 8. Let suppose t ∈ {EK + 1, . . . , EK + K}. The upper-bound on the difference of the gradient with respect to the input node embedding matrix at the `th graph convolutional layer given the same input D(`+1)t but different input H̃ (`−1) t ,H (`−1) t is defined as
  E[‖∇H f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇Hf (`)(D (`+1) t ,H (`−1) t ,W (`) t )‖2F]
  ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (178)
  Proof. For the gradient w.r.t. the node embedding matrices, we have
  E[‖∇H f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇Hf (`)(D (`+1) t ,H (`−1) t ,W (`) t )‖2F] = E[‖ ( D̃
  (`) t−1 + [L̃ (`)]>(D (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t − [L̃(`)]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1] ) − ( L>(D
  (`+1) t ◦ σ′(Z (`) t ))W (`) t ) ‖2F]
  ≤ 2E[‖ ( D̃
  (`) t−1 + [L̃ (`)]>(D (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t − [L̃(`)]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1] ) ︸ ︷︷ ︸
  (A) − ( L>(D
  (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t ) ‖2F]︸ ︷︷ ︸
  (A) + 2E[‖ ( L>(D
  (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t
  ) − ( L>(D
  (`+1) t ◦ σ′(Z (`) t ))W (`) t ) ‖2F]︸ ︷︷ ︸
  (B)
  (179) Let first take a closer look at term (A). Let suppose t ∈ {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
  (A) ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (180)
  Then we take a closer look at term (B). E[‖ ( L>(D
  (`+1) t ◦ σ′(Z̃ (`) t ))W (`) t
  ) − ( L>(D
  (`+1) t ◦ σ′(Z (`) t ))W (`) t ) ‖2F]
  ≤ B2LAB2DB2WL2σ E[‖Z̃ (`) t − Z (`) t ‖2F]︸ ︷︷ ︸
  (C)
  (181)
  The term (C) can be decomposed as
  E[‖Z̃(`)t − Z (`) t ‖2F]
  = E[‖Z̃(`)t − LH (`−1) t W (`) t ‖2F] ≤ 2E[‖Z̃(`)t − LH̃ (`−1) t W (`) t ‖2F] + 2E[‖LH̃ (`−1) t W (`) t − LH (`−1) t W (`) t ‖2F] ≤ 2E[‖Z̃(`)t − LH̃ (`−1) t W (`) t ‖2F] + 2B2LAB2WC2σE[‖Z̃ (`−1) t − Z (`−1) t ‖2F]
  (182)
  By induction, we have
  E[‖Z̃(`)t −Z (`) t ‖2F] ≤ O ( E[‖Z̃(`)t − LH̃ (`−1) t W
  (`) t ‖2F]︸ ︷︷ ︸
  (D)
  ) +. . .+O ( E[‖Z̃(1)t −LXW (1) t ‖2F] ) (183)
  The upper-bound for term (D) is similar to one we have in the proof of Lemma 4.
  ‖LH̃(`−1)t W (`) t − Z̃ (`) t ‖2F
  = ‖[LH̃(`−1)t W (`) t − LH̃ (`−1) t−1 W (`) t−1] + [LH̃ (`−1) t−1 W (`) t−1 − Z̃ (`) t−1]− [Z̃ (`) t − Z̃ (`) t−1]‖2F = ‖LH̃(`−1)t W (`) t − LH̃ (`−1) t−1 W (`) t−1‖2F + ‖LH̃ (`−1) t−1 W (`) t−1 − Z̃ (`) t−1‖2F + ‖Z̃ (`) t − Z̃ (`) t−1‖2F
  + 2〈LH̃(`−1)t W (`) t ,LH̃ (`−1) t−1 W (`) t−1〉 − 2〈LH̃ (`−1) t W (`) t , Z̃ (`) t − Z̃ (`) t−1〉
  − 2〈Z̃(`)t − Z̃ (`) t−1,LH̃ (`−1) t−1 W (`) t−1〉
  (184)
  Taking expectation condition on Ft and using
  E[Z̃(`)t − Z̃ (`) t−1|Ft] = LH̃ (`−1) t W (`) t − LH̃ (`−1) t−1 W (`) t−1 (185)
  we have
  E[‖LH̃(`−1)t W (`) t − Z̃ (`) t ‖2F|Ft]
  = ‖LH̃(`−1)t−1 W (`) t−1 − Z̃ (`) t−1‖2F + E[‖Z̃ (`) t − Z̃ (`) t−1‖2F|Ft]− ‖LH̃ (`−1) t W (`) t − LH̃ (`−1) t−1 W (`) t−1‖2F
  (186)
  Take expectation over Ft we have
  E[‖LH̃(`−1)t W (`) t − Z̃ (`) t ‖2F]
  = ‖LH̃(`−1)t−1 W (`) t−1 − Z̃ (`) t−1‖2F + E[‖Z̃ (`) t − Z̃ (`) t−1‖2F]− ‖LH̃ (`−1) t W (`) t − LH̃ (`−1) t−1 W (`) t−1‖2F
  (187)
  Let suppose t ∈ {EK + 1, . . . , EK +K}. Then we can denote t = EK + k for some k ≤ K such that
  E[‖LH̃(`−1)t W (`) t − Z̃ (`) t ‖2F]
  = E[‖LH̃(`−1)EK+kW (`) EK+k − Z̃ (`) EK+k‖ 2 F]
  = ‖LH̃(`−1)EK W (`) EK − Z̃ (`) EK‖ 2 F + EK+K∑ t=EK+1 E[‖Z̃(`)t − Z̃ (`) t−1‖2F]− ‖LH̃ (`−1) t W (`) t − LH̃ (`−1) t−1 W (`) t−1‖2F
  ≤ EK+K∑ t=EK+1 O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖W (`) t −W (`) t−1‖2F] ) (188)
  Plugging (D) to (C) and (C) to (B), we have
  (B) ≤ ∑̀ j=1 EK+K∑ t=EK+1 O ( |E[‖L̃(j)‖2F]− ‖L‖2F| × E[‖W (j) t −W (j) t−1‖2F] ) (189)
  Combing with term (A) we have
  E[‖∇H f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇Hf (`)(D (`+1) t ,H (`−1) t ,W (`) t )‖2F]
  ≤ L∑ `=1 EK+K∑ t=EK+1 O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖W (`) t −W (`) t−1‖2F] )
  = EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (190)
  In the following lemma, we derive the upper-bound on the difference of the gradient with respect to the weight matrix at each graph convolutional layer. Suppose the input node embedding matrix for the `th GCN layer is defined as H̃(`−1)t , the gradient calculated for the `th weight matrix in SGCN++ is defined as
  ∇W f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )
  = G̃ (`) t−1 + [L̃ (`)H̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t ))− [L̃(`)H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))
  (191)
  and the backward propagation for the `th layer in FullGCN is defined as
  ∇W f (`)(D(`+1)t , H̃ (`−1) t ,W (`) t ) = [LH̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t )) (192)
  Lemma 9. Let suppose t ∈ {EK + 1, . . . , EK + K}. The upper-bound on the difference of the gradient with respect to the `th graph convolutional layer given the same input D(`+1)t and H̃ (`−1) t
  is defined as
  ‖∇W f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇W f (`)(D (`+1) t , H̃ (`−1) t ,W (`) t )‖2F
  ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (193)
  Proof. To simplify the presentation, let us denote G̃(`)t = ∇W f̃ (`)(D (`+1) t , H̃ (`−1) t ,W (`) t ). Then, by definition, we have
  G̃ (`) t = G̃ (`) t−1 + [L̃ (`)H̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t ))− [L̃(`)H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1)) (194)
  Therefore, we know that
  ‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− G̃ (`) t ‖2F = ∥∥∥[[LH̃(`−1)t ]>(D(`+1)t ◦ σ′(Z̃(`)t ))− [LH̃(`−1)t−1 ]>(D(`+1)t−1 ◦ σ′(Z̃(`)t−1))]
  + [ [LH̃
  (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))− G̃ (`) t−1
  ] − [ G̃
  (`) t − G̃ (`) t−1 ]∥∥∥2 F
  ≤ ∥∥∥ [LH̃(`−1)t ]>(D(`+1)t ◦ σ′(Z̃(`)t ))− [LH̃(`−1)t−1 ]>(D(`+1)t−1 ◦ σ′(Z̃(`)t−1))︸ ︷︷ ︸
  (A1)
  ∥∥∥2 F
  + ‖ [LH̃(`−1)t−1 ]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))− G̃ (`) t−1︸ ︷︷ ︸
  A2
  ‖2F + ‖ G̃ (`) t − G̃ (`) t−1︸ ︷︷ ︸
  A3
  ‖2F
  + 2〈A1, A2〉 − 2〈A1, A3〉 − 2〈A2, A3〉
  (195)
  Taking expectation condition on Ft on both side, and using the fact that
  E[G̃(`)t − G̃ (`) t−1|Ft] = E[[L̃(`)H̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t ))− [L̃(`)H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))|Ft]
  = [LH̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t ))− [LH̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))
  (196) we have the following inequality holds
  E[‖[LH̃(`)t ]>(D (`+1) t ◦ σ′(Z̃ (`−1) t ))− G̃ (`) t ‖2F|Ft]
  ≤ ‖[LH̃(`−1)t−1 ]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))− G̃ (`) t−1‖2F + E[‖G̃ (`) t − G̃ (`) t−1‖2F|Ft]
  − ‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [LH̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1‖2F
  (197)
  Then, taking expectation over Ft, we have
  E[‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− G̃ (`) t ‖2F]
  ≤ E[‖[LH̃(`−1)t−1 ]>(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))W (`) t−1 − G̃ (`) t−1‖2F] + E[‖G̃ (`) t − G̃ (`) t−1‖2F]
  − [‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [LH̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))‖2F]
  (198)
  Let suppose t ∈ {EK + 1, . . . , EK +K}. Then we can denote t = EK + k for some k ≤ K such that
  E[‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− G̃ (`) t ‖2F] E[‖[LH̃(`−1)EK+k] >(D (`+1) EK+k ◦ σ ′(Z̃ (`) EK+k))− G̃ (`) EK+k‖ 2 F]
  ≤ E[‖[LH̃(`−1)EK ] >(D (`+1) EK ◦ σ ′(Z̃ (`) EK))− G̃ (`) EK‖ 2 F] + EK+K∑ t=EK+1 ( E[‖G̃(`)t − G̃ (`) t−1‖2F]
  − [‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [LH̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))‖2F] ) ≤ E[‖L>(D(`+1)EK ◦ σ ′(Z̃ (`) EK))W (`) EK − D̃ (`) EK‖ 2 F]
  + EK+K∑ t=EK+1 ( E[‖[L̃(`)H̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [L̃(`)H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))‖2F]
  − E[‖[LH̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [LH̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))‖2F] ) (199)
  Knowing that we are taking full-batch gradient descent when (t mod K) = 0, we have
  E[‖L>(D(`+1)t ◦ σ′(Z̃ (`) t ))W (`) t − D̃ (`) t ‖2F]
  ≤ EK+K∑ t=EK+1 ( |E[‖L̃(`)‖2F]− ‖L‖2F| ) × E[‖[H̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))‖2F]︸ ︷︷ ︸
  (B)
  (200)
  Let take closer look at term (B).
  E[‖[H̃(`−1)t ]>(D (`+1) t ◦ σ′(Z̃ (`) t ))− [H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1))‖2F]
  ≤ 3(B2DC4σ +B2HB2DL2σ)E[‖Z̃ (`−1) t − Z̃ (`−1) t−1 ‖2F]︸ ︷︷ ︸
  (C1)
  +3B2HC 2 σ E[‖D (`+1) t −D (`+1) t−1 ‖2F]︸ ︷︷ ︸
  (C2)
  (201)
  For term (C1) by definition we know
  ‖D(`+1)t −D (`+1) t−1 ‖2F = ‖
  ( L>(D
  (`+2) t ◦ σ′(Z (`+1) t ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t−1 ) ‖2F
  ≤ 3‖ ( L>(D
  (`+2) t ◦ σ′(Z (`+1) t ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t ))W (`+1) t ) ‖2F
  + 3‖ ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t ) ‖2F
  + 3‖ ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t
  ) − ( L>(D
  (`+2) t−1 ◦ σ′(Z (`+1) t−1 ))W (`+1) t−1 ) ‖2F
  ≤ O(‖D(`+2)t −D (`+2) t−1 ‖2F) +O(‖Z (`+1) t − Z (`+1) t−1 ‖2F) +O(‖W (`+1) t −W (`+1) t−1 ‖2F)
  (202) By induction, we have
  ‖D(`+1)t −D (`+1) t−1 ‖2F ≤ O(‖D (L+1) t −D (L+1) t−1 ‖2F)︸ ︷︷ ︸
  (D1)
  +O(‖Z(`+1)t − Z (`+1) t−1 ‖2F)︸ ︷︷ ︸
  (D2)
  + . . .+O(‖Z(L)t − Z (L) t−1‖2F)
  +O(‖W(`+1)t −W (`+1) t−1 ‖2F) + . . .+O(‖W (L) t −W (L) t−1‖2F)
  (203)
  For term (D1) we have
  ‖D(L+1)t −D (L+1) t−1 ‖2F = ‖ ∂L(θt) ∂W
  (L) t
  − ∂L(θt−1) ∂W
  (L) t−1
  ‖2F
  ≤ L2lossC2σ‖Z (L) t − Z (L) t−1‖2F
  (204)
  For term (D2) we have
  ‖Z(`+1)t − Z (`+1) t−1 ‖2F ≤ C2σ‖LH (`) t W (`+1) t − LH (`) t−1W (`+1) t−1 ‖2F
  ≤ C2σB2LA‖H (`) t W (`+1) t −H (`) t−1W (`+1) t + H (`) t−1W (`+1) t −H (`) t−1W (`+1) t−1 ‖2F
  ≤ 2C2σB2LA‖H (`) t W (`+1) t −H (`) t−1W (`+1) t ‖2F + 2C2σB2LA‖H (`) t−1W (`+1) t −H (`) t−1W (`+1) t−1 ‖2F
  ≤ 2C4σB2LAB2W ‖Z (`) t − Z (`) t−1‖2F + 2C2σB2LAB2H‖W (`+1) t −W (`+1) t−1 ‖2F
  (205) By induction we have
  ‖Z(`+1)t − Z (`+1) t−1 ‖2F ≤ O(‖W (`+1) t −W (`+1) t−1 ‖2F) + . . .+O(‖W (1) t −W (1) t−1‖2F) (206)
  For term (C2) by definition we have
  E[‖Z̃(`)t − Z̃ (`) t−1‖2F] = E[‖L̃(`)H̃ (`−1) t W (`) t − L̃(`)H̃ (`−1) t−1 W (`) t−1‖2F]
  ≤ 2B2LAB2W`C2σE[‖Z̃ (`−1) t − Z̃ (`−1) t−1 ‖2F] + 2B2LAB2HE[‖W (`) t −W (`) t−1‖2F]
  (207) By induction we have
  E[‖Z̃(`)t − Z̃ (`) t−1‖2F] ≤ O(E[‖W (`) t −W (`) t−1‖2F]) + . . .+O(E[‖W (1) t −W (1) t−1‖2F]) (208)
  Plugging (D1), (D2) to (C2) and C1, (C2) to (B), we have
  E[‖[LH̃(`−1)t ]>(D (`+1) t ◦σ′(Z̃ (`) t ))−G̃ (`) t ‖2F] ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]−‖L‖2F|×E[‖∇L̃(θt−1)‖2F] ) (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 t ∈ {EK + 1, . . . , EK + 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(`+1)t but different input H̃(`−1)t ,H (`−1) t is defined as
  E[‖∇W f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇W f (`)(D (`+1) t ,H (`−1) t ,W (`) t )‖2F]
  ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (210)
  Proof. For the gradient w.r.t. the weight matrices, we have
  E[‖∇W f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇W f (`)(D (`+1) t ,H (`−1) t ,W (`) t )‖2F] = E[‖ ( G̃
  (`) t−1 + [L̃ (`)H̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t ))− [L̃(`)H̃ (`−1) t−1 ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1)) ) − ( [LH (`−1) t ] >(D (`+1) t ◦ σ′(Z (`) t )) ) ‖2F]
  ≤ 2E[‖ ( G̃
  (`) t−1 + [L̃ (`)H̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t ))− [L̃(`)H̃ (`−1) t ] >(D (`+1) t−1 ◦ σ′(Z̃ (`) t−1)) ) ︸ ︷︷ ︸
  (A) − ( [LH̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t )) ) ‖2F]︸ ︷︷ ︸
  (A) + 2E[‖ ( [LH̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t )) ) − ( [LH (`−1) t ] >(D (`+1) t ◦ σ′(Z (`) t )) ) ‖2F]︸ ︷︷ ︸
  (B)
  (211) Let first take a closer look at term (A). Let suppose t ∈ {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
  (A) ≤ EK+K∑ t=EK+1 ( O(‖W(1)t −W (1) t−1‖2F) + . . .+O(‖W (L) t −W (L) t−1‖2F) ) (212)
  Then we take a closer look at term (B). E[‖ ( [LH̃ (`−1) t ] >(D (`+1) t ◦ σ′(Z̃ (`) t )) ) − ( [LH (`−1) t ] >(D (`+1) t ◦ σ′(Z (`) t )) ) ‖2F]
  ≤ 2(B2LAB2DC4σ +B2LAB2HB2D)E[‖Z̃ (`) t − Z (`) t ‖2F]︸ ︷︷ ︸
  (C)
  (213)
  Let suppose t ∈ {EK + 1, . . . , EK + K}. Then we can denote t = EK + k for some k ≤ K. From Eq. 220, we know that
  (C) = E[‖LH̃(`−1)t W (`) t − Z̃ (`) t ‖2F]
  ≤ EK+K∑ t=EK+1 O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖W (`) t −W (`) t−1‖2F] ) (214) Plugging (C) back to (B), combing with (A) we have
  E[‖∇W f̃ (`)(D(`+1)t , H̃ (`−1) t ,W (`) t )−∇W f (`)(D (`+1) t ,H (`−1) t ,W (`) t )‖2F]
  ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt)‖2F] ) (215)
  In the following lemma, we provide the upper-bound of Eq. 153, which is one of the three key factors that affect the mean-square error of stochastic gradient at the `th layer.
  Lemma 11. Let suppose t ∈ {EK + 1, . . . , EK +K}. Then the upper bound of
  E[‖D̃(L+1)t −D (L+1) t ‖2F] ≤ EK+K∑ t=EK+1 η2 ×O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (216)
  Proof. By definition we have
  ‖D̃(L+1)t −D (L+1) t ‖2F ≤ ∥∥∥∂L̃(θt) ∂H̃
  (L) t
  − ∂L(θt) ∂H
  (L) t−1
  ∥∥∥2 F
  ≤ L2loss‖H̃ (L) t −H (L) t ‖2F
  ≤ L2lossC2σ ‖Z̃ (L) t − Z (L) t ‖2F︸ ︷︷ ︸
  (A)
  (217)
  Let take closer look at term (A):
  ‖Z̃(L)t − Z (L) t ‖2F = ‖f̃ (L)(H̃ (L−1) t ,W (L))− f (L)(H(L−1)t ,W(L))‖2F ≤ 2‖f̃ (L)(H̃(L−1)t ,W(L))− f (L)(H̃ (L−1) t ,W
  (L))‖2F + 2‖f (L)(H̃(L−1)t ,W(L))− f (L)(H (L−1) t ,W
  (L))‖2F = 2‖f̃ (L)(H̃(L−1)t ,W(L))− f (L)(H̃ (L−1) t ,W
  (L))‖2F + 2‖σ(LH̃(L−1)t W(L))− σ(LH (L−1) t W
  (L))‖2F ≤ 2‖f̃ (L)(H̃(L−1)t ,W(L))− f (L)(H̃ (L−1) t ,W
  (L))‖2F + 2C4σB 2 LAB 2 W ‖Z̃ (L−1) t − Z (L−1) t ‖2F
  (218)
  By induction, we have
  ‖Z̃(L)t − Z (L) t ‖2F ≤ L∑ `=1 O(‖f̃ (`)(H̃(`−1)t ,W(`))− f (`)(H̃ (`−1) t ,W (`))‖2F) (219)
  Let suppose t ∈ {EK + 1, . . . , EK + K}. Then we can denote t = EK + k for some k ≤ K. From Eq. 220, we know that
  ‖f̃ (`)(H̃(`−1)t ,W(`))− f (`)(H̃ (`−1) t ,W (`))‖2F
  ≤ EK+K∑ t=EK+1 O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖W (`) t −W (`) t−1‖2F] ) (220) Therefore, we know
  ‖D̃(L+1)t −D (L+1) t ‖2F ≤ L2lossC2σ‖Z̃ (L) t − Z (L) t ‖2F
  ≤ L∑ `=1 EK+K∑ t=EK+1 O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖W (`) t −W (`) t−1‖2F] )
  = EK+K∑ t=EK+1 η2 ×O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (221)
  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 + 1, . . . , EK +K}. Then we can denote t = EK + k for some k ≤ K such that
  E[‖∇L̃(θt)−∇L(θt)‖2F] ≤ EK+K∑ t=EK+1 η2O ( |E[‖L̃(`)‖2F]− ‖L‖2F| × E[‖∇L̃(θt−1)‖2F] ) (222)
  Proof. From Lemma 6, we have
  E[‖G̃(`) −G(`)‖F 2]
  ≤ O(E[‖D̃(L+1) −D(L+1)‖2F])
  +O(E[‖∇H f̃ (L)(D(L+1), H̃(L−1)W(L))−∇Hf (L)(D(L+1),H(L−1)W(L))‖2F]) + . . .
  +O(E[‖∇H f̃ (`+1)(D(`+2), H̃(`),W(`+1))−∇Hf (`+1)(D(`+2),H(`),W(`+1))‖2F])
  +O(E[‖∇W f̃ (`)(D(`+1), H̃(`−1),W(`))−∇W f (`)(D(`+1),H(`−1),W(`))‖2F]) (223)
  Plugging the result from support lemmas, i.e., Lemma 8, 10, 11 and using the definition of stochastic gradient for all model parameters∇L̃(θt) = {G̃(`)t }L`=1 we have
  E[‖∇L̃(θt)−∇L(θt)‖2F] ≤ EK+K∑ t=EK+1 η2O ( L∑ `=1 |E[‖L̃(`)‖2F]− ‖L‖2F| ×E[‖∇L̃(θt−1)‖2F] ) (224)
  I.2 REMAINING STEPS TOWARD THEOREM 3
  By the smoothness of L(θt), we have
  L(θT+1) ≤ L(θt) + 〈∇L(θt),θt+1 − θt〉+ Lf 2 ‖θt+1 − θt‖2F
  = L(θt)− η〈∇L(θt),∇L̃(θt)〉+ η2Lf
  2 ‖∇L̃(θt)‖2F
  = (a) L(θt)−
  η 2 ‖∇L(θt)‖2F + η 2 ‖∇L(θt)−∇L̃(θt)‖2F − (η 2 − Lfη 2 2 ) ‖∇L̃(θt)‖2F
  (225) where equality (a) is due to the fact 2〈x,y〉 = ‖x‖2F + ‖y‖2F − ‖x− y‖2F for any x,y. Take expectation on both sides, we have
  E[L(θT+1)] ≤ E[L(θt)]− η
  2 E[‖∇L(θt)‖2F]+
  η 2 E[‖∇L(θt)−∇L̃(θt)‖2F]− (η 2 −Lfη 2 2 ) E[‖∇L̃(θt)‖2F]
  (226)
  By summing over t = 1, . . . , T where T is the inner-loop size, we have T∑ t=1 E[‖∇L(θt)‖2F]
  ≤ 2 η
  ( E[L(θ1)]− E[L(θT+1)] ) + T∑ t=1 [E[‖∇L(θt)−∇L̃(θt)‖2F]− ( 1− Lfη ) E[‖∇L̃(θt)‖2F]]
  ≤ 2 η
  ( E[L(θ1)]− E[L(θ?)] ) + T∑ t=1 [E[‖∇L(θt)−∇L̃(θt)‖2F]− ( 1− Lfη ) E[‖∇L̃(θt)‖2F]]
  (227) where L(θ?) is the global optimal solution.
  Let us consider each inner-loop with E ∈ [0, T/K − 1] and t ∈ {EK + 1, . . . , EK + K}. Using Lemma 12 we have
  T/K−1∑ E=0 EK+K∑ t=EK+1 E[‖L(θt)−∇L̃(θt)‖2F]− ( 1− Lfη ) T/K−1∑ E=0 EK+K∑ t=EK+1 E[‖∇L̃(θt)‖2F]
  ≤ O ( L∑ `=1 |E[‖L̃(`)‖2F]− ‖L‖2F| )( T/K−1∑ E=0 EK+K∑ t=EK+1 t−EK∑ j=2 η2E[‖∇L̃(θ)j−1‖2F] )
  − ( 1− Lfη ) T/K−1∑
  E=0 EK+K∑ t=EK+1 E[‖∇L̃(θt)‖2F]
  ≤ [η2K ×O ( L∑ `=1 |E[‖L̃(`)‖2F]− ‖L‖2F| ) − ( 1− Lfη ) ] T/K−1∑ E=0 EK+K∑ t=EK+1 E[‖∇L̃(θt)‖2F] = [η2∆++′b+n − (1− Lfη)] T/K−1∑ E=0 EK+K∑ t=EK+1 E[‖∇L̃(θt)‖2F]
  (228)
  Notice that η = 2 Lf+ √ L2f+4∆ ++′ b+n is a root of equation η2∆++′b+n − (1 − Lfη) = 0. Therefore we have T∑ t=1 E[‖∇L(θt)‖2F] ≤ 2 η ( E[L(θ1)]− E[L(θ?)] ) (229)
  which implies
  1
  T T∑ t=1 E[‖∇L(θt)‖2F] ≤ 1 T ( Lf + √ L2f + 4∆ ++′ b+n )( E[L(θ1)]− E[L(θ?)] ) (230)
  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 a jupyter notebook file for a 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 to find that the overall trend of loss curves and conclusions will remains the same.
  This implementation is based on 4PyTorch 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 requirements.txt
  Experiments are produced on PPI, PPI-Large, Flickr, Reddit, and Yelp datasets. The utilized datasets can be downloaded from 5Google drive.
  # create folders that save experiment results and datasets $ mkdir ./results $ mkdir ./data # please download the dataset and put them inside this folder
  4https://pytorch.org/ 5https://drive.google.com/drive/folders/15eP7OHiHQUnDrHKYh1YPxXkiqGoJhbis?
  usp=sharing
  To reproduce the results, please run the following commands:
  $ python train.py --sample_method 'ladies' --dataset 'reddit' $ python train.py --sample_method 'fastgcn' --dataset 'reddit' $ python train.py --sample_method 'graphsage' --dataset 'reddit' $ python train.py --sample_method 'vrgcn' --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' $ python train.py --sample_method 'fastgcn' --dataset 'ppi' $ python train.py --sample_method 'graphsage' --dataset 'ppi' $ python train.py --sample_method 'vrgcn' --dataset 'ppi' $ python train.py --sample_method 'graphsaint' --dataset 'ppi' $ python train.py --sample_method 'exact' --dataset 'ppi'
  $ python train.py --sample_method 'ladies' --dataset 'flickr' $ python train.py --sample_method 'fastgcn' --dataset 'flickr' $ python train.py --sample_method 'graphsage' --dataset 'flickr' $ python train.py --sample_method 'vrgcn' --dataset 'flickr' $ python train.py --sample_method 'graphsaint' --dataset 'flickr' $ python train.py --sample_method '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_method 'graphsage' --dataset 'ppi-large' $ python train.py --sample_method 'vrgcn' --dataset 'ppi-large' $ python train.py --sample_method 'graphsaint' --dataset 'ppi-large' $ python train.py --sample_method 'exact' --dataset 'ppi-large'
  $ python train.py --sample_method 'ladies' --dataset 'yelp' $ python train.py --sample_method 'fastgcn' --dataset 'yelp' $ python train.py --sample_method 'graphsage' --dataset 'yelp' $ python train.py --sample_method 'vrgcn' --dataset 'yelp' $ python train.py --sample_method 'graphsaint' --dataset 'yelp' $ python train.py --sample_method 'exact' --dataset 'yelp'