Year: 2021

Abstract:
Graph Convolution Networks (GCN) are used in numerous settings involving a large underlying graph as well as several layers. Standard SGD-based training scales poorly here since each descent step ends up updating node embeddings for a large portion of the graph. Recent methods attempt to remedy this by subsampling the graph which does reduce the compute load, but at the cost of biased gradients which may offer suboptimal performance. In this work we introduce a new method IGLU that caches intermediate computations at various GCN layers. This enables IGLU to perform lazy updates that do not require updating a large number of node embeddings during descent and offers much faster convergence without significantly biasing the gradients. Under standard assumptions such as objective smoothness, IGLU provably converges to a first-order saddle point. We validate IGLU extensively on a variety of benchmarks, where it offers up to 1.2% better accuracy despite requiring up to 88% less wall-clock time.

heading: 1 INTRODUCTION
text: The Graph Convolution Network (GCN) model has received much attention as an effective graph representation learning technique. It can exploit network topology while embedding data points enabling superior performance in several applications such as node classification on graphs (Kipf & Welling, 2017), recommendation systems (Ying et al., 2020) and program repair (Yasunaga & Liang, 2020). Their success notwithstanding, training GCNs at scale remains challenging especially on large and dense graphs and multiple convolution layers. This is mainly due to the aggregation operation that enables GCNs to adapt to graph topology – a node’s output layer embedding in a GCN is influenced by embeddings of its neighbors in the previous layer which recursively depend on embeddings of their own neighbors in the still previous layer, and so on.
Even in GCNs with 2-3 layers, when processing loss terms corresponding to a small sample of nodes in a mini-batch, this causes back propagation to update a large multi-hop neighborhood, causing mini-batch SGD-based techniques to scale poorly especially in dense graphs. Efforts to overcome this problem try to limit the number of nodes that receive updates as a result of a back-propagation step Chiang et al. (2019); Hamilton et al. (2018); Zeng et al. (2020). This is done either by sub-sampling the neighborhood or clustering (note the distinction between nodes sampled to create a mini-batch and neighborhood sampling done to limit the neighborhood of the mini-batch that receives updates). Variance reduction techniques have also been studied Chen et al. (2018a) to reduce additional variance introduced by neighborhood sampling. However, these techniques often require large graphs to be heavily subsampled resulting in poor accuracy due to insufficient aggregation. They also do not guarantee unbiased gradients or rigorous convergence guarantees. See Section 2 for a more detailed discussion about state-of-the-art methods for GCN training.
Our Contributions: This paper presents IGLU, an efficient technique for training GCNs based on lazy updates. An analysis of the gradient structure in GCNs reveals the most expensive component of the back-propagation step initiated at a node to be (re-)computation of forward-pass embeddings for its vast multi-hop neighborhood. Based on this observation, IGLU performs back-propagation with significantly reduced complexity using intermediate computations that are cached at regular intervals. This completely avoids neighborhood sampling and is a stark departure from the state-of-the-art. IGLU is architecture-agnostic and can be readily implemented on a wide range of GCN architectures. Avoiding neighborhood sampling also allows IGLU to completely avoid variance artifacts and offer
provable convergence to a first-order stationary point under standard assumptions. In experiments, IGLU offered superior accuracies and accelerated convergence on a range of benchmark datasets.

heading: 2 RELATED WORKS
text: (Bruna et al., 2013; Defferrard et al., 2016; Kipf & Welling, 2017) introduced the GCN architecture for transductive learning on graphs. Later works extended to inductive settings and explored architectural variants such as the GIN (Xu et al., 2019). Much effort has focused on speeding-up GCN training.
Sampling Based Approaches: The simplest neighborhood sampling strategy adopted by GraphSAGE (Hamilton et al., 2018) iteratively sub-samples the multi-hop neighborhood of a node before proceeding to initiate back-propagation at that node. Only the sub-sampled neighbors participate in the back-propagation step limiting the amount of computation. Layer sampling strategies such as FastGCN (Chen et al., 2018b), LADIES (Zou et al., 2019) and ASGCN (Huang et al., 2018) decrease computation by sampling nodes at each GCN layer, using importance sampling to reduce variance and improve connectivity among sampled nodes. FastGCN used the same sampling distribution for all layers and struggled to maintain connectivity unless large batch-sizes were used. LADIES used a per-layer distribution conditioned on nodes sampled for the layer above allowing more efficient use of samples under a budget. ASGCN used a linear model to jointly infer node importance weights. Recent works such as Cluster-GCN (Chiang et al., 2019) and GraphSAINT (Zeng et al., 2020) propose subgraph sampling creating mini-batches out of subgraphs and restricting back-propagation to nodes within the subgraph. To avoid losing too many edges, large mini-batch sizes are used. Cluster-GCN performs graph clustering and chooses multiple clusters per mini-batch while reinserting cross-cluster edges whereas GraphSAINT samples large subgraphs directly using random walks.
Bias and Variance: Sampling techniques face bias as non-linear activations in the GCN architecture make it difficult to offer unbiased estimates of the loss function e.g. Zeng et al. (2020) offer un-biased estimates only if non-linearities are discarded. Sampling techniques also face increased variance for which variance-reduction techniques have been proposed such as VR-GCN (Chen et al., 2018a), MVS-GNN (Cong et al., 2020) and AS-GCN (Huang et al., 2018). VR-GCN samples nodes whose embeddings are to be updated at each layer and uses stale embeddings for the rest and offers variance elimination in the limit under suitable conditions. MVS-GNN handles variance due to mini-batch creation by performing importance weighted sampling to construct mini-batches. Bandit Sampler (Liu et al., 2020) formulates variance reduction as an adversarial bandit problem.
Other Approaches: Recent approaches decouple propagation from prediction as a pre-processing step e.g. PPRGo (Bojchevski et al., 2020), APPNP (Klicpera et al., 2018) and SIGN (Frasca et al., 2020). APPNP makes use of the relationship between the GCNs and PageRank to construct improved propagation schemes via personalized PageRank. PPRGo extends APPNP by approximating the dense propagation matrix via the push-flow algorithm. SIGN proposes inception style pre-computation of graph convolutional filters to speed up training and inference. GNNAutoScale (Fey et al., 2021) builds on VR-GCN and makes use of historical embeddings for scaling GNN training to large graphs.
IGLU in Context of Related Work: IGLU avoids neighborhood sampling entirely and instead speeds-up learning using stale computations. Intermediate computations are cached and lazily updated at regular intervals e.g. once per epoch. We note that IGLU’s caching is distinct and much more aggressive (e.g. lasting an entire epoch) than the internal caching performed by popular frameworks such as TensorFlow and PyTorch (where caches last only a single iteration). Recomputing these values in bulk offers IGLU economies of scale. IGLU faces zero sampling variance but incurs bias due to the use of stale computations. Fortunately, this bias is provably bounded, and can be made arbitrarily small by adjusting the step length and refresh frequency of the stale computations.

