Title: FASTGCN: FAST LEARNING WITH GRAPH CONVOLU-TIONAL NETWORKS VIA IMPORTANCE SAMPLING

Abstract: The graph convolutional networks (GCN) recently proposed by Kipf and Welling are an effective graph model for semi-supervised learning. This model, however, was originally designed to be learned with the presence of both training and test data. Moreover, the recursive neighborhood expansion across layers poses time and memory challenges for training with large, dense graphs. To relax the requirement of simultaneous availability of test data, we interpret graph convolutions as integral transforms of embedding functions under probability measures. Such an interpretation allows for the use of Monte Carlo approaches to consistently estimate the integrals, which in turn leads to a batched training scheme as we propose in this work-FastGCN. Enhanced with importance sampling, FastGCN not only is efficient for training but also generalizes well for inference. We show a comprehensive set of experiments to demonstrate its effectiveness compared with GCN and related models. In particular, training is orders of magnitude more efficient while predictions remain comparably accurate. * These two authors contribute equally.

1. INTRODUCTION:

Graphs are universal representations of pairwise relationship. Many real world data come naturally in the form of graphs; e.g., social networks, gene expression networks, and knowledge graphs. To improve the performance of graph-based learning tasks, such as node classification and link prediction, recently much effort is made to extend well-established network architectures, including recurrent neural networks (RNN) and convolutional neural networks (CNN), to graph data; see, e.g., Bruna et al. (2013) ; Duvenaud et al. (2015) ; Li et al. (2015) ; Jain et al. (2015) ; Henaff et al. (2015) ; Niepert et al. (2016) ; Kipf & Welling (2016a; b) .

Whereas learning feature representations for graphs is an important subject among this effort, here, we focus on the feature representations for graph vertices. In this vein, the closest work that applies a convolution architecture is the graph convolutional network (GCN) (Kipf & Welling, 2016a; b) . Borrowing the concept of a convolution filter for image pixels or a linear array of signals, GCN uses the connectivity structure of the graph as the filter to perform neighborhood mixing. The architecture may be elegantly summarized by the following expression:

H (l+1) = σ( ÂH (l) W (l) ),

where Â is some normalization of the graph adjacency matrix, H (l) contains the embedding (rowwise) of the graph vertices in the lth layer, W (l) is a parameter matrix, and σ is nonlinearity.

As with many graph algorithms, the adjacency matrix encodes the pairwise relationship for both training and test data. The learning of the model as well as the embedding is performed for both data simultaneously, at least as the authors proposed. For many applications, however, test data may not be readily available, because the graph may be constantly expanding with new vertices (e.g. new members of a social network, new products to a recommender system, and new drugs for functionality tests). Such scenarios require an inductive scheme that learns a model from only a training set of vertices and that generalizes well to any augmentation of the graph.

A more severe challenge for GCN is that the recursive expansion of neighborhoods across layers incurs expensive computations in batched training. Particularly for dense graphs and powerlaw graphs, the expansion of the neighborhood for a single vertex quickly fills up a large portion of the graph. Then, a usual mini-batch training will involve a large amount of data for every batch, even with a small batch size. Hence, scalability is a pressing issue to resolve for GCN to be applicable to large, dense graphs.

To address both challenges, we propose to view graph convolutions from a different angle and interpret them as integral transforms of embedding functions under probability measures. Such a view provides a principled mechanism for inductive learning, starting from the formulation of the loss to the stochastic version of the gradient. Specifically, we interpret that graph vertices are iid samples of some probability distribution and write the loss and each convolution layer as integrals with respect to vertex embedding functions. Then, the integrals are evaluated through Monte Carlo approximation that defines the sample loss and the sample gradient. One may further alter the sampling distribution (as in importance sampling) to reduce the approximation variance.

The proposed approach, coined FastGCN, not only rids the reliance on the test data but also yields a controllable cost for per-batch computation. At the time of writing, we notice a newly published work GraphSAGE (Hamilton et al., 2017) that proposes also the use of sampling to reduce the computational footprint of GCN. Our sampling scheme is more economic, resulting in a substantial saving in the gradient computation, as will be analyzed in more detail in Section 3.3. Experimental results in Section 4 indicate that the per-batch computation of FastGCN is more than an order of magnitude faster than that of GraphSAGE, while classification accuracies are highly comparable.

2. RELATED WORK:

Over the past few years, several graph-based convolution network models emerged for addressing applications of graph-structured data, such as the representation of molecules (Duvenaud et al., 2015) . An important stream of work is built on spectral graph theory (Bruna et al., 2013; Henaff et al., 2015; Defferrard et al., 2016) . They define parameterized filters in the spectral domain, inspired by graph Fourier transform. These approaches learn a feature representation for the whole graph and may be used for graph classification.

