Title: STOCHASTIC TRAINING OF GRAPH CONVOLUTIONAL NETWORKS

Abstract: Graph convolutional networks (GCNs) are powerful deep neural networks for graph-structured data. However, GCN computes nodes' representation recursively from their neighbors, making the receptive field size grow exponentially with the number of layers. Previous attempts on reducing the receptive field size by subsampling neighbors do not have any convergence guarantee, and their receptive field size per node is still in the order of hundreds. In this paper, we develop a preprocessing strategy and two control variate based algorithms to further reduce the receptive field size. Our algorithms are guaranteed to converge to GCN's local optimum regardless of the neighbor sampling size. Empirical results show that our algorithms have a similar convergence speed per epoch with the exact algorithm even using only two neighbors per node. The time consumption of our algorithm on the Reddit dataset is only one fifth of previous neighbor sampling algorithms.

1. INTRODUCTION:

Graph convolution networks (GCNs) (Kipf & Welling, 2017) generalize convolutional neural networks (CNNs) (LeCun et al., 1995) to graph structured data. The "graph convolution" operation applies same linear transformation to all the neighbors of a node, followed by mean pooling. By stacking multiple graph convolution layers, GCNs can learn nodes' representation by utilizing information from distant neighbors. GCNs have been applied to semi-supervised node classification (Kipf & Welling, 2017) , inductive node embedding (Hamilton et al., 2017a) , link prediction (Kipf & Welling, 2016; Berg et al., 2017) and knowledge graphs (Schlichtkrull et al., 2017) , outperforming multi-layer perceptron (MLP) models that do not use the graph structure and graph embedding approaches (Perozzi et al., 2014; Tang et al., 2015; Grover & Leskovec, 2016) that do not use node features.

However, the graph convolution operation makes it difficult to train GCN efficiently. A node's representation at layer L is computed recursively by all its neighbors' representations at layer L -1. Therefore, the receptive field of a single node grows exponentially with respect to the number of layers, as illustrated in Fig. 1(a) . Due to the large receptive field size, Kipf & Welling (2017) proposed training GCN by a batch algorithm, which computes the representation for all the nodes altogether. However, batch algorithms cannot handle large scale datasets because of their slow convergence and the requirement to fit the entire dataset in GPU memory. Hamilton et al. (2017a) made an initial attempt on developing stochastic algorithms to train GCNs, which is referred as neighbor sampling (NS) in this paper. Instead of considering all the neighbors, they randomly subsample D (l) neighbors at the l-th layer. Therefore, they reduce the receptive field size to l D (l) , as shown in Fig. 1(b) . They found that for two layer GCNs, keeping D (1) = 10 and D (2) = 25 neighbors can achieve comparable performance with the original model. However, there is no theoretical guarantee on the predictive performance of the model learnt by NS comparing with the original algorithm. Moreover, the time complexity of NS is still D (1) D (2) = 250 times larger than training an MLP, which is unsatisfactory.

In this paper, we develop novel stochastic training algorithms for GCNs such that D (l) can be as low as two, so that the time complexity of training GCN is comparable with training MLPs. Our methods are built on two techniques. First, we propose a strategy which preprocesses the first graph convolution layer, so that we only need to consider all neighbors within L-1 hops instead of L hops. This is significant because most GCNs only have L = 2 layers (Kipf & Welling, 2017; Hamilton et al., 2017a) . Second, we develop two control variate (CV) based stochastic training algorithms. We show that our CV-based algorithms have lower variance than NS, and for GCNs without dropout, our algorithm provably converges to a local optimum of the model regardless of D (l) .

We empirically test on six graph datasets, and show that our techniques significantly reduce the bias and variance of the gradient from NS with the same receptive field size. Our algorithm with D (l) = 2 achieves the same predictive performance with the exact algorithm in comparable number of epochs on all the datasets, while the training time is 5 times shorter on our largest dataset.

2. BACKGROUNDS:

We now briefly review graph convolutional networks (GCNs) (Kipf & Welling, 2017) and the neighbor sampling (NS) algorithm (Hamilton et al., 2017a) .

2.1. GRAPH CONVOLUTIONAL NETWORKS:

The original GCN was presented in a semi-supervised node classification task (Kipf & Welling, 2017) . We follow this setting throughout this paper. Generalization of GCN to other tasks can be found in Kipf & Welling (2016) ; Berg et al. (2017) ; Schlichtkrull et al. (2017) and Hamilton et al. (2017b) . In the node classification task, we have an undirected graph G = (V, E) with V = |V| vertices and E = |E| edges, where each vertex v consists of a feature vector x v and a label y v . The label is only observed for some vertices V L and we want to predict the label for the rest vertices V U := V\V L . The edges are represented as a symmetric V × V adjacency matrix A, where A v,v is the weight of the edge between v and v , and the propagation matrix P is a normalized version of A: Ã = A + I, Dvv = v Ãvv , and P = D-1 2 Ã D-1 2 . A graph convolution layer is defined as

EQUATION

where H (l) is the activation matrix in the l-th layer, whose each row is the activation of a graph node. H (0) = X is the input feature matrix, W (l) is a trainable weight matrix, σ(•) is an activation function, and Dropout p (•) is the dropout operation (Srivastava et al., 2014) with keep probability p.

Finally, the loss is defined as

L = 1 |V L | v∈V L f (y v , Z (L) v )

, where f (•, •) can be the square loss, cross entropy loss, etc., depending on the type of the label.

When P = I, GCN reduces to a multi-layer perceptron (MLP) model which does not use the graph structure. Comparing with MLP, GCN is able to utilize neighbor information for node classification. We define n(v, L) as the set of all the L-neighbors of node v, i.e., the nodes that are reachable from v within L hops. It is easy to see from Fig. 1 (a) that in an L-layer GCN, a node uses the information from all its L-neighbors. This makes GCN more powerful than MLP, but also complicates the stochastic training, which utilizes an approximated gradient ∇L ≈ 1

|V B | v∈V B ∇f (y v , Z (L) v ),