heading: 3 IGLU: EFFICIENT GCN TRAINING VIA LAZY UPDATES
text: Problem Statement: Consider the problem of learning a GCN architecture on an undirected graph G(V, E) with each of the N nodes endowed with an initial feature vector x0i ∈ Rd0 , i ∈ V . X0 ∈ Rn×d0 denotes the matrix of these initial features stacked together. N (i) ⊂ V denotes the set of neighbors of node i. A denotes the (normalized) adjacency matrix of the graph. A multi-layer GCN architecture uses a parameterized function at each layer to construct a node’s embedding for the next
layer using embeddings of that node as well as those of its neighbors. Specifically
xki = f(x k−1 j , j ∈ {i} ∪ N (i);E k),
where Ek denotes the parameters of k-th layer. For example, a standard GCN layer is given by,
xki = σ ∑ j∈V Aij(W k)>xk−1j  , where Ek is simply the matrix W k ∈ Rdk−1×dk and dk is the embedding dimensionality after the kth layer. More involved architectures exist that incorporate operations such as layer normalization, batch normalization, etc. In this paper, Ek will always denote the collection of all parameters of the kth layer e.g. offset and scale parameters used in a layer norm operation. Xk ∈ Rn×dk will denote the matrix of kth layer embeddings stacked together, giving us the handy shorthand Xk = f(Xk−1;Ek). Given aK-layer GCN and a multi-label/multi-class task with C labels/classes, a fully-connected layer WK+1 ∈ RdK×C and appropriate activation function such as softmax is applied to get predictions which are then fed into the task loss. We note that IGLU does not require the task loss to decompose over the classes. The convergence proofs only require a smooth training objective function.
Neighborhood Explosion: To understand the reason behind neighborhood explosion and cost of mini-batch based SGD training, consider a univariate regression problem with a no-frills 2-layer GCN with unidimensional features and sigmoidal activation within the hidden layers i.e. K = 2 and C = 1 = d0 = . . . = dk. This GCN is parameterized by w1, w2, w3 ∈ R and offers the output ŷi = w3σ ( z2i ) where z2i = ∑ j∈V Aijw 2x1j ∈ R. In turn, we have x1j = σ ( z1i )
where z1i = ∑ j′∈V Ajj′w
1x0j′ ∈ R and x0j′ ∈ R are the initial features of the nodes. Given a task loss ` : RN × RN e.g. least squares, denoting `′i = `′(ŷi, yi) gives us
∂`(ŷi, yi)
∂w1 = `′i · ∂ŷi ∂z2i · ∂z
2 i
∂w1 = `′i · w3σ′(z2i ) · ∑ j∈V Aijw 2 ∂x1j ∂w1
= `′i · w3σ′(z2i ) · ∑ j∈V Aijw 2σ′(z1j ) · ∑ j′∈V Ajj′x 0 j′ .
The nesting of the summations is conspicuous and establishes the neighborhood explosion: when seeking gradients in a K-layer GCN on a graph with average degree m, up to an mK−1-sized neighborhood of a node may be involved in the back-propagation update initiated at that node. Note that the above expression involves terms such as σ′(z2i ), σ ′(z1j ) and that the values of z 2 i , z 1 j etc change whenever the model receives updates. Consequently, for a fresh mini-batch of nodes chosen for a particular iteration, terms such as σ′(z2i ), σ
′(z1j ) need to be computed afresh if the gradient is to be computed exactly. Performing these computations amounts to doing forward pass operations that frequently involve a large neighborhood of the nodes of the mini-batch. Sampling strategies try to limit this cost by directly restricting the neighborhood over which such forward passes are computed. However, this introduces both bias and variance into the gradient updates as discussed in Section 2. IGLU instead lazily updates various incomplete gradient (defined below) and node embedding terms that participate in the above expression. This completely eliminates sampling variance but introduce a bias due to the use of stale terms. However, this bias is not just provably bounded, but can be made arbitrarily small by adjusting the step length and frequency of refreshing these terms.
Lazy Updates for GCN Training: Consider an arbitrary GCN architecture with the following structure: for some parameterized layer functions we have Xk = f(Xk−1;Ek) where Ek denotes the collection of all parameters of the kth layer e.g. weight matrices, offset and scale parameters used in layer norm operations, etc. Xk ∈ RN×dk denotes the matrix of kth layer embeddings stacked together and X0 ∈ RN×d0 are the initial features. For a K-layer GCN on a multi-label/multi-class task with C labels/classes, a fully-connected layer WK+1 ∈ RdK×C is additionally included to offer ŷi = (W
K+1)>xKi ∈ RC . Let us use the shorthand Ŷ ∈ RN×C to denote the matrix where the predicted outputs ŷi over all the nodes are stacked. We assume a task loss function ` : RC×RC×R+ that need not decompose over the classes. We will use the abbreviation `i := `(ŷi, yi). For sake of simplicity, we assume that the loss function itself includes any activation such as softmax that needs to be applied over the predictions ŷi. Let L = ∑ i∈V `i denote the total loss function.
Motivation: We define the the loss derivative matrix G = [gic] = RN×C as gic := ∂`i∂ŷic . As the proof of Lemma 1 (See Appendix C) shows, the loss derivative with respect to parameters Ek at
any layer has the form ∂L ∂Ek = ∑N j=1 ∑dk p=1 (∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xkjp ) ∂Xkjp ∂Ek . Note that the loss derivative is expressed in terms of partial derivatives ∂Xkjp ∂Ek
that can be computed for any node using only embeddings of its neighbors in the (k − 1)th layer i.e. Xk−1 thus avoiding any neighborhood explosion. This means that neighborhood explosion must be encountered while computing the terms encapsulated in the round brackets. Let us first formally recognize these terms as incomplete gradients. The notation ∂P∂Q ∣∣∣ R
denotes the partial derivative of P w.r.t Q while keeping R fixed i.e. treated as a constant.
Definition 1. For any layer k ≤ K, define its incomplete task gradient to be αk = [αkjp] ∈ RN×dk ,
αkjp := ∂(G Ŷ) ∂Xkjp ∣∣∣∣∣ G = ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xkjp
The following lemma completely characterizes the loss gradients of the GCN but also shows that the incomplete gradient terms αk, k ∈ [K] can be efficiently computed using a recursive formulation that also does not involve any neighborhood explosion.
Lemma 1. The following results hold whenever the task loss L is differentiable:
1. For the final fully-connected layer we have ∂L ∂WK+1 = (XK)>G as well as for any k ∈ [K]
and any parameter Ek in the kth layer, ∂L ∂Ek = ∂(α k Xk) ∂Ek ∣∣∣ αk = ∑ i∈V ∑dk p=1 α k ip · ∂Xkip ∂Ek .
2. For the final layer, we have αK = G(WK+1)> as well as for any k < K, we have
αk = ∂(α k+1 Xk+1) ∂Xk ∣∣∣ αk+1 i.e. αkjp = ∑ i∈V ∑dk+1 q=1 α k+1 iq · ∂Xk+1iq ∂Xkjp .
Lemma 1 establishes a recursive definition of the incomplete gradients using terms such as ∂Xk+1iq ∂Xkjp that only concern a single layer. Thus, computing αk for any k ∈ [K] does not involve any neighborhood explosion since only the immediate neighbors of a node need be consulted. Lemma 1 also shows that if αk are computed and frozen, the loss derivatives ∂L
∂Ek only involve additional computation
of terms such as ∂Xkip ∂Ek
which yet again involve a single layer and do not cause neighborhood explosion. This motivates lazy updates to αk, Xk values in order to accelerate back-propagation. However, performing lazy updates to both αk, Xk offers suboptimal performance. Hence IGLU
adopts two variants described in Algorithms 1 and 2. The backprop variant* keeps embeddings Xk stale for an entire epoch but performs eager updates to αk. The inverted variant on the other hand keeps the incomplete gradients αk stale for an entire epoch but performs eager updates to Xk. Algorithm 1 IGLU: backprop order Input: GCN G, initial features X0, task loss L 1: Initialize model parameters Ek, k ∈ [K],WK+1 2: while not converged do 3: Do a forward pass to compute Xk for all k ∈
[K] as well as Ŷ 4: Compute G then ∂L
∂WK+1 using Lemma 1 (1)
and update WK+1 ←WK+1 − η · ∂L ∂WK+1
5: Compute αK using G,WK+1, Lemma 1 (2) 6: for k = K . . . 2 do 7: Compute ∂L
∂Ek using αk, Xk, Lemma 1 (1)
8: Update Ek ← Ek − η · ∂L ∂Ek
9: Update αk using αk+1 using Lemma 1 (2) 10: end for 11: end while
Algorithm 2 IGLU: inverted order Input: GCN G, initial features X0, task loss L 1: Initialize model parameters Ek, k ∈ [K],WK+1 2: Do an initial forward pass to computeXk, k ∈ [K] 3: while not converged do 4: Compute Ŷ,G and αk for all k ∈ [K] using
Lemma 1 (2) 5: for k = 1 . . .K do 6: Compute ∂L
∂Ek using αk, Xk, Lemma 1 (1)
7: Update Ek ← Ek − η · ∂L ∂Ek 8: Update Xk ← f(Xk−1;Ek) 9: end for
10: Compute ∂L ∂WK+1
using Lemma 1 (1) and use it to update WK+1 ←WK+1 − η · ∂L
∂WK+1
11: end while
SGD Implementation: Update steps in the algorithms (steps 4, 8 in Algorithm 1 and steps 7, 10 in Algorithm 2) are described as a single gradient step over the entire graph to simplify exposition – in practice these steps are implemented using mini-batch SGD. This is done as usual by sampling a mini-batch of nodes S and computing task gradients only w.r.t L̂S = ∑ i∈S `i instead of L.
Contribution: As noted in Section 2, IGLU uses caching in a manner fundamentally different from popular frameworks such as PyTorch or TensorFlow. Their caches are short-lived and always seek to compute exact gradients unlike IGLU that computes gradients much faster but with bounded bias. IGLU also extends other caching-based techniques such as VR-GCN which choose to cache node embeddings. In contrast, IGLU offers two variants and the variant that uses inverted order of updates (Algorithm 2) and caches incomplete gradients outperforms the backprop variant of IGLU (Algorithm 1) that caches node embeddings instead.
Theoretical Analysis: Note that, conditioned on the stale parameters (either αk or Xk depending on which variant is being executed), the gradients used by IGLU to perform model updates (steps 4, 8 in Algorithm 1 and steps 7, 10 in Algorithm 2) do not exhibit any sampling bias. However, the staleness itself introduces bias in the training process. By controlling the step length η and the frequency with which the stale parameters are updated, this bias can be provably controlled resulting in guaranteed convergence to a first-order stationary point. Due to lack of space, we postpone the detailed statement and proof of the convergence guarantee to Appendix C. Theorem 2 (IGLU Convergence (Informal)). Suppose the task loss function L has O (1)-Lipschitz gradients and IGLU is executed with small enough step lengths η with model updates (steps 4, 8 in Algorithm 1 and steps 7, 10 in Algorithm 2) being carried out either in a full-batch or mini-batch SGD manner, then within T iterations, IGLU converges to a model iterate satisfying:
1. ‖∇L‖22 ≤ O ( 1/T 2 3 ) if update steps are carried out on the entire graph in a full-batch.
2. ‖∇L‖22 ≤ O ( 1/ √ T ) if update steps are carried out using mini-batch SGD.
This result holds under minimal assumptions of objective smoothness and boundedness that are standard Chen et al. (2018a); Cong et al. (2020), yet offers convergence rates comparable to those offered by standard mini-batch SGD. However, whereas works such as (Chen et al., 2018a) assume bounds on the sup-norm i.e. L∞ norm of the gradients, Theorem 2 only requires an L2 norm bound. Note that objective smoothness requires the architecture to use smooth activation functions. However, IGLU offers similar performance whether using non-smooth activations e.g. ReLU or smooth ones e.g. GELU (see Appendix B.7) as is also observed by other works (Hendrycks & Gimpel, 2020).
*The backprop variant is named so since it updates model parameters in the order back-propagation would have updated them i.e. WK+1 followed by EK , EK−1, . . . whereas the inverted variant performs updates in the reverse order i.e. starting from E1, E2 all the way to WK+1.

heading: 4 EMPIRICAL EVALUATION
text: IGLU was compared to state-of-the-art (SOTA) baselines on several node classification benchmarks across multiple dimensions associated with GCN training such as accuracy and convergence rate. The inverted order of updates was used for IGLU that is justified empirically.
Datasets: The following five benchmark datasets were used: predicting the communities to which different posts belong in Reddit, classifying protein functions across various biological proteinprotein interaction graphs in PPI-Large (Hamilton et al., 2018), categorizing types of images based on descriptions and common properties in Flickr (Zeng et al., 2020), predicting paper-paper associations in OGBN-Arxiv (Hu et al., 2021) and categorizing meaningful associations between proteins in OGBN-Proteins (Hu et al., 2021). Training, validation and test splits were used for each dataset in a manner consistent with the original release of the dataset as were metrics: specifically ROC AUC for OGBN-Proteins and Micro-F1 for all the other datasets. Dataset descriptions and statistics are presented in the Appendix. The graphs considered are varied in size (56K nodes in PPI-Large to 232K nodes in Reddit), density (average degree of 13 in OGBN-Arxiv to 597 in OGBN-Proteins) and number of edges (800K to 39 Million). They are fairly diverse in the nature of information they capture and the type of node classification tasks they solve, covering both multi-label and multi-class problems, thus enabling an extensive evaluation of IGLU.
Baselines: IGLU was compared to state-of-the-art algorithms, namely - GCN (Kipf & Welling, 2017), GraphSAGE (Hamilton et al., 2018), VR-GCN (Chen et al., 2018a), Cluster-GCN (Chiang et al., 2019) and GraphSAINT (Zeng et al., 2020) (using the Random Walk Sampler which was reported by the authors to have the best performance). The mini-batched implementation of GCN provided by GraphSAGE authors was used since the version released by (Kipf & Welling, 2017) ran into run time error on all datasets. GraphSAGE and VRGCN alleviate the neighborhood explosion problem by sampling subsets of neighbors. ClusterGCN and GraphSAINT are subgraph sampling techniques using sampled subgraphs for constructing convolution layers. Thus, our baselines include neighbor sampling, layer-wise sampling, subgraph sampling and no-sampling methods. IGLU does not require any sampling and hence does not have additional sampling overheads. IGLU was implemented in TensorFlow and compared with the authors’ TensorFlow implementations of the baselines. Due to lack of space, comparisons with the following additional baselines is provided in the Appendix B.2: LADIES (Zou et al., 2019), L2-GCN (You et al., 2020), AS-GCN (Huang et al., 2018), MVS-GNN (Cong et al., 2020), FastGCN (Chen et al., 2018b), SIGN (Frasca et al., 2020), PPRGo (Bojchevski et al., 2020) and Bandit Sampler (Liu et al., 2020).
Architecture: We note that all these baselines propose a specific network architecture along with their proposed sampling strategy. We find that these architectures augment the standard GCN architecture (Kipf & Welling, 2017) e.g. multiple intermediate non-linear layers, normalization layers and concatenation layers, all of which can help improve performance. As IGLU is architecture-agnostic, it can be readily used with any such architectures. For the experiments, the network architectures proposed by VR-GCN and GraphSAINT was used with IGLU owing to their consistent performance across datasets as demonstrated by (Chen et al., 2018a; Zeng et al., 2020). Both architectures were considered and results for the best performing architecture are reported for each dataset.
Detailed Setup. The inductive supervised learning setting was considered for all five datasets, motivated by the fact that this setting is more general than the transductive setting which assumes availability of the entire graph during training. This setting also generalizes to the real-world where graphs can grow over time and information for all the nodes may not available (Hamilton et al., 2018). Experiments on PPI-Large, Reddit and Flickr used 2 Layer GCNs for both IGLU and baselines whereas OGBN-Proteins and OGBN-Arxiv used 3 Layer GCNs for all methods as prescribed in the original benchmark (Hu et al., 2021). Further implementation details are provided in Appendix A.1.
Model Selection and Hyperparameter Tuning. Model selection was done for all methods based on their validation set performance. Each experiment was run five times and mean and standard deviation is reported on the test performance in Table 1 along with training times. Although the embedding dimension varies across datasets, they are same for all methods for any dataset. For IGLU, GraphSAGE and VR-GCN an exhaustive grid search was done over general hyperparameters such as batch size, learning rate and dropout rate (Srivastava et al., 2014). In addition method-specific hyperparameter sweeps were also carried out that are detailed in Appendix A.4 due to lack of space.

heading: 4.1 RESULTS
text: For IGLU, the VR-GCN architecture performed better for Reddit and OGBN-Arxiv datasets while the GraphSAINT architecture performed better on the remaining three datasets. All baselines were extensively tuned on the 5 datasets and their reported performance was either replicated closely or else improved in some cases. Test accuracies are reported in Table 1 and convergence plots are shown in Figure 2. Additionally, Table 1 also reports the absolute accuracy gain of IGLU over the best baseline. The % speedup offered by IGLU is computed as follows: let the highest validation score obtained by the best baseline be v1 and t1 be the time taken to reach that score. Let the time taken by IGLU to reach v1 be t2. Then,
% Speedup := t1 − t2 t1 × 100 (1)
The wall clock training time in Figure 2 is strictly the optimization time for each method and excludes method-specific overheads such as pre-processing, sampling and sub-graph creation that other baselines incur. This is actually a disadvantage to IGLU since its overheads are much smaller. The various time and memory overheads incurred by all methods are summarized in Appendix A.2 and memory requirements for IGLU in Appendix A.3.
Performance on Test Set and Speedup Obtained: It is evident from Table 1 that IGLU significantly outperforms the baselines on PPI-Large, Flickr, OGBN-Proteins and OGBN-Arxiv and is competitive with best baselines on Reddit. On PPI-Large, IGLU improves accuracy upon the best baseline (VRGCN) by over 1.2% while providing speedups of upto 88%, i.e., IGLU is about 8× faster to train than VR-GCN. On OGBN-Proteins, IGLU improves accuracy upon the best baseline (GraphSAINT) by over 2.6% while providing speedup of 11%. On Flickr, IGLU offers 0.98% improvement in accuracy while simultaneously offering upto 45% speedup over the previous state-of-the-art method GraphSAINT. Similarly on Arxiv, IGLU provides 0.84% accuracy improvement over the best baseline GraphSAINT while offering nearly 14% speedup. On Reddit, an 8.1% speedup was observed in convergence to a high validation score of 0.955 while the final performance is within a standard deviation of the best baseline. The OGBN-Proteins and OGBN-Arxiv datasets are originally benchmarked in the transductive setting, with the entire graph information made available during training. However, we consider the more challenging inductive task yet IGLU outperforms the best transductive baseline by over 0.7% for Proteins while matching the best transductive baseline for Arxiv (Hu et al., 2021). It is important to note that OGBN-Proteins is an atypical dataset because the graph is highly dense. Because of this, baselines such as ClusterGCN and GraphSAINT that drop a lot of edges while creating subgraphs show a deterioration in performance. ClusterGCN runs into a runtime error on this dataset (denoted by || in the table), while GraphSAINT requires a very large subgraph size to achieve reasonable performance.
Convergence Analysis: Figure 2 shows that IGLU converges faster to a higher validation score than other baselines. For PPI-Large, while Cluster-GCN and VR-GCN show promising convergence in the initial stages of training, they stagnate at a much lower validation score in the end whereas IGLU is able to improve consistently and converge to a much higher validation score. For Reddit, IGLU’s final validation score is only marginally lower than GraphSAINT but IGLU offers rapid convergence in the early stages of training. For OGBN-Proteins, Flickr and OGBN-Arxiv, IGLU demonstrates a substantial improvement in both convergence and the final performance on the test set.

heading: 4.2 ABLATION STUDIES
text: Effect of the Order of Updates. IGLU offers the flexibility of using either the backprop order of updates or the inverted order of updates as mentioned in Section 3. Ablations were performed on the Reddit and OGBN-Proteins datasets to analyse the effect of the different orders of updates. Figure 3 offers epoch-wise convergence for the same and shows that the inverted order of updates provides more stability to IGLU as well as improved convergence in the early stages of training. We also observe that both of the orders eventually converge to similar validation scores. These ablations highlight the fact that the inverted order of updates accelerates convergence and justifies our choice of the inverted order of updates for IGLU for all ablation experiments hereon.
Analysis of Degrees of Staleness. The effect of frequency of updates to the incomplete gradients (αk) on the performance and convergence of IGLU was also analyzed. The ablation was conducted keeping all the other hyperparameters fixed. As the default setting all the αk’s were updated only once every epoch and referred to as update with frequency 1. Apart from this scheme, two other
update schemes are also considered: a) update the αk’s every two epochs, and b) update the αk’s twice within an epoch. We refer to these settings as update with frequency 0.5 and frequency 2 respectively. To summarize, the αk’s are the most fresh with update frequency 2 and most stale update frequency 0.5. This ablation study was performed on the PPI-Large dataset and each variant was trained for 200 epochs. Table 2 summarizes the results doing model selection and Figure 4 plots the convergence of these variants. Figure 4b shows that on PPI-Large, the default update frequency 1 has the best convergence on the validation set, followed by update frequency 0.5. Both update frequency 1 and 0.5 massively outperform update frequency 2. Figure 4a shows that IGLU with update frequency 2 has the lowest training loss but poor validation performance, indicating overfitting to the training dataset. We observe from this ablation that updating embeddings prematurely can cause unstable training resulting in convergence to a suboptimal solution. However, not refreshing the αk’s frequently enough can cause the model to take longer to converge to a good solution.
Additional Ablations and Experiments: Due to lack of space, we provide the following additional ablations and experiments in the Appendix: a) Ablation on degrees of staleness for the backprop order of updates in Appendix B.8, b) Experiments demonstrating IGLU’s scalability to deeper networks and larger datasets in Appendices B.3 and B.4 respectively, and c) Experiments demonstrating IGLU’s applicability and architecture-agnostic nature in Appendices B.5 and B.6.

heading: 5 DISCUSSION AND FUTURE WORK
text: We introduced IGLU, a novel method for training GCN architectures that uses biased updates based on cached intermediate computations to speed up gradient descent updates significantly. This bias is shown to be provably bounded so overall convergence is still effective (see Theorem 2). We validated IGLU’s performance on several datasets where it significantly outperformed SOTA methods in terms of accuracy and convergence speed. Ablation studies confirmed that IGLU is robust to its few hyperparameters to enable a near-optimal choice. Exploring other possible variants of IGLU, in particular reducing variance due to mini-batch SGD, sampling nodes to further speed-up updates, and exploring alternate staleness schedules are interesting future directions. On a theoretical side, it would be interesting to characterize properties of datasets and loss functions that influence the effect of lazy updates on convergence. Having such a property would allow practitioners to decide whether to execute IGLU with lazier updates or else reduce the amount of staleness. Exploring applicationand architecture-specific variants of IGLU is also an interesting direction.

heading: 6 REPRODUCIBILITY STATEMENT
text: We have made efforts to ensure that the results reported in this paper and the submitted material are reproducible and the methodology is easy to follow.
Theoretical Clarity: We discuss the problem setup and preliminaries in detail in Section 3, provide a detailed algorithm for our proposed methodology in Section 4. The discussion seeks to promote a sound theoretical understanding of the problem and our solution. Detailed proofs for all lemmata and theorems are provided in the appendices due to lack of space.
Experimental Reproducibility: Section 5 of the main paper and Appendix B contain information needed to reproduce the empirical results, namely: Datasets, their pre-processing and their source; implementation details for our method and the baselines, specifics regarding the architectures, hyperparameter search spaces and the best hyperparameters corresponding to the results in the paper.
Code Release: We will release the code for IGLU in the discussion forum to the reviewers and AC’s during the discussion period. Code for IGLU will also be released publicly upon acceptance of the paper in the spirit of reproducible research.

heading: 7 ETHICS STATEMENT
text: Our work presents a novel and effective mechanism to train GCN architectures on large graphs that outperforms the state of the art in terms of prediction accuracy and convergence speed. We do not consider any sensitive applications in this paper and our experiments focus primarily on publicly available benchmarks of scientific (e.g. PPI) and bibliographic (e.g. ArXiv) nature do not involve any user studies or human experimentation.

heading: APPENDIX
text: This appendix is segmented into three key parts.
1. Section A discusses additional implementation details. In particular, method-specific overheads are discussed in detail and detailed hyper-parameter settings for IGLU and the main baselines reported in Table 1 are provided. A key outcome of this analysis is that IGLU has the least overheads as compared to the other methods. We also provide the methodology for incorporating architectural modifications into IGLU, provide a detailed comparison with caching based methods and provide a more detailed descriptions of the algorithms mentioned in the main paper.
2. Section B reports dataset statistics, provides comparisons with additional baselines (not included in the main paper due to lack of space) and provides experiments on scaling to deeper models and larger graphs as mentioned in Section 4.2 in the main paper. It also provides experiments for the applicability of IGLU across settings and architectures, using smooth activation functions, continued ablation on degrees of staleness and optimization on the train set. The key outcome of this discussion is that IGLU scales well to deeper models and larger datasets, and continues to give performance boosts with state-of-the-art even when compared to these additional baselines. Experimental evidence also demonstrates IGLU’s ability to generalize easily to a different setting and architecture, perform equivalently with smooth activation functions and achieve lower training losses faster compared to all the baselines across datasets.
3. Section C gives detailed proofs of Lemma 1 and the convergence rates offered by IGLU. It is shown that under standard assumptions such as objective smoothness, IGLU is able to offer both the standard rate of convergence common for SGD-style algorithms, as well as a fast rate if full-batch GD is performed.

heading: A ADDITIONAL IMPLEMENTATION DETAILS
text: We recall from Section 4.1 that the wall-clock time reported in Figure 2 consists of strictly the optimization time for each method and excludes method-specific overheads. This was actually a disadvantage for IGLU since its overheads are relatively mild. This section demonstrates that other baselines incur much more significant overheads whereas IGLU does not suffer from these large overheads. When included, it further improves the speedups that IGLU provides over baseline methods.

heading: A.1 HARDWARE
text: We implement IGLU in TensorFlow 1.15.2 and perform all experiments on an NVIDIA V100 GPU (32 GB Memory) and Intel Xeon CPU processor (2.6 GHz). We ensure that all baselines are experimented with the exact same hardware.

heading: A.2 TIME AND MEMORY OVERHEADS FOR VARIOUS METHODS
text: We consider three main types of overheads that are incurred by different methods. This includes pre-processing overhead that is one-time, recurring overheads and additional memory overhead. We describe each of the overheads in the context of the respective methods below.
1. GraphSAGE - GraphSAGE recursively samples neighbors at each layer for every minibatch. This is done on-the-fly and contributes to a significant sampling overhead. Since this overhead is incurred for every minibatch, we categorize this under recurring overhead. We aggregate this overhead across all minibatches during training. GraphSAGE does not incur preprocessing or additional memory overheads.
2. VRGCN - Similar to GraphSAGE, VRGCN also recursively samples neighbors at each layer for every minibatch on-the-fly. We again aggregate this overhead across all minibatches during training. VRGCN also stores the stale/historical embeddings that are learnt for every node at every layer. This is an additional overhead of O (NKd), where K is the number of layers, N is the number of nodes in the graph and d = 1K ∑ k∈[K] dk is the average
embedding dimensionality across layers. 3. ClusterGCN - ClusterGCN creates subgraphs and uses them as minibatches for training.
For the creation of these subgraphs, ClusterGCN performs graph clustering using the highly optimized METIS tool†. This overhead is a one-time overhead since graph clustering is done before training and the same subgraphs are (re-)used during the whole training process. We categorize this under preprocessing overhead. ClusterGCN does not incur any recurring or additional memory overheads.
4. GraphSAINT - GraphSAINT, similar to ClusterGCN creates subgraphs to be used as minibatches for training. We categorize this minibatch creation as the preprocessing overhead for GraphSAINT. However, unlike ClusterGCN, GraphSAINT also periodically creates new subgraphs on-the-fly. We categorize this overhead incurred in creating new subgraphs as recurring overhead. GraphSAINT does not incur any additional memory overheads.
5. IGLU - IGLU creates mini-batches only once using subsets of nodes with their full neighborhood information which is then reused throughout the training process. In addition to this, IGLU requires initial values of both the incomplete gradients αk and the Xk embeddings (Step 3 and first part of Step 4 in Algorithm 1 and Step 2 and Step 4 in Algorithm 2) before optimization can commence. We categorize these two overheads - mini-batch creation and initializations of αk’s and Xk embeddings as IGLU’s preprocessing overhead and note that IGLU does not have any recurring overheads. IGLU does incur an additional memory overhead since it needs to store the incomplete gradients αk’s in the inverted variant and the embeddings Xk in the backprop variant. However, note that the memory occupied by Xk for a layer k is the same as that occupied by αk for that layer (see Definition 1). Thus, for both its variants, IGLU incurs an additional memory overhead of O (NKd), where K is the number of layers, N is the number of nodes in the graph and d = 1K ∑ k∈[K] dk is the
average embedding dimensionality across layers.
Tables 3 and 4 report the overheads incurred by different methods on the OGBN-Proteins and Reddit datasets (the largest datasets in terms of edges and nodes respectively). ClusterGCN runs into a runtime error on the Proteins dataset (denoted by || in the table). N/A stands for Not Applicable in the tables. In the tables, specifically for IGLU, pre-processing time is the sum of initialization time required to pre-computeαk, Xk and mini-batch creation time. We also report the individual overhead for both initialization and mini-batch creation for IGLU. The total pre-processing time for IGLU is denoted by T, overheads incurred for initialization by I and overheads incurred for mini-batch creation by MB.
†http://glaros.dtc.umn.edu/gkhome/metis/metis/download
For IGLU, the minibatch creation code is currently implemented in Python, while GraphSAINT uses a highly optimized C++ implementation. Specifically for the Reddit dataset, the number of subgraphs that GraphSAINT samples initially is sufficient and it does not incur any recurring overhead. However, on Proteins, GraphSAINT samples 200 subgraphs every 18 epochs, once the initially sampled subgraphs are used, leading to a sizeable recurring overhead.
Discussion on OGBN - Proteins: Figure 5 summarizes the total overheads for all the methods. On the OGBN-Proteins dataset, IGLU incurs ∼ 2× less overhead than GraphSAINT, the best baseline, while incurring as much as ∼ 7.8× less overhead than VRGCN and ∼ 4.7× less overhead than GraphSAGE. It is also important to note that these overheads are often quite significant as compared to the optimization time for the methods and can add to the overall experimentation time. For experiments on the OGBN-Proteins dataset, VRGCN’s total overheads equal 46.25% of its optimization time, GraphSAINT’s overheads equal 19.63% of its optimization time, GraphSAGE’s overhead equal 9.52% of its optimization time. However IGLU’s total overheads equal only 5.59% of its optimization time which is the lowest out of all methods. Upon re-computing the speedup provided by IGLU using the formula defined in equation (1), but this time with overheads included, it was observed that IGLU offered an improved speedup of 16.88% (the speedup was only 11.05% when considering only optimization time).
Discussion on Reddit: On the Reddit dataset, IGLU incurs upto ∼ 2.1× less overhead than VRGCN and upto ∼ 4.5× less overhead than GraphSAGE. However, IGLU incurs marginally higher overhead (∼ 1.4×) than GraphSAINT. This can be attributed to the non-optimized minibatch creation
routine currently used by IGLU compared to a highly optimized and parallelized implementation in C++ used by GraphSAINT. This is an immediate avenue for future work. Nevertheless, VRGCN’s total overhead equals 15.41% of its optimization time, GraphSAINT’s overhead equals 43.41% of its optimization time, GraphSAGE’s overhead equals 31.25% of its optimization time, ClusterGCN’s overhead equals 41.02% of its optimization time while IGLU’s overhead equals 44.09% of its optimization time. Whereas the relative overheads incurred by IGLU and GraphSAINT in comparison to optimization time may seem high for this dataset, this is because the actual optimization times for these methods are rather small, being just 15.43 and 20.86 seconds for GraphSAINT and IGLU respectively in comparison to the other methods such as VRGCN, ClusterGCN and GraphSAGE whose optimization times are 124s, 131s and 133s respectively, almost an order of magnitude larger than that of IGLU.

heading: A.3 MEMORY ANALYSIS FOR IGLU
text: While IGLU requires storing stale variables which can have additional memory costs, for most scenarios with real world graphs, saving these stale representations on modern GPUs are quite reasonable. We provide examples of additional memory usage required for two of the large datasets - Reddit and Proteins in Table 6 and we observe that IGLU requires only 150MB and 260MB of additional GPU memory. Even for a graph with 1 million nodes, the additional memory required would only be ∼ 2.86GB which easily fit on modern GPUs. For even larger graphs, CPU-GPU
interfacing can be used. CPU and GPU interfacing for data movement is a common practice in training machine learning models and hence a potential method to mitigate the issue of limited memory availability in settings with large datasets. This has been explored by many works for dealing with large datasets in the context of GCNs, such as VR-GCN (Chen et al., 2018a) for storing historical activations in main memory (CPU). Such an interfacing is an immediate avenue of future work for IGLU.
We observe that IGLU enjoys significant speedups and improvements in training cost across datasets as compared to the baselines as a result of using stale variables. Additionally, since IGLU requires training just a single layer at a time, there is score for further reduction in memory usage by using only the variables required for the current layer and by re-using the computation graphs across layers, and therefore making the IGLU even less memory expensive.

heading: A.4 HYPERPARAMETER CONFIGURATIONS FOR IGLU AND BASELINES
text: Table 5 summarizes the source URLs and commit stamps using which baseline methods were obtained for experiments. The Adam optimizer Kingma & Ba (2014) was used to train IGLU and all the baselines until convergence for each of the datasets. A grid search was performed over the other hyper-parameters for each baseline which are summarized below:
1. GraphSAGE: Learning Rate - {0.01, 0.001}, Batch Size - {512, 1024, 2048}, Neighborhood Sampling Size - {25, 10, 5}, Aggregator - {Mean, Concat}
2. VRGCN : Batch Size - {512, 1000, 2048}, Degree - {1, 5, 10}, Method - {CV, CVD}, Dropout - {0, 0.2, 0.5}
3. ClusterGCN : Learning Rate - {0.01, 0.001, 0.005}, Lambda - {-1, 1, 1e-4}, Number of Clusters - {5, 50, 500, 1000, 1500, 5000}, Batch Size - {1, 5, 50, 100, 500}, Dropout - {0, 0.1, 0.3, 0.5}
4. GraphSAINT-RW : Aggregator - {Mean, Concat}, Normalization - {Norm, Bias}, Depth - {2, 3, 4}, Root Nodes - {1250, 2000, 3000, 4500, 6000}, Dropout - {0.0, 0.2, 0.5}
5. IGLU : Learning Rate - {0.01, 0.001} with learning rate decay schemes, Batch Size - {512, 2048, 4096, 10000}, Dropout - {0.0, 0.2, 0.5, 0.7}
A.5 INCORPORATING RESIDUAL CONNECTIONS, BATCH NORMALIZATION AND VIRTUAL NODES IN IGLU
General-purpose techniques such as BatchNorm and skip/residual connections, and GCN-specific advancements such as bi-level aggregation using virtual nodes offer performance boosts. The current implementation of IGLU already incorporates normalizations as described in Sections 3 and 4. Below we demonstrate how all these aforementioned architectural variations can be incorporated into IGLU with minimal changes to Lemma 1 part 2. Remaining guarantees such as those offered by Lemma 1 part 1 and Theorem 2 remain unaltered but for changes to constants.
Incorporating BatchNorm into IGLU: As pointed out in Section 3, IGLU assumes a general form of the architecture, specifically
xki = f(x k−1 j , j ∈ {i} ∪ N (i);E k),
where f includes the aggregation operation such as using the graph Laplacian, weight matrices and any non-linearity. This naturally allows operations such as normalizations (LayerNorm, BatchNorm etc) to be carried out. For instance, the f function for BatchNorm with a standard GCN would look
like the following
zki = σ ∑ j∈V Aij(W k)>xk−1j  ẑki =
zki −µkB√ νkB +
xki = Γ k · ẑki + βk,
where σ denotes a non-linearity like the sigmoid, Γk ∈ Rdk×dk is a diagonal matrix and βk ∈ Rdk is a vector, andµkB , ν k B ∈ Rdk are vectors containing the dimension-wise mean and variance values over a mini-batch B (division while computing ẑki is performed element-wise). The parameter E k is taken
to collect parameters contained in all the above operations i.e. Ek = { W k,Γk,βk } in the above
example ( µkB , ν k B are computed using samples in a mini-batch itself). Downstream calculations in Definition 1 and Lemma 1 continue to hold with no changes.
Incorporating Virtual Nodes into IGLU: Virtual nodes can also be seamlessly incorporated simply by re-parameterizing the layer-wise parameter Ek. We are referring to (Pei et al., 2020) [Pei et al ICLR 2020] for this discussion. Let R be the set of relations and let Ng(·),Ns(·) denote graph neighborhood and latent-space neighborhood functions respectively. A concrete example is given below to illustrate how virtual nodes can be incorporated. We note that operations like BatchNorm etc can be additionally incorporated and alternate aggregation operations e.g. concatenation can be used instead.
g (k,r) i = σ  ∑ j∈Ng(i) τ(i,j)=r Aij(L k g) >xk−1j 
s (k,r) i = σ  ∑ j∈Ns(i) τ(i,j)=r Tij(L k s) >xk−1j  (Low-level aggregation)
mki = 1 |R| ∑ r∈R ( g (k,r) i + s (k,r) i ) (High-level aggregation)
xki = σ ( (W k)>mki ) (Non-linear transform)
where τ denotes the relationship indicator, Aij denotes the graph edge weight between nodes i and j and Tij denotes their geometric similarity in latent space. Note that g (k,r) i , s (k,r) i corresponds to embeddings of the virtual nodes in the above example. To implement the above, the parameter Ek can be taken to collect the learnable parameters contained in all the above operations i.e. Ek ={ Lkg , L k s ,W k }
in the above example. Downstream calculations in Definition 1 and Lemma 1 continue to hold as is with no changes.
Incorporating Skip/Residual Connections into IGLU: The architecture style presented in Section 3 does not directly allow skip connections but they can be incorporated readily with no change to Definition 1 and minor changes to Lemma 1. Let us introduce the notation k → m to denote a direct forward (skip) connection directed from layer k to some layer m > k. In a purely feed-forward style architecture, we would only have connections of the form k → k + 1. The following gives a simple example of a GCN with a connection that skips two layers, specifically (k − 2)→ k.
xki = f ( xk−1j , j ∈ {i} ∪ N (i);E k) ) + xk−2i ,
where f includes the aggregation operation with the graph Laplacian, the transformation weight matrix and any non-linearity. Definition 1 needs no changes to incorporate such architectures. Part 1 of Lemma 1 also needs no changes to address such cases. Part 2 needs a simple modification as shown below
Lemma 3 (Lemma 1.2 adapted to skip connections). For the final layer, we continue to have (i.e. no change) αK = G(WK+1)>. For any k < K, we have αk = ∑ m:k→m ∂(αm Xm) ∂Xk ∣∣∣ allαm s.t. k→m
i.e. αkjp = ∑ m:k→m ∑ i∈V ∑dm q=1 α m iq · ∂Xmiq ∂Xkjp . Note that as per the convention established in the paper, ∑ m:k→m ∂(αm Xm) ∂Xk ∣∣∣ allαm s.t. k→m
implies that while taking the derivatives, αm values are fixed (treated as a constant) for all m > k such that k → m. This “conditioning” is important since αm also indirectly depends on Xk if m > k. Proof of Lemma 1.2 adapted to skip connections: We consider two cases yet again and use Definition 1 that tells us that
αkjp = ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂XKjp
Case 1 (k = K): Since this is the top-most layer and there are no connections going ahead let alone skipping ahead, the analysis of this case remains unchanged and continues to yield αK = G(WK+1)>.
Case 2 (k < K): Using Definition 1 and incorporating all layers to which layer k has a direct or skip connection gives us
αkjp = ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xkjp = ∑ m:k→m ∑ i∈V ∑ c∈[C] gic · ∑ l∈V dm∑ q=1 ∂ŷic ∂Xmlq ∂Xmlq ∂Xkjp
Rearranging the terms gives us
αkjp = ∑
m:k→m ∑ l∈V dm∑ q=1 ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xmlq  · ∂Xmlq ∂Xkjp = ∑ m:k→m ∑ l∈V dm∑ q=1 αmlq · ∂Xmlq ∂Xkjp ,
where we simply used Definition 1 in the second step. However, the resulting term simply gives us αkjp = ∑ m:k→m ∂(αm Xm) ∂Xkjp ∣∣∣ allαm s.t. k→m which conditions on, or treats as a constant, the term αm for all m > k such that k → m according to our notation convention. This finishes the proof of part 2 adapted to skip connections.

heading: A.6 COMPARISON OF IGLU WITH VR-GCN, MVS-GNN AND GNNAUTOSCALE
text: We highlight the key difference between IGLU and earlier works that cache intermediate results for speeding up GNN training below.
A.6.1 VRGCN V/S IGLU
1. Update of Cached Variables: VR-GCN (Chen et al., 2018a) caches only historical embeddings, and while processing a single mini-batch these historical embeddings are updated for a sampled subset of the nodes. In contrast IGLU does not update any intermediate results after processing each mini-batch. These are updated only once per epoch, after all parameters for individual layers have been updated.
2. Update of Model Parameters: VR-GCN’s backpropagation step involves update of model parameters of all layers after each mini-batch. In contrast IGLU updates parameters of only a single layer at a time.
3. Variance due to Sampling: VR-GCN incurs additional variance due to neighborhood sampling which is then reduced by utilizing historical embeddings for some nodes and by computing exact embeddings for the others. IGLU does not incur such variance since IGLU uses all the neighbors.

heading: A.6.2 MVS-GNN V/S IGLU
text: MVS-GNN (Cong et al., 2020) is another work that caches historical embeddings. It follows a nested training strategy wherein firstly a large batch of nodes are sampled and mini-batches are further
created from this large batch for training. MVS-GNN handles variance due to this mini-batch creation by performing importance weighted sampling to construct mini-batches.
1. Update of Cached Variables and Variance due to Sampling: Building upon VR-GCN, to reduce the variance in embeddings due to its sampling of nodes at different layers, MVS-GNN caches only embeddings and uses historical embeddings for some nodes and recompute the embeddings for the others. Similar to VR-GCN, these historical embeddings are updated as and when they are part of the mini-batch used for training. As discussed above, IGLU does not incur such variance since IGLU uses all the neighbors.
2. Update of Model Parameters: Update of model parameters in MVS-GNN is similar to that of VR-GCN, where backpropagation step involves update of model parameters of all layers for each mini-batch. As described already, IGLU updates parameters of only a single layer at a time.

heading: A.6.3 GNNAUTOSCALE V/S IGLU
text: GNNAutoScale (Fey et al., 2021) extends the idea of caching historical embeddings from VR-GCN and provides a scalable solution.
1. Update of intermediate representations and model parameters: While processing a minibatch of nodes, GNNAutoScale computes the embeddings for these nodes at each layer while using historical embeddings for the immediate neighbors outside the current minibatch. After processing each mini-batch, GNNAutoScale updates the historical embeddings for nodes considered in the mini-batch. Similar to VR-GCN and MVS-GNN, GNNAutoScale updates all parameters at all layers while processing a mini-batch of nodes. In contrast IGLU does not update intermediate results (intermediate representations in Algorithm 1 and incomplete gradients in Algorithm 2) after processing each minibatch. In fact, these are updated only once per epoch, after all parameters for individual layers have been updated.
2. Partitioning: GNNAutoScale relies on the METIS clustering algorithm for creating minibatches that minimize inter-connectivity across batches. This is done to minimize access to historical embeddings and reduce staleness. This algorithm tends to bring similar nodes together, potentially resulting in the distributions of clusters being different from the original dataset. This may lead to biased estimates of the full gradients while training using minibatch SGD as discussed in Section 3.2, Page 5 of Cluster-GCN (Chiang et al., 2019). IGLU does not rely on such algorithms since it’s parameter updates are concerned with only a single layer and also avoids potential additional bias.
Similarity of IGLU with GNNAutoScale: Both of the methods avoid a neighborhood sampling step, thereby avoiding additional variance due to neighborhood sampling and making use of all the edges in the graph. Both IGLU and GNNAutoScale propose methods to reduce the neighborhood explosion problem, although in fundamentally different manners. GNNAutoScale does so by pruning the computation graph by using historical embeddings for neighbors across different layers. IGLU on the other hand restricts the parameter updates to a single layer at a time by analyzing the gradient structure of GNNs therefore alleviating the neighborhood explosion problem.

heading: A.6.4 SUMMARY OF IGLU’S TECHNICAL NOVELTY AND CONTRAST WITH CACHING BASED RELATED WORKS
text: To summarize, IGLU is fundamentally different from these methods that cache historical embeddings in that it changes the entire training procedure of GCNs in contrast with the aforementioned caching based methods as follows:
• The above methods still follow standard SGD style training of GCNs in that they update the model parameters at all the layers after each mini-batch. This is very different from IGLU’s parameter updates that concern only a single layer while processing a mini-batch.
• IGLU can cache either incomplete gradients or embeddings which is different from the other approaches that cache only embeddings. This provides alternate approaches for training GCNs and we demonstrate empirically that caching incomplete gradients, in fact, offers superior performance and convergence.
• Unlike GNNAutoScale and VR-GCN that update some of the historical embeddings after each mini-batch is processed, IGLU’s caching is much more aggressive and the stale variables are updated only once per epoch, after all parameters for all layers have been updated.
• Theoretically, we provide good convergence rates and bounded bias even while using stale gradients, which has not been discussed in any prior works.
These are the key technical novelties of our proposed method and they are a consequence of a careful understanding of the gradient structure of GCN’s themselves.

heading: A.6.5 EMPIRICAL COMPARISON WITH GNNAUTOSCALE
text: We provide an empirical comparison of IGLU with GNNAutoScale and summarize the results in Table 7 and Figure 6. It is important to note that the best results for GNNAutoScale as reported by the authors in the paper, correspond to varying hyperparameters such as number of GNN layers and different embedding dimensions across methods, datasets and architectures. However, for the experiments covered in the main paper, we use 2 layer settings for PPI-Large, Flickr and Reddit and 3 layer settings for OGBN-Arxiv and OGBN-Proteins datasets consistently for IGLU and the baseline methods, as motivated by literature. We also ensure that the embedding dimensions are uniform across IGLU and the baselines. Therefore, to ensure a fair comparison, we perform additional experiments with these parameters for GNNAutoScale set to values that are consistent with our experiments for IGLU and the baselines. We train GNNAutoScale with three variants, namely GCN, GCNII (Chen et al., 2020) and PNA (Corso et al., 2020) and report the results for each of the variant. We also note here that GNNAutoScale was implemented in PyTorch (Paszke et al., 2019) while IGLU was implemented in TensorFlow (Abadi et al., 2016). While this makes a wall-clock time comparison unsuitable as discussed in Appendix B.2, we still provide a wall-clock time comparison for completeness. We also include the best performance numbers for GNNAutoScale on these datasets (as reported by the authors in Table 5, Page 9 of the GNNAutoScale paper) across different architectures. Note that we do not provide comparisons on the OGBN-Proteins dataset since we ran into errors while trying to incorporate the dataset into the official implementation of GNNAutoScale.
Results: Figure 6 provides convergence plots comparing IGLU with the different architectures of GNNAutoScale and Table 7 summarizes the test performance on PPI-Large, Flickr, Reddit and OGBNArxiv (transductive) datasets. From the table, we observe that IGLU offers competitive performance compared to the GCN variant of GAS for the majority of the datasets. We also observe from Figure 6 that IGLU offers significant improvements in training time with rapid early convergence on the validation set. We note that more complex architectures such as GCNII and PNA offer improvements in performance to GNNAutoScale. IGLU being architecture agnostic can be incorporated with these architectures for further improvements in performance. We leave this as an avenue for future work.

heading: A.7 DETAILED DESCRIPTION OF ALGORITHMS 1 AND 2
text: We present Algorithms 1 and 2 again below (as Algorithms 3 and 4 respectively) with details of each step.
IGLU: backprop order Algorithm 3 implements the IGLU algorithm in its backprop variant. Node embeddings Xk, k ∈ [K] are calculated and kept stale. They are not updated even when model parameters Ek, k ∈ [K] get updated during the epoch. On the other hand, the incomplete task gradients αk, k ∈ [K] are kept refreshed using the recursive formulae given in Lemma 1. For sake of simplicity, the algorithm been presented with staleness duration of one epoch i.e. Xk are refreshed at the beginning of each epoch. Variants employing shorter or longer duration of staleness can be also explored simply by updating Xk, k ∈ [K] say twice in an epoch or else once every two epochs. IGLU: inverted order Algorithm 4 implements the IGLU algorithm in its inverted variant. Incomplete task gradients αk, k ∈ [K] are calculated once at the beginning of every epoch and kept stale. They are not updated even when node embeddings Xk, k ∈ [K] get updated during the epoch. On the other hand, the node embeddings Xk, k ∈ [K] are kept refreshed. For sake of simplicity, the algorithm been presented with staleness duration of one epoch i.e. αk are refreshed at the beginning of each epoch. Variants employing shorter or longer durations of staleness can be also explored

heading: Algorithm PPI-Large Reddit Flickr Arxiv (Trans)
text: simply by updating αk say twice in an epoch or else once every two epochs. This has been explored in Section 4.2 (see paragraph on "Analysis of Degrees of Staleness.").
Algorithm 3 IGLU: backprop order
Input: GCN G, initial features X0, task loss L 1: Initialize model parameters Ek, k ∈ [K],WK+1 2: for epoch = 1, 2, . . . do 3: for k = 1 . . .K do 4: Refresh Xk ← f(Xk−1;Ek) //Xk will be kept stale till next epoch 5: end for 6: Ŷ ← XKWK+1 //Predictions 7: G← [ ∂`i ∂ŷic ] N×C //The loss derivative matrix
8: Compute ∂L ∂WK+1 ← (XK)>G //Using Lemma 1.1 here 9: Update WK+1 ←WK+1 − η · ∂L
∂WK+1
10: Refresh αK ← G(WK+1)> //Using Lemma 1.2 here 11: for k = K . . . 2 do 12: Compute ∂L
∂Ek ← ∂(α k Xk) ∂Ek ∣∣∣ αk , //Using Lemma 1.1 here
13: Update Ek ← Ek − η · ∂L ∂Ek 14: Refresh αk−1 ← ∂(α k Xk)
∂Xk−1 ∣∣∣ αk
//Using Lemma 1.2 here 15: end for 16: end for
Algorithm 4 IGLU: inverted order
Input: GCN G, initial features X0, task loss L 1: Initialize model parameters Ek, k ∈ [K],WK+1 2: for k = 1 . . .K do 3: Xk ← f(Xk−1;Ek) //Do an initial forward pass 4: end for 5: for epoch = 1, 2, . . . do 6: Ŷ ← XKWK+1 //Predictions 7: G← [ ∂`i ∂ŷic ] N×C //The loss derivative matrix
//Use Lemma 1.2 to refresh αk, k ∈ [K]. 8: Refresh αK ← G(WK+1)> 9: for k = K . . . 2 do
10: Refresh αk−1 ← ∂(α k Xk)
∂Xk−1 ∣∣∣ αk
11: end for //These αk, k ∈ [K] will now be kept stale till next epoch 12: for k = 1 . . .K do 13: Compute ∂L
∂Ek ← ∂(α k Xk) ∂Ek ∣∣∣ αk , //Using Lemma 1.1 here
14: Update Ek ← Ek − η · ∂L ∂Ek 15: Refresh Xk ← f(Xk−1;Ek) 16: end for 17: Compute ∂L
∂WK+1 ← (XK)>G //Using Lemma 1.1 here
18: Update WK+1 ←WK+1 − η · ∂L ∂WK+1 19: end for

