Title: ON THE APPROXIMATION ERRORS OF NODE SAM-PLING FOR GRAPH NEURAL NETWORKS Anonymous authors Paper under double-blind review

Abstract: The recent advancements in graph neural networks (GNNs) have led to state-of-theart performances in various applications, including chemo-informatics, questionanswering systems, and recommender systems. However, scaling up these methods to huge graphs such as social network graphs and web graphs still remains a challenge. In particular, the existing methods for accelerating GNNs are either not theoretically guaranteed in terms of approximation error, or they require at least a linear time computation cost. In this study, we reveal the query complexity of the uniform node sampling scheme for GraphSAGE, the graph attention networks (GAT), and the graph convolutional networks (GCN). The key advantage of our analysis is that the complexity is completely independent of the numbers of the nodes, edges, and neighbors of the input and depends only on the error tolerance and confidence probability while providing theoretical guarantee for the approximation error. To the best of our knowledge, this is the first work to give a theoretical guarantee of approximation for GNNs within a constant time complexity. Through experiments using synthetic and real-world datasets, we demonstrate the speed and precision of the node sampling scheme and validate our theoretical results.

1. INTRODUCTION:

Machine learning on graph structures has various applications such as chemo-informatics (Gilmer et al., 2017) , question answering systems (Schlichtkrull et al., 2018) , and recommender systems (Fan et al., 2019) . Recently, a novel machine learning model for graph data called graph neural networks (GNNs) (Gori et al., 2005; Scarselli et al., 2009; Kipf & Welling, 2017; Hamilton et al., 2017) demonstrated state-of-the-art performances in various graph learning tasks. However, large scale graphs such as social network graphs and web graphs contain billions of nodes, and even a linear time computation cost per iteration is prohibited. Therefore, applying GNNs to huge graphs is challenging. Although Ying et al. (2018) succeeded in applying GNNs to a web-scale network using MapReduce, it still requires massive computational resources.

There are several node sampling techniques to reduce GNN computation. For example, an empirical neighbor sampling scheme is used to speed up GraphSAGE (Hamilton et al., 2017) . FastGCN employs a random layer-wise node sampling (Chen et al., 2018b) . Huang et al. (2018) further improved FastGCN by using an adaptive sampling technique to reduce the variance of estimators. Chen et al. (2018a) proposed a variant of neighbor sampling, which used historical activations to reduce the estimator variance. Overall, the existing sampling techniques for GNNs work well in practice. However, these techniques are either not theoretically guaranteed in terms of approximation error, or thy require at least a linear time computation cost.

In this study, we consider the problem of approximating the embedding of one node using GNNs in constant time with maximum precision1 . We analyze the neighbor sampling technique (Hamilton et al., 2017) to show that a constant number of samples are needed to guarantee the approximation error. It should be noted that the neighbor sampling was introduced as a heuristic method originally, and they did not provide any theoretical guarantees. Specifically, given an error tolerance ε and Table 1 : indicates neighbor sampling approximates the network in constant time. indicates any algorithm cannot approximate the network in constant time. in the Gradient column indicates the error between the gradient of the approximated embedding and that of the exact embedding is also theoretically bounded. * needs an additional condition to approximate it in constant time. confidence probability 1 -δ, our analysis shows that the estimate ẑv of the exact embedding z v of a node v such that Pr[ ẑvz v 2 ≥ ε] ≤ δ and the estimate ∂zv ∂θ of the exact gradient ∂zv ∂θ of the embedding z v with respect to the network parameters θ, such that Pr[ ∂zv ∂θ -∂zv ∂θ F ≥ ε] ≤ δ can be computed in a constant time. Especially, the uniform node sampling can approximate the exact embedding and its gradients within O( 1 ε 2L (log 1 ε + log 1 δ ) L-1 log 1 δ ) time, where L denotes the number of layers. This complexity is completely independent of the number of nodes, edges, and neighbors of the input, which enables us to deal with graphs irrespective of their size. Moreover, the complexity is a polynomial with respect to 1 ε and log 1 δ . We demonstrate that the time complexity is optimal when L = 1 with respect to the error tolerance ε.

Through experiments, we show that the approximation error between the exact computation and its approximation rapidly converges to zero. To the best of our knowledge, this is the first constant time approximation algorithm for GNNs with a theoretical guarantee in terms of approximation error.

Contributions::

The contributions of this paper are summarized as follows:

• We analyze the neighbor sampling technique for GraphSAGE, GAT, and GCN to provide theoretical justification. Especially, our analysis shows that the complexity is completely independent of the number of nodes, edges, and neighbors of the input. • We show that some existing GNNs, including the original GraphSAGE (Hamilton et al., 2017) , cannot be approximated in constant time by any algorithm (see Table 1 for details). • We empirically validate our theorems using synthetic and real-world datasets.

