Year: 2021

Abstract:
To accelerate the training of graph convolutional networks (GCNs), many sampling-based methods have been developed for approximating the embedding aggregation. Among them, a layer-wise approach recursively performs importance sampling to select neighbors jointly for existing nodes in each layer. This paper revisits the approach from a matrix approximation perspective. We identify two issues in the existing layer-wise sampling methods: sub-optimal sampling probabilities and the approximation bias induced by sampling without replacement. We thus propose remedies to address these issues. The improvements are demonstrated by extensive analyses and experiments on common benchmarks.

text: To accelerate the training of graph convolutional networks (GCNs), many sampling-based methods have been developed for approximating the embedding aggregation. Among them, a layer-wise approach recursively performs importance sampling to select neighbors jointly for existing nodes in each layer. This paper revisits the approach from a matrix approximation perspective. We identify two issues in the existing layer-wise sampling methods: sub-optimal sampling probabilities and the approximation bias induced by sampling without replacement. We thus propose remedies to address these issues. The improvements are demonstrated by extensive analyses and experiments on common benchmarks.

heading: 1 INTRODUCTION
text: Graph Convolutional Networks (Kipf & Welling, 2017) are popular methods for learning the representation of nodes. However, it is computationally challenging to train a GCN over large-scale graphs due to the inter-dependence of nodes in a graph. In the mini-batch training for an L-layer GCN, the computation of embeddings involves not only the batch nodes but also batch nodes’ Lhop neighbors, which is known as the phenomenon of “neighbor explosion” (Zeng et al., 2019) or “neighbor expansion” (Chen et al., 2018a; Huang et al., 2018). To alleviate such a computation issue for large graphs, sampling-based methods are proposed to accelerate the training and reduce the memory cost. These approaches can be categorized as node-wise sampling approaches (Hamilton et al., 2017; Chen et al., 2018a), subgraph sampling approaches (Zeng et al., 2019; Chiang et al., 2019; Cong et al., 2020), and layer-wise sampling approaches (Chen et al., 2018b; Huang et al., 2018; Zou et al., 2019). We focus on layer-wise sampling in this work, which enjoys the efficiency and variance reduction by sampling columns of renormalized Laplacian matrix in each layer.
This paper is a study of the existing sampling schemes in layer-wise sampling methods. We identify two potential drawbacks in the common practice for layer-wise sampling (Chen et al., 2018b; Zou et al., 2019). First, the current sampling probabilities used are sub-optimal since a core assumption in FastGCN and LADIES does not hold in many common graph benchmarks, such as Reddit (Hamilton et al., 2017) and OGB(Hu et al., 2020). Secondly, the previous implementations of the layer-wise sampling methods slightly deviate from their theoretical results, and introduce bias in the estimation due to the usage of sampling without replacement. Realizing the two issues, we accordingly propose the remedies with new sampling probabilities and a debiasing algorithm. The improvements of the proposed methods are demonstrated by extensive experiments on evaluating the matrix approximation error and prediction accuracy on large-scale benchmarks along with some theoretical analyses.
To the best of our knowledge, our result is the first to recognize and resolve the issues with the default assumption and the practical implementation of layer-wise sampling methods for GCN. Once these sub-optimal practices are addressed, we observe that the GCN models consistently converge faster in training and usually enjoy a higher prediction accuracy. We believe the proposed methods can more generally improve the training for GCN as well, e.g., the same strategy can allow node-wise sampling methods to adopt sampling without replacement and further improve the approximation accuracy. Moreover, our discussion on the bias induced by sampling without replacement is not limited to GCN, and the debiasing algorithm we develop can contribute to other sampling-based machine learning models beyond layer-wise sampling.

heading: 1.1 BACKGROUND AND RELATED WORK
text: GCN Graph Convolutional Networks (GCNs, Kipf & Welling (2017)) effectively incorporate the technique of convolution filter into the graph domain (Wu et al., 2020; Bronstein et al., 2017). Viewed as an approximation for the spectral graph convolutions (Bruna et al., 2014; Defferrard et al., 2016), GCN has achieved great success in learning tasks such as node classification and link prediction, with applications ranging from recommender systems (Ying et al., 2018), traffic prediction (Cui et al., 2019; Rahimi et al., 2018), and knowledge graphs (Schlichtkrull et al., 2018).
Sampling Based GCN Training To name a few of sampling schemes, GraphSAGE (Hamilton et al., 2017) first introduces the “node-wise” neighbor sampling scheme, where a fixed number of neighbors are uniformly and independently sampled for each node in every layer. To reduce variance in node-wise sampling training, VR-GCN (Chen et al., 2018a) applies a control variate approach with historical activation. Instead of sampling for each node, “layer-wise” sampling is a more efficient approach: joint sampling scheme for all the existing nodes in each layer so that these nodes can share the sampled neighboring node. FastGCN (Chen et al., 2018b) first introduces this scheme with importance sampling. AS-GCN (Huang et al., 2018) proposes an alternative sampling probability for layer-wise sampling by approximating the hidden layer in the sampling procedure. Then Zou et al. (2019) propose a layer-dependent importance sampling scheme (LADIES) to further reduce the variance in training. This alleviates the issue of empty rows in sampled adjacency matrix for FastGCN. For the “subgraph” approach, ClusterGCN (Chiang et al., 2019) samples a dense subgraph associated with the batch nodes by graph clustering algorithm; GraphSAINT (Zeng et al., 2019) introduces normalization and variance reduction in subgraph sampling.

heading: 2 NOTATIONS AND PRELIMINARIES