Another line of work learns embeddings for graph vertices, for which Goyal & Ferrara (2017) is a recent survey that covers comprehensively several categories of methods. A major category consists of factorization based algorithms that yield the embedding through matrix factorizations; see, e.g., Roweis & Saul (2000) ; Belkin & Niyogi (2001) ; Ahmed et al. (2013) ; Cao et al. (2015) ; Ou et al. (2016) . These methods learn the representations of training and test data jointly. Another category is random walk based methods (Perozzi et al., 2014; Grover & Leskovec, 2016 ) that compute node representations through exploration of neighborhoods. LINE (Tang et al., 2015) is also such a technique that is motivated by the preservation of the first and second-order proximities. Meanwhile, there appear a few deep neural network architectures, which better capture the nonlinearity within graphs, such as SDNE (Wang et al., 2016) . As motivated earlier, GCN (Kipf & Welling, 2016a) is the model on which our work is based.

The most relevant work to our approach is GraphSAGE (Hamilton et al., 2017) , which learns node representations through aggregation of neighborhood information. One of the proposed aggregators employs the GCN architecture. The authors also acknowledge the memory bottleneck of GCN and hence propose an ad hoc sampling scheme to restrict the neighborhood size. Our sampling approach is based on a different and more principled formulation. The major distinction is that we sample vertices rather than neighbors. The resulting computational savings are analyzed in Section 3.3.

3. TRAINING AND INFERENCE THROUGH SAMPLING:

One striking difference between GCN and many standard neural network architectures is the lack of independence in the sample loss. Training algorithms such as SGD and its batch generalization are designed based on the additive nature of the loss function with respect to independent data samples. For graphs, on the other hand, each vertex is convolved with all its neighbors and hence defining a sample gradient that is efficient to compute is beyond straightforward.

Concretely, consider the standard SGD scenario where the loss is the expectation of some function g with respect to a data distribution D:

L = E x∼D [g(W ; x)].

Here, W denotes the model parameter to be optimized. Of course, the data distribution is generally unknown and one instead minimizes the empirical loss through accessing n iid samples x 1 , . . . , x n :

L emp = 1 n n i=1 g(W ; x i ), x i ∼ D, ∀ i.

In each step of SGD, the gradient is approximated by ∇g(W ; x i ), an (assumed) unbiased sample of ∇L. One may interpret that each gradient step makes progress toward the sample loss g(W ; x i ).

The sample loss and the sample gradient involve only one single sample x i .

For graphs, one may no longer leverage the independence and compute the sample gradient ∇g(W ; x i ) by discarding the information of i's neighboring vertices and their neighbors, recursively. We therefore seek an alternative formulation. In order to cast the learning problem under the same sampling framework, let us assume that there is a (possibly infinite) graph G with the vertex set V associated with a probability space (V , F, P ), such that for the given graph G, it is an induced subgraph of G and its vertices are iid samples of V according to the probability measure P . For the probability space, V serves as the sample space and F may be any event space (e.g., the power set F = 2 V ). The probability measure P defines a sampling distribution.

To resolve the problem of lack of independence caused by convolution, we interpret that each layer of the network defines an embedding function of the vertices (random variable) that are tied to the same probability measure but are independent. See Figure 1 . Specifically, recall the architecture of GCN

H(l+1) = ÂH (l) W (l) , H (l+1) = σ( H(l+1) ), l = 0, . . . , M -1, L = 1 n n i=1

g(H (M ) (i, :)).

(1) For the functional generalization, we write

EQUATION

EQUATION

Here, u and v are independent random variables, both of which have the same probability measure P . The function h (l) is interpreted as the embedding function from the lth layer. The embedding functions from two consecutive layers are related through convolution, expressed as an integral transform, where the kernel Â(v, u) corresponds to the (v, u) element of the matrix Â. The loss is the expectation of g(h (M ) ) for the final embedding h (M ) . Note that the integrals are not the usual Riemann-Stieltjes integrals, because the variables u and v are graph vertices but not real numbers; however, this distinction is only a matter of formalism.

Writing GCN in the functional form allows for evaluating the integrals in the Monte Carlo manner, which leads to a batched training algorithm and also to a natural separation of training and test data, as in inductive learning. For each layer l, we use t l iid samples u

EQUATION

t l ∼ P to approximately evaluate the integral transform (2); that is, 0) . Then, the loss L in (3) admits an estimator

EQUATION