where V B ⊂ V L is a minibatch of training data. The large receptive field size | ∪ v∈V B n(v, L)| per minibatch leads to high time complexity, space complexity and amount of IO. See Table 1 for the average number of 1-and 2-neighbors of our datasets.

2.2. ALTERNATIVE NOTATION:

We introduce alternative notations to help compare different algorithms. Let U (l) = P H(l) , or u l) neighbors. In this way, they reduce the receptive field size from l) ). Neighbor sampling can also be written in a matrix form as

(l) v = v ∈n(v,1) P v,v h(l) v ,

v∈D (l) p v h v , where D (l) ⊂ [D] + is a subset of D (

| ∪ v∈V B n(v, L)| to O(|V B | L l=1 D (

EQUATION

where P (l) is a sparser unbiased estimator of P , i.e., E P (l) = P . The approximate prediction Z

(L) N S

used for testing and for computing stochastic gradient1 

EQUATION

CV,v ) during training. The NS estimator u N S is unbiased. However it has a large variance, which leads to biased prediction and gradients after the non-linearity in subsequent layers. Due to the biased gradients, training with NS does not converge to the local optimum of GCN. When D (l) is moderate, NS may has some regularization effect like dropout (Srivastava et al., 2014) , where it drops neighbors instead of features. However, for the extreme ease D (l) = 2, the neighbor dropout rate is too high to reach high predictive performance, as we will see in Sec. 5.4. Intuitively, making prediction solely depends on one neighbor is inferior to using all the neighbors. To keep comparable prediction performance with the original GCN, Hamilton et al. (2017a) use relatively large D (1) = 10 and D (2) = 25. Their receptive field size D (1) × D (2) = 250 is still much larger than MLP, which is 1.

3. PREPROCESSING FIRST LAYER:

We first present a technique to preprocess the first graph convolution layer, by approximating ADropout p (X) with Dropout p (AX). The model becomes

Z (l+1) = Dropout p (P H (l) )W (l) , H (l+1) = σ(Z l+1 ).

(3)

This approximation does not change the expectation because E ADropout p (X) = E Dropout p (AX) , and it does not affect the predictive performance, as we shall see in Sec. 5.1.

The advantage of this modification is that we can preprocess U (0) = P H (0) = P X and takes U (0) as the new input. In this way, the actual number of graph convolution layers is reduced by one -the first layer is merely a fully connected layer instead of a graph convolution one. Since most GCNs only have two graph convolution layers (Kipf & Welling, 2017; Hamilton et al., 2017a) , this gives a significant reduction of the receptive field size from the number of L-neighbors

| ∪ v∈V B n(v, L)| to the number of L -1-neighbors | ∪ v∈V B n(v, L -1)|.

The numbers are reported in Table 1 .

4. CONTROL VARIATE BASED STOCHASTIC APPROXIMATION:

We now present two novel control variate based estimators that have smaller variance as well as stronger theoretical guarantees than NS.

4.1. CONTROL VARIATE BASED ESTIMATOR:

We assume that the model does not have dropout for now and will address dropout in Sec. 4.2. The idea is that we can approximate u = D v=1 p v h v better if we know the latest historical activations hv of the neighbors, where we expect hv and h v are similar if the model weights do not change too fast during the training. With the historical activations, we approximate

EQUATION

)

where v is a random neighbor, and ∆h v = h v -hv . For the ease of presentation, we assume that we only use the latest activation of one neighbor, while the implementation also include the node itself besides the random neighbor, so D (l) = 2. Using historical activations is cheap because they need not to be computed recursively using their neighbors' activations, as shown in Fig. 1(c ).

Unlike NS, we apply Monte-Carlo approximation on v p v ∆h v instead of v p v h v . Since we expect h v and hv to be close, ∆h v will be small and u CV should have a smaller variance than u N S . Particularly, if the model weight is kept fixed, hv should be eventually equal with h v , so that

u CV = 0 + D v=1 p v hv = D v=1 p v h v = u, i.e.

, the estimator has zero variance. The term (Ripley, 2009, Chapter 5) , which has zero mean and large correlation with u N S , to reduce its variance. We refer this stochastic approximation algorithm as CV, and we will formally analyze the variance and prove the convergence of the training algorithm using CV for stochastic gradient in subsequent sections.

CV = u CV -u N S = -Dp v hv + D v=1 p v hv is a control variate

In matrix form, CV computes the approximate predictions as follows, where we explicitly write down the iteration number i and add the subscript CV to the approximate activations2 Z (l+1)

EQUATION

EQUATION

where h(l) CV,i,v stores the latest activation of node v on layer l computed before time i. Formally, let m (l) i ∈ R V ×V be a diagonal matrix, and (m

(l) i ) vv = 1 if ( P (l) i ) v v > 0 for any v .

After finishing one iteration we update history H with the activations computed in that iteration as Eq. ( 6).

4.2. CONTROL VARIATE FOR DROPOUT:

With dropout, the activations H are no longer deterministic. They become random variables whose randomness come from different dropout configurations. Therefore, ∆h v = h v -hv is not necessarily small even if h v and hv have the same distribution. We develop another stochastic approximation algorithm, control variate for dropout (CVD), that works well with dropout.

Our method is based on the weight scaling procedure (Srivastava et al., 2014) to approximately compute the mean

µ v := E [h v ].

That is, along with the dropout model, we can run a copy of the model with no dropout to obtain the mean µ v , as illustrated in Fig. 1(d) . With the mean, we can obtain a better stochastic approximation by separating the mean and variance

u = D v=1 p v [(h v -µ v ) + (µ v -μv ) + μv ] ≈ √ Dp v (h v -µ v )+Dp v ∆µ v + D v=1 p v μv := u CV D ,

where μv is the historical mean activation, obtained by storing µ v instead of h v , and ∆µ = µ v -μv . u CV D an unbiased estimator of u because the term √ Dp v (h v -µ v ) has zero mean, and the Monte-Carlo approximation

D v=1 p v (µ v -μv ) ≈ Dp v ∆µ v does not change the mean. The approximation D v=1 p v (h v -µ v ) ≈ √ Dp v (h v -µ v

) is made by assuming h v 's to be independent Gaussians, which we will soon clarify. The pseudocodes for CV and CVD are in Appendix E.

4.3. VARIANCE ANALYSIS:

NS, CV and CVD are all unbiased estimators of u = v p v h v . We analyze their variance in a simple independent Gaussian case, where we assume that activations are Gaussian random variables Alg.

Estimator Var. from MC. approx.

Var. from dropout:

Exact u = v pvhv 0 s 2 NS uNS = Dp v h v 1 2 v,v (pvµv -p v µ v ) 2 Ds 2 CV uCV = Dp v ∆h v + v pv hv 1 2 v,v (pv∆µv -p v ∆µ v ) 2 Ds 2 + (D -1)s 2 CVD uCV D = √ Dp v (h v -µ v ) 1 2 v,v (pv∆µv -p v ∆µ v ) 2 s 2 +Dp v ∆µ v + v pv μv

Table 2 : Variance of different algorithms in the independent Gaussian case. Wang & Manning (2013) . Without loss of generality, we assume that all the activations h v are one dimensional. We also assume that all the activations h 1 , . . . , h D and historical activations h1 , . . . , hD are independent, where the historical activations hv ∼ N (μ v , s2 v ). We introduce a few more notations. ∆µ v and ∆s 2 v are the mean and variance of

h v ∼ N (µ v , s 2 v ) following

∆h v = h v -hv , where ∆µ v = µ v -μv and ∆s 2 v = s 2 v + s2 v .

µ and s 2 are the mean and variance of v p v h v , where

µ = v p v µ v and s 2 = v p 2 v s 2 v .

Similarly, ∆µ, ∆s 2 , μ and s2 are the mean and variance of v p v ∆h v and v p v hv , respectively. With these assumptions and notations, we list the estimators and variances in Table 2 , where the derivations can be found in Appendix C. We decompose the variance as two terms: variance from Monte-Carlo approximation (VMCA) and variance from dropout (VD).

If the model has no dropout, the activations have zero variance, i.e., s v = sv = 0, and the only source of variance is VMCA. We want VMCA to be small. As in Table 2 , the VMCA for the exact estimator is 0. For the NS estimator, VMCA is

1 2 v,v (p v µ v -p v µ v ) 2 , whose magnitude depends on the pairwise difference (p v µ v -p v µ v ) 2 , and VMCA is zero if and only if p v µ v = p v µ v for all v, v . Similarly, VMCA for both CV and CVD estimators is 1 2 v,v (p v ∆µ v -p v ∆µ v ) 2 , which should be smaller than NS estimator's VMCA if (p v ∆µ v -p v ∆µ v ) 2 < (p v µ v -p v µ v ) 2 ,

which is likely because ∆µ v should be smaller than µ v . Since CV and CVD estimators have the same VMCA, we adopt the CV estimator for models without dropout, due to its simplicity.

The VD of the exact estimator is s 2 , which is overestimated by both NS and CV. NS overestimates VD by D times, and CV has even larger VD. Meanwhile, the VD of the CVD estimator is the same as the exact estimator, indicating CVD to be the best estimator for models with dropout.

4.4. EXACT TESTING:

Besides smaller variance, CV also has stronger theoretical guarantees than NS. We can show that during testing, CV's prediction becomes exact after a few testing epochs. For models without dropout, we can further show that training using the stochastic gradients obtained by CV converges to GCN's local optimum. We present these results in this section and Sec. 4.5. Note that the analysis does not need the independent Gaussian assumption.

Given a model W , we compare the exact predictions (Eq. 1) and CV's approximate predictions (Eq. 5,6) during testing, which uses the deterministic weight scaling procedure. To make predictions, we run forward propagation by epochs. In each epoch, we randomly partition the vertex set V as I minibatches V 1 , . . . , V I and in the i-th iteration, we run a forward pass to compute the prediction for nodes in V i . Note that in each epoch we scan all the nodes instead of just testing nodes, to ensure that the activation of each node is computed at least once per epoch. The following theorem reveals the connection of the exact predictions and gradients, and their approximate versions by CV.

Theorem 1. For a fixed W and any i > LI we have: (1) (Exact Prediction) The activations computed by CV are exact, i.e., Z (l)

CV,i = Z (l) for each l ∈ [L] and H (l) CV,i = H (l) for each l ∈ [L-1]. (2) (Unbiased Gradient) The stochastic gradient g CV,i (W ) := 1 |V B | v∈V B ∇ W f (y v , z (L) CV,i,v ) is an unbiased estimator of GCN's gradient, i.e., E P ,V B g CV,i (W ) = ∇ W 1 |V | v∈V f (y v , z (L) v ).

Theorem 1 shows that at testing time, we can run forward propagation with CV for L epoches and get the exact prediction. This outperforms NS, which cannot recover the exact prediction. Comparing with directly making exact predictions by a batch algorithm, CV is more scalable because it does not need to load the entire graph into memory. The proof can be found in Appendix A.

4.5. CONVERGENCE GUARANTEE:

The following theorem shows that for a model without dropout, training using CV's approximated gradients converges to a local optimum of GCN, regardless of the neighbor sampling size D (l) . Therefore, we can choose arbitrarily small D (l) without worrying about the convergence. Theorem 2. Assume that (1) all the activations are ρ-Lipschitz, (2) the gradient of the cost function ∇ z f (y, z) is ρ-Lipschitz and bounded, (3) g CV (W ) ∞ and g(W ) ∞ = ∇L(W ) ∞ are bounded by G > 0 for all P , V B and W . (4) The loss

L(W ) is ρ-smooth, i.e., |L(W 2 ) - L(W 1 ) -∇L(W 1 ), W 2 -W 1 | ≤ ρ 2 W 2 -W 1 2 ∀W 1 , W 2 , where A, B = tr(A B)

is the inner product of matrix A and matrix B. We randomly run SGD for R ≤ N iterations, where

P R (R = i) = 2γi-ργ 2 i N i=1 (2γi-ργ 2 i ) . Then, for the updates W i+1 = W i -γ i g CV (W i ) and step sizes γ i = min{ 1 ρ , 1 √ N }

, there exists constants K 1 and K 2 which are irrelevant with N , s.t. ∀N > LI,

E R∼P R E P ,V B ∇L(W R ) 2 ≤ 2ρ(L(W 1 ) -L(W * )) + K 2 N + 2(L(W 1 ) -L(W * )) + K 1 √ N .

The proof can be found in Appendix B. Particularly,

lim N →∞ E R∼P R E P ,V B ∇L(W R ) 2 = 0.

Therefore, our algorithm converges to a local optimum.

4.6. TIME COMPLEXITY AND IMPLEMENTATION DETAILS:

Finally we discuss the time complexity of different algorithms. We decompose the time complexity as sparse time complexity for sparse-dense matrix multiplication such as P H(l) , and dense time complexity for dense-dense matrix multiplication such as U (l) W (l) . Assume that the node feature is K-dimensional and the first hidden layer is A-dimensional, Our implementation is similar as Kipf & Welling (2017) . We store the node features in the main memory, without assuming that they fit in GPU memory as Hamilton et al. (2017a) , which makes our implementation about 2 times slower than theirs. We keep the histories in GPU memory for efficiency since they are only LH < K dimensional.

5. EXPERIMENTS:

We examine the variance and convergence of our algorithms empirically on six datasets, including Citeseer, Cora, PubMed and NELL from Kipf & Welling (2017) and Reddit, PPI from Hamilton et al. (2017a), as summarized in Table 1 . To measure the predictive performance, we report Micro-F1 for the multi-label PPI dataset, and accuracy for all the other multi-class datasets. We use the same model architectures with previous papers but slightly different hyperparameters (see Appendix D for the details). We repeat the convergence experiments 10 times on Citeseer, Cora, PubMed and NELL, and 5 times on Reddit and PPI. The experiments are done on a Titan X (Maxwell) GPU.

5.1. IMPACT OF PREPROCESSING:

We first examine the approximation in Sec. 3 that switches the order of dropout and aggregating the neighbors. Let M0 be the original model (Eq. 1) and M1 be our approximated model (Eq. 3), we compare three settings: (1) M0, D (l) = ∞ is the exact algorithm without any neighbor sampling. (2) M1+PP, D (l) = ∞ changes the model from M0 to M1. Preprocessing does not affect the training for D (l) = ∞. (3) M1+PP, D (l) = 20 uses NS with a relatively large number of neighbors. In Table 3 we can see that all the three settings performs similarly, i.e., our approximation does not affect the predictive performance. Therefore, we use M1+PP, D (l) = 20 as the exact baseline in following convergence experiments because it is the fastest among these three settings. 70.8 ± .1 81.7 ± .5 79.0 ± .4 -97.9 ± .04 96.2 ± .04 M1+PP, D (l) = ∞ 70.9 ± .2 82.0 ± .8 78.7 ± .3 64.9 ± 1.7 97.8 ± .05 96.3 ± .07 M1+PP, D (l) = 20 70.9 ± .2 81.9 ± .7 78.9 ± .5 64.2 ± 4.6 97.6 ± .09 96.3 ± .04

Table 3 : Testing accuracy of different algorithms and models after fixed number of epochs. Our implementation does not support M0, D (l) = ∞ on NELL so the result is not reported 

5.2. CONVERGENCE WITH NO DROPOUT:

We now study how fast our algorithms converge with a very small neighbor sampling size D (l) = 2. We compare the following algorithms: (1) Exact, which is M1+PP, D (l) = 20 in Sec. 5.1 as a surrogate of the exact algorithm. (2) NS, which is the NS algorithm with no preprocessing and D (l) = 2. (3) NS+PP, which is same with NS but uses preprocessing. (4) CV+PP, which replaces the NS estimator in NS+PP with the CV estimator. ( 5) CVD+PP, which uses the CVD estimator.

We first validate Theorem 2, which states that CV+PP converges to a local optimum of Exact, for models without dropout, regardless of D (l) . We disable dropout and plot the training loss with respect to number of epochs as Fig. 2 . We can see that CV+PP can always reach the same training loss with Exact, which matches the conclusion of Theorem 2. Meanwhile, NS and NS+PP have a higher training loss because their gradients are biased.

5.3. CONVERGENCE WITH DROPOUT:

Next, we compare the predictive accuracy obtained by the model trained by different algorithms, with dropout turned on. We use different algorithms for training and the same Exact algorithm for testing, and report the validation accuracy at each training epoch. The result is shown in Fig. 3 . We find that CVD+PP is the only algorithm that is able to reach comparable validation accuracy with Exact on all datasets. Furthermore, its convergence speed with respect to the number of epochs is comparable with Exact despite its D (l) is 10 times smaller. Note that CVD+PP performs much better than Exact on the PubMed dataset; we suspect it finds a better local optimum. Meanwhile, simper algorithms CV+PP and NS+PP work acceptably on most of the datasets. CV+PP reaches a comparable accuracy with Exact for all datasets except PPI. NS+PP works slightly worse but the final validation accuracy is still within 2%. These algorithms can be adopted if there is no strong need for predictive performance. We however emphasize that exact algorithms must be used for making predictions, as we will show in Sec. 5.4. Finally, the algorithm NS without preprocessing works much worse than others, indicating the significance of our preprocessing strategy.

5.4. FURTHER ANALYSIS ON TIME COMPLEXITY, TESTING ACCURACY AND VARIANCE:

Table 4 reports the average number of epochs, time, and total number of floating point operations to reach a given 96% validation accuracy on the largest Reddit dataset. Sparse and dense computations are defined in Sec. 4.6. We found that CVD+PP is about 5 times faster than Exact due to the significantly reduced receptive field size. Meanwhile, simply setting D (l) = 2 for NS does not converge to the given accuracy.

We compare the quality of the predictions made by different algorithms, using the same model trained by Exact in Fig. 4 . As Thm. 1 states, CV reaches the same testing accuracy as Exact, while NS and NS+PP perform much worse. Testing using exact algorithms (CV or Exact) corresponds to the weight scaling algorithm for dropout (Srivastava et al., 2014) .

Finally, we compare the average bias and variance of the gradients per dimension for first layer weights relative to the weights' magnitude in Fig. 5 . For models without dropout, the gradient of CV+PP is almost unbiased. For models with dropout, the bias and variance of CV+PP and CVD+PP are ususally smaller than NS and NS+PP, as we analyzed in Sec. 4.3.

6. CONCLUSIONS:

The large receptive field size of GCN hinders its fast stochastic training. In this paper, we present a preprocessing strategy and two control variate based algorithms to reduce the receptive field size.

Our algorithms can achieve comparable convergence speed with the exact algorithm even the neighbor sampling size D (l) = 2, so that the per-epoch cost of training GCN is comparable with training MLPs. We also present strong theoretical guarantees, including exact prediction and convergence to GCN's local optimum, for our control variate based algorithm.

A PROOF OF THEOREM 1:

Proof. 1. We prove by induction. After the first epoch the activation h

(0) i,v is at least computed once for each node v, so H(0) CV,i = H (0)

CV,i = H (0) for all i > I. Assume that we have

H(l) CV,i = H (l)

CV,i = H (l) for all i > (l + 1)I. Then for all i > (l + 1)I Z (l+1)

EQUATION

H (l+1) CV,i = σ(Z (l+1) CV,i ) = H (l+1)

After one more epoch, all the activations h (l+1)

CV,i,v are computed at least once for each v, so H(l+1)

CV,i = H (l+1)

CV,i = H (l+1) for all i > (l + 2)I. By induction, we know that after LI steps, we have -1) . By Eq. 7 we have