heading: 2.1 GRAPH CONVOLUTIONAL NETWORKS
text: The GCN architecture for semi-supervised node classification is introduced by Kipf & Welling (2017). Suppose we have an undirected graph G = (V, E), where V is the set of n nodes and E is the set of E edges. Denote node i in V as vi, where i ∈ [n] is the index of nodes in the graph and [n] denotes the set {1, 2, ..., n}. Each node vi ∈ V is associated with a feature vector xi ∈ Rp and a label vector yi ∈ Rq . Though we can observe the feature of every node in V and every edge in E , i.e. the n × n adjacency matrix A, we are only able to observe the label of partial nodes Vtrain, satisfying Vtrain ⊂ V . Thus, we need to predict the labels for the rest nodes in V\Vtrain and it becomes a semi-supervised learning task. A graph convolution layer is defined as:
Z(l+1) = PH(l)W (l), H(l) = σ(Z(l)), (1)
where σ is an activation function and P is obtained from applying normalization to the graph adjacency matrix A; H(l) is the embedding matrix of the graph nodes in the l-th layer, and W (l) is the parameter matrix of the same layer. In particular, H(0) is the n × p feature matrix. For mini-batch training, the training loss for an L-layer GCN is defined as 1|Vbatch| ∑ vi∈Vbatch `(yi, z (L) i ), where ` is the loss function, batch nodes Vbatch is a subset of Vtrain at each iteration. z(L)i is the i-th row in Z(L), | · | denotes the cardinality of a set.
In this paper, we set P = D̃−1/2(A + I)D̃−1/2, where D̃ is a diagonal matrix with Dii = 1 + ∑ iAij . The matrix P is constructed as a renormalized Laplacian matrix to help alleviate overfitting and exploding/vanishing gradients issues (Kipf & Welling, 2017), which is used by Kipf & Welling (2017); Chen et al. (2018a); Cong et al. (2020).