heading: A.8 OGBN-PROTEINS: VALIDATION PERFORMANCE AT A MINI-BATCH LEVEL
text: The performance analysis in the main paper was plotted at a coarse granularity of an epoch. We refer to an epoch as one iteration of steps 6 to 18 (both inclusive) in Algorithm 4. For a finer analysis of IGLU’s performance on the OGBN-Proteins dataset, we measure the Validation ROC-AUC at the granularity of a mini-batch. As mentioned in the “SGD Implementation" paragraph in Page 5 below the algorithm description in the main paper, steps 14 and 18 in Algorithm 4 are implemented using mini-batch SGD. Recall that OGBN-Proteins used a 3 layer GCN. For the proteins dataset we have ∼ 170 mini-batches for training. We update the parameters for each layer using all the mini-batches
from the layer closest to the input to the layer closest to the output as detailed in Algorithm 4. To generate predictions, we compute partial forward passes after the parameters of each layer is updated. We plot the validation ROC-AUC as the first epoch progresses in Figure 7. We observe that when the layer closest to the input is trained (Layer 1 in the figure), IGLU has an ROC-AUC close to 0.5 on the validation set. Subsequently once the second layer (Layer 2 in the figure) is trained, we observe that the validation ROC-AUC improves from around ∼ 0.51 to ∼ 0.57 and finally once the layer closest to the output is trained (Layer 3 in the figure), the ROC-AUC progresses quickly to a high validation ROC-AUC of ∼ 0.81. In the figure in the main paper the high validation score reflects the result at the end of the first epoch. Training the GCN using a total of ∼ 510 minibatches (∼ 170 per layer) approximately takes 5 seconds as depicted in Figure 2 in the main paper.