H(L-1) CV,i = H (L-1) CV,i = H (L

Z(L) CV,i = Z (L) .

2. We omit the time subscript i and denote

EQUATION

CV,v ). By back propagation, the approximated gradients by CV can be computed as follows

EQUATION

where • is the element wise product and σ (Z

CV ) is the element-wise derivative. Similarly, denote

f v := f (y v , Z (l)

v ), the exact gradients can be computed as follows

EQUATION

Applying E P = E P (1) ,..., P (L) to both sides of Eq. 8, and utilizing

EQUATION

CV = Z (l) , so ∇ Z (L) CV f CV,v is also determinstic. • E P [∇ Z (l) f CV,v ] = E P (l) ,..., P (L) [∇ Z (l) f CV,v ]

.

• E P [∇ H (l) f CV,v ] = E P (l) ,..., P (L) [∇ H (l) f CV,v ].

we have

EQUATION

Comparing Eq. 10 and Eq. 9 we get

E P ∇ W (l) f CV,v = ∇ W (l) f v , l = 0, . . . , L -1, so E P ,V B g CV (W ) = E V B 1 |V B | v∈V B E P ∇ W f CV,v = 1 V v∈V ∇ W f v .

B PROOF OF THEOREM 2:

We proof Theorem 2 in 3 steps:

1. Lemma 1: For a sequence of weights W (1) , . . . , W (N ) which are close to each other, CV's approximate activations are close to the exact activations.

2. Lemma 2: For a sequence of weights W (1) , . . . , W (N ) which are close to each other, CV's gradients are close to be unbiased.

3. Theorem 2: An SGD algorithm generates the weights that changes slow enough for the gradient bias goes to zero, so the algorithm converges.

The following proposition is needed in our proof Proposition 1.

Let A ∞ = max ij |A ij |, then • AB ∞ ≤ col(A) A ∞ B ∞ , where col(A)

is the number of columns of the matrix A.

• A • B ∞ ≤ A ∞ B ∞ . • A + B ∞ ≤ A ∞ + B ∞ .

Proof.

AB ∞ = max ij k A ik B ik ≤ max ij k A ∞ B ∞ = col(A) A ∞ B ∞ . A • B ∞ = max ij |A ij B ij | ≤ max ij A ∞ B ∞ = A ∞ B ∞ . A + B ∞ = max ij |A ij + B ij | ≤ max ij {|A ij | + |B ij |} ≤ max ij |A ij | + max ij |B ij | = A ∞ + B ∞ .

We define C := max{col(P ), col(H (0) ), . . . , col(H (L) )}.

B.1 PROOF OF LEMMA 1:

Proposition 2. There are a series of T inputs X 1 , . . . , X T , X CV,1 , . . . , X CV,T and weights W 1 , . . . , W T feed to an one-layer GCN with CV

Z CV,i = Pi (X i -Xi ) + P Xi W i , H CV,i = σ(Z CV,i ), HCV,i+1 = s i H CV,i +(1-s i ) HCV,i .

and an one-layer exact GCN

Z i = P X i W i , H i = σ(Z i ). If 1. The activation σ(•) is ρ-Lipschitz;

2. X CV,i -X CV,j ∞ < and X CV,i -X i ∞ < for all i, j ≤ T and > 0.

Then there exists some K > 0, s.t., H CV,i -H CV,j ∞ < K and H CV,i -H i ∞ < K for all I < i, j ≤ T , where I is the number of iterations per epoch.

Proof.:

Because for all i > I, the elements of XCV,i are all taken from previous epochs, i.e., X CV,1 , . . . , X CV,i-1 , we know that

EQUATION

By triangular inequality, we also know

EQUATION

)

Since X CV,1 ∞ , . . . , X CV,T ∞ are bounded, XCV,i ∞ is also bounded for i > I. Then,

EQUATION

H CV,i -H i ∞ ≤ ρ Z CV,i -Z i ∞ ≤ ρ Pi (X CV,i -XCV,i ) + P ( XCV,i -X i ) ∞ W i ≤ ρC( Pi ∞ + 2 P ∞ ) W i ∞ ≤ K 2 .