2. RELATED WORK:



2.1. GRAPH NEURAL NETWORKS:

Graph neural networks (GNNs) were first introduced by Gori et al. (2005) and Scarselli et al. (2009) . They obtained node embedding by recursively applying the propagation function until convergence. Kipf & Welling (2017) proposed graph convolutional networks (GCN), which significantly outperformed the existing methods, including non-neural network based approaches. Gilmer et al. (2017) proposed the message passing neural networks (MPNNs), a general framework of GNNs using the message passing mechanism. Veličković et al. (2018) proposed the graph attention networks (GAT), which incorporate the attention mechanism into GNNs. With the advent of GAT, various GNN models with the attention mechanism have been recently proposed (Wang et al., 2019; Park et al., 2019) .

GraphSAGE (Hamilton et al., 2017 ) is another GNN model, which employs neighbor sampling to reduce the computational costs of training and inference. Owing to neighbor sampling, GraphSAGE can deal with large graphs. However, neighbor sampling was introduced without any theoretical guarantee, and the number of samples is chosen empirically. An alternative computationally efficient GNN would be FastGCN (Chen et al., 2018b) , which employs layer-wise random node sampling to speed up training and inference. Huang et al. (2018) further improved FastGCN by using an adaptive node sampling technique to reduce the variance of estimators. Thanks to the adaptive sampling technique, it reduces the computational costs and outperforms neighbor sampling in terms of classification accuracy and convergence speed. Chen et al. (2018a) proposed an alternative neighbor sampling technique, which uses historical activations to reduce the estimator variance. Additionally, it could achieve zero variance after a certain number of iterations. However, because it used the same sampling technique of GraphSAGE to obtain the initial solution, the approximation error was not theoretically bounded until the Ω(n)-th iteration. Overall, the existing sampling techniques work well in practice. However, these techniques are either not theoretically guaranteed in terms of approximation error, or they require at least a linear time computation cost to calculate the embedding of a node and its gradient of GNN models. Moreover, it is not clear whether we can apply the sampling techniques to state-of-the-art GNN models such as GAT.

2.2. SUBLINEAR TIME ALGORITHMS:

The sublinear time algorithms were originally proposed for property testing (Rubinfeld & Sudan, 1996) . Sublinear property testing algorithms check whether the input has some property π or the input is sufficiently far from the property π with high probability in sublinear time with respect to the input size. Sublinear time approximation algorithms are another type of sublinear time algorithms. More specifically, they calculate a value sufficiently close to the exact value with high probability in sublinear time. Constant time algorithms are a subclass of sublinear time algorithms. They work not only in sublinear time with respect to the input size but also in constant time. The proposed algorithm is classified as a constant time approximation algorithm.

The examples of sublinear time approximation algorithms include minimum spanning tree in metric space (Czumaj & Sohler, 2004) and minimum spanning tree with integer weights (Chazelle et al., 2005) . Parnas & Ron (2007) proposed a method to convert distributed local algorithms into constant time approximation algorithms. In their study, they proposed a method to construct constant time algorithms for the minimum vertex cover problem and dominating set problem. A classic example of sublinear time algorithms related to machine learning includes clustering (Indyk, 1999; Mishra et al., 2001) . Examples of recent work in this stream include constant time approximation of the minimum value of quadratic functions (Hayashi & Yoshida, 2016) and constant time approximation of the residual error of the Tucker decomposition (Hayashi & Yoshida, 2017) . They adopted simple sampling strategies to obtain theoretical guarantee similar to our work.

In this paper, we provide theoretical guarantee for approximation of GNNs within a constant time for the first time.

3. BACKGROUND:



3.1. NOTATIONS:

Let G be the input graph, V = {1, 2, . . . , n} be the set of nodes, n = |V| be the number of nodes, E be the set of edges, m = |E| be the number of edges, deg(v) be the degree of a node v, N (v) be the set of neighbors of a node v, x v ∈ R d0 be the feature vector associated to a node v ∈ V, X = (x 1 , x 2 . . . , x n ) ∈ R n×d0 be the stacked feature vectors, and denotes the matrix transpose.

3.2. NODE EMBEDDING MODEL:

Algorithm 1 O z : Exact embedding

Require: Graph G = (V, E); Features X ∈ R n×d0 ; Node index v ∈ V; Model parameters θ. Ensure: Exact embedding z v 1: z (0) i ← x i (∀i ∈ V) 2: for l ∈ {1, . . . , L} do 3: for i ∈ V do 4: h (l) i ← u∈N (i) M liu (z (l-1) i , z (l-1) u , e iu , θ) 5: z (l) i ← U l (z (l-1) i , h (l) i , θ) 6: end for 7: end for 8: return z (L) v