heading: B DATASET STATISTICS AND ADDITIONAL EXPERIMENTAL RESULTS

heading: B.1 DATASET STATISTICS
text: Table 8 provides details on the benchmark node classification datasets used in the experiments. The following five benchmark datasets were used to empirically demonstrate the effectiveness of IGLU: predicting the communities to which different posts belong in Reddit‡ (Hamilton et al., 2018), classifying protein functions across various biological protein-protein interaction graphs in PPI-Large§ (Hamilton et al., 2018), categorizing types of images based on descriptions and common properties in Flickr¶ (Zeng et al., 2020), predicting paper-paper associations in OGBN-Arxiv|| (Hu et al., 2021) and categorizing meaningful associations between proteins in OGBN-Proteins** (Hu et al., 2021).

heading: B.2 COMPARISON WITH ADDITIONAL BASELINES
text: In addition to the baselines mentioned in Table 1, Table 9 compares IGLU to LADIES (Zou et al., 2019), L2-GCN (You et al., 2020), AS-GCN (Huang et al., 2018), MVS-GNN (Cong et al., 2020), FastGCN (Chen et al., 2018b), SIGN (Frasca et al., 2020), PPRGo (Bojchevski et al., 2020) and Bandit Sampler’s (Liu et al., 2020) performance on the test set. However, a wall-clock time comparison with these methods is not provided since the author implementations of LADIES, L2GCN, MVS-GNN and
‡http://snap.stanford.edu/graphsage/reddit.zip §http://snap.stanford.edu/graphsage/ppi.zip ¶https://github.com/GraphSAINT/GraphSAINT - Google Drive Link ||https://ogb.stanford.edu/docs/nodeprop/#ogbn-arxiv
**https://ogb.stanford.edu/docs/nodeprop/#ogbn-proteins
SIGN are in PyTorch (Paszke et al., 2019) which has been shown to be less efficient than Tensorflow (Chiang et al., 2019; Abadi et al., 2016) for GCN applications. Also, the official AS-GCN, FastGCN and Bandit Sampler implementations released by the authors were for 2 layer models only, whereas some datasets such as Proteins and Arxiv require 3 layer models for experimentation. Attempts to generalize the code to a 3 layer model ran into runtime errors, hence the missing results are denoted by ** in Table 9 and these methods are not considered for a timing analysis. MVS-GNN also runs into a runtime error on the Proteins dataset denoted by ||. IGLU continues to significantly outperform all additional baselines on all the datasets.