heading: 2.2 LAYER-WISE SAMPLING
text: To address “neighbor explosion” issue for graph neural networks, sampling methods are integrated into the stochastic training. Motivated by the idea to approximate the matrix PH(l) in (1), FastGCN (Chen et al., 2018b) applies an importance-sampling-based strategy. Instead of individually sampling neighbors for each node in the l-th layer, they sample a set of s neighbors S(l) from V with importance sampling probability pi, where pi ∝ ∑n j=1P 2 ji and ∑ i pi = 1. For the (l− 1)-th layer,
they naturally set V(l−1) = S(l). LADIES (Zou et al., 2019) improves the importance sampling probability pi as
p (l) i ∝ ∑ vj∈N (l)P 2 ji,∀i ∈ [n] (2)
where N (l) = ∪vi∈V(l)N (vi) and ∑ j p (l) j = 1. In this case, S(l) the nodes sampled for the l-th layer are guaranteed to be within the neighborhood of V(l). The whole procedure can be concluded by a diagonal matrix S(l) ∈ Rn×n and a row selection matrixQ(l) ∈ Rsl×n, which are defined as
Q (l) k,j =
{ 1, j = i
(l) k
0, else , S
(l) j,j =
{ (slp (l)
i (l) k
)−1, j = i (l) k
0, else,
where {i(l)k } sl k=1 are the indices of rows selected in the l-th layer. The forward propagation with layer-wise sampling can thus be equivalently represented as Z̃(l+1) = Q(l+1)PS(l)H(l)W (l),H(l) = (Q(l))Tσ(Z̃(l)), where Z̃(l+1) is the approximation of the embedding matrix for layer l.

heading: 3 EXPERIMENTAL SETUP
text: In advance of the formal introduction to the the issues and the corresponding remedies in Section 4 and Section 5, we state the basic setups of the main experiments and datasets as they appear multiple times across the paper. Details about GCN model training are deferred to the according sections.
Main experiments. To study the influence of the aforementioned issues we evaluate the matrix approximation error (c.f. Figure 1) of different methods and consider it as a new metric to reflect the performance of the sampling strategy on approximating the original mini-batch training in one-step propagation. Since the updates of parameters in the training is not involved in the simple metric above, in Section 6 we further evaluate the prediction accuracy on testing sets of both intermediate models during training and final outputs, using the metrics in Table 2.
Benchmarks. Empirical experiments are conducted on 5 datasets (see details at Table 2 in Appendix B): Reddit (Hamilton et al., 2017), ogbn-arxiv, ogbn-proteins, ogbn-mag and ogbn-products (Hu et al., 2020). Reddit is a traditional large graph dataset used by Chen et al. (2018b); Zou et al. (2019); Chen et al. (2018a); Cong et al. (2020); Zeng et al. (2019). Ogbn-arxiv, ogbn-proteins and ogbn-products are Open Graph Benchmarks (OGB) proposed by Hu et al. (2020). Compared to traditional datasets, our selected OGB data have larger volume (up to million-node scale) with more challenging data split. The metrics in Table 2 follow the choices of recent works and the recommendation by (Hu et al., 2020).

heading: 4 RECONSIDER IMPORTANCE SAMPLING PROBABILITIES
text: The efficiency of layer-wise sampling comes from sampling, and the choice of sampling probabilities impacts the prediction accuracy of GCNs. To minimize the variance (for the sake of notational
brevity, from now on we omit the superscript (l) when the objects are from the same layer)
E ‖QPSHW −QPHW ‖2F ,
Zou et al. (2019) show that the optimal sampling probability pi for node i satisfies (also see Appendix E for a derivation from a perspective of approximate matrix multiplication)
pi ∝ ‖(HW )[i]‖ · ‖QP [i]‖,
where for a matrixA,A[i] andA[i] respectively represent the i-th row / column of matrixA.

heading: 4.1 CURRENT APPROXIMATION STRATEGIES AND THEIR LIMITATIONS
text: The optimal sampling probabilities discussed above are usually unavailable during the mini-batch training. W andH keep changing in the training and even worse we only have access to part of the matrix HW in each batch. To approximate the optimal sampling probability, previous works develop two different strategies, 1) approximating the hidden activationHW , or 2) sampling without the information fromHW .
A representatives of the former strategy is AS-GCN (Huang et al., 2018), which linearly approximate HW by node features. This strategy improves the accuracy of efficient GCN, while the estimation of hidden activation causes considerable overall training time, which in practice is “even longer than vanilla GCN” (Zeng et al., 2019). Another issue of this strategy is that the approximation must be updated during the training procedure. This dependence between sampling and training renders it impossible to save training time by preparing sampling results in advance of the training.
Instead, the second strategy, the theme of our paper, allows the decoupling of sampling and training, and has been adopted by FastGCN (Chen et al., 2018b) and LADIES (Zou et al., 2019). To proceed without the information from HW , FastGCN (resp. LADIES) assumes ‖(HW )[i]‖ ∝ ‖P [i]‖ (resp. ‖QP [i]‖), and sets their sampling probabilities as pi ∝ ‖P [i]‖2 (resp. ‖QP [i]‖2), ∀i ∈ [n]. However, we find this assumption is too strong. To test it, we conduct linear regression y(l) ∼ β0 + β1x for different layers separately (l = 1, 2, 3), where each y(l) is the `2 norm of one row in H(l)W (l) and each x is the norm of corresponding column in P . The result of ogbn-arxiv is presented in Figure 2. For layer 1 to layer 3, the intercepts (with stand error) of the regression lines are 2.364(±0.002), 4.014(±0.008) and 19.722(±0.024) respectively, which are much larger than zero with the statistical significance; however, the according slopes β1’s are at a smaller scale: −0.154(±0.005), 4.468(±0.021) and 2.158(±0.069). Note that β1 for layer 1 is even negative with statistical significance. Such evidence shows that the relation between ‖(HW )[i]‖ and ‖P [i]‖ is pretty unstable and might be even negatively correlated. This violates the assumption that ‖(HW )[i]‖ ∝ ‖P [i]‖. Similar patterns are also observed on the other datasets (See Appendix C.3).

heading: 4.2 PROPOSED SAMPLING PROBABILITIES
text: As analyzed above, the sampling probabilities in LADIES are sub-optimal due to the failure of their assumption. To address this issue, we instead admit that we have no prior knowledge of HW , and
tend to assume a uniform distribution of ‖(HW )[i]‖’s. With this belief we propose the following sampling probabilities:
pi ∝ ‖QP [i]‖, ∀i ∈ [n]. (3) Compared to the sampling probabilities of LADIES in Equation (2), our proposed sampling probabilities pi’s are more conservative. From a matrix approximation perspective, we rewrite the target matrix product as QPIHW , and only aim to approximate the known part QPI . It turns out that assuming the uniform distribution of the norms of rows inHW can help improve both the variance of the matrix approximation and the prediction accuracy of GCNs. We empirically justify our sampling probabilities in the next section and Section 6. The theoretical analysis of the new probability is presented in Section 4.4.

heading: 4.3 EVALUATION ON MATRIX APPROXIMATION ERROR
text: To further justify our proposed sampling probability, we consider the following 1-layer matrix approximation error, which evaluates the propagation approximation to the embedding aggregation of full-batch.
‖Z̃(1)batch − Z̃ (1) sampling‖F = ‖QbatchPH (0)W (0) −QbatchPSH(0)W (0)‖F ,
where Z̃(1)batch and Z̃ (1) sampling are the embedding at the bottom layer for the whole batch (using all avaiable neighbors without sampling) and for a certain sampling method; S is the sampling matrix; Qbatch’s 0, 1 diagonal entries indicate if a node is in the batch; ‖ · ‖F denotes the Frobenius norm of a matrix. The experiments are repeated 200 times, in which we regenerate the batch nodes (shared by all sampling methods) and the sampling matrix for each method. The batch size is fixed as 512 and the numbers of sampled neighbors are 256, 512, 768, 1024, 2536, and 2048. W (0) is fixed and inherited from the trained model reported in Section 6.
It is not surprising to observe that in Figure 1, the result of our proposed sampling probabilities (denoted as “LAIDES+flat”, blue solid line) is consistently better than that of the original LADIES method (black dashed line) and of FastGCN (black solid line) on every dataset. This supports the new sampling probability in Equation (3). The discussion of the other (debiased) methods in this figure is deferred to Section 5.

heading: 4.4 THEORETICAL ANALYSIS OF THE NEW PROBABILITIES
text: Whether our choice of probabilities can outperform the previous ones depends on the distribution of the norms of rows in HW . When ‖HW(i)‖ is not proportional to the corresponding `2 norm of column (QP )(i), our proposed probabilities can benefit the approximate matrix multiplication task more than the ones assuming a relation of proportionality. We find the common long-tail distribution of numbers suffices to exert the strengths of the new probabilities, which can be concluded as the following assumption: Assumption 1. To simplify the notation we denote B := QP and C := HW , where P is an n-by-n matrix as defined above. Let m be the number of non-zero columns in B, and define C1 :=
‖B‖2F /m ( ∑n i=1 ‖B[i]‖/m) 2 ≥ 1. There also exists a constant C2 ≥ 1 such that 1C2 ‖C‖ 2 F /n ≤ ‖C[i]‖2 ≤ C2‖C‖2F /n. Assume C1/C22 ≥ 1.
With the assumption above, we show the variance of the approximation with our proposed probabilities is smaller than the variance of LADIES by the following lemma. Lemma 1. We denote the sampling matrix with our probabilities in Equation (3) as S1, and denote the sampling matrix with probabilities of LADIES in Equation (2) as S0. If Assumption 1 holds, then we have
E ‖BS1C −BC‖2F ≤ E ‖BS0C −BC‖2F .
The proof is provided in Appendix F.
Remark. Assumption 1 is related to the uniformity in the distributions of ‖B[i]‖’s and ‖C[i]‖’s. We tentatively discuss the implication of the assumption in Appendix F. We remark the assumption
indicates it is unrealistic that the new probabilities can outperform the ones in LADIES, as distributions of datasets can vary. Nevertheless, as shown in Section 6 it can be an effective attempt to improve the prediction accuracy of LADIES by simply adopting the conservative sampling scheme.