We consider the node embedding problem using GNNs. Especially, we employ the message passing neural networks (MPNNs) framework (Gilmer et al., 2017) . This framework includes many GNN models, such as GraphSAGE and GCN. Algorithm 1 shows the algorithm of MPNNs. We refer to z (L) v as z v simply. The aim of this study is to develop a constant time approximation algorithm for calculating the embedding vector z v and gradients ∂zv ∂θ with the given model parameters θ and node v.

3.3. COMPUTATIONAL MODEL ASSUMPTIONS:

We have to specify how to access the input to design constant time algorithms because the constant time algorithms cannot read the entire input. We follow the standard convention of sublinear time algorithms (Parnas & Ron, 2007; Nguyen & Onak, 2008) . We model our algorithm as an oracle machine that can generate queries regarding the input, and we measure the complexity by query complexity. Algorithms can access the input only by querying the following oracles:  and (3) O feature (v) : the feature of node v. We assume that our algorithm can query the oracles in constant time per query.

(1) O deg (v): the degree of node v, (2) O G (v, i): the i-th neighbor of node v,

3.4. PROBLEM FORMULATION:

Given a node v, we calculate the following functions with the least number of oracle accesses:

(1) O z (v): the embedding z v and (2) O g (v): the gradients of parameters ∂zv ∂θ . However, the exact computation of O z and O g needs at least deg(v) queries to aggregate the features from the neighbor nodes, which is not constant time with respect to the input size. Thus, it is computationally expensive to execute the algorithm for a huge network. Therefore, we consider making the following approximations:

• Ôz (v, ε, δ): an estimate ẑv of z v such that Pr[ ẑi -z i 2 ≥ ε] ≤ δ, • Ôg (v, ε, δ): an estimate ∂zv ∂θ of ∂zv ∂θ such that Pr[ ∂zv ∂θ -∂zv ∂θ F ≥ ε] ≤ δ,

where ε > 0 is the error tolerance, 1 -δ is the confidence probability, and • 2 and • F are the Euclidean and Frobenius norm, respectively.

Under the fixed model structure (i.e., the number of layers L, the message passing functions, and the update functions), we construct an algorithm that calculates Ôz and Ôg in constant time irrespective of the number of nodes, edges, and neighbors of the input.

However, it is impossible to construct a constant time algorithm without any assumption about the inputs. Therefore, we make the following mild assumptions:

Assumption 1 ∃B ∈ R s.t. x i 2 ≤ B, e iu 2 ≤ B, and θ 2 ≤ B.

Assumption 2 deg(i)M liu and U l are uniformly continuous in any bounded domain.

Assumption 3 (Only for gradient computation) deg(i)DM liu and DU l are uniformly continuous in any bounded domain, where D denotes the Jacobian operator.

Building a constant time algorithm is impossible without these assumptions as shown in Section 5. Furthermore, we derive the query complexity of our algorithm when the message functions and update functions satisfy the following assumptions. 

4. PROPOSED METHOD:



4.1. CONSTANT TIME EMBEDDING APPROXIMATION:

Here, we build a constant time approximation algorithm based on neighbor sampling, which approximates the embedding z v with an absolute error of at most ε and probability 1 -δ. We recursively construct the algorithm layer by layer by sampling r (l) neighboring nodes in layer l. We refer to the algorithm that calculates the estimate of embeddings in the l-th layer z (l) as Ô(l) z (l = 1, . . . , L). The pseudo code is presented in Algorithm 2. Here, Ô(l-1) z (v ← u) means calling the function Ô(l-1) z with the same parameters as the current function except for v, which is replaced by u. In the following, we demonstrate the theoretical properties of Algorithm 2.

Algorithm 2 Ô(l) z : Estimate the embedding z (l) v

Require: Graph G = (V, E) (as oracle); Features X ∈ R n×d0 (as oracle); Node index v ∈ V; Model parameters θ; Error tolerance ε; Confidence probability 1 -δ. Ensure: Approximation of the embedding z (l) v 1: S (l) ← sample r (l) (ε, δ) neighbors of v with uniform random with replacement.

2:

ĥ(l) v ← Odeg(v) r (l) u∈S (l) M lvu (O feature (v), O feature (u), e iu , θ) (l = 1) Odeg(v) r (l) u∈S (l) M lvu ( Ô(l-1) z (v ← v), Ô(l-1) z (v ← u), e iu , θ) (l > 1) 3: ẑ(l) v ← U l ( Ô(l-1) z (v ← v), ĥ(l) v , θ) if l > 1 otherwise U l (O feature (v), ĥ(l) v , θ) 4: return ẑi

Theorem 1. For all ε > 0, 1 > δ > 0, there exists r (l) (ε, δ) (l = 1, . . . , L) such that for all inputs satisfying Assumptions 1 and 2, the following property holds true:

Pr[ Ôz (v, ε, δ) -z v 2 ≥ ε] ≤ δ.

Theorem 1 shows that the approximation error of Algorithms 2 is bounded by ε with the probability 1 -δ. It is proved by Hoeffding's inequality (Hoeffding, 1963) . The complete proof is available in the appendices. Because the number of sampled nodes depends only on ε and δ and independent of the number of nodes, edges, and neighbors of the input, Algorithm 2 works in constant time. We provide the complexity when the functions are Lipschitz continuous.

Theorem 2. Under Assumption 1 and 4, r

(L) = O( 1 ε 2 log 1 δ ) and r (1) , . . . , r (L-1) = O( 1 ε 2 (log 1 ε + log 1 δ )) are sufficient, and the query complexity of Algorithms 2 is O( 1 ε 2L (log 1 ε + log 1 δ ) L-1 log 1 δ ).

We show that the query complexity of Algorithm 2 is optimal with respect to ε if the number of layers is one. In other words, a one-layer model cannot be approximated in o( 1 ε 2 ) time by any algorithm. Theorem 3. Under Assumptions 1 and 4 and L = 1, the time complexity of Algorithm 2 in Theorem 2 is optimal with respect to the error tolerance ε.

The proof is based on Chazelle et al.'s lemma (Chazelle et al., 2005) . The optimality when L ≥ 2 is an open problem.

4.2. CONSTANT TIME GRADIENT APPROXIMATION:

Next, we propose a constant time algorithm that approximates the gradient of embeddings with respect to the model parameters with an absolute error of at most ε and probability 1 -δ. The basic strategy is to execute Algorithm 2 and calculate the gradients of the embedding z v . Let ∂zv ∂θ be the gradient of the embedding z v with respect to the model parameter θ (i.e., ( ∂zv ∂θ ) ijk = ∂zvi ∂θ jk ). Theorem 4. For all ε > 0, 1 > δ > 0, there exists r (l) (ε, δ) (l = 1, . . . , L) such that for all inputs satisfying Assumptions 1, 2, and 3, the following property holds true:

Pr[ ∂z (L) v ∂θ - ∂z (L) v ∂θ F ≥ ε] ≤ δ, where ∂z (L) v ∂θ is the gradient of ẑ(L) v , which is obtained by Ô(L) z (v, ε, δ), with respect to θ.

Therefore, we can calculate an estimate of the gradient of the embedding with respect to parameters with an absolute error of at most ε and probability 1 -δ by running Ô(L) z (v, ε, δ) and calculating the gradient of the obtained estimate of the embedding. We provide the complexity when the functions are Lipschitz continuous.

Theorem 5. Under Assumptions 1, 4, and 5, r

(L) = O( 1 ε 2 log 1 δ ) and r (1) , . . . , r (L-1) = O( 1 ε 2 (log 1 ε + log 1 δ ))

are sufficient, and the gradient of the embedding with respect to parameters can be approximated with an absolute error of at most ε and probability

1 -δ in O( 1 ε 2L (log 1 ε + log 1 δ ) L-1 log 1 δ ) time.

Note that, technically, MPNNs do not include GAT because GAT uses embeddings of other neighboring nodes to calculate the attention value. However, Theorem 1 and 4 can be naturally extended to GAT, and we can approximate GAT in constant time by neighbor sampling. The details are described in the appendices.

5. INAPPROXIMABILITY:

In this section, we show that some existing GNNs cannot be approximated in constant time. These theorems state that these models cannot be approximated in constant time not only by our algorithm but also by any other algorithm. In other words, for any algorithm that works in constant time, there exists an error tolerance ε, a confidence probability 1 -δ, and a counter example input such that the approximation error for the input is more than ε with probability δ. It indicates that the application of an approximation method to these models requires great supervision because the obtained embedding may be significantly different from the exact embedding. Theorem 6. If x i 2 or θ F is not bounded, even under Assumption 1 and 2, the embeddings of GraphSAGE-GCN cannot be approximated with arbitrary precision and probability in constant time.

Theorem 7. Even under Assumption 1, the embeddings and gradients of GraphSAGE-GCN with ReLU activation and normalization cannot be approximated with arbitrary precision and probability in constant time.

We confirm Theorem 7 through computational experiments in Section 6. Theorem 8. Even under Assumptions 1 and 2, the gradients of GraphSAGE-GCN with ReLU activation cannot be approximated with arbitrary precision and probability in constant time.

However, it should be noted that the embeddings of GraphSAGE-GCN with ReLU activation (without normalization layer) can be approximated in constant time using our algorithm by Theorem 1. The following two theorems state that these models cannot be approximated in constant time even under Assumptions 1, 2, and 3. Theorem 9. Even under assumptions 1, 2, and 3, the embeddings of GraphSAGE-pool cannot be approximated with arbitrary precision and probability in constant time.