The following lemma bounds CV's approximation error of activations Lemma 1. Given a sequence of model weights W 1 , . . . , W T . If W i -W j ∞ < , ∀i, j, and all the activations are ρ-Lipschitz, there exists K > 0, s.t.,

• H l i -H l CV,i ∞ < K , ∀i > LI, l = 1, . . . , L -1, • Z l i -Z l CV,i ∞ < K , ∀i > LI, l = 1, . . . , L.

Proof. We prove by induction. Because H 0 = X is constant, H0 CV,i = H 0 i after I iterations. So H 1 CV,i = σ( Pi (H 0 CV,i -H0 CV,i ) + P H0 CV,i W 0 i = σ(P XW 0 i ) = H 1 i , and

H 1 CV,i -H 1 CV,j ∞ = σ(P XW 0 i ) -σ(P XW 0 j ) ∞ ≤ ρC P ∞ X ∞ .

Repeatedly apply Proposition B.1 for L -1 times, we get the intended results.

B.2 PROOF OF LEMMA 2:

The following lemma bounds the bias of CV's approximate gradient Lemma 2. Given a sequence of model weights W 1 , . . . , W T , if

1. W i -W j ∞ < , ∀i, j,

2. all the activations are ρ-Lipschitz, 3. the gradient of the cost function ∇ z f (y, z) is ρ-Lipschitz and bounded, then there exists K > 0, s.t.,

E P ,V B g CV (W i ) -g(W i ) ∞ < K , ∀i > LI.

Proof. By Lipschitz continuity of ∇ z f (y, z) and Lemma 1, there exists K > 0, s.t.,

∇ Z (l) CV f CV,v -∇ Z (l) f v ∞ ≤ ρ Z (l) CV -Z (l) ∞ ≤ ρK .

Assume that

E P ∇ Z (l+1) CV f CV,v -∇ Z (l+1) f v ∞

< K 1 , we now prove that there exists K > 0, s.t.,

E P ∇ Z (l) CV f CV,v -∇ Z (l) f v ∞ < K

. By Eq. 9, Eq. 10 and Lemma 1, we have

EQUATION

E P ∇ Z (l) CV f CV,v -∇ Z (l) f v ∞ = E P σ (Z (l) CV ) • ∇ H (l) CV f CV,v -σ (Z (l) ) • ∇ H (l) f v ∞ ≤ E P (σ (Z (l) CV ) -σ (Z (l) )) • ∇ H (l) CV f CV,v ∞ + σ (Z (l) )(E P [∇ H (l) CV f CV,v ] -∇ H (l) f v ) ∞ ≤ E P ρKC • ∇ H (l) CV f CV,v ] ∞ + σ (Z (l) ) ∞ P ∞ K 1 C 3 W (l) ∞ ≤ρKC 2 E P ∇ H (l) CV f CV,v ∞ + σ (Z (l) ) ∞ P ∞ K 1 C 3 W (l) ∞ ≤ K 2

By induction we know that for l = 1, . . . , L there exists K, s.t.,

E P ∇ Z (l) CV f CV,v -∇ Z (l) f v ∞ ≤ K .

Again by Eq. 9, Eq. 10, and Lemma 1,

E P ∇ W (l) f CV,v -∇ W (l) f v ∞ = E P H (l) CV P ∇ Z (l) CV f CV,v -H (l) P ∇ Z (l) f v ∞ ≤ E P H (l) CV -H (l) P ∇ Z (l) CV f CV,v ∞ + H (l) P E P ∇ Z (l) CV f CV,v -∇ Z (l) f v ∞ ≤ E P KC P ∇ Z (l) CV f CV,v ∞ + C 2 H (l) ∞ P ∞ K ≤KC 2 P ∞ E P ∇ Z (l) CV f CV,v ∞ + H (l) ∞ P ∞ KC 2 ≤ K 3 Finally, E P ,V B g CV (W i ) -g(W i ) ∞ = E V B 1 |V B | v∈V B E P [∇ W (l) f CV,v ] - 1 V v∈V ∇ W (l) f v ∞ = 1 V v∈V E P [∇ W (l) f CV,v ] -∇ W (l) f v ∞ ≤ K 3 .

B.3 PROOF OF THEOREM 2:

Proof. This proof is a modification of Ghadimi & Lan (2013) , but using biased stochastic gradients instead. We assume the algorithm is already warmed-up for LI steps with the initial weights W 0 , so that Lemma 2 holds for step i > 0. Denote δ i = g CV (W i ) -∇L(W i ). By smoothness we have

EQUATION

)

Consider the sequence of LI + 1 weights W i-LI , . . . , W i .

max i-LI≤j,k≤i W j -W k ∞ ≤ i-1 j=i-LI W j -W j+1 ∞ = i-1 j=i-LI γ j g CV (W j ) ∞ ≤ i-1 j=i-LI γ j G ≤ LIGγ i-LI .

By Lemma 2, there exists K > 0, s.t.

E P ,V B δ i ∞ = E P ,V B g CV (W i ) -∇L(W i ) ∞ ≤ KLIGγ i-LI , ∀i > 0. Assume that W is D-dimensional, E P ,V B ∇L(W i ), δ i ≥ -E P ,V B D ∇L(W i ) ∞ δ i ∞ ≥ -KLIDG 2 γ i-LI = K 1 γ i-LI , E P ,V B δ i 2 ≤ D E P ,V B δ i ∞ 2 ≤ DK 2 L 2 B 2 G 2 γ i-LI = K 2 γ i-LI ,

where K 1 = KLIDG 2 and K 2 = DK 2 L 2 B 2 G 2 . Taking E P ,V B to both sides of Eq. 14 we have

L(W i+1 ) ≤ L(W i ) + (γ i -ργ 2 i )K 1 γ i-LI -(γ i - ργ 2 i 2 )E P ,V B ∇L(W i ) 2 + ρ 2 γ 2 i K 2 γ i-LI .

Summing up the above inequalities and re-arranging the terms, we obtain,