t l+1 (v) := σ( h(l+1) t l+1 (v)), l = 0, . . . , M -1, with the convention h (0) t0 ≡ h (

L t0,t1,...,t M := 1 t M t M i=1 g(h (M ) t M (u (M ) i )).

The follow result establishes that the estimator is consistent. The proof is a recursive application of the law of large numbers and the continuous mapping theorem; it is given in the appendix.

Graph convolution view:

Integral transform view batch

H (2)

H ( 1)

H (0) h (2) (v) h (1) (v) h (0) (v)

Figure 1 : Two views of GCN. On the left (graph convolution view), each circle represents a graph vertex. On two consecutive rows, a circle i is connected (in gray line) with circle j if the two corresponding vertices in the graph are connected. A convolution layer uses the graph connectivity structure to mix the vertex features/embeddings. On the right (integral transform view), the embedding function in the next layer is an integral transform (illustrated by the orange fanout shape) of the one in the previous layer. For the proposed method, all integrals (including the loss function) are evaluated by using Monte Carlo sampling. Correspondingly in the graph view, vertices are subsampled in a bootstrapping manner in each layer to approximate the convolution. The sampled portions are collectively denoted by the solid blue circles and the orange lines.

Theorem 1. If g and σ are continuous, then lim t0,t1,...,t M →∞ L t0,t1,...,t M = L with probability one.

In practical use, we are given a graph whose vertices are already assumed to be samples. Hence, we will need bootstrapping to obtain a consistent estimate. In particular, for the network architecture (1), the output H (M ) is split into batches as usual. We will still use u

(M ) 1 , . . . , u (M )

t M to denote a batch of vertices, which come from the given graph. For each batch, we sample (with replacement) uniformly each layer and obtain samples u (l) i , i = 1, . . . , t l , l = 0, . . . , M -1. Such a procedure is equivalent to uniformly sampling the rows of H (l) for each l. Then, we obtain the batch loss

EQUATION

where, recursively,

EQUATION

Here, the n inside the activation function σ is the number of vertices in the given graph and is used to account for the normalization difference between the matrix form (1) and the integral form (2). The corresponding batch gradient may be straightforwardly obtained through applying the chain rule on each H (l) . See Algorithm 1.

3.1. VARIANCE REDUCTION:

As for any estimator, one is interested in improving its variance. Whereas computing the full variance is highly challenging because of nonlinearity in all the layers, it is possible to consider each single layer and aim at improving the variance of the embedding function before nonlinearity. Specifically, consider for the lth layer, the function h(l+1) t l+1 (v) as an approximation to the convolution Â(v, u)h

(l) t l (u)W (l) dP (u). When taking t l+1 samples v = u (l+1) 1 , . . . , u (l+1)

t l+1 , the sample average of h(l+1) t l+1 (v) admits a variance that captures the deviation from the eventual loss contributed by this layer. Hence, we seek an improvement of this variance. Now that we consider each layer separately, we will do the following change of notation to keep the expressions less cumbersome: If v is sampled in the next layer,

∇ H(l+1) (v, :) ← n t l t l j=1 Â(v, u (l) j )∇ H (l) (u (l) j , :)W (l)

5:

end for 6:

W ← W -η∇L batch SGD step 7: end for Function Samples Num. samples Layer l + 1; random variable v h(l+1) t l+1 (v) → y(v) u (l+1) i → v i t l+1 → s Layer l; random variable u h (l) t l (u)W (l) → x(u) u (l) j → u j t l → t

Under the joint distribution of v and u, the aforementioned sample average is

G := 1 s s i=1 y(v i ) = 1 s s i=1   1 t t j=1 Â(v i , u j )x(u j )   .

First, we have the following result. Proposition 2. The variance of G admits

EQUATION

where

R = 1 s 1 - 1 t e(v) 2 dP (v) - 1 s e(v) dP (v)

2 and e(v) = Â(v, u)x(u) dP (u).

The variance (6) consists of two parts. The first part R leaves little room for improvement, because the sampling in the v space is not done in this layer. The second part (the double integral), on the other hand, depends on how the u j 's in this layer are sampled. The current result ( 6) is the consequence of sampling u j 's by using the probability measure P . One may perform importance sampling, altering the sampling distribution to reduce variance. Specifically, let Q(u) be the new probability measure, where the u j 's are drawn from. We hence define the new sample average approximation

y Q (v) := 1 t t j=1 Â(v, u j )x(u j ) dP (u) dQ(u) uj , u 1 , . . . , u t ∼ Q,

and the quantity of interest

G Q := 1 s s i=1 y Q (v i ) = 1 s s i=1   1 t t j=1 Â(v i , u j )x(u j ) dP (u) dQ(u) uj   .

Clearly, the expectation of G Q is the same as that of G, regardless of the new measure Q. The following result gives the optimal Q. Theorem 3. If

EQUATION

)

then the variance of G Q admits

EQUATION

)

where R is defined in Proposition 2. The variance is minimum among all choices of Q.

A drawback of defining the sampling distribution Q in this manner is that it involves |x(u)|, which constantly changes during training. It corresponds to the product of the embedding matrix H (l) and the parameter matrix W (l) . The parameter matrix is updated in every iteration; and the matrix product is expensive to compute. Hence, the cost of computing the optimal measure Q is quite high.

As a compromise, we consider a different choice of Q, which involves only b(u). The following proposition gives the precise definition. The resulting variance may or may not be smaller than (6).

In practice, however, we find that it is almost always helpful. Proposition 4. If

dQ(u) = b(u) 2 dP (u) b(u) 2 dP (u)

where b(u) is defined in (7), then the variance of G Q admits

EQUATION

)

where R is defined in Proposition 2.

With this choice of the probability measure Q, the ratio dQ(u)/dP (u) is proportional to b(u) 2 , which is simply the integral of Â(v, u) 2 with respect to v. In practical use, for the network architecture (1), we define a probability mass function for all the vertices in the given graph:

q(u) = Â(:, u) 2 / u ∈V Â(:, u ) 2 , u ∈ V

and sample t vertices u 1 , . . . , u t according to this distribution. From the expression of q, we see that it has no dependency on l; that is, the sampling distribution is the same for all layers. To summarize, the batch loss L batch in (4) now is recursively expanded as

H (l+1) (v, :) = σ   1 t l t l j=1 Â(v, u (l) j )H (l) (u (l) j , :)W (l) q(u (l) j )   , u (l) j ∼ q, l = 0, . . . , M -1. (10)