heading: 5 DEBIAS SAMPLING WITHOUT REPLACEMENT
text: The theoretical results derived by Chen et al. (2018b); Zou et al. (2019) guarantee their methods are unbiased approximations of the embedding aggregations in GCN training, with an independent and identically distributed (i.i.d.) sub-sample assumption. However, biases are induced in the implementations due to the usage of sampling without replacement. As illustrated by the U-shape curves of ogbn-arxiv and ogbn-mag datasets in Figure 1, in the long run the matrix approximation error of LADIES even increases with the number of sub-samples. The trend indicates the existence of biases since given an unbiased estimation, the matrix approximation error is expected to decrease with more sub-samples used. In particular, such bias is more significant in sparse graphs with small average degrees. In the following subsections, we analyze the implementation of sampling without replacement in LADIES to demonstrate the origin of the bias, and propose a new correction method to debias the estimation while keep using sampling without replacement.
Remark. We insist on the usage of sampling without replacement because it does diminish the variance of the estimation. This practice indeed works for node-wise sampling methods adopting simple random sampling (SRS, sampling with all-equal probabilities). Specifically, assuming m items are uniformly sampled without replacement from a population of size n, using the same formula for estimation as in sampling with replacement is still unbiased and enjoys the reduction of variance by a finite population correction (FPC) factor n−mn−1 (Lohr, 2019, Section 2.3).

heading: 5.1 WEIGHTED RANDOM SAMPLING WITHOUT REPLACEMENT
text: We first introduce the implementation and necessary notations of non-uniform sampling without replacement applied in FastGCN and LADIES. Given a set V = [n] representing the indices of n items (matrices in the layer-wise sampling setting) {Xi}ni=1 and associated sampling probabilities {pi}ni=1, the algorithm can be stated as a sequential procedure named as weighted random sampling (WRS) (Efraimidis & Spirakis, 2006, Algorithm D). Specifically, aside from the set Sk of k previously sampled indices (0 ≤ k ≤ m− 1, S0 := ∅), the k + 1-th random index Ik+1 is sampled from the set V − Sk of the rest n− k indices with probabilities
p (0) i := pi, ∀i ∈ V = [n],
p (k) i := pi∑ j∈V−Sk pj , ∀k ∈ [m− 1], i ∈ V − Sk.
In the implementation of previous layer-wise sampling methods, the formula used to estimate the target ∑n i=1Xi is displayed (after adapted to the notations in this subsection) as follows:
1
m m−1∑ k=0 XIk+1/pIk+1 , (4)
which is derived under the sampling with replacement assumption. With the notations introduced, we are now able to analyze the effect of directly using Equation (4) while the WRS algorithm is performed. The expectation of a certain summandXIk+1/pIk+1 will be
E XIk+1 pIk+1 = E E XIk+1 p (k) Ik+1 p (k) Ik+1 pIk+1 | Fk  = E[ 1∑ i∈V−Sk pi ∑ i∈V−Sk Xi ] , (5)
where Fk is the σ-algebra generated by the random indices inside the corresponding set Sk,∀k = 0, 1, · · · ,m − 1, and the second equation holds because p (k) Ik+1
pIk+1 = 1∑ i∈V−Sk pi
is Fk-measurable. The expectation is in general unequal to the target ∑n i=1Xi for k > 0, except for some extreme conditions such as all-equal pi’s. The bias in each summand (except for the first term with k = 0) accumulates and results in the ultimate bias in the given estimation.

heading: 5.2 A REMEDY TO ATTAIN UNBIASED ESTIMATION
text: The bias induced by the sequential WRS algorithm has been extensively analyzed by many other studies, especially the ones on stochastic gradient estimators (Liang et al., 2018; Liu et al., 2019; Kool et al., 2019). Given a sequence of random indices sampled through WRS, as far as we know there are two genres to assign coefficients to summands in Equation (4). Both of the two methods relate to the stochastic sum-and-sample estimator (Liang et al., 2018; Liu et al., 2019), which can be derived from Equation (5). Using the fact E XIk+1pIk+1 ∑ i∈V−Sk pi = E [∑ i∈V−Sk Xi ] , a stochastic
sum-and-sample estimator of ∑n
i=1Xi can be immediately constructed as
Πk+1 = ∑ j∈Sk Xi + XIk+1 p (k) Ik+1 ,∀k = 0, 1, · · ·m− 1. (6)
(The proof of unbiasedness is brief and provided by Kool et al. (2019, Appendix C.1).) To minimize the variance, Liang et al. (2018); Liu et al. (2019) develop the first genre to focus on the last estimator Πm and propose methods to pick the initial m − 1 random indices. Kool et al. (2019, Theorem 4) turn to the second genre which utilize Rao-Blackwellization (Casella & Robert, 1996) of Πm.
In fast training for GCN, both of the two genres are somewhat inefficient from a practitioner’s perspective. The first genre works well when ∑ i∈Sm−1 pi is close to 1, otherwise the last term in Πm, XIk+1
p (k) Ik+1
, will brings in large variance and reduce the sample efficiency; for the second genre,
the time cost to perform Rao-Blackwellization (Kool et al., 2019) is extremely high (O(2m) even with approximation by numerical integration) and conflicts with the purpose of fast training. To overcome the issues of the two existing genres, we propose an iterative method to fully utilize each estimator Πk+1 with acceptable runtime to decide the coefficients for each term in Equation (4).
Conceptually, we recursively perform the weighted averaging below:
Z0 := 0, Zk+1 := (1− αk+1)Zk + αk+1Πk+1,∀k = 0, 1, · · · ,m− 1,
where αk+1 is a constant decided by k. When α1 = 1, Z1 = Π1 = XI1 pI1
is unbiased and the unbiasedness of Zm can be obtained by induction as each Πk+1 is unbiased as well. For the choice of αk+1’s, we intentionally specify αk+1 = n(n−k)(k+1) , motivated by the preference that if all pi’s are 1/n, the output coefficients of the algorithm will be all 1/m, the same as the ones in an SRS setting.
In practice, it suffices for the debiasing algorithm to find the coefficients for each Xk+1, which can be obtained in advance of the final computation of Zm. We, therefore, adapt the preceding recursive estimation to iterative updates of coefficients, and conclude the procedure in Algorithm 1 in Appendix A. The time complexity of Algorithm 1 is O(m2), as in iteration k + 1 we update the coefficients for the first k + 1 random indices sampled. The time complexity is comparable to the one of embedding aggregation in layer-wise sampling, as shown in Appendix D.