Theorem 10. Even under Assumptions 1, 2, and 3, the embeddings of GCN cannot be approximated with arbitrary precision and probability in constant time.

Constant Time Approximaiton for GCN: Due to the inapproximability theorems, we cannot approximate GCN in constant time. However, we can approximate GCN in constant time if the input graph satisfies the following property:

Assumption 6 There exists a constant C ∈ R such that for any input graph G = (V, E) and node v, u ∈ V, the ratio of deg(v

) to deg(u) is at most C (i.e., deg(v)/deg(u) ≤ C).

Assumption 6 prohibits input graphs that have a skewed degree distribution. GCN needs Assumption 6 because (1) the norm of the embedding is not bounded and (2) the influence of anomaly nodes with low degrees is significant without Assumption 6. It should be noted that the GraphSAGE-pool cannot be approximated in constant time even under Assumption 6.

6. EXPERIMENTS:

We will answer the following questions through experiments:

Q1: How fast is the constant time approximation algorithm (Algorithm 2)? Q2: Does Algorithm 2 accurately approximate the embeddings of GraphSAGE-GCN without normalization (Theorem 1), whereas it cannot approximate the original one (Theorem 7)? Q3: Does Algorithm 2 accurately approximate the gradients of GraphSAGE-GCN with sigmoid (Theorem 4), whereas it cannot approximate that with ReLU activation (Theorem 8)? Q4: Is the theoretical rate of the approximation error of Algorithm 2 tight? It should be noted that we focus on showing the approximation error of neighbor sampling in this study as Hamilton et al. (2017) have already reported the effect of neighbor sampling for downstream machine learning tasks (e.g., classification). The experimental details are described in the appendices.

Experiments for Q1: We measure the speed of exact computation and the constant time approximation of the two-layer GraphSAGE-GCN and two-layer GAT. We initialize parameters using the i.i.d. standard multivariate normal distribution. The input graph is a clique K n . We use ten-dimensional vectors from the i.i.d. standard multivariate normal distribution as the node features. We take r (1) = r (2) = 100 samples. Figure 1 (a) shows the speed of these methods as the number of nodes increases. This shows that the constant time approximation is several orders of magnitude faster than the exact computation when the number of nodes is large.

Experiments for Q2: We use the original one-layer GraphSAGE-GCN (with ReLU activation and normalization) and one-layer GraphSAGE-GCN with ReLU activation. The input graph is a clique K n , whose features are x 1 = (1, 0) and x i = (0, 1/n) (i = 1), and the weight matrix is an identity matrix I 2 . We use r (1) = 5, 30, and 100 as the sample size. If a model can be approximated in constant time, the approximation error goes to zero as the sample size increases even if the graph size reaches infinity. The approximation errors of both models are illustrated in Figure 1 (b) .

The approximation error of the original GraphSAGE-GCN converges to approximately 0.75 even if the sample size increases. In contrast, the approximation error without normalization becomes increasingly bounded as the sample size increases. This is consistent with Theorems 1 and 7.

Experiments for Q3: We study the approximation errors of the gradients. We use the one-layer GraphSAGE-GCN with ReLU activation and sigmoid activation. The input graph is a clique K n , whose features are x 1 = (1, 2) and x i = (1, 1) (i = 1), and the weight matrix is ((-1, 1)). We use r (1) = 5, 30, and 100 as the sample size. The approximation error of both models are illustrated in Figure 1 (c ). The approximation error with ReLU activation converges to approximately 1.0 even if the sample size increases. In contrast, the approximation error with sigmoid activation becomes increasingly bounded as the sample size increases. This is consistent with Theorems 4 and 8.

Experiments for Q4: We use the one-layer GraphSAGE-GCN with sigmoid activation. The input graph is a clique K n , where the number of nodes is n = 40000. We set the dimensions of intermediate embeddings as 2, and each feature value is set to 1 with probability 0.5 and -1 otherwise. This satisfies Assumption 1 (i.e., x i 2 ≤ √ 2). We compute the approximation errors of Algorithm 2 with different numbers of samples. Figure 1 (d) illustrates the 99-th percentile point of empirical approximation errors and the theoretical bound by Theorem 2 (i.e., ε = O(r -1/2 )). It shows that the approximation error decreases along with the theoretical rate. It indicates the theoretical rate is tight.