The major difference between (5) and ( 10) is that the former obtains samples uniformly whereas the latter according to q. Accordingly, the scaling inside the summation changes. The corresponding batch gradient may be straightforwardly obtained through applying the chain rule on each H (l) . See Algorithm 2.

Algorithm 2 FastGCN batched training (one epoch), improved version:

1: For each vertex u, compute sampling probability q(u) ∝ Â(:, u) 2 2: for each batch do If v is sampled in the next layer,

∇ H(l+1) (v, :) ← 1 t l t l j=1 Â(v, u (l) j ) q(u (l) j ) ∇ H (l) (u (l) j , :)W (l)

6: end for 7:

W ← W -η∇L batch SGD step 8: end for

3.2. INFERENCE:

The sampling approach described in the preceding subsection clearly separates out test data from training. Such an approach is inductive, as opposed to transductive that is common for many graph algorithms. The essence is to cast the set of graph vertices as iid samples of a probability distribution, so that the learning algorithm may use the gradient of a consistent estimator of the loss to perform parameter update. Then, for inference, the embedding of a new vertex may be either computed by using the full GCN architecture (1), or approximated through sampling as is done in parameter learning. Generally, using the full architecture is more straightforward and easier to implement.

3.3. COMPARISON WITH GRAPHSAGE:

GraphSAGE (Hamilton et al., 2017 ) is a newly proposed architecture for generating vertex embeddings through aggregating neighborhood information. It shares the same memory bottleneck with GCN, caused by recursive neighborhood expansion. To reduce the computational footprint, the authors propose restricting the immediate neighborhood size for each layer. Using our notation for the sample size, if one samples t l neighbors for each vertex in the lth layer, then the size of the expanded neighborhood is, in the worst case, the product of the t l 's. On the other hand, FastGCN samples vertices rather than neighbors in each layer. Then, the total number of involved vertices is at most the sum of the t l 's, rather than the product. See experimental results in Section 4 for the order-of-magnitude saving in actual computation time.

4. EXPERIMENTS:

We follow the experiment setup in Kipf & Welling (2016a) and Hamilton et al. (2017) to demonstrate the effective use of FastGCN, comparing with the original GCN model as well as Graph-SAGE, on the following benchmark tasks: (1) classifying research topics using the Cora citation data set (McCallum et al., 2000) ; (2) categorizing academic papers with the Pubmed database; and

(3) predicting the community structure of a social network modeled with Reddit posts. These data sets are downloaded from the accompany websites of the aforementioned references. The graphs have increasingly more nodes and higher node degrees, representative of the large and dense setting under which our method is motivated. Statistics are summarized in Table 1 . We adjusted the training/validation/test split of Cora and Pubmed to align with the supervised learning scenario. Specifically, all labels of the training examples are used for training, as opposed to only a small portion in the semi-supervised setting (Kipf & Welling, 2016a) . Such a split is coherent with that of the other data set, Reddit, used in the work of GraphSAGE. Additional experiments using the original split of Cora and Pubmed are reported in the appendix. We first consider the use of sampling in FastGCN. The left part of Table 2 (columns under "Sampling") lists the time and classification accuracy as the number of samples increases. For illustration purpose, we equalize the sample size on both layers. Clearly, with more samples, the per-epoch training time increases, but the accuracy (as measured by using micro F1 scores) also improves generally.

An interesting observation is that given input features H (0) , the product ÂH (0) in the bottom layer does not change, which means that the chained expansion of the gradient with respect to W (0) in Next, we compare the sampling approaches for FastGCN: uniform and importance sampling. Figure 2 summarizes the prediction accuracy under both approaches. It shows that importance sampling consistently yields higher accuracy than does uniform sampling. Since the altered sampling distribution (see Proposition 4 and Algorithm 2) is a compromise alternative of the optimal distribution that is impractical to use, this result suggests that the variance of the used sampling indeed is smaller than that of uniform sampling; i.e., the term (9) stays closer to (8) than does (6). A possible reason is that b(u) correlates with |x(u)|. Hence, later experiments will apply importance sampling.

We now demonstrate that the proposed method is significantly faster than the original GCN as well as GraphSAGE, while maintaining comparable prediction performance. See Figure 3 Hamilton et al. (2017) . The timings of using other aggregators, such as GraphSAGE-mean, are similar. GCN refers to using batched learning, as opposed to the original version that is nonbatched; for more details of the implementation, see the appendix. The nonbatched version of GCN runs out of memory on the large graph Reddit. The sample sizes for FastGCN are 400, 100, and 400, respectively for the three data sets.

In the discussion period, the authors of GraphSAGE offered an improved implementation of their codes and alerted that GraphSAGE was better suited for massive graphs. The reason is that for small graphs, the sample size (recalling that it is the product across layers) is comparable to the graph size and hence improvement is marginal; moreover, sampling overhead might then adversely affect the timing. For fair comparison, the authors of GraphSAGE kept the sampling strategy but improved the implementation of their original codes by eliminating redundant calculations of the sampled nodes. Now the per-batch training time of GraphSAGE compares more favorably on the smallest graph Cora; see Table 3 . Note that this implementation does not affect large graphs (e.g., Reddit) and our observation of orders of magnitude faster training remains valid. 