heading: 5.3 EFFECTS OF DEBIASING
text: As an immediate result, our proposed debiasing method can significantly improve the matrix approximation error in a single step of embedding aggregation. In Figure 1, we observe in all the five datasets, the debiased variant of LADIES has a smaller error, and its error ordinarily drops with the number of sub-samples. When combining the debiasing with the new probabilities raised in Section 2.2, the effect of debiasing is significant only when the graphs (ogbn-arxiv and ogbn-mag) are sparse and when the sub-sample sizes are large enough. We have the following conjectures for this phenomenon. The insignificant improvement is caused by the flat distribution of our proposed probabilities, as the debiasing procedure will maintain the estimation if the sampling probabilities are all-equal. For sparse graphs, however, the number of sampled columns (m) is comparable to the neighborhood size of the batch nodes (n), and as the biases accumulate the effect of the debiasing algorithm emerges.
We further comment that debiasing is not a panacea for GCN accuracy, and the results in Figure 1 cannot fully predict the corresponding GCN accuracy. Conceptually, the columns sampled are
mainly decided by the sampling probabilities, and a debiasing procedure solely adjusts the weights for different columns. This adjustment partially overlaps with the updates of the parametersW (l)’s in a GCN, and therefore the improvement of debiasing to the GCN accuracy will be limited. We claim debiasing should be considered as a means to save some parameter updates, more than a technique designed to improve the GCN accuracy. We will empirically verify the effect of debiasing in the next section.

heading: 6 EXPERIMENTS
text: In this section, we empirically evaluate the performance of each method on five node prediction datasets: Reddit, ogbn-arxiv, ogbn-proteins, ogbn-mag, ogbn-products (See Table 2). We denote “LADIES+flat”, “LADIES+debiased”, and “LADIES+flat+debiased” respectively as the variants of LADIES with the improvements from Section 4, from Section 5, and from both. We compare our methods to the original GCN with mini-batch stochastic training (denoted by full-batch), two layer-wise sampling methods: FastGCN and LADIES. Apart from that, we also implement the other common fast GCN training methods, including GraphSAGE (Hamilton et al., 2017) (vanilla nodewise sampling while keep using the GCN architecture), VR-GCN (Chen et al., 2018a), and subgraph sampling method GraphSAINT (Zeng et al., 2019).
For the training settings, we use a 2-layer GCN for each task with an ADAM optimizer of learning rate 0.001. (Due to limited computational resources, we have to use the shallow GCN since the full-batch method and node-wise sampling methods require much more GPU memory even when L = 3.) The number of hidden variables is 256 and the batch size is 512. For layer-wise sampling methods, we adopt both the classical setting that in each layer the numbers of columns sampled are equal to the batch size 512, and an “increasing” setting (denoted with (2)) that twice nodes will be sampled in the next layer. For the node-wise sampling methods, we combine the settings in the original papers of GraphSAGE (Hamilton et al., 2017) VR-GCN (Chen et al., 2018a) that we separately set the number of neighbors per node as 2 (denoted with (2)) and 10 (denoted with (10)). For the subgraph sampling method GraphSAINT (Zeng et al., 2019), the subgraph size is by default equal to the batch size. The experiment results are reported in Table 1, with the means and associated standard deviations based on 5 runs. More details of the settings are deferred to Appendix C.1.

heading: 6.1 TRAJECTORIES OF MODEL CONVERGENCE
text: We first compare the convergence rates of different methods. The results on Reddit, ogbn-proteins, and ogbn-products are shown in Figure 3, and the results on other datasets are deferred to Figure 4 in Appendix C.2 due to limited space.
In Figure 3 (and Figure 4), we observe LADIES (the green curve) tends to move fast at the beginning of training while slow down after the first stage. In contrast, the debiased variants of LADIES have convergence patterns similar to the outstanding while time-consuming methods Full-batch and VRGCN, and seem to accelerate the training of LADIES at the local convergence stage. Figure 3 also demonstrates that in most cases, “LADIES+debiased” converges slower than “LADIES+flat”, which indicates that regarding convergence, the effect of debiasing is not as important as choosing a proper sampling scheme. This observation matches our analysis in Section 5.3 that debiasing can save some turns of parameter updates and sampling probabilities decide the overall quality of the approximation to embedding aggregations.

heading: 6.2 PREDICTION ACCURACY
text: The prediction accuracy on test sets of different datasets is reported in Table 1. We observe the new sampling probabilities and the debiasing algorithm show consistent improvement over existing layer-wise sampling methods, FastGCN and LADIES. In addition, the experiments further verify the previous observation that the two remedies we apply have the overlapping effect, which is implied in Figure 1 and Figure 3 that (when sub-sample size is small) the improvement of “LADIES+flat+debiased” over its ablation method “LADIES+flat” is insignificant. As for the four benchmark methods, the full-batch method and two node-wise sampling methods have outstanding performance on the accuracy, as they benefit from the vast amount of computation. The graphSAINT method demonstrates a tendency to perform well on the dataset with a high average degree, such as ogbn-proteins. We comment that when the amount of computation is controlled, layer-wise sampling can enjoy similar prediction accuracy to node-wise sampling, as shown in Table 1.
We close this section with a remark on the seemingly strange phenomenon that some efficient GCNs aim to mimic the original full-batch GCN while having a higher prediction accuracy on some datasets. We speculate the reason behind this phenomenon is that a good approximation can recover the principal components in the original embedding matrix, and also restrain the noise via the sparse / low-rank structure. A similar observation (Sanyal et al., 2018) is found in Convolutional Neural Networks (CNN) as well, that applying a low-rank regularizer, such as SVD, to the representation of the intermediate layers can improve the prediction errors of the CNN models.

heading: 7 CONCLUSION
text: We carefully analyze layer-wise sampling in this work and make two improvements to current layerwise sampling models. We first show that a conservative choice of sampling probabilities outperforms the existing ones. The latter probabilities assume the L2 norm of embedding is proportional to the norm of the corresponding column of the renormalized Laplacian matrix, which is usually too strong in practice. We further propose a new correction method to debias the estimation in layer-wise sampling through iterative updates of coefficients for columns sampled. The empirical experiments justify our proposed method and show that our method achieves high accuracy close to the SOTA node-wise sampling method, VR-GCN.

