FastGCN: Fast Learning with Graph Convolutional Networks via Importance Sampling

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:

where $\hat{A}$ is some normalization of the graph adjacency matrix, $H^{(l)}$ contains the embedding (rowwise) of the graph vertices in the $l$th layer, $W^{(l)}$ is a parameter matrix, and $\sigma$ 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.

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.

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$:

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},\ldots,x_{n}$:

In each step of SGD, the gradient is approximated by $\nabla g(W;x_{i})$, an (assumed) unbiased sample of $\nabla 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 $\nabla 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^{\prime}$ with the vertex set $V^{\prime}$ associated with a probability space $(V^{\prime},F,P)$, such that for the given graph $G$, it is an induced subgraph of $G^{\prime}$ and its vertices are iid samples of $V^{\prime}$ according to the probability measure $P$. For the probability space, $V^{\prime}$ serves as the sample space and $F$ may be any event space (e.g., the power set $F=2^{V^{\prime}}$). 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

For the functional generalization, we write

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 $l$th layer. The embedding functions from two consecutive layers are related through convolution, expressed as an integral transform, where the kernel $\hat{A}(v,u)$ corresponds to the $(v,u)$ element of the matrix $\hat{A}$. 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^{(l)}_{1},\ldots,u^{(l)}_{t_{l}}\sim P$ to approximately evaluate the integral transform (2); that is,

with the convention $h^{(0)}_{t_{0}}\equiv h^{(0)}$. Then, the loss $L$ in (3) admits an estimator

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.

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},\ldots,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,\ldots,t_{l}$, $l=0,\ldots,M-1$. Such a procedure is equivalent to uniformly sampling the rows of $H^{(l)}$ for each $l$. Then, we obtain the batch loss

Here, the $n$ inside the activation function $\sigma$ 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.

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

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

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$.

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.

where $b(u)$ is defined in (7), then the variance of $G_{Q}$ admits

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 $\hat{A}(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:

and sample $t$ vertices $u_{1},\ldots,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_{\text{batch}}$ in (4) now is recursively expanded as

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.

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 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 $l$th 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.

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 GraphSAGE, 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.

Implementation details are as following. All networks (including those under comparison) contain two layers as usual. The codes of GraphSAGE and GCN are downloaded from the accompany websites and the latter is adapted for FastGCN. Inference with FastGCN is done with the full GCN network, as mentioned in Section 3.2. Further details are contained 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 $\hat{A}H^{(0)}$ in the bottom layer does not change, which means that the chained expansion of the gradient with respect to $W^{(0)}$ in 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 Table 2 (columns under “Precompute”). One sees that the training time substantially decreases while the accuracy is comparable. Hence, all the experiments that follow use precomputation.

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. 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.

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.

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.

References

Appendix A Proofs

Because the samples $u^{(0)}_{j}$ are iid, by the strong law of large numbers,

Conditioned on $v$, the expectation of $y(v)$ is

and the variance is $1/t$ times that of $\hat{A}(v,u)x(u)$, i.e.,

Instantiating (11) and (12) with iid samples $v_{1},\ldots,v_{s}\sim P$ and taking variance and expectation in the front, respectively, we obtain

Conditioned on $v$, the variance of $y_{Q}(v)$ is $1/t$ times that of

Then, following the proof of Proposition 2, the overall variance is

Hence, the optimal $dQ(u)$ must be proportional to $b(u)|x(u)|\,dP(u)$. Because it also must integrate to unity, we have

Conditioned on $v$, the variance of $y_{Q}(v)$ is $1/t$ times that of

The rest of the proof follows that of Proposition 2. ∎

Appendix B Additional Experiment Details

GCN: The original GCN cannot work on very large graphs (e.g., Reddit). So we modified it into a batched version by simply removing the sampling in our FastGCN (i.e., using all the nodes instead of sampling a few in each batch). For relatively small graphs (Cora and Pubmed), we also compared the results with the original GCN.

GraphSAGE: For training time comparison, we use GraphSAGE-GCN that employs GCN as the aggregator. It is also the fastest version among all choices of the aggregators. For accuracy comparison, we also compared with GraphSAGE-mean. We used the codes from https://github.com/williamleif/GraphSAGE. Following the setting of Hamilton et al. (2017), we use two layers with neighborhood sample sizes $S_{1}=25$ and $S_{2}=10$. For fair comparison with our method, the batch size is set to be the same as FastGCN, and the hidden dimension is 128.

B.2 Experiment Setup

Datasets: The Cora and Pubmed data sets are from https://github.com/tkipf/gcn. As we explained in the paper, we kept the validation index and test index unchanged but changed the training index to use all the remaining nodes in the graph. The Reddit data is from http://snap.stanford.edu/graphsage/.

Experiment Setting: We preformed hyperparameter selection for the learning rate and model dimension. We swept learning rate in the set {0.01, 0.001, 0.0001}. The hidden dimension of FastGCN for Reddit is set as 128, and for the other two data sets, it is 16. The batch size is 256 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.

Appendix C Additional Experiments

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 Table 4. Although the orders-of-magnitude 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 Figure 4 the evolution of prediction accuracy as training progresses.

C.2 Original Data Split for Cora and Pubmed

As explained in Section 4, we increased the number of labels used for training in Cora and Pubmed, to align with the supervised learning setting of Reddit. For reference, here we present results by using the original data split with substantially fewer training labels. We also fork a separate version of FastGCN, called FastGCN-transductive, that uses both training and test data for learning. See Table 5.

The results for GCN are consistent with those reported by Kipf & 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 of GCN (original) is slower than what is reported in Table 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.

Appendix 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

for all $x$ and $y$. We use $g$ to denote the gradient estimator. Specifically, denote by $g(x;\xi_{N})$, with $\xi_{N}$ being a random variable, a strongly consistent estimator of $\nabla f(x)$; that is,

Moreover, we consider the SGD update rule

where $\xi_{N}^{(k)}$ is an indepedent sample of $\xi_{N}$ for the $k$th update. The following result states that the update converges on the order of $O(1/k)$.

Let $x^{*}$ be the (global) minimum of $f$ and assume that $\|\nabla f(x)\|$ is uniformly bounded by some constant $G>0$. If $\gamma_{k}=(lk)^{-1}$, then there exists a sequence $B_{k}$ with

such that $\|x_{k}-x^{*}\|^{2}\to B_{k}$ w.p.$1$.

Expanding $\|x_{k+1}-x^{*}\|^{2}$ by using the update rule (14), we obtain

where $g_{k}\equiv g(x_{k};\xi_{N}^{(k)})$. Because for a given $x_{k}$, $g_{k}$ converges to $\nabla f(x_{k})$ w.p.$1$, we have that conditioned on $x_{k}$,

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

Substituting (16) to (15), we have that conditioned on $x_{k}$,

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})\leq L/(T+1)$. Letting this limit be $B_{T+1}$, we complete the induction proof. ∎