5. CONCLUSIONS:

We have presented FastGCN, a fast improvement of the GCN model recently proposed by Kipf & Welling (2016a) for learning graph embeddings. It generalizes transductive training to an inductive manner and also addresses the memory bottleneck issue of GCN caused by recursive expansion of neighborhoods. The crucial ingredient is a sampling scheme in the reformulation of the loss and the gradient, well justified through an alternative view of graph convoluntions in the form of integral transforms of embedding functions. We have compared the proposed method with additionally GraphSAGE (Hamilton et al., 2017) , a newly published work that also proposes using sampling to restrict the neighborhood size, although the two sampling schemes substantially differ in both algorithm and cost. Experimental results indicate that our approach is orders of magnitude faster than GCN and GraphSAGE, while maintaining highly comparable prediction performance with the two.

The simplicity of the GCN architecture allows for a natural interpretation of graph convolutions in terms of integral transforms. Such a view, yet, generalizes to many graph models whose formulations are based on first-order neighborhoods, examples of which include MoNet that applies to (meshed) manifolds (Monti et al., 2017) , as well as many message-passing neural networks (see e.g., Scarselli et al. (2009) ; Gilmer et al. (2017) ). The proposed work elucidates the basic Monte Carlo ingredients for consistently estimating the integrals. When generalizing to other networks aforementioned, an additional effort is to investigate whether and how variance reduction may improve the estimator, a possibly rewarding avenue of future research. KDD '16, pp. 1225 -1234 , 2016. ISBN 978-1-4503-4232-2 .

A PROOFS:

Proof of Theorem 1. Because the samples u (0) j are iid, by the strong law of large numbers,

h(1) t1 (v) = 1 t 0 t0 j=1 Â(v, u (0) j )h (0) (u (0) j )W (0)

converges almost surely to h(1) (v). Then, because the activation function σ is continuous, the continuous mapping theorem implies that h

(1) 1) dP (u), where note that the probability space is with respect to the 0th layer and hence has nothing to do with that of the variable u or v in this statement. Similarly,

t1 (v) = σ( h(1) t1 (v)) converges almost surely to h (1) (v) = σ( h(1) (v)). Thus, Â(v, u)h (1) t1 (u)W (1) dP (u) converges almost surely to h(2) (v) = Â(v, u)h (1) (u)W (

EQUATION

(1)

j )W (1)

converges almost surely to Â(v, u)h

(1) t1 (u)W (1) dP (u) and thus to h(2) (v). A simple induction completes the rest of the proof.

Proof of Proposition 2. Conditioned on v, the expectation of y(v) is

EQUATION

and the variance is 1/t times that of Â(v, u)x(u), i.e.,

EQUATION

)

Instantiating ( 11) and ( 12) with iid samples v 1 , . . . , v s ∼ P and taking variance and expectation in the front, respectively, we obtain

EQUATION

E Var 1 s s i=1 y(v i ) v 1 , . . . , v s = 1 st Â(v, u) 2 x(u) 2 dP (u) dP (v) - 1 st e(v) 2 dP (v).

Then, applying the law of total variance

Var 1 s s i=1 y(v i ) = Var E 1 s s i=1 y(v i ) v 1 , . . . , v s +E Var 1 s s i=1 y(v i ) v 1 , . . . , v s ,

we conclude the proof.

for Cora and Reddit, and 1024 for Pubmed. Dropout rate is set as 0. We use Adam as the optimization method for training. In the test phase, we use the trained parameters and all the graph nodes instead of sampling. For more details please check our codes in a temporary git repository https://github.com/matenure/FastGCN.

Hardware: Running time is compared on a single machine with 4-core 2.5 GHz Intel Core i7, and 16G RAM. 4 , even though fewer labels are used here. The reason is that we adopt the same hyperparameters as in Kipf & Welling (2016a) to reproduce the F1 scores of their work, whereas for Table 4 , a better learning rate is found that boosts the performance on the new split of the data, in which case GCN (original) converges faster. 

D CONVERGENCE:

Strictly speaking, the training algorithms proposed in Section 3 do not precisely follow the existing theory of SGD, because the gradient estimator, though consistent, is biased. In this section, we fill the gap by deriving a convergence result. Similar to the case of standard SGD where the convergence rate depends on the properties of the objective function, here we analyze only a simple case; a comprehensive treatment is out of the scope of the present work. For convenience, we will need a separate system of notations and the same notations appearing in the main text may bear a different meaning here. We abbreviate "with probability one" to "w.p.1" for short.

We use f (x) to denote the objective function and assume that it is differentiable. Differentiability is not a restriction because for the nondifferentiable case, the analysis that follows needs simply change the gradient to the subgradient. The key assumption made on f is that it is l-strictly convex; that is, there exists a positive real number l such that f (x) -f (y) ≥ ∇f (y), x -y + l 2 x -y 2 , (13) for all x and y. We use g to denote the gradient estimator. Specifically, denote by g(x; ξ N ), with ξ N being a random variable, a strongly consistent estimator of ∇f (x); that is, Moreover, we consider the SGD update rule

EQUATION

where ξ (k)