heading: A ALGORITHM TO SPECIFY THE COEFFICIENTS FOR DEBIASING
text: The algorithm described in Section 5.2 is provided in this section.
Algorithm 1: Iterative updates of coefficients to construct the ultimate debiased estimator Zm Input: probabilities {pi}ni=1, random indices {Ik+1} m−1 k=0 generated by WRS with {pi}ni=1 Output: a length m coefficient vector β Initialize the vector β = 0 ∈ Rm, sum of probabilities pS = 0; for k ← 0 to m− 1 do
α = n(n−k)(k+1) ; β[k] = α(1− pS)/pIk+1 ; for j ← 0 to k − 1 do
β[j] = (1− α)β[j] + α; end pS = pS + pIk+1 ;
end return β;
B INFORMATION OF DATASETS
Table 2 summaries the information of datasets used in our experiments.

heading: C ADDITIONAL EXPERIMENT DETAILS
text: C.1 ADDITIONAL DETAILS ON EXPERIMENTAL SETUPS
In this section, we describe the additional details of experiment setups for Section 6. All the models are implemented by Pytorch. We use one Tesla V100 SXM2 32GB GPU with 10 CPU threads to train the models in Section 6. Our implementation of Full-batch method, FastGCN, LADIES are adapted from the codes by Zou et al. (2019); our implementation of vanilla node-wise sampling, VR-GCN, GraphSAINT is adapted from the codes by Cong et al. (2020). For the vanilla node-wise sampling method, there are several variants of structures Ying et al. (2018) while we fix the model structure as GCN in our experiments for fair comparison. We use ELU as the activation function in the convolutional layer for all the models: ELU(x) = x for x > 0, ELU(x) = exp(x) − 1 for x ≤ 0. We choose dropout rate as 0.2, which means 20 percents of units are randomly dropped during the training. Validation and testing are performed with Full-batch inference (using all possible neighbors) on validation and testing nodes. Note that some existing Pytorch implementations of GCNs involve several ad-hoc tricks, such as row-normalizing sampled Laplacian matrix. For the accuracy evaluation experiments in Section 6, we stop training when the validation F1 score does not increase for 200 batches. For a fair comparison, we remove certain tricks in our experiments, such as normalization of each rows in the sampled Laplacian matrix in layer-wise sampling. Such a trick may help in the practice but it might not be compatible to some other methods and is out of the focus of our study. We use the metrics in Table 2 to evaluate the accuracy of each method. Concretely, Reddit is a multi-class classification task and we use the Micro-F1 score with function “sklearn.metrics.f1 score”. For OGB data, we use the built-in evaluator function in module ogb by Hu et al. (2020).
C.2 ADDITIONAL RESULTS ON MODEL CONVERGENCE
The convergence results on ogbn-arxiv and ogbn-mag are provided in Figure 4. The setting of each model is the same as in Figure 3.
C.3 ADDITIONAL REGRESSION EXPERIMENTS
For the regression experiments in Section 4.1, we train a 3-layer GCN without sampling, i.e., fullbatch SGD training, with 256 hidden variables per layer. The batch size is 512. Early stopping training policy is applied. We also remark that these experiments are conducted to check the assumption of importance sampling, rather than pursuing SOTA performance. When we finish training the model, the norms of rows in HW are extracted through a full-batch inference with all training nodes. Hence, we only record the pairs (x(l)i , yi)’s for sampled nodes in regression, where i indicates the index of the sampled nodes in the l-th layer.
Figure 5 is a supplementary to Figure 2. The assumption: ‖(HW )[i]‖ ∝ ‖P [i]‖ does not hold on all of these datasets. We do not have the regression result on the ogbn-product dataset since the training of 3-layer GCN fails due to the memory limitation.
C.4 DISTRIBUTIONS OF THE NORM OF ROWS / COLUMNS IN RELATED MATRICES
Figure 6 demonstrates the distribution of ‖P [i]‖ and ‖(HW )[i]‖ (layer 1, 2, 3) for Reddit, ogbnarxiv, ogbn-protein and ogbn-mag datasets. The ‖(HW )[i]‖’s are obtained from the experiment in Section C.3. The outliers larger than the 99.9% quantile or small than the 0.1% quantile are removed.
As shown in the histograms, our analysis regarding Assumption 1 tends to hold generally on these datasets. For the norms of columns in P (as a replacement for QP for clarity), we observe there are some columns with large norms far beyond the average. Those columns contribute a lot to the quadratic mean, which results in a huge C1 in Assumption 1. In contrast, the norms of rows inHW concentrate around their average, inducing a small C2. Those facts together with Assumption 1 and Lemma 1 in the main paper explain why our proposed sampling probabilities are more proper for some real datasets.

heading: D TIME COMPLEXITY ANALYSIS
text: We analyze the complexity of vanilla node-wise sampling and layer-wise sampling in this section. The analysis is adapted from the work by Zou et al. (2019), but we show a lighter bound for layerwise sampling. For l such that 0 ≤ l ≤ L − 1, the propagation formulas for sampling based GCN can be formulated as:
Z̃(l+1) = P̄ (l)H̃(l)W (l),
where H̃(l) = Z̃(l) ∈ Rsl×p, P̂ (l) ∈ Rsl+1×sl , W (l) ∈ Rp×p. In particular, for LADIES, P̄
(l) LADIES = Q (l+1)PS(l).
For simplicity, we suppose that the number of hidden variables in each layer is fixed as p, the same as the dimension of H(0). The batch size and the numbers of nodes sampled in each layer equal are set all equal to a fixed constant s. We assume the number of sampled neighbors per node in node-wise sampling is b. We denote the maximal degree of all the nodes in the graph as c. The computational cost of the propagation comes from two parts: the linear transformation, a dense matrix product, H̃(l)W (l) and the node aggregation, a sparse matrix product, P̄ (l)(H̃(l)W (l)). The time complexity is summarized in Table 3. We additionally comment although the time cost of two parts both linearly depend on the number of nodes involved (number of non-zero elements in Q(l+1)), the node aggregation part usually dominates since the sparse matrix product involved is less efficient than the dense matrix product involved in a modern computer.
The linear transformation, H̃(l)W (l) is dense matrix production. The cost depends on the shape of two matrices, and is given as O(slp2). LADIES fixes sl as s for each layer, so O(slp2) = O(sp2). For node-wise sampling sl = sbL−l, since number of node exponentially grows. Thus by summation over all the layers, we have the results in the second row in Table 3.
The node aggregation, P̄ (l)(H̃(l)W (l)) is a sparse matrix production, since P̄ (l) is sparse. For simplicity, we denote H̃(l)W (l) as C(l) ∈ Rsl×p. Thus the time complexity of this sparse matrix production becomes O(nnz(Pl)p), where nnz(Pl) is the number of non-zero entries in P̄ (l). For layer-wise sampling, since we sample s nodes for each layer and each node has at most c neighbors, so nnz(Pl) ≤ sc. For node-wise sampling, since each node has b neighbors and the neighbors are not shared by all the nodes in each layer, nnz(Pl) = bsl = sbL+1−l. By summation over all the layers, we attain the results in the first row of Table 3.