N i=1 (γ i - ργ 2 i 2 )E P ,V B ∇L(W i ) 2 ≤L(W 1 ) -L(W * ) + K 1 N i=1 (γ i -ργ 2 i )γ i-LI + ρK 2 2 N i=1 γ 2 i γ i-LI .

Dividing both sides by

N i=1 (γ i - ργ 2 i 2 ), E R∼P R E P ,V B ∇L(W R ) 2 ≤2 L(W 1 ) -L(W * ) + K 1 N i=1 (γ i -ργ 2 i )γ i-LI + ρK2 2 N i=1 γ 2 i γ i-LI N i=1 γ i (2 -ργ i )

.

Taking γ i = γ := min{ 1 ρ , 1 √ N }, for all i = 1, . . . , N , we have

E R∼P R E P ,V B ∇L(W R ) 2 ≤2 L(W 1 ) -L(W * ) + K 1 N (γ -ργ 2 )γ + ρK2 2 N γ 3 N γ(2 -ργ) ≤2 L(W 1 ) -L(W * ) + K 1 N (γ -ργ 2 )γ + ρK2 2 N γ 3 N γ ≤2 L(W 1 ) -L(W * ) N γ + K 1 γ(1 -ργ) + ρK 2 γ 2 ≤ 2ρ(L(W 1 ) -L(W * )) + K 2 N + 2(L(W 1 ) -L(W * )) + K 1 √ N .

Particularly, when N → ∞, we have E R∼P R E P ,V B ∇L(W R ) 2 = 0, which implies that the gradient is asymptotically unbiased. Training with the CVD estimator is similar with the CV estimator, except it runs two versions of the network, with and without dropout, to compute the samples H and their mean µ of the activation. The matrix Pv,v = Pv,v / |n(v, 1)| , where |n(v, 1)| is the degree of node v.

C DERIVATION OF THE VARIANCE Var:

[u] =E v p v (h v -µ v ) 2 = v p 2 v E [(h v -µ v )] 2 = v p 2 v s 2 v =s 2 . Var[u N S ] =E [Dp v h v -µ] 2 =E v {D 2 p 2 v (µ 2 v + s 2 v ) + µ 2 -Dµp v µ vp } =Ds 2 + (D v p 2 v µ 2 v -µ 2 ) =Ds 2 + 1 2 v,v (p v µ v -p v µ v ) 2 .

F EXPERIMENT FOR 3-LAYER GCNS:

We test 3-layer GCNs on the Reddit dataset. The settings are the same with 2-layer GCNs in Sec. 5.3.

To ensure the exact algorithm can run in a reasonable amount of time, we subsample the graph so that the maximum degree is 10. The convergence result is shown as Fig. 6 , which is similar with the two-layer models. The time consumption to reach 0.94 testing accuracy is shown in Table 5 .

G JUSTIFICATION OF THE INDEPENDENT GAUSSIAN ASSUMPTION G.1 RESULTS FOR 2-LAYER GCNS:

We justify the independent Gaussian assumption in Sec. 4.3 by showing that for a 2-layer GCN with the first layer pre-processed, the neighbor's activations are independent. Without loss of generality, suppose that we want to compute z

(2)

1 , and the neighbors of node 1 are 1, . . . , D. By Eq. (1), h where the indices i, j ∈ n(l) (v). Then, we compute the average correlation of all pairs of neighbors i = j.

(1) v = σ (φ v • u (0) v )W (0)

AvgCorr (l,v,d) := 1 n(l) (v) ( n(l) (v) -1) i =j Corr (l,v,d) ij

, and define the average neighbor correlation on layer l as AvgCorr (l,v,d) averaged over all the nodes v and dimensions d.

We report the average feature correlation and the average neighbor correlation per layer, on the Citeseer, Cora, PubMed and PPI datasets. These quantities are too expensive to compute for NELL and Reddit. On each dataset, we train a GCN with 10 graph convoluation layers until early stopping criteria is met, and compute the average feature correlation and the average neighbor correlation for layer 1 to 9. We are not interested in the correlation on layer 10 because there are no more graph convolutional layers after it. The result is shown as Fig. 7 . As analyzed in Sec. G.1, the average neighbor correlation is close to zero on the first layer, but it is not exactly zero due to the finite sample size for computing the empirical covariance. There is no strong tendency of increased correlation as the number of layers increases, after the third layer. The average neighbor correlation and the average feature correlation remain on the same order of magnitude, so bringing correlated neighbors does not make the activations much more correlated than the MLP case (Wang & Manning, 2013) . Finally, both correlations are much smaller than one.

:

For an integer N , we define [N ] = {0, . . . , N } and [N ]+ = {1, . . . , N }.

We will omit the subscripts CV and i in subsequent sections when there is no confusion.

Figures and Tables:

Figure 1: Two-layer graph convolutional networks, and the receptive field of a single vertex.

the batch GCN has O(EK) sparse and O(V KA) dense time complexity per epoch. NS has O(V L l=1 D (l) K) sparse and O(V L l=2 D (l) KA) dense time complexity per epoch. The dense time complexity of CV is the same as NS. The sparse time complexity depends on the cost of computing the sum v p v μv . There are V L l=2 D (l) such sums to compute on the first graph convolution layer, and overall cost is not larger than O(V D L l=2 D (l) K), if we subsample the graph such that the max degree is D, following Hamilton et al. (2017a). The sparse time complexity is D/D (1) times higher than NS.

Figure 2: Comparison of training loss with respect to number of epochs without dropout. The CV+PP curve overlaps with the Exact curve in the first four datasets.

Figure 3: Comparison of validation accuracy with respect to number of epochs. NS converges to 0.94 on the Reddit dataset and 0.6 on the PPI dataset.

Figure 4: Comparison of the accuracy of different testing algorithms. The yaxis is Micro-F1 for PPI and accuracy otherwise.

Figure 6: Comparison of validation accuracy with respect to number of epochs for 3-layer GCNs.