N is an indepedent sample of ξ N for the kth update. The following result states that the update converges on the order of O(1/k).

Theorem 5. Let x * be the (global) minimum of f and assume that ∇f (x) is uniformly bounded by some constant G > 0. If γ k = (lk) -1 , then there exists a sequence B k with

B k ≤ max{ x 1 -x * 2 , G 2 /l 2 } k such that x k -x * 2 → B k w.p.1.

Proof. Expanding x k+1 -x * 2 by using the update rule (14), we obtain

x k+1 -x * 2 = x k -x * 2 -2γ k g k , x k -x * + γ 2 k g k 2 ,

where g k ≡ g(x k ; ξ (k) N ). Because for a given x k , g k converges to ∇f (x k ) w.p.1, we have that conditioned on x k ,

EQUATION

)

On the other hand, applying the strict convexity (13), by first taking x = x k , y = x * and then taking x = x * , y = x k , we obtain

EQUATION

)

Substituting ( 16) to ( 15), we have that conditioned on x k ,

x k+1 -x * 2 → C k w.p.1 for some

EQUATION

Now consider the randomness of x k and apply induction. For the base case k = 2, the theorem clearly holds with B 2 = C 1 . If the theorem holds for k = T , let L = max{ x 1 -x * 2 , G 2 /l 2 }. Then, taking the probabilistic limit of x T on both sides of (17), we have that C T converges w.p.1 to some limit that is less than or equal to (1 -2/T )(L/T ) + G 2 /(l 2 T 2 ) ≤ L/(T + 1). Letting this limit be B T +1 , we complete the induction proof.

:



Figures and Tables:

FastGCN batched training (one epoch) 1: for each batch do 2:For each layer l, sample uniformly t l vertices u each layer l do Compute batch gradient ∇L batch 4:

Figure 3 in the main text compares the per-batch training time for different methods. Here, we list the total training time for reference. It is impacted by the convergence of SGD, whose contributing factors include learning rate, batch size, and sample size. See Table4. Although the orders-ofmagnitude speedup of per-batch time is slightly weakened by the convergence speed, one still sees a substantial advantage of the proposed method in the overall training time. Note that even though the original GCN trains faster than the batched version, it does not scale because of memory limitation. Hence, a fair comparison should be gauged with the batched version. We additionally show in Figure4the evolution of prediction accuracy as training progresses.

ξ N ) = ∇f (x) w.p.1.

Dataset Statistics

Table content: <table><tr><td>Dataset</td><td>Nodes</td><td>Edges</td><td colspan="2">Classes Features</td><td>Training/Validation/Test</td></tr><tr><td>Cora</td><td>2, 708</td><td>5, 429</td><td>7</td><td>1, 433</td><td>1, 208/500/1, 000</td></tr><tr><td colspan="2">Pubmed 19, 717</td><td>44, 338</td><td>3</td><td>500</td><td>18, 217/500/1, 000</td></tr><tr><td>Reddit</td><td colspan="2">232, 965 11, 606, 919</td><td>41</td><td>602</td><td>152, 410/23, 699/55, 334</td></tr><tr><td colspan="6">Implementation details are as following. All networks (including those under comparison) contain</td></tr><tr><td colspan="6">two layers as usual. The codes of GraphSAGE and GCN are downloaded from the accompany</td></tr><tr><td colspan="6">websites and the latter is adapted for FastGCN. Inference with FastGCN is done with the full GCN</td></tr><tr><td colspan="6">network, as mentioned in Section 3.2. Further details are contained in the appendix.</td></tr></table>

Benefit of precomputing ÂH(0) for the input layer. Data set: Pubmed. The three data sets from top to bottom are ordered the same as Table1.the last step is a constant throughout training. Hence, one may precompute the product rather than sampling this layer to gain efficiency. The compared results are listed on the right part of Table2(columns under "Precompute"). One sees that the training time substantially decreases while the accuracy is comparable. Hence, all the experiments that follow use precomputation.

Table content: <table><tr><td colspan="2">Train-ing time is in seconds, per-epoch (batch size 1024). Accuracy is measured by using micro F1 score.</td><td>F1 Score F1 Score</td><td>0.8 0.85 0.75 0.8 0.85 0.9</td><td>10 10</td><td>25 25</td><td>50 50</td><td>Uniform Importance Uniform Importance</td></tr><tr><td>Sampling Time F1 25 0.760 0.873 0.144 0.879 Precompute t 1 Time F1 5 10 0.755 0.863 0.141 0.870 0.737 0.859 0.139 0.849</td><td/><td>F1 Score</td><td>0.9 0.92 0.94</td><td colspan="3">Sample size 25 50 100</td><td>Uniform Importance</td></tr><tr><td>50 0.774 0.864 0.142 0.880</td><td colspan="7">Figure 2: Prediction accuracy: uniform versus impor-</td></tr><tr><td/><td colspan="3">tance sampling.</td><td/><td/><td/></tr></table>

. The bar heights indicate the per-batch training time, in the log scale. One sees that GraphSAGE is a substantial improvement of GCN for large and dense graphs (e.g., Reddit), although for smaller ones (Cora and Pubmed), GCN trains faster. FastGCN is the fastest, with at least an order of magnitude improvement compared with the runner up (except for Cora), and approximately two orders of magnitude speed up compared with the slowest. Here, the training time of FastGCN is with respect to the sample size that achieves the best prediction accuracy. As seen from the table on the right, this accuracy is highly comparable with the best of the other two methods.