heading: E RESULTS IN APPROXIMATE MATRIX MULTIPLICATION
text: In this section, we revisit approximate matrix multiplication to derive the previous layer-wise sampling methods. Specifically, the sampling matrix S used in FastGCN and LADIES can be decomposed as S = ΠΠT , where Π ∈ Rn×d is a sub-sampling sketching matrix defined as follows: Definition E.1 (Sub-sampling sketching matrix). Consider a discrete distribution which draws i with probability pi > 0,∀i ∈ [n]. For a random matrix Π ∈ Rn×d, if Π has i.i.d. columns and each column Π(j) can randomly be 1√
dpi ei with probability pi, where ei is the i-th column of the n-by-n
identity matrix In, then Π is called a sub-sampling sketching matrix with sub-sampling probabilities {pi}ni=1.
With this definition, we introduce a result in AMM to construct the sub-sampling sketching matrix, which coincides with the conclusion in FastGCN and LADIES.
Theorem E.1 (Theorem 1 (Drineas et al., 2006)). Suppose B ∈ RnB×n, C ∈ Rn×nC , the number of sub-sampled columns d ∈ Z+ such that 1 ≤ d ≤ n, and the sub-sampling probabilities {pi}ni=1 are such that ∑n i=1 pi = 1 and such that for a quality coefficient β ∈ (0, 1]
pi ≥ β ‖B[i]‖‖C[i]‖∑n
i′=1 ‖B[i ′]‖‖C[i′]‖
,∀i ∈ [n]. (7)
Construct a sub-sampling sketching matrix Π ∈ Rn×d with sub-sampling probabilities {pi}ni=1 as in Definition E.1, and let BΠΠTC be an approximation to BC. Let δ ∈ (0, 1) and η =
1 + √ (8/β) log(1/δ). Then with probability at least 1− δ,
‖BC −BΠΠTC‖2F ≤ η2
βd ‖B‖2F ‖C‖2F . (8)
Remark. The theorem is closely related to Lemma 1 in Appendix B of LADIES, which studies the variance E ‖BC −BSC‖2F . For the choice of sub-sampling probabilities, Equation (7) reproduces the conclusion in FastGCN and LADIES, when we respectively takeB as P andQP .

heading: F PROOF OF LEMMA 1 AND SOME REMARKS
text: To prove Lemma 1 in the main paper, we first adapt a technical lemma (Zou et al., 2019, Lemma 1), which relates the sampling matrix to the variance (expectation of squared Frobenius norm) of the approximate matrix multiplication.
Lemma F.1 (Adapted from Lemma 1 (Zou et al., 2019)). Given two matrices B ∈ RnB×n and C ∈ Rn×nC , for any i ∈ [n] define the positive probabilities pi’s such that ∑n i=1 pi = 1. We further require the probability pi = 0 if and only if the corresponding column B[i] or row C[i] is all-zero. The sub-sampling sketching matrix Π ∈ Rn×d is generated accordingly. Let S := ΠΠT , it holds that
ES [ ‖BSC −BC‖2F ] = 1
d  ∑ i:pi>0 1 pi ∥∥∥B[i]∥∥∥2 · ∥∥C[i]∥∥2 − ‖BC‖2F 
where d is the number of samples.
With the lemma above, the proof of Lemma 1 in the main paper is provided as follows.
Proof. Recall the notation in Lemma 1 in the main paper is simplified as B := QP ,C := HW . As the union of neighbors of nodes inQ cannot cover all the nodes, some columns inB are all-zero, and we accordingly define anQ-measurable matrix L as in Lemma F.1. We have
E [ ‖BS1C −BC‖2F ] =EQ [ ES1 ( ‖BS1C −BC‖2F ∣∣Q)] = 1
d EQ  ∑ i:pi>0 1 pi ∥∥∥B[i]∥∥∥2 · ∥∥C[i]∥∥2 − ‖BC‖2F  .
where the second equation holds as we apply Lemma F.1 to the inner expectation in the right-hand side of the first line. Plugging pi ∝ ‖B[i]‖ (Equation (3) in the main paper) into the preceding probabilities pi’s, we reach
E [ ‖BS1C −BC‖2F ] = EQ [(∑ i:pi>0 ∥∥B[i]∥∥)(∑i:pi>0 ∥∥B[i]∥∥ ∥∥C[i]∥∥2)] d − EQ [ ‖BC‖2F ] d
= 1
d EQ
[( n∑
i=1 ∥∥∥B[i]∥∥∥)( n∑ i=1 ∥∥∥B[i]∥∥∥∥∥C[i]∥∥2 )] − 1 d EQ [ ‖BC‖2F ] .
As computed by Zou et al. (2019), the variance of LADIES is similarly given as
E [ ‖BS0C −BC‖2F ] = 1
d EQ  ∑ i:pi>0 ∥∥∥B[i]∥∥∥2  ∑ i:pi>0 ∥∥C[i]∥∥2 − 1 d EQ [ ‖BC‖2F ] .
Consequently, to prove the lemma it suffices to show that ∑ i:pi>0 ∥∥∥B[i]∥∥∥  ∑ i:pi>0 ∥∥∥B[i]∥∥∥∥∥C[i]∥∥2  ≤  ∑ i:pi>0 ∥∥∥B[i]∥∥∥2  ∑ i:pi>0 ∥∥C[i]∥∥2  , (9)
and the inequality above follows with Assumption 1. Specifically, plugging the inequality ‖C[i]‖2 ≤ C2‖C‖2F /n,∀i ∈ [n] in the left-hand-side above, we have ∑
i:pi>0
∥∥∥B[i]∥∥∥  ∑
i:pi>0
∥∥∥B[i]∥∥∥∥∥C[i]∥∥2  ≤ ( n∑
i=1
‖B[i]‖ )2 C2 n ‖C‖2F = m C1 ‖B‖2F C2 nm m‖C‖2F ,
in which the last equation comes from the definitionC1 := ‖B‖2F /m ( ∑n i=1 ‖B[i]‖/m) 2 . To close the proof, we utilize the inequality 1C2 ‖C‖ 2 F /n ≤ ‖C[i]‖2 and bound m‖C‖2F by nC2 ∑ i:pi>0 ∥∥C[i]∥∥2. Finally we attain Equation (9) with the core assumption C1
C22 ≥ 1. ♦
Remark. In Assumption 1 we indeed implicitly assume ‖B[i]‖’s follow a long-tail distribution that most norms are around the average while a few columns have large norms. The high nonuniformity makes the average of squared norms much larger than the square of averaged norms. For ‖C[i]‖’s, considering the normalization techniques (such as batch or layer normalization) to stabilize the scale of the parameters, they tend to not vary widely, which implies a small C2. The numerical experiments on the comparison of approximation error (see Figure 1) and the histograms of the norms in trained models shown in Figure 6 further validate the assumption. Based on the empirical analysis above, we claim the assumption is mild and tends to hold at least for some datasets.

