Year: 2019

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.

heading: 1 INTRODUCTION
text: 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
1e.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.
confidence probability 1 − δ, our analysis shows that the estimate ẑv of the exact embedding zv of a node v such that Pr[‖ẑv − zv‖2 ≥ ε] ≤ δ and the estimate ∂̂zv∂θ of the exact gradient ∂zv ∂θ of
the embedding zv 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.

heading: 2 RELATED WORK

heading: 2.1 GRAPH NEURAL NETWORKS
text: 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. Veličković 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.

heading: 2.2 SUBLINEAR TIME ALGORITHMS
text: 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.

heading: 3 BACKGROUND

heading: 3.1 NOTATIONS
text: 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, xv ∈ Rd0 be the feature vector associated to a node v ∈ V , X = (x1,x2 . . . ,xn) > ∈ Rn×d0 be the stacked feature vectors, and > denotes the matrix transpose.
3.2 NODE EMBEDDING MODEL Algorithm 1 Oz: Exact embedding
Require: Graph G = (V, E); Features X ∈ Rn×d0 ; Node index v ∈ V; Model parameters θ. Ensure: Exact embedding zv 1: z(0)i ← xi (∀i ∈ V) 2: for l ∈ {1, . . . , L} do 3: for i ∈ V do 4: h(l)i ← ∑ u∈N (i)Mliu(z (l−1) i , z (l−1) u , eiu,θ)
5: z(l)i ← Ul(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 zv simply. The aim of this study is to develop a constant time approximation algorithm for calculating the embedding vector zv and gradients ∂zv∂θ with the given model parameters θ and node v.

heading: 3.3 COMPUTATIONAL MODEL ASSUMPTIONS
text: 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: (1) Odeg(v): the degree of node v, (2) OG(v, i): the i-th neighbor of node v, and (3) Ofeature(v): the feature of node v. We assume that our algorithm can query the oracles in constant time per query.

heading: 3.4 PROBLEM FORMULATION
text: Given a node v, we calculate the following functions with the least number of oracle accesses: (1) Oz(v): the embedding zv and (2) Og(v): the gradients of parameters ∂zv∂θ . However, the exact computation of Oz and Og 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:
• Ôz(v, ε, δ): an estimate ẑv of zv such that Pr[‖ẑi − zi‖2 ≥ ε] ≤ δ,
• Ô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 Ôz and Ô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. ‖xi‖2 ≤ B, ‖eiu‖2 ≤ B, and ‖θ‖2 ≤ B. Assumption 2 deg(i)Mliu and Ul are uniformly continuous in any bounded domain. Assumption 3 (Only for gradient computation) deg(i)DMliu and DUl 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.
Assumption 4 ∃K ∈ R s.t. deg(i)Mliu and Ul are K-Lipschitz continuous in any bounded domain. Assumption 5 (Only for gradient computation) ∃K ′ ∈ R s.t. deg(i)DMliu and DUl are K ′-
Lipschitz continuous in any bounded domain.

heading: 4 PROPOSED METHOD

heading: 4.1 CONSTANT TIME EMBEDDING APPROXIMATION
text: Here, we build a constant time approximation algorithm based on neighbor sampling, which approximates the embedding zv 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 Ô(l)z (l = 1, . . . , L). The pseudo code is presented in Algorithm 2. Here, Ô(l−1)z (v ← u) means calling the function Ô(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 Ô(l)z : Estimate the embedding z(l)v Require: Graph G = (V, E) (as oracle); Features X ∈ Rn×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: ĥ(l)v ←
{ Odeg(v) r(l) ∑ u∈S(l) Mlvu(Ofeature(v),Ofeature(u), eiu,θ) (l = 1)
Odeg(v) r(l) ∑ u∈S(l) Mlvu(Ô (l−1) z (v ← v), Ô(l−1)z (v ← u), eiu,θ) (l > 1)
3: ẑ(l)v ← Ul(Ô(l−1)z (v ← v), ĥ(l)v ,θ) if l > 1 otherwise Ul(Ofeature(v), ĥ(l)v ,θ) 4: return ẑ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[‖Ôz(v, ε, δ)− zv‖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.

heading: 4.2 CONSTANT TIME GRADIENT APPROXIMATION
text: 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 zv. Let ∂zv∂θ be the gradient of the embedding zv 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 ẑ (L) v , which is obtained by Ô(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 Ô(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.

heading: 5 INAPPROXIMABILITY
text: 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 ‖xi‖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.

heading: 6 EXPERIMENTS
text: 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?
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)?
Q6: Does neighbor sampling efficiently work for real data?
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 Kn. 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 Kn, whose features are x1 = (1, 0)> and xi = (0, 1/n)> (i 6= 1), and the weight matrix is an identity matrix I2. 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 Kn, whose features are x1 = (1, 2)> and xi = (1, 1)> (i 6= 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 cliqueKn, 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., ‖xi‖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 = r2 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 Erdős-Rényi model (ER model) (Erdő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 matricesW (1) andW (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.

heading: 7 CONCLUSION
text: 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.

heading: A MODELS
text: We introduce GraphSAGE-GCN, GraphSAGE-mean, GraphSAGE-pool, the graph convolutional networks (GCN), and the graph attention networks (GAT) for completeness of our paper.
GraphSAGE-GCN (Hamilton et al., 2017): The message function and the update function of this model is
Mlvu(zv, zu, evu,θ) = zu
deg(v) ,
Ul(zv,hv,θ) = σ(W (l)hv),
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
Mlvu(zv, zu, evu,θ) = zv
deg(v) ,
Ul(zv,hv,θ) = σ(W (l)[zv,hv]),
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
Mlvu(zv, zu, evu,θ) = z(l)√
deg(v)deg(u) ,
Ul(zv,hv,θ) = σ(W (l)hv).
GAT (Veličković et al., 2018): The message function and the update function of this model is
α(l)vu = exp(LEAKYRELU(a(l)>[W (l)zv,W (l)zu]))∑
u′∈N (v) exp(LEAKYRELU(a(l)>[W (l)zv,W (l)zu′ ])) ,
Mlvu(zv, zu, evu,θ) = α (l) vuzu,
Ul(zv,hv,θ) = σ(W (l)hv).
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
ẑ(0)v = z (0) v
α̂(l)vu = exp(LEAKYRELU(a(l)>[W (l)ẑ(l−1)v ,W (l)ẑ (l−1) u ]))∑
u′∈S(l) exp(LEAKYRELU(a(l)>[W (l)ẑ (l−1) v ,W (l)ẑ (l−1) u′ ]))
,
ĥ(l)v = ∑ u∈S(l) α̂(l)vuẑ (l−1) u ,
ẑ(l)v = Ul(ẑ (l−1) v , ĥ (l) v ,θ).
Then, the following property holds true.
Pr[‖z(L)v − ẑ(L)v ‖2 ≥ ε] < δ.
We prove Proposition 11 in Section C.

heading: B EXPERIMENTAL SETUP
text: Experiments for Q1: For each algorithm and for n = 27, 28, . . . , 215, we (1) run the algorithm with an input graph, (2) we calculate the exact embedding of one node and its estimate using 100 samples (i.e., r(1) = r(2) = 100), (3) calculate the approximation error. We run the process above 10 times, and we report the mean value of time consumption.
Experiments for Q4: We initialize the weight matrix W (1) with normal distribution and then normalize it so that the operator norm ‖W (1)‖op of the matrix equal to 1. This satisfies Assumption 1 (i.e., ‖θ‖2 ≤ √ 2). For each r = 1, . . . , 10000, we (1) initialize the weight matrix, (2) choose 400 nodes, (3) calculate the exact embedding of each chosen node, (4) calculate the estimate for each chosen node with r samples (i.e., r(1) = r), and (5) calculate the approximation error of each chosen node.
Experiments for Q5: In the experiment for the AB model, we (1) iterate r from 8 to 1000, (2) set n = r2, (3) generate 10 graphs with n nodes using the BA model, (4) choose the node that has the maximum degree for each generated graph, (5) calculate the exact embeddings and its estimate for each chosen node with r samples (i.e., r(1) = r(2) = r), and (6) calculate the approximation error.
The experimental process for the ER model is similar to that for the ER model, but we (1) use the ER model instead of BA model and (2) set n = floor(r1.5) instead of n = r2 to reduce the computational cost.
Experiments for Q6: To calculate the approximation errors of embeddings, for each trained model, for r = 1 . . . 100, we (1) calculate the exact embedding of each test node, (2) calculate an estimate of embedding of each test node with r samples (i.e., r(1) = r(2) = r), and (3) calculate the approximation error of each test node. To calculate the approximation errors of gradients, for each dataset, we (1) initialize ten models with Xavier initializer (Glorot & Bengio, 2010), (2) choose random 1000 nodes, and (3) for each model, for each chosen node, and for r = 1 . . . 100, calculate the exact and approximation gradients of the classification loss with respect to the parameters, and we calculate their approximation error.

heading: C PROOFS
text: We introduce the following multivariate version of the Hoeffding’s inequality (Hoeffding, 1963) to prove the theoretical bound (Theorem 1 and 4).
Lemma 12 (multivariate Hoeffding’s inequality). Let x1,x2, . . . ,xn be independent d-dimensional random variables whose two-norms are bounded ‖xi‖2 ≤ B, and let x̄ be the empirical mean of these variables x̄ = 1n ∑n i=1 xi. Then, for any ε > 0,
Pr[‖x̄− E[x̄]‖2 ≥ ε] ≤ 2d exp ( − nε 2
2B2d ) holds true.
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 X1, X2, . . . , Xn be independent random variables bounded by the intervals [−B,B] and let X̄ be the empirical mean of these variables X̄ = 1n ∑n i=1Xi. Then, for any ε > 0,
Pr[|X̄ − E[X̄]| ≥ ε] ≤ 2 exp ( − nε 2
2B2 ) holds true.
Proof of Lemma 12. Apply Lemma 13 to each dimension k of Xi. Then,
Pr[|X̄k − E[X̄]k| ≥ ε√ d
] ≤ 2 exp ( − nε 2
2B2d
) .
It should be noted that |Xik| < B because ‖Xi‖2 < B. Therefore,
Pr[∃k ∈ {1, 2, . . . , d} |X̄k − E[X̄]k| ≥ ε√ d
] ≤ 2d exp ( − nε 2
2B2d
) .
If |X̄k − E[X̄]k| < ε√d holds true for all dimension k, then
‖X̄ − E[X̄]‖2 = √√√√ d∑ k=1 (X̄k − E[X̄]k)2 < √ d · ε 2 d = ε.
Therefore,
Pr[‖X̄ − E[X̄]‖2 ≥ ε] ≤ 2d exp ( − nε 2
2B2d
) .
First, we prove that the embedding in each layer is bounded.
Lemma 14. Under Assumptions 1 and 2, the norms of the embeddings ‖z(l)v ‖2, ‖ẑ(l)v ‖2, ‖h(l)v ‖2, and ‖ĥ(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)Mliu and the update function Ul is continuous by Assumption 2. Since the image f(X) of a compact set X ∈ Rd is compact if f is continuous, the images of deg(i)Mliu and Ul 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.
Because UL is uniform continuous, ∃ε′ > 0,∀zv,hv,h′v,θ, ‖hv − h′v‖2 < ε′ ⇒ ‖UL(zv,hv,θ)− UL(zv,h′v,θ)‖2 < ε. (1)
Let xk be the k-th sample in S(L) and Xk = deg(v)MLvu(z(0)v , z(0)xk , evxk ,θ). Then, E[Xk] = ∑
u∈N (v)
MLvu(z (0) v , z (0) u , evu,θ) = h (L) v . (2)
There exists a constant C ∈ R such that for any input satisfying Assumption 1, ‖Xk‖2 < C (3)
holds true because ‖z(0)v ‖2, ‖z(0)xk ‖2, ‖evxk‖2, and ‖θ‖2 are bounded by Assumption 1 and deg(v)MLvu is continuous. Therefore, if we take r(L) = O( 1ε′2 log 1 δ ) samples, Pr[‖ĥ (L) v −h(L)v ‖2 ≥ ε′] ≤ δ by the Hoeffding’s inequality and equations (2) and (3). Therefore, Pr[‖ẑ(L)v − z(L)v ‖2 ≥ ε] ≤ δ.
Inductive step: It is shown that the statement holds true for L = l+ 1 if it holds true for L = l. The induction hypothesis is ∀ε > 0, 1 > δ > 0, ∃r(1)(ε, δ), . . . , r(L−1)(ε, δ) such that ∀v ∈ V , Pr[‖Ô(L−1)z (v, ε, δ)− zv‖2 ≥ ε] ≤ δ.
Because UL is uniform continuous,
∃ε′ > 0,∀zv,hv, z′v,h′v,θ, ‖[zv,hv]− [z′v,h′v]‖2 < ε′ ⇒ ‖UL(zv,hv,θ)− UL(z′v,h′v,θ)‖2 < ε, (4)
where [·] denotes concatenation. By the induction hypothesis,
∃r′(1), . . . , r′(L−1) such that Pr[‖Ô(L−1)z (v)− z(L−1)v ‖2 ≥ ε′/ √
2] ≤ δ/2. (5) holds true. Let
h̃(L)v = deg(v) r(L) ∑ u∈S(L) MLvu(z (L−1) v , z (L−1) u , evu,θ).
Let xk be the k-th sample in S(L) and Xk = deg(v)MLvu(z(L−1)v , z(L−1)xk , evxk ,θ). Then, E[Xk] = ∑
u∈N (v)
MLvu(z (L−1) v , z (L−1) xk , evxk ,θ) = h (L) v . (6)
There exists a constant C ∈ R such that for any input satisfying Assumption 1, ‖Xk‖2 < C, (7)
because ‖z(L−1)v ‖2, ‖z(L−1)xk ‖2, ‖evxk‖2, and ‖θ‖2 are bounded by Assumption 1 and Theorem 14, and deg(v)MLvu is continuous. If we take r(L) = O( 1ε′2 log 1 δ ), then
Pr[‖h̃(l)v − h(l)v ‖2 ≥ ε′/(2 √
2)] ≤ δ/4, (8) by the Hoeffding’s inequality and equations (6) and (7). Because deg(v)MLvu is uniform continuous,
∃ε′′ > 0 such that ‖[z(L−1)v , z(L−1)u ]− [z′(L−1)v , z′(L−1)u ]‖2 ≤ ε′′
⇒ deg(v)‖MLvu(z(L−1)v , z(L−1)u , evu,θ)−MLvu(z′(L−1)v , z′(L−1)u , evu,θ)‖2 ≤ ε′/(2 √
2). (9)
By the induction hypothesis,
∃r′′(1), . . . , r′′(l) such that Pr[‖Ô(L−1)z (v)− z(L−1)v ‖2 ≥ ε′′/ √
2] ≤ δ/(8r(L)). (10) Therefore, the probability that the errors of all oracle calls are bounded is
Pr[∃u ∈ S(L), ‖[Ô(L−1)z (v), Ô(L−1)z (u)]− [z(L−1)v , z(L−1)u ]‖2 ≥ ε′′] ≤ δ/4. (11) By equations (9) and (11),
Pr[∃u ∈ S(L), deg(v)‖MLvu(z(L−1)v , z(L−1)u , evu,θ)−MLvu(ẑ(L−1)v , ẑ(L−1)u , evu,θ)‖2 ≥ ε′/(2 √ 2)] ≤ δ/4.
Pr[‖h̃(L)v − ĥ(L)v ‖2 ≥ ε′/(2 √
2)] ≤ δ/4. (12) By the triangular inequality and equations (8) and (12),
Pr[‖ĥ(L)v − h(L)v ‖2 ≥ ε′/ √ 2] ≤ δ/2. (13)
Therefore, if we take r(1) = max(r′(1), r′′(1)), . . . , r(L−1) = max(r′(L−1), r′′(L−1)), by equations (5) and (13),
Pr[‖[z(L−1)v ,h(L)v ]− [Ô(L−1)z (v), ĥ(L)v ]‖2 ≥ ε′] ≤ δ. (14) Therefore, by equations (9) and (14),
Pr[‖ẑ(L)v − z(L)v ‖2 ≥ ε] ≤ δ.
Proof of Theorem 2. We prove this by performing mathematical induction on the number of layers.
Base case: It is shown that the statement holds true for L = 1.
If UL is K-Lipschitz continuous, ε′ = O(ε) in equation (1). Therefore, r(L) = O( 1ε2 log 1 δ ).
Inductive step: It is shown that the statement holds true for L = l + 1 if it holds true for L = l.
If UL and MLvu are K-Lipschitz continuous, ε′ = O(ε) in equation (4) and ε′′ = O(ε) in equation (9). Therefore, r(L) = O( 1ε2 log 1 δ ). We call Ô (L−1) z (v) such that Pr[‖Ô(L−1)z (v) − z(L−1)v ‖2 ≥
ε′/ √
2] ≤ δ/2 in equation (5). Therefore, r′(1), . . . , r′(L−1) = O( 1ε2 (log 1 ε + log 1 δ ) are sufficient
by the induction hypothesis. We call Ô(L−1)z (v) such that Pr[‖Ô(L−1)z (v)− z(L−1)v ‖2 ≥ ε′′/ √
2] ≤ δ/(8r(L)) in equation (10). Therefore, r′′(1), . . . , r′′(L−1) = O( 1ε2 (log 1 ε + log 1 δ ) are sufficient by the induction hypothesis because log 1 δ/(8r(L)) = O(log 1ε + log 1 δ ). In total, the complexity is O( 1 ε2L (log 1ε + log 1 ε ) L−1 log 1δ )
Lemma 15 (Chazelle et al. (2005)). Let Ds be Bernoulli( 1+sε2 ). Let n-dimentional distribution D be (1) pick s = 1 with probability 1/2 and s = −1 otherwise; (2) then draw n values from Ds. Any probabilistic algorithm that can guess the value of s with a probability error below 1/4 requires Ω( 1ε2 ) bit lookup on average.
Proof of Theorem 3. We prove there is a counter example in the GraphSAGE-GCN models. Suppose there is an algorithm that approximates the one-layer GraphSAGE-GCN within o(ε2) queries. We prove that this algorithm can distinguish D in Lemma 15 within o(ε2) queries and derive a contradiction.
Let σ be any non-constant K-Lipschitz activation function. There exists a, b ∈ R (a > b) such that σ(a) 6= σ(b) because σ is not constant. Let S = |σ(a)−σ(b)|a−b > 0. Let ε > 0 be any sufficiently small positive value and t ∈ {0, 1}n be a random variable drawn from D. We prove that we can determine s with high provability within o(ε2) queries using the algorithm. Let G be a clique Kn andW (1) = 1. Let us calculate aε and bε using the following steps: (1) set aε = a and bε = b; (2) if aε − bε < ε, return aε and bε; (3) m = aε+bε2 ; (4) if |σ(aε) − σ(m)| > |σ(m) − σ(bε)|, then set bε = m, otherwise aε = m; and (5) go back to (2). Here, ε/2 ≤ aε − bε < ε, a ≤ aε+bε2 ≤ b, and |σ(aε) − σ(bε)| ≥ S2 ε hold true. Let xv = aε+bε 2 + (2tv − 1) aε−bε 2ε for all v ∈ V . Then, E[hv | s = 1] = aε and E[hv | s = −1] = bε. Therefore, Pr[|zv − σ(aε)| < S8 ε | s = 1] → 1 as n→∞ and Pr[|zv − σ(bε)| < S8 ε | s = −1]→ 1 as n→∞ because σ is K-Lipschitz. We set the error tolerance to S8 ε and n to a sufficiently large number. Then s = 1 if |ẑv − σ(aε)| < S 4 ε and s = −1 otherwise with high probability. However, the algorithm accesses t (i.e., accesses Ofeature) o(ε2) times. This contradicts with Lemma 15.
Lemma 16. Under Assumptions 1, 2 and 3, the norms of the gradients of the message functions and the update functions ‖DUl(z(l−1)v ,h(l)v ,θ)‖F , ‖DUl(ẑ(l−1)v , ĥ(l)v ,θ)‖F , ‖deg(v)DMlvu(z(l−1)v ,v(l)u , evu,θ)‖F , and ‖deg(v)DMlvu(ẑ(l−1)v , v̂(l)u , evu,θ)‖F are bounded by a constant B′ ∈ R.
Proof of Lemma 16. The input of each function is bounded by Lemma 14. Because DUl and deg(v)DMlvu is uniform continuous, these images are bounded.
Proof of Theorem 4. 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.
When the number of layers is one,
∂z (L) v
∂θ = ∂UL ∂θ (z(0)v ,h (L) v ,θ) + ∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)
∂h (L) v
∂θ
= ∂UL ∂θ (z(0)v ,h (L) v ,θ) + ∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ) ∑ u∈N (v) ∂MLvu ∂θ (z(0)v , z (0) u , evu,θ).
∂ẑ (L) v
∂θ = ∂UL ∂θ (z(0)v , ĥ (L) v ,θ) + ∂UL
∂ĥ (L) v
(z(0)v , ĥ (L) v ,θ)
∂ĥ (L) v
∂θ
= ∂UL ∂θ (z(0)v , ĥ (L) v ,θ) + ∂UL
∂ĥ (L) v
(z(0)v , ĥ (L) v ,θ) deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂θ (z(0)v , z (0) u , evu,θ).
Because DUL is uniform continuous,
∃ε′ > 0 s.t. for all input,‖h(L)v − ĥ(L)‖2 < ε′ ⇒
(‖∂UL ∂θ (z(0)v ,h (L) v ,θ)− ∂UL ∂θ (z(0)v , ĥ (L) v ,θ)‖F < ε/2 ∧
‖ ∂UL ∂h (L) v (z(0)v ,h (L) v ,θ)− ∂UL ∂ĥ (L) v (z(0)v , ĥ (L) v ,θ)‖F < ε/(4B′)). (15)
If we take r(L) = O( 1ε′2 log 1 δ ),
Pr[‖h(L)v − ĥ(L)‖2 ≥ ε′] ≤ δ/2 (16)
holds true for any input by the argument of the proof of Theorem 1.
Let xk be the k-th sample in S(L) and Xk = deg(v)∂MLvu∂θ (z (0) v , z (0) u , evu,θ). Then,
E[Xk] = ∑
u∈N (v)
∂MLvu ∂θ (z(0)v , z (0) u , evu,θ) =
∂h (L) v
∂θ . (17)
There exists a constant C ∈ R such that for any input satisfying Assumption 1,
‖Xk‖2 < C, (18)
because ‖z(0)v ‖2, ‖z(0)xk ‖2, ‖evxk‖2, and ‖θ‖2 are bounded by Assumption 1 and deg(v)DMLvu is continuous. Therefore, if we take r(L) = O( 1ε2 log 1 δ ) samples,
Pr[‖∂ĥ (L) v
∂θ − ∂h
(L) v
∂θ ‖F ≥ ε/(4B′)] ≤ δ/2 (19)
holds true by the Hoeffding’s inequality and equations (17) and (18). If
‖∂ĥ (L) v
∂θ − ∂h
(L) v
∂θ ‖F < ε/(4B′)
and
‖ ∂UL ∂h (L) v (z(0)v ,h (L) v ,θ)− ∂UL ∂ĥ (L) v (z(0)v , ĥ (L) v ,θ)‖F < ε/(4B′))
hold true, then
‖∂ĥ (L) v
∂θ
∂UL
∂ĥ (L) v
(z(0)v , ĥ (L) v ,θ)−
∂h (L) v
∂θ
∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)‖F
≤ ‖∂ĥ (L) v
∂θ
∂UL
∂ĥ (L) v
(z(0)v , ĥ (L) v ,θ)−
∂ĥ (L) v
∂θ
∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)‖F
+ ‖∂ĥ (L) v
∂θ
∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)−
∂h (L) v
∂θ
∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)‖F
= ‖∂ĥ (L) v
∂θ ‖F ‖
∂UL
∂ĥ (L) v
(z(0)v , ĥ (L) v ,θ)−
∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)‖F
+ ‖∂ĥ (L) v
∂θ − ∂h
(L) v
∂θ ‖F ‖
∂UL
∂h (L) v
(z(0)v ,h (L) v ,θ)‖F
≤ B′‖ ∂UL ∂ĥ (L) v (z(0)v , ĥ (L) v ,θ)− ∂UL ∂h (L) v (z(0)v ,h (L) v ,θ)‖F +B′‖
∂ĥ (L) v
∂θ − ∂h
(L) v
∂θ ‖F
< B′ ε
4B′ +B′
ε 4B′ = ε 2 .
Therefore, Pr[‖∂z (L) v ∂θ − ∂z(L)v ∂θ ‖F ≥ ε] ≤ δ by equations (15), (16), and (19).
Inductive step: It is shown that the statement holds true for L = l + 1 if it holds true for L = l.
∂z (L) v
∂θ = ∂UL ∂θ (z(L−1)v ,h (L) v ,θ) + ∂UL
∂z (L−1) v
(z(L−1)v ,h (L) v ,θ)
∂z (L−1) v
∂θ
+ ∂UL
∂h (L) v
(z(L−1)v ,h (L) v ,θ) ∑ u∈N (v) ∂MLvu ∂θ (z(L−1)v , z (L−1) u , evu,θ)
+ ∂UL
∂h (L) v
(z(L−1)v ,h (L) v ,θ) ∑ u∈N (v) ∂MLvu ∂z (L−1) v (z(L−1)v , z (L−1) u , evu,θ) ∂z (L−1) v ∂θ
+ ∂UL
∂h (L) v
(z(L−1)v ,h (L) v ,θ) ∑ u∈N (v) ∂MLvu ∂z (L−1) u (z(L−1)v , z (L−1) u , evu,θ) ∂z (L−1) u ∂θ .
∂ẑ (L) v
∂θ = ∂UL ∂θ (ẑ(L−1)v , ĥ (L) v ,θ) + ∂UL
∂ẑ (L−1) v
(ẑ(L−1)v , ĥ (L) v ,θ)
∂ẑ (L−1) v
∂θ
+ ∂UL
∂ĥ (L) v
(ẑ(L−1)v , ĥ (L) v ,θ) deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂θ (ẑ(L−1)v , ẑ (L−1) u , evu,θ)
+ ∂UL
∂ĥ (L) v
(ẑ(L−1)v , ĥ (L) v ,θ) deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂ẑ (L−1) v (ẑ(L−1)v , ẑ (L−1) u , evu,θ) ∂ẑ (L−1) v ∂θ
+ ∂UL
∂ĥ (L) v
(ẑ(L−1)v , ĥ (L) v ,θ) deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂ẑ (L−1) u (ẑ(L−1)v , ẑ (L−1) u , evu,θ) ∂ẑ (L−1) u ∂θ .
Because DUL is uniform continuous,
∃ε′ > 0 s.t. for all input,‖[z(L−1)v ,h(L)v ]− [ẑ(L−1)v , ĥ(L)v ]‖2 < ε′ ⇒
(‖∂UL ∂θ (z(L−1)v ,h (L) v ,θ)− ∂UL ∂θ (ẑ(L−1)v , ĥ (L) v ,θ)‖F < O(ε) ∧
‖ ∂UL ∂h (L) v (z(L−1)v ,h (L) v ,θ)− ∂UL ∂ĥ (L) v (ẑ(L−1)v , ĥ (L) v ,θ)‖F < O(ε) ∧
‖ ∂UL ∂z (L−1) v (z(L−1)v ,h (L) v ,θ)− ∂UL ∂ẑ (L−1) v (ẑ(L−1)v , ĥ (L) v ,θ)‖F < O(ε)). (20)
Because deg(v)DMLvu is uniform continuous,
∃ε′′ > 0 s.t. for all input,‖[z(L−1)v , z(L−1)u ]− [ẑ(L−1)v , ẑ(L−1)u ]‖2 < ε′′ ⇒
(deg(v)‖∂MLvu ∂θ (z(L−1)v , z (L−1) u , evu,θ)− ∂UL ∂θ (ẑ(L−1)v , ẑ (L−1) u , evu,θ)‖F < O(ε) ∧
deg(v)‖∂MLvu ∂h (L) v (z(L−1)v , z (L−1) u , evu,θ)− ∂UL ∂ĥ (L) v (ẑ(L−1)v , ẑ (L−1) u , evu,θ)‖F < O(ε) ∧
deg(v)‖ ∂MLvu ∂z (L−1) v (z(L−1)v , z (L−1) u , evu,θ)− ∂UL ∂ẑ (L−1) v (ẑ(L−1)v , ẑ (L−1) u , evu,θ)‖F < O(ε)). (21)
By the argument of the proof of Theorem 1, if we take sufficiently large number of samples,
Pr[‖z(L−1)v − ẑ(L−1)v ‖2 ≥ O(min(ε′, ε′′))] ≤ O(εδ), (22)
Pr[‖h(L)v − ĥ(L)v ‖2 ≥ O(ε′)] ≤ O(δ). (23)
By the induction hypothesis, there exists r(1), . . . , r(L−1) such that
Pr[‖∂z (L−1) u
∂θ − ∂ẑ
(L−1) u ∂θ ‖F ≥ O(ε)] ≤ O(εδ),
Pr[‖∂z (L−1) v
∂θ − ∂ẑ
(L−1) v ∂θ ‖F ≥ O(ε)] ≤ O(εδ),
If we take r(L) = O( 1ε2 log 1 δ ),
Pr[‖ ∑
u∈N (v)
∂MLvu ∂θ (z(L−1)v , z (L−1) u , evu,θ)
− deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂θ (z(L−1)v , z (L−1) u , evu,θ)‖ ≥ O(ε)] ≤ O(δ), (24)
Pr[‖ ∑
u∈N (v)
∂MLvu ∂z (L−1) v (z(L−1)v , z (L−1) u , evu,θ)
∂z (L−1) v
∂θ
− deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂z (L−1) v (z(L−1)v , z (L−1) u , evu,θ) ∂z (L−1) v ∂θ ‖ ≥ O(ε)] ≤ O(δ), (25)
Pr[‖ ∑
u∈N (v)
∂MLvu ∂z (L−1) u (z(L−1)v , z (L−1) u , evu,θ)
∂z (L−1) u
∂θ
− deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂z (L−1) u (z(L−1)v , z (L−1) u , evu,θ) ∂z (L−1) u ∂θ ‖ ≥ O(ε)] ≤ O(δ), (26)
holds true by the Hoeffding’s inequality. Therefore,
Pr[‖ ∑
u∈N (v)
∂MLvu ∂θ (z(L−1)v , z (L−1) u , evu,θ)
− deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂θ (ẑ(L−1)v , ẑ (L−1) u , evu,θ)‖ ≥ O(ε)] ≤ O(δ), (27)
Pr[‖ ∑
u∈N (v)
∂MLvu ∂z (L−1) v (z(L−1)v , z (L−1) u , evu,θ)
∂z (L−1) v
∂θ
− deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂ẑ (L−1) v (ẑ(L−1)v , ẑ (L−1) u , evu,θ) ∂ẑ (L−1) v ∂θ ‖ ≥ O(ε)] ≤ O(δ), (28)
Pr[‖ ∑
u∈N (v)
∂MLvu ∂z (L−1) u (z(L−1)v , z (L−1) u , evu,θ)
∂z (L−1) u
∂θ
− deg(v) r(L) ∑ u∈S(L) ∂MLvu ∂ẑ (L−1) u (ẑ(L−1)v , ẑ (L−1) u , evu,θ) ∂ẑ (L−1) u ∂θ ‖ ≥ O(ε)] ≤ O(δ), (29)
holds true by equations (20), (21), (22), (23), (24), (25), and (26).
Therefore, if we take r(1), . . . r(L) sufficiently large, Pr[‖∂z (L) v ∂θ − ∂z(L)v ∂θ ‖F ≥ ε] ≤ δ holds true by equations (20), (21), (22), (23), (27), (28), and (29).
Proof of Theorem 5. We prove this by performing mathematical induction on the number of layers.
Base case: It is shown that the statement holds true for L = 1.
If DUL is K ′-Lipschitz continuous, ε′ = O(ε) in equation (15). Therefore, r(L) = O( 1ε2 log 1 δ ) is sufficient.
Inductive step: It is shown that the statement holds true for L = l + 1 if it holds true for L = l.
If DUL and deg(v)DULvu is K ′-Lipschitz continuous, ε′ = O(ε) in equation (20) and ε′′ = O(ε) in equation (21). Therefore, r(L) = O( 1ε2 log 1 δ ) and r (1), . . . , r(L−1) = O( 1ε2 (log 1 ε + log 1 δ )) are sufficient.
Proof of Theorem 6. We show that one-layer GraphSAGE-GCN whose activation function is not constant cannot be approximated in constant time if ‖xv‖2 or ‖θ‖2 are not bounded. There exists a ∈ R such that σ(a) 6= σ(0) because σ is not constant. We consider the following two types of inputs:
• G is the clique Kn,W (1) = 1, and xi = 0 for all nodes i ∈ V .
• G is the clique Kn,W (1) = 1, xi = 0(i 6= v) for some v ∈ V , and xv = an.
Then, for the former input, z(1)v = σ(0). For the latter type of inputs, z (1) v = σ(a). Let A be an arbitrary constant algorithm and C be the number of queries A makes when we set ε = |σ(a)− σ(0)|/3. When A calculates the embedding of u 6= v ∈ V , the states of all nodes but u are symmetrical until A makes a query about that node. Therefore, if n is sufficiently large, A does not make any query about v with high probability (i.e., at least (1 − 1n−1 )
C). If A does not make any query about v, the state of A is the same for both types of inputs. If the approximation error is less than ε for the first type of inputs, the approximation error is larger than ε for the second type of inputs by the triangle inequality and vice versa. Therefore, A fails to approximate the embeddings of either type of inputs with the absolute error of at most ε. As for θ, we setW (1) = an and xv = 1 for the second type of inputs. Then, the same argument follows.
Proof of Theorem 7. We consider the one-layer GraphSAGE-GCN with ReLU and normalization (i.e., σ(x) = RELU(x)/‖RELU(x)‖2). We use the following two types of inputs:
• G is the clique Kn, W (1) is the identity matrix I2, xi = (0, 0)>(i 6= v) for some node v ∈ V , and xv = (1, 0)>.
• G is the clique Kn, W (1) is the identity matrix I2, xi = (0, 0)>(i 6= v) for some node v ∈ V , and xv = (0, 1)>.
Then, for the former type of inputs, hi = (1/n, 0)>, zi = (1, 0)>, and ∂zi2∂W21 = 1 for all i ∈ V . For the latter type of inputs, hi = (0, 1/n)>, zi = (0, 1)>, and ∂zi2∂W21 = 0 for all i ∈ V . Let A be an arbitrary constant algorithm and C be the number of queries A makes when we set ε = 1/3. When A calculates the embedding or gradient of u 6= v ∈ V , the states of all nodes but u are symmetrical until A makes a query about that node. Therefore, if n is sufficiently large, A does not make any query about v with high probability (i.e., at least (1− 1n−1 )
C). If A does not make any query about v, the state of A is the same for both types of inputs. If the approximation error is less than ε for the first type of inputs, the approximation error is larger than ε for the second type of inputs by the triangle inequality and vice versa. Therefore, A fails to approximate the embeddings and gradients of either type of inputs with the absolute error of at most ε.
Proof of Theorem 8. We consider the one-layer GraphSAGE-GCN with ReLU (i.e., σ(x) = RELU(x)). We use the following two types of inputs:
• G is the clique Kn, W (1) = (−1, 1), xi = (1, 1)>(i 6= v) for some node v ∈ V , and xv = (1, 2) >.
• G is the clique Kn, W (1) = (−1, 1), xi = (1, 1)>(i 6= v) for some node v ∈ V , and xv = (1, 0) >.
Then, for the former type of inputs, MEAN({xu | u ∈ N (v)}) = (1, 1 + 1n ) >, hv = zv = 1n , and ∂zv∂W = (1, 1 + 1 n ) for all i ∈ V . For the latter type of inputs, MEAN({xu | u ∈ N (v)}) = (1, 1 − 1n ) >, hv = − 1n , zv = 0, and ∂zv ∂W = (0, 0) for all i ∈ V . Let A be an arbitrary constant algorithm and C be the number of queries A makes when we set ε = 1/3. When A calculates the gradient of u 6= v ∈ V , the states of all nodes but u are symmetrical until A makes a query about that node. Therefore, if n is sufficiently large, A does not make any query about v with high probability (i.e., at least (1− 1n−1 )
C ). If A does not make any query about v, the state of A is the same for both types of inputs. If the approximation error is less than ε for the first type of inputs, the approximation error is larger than ε for the second type of inputs by the triangle inequality and vice versa. Therefore, A fails to approximate the gradients of either type of inputs with the absolute error of at most ε.
Proof of Theorem 9. We consider the one-layer GraphSAGE-pool whose activation function satisfies σ(1) 6= σ(0) and the following two types of inputs:
• G is the clique Kn,W (1) = 1, b = 0, and xi = 0 for all nodes v ∈ V .
• G is the clique Kn,W (1) = 1, b = 0, xi = 0 (i 6= v) for some node v ∈ V , and xv = 1.
Then, for the former type of inputs, zi = σ(0) for all i ∈ V . For the latter type of inputs, zi = σ(1) for all i ∈ V . LetA be an arbitrary constant algorithm and C be the number of queriesAmakes when we set ε = |σ(1)− σ(0)|/3. When A calculates the embedding of u 6= v ∈ V , the states of all nodes but u are symmetrical until A makes a query about that node. Therefore, if n is sufficiently large, A does not make any query about v with high probability (i.e., at least (1− 1n−1 )
C ). IfA does not make any query about v, the state of A is the same for both types of inputs. If the approximation error is less than ε for the first type of inputs, the approximation error is larger than ε for the second type of inputs by the triangle inequality and vice versa. Therefore, A fails to approximate the embeddings of either type of inputs with the absolute error of at most ε.
Proof of Lemma 10. We consider the one-layer GCN whose activation function satisfies σ(1) 6= σ(0). We use the following two types of inputs:
• G is a star graph, where v ∈ V is the center of G,W (1) = 1, and all features are 0. • G is a star graph, where v ∈ V is the center of G,W (1) = 1, and the features of √
2n leafs are 1 and the features of other nodes are 0.
Then, for the former type of inputs, zv = σ(0). For the latter type of inputs, zv = σ(1). Let A be an arbitrary constant algorithm and C be the number of queries A makes when we set ε = |σ(1)− σ(0)|/3. When A calculates the embedding of u ∈ V that xu = 0, the states of all nodes but u are symmetrical until A makes a query about that node. Therefore, if n is sufficiently large, A does not make any query about v with high probability (i.e., at least (1 − √ 2n
n−1 ) C). If A does not make
any query about v, the state of A is the same for both types of inputs. If the approximation error is less than ε for the first type of inputs, the approximation error is larger than ε for the second type of inputs by the triangle inequality and vice versa. Therefore, A fails to approximate the embeddings of either type of inputs with the absolute error of at most ε.
We prove the following lemma to prove Proposition 11.
Lemma 17. If Assumptions 1 holds true and σ is Lipschitz continuous, ‖z(l)v ‖2 and ‖z(l)v ‖2 (l = 1, . . . , L) of the GAT model are bounded by a constant
Proof of Lemma 17. We prove this by performing mathematical induction on the number of layers. The norm of the input of the first layer is bounded by Assumption 1. If ‖z(l−1)u ‖2 and ‖ẑ(l−1)u ‖2 are bounded for all u ∈ V , ‖h(l)v ‖2 and ‖ĥ(l)v ‖2 are bounded because h(l)v and ĥ(l)v are the weighted sum of z(l−1)u and ẑ (l−1) u . Therefore, ‖z(l)u ‖2 and ‖ẑ(l)u ‖2 are bounded because Ul is continuous.
Proof of Proposition 11. 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.
Because UL is Lipschitz continuous,
∀zv,hv,h′v,θ, ‖hv − h′v‖2 < O(ε)⇒ ‖UL(zv,hv,θ)− UL(zv,h′v,θ)‖2 < ε. (30)
Let eu = exp(LEAKYRELU(a(l)>[W (0)z (0) v ,W (0)z (0) u u])). Then,
ĥ(L)v − h(L)v = ∑ u∈S(L) α̂vuz (0) u − ∑ u∈N (v) αvuz (0) u
= 1
r(L) ∑ u∈S(L) eu 1 r(L) ∑ u′∈S(L) eu′ z(0)u − 1 deg(v) ∑ u∈N (v) eu 1 deg(v) ∑ u′∈N (v) eu′ z(0)u .
Let xk be the k-th sample in S(L) and Xk = exk . Then,
E[Xk] = 1 N (v) ∑
u∈N (v)
eu. (31)
There exists a constant c > 0, C > 0 such that for any input satisfying Assumption 1,
c < |Xk| < C, (32)
because ‖z(0)v ‖2, ‖z(0)xk ‖2, ‖W (0)‖F , and ‖a(0)‖2 are bounded by Assumption 1. Therefore, if we take r(L) = O( 1ε2 log 1 δ ) samples,
Pr[| 1 r(L) ∑ u∈S(L) eu − 1 N (v) ∑ u∈N (v) eu| ≥ O(ε)] ≤ O(δ)
by the Hoeffding’s inequality and equations (31) and (32). Because f(x) = 1/x is Lipschitz continuous in x > c > 0,
Pr[| 11 r(L) ∑ u∈S(L) eu − 11 N (v) ∑ u∈N (v) eu | ≥ O(ε)] ≤ O(δ) (33)
Let Yk =
exk 1 N (v) ∑ u′∈N (v) eu′ z(0)xk .
Then,
E[Yk] = 1 N (v) ∑
u∈N (v)
eu 1 N (v) ∑ u′∈N (v) eu′ z(0)u . (34)
There exists a constant C ′ ∈ R such that for any input satisfying Assumption 1, ‖Yk‖2 < C ′ (35)
holds true because ‖z(0)u ‖2 are bounded, and c < |eu| < C. Therefore, if we take r(L) = O( 1ε2 log 1 δ ) samples,
Pr[‖ 1 r(L) ∑ u∈S(L) eu 1 deg(v) ∑ u′∈N (v) eu′ z(0)u
− 1 deg(v) ∑ u∈N (v) eu 1 deg(v) ∑ u′∈N (v) eu′ z(0)u ‖2 ≥ O(ε)] ≤ O(δ) (36)
holds true by the Hoeffding’s inequality and equations (34) and (35). Therefore,
Pr[‖ 1 r(L) ∑ u∈S(L) eu 1 r(L) ∑ u′∈S(L) eu′ z(0)u
− 1 deg(v) ∑ u∈N (v) eu 1 deg(v) ∑ u′∈N (v) eu′ z(0)u ‖2 ≥ O(ε)] ≤ O(δ) (37)
holds true by the triangle inequality and equations (33) and (36), and Pr[‖ẑ(L)v − z(L)v ‖2 ≥ ε] ≤ δ holds true by equations (30) and (37).
Inductive step: It is shown that the statement holds true for L = l + 1 if it holds true for L = l.
Because UL is Lipschitz continuous,
∀zv,hv,h′v,θ, ‖hv − h′v‖2 < O(ε)⇒ ‖UL(zv,hv,θ)− UL(zv,h′v,θ)‖2 < ε (38)
holds true.
Pr[‖ 1 r(L) ∑ u∈S(L) eu 1 r(L) ∑ u′∈S(L) eu′ z(L−1)u
− 1 deg(v) ∑ u∈N (v) eu 1 deg(v) ∑ u′∈N (v) eu′ z(L−1)u ‖2 ≥ O(ε)] ≤ O(δ) (39)
holds true by the same argument as the base step. If we take r(1), . . . , r(L−1) = O( 1ε2 (log 1 ε +log 1 δ )) samples,
Pr[‖ẑ(L−1)u − z(L−1)u ‖2 ≥ O(ε)] ≤ O(εδ) (40)
holds true by the induction hypothesis. Therefore, Pr[‖ẑ(L)v − z(L)v ‖2 ≥ ε] ≤ δ holds true by equations (38), (39), and (40).

heading: D COMPUTATIONAL MODEL ASSUMPTIONS
text: In this study, we model our algorithm as an oracle machine that can make queries about the input and measure the complexity by query complexity. Modeling our algorithm as an oracle machine and measuring the complexity by query complexity are reasonable owing to the following reasons:
• In a realistic setting, data is stored in a storage or cloud and we may not be able to load all the information of a huge network on to the main memory. Sometimes the network is constructed for access to the information on demand (e.g., in web graph mining, the edge information is retrieved when queried). In such cases, reducing the number of queries is crucial because accessing storage or cloud is very expensive. • Our algorithm executes a constant number of elementary operations of O(log n) bits (e.g.,
accessing the O(n)-th address, sampling one element from O(n) elements). Therefore, if we assume that these operations can be done in constant time, the total computational complexity of our algorithms will be constant. This assumption is natural because most computers in the real-world can handle 64 bit integers at once and most of the network data contain less than 264 ≈ 1019 nodes. • Even if the above assumption is not satisfied, our algorithm can be executed inO(log n) time
in terms of the strict meaning of computational complexity. This indicates that our algorithm is still sub-linear and therefore it scales well. It should be noted that it is impossible to access even a single node in o(log n) time in the strict meaning of computational complexity because we cannot distinguish n nodes with o(log n) bits. Therefore, our algorithm has optimal complexity with respect to the number of nodes n.

heading: E GRAPH EMBEDDING
text: In this study, we model our algorithm as an oracle machine that can make queries about the input and measure the complexity by query complexity. Modeling our algorithm as an oracle machine and measuring the complexity by query complexity are reasonable owing to the following reasons:
• In a realistic setting, data is stored in a storage or cloud, and we may not be able to load all the information of a huge network on to the main memory. Sometimes the network is constructed for access to the information on demand (e.g., in web graph mining, the edge information is retrieved when queried). In such cases, reducing the number of queries is crucial because accessing storage or cloud is very expensive. • Our algorithm executes a constant number of elementary operations of O(log n) bits (e.g.,
accessing the O(n)-th address, sampling one element from O(n) elements). Therefore, if we assume that these operations can be done in constant time, the total computational complexity of our algorithms will be constant. This assumption is natural because most computers in the real-world can handle 64 bit integers at once, and most of the network data contain less than 264 ≈ 1019 nodes. • Even if the above assumption is not satisfied, our algorithm can be executed inO(log n) time
in terms of the strict meaning of computational complexity. This indicates that our algorithm is still sub-linear and therefore, it scales well. It should be noted that it is impossible to access even a single node in o(log n) time in the strict meaning of computational complexity because we cannot distinguish n nodes with o(log n) bits. Therefore, our algorithm has optimal complexity with respect to the number of nodes n.

heading: F TIME COMPLEXITY
text: We summarize the time complexity of approximation and exact algorithms in Table 2. Let BL = {v} and Bl = ⋃ u∈Bl+1 N (u) (l = 0, . . . , L− 1). In other words, BL−k is the k-hop neighbors of node v. We need all features of B0 to calculate the exact embedding of node v. If the graph is sparse, the size of BL−k grows exponentially with respect to k because deg(u) nodes are added to BL−k+1 for each node u ∈ BL−k. Namely, it is bounded by ∆k. If the graph is dense, BL−k ≈ V (1 ≤ k ≤ L) Therefore, complexity is linear with respect to the number of layers L and edges m. In contrast, the approximation algorithm runs in constant time irrespective of the density of the input graph.

Review summaries:
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.

In this paper, the authors provide a theoretical framework for characterizing the approximation guarantees provided by node sampling to estimate embeddings in various GNN architectures. In particular, they prove several PAC learning-style bounds on the embedding and gradient estimation when using node sampling approaches. They also observe that since the number of nodes selected for sampling is not dependent on the size of the graph, this amounts to a constant time operation for determining embedding and gradient estimates.