Table content: <table><tr><td/><td/><td/><td/><td/><td>FastGCN</td></tr><tr><td/><td colspan="2">10</td><td>0</td><td/><td>GraphSAGE GCN</td><td>Micro F1 Score</td></tr><tr><td>Time (seconds)</td><td colspan="3">10 10 -1 -2</td><td/><td>Cora Pubmed Reddit 0.850 0.880 0.937 GraphSAGE-GCN 0.829 FastGCN 0.849 0.923 GraphSAGE-mean 0.822 0.888 0.946 GCN (batched) 0.851 0.867 0.930</td></tr><tr><td/><td/><td/><td/><td/><td>GCN (original)</td><td>0.865</td><td>0.875</td><td>NA</td></tr><tr><td/><td>10</td><td colspan="2">-3</td><td>Cora</td><td>Pubmed</td><td>Reddit</td></tr></table>

Further comparison of per-batch training time (in seconds) with new implementation of GraphSAGE for small graphs. The new implementation is in PyTorch whereas the rest are in Ten-sorFlow.

Table content: <table><tr><td/><td>Cora</td><td>Pubmed Reddit</td></tr><tr><td>FastGCN</td><td colspan="2">0.0084 0.0047 0.0129</td></tr><tr><td colspan="3">GraphSAGE-GCN (old impl) 1.1630 0.3579 0.4260</td></tr><tr><td colspan="3">GraphSAGE-GCN (new impl) 0.0380 0.3989</td><td>NA</td></tr><tr><td>GCN (batched)</td><td colspan="2">0.0166 0.0815 2.1731</td></tr></table>

Jian Tang, Meng Qu, Mingzhe Wang, Ming Zhang, Jun Yan, and Qiaozhu Mei. Line: Large-scale information network embedding. In Proceedings of the 24th International Conference on World Wide Web, WWW '15, pp. 1067-1077, 2015. ISBN 978-1-4503-3469-3. Daixin Wang, Peng Cui, and Wenwu Zhu. Structural deep network embedding. In Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining,

Table content: <table/>

Total training time (in seconds).

Table content: <table><tr><td/><td/><td/><td/><td/><td colspan="4">Cora Pubmed Reddit</td></tr><tr><td/><td/><td/><td colspan="2">FastGCN</td><td>2.7</td><td>15.5</td><td/><td>638.6</td></tr><tr><td/><td/><td/><td colspan="3">GraphSAGE-GCN 72.4</td><td>259.6</td><td colspan="2">3318.5</td></tr><tr><td/><td/><td/><td colspan="2">GCN (batched)</td><td>6.9</td><td colspan="3">210.8 58346.6</td></tr><tr><td/><td/><td/><td colspan="2">GCN (original)</td><td>1.7</td><td>21.4</td><td/><td>NA</td></tr><tr><td/><td>FastGCN</td><td>GraphSAGE</td><td>GCN (batched)</td><td>FastGCN</td><td>GraphSAGE</td><td colspan="2">GCN (batched)</td><td>FastGCN</td><td>GraphSAGE</td><td>GCN (batched)</td></tr><tr><td/><td>1</td><td/><td/><td>0.9</td><td/><td/><td/><td>0.95</td></tr><tr><td>Training accuracy</td><td>0.2 0.4 0.6 0.8</td><td/><td>Training accuracy</td><td>0.82 0.84 0.86 0.88</td><td/><td/><td>Training accuracy</td><td>0.75 0.8 0.85 0.9</td></tr><tr><td/><td>10 -2 0</td><td>10 0</td><td>10 2</td><td>10 -2 0.8</td><td>10 0</td><td>10 2</td><td>10 4</td><td>10 0 0.7</td><td>10 5</td></tr><tr><td/><td colspan="3">Training time (seconds)</td><td colspan="3">Training time (seconds)</td><td/><td colspan="2">Training time (seconds)</td></tr><tr><td/><td>FastGCN</td><td>GraphSAGE</td><td>GCN (batched)</td><td>FastGCN</td><td>GraphSAGE</td><td colspan="2">GCN (batched)</td><td>FastGCN</td><td>GraphSAGE</td><td>GCN (batched)</td></tr><tr><td/><td>1</td><td/><td/><td>0.9</td><td/><td/><td/><td>0.95</td></tr><tr><td>Test accuracy</td><td>0.4 0.6 0.8</td><td/><td>Test accuracy</td><td>0.84 0.86 0.88</td><td/><td/><td>Test accuracy</td><td>0.8 0.85 0.9</td></tr><tr><td/><td>10 -2 0.2</td><td>10 0</td><td>10 2</td><td>10 -2 0.82</td><td>10 0</td><td>10 2</td><td>10 4</td><td>10 0 0.75</td><td>10 5</td></tr><tr><td/><td colspan="3">Training time (seconds)</td><td colspan="3">Training time (seconds)</td><td/><td colspan="2">Training time (seconds)</td></tr><tr><td colspan="10">Figure 4: Training/test accuracy versus training time. From left to right, the data sets are Cora,</td></tr><tr><td colspan="4">Pubmed, and Reddit, respectively.</td><td/><td/><td/><td/><td/></tr><tr><td colspan="5">C.2 ORIGINAL DATA SPLIT FOR CORA</td><td/><td/><td/><td/></tr></table>