Experiments for Q5: We analyze the instances when neighbor sampling succeeds and fails for a variety of models. First, we use the Barabasi-Albert model (BA model) (Barabasi & Albert, 1999) to generate input graphs. The degree distribution of the BA model follows a power-law. Therefore, neighbor sampling will fail to approximate GCN (Theorem 10, Assumption 6). We use two-layer GraphSAGE-GCN, GAT, and GCN with ReLU activation. We use ten-dimensional vectors from the i.i.d. standard multivariate normal distribution as the node features. In this experiment, we use the same number r of samples in the first and second layer (i.e., r = r (1) = r (2) ) and we change r from 8 to 1000. We use graphs with n = r 2 nodes. Note that if the constant time approximation is possible, the error is bounded as the number of samples increases even if the number of nodes increases. Figure 1 (e) shows the approximation error. It shows that the error of GCN linearly increases even if the number of samples increases. However, the errors of GraphSAGE-GCN and GAT gradually decrease as the number of samples increases. It indicates that we cannot bound the approximation error of GCN however large number of examples we use. This is consistent with Theorem 10. This result indicates that approximating GCN requires great supervision when the input graph is a social network because the degree distribution of a social network presents the power-law as the BA-model.

Next, we use the Erdős-Rényi model (ER model) (Erdős & Rényi, 1959) . It generates graphs with flat degree distribution. We use the two-layer GraphSAGE-GCN, GAT, GCN, and GraphSAGE-pool. Figure 1 (f) shows the approximation error. It shows that the errors of GraphSAGE-GCN, GAT, and GCN gradually decrease as the number of samples increases. This is consistent with Theorem 1 and Assumption 6. In contrast, the approximation error of GraphSAGE-pool does not decrease even if the input graphs are generated by the ER model. This is consistent with Theorem 9. Experiments for Q6: We use three real-world datasets: Cora, PubMed, and Reddit. They contain 2708, 19717, and 232965 nodes, respectively. We randomly choose 500 nodes for validation and 1000 nodes for test and use the remaining nodes for training. We use the two-layer GraphSAGE-GCN with sigmoid activation in this experiment. The dimensions of the hidden layers are set to 128, and we use an additional fully connected layer to predict the labels of the nodes from the embeddings. We train the models with Adam (Kingma & Ba, 2015) with a learning rate of 0.001. We first train ten models with training nodes for each dataset. The micro-F1 scores of Cora, PubMed, and Reddit are 0.877, 0.839, and 0.901, respectively. It should be noted that we do not aim to obtain high classification accuracy here but intend to sanity check the models. In this experiment, we use the same number r of samples in the first and second layer (i.e., r = r (1) = r (2) ) and we change r from 1 to 100. We also calculate the gradient of parameter matrices W (1) and W (2) with respect to the embedding obtained by Algorithm 2. We compare these values with exact embeddings and exact gradients. Figure 6 illustrates the 99-th percentile point of the empirical approximation errors. We normalize them to ensure that the value is 1.0 at r = 1 to demonstrate the decreasing rate of each error. It shows that the approximation errors of embeddings and gradients rapidly decrease for the real-world data.

) Cora PubMed Reddit Cora W (1) Cora W (2) PubMed W (1) PubMed W (2) Reddit W (1) Reddit W (2)

7. CONCLUSION:

We proposed a constant time approximation algorithm for the embedding and gradient computation of GNNs, where the complexity is completely independent of the number of nodes, edges, and neighbors of the input. We proved its theoretical guarantee in terms of the approximation error. This is the first constant time approximation algorithm for GNNs in the literature. We further demonstrated that some existing GNNs cannot be approximated in constant time by any algorithm. Lastly, we validate the theory through experiments using synthetic and real-world datasets.

A MODELS:

We introduce GraphSAGE-GCN, GraphSAGE-mean, GraphSAGE-pool, the graph convolutional networks (GCN), and the graph attention networks (GAT) for completeness of our paper. (Hamilton et al., 2017) : The message function and the update function of this model is

GraphSAGE-GCN:

M lvu (z v , z u , e vu , θ) = z u deg(v) , U l (z v , h v , θ) = σ(W (l) h v ),

where W (l) is a parameter matrix and σ is an activation function such as sigmoid and ReLU. GraphSAGE-GCN includes the center node itself in the set of adjacent nodes (i.e., N (v) ← N (v) ∪ {v}).

GraphSAGE-mean (Hamilton et al., 2017) : The message function and the update function of this model is

M lvu (z v , z u , e vu , θ) = z v deg(v) , U l (z v , h v , θ) = σ(W (l) [z v , h v ]),

where [•] denotes concatenation GraphSAGE-pool (Hamilton et al., 2017) : We do not formulate GraphSAGE-pool using the message function and the update function because it takes maximum instead of summation. The model of GraphSAGE-pool is

z (l) v = max({σ(W (l) z (l-1) u + b) | u ∈ N (v)}).

GCN (Kipf & Welling, 2017) : The message function and the update function of this model is

M lvu (z v , z u , e vu , θ) = z (l) deg(v)deg(u) , U l (z v , h v , θ) = σ(W (l) h v ).

GAT (Veličković et al., 2018) : The message function and the update function of this model is