heading: B.3 TIMING ANALYSIS FOR SCALING TO MORE LAYERS
text: To compare the scalability and performance of different algorithms for deeper models, models with 2, 3 and 4 layers were trained for IGLU and the baseline methods. IGLU was observed to offer a per-epoch time that scaled roughly linearly with the number of layers as well as offer the highest gain in test performance as the number of layers was increased.

heading: B.3.1 PER EPOCH TRAINING TIME
text: Unlike neighbor sampling methods like VRGCN and GraphSAGE, IGLU does not suffer from the neighborhood explosion problem as the number of layers increases, since IGLU updates involve only a single layer at any given time. We note that GraphSAINT and ClusterGCN also do not suffer from the neighborhood explosion problem directly since they both operate by creating GCNs on subgraphs. However, these methods may be compelled to select large subgraphs or else suffer from poor convergence. To demonstrate IGLU’s effectiveness in solving the neighborhood explosion problem, the per-epoch training times are summarized as a function of the number of layers in Table
10 on the Proteins dataset. A comparison with ClusterGCN could not be provided as since it ran into runtime errors on this dataset. In Table 10, while going from 2 to 4 layers, GraphSAINT was observed to require ∼ 1.6× more time per epoch while IGLU required ∼ 2.3× more time per epoch, with both methods scaling almost linearly with respect to number of layers as expected. However, GraphSAGE suffered a ∼ 62× increase in time taken per epoch in this case, suffering from the neighborhood explosion problem. VRGCN ran into a run-time error for the 4 layer setting (denoted by || in Table 10). Nevertheless, even while going from 2 layers to 3 layers, VRGCN and GraphSAGE are clearly seen to suffer from the neighborhood explosion problem, resulting in an increase in training time per epoch of almost ∼ 9.4× and ∼ 5.6× respectively. We note that the times for GraphSAINT in Table 10 are significantly smaller than those of IGLU even though earlier discussion reported IGLU as having the fastest convergence. However, there is no contradiction – GraphSAINT operates with very large subgraphs, with each subgraph having almost 10 % of the nodes of the entire training graph (∼ 10000 nodes in a minibatch), while IGLU operates with minibatches of size 512, resulting in IGLU performing a lot more gradient updates within an epoch, as compared with GraphSAINT. Consequently, IGLU also takes fewer epochs to converge to a better solution than GraphSAINT, thus compensating for the differences in time taken per epoch.