The results for GCN are consistent with those reported byKipf & Welling (2016a). Because labeled data are scarce, the training of GCN is quite fast. FastGCN beats it only on Pubmed. The accuracy results of FastGCN are inferior to GCN, also because of the limited number of training labels. The transductive version FastGCN-transductive matches the accuracy of that of GCN. The results for GraphSAGE are curious. We suspect that the model significantly overfits the data, because perfect training accuracy (i.e., 1) is attained.One may note a subtlety that the training ofGCN (original)  is slower than what is reported in Table

Table content: <table/>

Total training time and test accuracy for Cora and Pubmed, original data split. Time is in seconds.

Table content: <table><tr><td/><td>Cora</td><td/><td colspan="2">Pubmed</td></tr><tr><td/><td>Time</td><td>F1</td><td>Time</td><td>F1</td></tr><tr><td>FastGCN</td><td colspan="2">2.52 0.723</td><td colspan="2">0.97 0.721</td></tr><tr><td>FastGCN-transductive</td><td colspan="2">5.88 0.818</td><td colspan="2">8.97 0.776</td></tr><tr><td>GraphSAGE-GCN</td><td colspan="4">107.95 0.334 39.34 0.386</td></tr><tr><td>GCN (original)</td><td colspan="4">2.18 0.814 32.65 0.795</td></tr></table>

Peer review (iclr2018):
 Update: I have read the rebuttal and the revised manuscript. Additionally I had a brief discussion with the authors regarding some aspects of their probabilistic framework. I think that batch training of GCN is an important problem and authors have proposed an interesting solution to this problem. I appreciated all the work authors put into the revision. In this regard, I have updated my rating. However, I am not satisfied with how the probabilistic problem formulation was presented in the paper. I would appreciate if authors were more upfront about the challenges of the problem they formulated and limitations of their results. I briefly summarize the key missing points below, although I acknowledge that solution to such questions is out of scope of this work. 1. Sampling of graph nodes from P is not iid. Every subsequent node can not be equal to any of the previous nodes. Hence, the distribution changes and subsequent nodes are dependent on previous ones. However, exchangeability could be a reasonable assumption to make as order (in the joint distribution) does not matter for simple choices of P. Example: let V be {1,2,3} and P a uniform distribution. First node can be any of the {1,2,3}, second node given first (suppose first node is *2*) is restricted to {1,3}. There is clearly a dependency and change of distribution. 2. Theorem 1 is proven under the assumption that it is possible to sample from P and utilize Monte Carlo type argument. However, in practice, sampling is done from a uniform distribution over observed samples. Also, authors suggest that V may be infinite. Recall that for Monte Carlo type approaches to work, sampling distribution is ought to contain support of the true distribution. Observed samples (even as sample size goes to infinity) will never be able to cover an infinite V. Hence, Theorem 1 will never be applicable (for the purposes of evaluating population loss). Also note that this is different from a more classical case of continuous distributions, where sampling from a Gaussian, for instance, will cover any domain of true distribution. In the probabilistic framework defined by the authors it is impossible to cover domain of P, unless whole V is observed. ---------------------------------------------------------------------- This work addresses a major shortcoming of recently popularized GCN. That is, when the data is equipped with the graph structure, classic SGD based methods are not straightforward to apply. Hence it is not clear how to deal with large datasets (e.g., Reddit). Proposed approach uses an adjacency based importance sampling distribution to select only a subset of nodes on each GCN layer. Resulting loss estimate is shown to be consistent and its gradient is used to perform the weight updates. Proposed approach is interesting and the direction of the work is important given recent popularity of the GCN. Nonetheless I have two major question and would be happy to revisit my score if at least one is addressed. Theory: SGD requires an unbiased estimate of the gradient to converge to the global optima in the convex loss case. Here, the loss estimate is shown to be consistent, but not guaranteed to be unbiased and nothing is said about the gradient in Algorithm 1. Could you please provide some intuition about the gradient estimate? I might not be familiar with some relevant results, but it appears to me that Algorithm 1 will not converge to the same solution as full data GD would. Practice: Per batch timings in Fig. 3 are not enough to argue that the method is faster as it might have poor convergence properties overall. Could you please show the train/test accuracies against training time for all compared methods? Some other concerns and questions: - It is not quite cleat what P is. You defined it as distribution over vertices of some (potentially infinite) population graph. Later on, sampling from P becomes equivalent to uniform sampling over the observed nodes. I don*t see how you can define P over anything outside of the training nodes (without defining loss on the unobserved data), as then you would be sampling from a distribution with 0 mass on the parts of the support of P, and this would break the Monte Carlo assumptions. - Weights disappeared in the majority of the analysis. Could you please make the representation more consistent. - a(v,u) in Eq. 2 and A(v,u) in Eq. 5 are not defined. Do they both correspond to entries of the (normalized) adjacency?