α (l) vu = exp(LEAKYRELU(a (l) [W (l) z v , W (l) z u ])) u ∈N (v) exp(LEAKYRELU(a (l) [W (l) z v , W (l) z u ])) , M lvu (z v , z u , e vu , θ) = α (l) vu z u , U l (z v , h v , θ) = σ(W (l) h v ).

Note that, technically, MPNNs do not include GAT because GAT uses embeddings of other neighboring nodes to calculate the attention value α vu . However, we can apply the same argument as MPNNs to GAT, and neighbor sampling can approximate GAT in constant time as other MPNNs. To be precise, the following proposition holds true. Proposition 11. If Assumptions 1 holds true and σ is Lipschitz continuous, if we take r (L) = O( 1 ε 2 log 1 δ ) and r (1) , . . . , r (L-1) = O( 1 ε 2 (log 1 ε + log 1 δ )) samples, and let ẑ(0

) v = z (0)

v Lemma 12 states that the empirical mean of n = O( 1 ε 2 log 1 δ ) samples independently sampled from the same distribution is the approximation of the exact mean with an absolute error of at most ε and probability 1 -δ. Lemma 13 (Hoeffding's inequality (Hoeffding, 1963) ). Let X 1 , X 2 , . . . , X n be independent random variables bounded by the intervals [-B, B] and let X be the empirical mean of these variables X = 1 n n i=1 X i . Then, for any ε > 0,

Pr[| X -E[ X]| ≥ ε] ≤ 2 exp - nε 2 2B 2 holds true.

Proof of Lemma 12. Apply Lemma 13 to each dimension k of X i . Then,

Pr[| Xk -E[ X] k | ≥ ε √ d ] ≤ 2 exp - nε 2 2B 2 d . It should be noted that |X ik | < B because X i 2 < B. Therefore, Pr[∃k ∈ {1, 2, . . . , d} | Xk -E[ X] k | ≥ ε √ d ] ≤ 2d exp - nε 2 2B 2 d . If | Xk -E[ X] k | < ε √ d holds true for all dimension k, then X -E[ X] 2 = d k=1 ( Xk -E[ X] k ) 2 < d • ε 2 d = ε. Therefore, Pr[ X -E[ X] 2 ≥ ε] ≤ 2d exp - nε 2 2B 2 d .

First, we prove that the embedding in each layer is bounded.

Lemma 14. Under Assumptions 1 and 2, the norms of the embeddings z

EQUATION

v 2 , and ĥ(l) v 2 (l = 1, . . . , L) are bounded by a constant B ∈ R.

Proof of Lemma 14. We prove the theorem by performing mathematical induction. The norm of the input to the first layer is bounded by Assumption 1. The message function deg(i)M liu and the update function U l is continuous by Assumption 2. Since the image f (X) of a compact set X ∈ R d is compact if f is continuous, the images of deg(i)M liu and U l are bounded by induction.

Proof of Theorem 1. We prove the theorem by performing mathematical induction on the number of layers L.

Base case: It is shown that the statement holds true for L = 1.

EQUATION

x k , e vx k , θ). Then,

EQUATION

There exists a constant C ∈ R such that for any input satisfying Assumption 1,

EQUATION

x k 2 , e vx k 2 , and θ 2 are bounded by Assumption 1 and deg(v)M Lvu is continuous. Therefore, if we take r

(L) = O( 1 ε 2 log 1 δ ) samples, Pr[ ĥ(L) v -h (L) v 2 ≥ ε ] ≤ δ

by the Hoeffding's inequality and equations ( 2) and (3). Therefore, Pr[

ẑ(L) v -z (L) v 2 ≥ ε] ≤ δ.

When the number of layers is one,

EQUATION

If we take r

EQUATION

holds true for any input by the argument of the proof of Theorem 1.

Let x k be the k-th sample in S (L) and

EQUATION

u , e vu , θ). Then,

EQUATION

There exists a constant C ∈ R such that for any input satisfying Assumption 1,

EQUATION

because z

EQUATION

x k 2 , e vx k 2 , and θ 2 are bounded by Assumption 1 and deg(v)DM Lvu is continuous. Therefore, if we take r

(L) = O( 1 ε 2 log 1 δ ) samples, Pr[ ∂ ĥ(L) v ∂θ - ∂h (L) v ∂θ F ≥ ε/(4B )] ≤ δ/2 (19)

holds true by the Hoeffding's inequality and equations ( 17) and ( 18). If ∂ ĥ(L) 