is a random variable with respect to φ v , where φ v ∼ Bernoulli(p) is the dropout mask and u(0) v = (P H (0) ) v . The indepedent Gaussian assumption states that 1are independent, for v = v .Assumption 1 is not GCN-specific and is discussed inWang & Manning (2013), we now prove assumption 2 by the following lemma.Lemma 3. If a and b are independent random variables, then their transformations f 1 (a) and f 2 (b) are independent.Because for any event A and B, P(f 1 (a) ∈ f 1 (A), f 2 (b) ∈ f 2 (B)) = P (a ∈ A, b ∈ B) = P (a ∈ A)P (b ∈ B) = P (f 1 (a) ∈ f 1 (A))P (f 2 (B) ∈ f 2 (B)),where f 1 (A) = {f 1 (a)|a ∈ A} and f 2 (B) = {f 2 (b)|b ∈ B}.

To keep notations simple, we omit all the subscripts and tildes, and exchange the ID of nodes such Number of vertexes, edges, average number of 1-and 2-neighbors per node for each dataset. Undirected edges are counted twice and self-loops are counted once. Reddit is already subsampled to have a max degree of 128 followingHamilton et al. (2017a).

Table content: <table><tr><td/><td>Dataset</td><td>V</td><td>E</td><td colspan="2">Degree Degree 2</td><td>Type</td></tr><tr><td/><td>Citeseer</td><td>3,327</td><td>12,431</td><td>4</td><td>15</td><td>Document network</td></tr><tr><td/><td>Cora</td><td>2,708</td><td>13,264</td><td>5</td><td>37</td><td>Document network</td></tr><tr><td/><td>PubMed</td><td>19,717</td><td>108,365</td><td>6</td><td>60</td><td>Document network</td></tr><tr><td/><td>NELL</td><td>65,755</td><td>318,135</td><td>5</td><td>1,597</td><td>Knowledge graph</td></tr><tr><td/><td>PPI</td><td>14,755</td><td>458,973</td><td>31</td><td>970</td><td>Protein-protein interaction</td></tr><tr><td/><td>Reddit</td><td colspan="2">232,965 23,446,803</td><td>101</td><td>10,858</td><td>Document network</td></tr><tr><td colspan="7">that n(v, 1) = [D] + , 1 where D = |n(v, 1)| is the number of neighbors. We get the propagation</td></tr><tr><td>rule u =</td><td colspan="6">D v=1 p v h v , which is used interchangeably with the matrix form U (l) = P H(l) .</td></tr><tr><td colspan="3">2.3 NEIGHBOR SAMPLING</td><td/><td/><td/></tr><tr><td colspan="7">To reduce the receptive field size, Hamilton et al. (2017a) propose a neighbor sampling (NS) al-</td></tr><tr><td colspan="7">gorithm. On the l-th layer, they randomly choose D (l) neighbors for each node, and develop an</td></tr><tr><td colspan="7">estimator u N S of u based on Monte-Carlo approximation u ≈ u N S = D D (l)</td></tr></table>

Time complexity comparison of different algorithms on the Reddit dataset.

Table content: <table><tr><td>Alg.</td><td colspan="2">Valid. Epochs acc.</td><td>Time (s)</td><td colspan="2">Sparse GFLOP TFLOP Dense</td></tr><tr><td>Exact</td><td>96.0</td><td>4.2</td><td>252</td><td>507</td><td>7.17</td></tr><tr><td>NS</td><td>94.4</td><td>102.0</td><td>577</td><td>76.5</td><td>21.4</td></tr><tr><td>NS+PP</td><td>96.0</td><td>35.0</td><td>195</td><td>2.53</td><td>7.36</td></tr><tr><td>CV+PP</td><td>96.0</td><td>7.8</td><td>56</td><td>40.6</td><td>1.64</td></tr><tr><td>CVD+PP</td><td>96.0</td><td>5.8</td><td>50</td><td>60.3</td><td>2.44</td></tr></table>

Time to reach 0.95 testing accuracy.

Table content: <table><tr><td>Alg.</td><td colspan="2">Valid. Epochs acc.</td><td colspan="3">Time Sparse (s) GFLOP TFLOP Dense</td></tr><tr><td>Exact</td><td>0.940</td><td>3.0</td><td>199</td><td>306</td><td>11.7</td></tr><tr><td>NS</td><td>0.940</td><td>24.0</td><td>148</td><td>33.6</td><td>9.79</td></tr><tr><td>NS+PP</td><td>0.940</td><td>12.0</td><td>68</td><td>2.53</td><td>4.89</td></tr><tr><td>CV+PP</td><td>0.940</td><td>5.0</td><td>32</td><td>8.06</td><td>2.04</td></tr><tr><td colspan="2">CVD+PP 0.940</td><td>5.0</td><td>36</td><td>16.1</td><td>4.08</td></tr><tr><td colspan="3">E.2 TRAINING WITH THE CVD ESTIMATOR</td><td/><td/><td/></tr></table>

Peer review (iclr2018):
 This paper proposes a new training method for graph convolutional networks. The experimental results look interesting. However, this paper has some issues. This paper is hard to read. There are some undefined or multi-used notations. For instance, sigma is used for two different meanings: an activation function and variance. Some details that need to be explained are omitted. For example, what kind of dropout is used to obtain the table and figures in Section 5? Forward and backward propagation processes are not clearly explained In section 4.2, it is not clear why we have to multiply sqrt{D}. Why should we make the variance from dropout sigma^2? Proposition 1 is wrong. First, .|A.|_.infty should be max_{ij} |A_ij| not A_{ij}. Second, there is no order between .|AB.|_.infty and .|A.|_.infty .|B.|_.infty. When A=[1 1] and B is the transpose matrix of A, .|AB.|_.infty =2 and .|A.|_.infty .|B.|_.infty = 1. When, A’=[1 -1] and B is the same matrix defined just before, .|A’ B .|_.infty = 0 and .|A’.|_.infty .|B.|_.infty =1. So, both .|AB.|_.infty .le .|A.|_.infty .|B.|_.infty and .|AB.|_.infty .ge .|A.|_.infty .|B.|_.infty are not true. I cannot believe the proof of Theorem 2.