heading: B.3.2 TEST CONVERGENCE AUC
text: This section explores the effect of increasing the number of layers on convergence rate, optimization time, and final accuracy, when scaling to larger number of layers. To jointly estimate the efficiency of a method in terms of wall clock time and test performance achieved, the area under the test convergence plots (AUTC) for various methods was computed. A method that converged rapidly and that too to better performance levels would have a higher AUTC than a method that converges to suboptimal values or else converges very slowly. To fairly time all methods, each method was offered time that was triple of the time it took the best method to reach its highest validation score. Defining the cut-off time this way ensures that methods that may not have rapid early convergence still get a fair chance to improve by having better final performance, while also simultaneously penalizing methods that may have rapid early convergence but poor final performance. We rescale the wall clock time to be between 0 and 1, where 0 refers to the start of training while 1 refers to the cut-off time.
Results: Figure 8 summarizes the AUTC values. IGLU consistently obtains higher AUC values than all methods for all number of layers demonstrating its stability, rapid convergence during early phases of training and ability to generalize better as compared to other baselines. GraphSAGE suffered from neighborhood explosion leading to increased training times and hence decreased AUC values as the number of layers increase. VR-GCN also suffered from the same issue, and additionally ran into a run-time error for 4 layer models.
Test Performance with Increasing Layers: Table 11 summarizes the final test performances for different methods, across different number of layers. The performance for some methods is inconsistent as the depth increases whereas IGLU consistently outperforms all the baselines in this case as well with gains in performance as we increase the number of layers hence making it an attractive technique to train deeper GCN models.

heading: B.4 SCALABILITY TO LARGER DATASETS: OGBN - PRODUCTS
text: To demonstrate the ability of IGLU to scale to very large datasets, we performed experiments on the OGBN-Products dataset (Hu et al., 2021), one of the largest datasets in the OGB collection, with its statistics summarized in Table 12.
To demonstrate scalability, we conducted experiments in the transductive setup since this setup involves using the full graph. In addition, this was the original setup in which the dataset was
benchmarked, therefore allowing for a direct comparison with the baselines (results taken from Table 4 in OGB (Hu et al., 2021) and Table 6 in SIGN (Frasca et al., 2020)). We summarize the performance results in the Table 13, reporting Micro-F1 as the metric. We however do not provide timing comparisons with the baseline methods since IGLU is implemented in TensorFlow while the baselines in the original benchmark Hu et al. (2021) are implemented in PyTorch. This therefore renders a direct comparison of wall-clock time unsuitable. Please refer to Appendix B.2 for more details.
We observe that IGLU is able to scale to the OGBN-Products dataset with over 2.4 million nodes and outperforms all of the baseline methods.

heading: B.5 APPLICABILITY OF IGLU IN THE TRANSDUCTIVE SETTING: OGBN-ARXIV AND OGBN-PROTEINS
text: In addition to the results in the inductive setting reported in the main paper, we perform additional experiments on the OGBN-Arxiv and OGBN-Proteins dataset in the transductive setting to demonstrate IGLU’s applicability across inductive and transductive tasks and compare the performance of IGLU to that of transductive baseline methods in Table 14 (results taken from OGB (Hu et al., 2021)).
We observe that even in the transductive setting, IGLU outperforms the baseline methods on both the OGBN-Arxiv and OGBN-Proteins datasets.

heading: B.6 ARCHITECTURE AGNOSTIC NATURE OF IGLU
text: To demonstrate the applicability of IGLU to a wide-variety of architectures, we perform experiments on IGLU with Graph Attention Networks (GAT) (Veličković et al., 2017), GCN (Kipf & Welling, 2017) and GraphSAGE (Hamilton et al., 2018) based architectures and summarize the results for the same in Table 15 and 16 respectively as compared to these baseline methods. We use the Cora, Citeseer and Pubmed datasets as originally used in Veličković et al. (2017) for comparison with the GAT based architecture and the OGBN-Arxiv dataset for comparison with GCN and GraphSAGE based architectures (baseline results as originally reported in (Hu et al., 2021)).
We observe that the IGLU+GAT, IGLU+GCN and IGLU+GraphSAGE outperforms the baseline methods across datasets, thereby demonstrating the architecture agnostic nature of IGLU.

heading: B.7 EXPERIMENTS USING A SMOOTH ACTIVATION FUNCTION: GELU
text: To understand the effect of using non-smooth vs smooth activation functions on IGLU, we perform experiments using the GELU (Hendrycks & Gimpel, 2020) activation function which is a smooth function and in-line with the objective smoothness assumptions made in Theorem 2.
GELU(x) = xP (X ≤ x) = xΦ(x) = x · 1 2
[1 + erf(x/ √ 2)]
We compare the performance of IGLU using ReLU and GELU on all the 5 datasets in the main paper and summarize the results in Table 17.
We observe that IGLU is able to enjoy very similar performance across both GELU and ReLU as the activation functions, thereby justifying the practicality of our smoothness assumptions.

heading: B.8 ANALYSIS ON DEGREE OF STALENESS: BACKPROP ORDER OF UPDATES
text: To understand the effect of more frequent updates in the backprop variant, we performed additional experiments using the backprop variant on the PPI-Large dataset and varied the frequency of updates, the results of which are reported in Table 18 for a single run. We fix the hyperparameters and train for 200 epochs.
We observe that more frequent updates help stabilize training better. We also observe that update frequencies 1 and 2 perform competitively, and both significantly outperform update frequency 0.5.
However we note that with higher update frequency, we incur an additional computational cost since we need to re-compute embeddings more frequently. We believe that both improved stability and competitive performance can be attributed to the fresher embeddings, which are otherwise kept stale
within an epoch in this order of updates. The experiment with frequency 0.5 has a slower convergence and comparatively poor performance as expected.

heading: B.9 CONVERGENCE - TRAIN LOSS VS WALL CLOCK TIME
text: Figure 9 provides train loss curves for all datasets and methods in Table 1. IGLU is able to achieve a lower training loss faster compared to all the baselines across datasets.

heading: C THEORETICAL PROOFS
text: In this section we provide a formal restatement of Theorem 2 as well as its proof. We also provide the proof of Lemma 1 below.

heading: C.1 PROOF FOR LEMMA 1
text: Proof of Lemma 1 We consider the two parts of Lemma 1 separately and consider two cases while proving each part. For each part, Case 1 pertains to the final layer and Case 2 considers intermediate layers.
Clarification about some Notation in the statements of Definition 1 and Lemma 1: The notation ∂(G Ŷ) ∂Xkjp ∣∣∣ G is meant to denote the partial derivative w.r.t Xkjp but while keeping G fixed i.e. treated as a constant or being “conditioned upon” (indeed both G and Ŷ depend on Xkjp but the definition keeps G fixed while taking derivatives). Similarly, in Lemma 1 part 1, αk is fixed (treated as a constant) in the derivative in the definition of ∂L/∂Ek and in part 2, αk+1 is fixed (treated as a constant) in the derivative in the definition of αk.
Recalling some Notation for sake of completeness: We recall that xki is the k-th layer embedding of node i. yi is the C-dimensional ground-truth label vector for node i and ŷi is the C-dimensional predicted score vector for node i. K is total number of layers in the GCN. Li denotes the loss on the i-th node and L denotes the total loss summed over all nodes. dk is the embedding dimensionality after the k-th layer.
Proof of Lemma 1.1 We analyze two cases
Case 1 (k = K i.e. final layer): Recall that the predictions for node i are obtained as ŷi = (WK+1)>xKi . Thus we have
∂L ∂WK+1 = ∑ i∈V ∑ c∈[C] ∂Li ∂ŷic ∂ŷic ∂WK+1 .
Now, ∂Li∂ŷic = gic by definition. If we let w K+1 c denote the c-th column of the dk × C matrix WK+1, then it is clear that ŷic depends only on wK+1c and x K i . Thus, we have
∂ŷic ∂WK+1 = ∂ŷic ∂wK+1c · e>c = (xKi )>e>c ,
where ec is the c-th canonical vector in C-dimensions with 1 at the c-th coordinate and 0 everywhere else. This gives us
∂L ∂WK+1 = ∑ i∈V (xKi ) > ∑ c∈[C] gic · e>c = (XK)>G
Case 2 (k < K i.e. intermediate layers): We recall that Xk stacks all k-th layer embeddings as an N × dk matrix and Xk = f(Xk−1;Ek) where Ek denotes the parameters (weights, offsets, scales) of the k-th layer. Thus we have
∂L ∂Ek = ∑ i∈V ∑ c∈[C] ∂L ∂ŷic ∂ŷic ∂Ek .
As before, ∂L∂ŷic = gic by definition and we have
∂ŷic ∂Ek = ∑ j∈V dk∑ p=1 ∂ŷic ∂Xkjp ∂Xkjp ∂Ek
This gives us
∂L ∂Ek = ∑ i∈V ∑ c∈[C] gic · ∑ j∈V dk∑ p=1 ∂ŷic ∂Xkjp ∂Xkjp ∂Ek = ∑ j∈V dk∑ p=1 αkjp · ∂Xkjp ∂Ek = ∂(αk Xk) ∂Ek ∣∣∣∣ αk ,
where we used the definition of αkjp in the second step and used a “conditional” notation to get the third step. We reiterate that ∂(α k Xk) ∂Ek ∣∣∣ αk implies that αk is fixed (treated as a constant) while taking the derivative. This “conditioning” is critical since αk also depends on Ek. This concludes the proof.
Proof of Lemma 1.2: We consider two cases yet again and use Definition 1 that tells us that αkjp = ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xkjp
Case 1 (k = K): Since ŷi = (WK+1)>xKi , we know that ŷic depends only on xKi and wK+1c where as before, wK+1c is the c-th column of the matrix W
K+1. This gives us ∂ŷic ∂XKjp = 0 if i 6= j and ∂ŷic ∂XKjp = wKpc if i = j where w K pc is the (p, c)-th entry of the matrix W K+1 (or in other words, the p-th coordinate of the vector wK+1c ). This tells us that
αKjp = ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂XKjp = ∑ c∈[C] gjc · wKpc,
which gives us αK = G(WK+1)>.
Case 2 (k < K): By Definition 1 we have
αkjp = ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xkjp = ∑ i∈V ∑ c∈[C] gic · ∑ l∈V dk+1∑ q=1 ∂ŷic ∂Xk+1lq ∂Xk+1lq ∂Xkjp
Rearranging the terms gives us
αkjp = ∑ l∈V dk+1∑ q=1 ∑ i∈V ∑ c∈[C] gic · ∂ŷic ∂Xk+1lq  · ∂Xk+1lq ∂Xkjp = ∑ l∈V dk+1∑ q=1 αk+1lq · ∂Xk+1lq ∂Xkjp ,
where we simply used Definition 1 in the second step. However, the resulting term is simply ∂(αk+1 Xk+1)
∂Xkjp
∣∣∣ αk+1 which conditions on, or treats as a constant, the term αk+1 according to our
notation convention. This finishes the proof of part 2.