heading: G SUPPLEMENTARY EXPERIMENTS OF REBUTTAL REVISION
text: G.1 SUPPLEMENTARY REGRESSION EXPERIMENTS
In addition to the regression experiments in Section 4.1 and Appendix C.3, we also present the regression results in Figure 7 and Figure 5, for a 3-layer GCN with our LADIES+flat+debias sampler and LADIES sampler respectively Both of them have 512 nodes sampled at each layer and keep other hyper-parameters the same as regression experiments. The regression results in the following figures show similar patterns to Figures 2 and 5.
G.2 SAMPLING TIME AND TRAINING TIME
We compare the sampling time per batch for 1-layer GCN with layer-wise sampling methods (FastGCN, LADIES, and our proposed methods) and GraphSAGE. The time is presented in milliseconds. The batch size is 512, and the number of sampled nodes is 512 or 1024. The average sampling time (followed by standard deviation) over 200 batches is presented in Table 4. We note that the sampling time may involve some overhead costs. For example, the input Laplacian matrix is Scipy-sparematrix on the CPU, while in sampling, it is converted to a PyTorch-sparse-matrix.
By Table 4, we conclude that the cost of debiasing algorithm is acceptable. Moreover, since the debiasing only depends on the number of nodes sampled, its time cost will be dwarfed by sampling on very large graphs. For example, sampling 512 nodes, the average batch sampling time for “LADIES”, “LADIES + debiased”, “LADIES + flat + debiased”, are 8.3 ± 0.1, 11.7 ± 0.2, 11.7± 0.5 respectively on ogbn-arxiv while 83.4± 1.3 and 83.7± 0.8 and 83.5± 0.6 respectively on ogbn-products data. We remark that the node-wise sampling takes a significantly longer time because individually sampling from each row in the re-normalized Laplacian matrix (stored as a sparse matrix in implementation) leads to a large overhead cost.
We present the training time per batch for 2-layer and 3-layer GCNs with different methods in Tables 5 and Table 6 respectively. The time is presented in millisecond and averaged over 110 batches, where we discard the first 20 and the last 20 out of 150 total batches to disregard potential warm-up time for GPU. The other settings are kept the same as our experiments of accuracy evaluation in Table 1. We note that the timing on GPU is sensitive to the hardware and has a relatively large standard deviation.
As presented in Tables 5 and 6, our proposed methods have similar training time with LADIES due to the same propagation scheme of GCN with layer-wise sampling strategy. The VR-GCN generally shows superiority in prediction accuracy (see Table 1). However, it also takes a significantly longer time in training since its propagation involves using and updating historical activation.
G.3 FINER VISUALIZATION OF FIGURE 3 AND FIGURE 4
To better differentiate the distinct sampling-based methods in Figure 3 and Figure 4, we provide Figure 9 to group the methods of the same type by the same color and demonstrate the trajectories of model convergence. We also provide Figure 10 to only compare the layer-wise sampling methods (FastGCN, LADIES, LADIES+flat, LADIES+debiasd, LADIES+flat+debiasd).

Review summaries:
This paper studies existing sampling schemes employed in training graph neural network architectures (e.g., FastGCN [Chen et al.] and LADIES [Zon et al.]) that are improvements to graph convolutional networks (GCN) [Kipf & Welling 2017] from a matrix approximation perspective (looking a Frobenius norm).  The paper focuses on layer-wise sampling (cf. node and subgraph sampling) and observes that there are two drawbacks in common practice layer-wise sampling. First drawback is current probabilities (probability distributions) used for sampling are sub-optimal. The reason the paper gives is that a core assumption made by certain GCN schemes such as FastGCN and LADIES do not hold in datasets such as Reddit and OGB. Second drawback is the implementations of these schemes slightly deviate from their respective theoretical results. The implementations use sampling without replacement, and thus introduces bias.  To this, the paper presents new sampling probabilities (uniform proposal distribution), and an algorithm to reduce the bias. With these adjustments, training of GCN converges faster and potentially leads to higher prediction accuracy (for node prediction). Theoretical analysis is presented and experiments are conducted on common benchmarks (Reddit, ogbn-{arxiv,proteins,products,mag}). 

This paper revisits layer-wise sampling methods of graph neural networks by addressing issues such as sub-optimal sampling probabilities and the approximation bias due to the usage of sampling without replacement. This paper also presents a new metric to evaluate the performance of the sampling strategy. Through comprehensive experiments, the authors validate the effectiveness of their proposed methods.