(z (0) v , h (L) v , θ) F ≤ ∂ ĥ(L) v ∂θ ∂U L ∂ ĥ(L) v (z (0) v , ĥ(L) v , θ) - ∂ ĥ(L) v ∂θ ∂U L ∂h (L) v (z (0) v , h (L) v , θ) F + ∂ ĥ(L) v ∂θ ∂U L ∂h (L) v (z (0) v , h (L) v , θ) - ∂h (L) v ∂θ ∂U L ∂h (L) v (z (0) v , h (L) v , θ) F = ∂ ĥ(L) v ∂θ F ∂U L ∂ ĥ(L) v (z (0) v , ĥ(L) v , θ) - ∂U L ∂h (L) v (z (0) v , h (L) v , θ) F + ∂ ĥ(L) v ∂θ - ∂h (L) v ∂θ F ∂U L ∂h (L) v (z (0) v , h (L) v , θ) F ≤ B ∂U L ∂ ĥ(L) v (z (0) v , ĥ(L) v , θ) - ∂U L ∂h (L) v (z (0) v , h (L) v , θ) F + B ∂ ĥ(L) v ∂θ - ∂h (L) v ∂θ F < B ε 4B + B ε 4B = ε 2 .

Therefore, Pr[ 15), ( 16), and (19).

∂z (L) v ∂θ - ∂z (L) v ∂θ F ≥ ε] ≤ δ by equations (

Inductive step: It is shown that the statement holds true for L = l + 1 if it holds true for L = l. ∂z 

(L) v ∂θ = ∂U L ∂θ (z (L-1) v , h (L) v , θ) + ∂U L ∂z (L-1) v (z (L-1) v , h (L) v , θ) ∂z (L-1) v ∂θ + ∂U L ∂h (L) v (z (L-1) v , h (L) v , θ) u∈N (v)

:

e.g., predicting whether a user of an SNS clicks an advertisement by GNNs in real-time (i.e., when the user accesses). A user may have many neighbors, but GNNs must respond in limited time, which prohibits exact computation. It motivates us to approximate the exact computation in limited time with theoretical guarantee.

Figures and Tables:

Figure 1: (a) The inference speed of each method. (b) The approximation error of the original GraphSAGE-GCN (i.e., with ReLU activation and normalization) and GraphSAGE-GCN with ReLU activation. (c) The approximation error of the gradient with ReLU and sigmoid activations. (d) The approximation error of Algorithm 2 and its theoretical bound. (e) The approximation error of GraphSAGE-GCN, GAT, and GCN with the Barabasi-Albert model. (f) The approximation error of GraphSAGE-GCN, GAT, GCN, and GraphSAGE-pool with the Erdős-Rényi model. Q5: Does neighbor sampling fail to approximate GCN when the degree distribution is skewed (Theorem 10)? Does it succeed when the node distribution is flat (Assumption 6)?

Figure 2: Real data.

v , h (L) v , θ) -∂U L ∂ ĥ(L) v (z (0) v , ĥ(L) v , θ) F < ε/(4B ))

Because DU L is uniform continuous, ∃ε > 0 s.t. for all input, [z (L-1) v , h (L) v ] -[ ẑ(L-1) v , ĥ(L) v ] 2 < ε ⇒ ( ∂U L ∂θ (z (L-1) v , h (L) v , θ) -∂U L ∂θ ( ẑ(L-1) v , ĥ(L) v , θ) F < O(ε) ∧ ∂U L ∂h (L) v (z (L-1) v , h (L) v , θ)v , θ) F < O(ε) ∧

Peer review (iclr2020):
In this paper, the authors propose a constant-time approximation for graph convolution operation via theoretical analysis on the number of sampling from each neighbor. The authors prove that both node embedding and gradient can be approximated via constant number of samples among the neighbors. Extensive experiments are carried out to verify the correctness of the proposed bounds. Strength: 1. The authors establish rigorous bounds on the number of required neighbors to sample to guarantee good approximation of the graph convolution operations. The bounds are established on both nodes embedding and gradients. 2. The authors corroborate the theoretical results with extensive experiments. The experiments are carried under the case both when the assumption holds and not. Weakness: 1. Though the bounds do not depend on number of nodes, it does have an exponential dependence on the number of layers. For common GCN with L=2, even epsilon = 0.1 requires more than 10^4 samples. As a result, the bound is less practical in providing guidance for real-world application. Moreover, the bound does not depend on the embedding dimension at all. Is this to be expected? 2. The authors only provide constant bound for convolution on individual nodes. For the final metrics like classification accuracy, we need to guarantee good approximation for all node’s convolution. In this case, the bounds will definitely depend on number of nodes at least through the union bounds among all nodes. 3. The proposed method is exactly the uniform neighbor sampling methods. There is little contribution on the methodology side from the paper. Also, it would be interesting if the authors could provide theoretical analysis for adaptive sampling methods as well. 4. The authors provide in-approximability results under certain conditionals. However, uniform sampling methods works well empirically under this case as well. It would be better if the authors could explain the gap between theory and empirical performance. Detailed Comments: 1. Algorithm 2, please define notation like O(v<-v).