heading: C.2 STATEMENT OF CONVERGENCE GUARANTEE
text: The rate for full-batch updates, as derived below, is O ( 1
T 2 3
) . This fast rate offered by full-batch
updates is asymptotically superior to the O (
1√ T
) rate offered by mini-batch SGD updates. This is
due to the additional variance due to mini-batch construction that mini-batch SGD variants have to incur.
Theorem 4 (IGLU Convergence (Final)). Suppose the task loss function L has H-smooth and an architecture that offers bounded gradients and Lipschitz gradients as quantified below, then if IGLU in its inverted variant (Algorithm 2) is executed with step length η and a staleness count of τ updates per layer as in steps 7, 10 in Algorithm 2, then within T iterations, we must have
1. ‖∇L‖22 ≤ O ( 1/T 2 3 ) if model update steps are carried out on the entire graph in a
full-batch with step length η = O ( 1/T 1 3 ) and τ = O (1).
2. ‖∇L‖22 ≤ O ( 1/ √ T ) if model update steps are carried out using mini-batch SGD with
step length η = O ( 1/ √ T ) and τ = O ( T 1 4 ) .
It is curious that the above result predicts that when using mini-batch SGD, a non-trivial amount of staleness (τ = O ( T 1 4 ) as per the above result) may be optimal and premature refreshing of
embeddings/incomplete gradients may be suboptimal as was also seen in experiments reported in Table 2 and Figure 4. Our overall proof strategy is the following
1. Step 1: Analyze how lazy updates in IGLU affect model gradients 2. Step 2: Bound the bias in model gradients in terms of staleness due to the lazy updates 3. Step 3: Using various properties such as smoothness and boundedness of gradients, obtain
an upper bound for the bias in the gradients in terms of number of iterations since last update
4. Step 4: Use the above to establish the convergence guarantee
We will show the results for the variant of IGLU that uses the inverted order of updates as given in Algorithm 2. A similar proof technique will also work for the variant that uses the backprop order of updates. Also, to avoid clutter, we will from hereon assume that a normalized total loss function is used for training i.e. L = 1N ∑ i∈V `i where N = |V| is the number of nodes in the training graph.

heading: C.3 STEP 1: PARTIAL STALENESS AND ITS EFFECT ON MODEL GRADIENTS
text: A peculiarity of the inverted order of updates is that the embeddings Xk, k ∈ [K] are never stale in this variant. To see this, we use a simple inductive argument. The base case of X0 is obvious – it is never stale since it is never meant to be updated. For the inductive case, notice how, the moment any parameter Ek is updated in step 7 of Algorithm 2 (whether by mini-batch SGD or by full-batch GD), immediately thereafter in step 8 of the algorithm, Xk is updated using the current value of Xk−1 and Ek. Since by induction Xk−1 never has a stale value, this shows that Xk is never stale either, completing the inductive argument.
This has an interesting consequence: by Lemma 1, we have ∂L ∂Ek = 1N ∂(αk Xk) ∂Ek ∣∣∣ αk (notice the additional 1/N term since we are using a normalized total loss function now). However, as ∂(α k Xk) ∂Ek ∣∣∣ αk is completely defined given Ek,αk and Xk−1 and by the above argument, Xk−1 never has a stale value. This shows that the only source of staleness in ∂L
∂Ek is the staleness in values of the incomplete
task gradient αk. Similarly, it is easy to see that the only source of staleness in ∂L ∂WK+1 = (XK)>G is the staleness in G.
The above argument is easily mirrored for the backprop order of updates where an inductive argument similar to the one used to argue above that Xk values are never stale in the inverted update variant, would show that the incomplete task gradient αk values are never stale in the backprop variant and the only source of staleness in ∂L
∂Ek would then be the staleness in the Xk values.

heading: C.4 STEP 2: RELATING THE BIAS IN MODEL GRADIENTS TO STALENESS
text: The above argument allows us to bound the bias in model gradients as a result of lazy updates. To avoid clutter, we will present the arguments with respect to the Ek parameter. Similar arguments would hold for the WK+1 parameter as well. Let α̃k,αk denote the stale and actual values of the incomplete task gradient relevant to Ek. Let ∂̃L
∂Ek = ∂(α̃ k Xk) ∂Ek ∣∣∣ α̃k be the stale gradients used by
IGLU in its inverted variant to update Ek and similarly let ∂L ∂Ek = ∂(α k Xk) ∂Ek ∣∣∣ αk
be the true gradient that could have been used had there been no staleness.
We will abuse notation and let the vectorized forms of these incomplete task gradients be denoted by the same symbols i.e. we stretch the matrix αk ∈ RN×dk into a long vector denoted also by αk ∈ RN ·dk . Let dim(Ek) denote the number of dimensions in the model parameter Ek (recall that Ek can be a stand-in for weight matrices, layer norm parameters etc used in layer k of the GCN).
Similarly, we let Zkjp ∈ Rdim(E k), j ∈ [N ], p ∈ [dk] denote the vectorized form of the gradient ∂Xkip ∂Ek and let Zk ∈ RN ·dk×dim(Ek) denote the matrix with all these vectors Zkjp stacked up. As per the above notation, it is easy to see that the vectorized form of the model gradient is given by
∂L ∂Ek = (Zk)>αk N ∈ Rdim(E k)
The 1/N term appears since we are using a normalized total loss function. This also tells us that∥∥∥∥∥ ∂̃L∂Ek − ∂L∂Ek ∥∥∥∥∥
2
=
√ (α̃k −αk)>(Zk(Zk)>)(α̃k −αk)
N ≤
∥∥∥α̃k −αk∥∥∥ 2 · σmax(Zk)
N , (2)
where σmax(Zk) is the largest singular value of the matrix Zk.

heading: C.5 STEP 3: SMOOTHNESS AND BOUNDED GRADIENTS TO BOUND GRADIENT BIAS
text: The above discussion shows how to bound the bias in gradients in terms of staleness in the incomplete task gradients. However, to utilize this relation, we assume that the loss function L is H-smooth which is a standard assumption in literature. We will also assume that the network offers bounded gradients. Specifically, for all values of model parameters Ek, k ∈ [K],WK+1 we have
‖G‖2 , ∥∥∥∥∂Xkip∂Ek ∥∥∥∥
2
, ∥∥ ∂L ∂Ek ∥∥ 2 , ∥∥∥∂Xk+1∂Xk ∥∥∥2 ,∥∥αk∥∥2 ≤ B. For sake of simplicity, we will also assume the
same bound on parameters e.g. ∥∥WK∥∥
2 ≤ B. Assuming bounded gradients and bounded parameters
is also standard in literature. However, whereas works such as Chen et al. (2018a) assume bounds on the sup-norm i.e. L∞ norm of the gradients, our proofs only require an L2 norm bound.
We will now show that if the model parameters Ẽk, k ∈ [K], W̃K+1 undergo gradient updates to their new valuesEk, k ∈ [K],WK+1 and the amount of travel is bounded by r > 0 i.e. ∥∥∥Ẽk − Ek∥∥∥ 2 ≤ r,
then we have ∥∥∥α̃k −αk∥∥∥
2 ≤ Ik · r where α̃k,αk are the incomplete task gradients corresponding to
respectively old and new model parameter values and Ik depends on various quantities such as the smoothness parameter and the number of layers in the network.
Lemma 1 (part 2) tells us that for the final layer, we have αK = G(WK+1)> as well as for any k < K, we have αk = ∂(α
k+1 Xk+1) ∂Xk ∣∣∣ αk+1 . We analyze this using an inductive argument.
1. Case 1: k = K (Base Case): In this case we have
α̃K −αK = G̃(W̃K+1)> −G(WK+1)>
= (G̃−G)(W̃K+1)> + G(W̃K+1 −WK+1)> Now, the travel condition tells us that ∥∥∥W̃K+1 −WK+1∥∥∥
2 ≤ r, boundedness tells us that ‖G‖2 , ∥∥∥W̃K+1∥∥∥
2 ≤ B. Also, since the loss function is H-smooth, the task gradients are H-Lipschitz which implies along with the travel condition that ∥∥∥G̃−G∥∥∥
2 ≤ H · r. Put
together this tells us that ∥∥∥α̃K −αK∥∥∥ 2 ≤ (H + 1)B · r,
telling us that IK ≤ (H + 1)B. 2. Case 2: k < K (Inductive Case): In this case, let X̃k and Xk denote the embeddings with
respect to the old and new parameters respectively. Then we have
α̃k −αk = ∂(α̃ k+1 X̃k+1) ∂X̃k ∣∣∣∣∣ α̃k+1 − ∂(α k+1 Xk+1) ∂Xk ∣∣∣∣ αk+1
= ∂(α̃k+1 X̃k+1)
∂X̃k
∣∣∣∣∣ α̃k+1 − ∂(α k+1 X̃k+1) ∂X̃k ∣∣∣∣∣ αk+1︸ ︷︷ ︸
(P )
+ ∂(αk+1 X̃k+1)
∂X̃k
∣∣∣∣∣ αk+1 − ∂(α k+1 Xk+1) ∂Xk ∣∣∣∣ αk+1︸ ︷︷ ︸
(Q) By induction we have ∥∥∥α̃k+1 −αk+1∥∥∥
2 ≤ Ik+1 · r and bounded gradients tells us∥∥∥∂X̃k+1
∂X̃k ∥∥∥ 2 ≤ B giving us ‖(P )‖2 ≤ Ik+1B · r. To analyze the term ‖(Q)‖2, recall that we have Xk+1 = f(Xk;Ek+1) and since the overall task loss is H smooth, so must be the function f . Abusing notation to let H denote the Lipschitz constant of the network gives us ∥∥∥X̃k −Xk∥∥∥
2 ≤ H · r. Then we have
∂(X̃k+1)
∂X̃k − ∂X
k+1
∂Xk = f ′(X̃k; Ẽk+1)− f ′(Xk;Ek+1)
= f ′(X̃k; Ẽk+1)− f ′(Xk; Ẽk+1)︸ ︷︷ ︸ (M) + f ′(Xk; Ẽk+1)− f ′(Xk;Ek+1)︸ ︷︷ ︸ (N)
Applying smoothness, we have ‖(M)‖2 ≤ H2 · r as well as ‖(N)‖2 ≤ H · r. Together with bounded gradients that gives us ∥∥αk+1∥∥ 2 ≤ B, we have ‖(Q)‖2 ≤ BH(H + 1) · r.
Together, we have ∥∥∥α̃k −αk∥∥∥ 2 ≤ B(H(H + 1) + Ik+1) · r, telling us that Ik ≤ B(H(H + 1) + Ik+1).
The above tells us that in Algorithm 2, suppose the parameter update steps i.e. step 7 and step 10 are executed by effecting τ mini-batch SGD steps or else τ full-batch GD steps each time, with
step length η, then the amount of travel in any model parameter is bounded above by τηB i.e.∥∥∥Ẽk − Ek∥∥∥ 2 ≤ τηB and so on. Now, incomplete task gradients αk are updated only after model parameters for all layers have been updated once and the algorithm loops back to step 6. Thus, Lipschitzness of the gradients tells us that the staleness in the incomplete task gradients, for any k ∈ [K], is upper bounded by ∥∥∥α̃k −αk∥∥∥
2 ≤ τη · IB, where we take I = maxk≤K Ik for sake of simplicity. Now, since ∥∥∥∥∂Xkip∂Ek ∥∥∥∥
2
≤ B as gradients are
bounded, we have σmax(Zk) ≤ N · dkB. Then, combining with the result in Equation (2) gives us∥∥∥∥∥ ∂̃L∂Ek − ∂L∂Ek ∥∥∥∥∥
2
≤ τη · IB2dk
Taking in contributions of gradients of all layers and the final classifier layer gives us the bias in the total gradient using triangle inequality as∥∥∥∇̃L −∇L∥∥∥
2 ≤ ∑ k∈[K] ∥∥∥∥∥ ∂̃L∂Ek − ∂L∂Ek ∥∥∥∥∥ 2 + ∥∥∥∥∥ ∂̃L∂WK+1 − ∂L∂WK+1 ∥∥∥∥∥ 2 ≤ τη · IB2dmax(K + 1),
where dmax = maxk∈K dk is the maximum embedding dimensionality of any layer.

heading: C.6 STEP 4 (I): CONVERGENCE GUARANTEE (MINI-BATCH SGD)
text: Let us analyze convergence in the case when updates are made using mini-batch SGD in steps 7, 10 of Algorithm 2. The discussion above establishes an upper bound on the absolute bias in the gradients. However, our proofs later require a relative bound which we tackle now. Let us decide to set the step length to η = 1
C √ T for some constant C > 0 that will be decided later and also set some value φ < 1. Then two cases arise
1. Case 1: The relative gradient bias is too large i.e. τη · IB2dmax(K + 1) > φ · ‖∇L‖2. In this case we are actually done since we get
‖∇L‖2 ≤ τ · IB2dmax(K + 1)
φ · C √ T
,
i.e. we are already at an approximate first-order stationary point.
2. Case 2: The relative gradient bias is small i.e. τη · IB2dmax(K + 1) ≤ φ · ‖∇L‖2. In this case we satisfy the relative bias bound required by Lemma 5 (part 2) with δ = φ.
This shows that either Case 1 happens in which case we are done or else Case 2 keeps applying which means that Lemma 5 (part 2) keeps getting its prerequisites satisfied. If Case 1 does not happen for T steps, then Lemma 5 (part 2) assures us that we will arrive at a point where
E [ ‖∇L‖22 ] ≤ 2C(L
0 − L∗) +Hσ2/C (1− φ) √ T
where L0,L∗ are respectively the initial and optimal values of the loss function which we recall is H-smooth and σ2 is the variance due to mini-batch creation and we set η = 1
C √ T . Setting C = √ Hσ2
2(L0−L∗) tells us that within T steps, either we will achieve
‖∇L‖22 ≤ 2τ2 · I2B4d2max(K + 1)2(L0 − L∗)
φ2Hσ2T ,
or else we will achieve
E [ ‖∇L‖22 ] ≤ 2σ √
2(L0 − L∗)H (1− φ) √ T
Setting τ = T 1 4 · ( √ σ3H √ H
(L0−L∗) 1 4 IB2dmax(K+1)
) balances the two quantities in terms of their dependence
on T (module absolute constants such as φ). Note that as expected, as the quantities I,B, dmax,K increase, the above limit on τ goes down i.e. we are able to perform fewer and fewer updates to the model parameters before a refresh is required. This concludes the proof of Theorem 4 for the second case.

heading: C.7 STEP 4 (II): CONVERGENCE GUARANTEE (FULL-BATCH GD)
text: In this case we similarly have either the relative gradient bias to be too large in which case we get
‖∇L‖2 ≤ τη · IB2dmax(K + 1)
φ ,
or else we satisfy the relative bias bound required by Lemma 5 (part 1) with δ = φ. This shows that either Case 1 happens in which case we are done or else Case 2 keeps applying which means that Lemma 5 (part 1) keeps getting its prerequisites satisfied. If Case 1 does not happen for T steps, then Lemma 5 (part 1) assures us that we will arrive at a point where
‖∇L‖22 ≤ 2(L0 − L∗) η(1− φ)T
In this case, for a given value of τ , setting η = (
2(L0−L∗) (1−φ)τ2I2B4d2max(K+1)2T
) 1 3
gives us that within T iterations, we must achieve
‖∇L‖22 ≤ ( 2τ(L0 − L∗)IB2dmax(K + 1) (1− φ) ) 2 3 · 1 T 2 3
In this case, it is prudent to set τ = O (1) so as to not deteriorate the convergence rate. This concludes the proof of Theorem 4 for the first case.

heading: C.8 GENERIC CONVERGENCE RESULTS
text: Lemma 5 (First-order Stationarity with a Smooth Objective). Let f : Θ → R be an H-smooth objective over model parameters θ ∈ Θ that is being optimized using a gradient oracle and the following update for step length η:
θt+1 = θt − η · gt
Let θ∗ be an optimal point i.e. θ∗ ∈ arg minθ∈Θ f(θ). Then, the following results hold depending on the nature of the gradient oracle:
1. If a non-stochastic gradient oracle with bounded bias is used i.e. for some δ ∈ (0, 1), for all t, we have gt = ∇f(θt) + ∆t where ‖∆t‖2 ≤ δ · ‖∇f(θt)‖2 and if the step length satisfies η ≤ (1−δ)2H(1+δ2) , then for any T > 0, for some t ≤ T we must have∥∥∇f(θt)∥∥2
2 ≤ 2(f(θ 0)− f(θ∗)) η(1− δ)T
2. If a stochastic gradient oracle is used with bounded bias i.e. for some δ ∈ (0, 1), for all t, we have E [ |gtθt] = ∇f(θt) + ∆t where ‖∆t‖2 ≤ δ · ‖∇f(θt)‖2, as well as bounded variance i.e. for all t, we have E [ | ‖gt −∇f(θt)−∆t‖22θt ] ≤ σ2 and if the step length
satisfies η ≤ (1−δ)2H(1+δ2) , as well as η = 1 C √ T for some C > 0, then for any T > 4H 2(1+δ2)2
C2(1−δ)2 , for some t ≤ T we must have
E [∥∥∇f(θt)∥∥2
2
] ≤ 2C(f(θ
0)− f(θ∗)) +Hσ2/C (1− δ) √ T
Proof (of Lemma 5). To prove the first part, we notice that smoothness of the objective gives us f(θt+1) ≤ f(θt) + 〈 ∇f(θt), θt+1 − θt 〉 + H
2
∥∥θt+1 − θt∥∥2 2
Since we used the update θt+1 = θt − η · gt and we have gt = ∇f(θt) + ∆t, the above gives us
f(θt+1) ≤ f(θt)− η · 〈 ∇f(θt),gt 〉 + Hη2
2
∥∥gt∥∥2 2
= f(θt)− η · 〈 ∇f(θt),∇f(θt) + ∆t 〉 + Hη2
2
(∥∥∇f(θt) + ∆t∥∥2 2 ) = f(θt)− η · ∥∥∇f(θt)∥∥2 2 − η · 〈 ∇f(θt),∆t 〉 + Hη2
2
(∥∥∇f(θt) + ∆t∥∥2 2 ) Now, the Cauchy-Schwartz inequality along with the bound on the bias gives us−η · 〈∇f(θt),∆t〉 ≤ ηδ · ‖∇f(θt)‖22 as well as ‖∇f(θt) + ∆t‖ 2 2 ≤ 2(1 + δ2) · ‖∇f(θt)‖ 2 2. Using these gives us
f(θt+1) ≤ f(θt)− η(1− δ − ηH(1 + δ2)) · ∥∥∇f(θt)∥∥2
2
≤ f(θt)− η(1− δ) 2
· ∥∥∇f(θt)∥∥2
2 ,
since we chose η ≤ (1−δ)2H(1+δ2) . Reorganizing, taking a telescopic sum over all t, using f(θ T+1) ≥ f(θ∗) and making an averaging argument tells us that since we set , for any T > 0, it must be the case that for some t ≤ T , we have∥∥∇f(θt)∥∥2
2 ≤ 2(f(θ 0)− f(θ∗)) η(1− δ)T
This proves the first part. For the second part, we yet again invoke smoothness to get f(θt+1) ≤ f(θt) + 〈 ∇f(θt), θt+1 − θt 〉 + H
2
∥∥θt+1 − θt∥∥2 2
Since we used the update θt+1 = θt − η · gt, the above gives us f(θt+1) ≤ f(θt)− η · 〈 ∇f(θt),gt 〉 + Hη2
2
∥∥gt∥∥2 2
= f(θt)− η · 〈 ∇f(θt),gt 〉 + Hη2
2
(∥∥∇f(θt) + ∆t∥∥2 2 + ∥∥gt −∇f(θt)−∆t∥∥2 2 ) −Hη2 〈 ∇f(θt) + ∆t,gt −∇f(θt)−∆t
〉 Taking conditional expectations on both sides gives us
E [ | f(θt+1)θt ] ≤ f(θt)− η · ∥∥∇f(θt)∥∥2 2 − η · 〈 ∇f(θt),∆t 〉 + Hη2
2
(∥∥∇f(θt) + ∆t∥∥2 2 + σ2 )
Now, the Cauchy-Schwartz inequality along with the bound on the bias gives us−η · 〈∇f(θt),∆t〉 ≤ ηδ · ‖∇f(θt)‖22 as well as ‖∇f(θt) + ∆t‖ 2 2 ≤ 2(1 + δ2) · ‖∇f(θt)‖ 2 2. Using these and applying a total expectation gives us
E [ f(θt+1) ] ≤ E [ f(θt) ] − η(1− δ − ηH(1 + δ2)) · E [∥∥∇f(θt)∥∥2 2 ] + Hη2 2 · σ2
≤ E [ f(θt) ] − η(1− δ) 2 · E [∥∥∇f(θt)∥∥2 2 ] + Hη2 2 · σ2
where the second step follows since we set η ≤ (1−δ)2H(1+δ2) . Reorganizing, taking a telescopic sum over all t, using f(θT+1) ≥ f(θ∗) and making an averaging argument tells us that for any T > 0, it must be the case that for some t ≤ T , we have
E [∥∥∇f(θt)∥∥2
2
] ≤ 2(f(θ
0)− f(θ∗)) η(1− δ)T + Hησ2
1− δ However, since we also took care to set η = 1
C √ T , we get E [∥∥∇f(θt)∥∥2
2
] ≤ 2C(f(θ
0)− f(θ∗)) +Hσ2/C (1− δ) √ T
which proves the second part and finishes the proof.

Review summaries:
This paper studies tries to tackle the scalability challenge of training GNNs on large graphs. The authors propose IGLU, an architecture-agnostic method. IGLU caches intermediate computations and uses a lazy update strategy. Convergence analysis on IGLU is provided and empirical results show that IGLU has better performance than baselines such as GraphSAINT, Cluster-GCN, VR-GCN. 

The work proposes IGLU, an algorithm to scale-up GNNs using stale computations instead of traditional neighborhood sampling. IGLU has bounded bias if the loss and activation functions are smooth. Results on large-scale benchmarks show IGLU achieves SOTA and scales better than previous methods.

This paper introduces a new method, IGLU, that caches intermediate computations at various GCN layers. This enables IGLU to perform lazy updates that do not require updating a large number of node embeddings during descent and offers much faster convergence without significantly biasing the gradients.  Overall, this paper represents a novel solution towards efficient GNN training.