LEARN LOCALLY, CORRECT GLOBALLY: A DISTRIBUTED ALGORITHM FOR TRAINING GRAPH NEURAL NETWORKS
Year: 2022

Abstract:
Despite the recent success of Graph Neural Networks (GNNs), training GNNs on large graphs remains challenging. The limited resource capacities of the existing servers, the dependency between nodes in a graph, and the privacy concern due to the centralized storage and model learning have spurred the need to design an effective distributed algorithm for GNN training. However, existing distributed GNN training methods impose either excessive communication costs or large memory overheads that hinders their scalability. To overcome these issues, we propose a communication-efficient distributed GNN training technique named Learn Locally, Correct Globally (LLCG). To reduce the communication and memory overhead, each local machine in LLCG first trains a GNN on its local data by ignoring the dependency between nodes among different machines, then sends the locally trained model to the server for periodic model averaging. However, ignoring node dependency could result in significant performance degradation. To solve the performance degradation, we propose to apply Global Server Corrections on the server to refine the locally learned models. We rigorously analyze the convergence of distributed methods with periodic model averaging for training GNNs and show that naively applying periodic model averaging but ignoring the dependency between nodes will suffer from an irreducible residual error. However, this residual error can be eliminated by utilizing the proposed global corrections to entail fast convergence rate. Extensive experiments on real-world datasets show that LLCG can significantly improve the efficiency without hurting the performance.

heading: 1 INTRODUCTION
text: In recent years, Graph Neural Networks (GNNs) have achieved impressive results across numerous graph-based applications, including social networks (Hamilton et al., 2017; Deng et al., 2019),
recommendation systems (Ying et al., 2018; Wang et al., 2018), and drug discovery (Fout et al., 2017; Do et al., 2019). Despite their recent success, effective training of GNNs on large-scale real-world graphs, such as Facebook social network (Boldi & Vigna, 2004), remains challenging. Although several attempts have been made to scale GNN training by sampling techniques (Hamilton et al., 2017; Zou et al., 2019; Zeng et al., 2020; Chiang et al., 2019; Chen et al., 2018; Zhang et al., 2021), they are still inefficient for training on extremely large graphs, due to the unique structure of GNNs and the limited memory capacity/bandwidth of current servers. One potential solution to tackle these limitations is employing distributed training with data parallelism, which have become almost a de facto standard for fast and accurate training for natural language processing (Lin
et al., 2021; Hard et al., 2018) and computer vision (Bonawitz et al., 2019; Konečnỳ et al., 2018). For
∗Equal Contribution.
example, as shown in Figure 1, moving from single machine to multiple machines reduces the training time and alleviates the memory burden on each machine. Besides, scaling the training of GNNs with sampling techniques can result in privacy concerns: existing sampling-based methods require centralized data storage and model learning, which could result in privacy concerns in real-world scenarios (Shin et al., 2018; Wu et al., 2021). Fortunately, the privacy in distributed learning can be preserved by avoiding mutual access to data between different local machines, and using only a trusted third party server to access the entire data.
Nonetheless, generalizing the existing data parallelism techniques of classical distributed training to the graph domain is non-trivial, which is mainly due to the dependency between nodes in a graph. For example, unlike solving image classification problems where images are mutually independent, such that we can divide the image dataset into several partitions without worrying about the dependency between images; GNNs are heavily relying on the information inherent to a node and its neighboring nodes. As a result, partitioning graph leads to subgraphs with edges spanning subgraphs (cut-edges), which will cause information loss and hinder the performance of the model (Angerd et al., 2020). To cope with this problem, (Md et al., 2021; Jiang & Rumi, 2021; Angerd et al., 2020) propose to transfer node features and (Zheng et al., 2020; Tripathy et al., 2020; Scardapane et al., 2020) propose to transfer both the node feature and its hidden embeddings between local machines, both of which can cause significant storage/communication overhead and privacy concerns (Shin et al., 2018; Wu et al., 2021).
To better understand the challenge of distributed GNN training, we compare the validation F1score in Figure 2 (a) and the average data communicated per round in Figure 2 (b) for two different distributed GNN training methods on the Reddit dataset. On the one hand, we can observe that when ignoring the cut-edges, Parallel SGD with Periodic Averaging (PSGD-PA (Dean et al., 2012; Li et al., 2020b)) suffers from significant accuracy drop and cannot achieve the same accuracy as the single machine training, even by increasing the number of communication. However, Global Graph Sampling (GGS) can successfully reach the baseline by considering the cut-edges and allowing feature transfer, at the cost of significant communication overhead, and potential violation of privacy.
In this paper, we propose a communication-efficient distributed GNN training method, called Learn Locally, Correct Globally (LLCG). To reduce the communication overhead, inspired by the recent success of the distributed optimization with periodic averaging (Stich, 2019; Yu et al., 2019), we propose Local Training with Periodic Averaging: where each local machine first locally trains a GNN model by ignoring the cut-edges, then sends the trained model to the server for periodic model averaging, and receive the averaged model from server to continue the training. By doing so we eliminate the features exchange phase between server and local machines, but it can result in a significant performance degradation due to the lack of the global graph structure and the dependency between nodes among different machines. To compensate for this error, we propose a Global Server Correction scheme to take advantage of the available global graph structure on the server and refine the averaged locally learned models before sending it back to each local machine. Notice that without Global Server Correction, LLCG is similar to PSGD-PA as introduced in Figure 2.
To get a deeper understanding on the necessity of Global Server Correction, we provide the first theoretical analysis on the convergence of distributed training for GNNs with periodic averaging. In particular, we show that solely averaging the local machine models and ignoring the global graph structure will suffer from an irreducible residual error, which provides sufficient explanation on why Parallel SGD with Periodic Averaging can never achieve the same performance as the model trained on a single machine in Figure 2 (a). Then, we theoretically analyze the convergence of our proposal LLCG . We show that by carefully choosing the number of global correction steps, LLCG can overcome the aforementioned residual error and enjoys O ( 1/ √ PT )
convergence rate with P local machines and T iterations of gradient updates, which matches the rate of (Yu et al.,
2019) on a general (not specific for GNN training) non-convex optimization setting. Finally, we conduct comprehensive evaluations on real-world graph datasets with ablation study to validate the effectiveness of LLCG and its improvements over the existing distributed methods.
Related works. Recently, several attempts have been made on distributed GNN training. According to how they deal with the input/hidden feature of nodes that are associated with the cut-edges (i.e., the edges spanning subgraphs of each local machine), existing methods can be classified into two main categories : (1) Input feature only communication-based methods: In these methods, each local machine receives the input features of all nodes required for the gradient computation from other machines, and trains individually. However, since the number of required nodes grows exponentially with the number of layers, these methods suffer from a significant communication and storage overhead. To alleviate these issues, (Md et al., 2021) proposes to split the original graph using a min-cut graph partition algorithm that can minimize the number of cut-edges. (Jiang & Rumi, 2021) proposes to use importance sampling to assign nodes on the local machine with a higher probability. (Angerd et al., 2020) proposes to sample and save a small subgraph from other local machines as an approximation of the original graph structure. Nonetheless, these methods are limited to a very shallow GNN structure and suffer from significant performance degradation when the original graph is dense. (2) Input and hidden feature communication-based methods: These methods propose to communicate hidden features in addition to the input node features. Although these methods reduce the number of transferred bytes during each communication round (due to the smaller size of hidden embedding and less required nodes features), the number of communication rounds grows linearly as the number of layers, and are prone to more communication delay. To address these issues, in addition to optimal partitioning of the graph, (Zheng et al., 2020) proposes to use sparse embedding to reduce the number of bytes to communicate and (Tripathy et al., 2020) proposes several graph partitioning techniques to diminish the communication overhead.

heading: 2 BACKGROUND AND PROBLEM FORMULATION
text: In this section, we start by describing Graph Convolutional Network (GCN) and its training algorithm on a single machine, then formulate the problem of distributed GCN training. Note that we use GCN with mean aggregation for simplicity, however, our discussion is also applicable to other GNN architectures, such as SAGE (Hamilton et al., 2017), GAT (Velickovic et al., 2018), ResGCN (Li et al., 2019) and APPNP (Klicpera et al., 2019).
Training GCN on a single machine. Here, we consider the semi-supervised node classification in an undirected graph G(V, E) with N = |V| nodes and |E| edges. Each node vi ∈ V is associated with a pair (xi,yi), where xi ∈ Rd is the input feature vector, yi ∈ R|C| is the ground truth label, and C is the candidate labels in the multi-class classifications. Besides, let X = [x1, . . . ,xN ] ∈ RN×d denote the input node feature matrix. Our goal is to find a set of parameters θ = {W(`)}L`=1 by minimizing the empirical loss L(θ) over all nodes in the training set, i.e.,
L(θ) = 1 N ∑ i∈V φ(h (L) i ,yi), h (`) i = σ ( 1 |N (vi)| ∑ j∈N (vi) h (`−1) j W (`) ) , (1)
where φ(·, ·) is the loss function (e.g., cross entropy loss), σ(·) is the activation function (e.g., ReLU), and N (vi) is the neighborhood of node vi. In practice, we can update the model parameters by the stochastic gradient computed on a sampled mini-batch (using full-neighbors) by
∇̃L(θ, ξ) = 1 B ∑ i∈ξ ∇φ(h(L)i ,yi), (2)
where ξ denotes an i.i.d. sampled mini-batch of size B and we have E[∇̃L(θ, ξ)] = ∇L(θ). Distributed GCN training with periodic averaging. In this paper, we consider the distributed learning setting with P local machines and a single parameter server. The original input graph G is partitioned into P subgraphs, where Gp(Vp, Ep) denotes the subgraph on the p-th local machine with Np = |Vp| nodes, and Xp ∈ RNp×d as the input feature of all nodes in Vp located on the p-th machine. Then, the full-batch local gradient∇Llocalp (θp) is computed as ∇Llocalp (θp) = 1
Np
∑ i∈Vp ∇φ(h(L)i ,yi), h (`) i = σ ( 1 |Np(vi)| ∑ j∈Np(vi) h (`−1) j W (`) p ) , (3)
where θp = {W(`)p }L`=1 is the model parameters on the p-th local machine, Np(vi) = {vj |(vi, vj) ∈ Ep} is the local neighbors of node vi on the p-th local machine. When the graph is large, the
Algorithm 1 Distributed GCN training with “Parallel SGD with Periodic Averaging” Input: Global parameters θ̄0, local parameters θ0p = θ̄0, time-step t = 0, learning rate η. 1: for r← 1 to R do 2: for p← 1 to P do in parallel . Parallel training on local machines 3: Local machine p receives the global parameters θtp ← θ̄t. . Communication 4: for k← 1 to K do 5: t← t+ 1. 6: Local machine p constructs the mini-batch ξtp with neighbor sampling. 7: Local machine p computes the stochastic gradients ∇̃Llocalp (θtp, ξtp). 8: Local machine p updates the local parameter by θt+1p = θtp− η∇̃Llocalp (θtp, ξtp). 9: end for 10: Local machine p sends the local parameters θt+1p to the server. . Communication 11: end for 12: Server updates the global parameters by parameter averaging θ̄t+1 = 1
P ∑P p=1 θ t+1 p .
13: end for Output: Server returns trained GCN model with mint E[‖∇L(θ̄t)‖2].
computational complexity of forward and backward propagation could be very high. One practical solution is to compute the stochastic gradient on a sampled mini-batch with neighbor sampling, i.e.,
∇̃Llocalp (θp, ξp) = 1
Bp
∑ i∈ξp ∇φ(h̃(L)i ,yi), h̃ (`) i = σ ( 1 |Ñp(vi)| ∑ j∈Ñp(vi) h̃ (`−1) j W (`) p ) , (4)
where ξp is an i.i.d. sampled mini-batch of Bp nodes, Ñp(vi) ⊂ N (vi) is the sampled neighbors. An illustration of distributed GCN training with Parallel SGD with Periodic Averaging (PSGD-PA) is summarized in Algorithm 1. Before training, the server maintains a global model θ̄0 and each local machine keeps a local copy of the same model θ0p. During training, the local machine first updates the local model θtp using the stochastic gradient ∇̃Llocalp (θtp, ξtp) computed by Eq. 4 for K iterations (line 8), then sends the local model θtp to the server (line 10). At each communication step, the server collects and averages the model parameters from the local machines (line 12) and send the averaged model θt+1p back to each local machine.
Limitations. Although PSGD-PA can significantly reduce the communication overhead by transferring the locally trained models instead of node feature/embeddings (refer to Figure 2 (b)), it suffers from performance degeneration due to ignorance of the cut-edges (refer to Figure 2 (a)). In the next section, we introduce a communication-efficient algorithm LLCG that does not suffer from this issue, and can achieve almost the same performance as training the model on a single machine.

heading: 3 PROPOSED ALGORITHM: LEARN LOCALLY CORRECT GLOBALLY
text: In this section, we describe Learn Locally, Correct Globally (LLCG) for distributed GNN training. LLCG includes two main phases, local training with periodic model averaging and global server correction, to help reduce both the number of required communications and size of transferred data, without compromising the predictive accuracy. We summarize the details of LLCG in Algorithm 2.

heading: 3.1 LOCAL TRAINING WITH PERIODIC MODEL AVERAGING
text: At the beginning of a local epoch, each local machine receives the latest global model parameters from the server (line 3). Next, each local machine runs Kρr iterations to update the local model (line 4 to 9), where K and ρ are the hyper-parameters that control the local epoch size. Note that instead of using a fixed local epoch size as Algorithm 1, we choose to use exponentially increasing local epoch size in LLCG with ρ > 1. The reasons are as follows.
At the beginning of the training phase, all local models θtp are far from the optimal solution and will receive a gradient ∇̃Llocalp (θtp, ξtp) computed by Eq. 4. Using a smaller local update step at the early stage guarantees each local model does not diverge too much from each other before the model averaging step at the server side (line 12). However, towards the end of the training, all local models θtp will receive relatively smaller gradient ∇̃Llocalp (θtp, ξtp), such that we can chose a larger local epoch size to reduce the number of communications, without worrying about the divergence of
Algorithm 2 Distributed GCN training by “Learn Locally, Correct Globally” Input: Global parameters θ̄0, local parameters θ0p, time-step t = 0, local step size hyper-
parameters K, ρ, and learning rate γ, η 1: for r← 1 to R do 2: for p← 1 to P do in parallel . Parallel training on local machine 3: Local machine p receives the global parameters θtp ← θ̄t . Communication 4: for k← 1 to Kρr do 5: t← t+ 1 6: Local machine p constructs the mini-batch ξtp with neighbor sampling 7: Local machine p computes stochastic gradients ∇̃Llocalp (θtp, ξtp) 8: Local machine p updates model parameter by θt+1p = θtp − η∇̃Llocalp (θtp, ξtp) 9: end for
10: Local machine p sends the local parameters θt+1p to the server . Communication 11: end for 12: Server updates the global parameters using parameter averaging θ̄t+1 = 1
P ∑P p=1 θ t+1 p
13: for s← 1 to S do . Server Correction 14: t← t+ 1 15: Server constructs a mini-batch ξt with full-neighbors 16: Server computes the stochastic gradient ∇̃L(θ̄t, ξt) 17: Server updates the global parameters by θ̄t+1 = θ̄t − γ∇̃L(θ̄t, ξt) 18: end for 19: end for Output: Server return GCN model with trained mint E[‖∇L(θ̄t)‖2]
local models. By doing so, after total number of T = ∑R r=1Kρ
r iterations, LLCG only requires R = logρ T K rounds of communications. Therefore, compared to the fully-synchronous method, we can significantly reduce the total number of communications from O(T ) to O(logρ TK ).

heading: 3.2 GLOBAL SERVER CORRECTION
text: The design of the global server correction is to ensure that the trained model not only learns from the data on each local machine, but also learns the global structure of the graph, thus reducing the information loss caused by graph partitioning and avoiding cut-edges. Before the correction, the server receives the locally trained models from all local machines (line 10) and applies model parameter averaging (line 12). Next, S server correction steps are applied on top of the averaged model (line 13 to 18). During the correction, the server first constructs a mini-batch ξt using fullneighbors1(line 15), compute the stochastic gradient ∇̃L(θ̄t, ξt) on the constructed mini-batch by Eq. 2 (line 16) and update the averaged model θ̄t for S iterations (line 17). The number of correction steps S 2 depends on the heterogeneity among the subgraphs on each local machine: the more heterogeneous the subgraphs are, the more correction steps are required to better refine the averaged model and reduce the divergence across the local models. Note that, the heterogeneity is minimized when employing GGS (Figure 2) with the local machines having access to the full graph, as a result. However, GGS requires sampling from the global graph and communication at every iteration, which results in additional overhead and lower efficiency. Instead, in LLCG we are trading computation on the server for the costly feature communication, and only requires periodic communication.

heading: 4 THEORETICAL ANALYSIS
text: In this section, we provide the convergence analysis on the distributed training of GCN under two different settings, i.e., with and without server correction. In the following, we first introduce the notations and assumptions for the analysis (Section 4.1). Then, we show that periodic averaging of local machine models alone and ignoring the global graph structure will suffer from an irreducible residual error (Section 4.2). Finally, we show that this residual error can be eliminated by running server correction steps after each periodic averaging step on the server (Section 4.3).
1Note that using full neighbors is required for the server correction but not the local machines 2In practice, we found S = 1 or S = 2 works well on most datasets.
Central
∇ℒ 𝜃 ≠ 0.5× ∇ℒ!"#$%" 𝜃 + ∇ℒ&"#$%"(𝜃)
𝔼 -∇ℒ! "#$%"(𝜃) ≠ ∇ℒ! "#$%" 𝜃 𝔼 -∇ℒ&"#$%"(𝜃) ≠ ∇ℒ&"#$%"(𝜃)
\mathcal{N}(v_i)

heading: 4.1 NOTATIONS AND ASSUMPTIONS
text: Let us first recall the notations defined in Section 2, where L(θ) denotes the global objective function computed using the all node features X and the original graph G, Lp(θ) denotes the local objective function computed using the local node features Xp and local graph Gp, θtp denotes the model parameters on the p-th local machine at the t-th step, and θ̄t = 1P ∑P p=1 θ t p denotes the virtual averaged model at the t-th step. In the non-convex optimization, our goal is to show the expected gradient of the global objective on the virtual averaged model parameters E[‖∇L(θ̄t)‖2] decreases as the number of local machines P and the number of training steps T increase. Besides, we introduce ∇Lfullp (θ) as the gradient computed on the p-th local machine but have access the full node features X and the original graph structure G as ∇Lfullp (θ) = 1
|Vp| ∑ i∈Vp ∇φ(h(L)i , yi), h (`) i = σ ( 1 |N (vi)| ∑ j∈N (vi) h (`−1) j W (`) p ) . (5)
Please refer to Figure 3 for an illustration of different gradient computations. Besides, we introduce local-global gradient discrepancy as κ2 = κ2A + κ 2 X, where κ 2 A = maxp∈[P ]{‖∇Llocalp (θ) − ∇Lfullp (θ)‖2} is the maximum difference between the gradient computed on the local machine with and without having access to the global graph structure, which is mainly due to fact that the local machines are oblivious to the full graph information; and κ2X = maxp∈[P ]{‖∇Lfullp (θ)−∇L(θ)‖2} is the maximum difference between the gradient computed using the local node and all nodes, which is mainly due to the heterogeneity of the node features on each local machine, and we have κ2X = 0 if the nodes are i.i.d. sampled to each local machine. Notice that local-global gradient discrepancy κ2 plays an important role in our theoretical results.
For the convergence analysis, we make the following standard assumptions.
Assumption 1 The stochastic gradient on the p-th local machine (with neighbor sampling) has stochastic gradient variance bounded by σ2var and stochastic gradient bias bounded by σ 2 bias, i.e., E[‖∇̃Llocalp (θ; ξ)− E[∇̃Llocalp (θ; ξ)]‖2] ≤ σ2var, E[‖E[∇̃Llocalp (θ; ξ)]−∇Llocalp (θ)‖2] ≤ σ2bias.
Assumption 2 The stochastic gradient for global server correction (with full neighbors) has stochastic gradient variance bounded by σ2global, i.e., E[‖∇̃Lfullp (θ; ξ)−∇Lfullp (θ)]‖2] ≤ σ2global. The existence of stochastic gradient bias and variance in sampling-based GNN training have been studied in (Cong et al., 2020; 2021), where (Cong et al., 2021) further quantify the stochastic gradient bias and variance as a function of the number of GCN layers. In particular, they show that the existence of σ2bias is due to neighbor sampling and non-linear activation, and we have σ 2 bias = 0 if all neighbors are used or the non-linear activation is removed. The existence of σ2var is because we are sampling mini-batches to compute the stochastic gradient on each local machine during training. As the mini-batch size increases, σ2var will be decreasing, and we have σ 2 var = 0 when using full-batch.

heading: 4.2 DISTRIBUTED GNN VIA PARAMETER AVERAGING
text: In the following, we provide the first convergence analysis on distributed training of GCN. We show that solely periodic averaging of the local machine models and ignoring the global graph structure suffers from an upper bound that is irreducible with the number of training steps. Comparing to the traditional distributed training (e.g., distributed training Convolutional Neural Network for image classification (Dean et al., 2012; Li et al., 2020b)), the key challenges in the distributed GCN training is the two different types of gradient bias: (1) The expectation of the local full-batch gradient is
a biased estimation of the global full-batch gradient, i.e., 1P ∑P p=1∇Llocalp (θ) 6= ∇L(θ). This is because each local machine does not have access to the original input graph and full node feature matrix. Note that the aforementioned equivalence is important for the classifcal distributed training analysis Dean et al. (2012); Yu et al. (2019). (2) The expectation of the local stochastic gradient is a biased estimation of the local full-batch gradient i.e., E[∇̃Llocalp (θ, ξp)] 6= ∇Llocalp (θ). This is because the stochastic gradient on each local machine is computed by using neighbor sampling, which has been studied in (Cong et al., 2021).
Theorem 1 (Distributed GCN via Parameter Averaging) Consider applying model averaging for GNN training under Assumption 1 and 2. If we choose learning rate η =
√ P√ T
and the local step size K ≤ √ 2T 1/4
8LP 3/4 , then for any T ≥ L2P steps of gradient updates we have
1
T ∑T−1 t=0
E[‖∇L(θ̄t)‖2] = O (
1√ PT
) +O(κ2 + σ2bias).
Theorem 1 implies that, by carefully choosing the learning rate η and the local step sizeK, the gradient norm computed on the virtual averaged model is bounded byO(1/ √ PT ) afterR = T/K = O(P 3/4
T 3/4 )
communication rounds, but suffers from an irreducible residual error upper bound O(κ2 + σ2bias). In the next section, we show that this residual error can be eliminated by applying server correction.

heading: 4.3 DISTRIBUTED GCN VIA SERVER CORRECTION
text: Before proceeding to our result, in order to simplify the presentation, let us first define the notation Grglobal = mint∈Tglobal(r) E[‖∇L(θ̄t)‖2] and Grlocal = mint∈Tlocal(r) E [∥∥ 1 P ∑P p=1∇Llocalp (θtp) ∥∥2] as the minimum gradient computed at the r-th round global and local step, where Tglobal(r) and Tlocal(r) are the number of iteration run after the r-th communication round on server and local machine, respectively. Please refer to Eq. 42 in Appendix C.2 for a formal definition.
Theorem 2 Consider applying model averaging for GCN training under Assumption 1 and 2. If we choose learning rate γ = η =
√ P√ T , the local step size K, ρ such that ∑R r=1K 2ρ2r ≤ RT 1/2 32L2P 3/2 , and
server correction step size S = maxr∈[R] ( κ2+2σ2bias 1−L( √ P/T ) −Grlocal ) Kρr Grlocal , then for any T ≥ L2P steps of gradient updates we have: 1T ∑T t=1 E[‖∇L(θ̄t)‖2] = O ( 1√ PT ) .
Theorem 2 implies that, by carefully choosing the learning rates γ and η, the local step size hyperparameters K, ρ, and the number of global correction steps S, after T steps (R rounds of communication), employing parameter averaging with Global Server Correction, we have the norm of gradient bounded by O(1/ √ PT ), without suffering the residual error that exists in the naive parameter averaging (in Theorem 1). Besides, the server correction step size is proportional to the scale of κ2 and local stochastic gradient bias σ2bias. The larger κ
2 and σ2bias, the more corrections are required to eliminate the residual error. However, in practice, we observe that a very small number of correction steps (e.g., S = 1) performs well, which minimizes the computation overhead on the server.

heading: 5 EXPERIMENTS
text: Real-world simulation. In a real-world distributed setting, the server and local machines are located on different machines, connected through the network (Li et al., 2020a). However, for our experiments, we only have access to a single machine with multiple GPUs. As a result, we simulate a real-world distributed learning scenario, such that each GPU is responsible for the computation of two local machines (8 in total) and the CPU acts as the server. For these reasons, in our evaluations, we opted to report the communication size and number of communication rounds, instead of the wall-clock time, which can show the benefit of distributed training. We argue that these are acceptable measures in real-world scenarios as well since the two main factors in distributed training are initializing connection overhead and bandwidth (Tripathy et al., 2020). Baselines. To illustrate the effectiveness of LLCG, we setup two general synchronized distributed training techniques as the our baseline methods, namely “Parallel SGD with Parameter Averaging” (PSGD-PA) and “Global Graph Sampling” (GGS), as introduced in Figure 2, where the cut-edges in PSGD-PA are ignored and only the model parameters are transferred, but the cut-edges in GGS are considered and the node features of the cut-edges are transferred to the corresponding machine.
Note that we choose GGS as a reasonable representation for most existing proposals (Md et al., 2021; Zheng et al., 2020; Tripathy et al., 2020) for distributed GNN training, since these methods have very close communication cost and also require a large cluster of machines to truly show their performance improvement. We also use PSGD-PA as a lower bound for communication size, which is widely used in traditional distributed training and similar to the one used in (Angerd et al., 2020; Jiang & Rumi, 2021). However, we did not specifically include these methods in our results since we could not reproduce their results in our settings. Please refer to Appendix A for a detailed description of implementation, hardware specification and link to our source code. Datasets and evaluation metric. We compare LLCG and other baselines on real-world semisupervised node classification datasets, detail of which are summarized in Table 2 in the Appendix. The input graphs are splitted into multiple subgraphs using METIS before training, then the same set of subgraphs are used for all baselines. For training, we use neighborhood sampling (Hamilton et al., 2017) with 10 neighbors sampled per node and ρ = 1.1 for LLCG. For a fair comparison, we chose the base local update step K such that LLCG has the same number of local update steps as PSGD-PA. During evaluation, we use full-batch without sampling, and report the performance on the full graph using AUC ROC and F1 Micro as the evaluation metric. Unless otherwise stated, we conduct each experiment five times and report the mean and standard deviation.

heading: 5.1 PRIMARY RESULTS
text: In this section, we compare our proposed LLCG algorithm with baselines on four datasets. Due to space limitations we defer the detailed discussion on additional datasets to the Appendix A.4.
LLCG requires same number of communications. Figure 4 (a) through 4 (d) illustrate the validation accuracy per communication rounds on four different datasets. We run a fixed number of communication rounds and plot the global validation score (the validation score computed using the full-graph on the server) at the end of each communication step. For PSGD-PA and GGS, the score is calculated on the averaged model, whereas for LLCG the validation is calculated after the correction
0 20 40 Number of Communication
0.5
0.6
0.7
0.8
Gl ob
al Va
lid ati
on F
1 Sc
or e Arxiv
K = 1 K = 4
K = 16 K = 64
K = 128
Figure 5: Effect of local epoch size (K) 20 40 60 Number of Communications
0.0
0.5
Gl ob
al Va
l. F1
S co
re
Reddit
w/o sampling ,S=1 w/ sampling 20%,S=1 w/ sampling 5%, S=1 w/ sampling 5%, S=2
20 40 60 Number of Communications
0.5
0.6
0.7
Gl ob
al Va
l. F1
S co
re
Arxiv
w/o sampling ,S=1 w/ sampling 20%,S=1 w/ sampling 5%, S=1 w/ sampling 5%, S=4
Figure 6: Effect of sampling on local machine and number of correction steps on the server
step. It can be seen that PSGD-PA suffers from performance drop compared to other two methods, due to the residual error we discussed in Section 4, while both GGS and LLCG perform well and can achieve the expected accuracy. Note that the performance drop of PSGD-PA can vary across different datasets; in some cases such as Reddit, PSGD-PA can significantly hurt the accuracy, while on other datasets the gap is smaller. Nevertheless, LLCG can always close the gap between PSGD-PA and GGS with minimal overhead. LLCG convergences as fast as GGS. To represent the effect of communication on the real-time convergence, in Figure 4 (e) and 4 (f), we plot the global training loss (training loss computed on the full-graph on the server) after each communication round. Similar to the accuracy score, the training loss is also computed on the server averaged (and corrected, in case of LLCG) global model. These results clearly indicate that LLCG can improve the convergence over PSGD-PA, while it shows a similar convergence speed to GGS. LLCG exchanges data as little as PSGD-PA. Figure 4 (g) and 4 (h) show the relation between global validation accuracy with the average size (volume) of communication in bytes. As expected, this figure clearly shows the effectiveness of LLCG . On the one hand, LLCG has a similar amount of communication volume as PSGD-PA but can achieve a higher accuracy. On the other hand, LLCG requires significantly less amount of communication volume than GGS to achieve the same accuracy, which leads to slower training time in real world settings. LLCG works with various GNN models and aggregations. We evaluate four popular GNN models, used in recent graph learning literature: GCN Kipf & Welling (2017), SAGE Hamilton et al. (2017), GAT Velickovic et al. (2018) and APPNP Klicpera et al. (2019). In Table 1, we summarize the test score and average communication size (in MB) on different datasets for a fixed number of communication rounds. Note that we only include the results for the aggregation methods (GCN or SAGE) that have higher accuracy for the specific datasets, details of which can be found in Appendix A.2. As shown here, LLCG can consistently improve the test accuracy for all different models compared to PSGD-PA, while the communication size is significantly lower than GGS, since LLCG only needs to exchange the model parameters. Effect of local epoch size. Figure 5 compares the effect of various values of local epoch size K ∈ {1, 4, 16, 64, 128} for fixed ρ and S on the OGB-Arxiv dataset. When using fully synchronous with K = 1, the model suffers from very slow convergence and needs more communications. Further increasing the K to larger values can speed up the training; however, we found a diminishing return point for K > 128 in this dataset and extremely large K in general. Effect of sampling in local machines. In Figure 6, we report the validation scores per round of communication to compare the effect of neighborhood sampling at local machines. We can observe that when the neighborhood sampling size is reasonably large (i.e., 20%), the performance is very similar to full neighborhood training. However, reducing the neighbor sampling ratio to 5% could result in a larger local stochastic gradient bias σ2bias, which requires using more correction steps (S).

heading: 6 CONCLUDING REMARKS
text: In this paper, we propose a novel distributed algorithm for training Graph Neural Networks (GNNs). We theoretically analyze various GNN models and discover that, unlike the traditional deep neural networks, due to inherent data samples dependency in GNNs, naively applying periodic parameter averaging leads to a residual error and current solutions to this issue impose huge communication overheads. Instead, our proposal tackles these problems by applying correction on top of locally learned models, to infuse the global structure of the graph back into the network and avoid any costly communication. In addition, through extensive empirical analysis, we support our theoretical findings and demonstrate that LLCG can achieve high accuracy without additional communication costs.

heading: ACKNOWLEDGEMENTS
text: This work was supported in part by CRISP, one of six centers in JUMP, a Semiconductor Research Corporation (SRC) program sponsored by DARPA and NSF grants 1909004, 1714389, 1912495, 1629915, 1629129, 1763681, 2008398.

heading: REPRODUCIBILITY STATEMENT
text: We provide a GitHub repository in Appendix A including all code and scripts used in our experimental studies. This repository includes a README.md file, explaining how to install and prepare the code and required packages. Detailed instruction on how to use the partitioning scripts is provided for various datasets. In addition, we provide several configuration files (under scripts/configs) folder for different hyper-parameters on each individual dataset, and a general script (scripts/run-config.py) to run and reproduce the results with these configurations. Details of various models and parameters used in our evaluation studies can also be found in Appendix A.

heading: SUPPLEMENTARY MATERIAL

heading: LEARN LOCALLY, CORRECT GLOBALLY: A DISTRIBUTED ALGORITHM FOR TRAINING GRAPH NEURAL NETWORKS

heading: Appendix

heading: Table of Contents

heading: A Additional experimental setup and studies 14
text: A.1 Hardware specification and environment . . . . . . . . . . . . . . . . . . . . . 14 A.2 Dataset and models details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A.3 Minibatch selection for correction step . . . . . . . . . . . . . . . . . . . . . . 15 A.4 Effectiveness of correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 A.5 Effectiveness in large-scale settings . . . . . . . . . . . . . . . . . . . . . . . . 16

heading: B Proof of Theorem 1 17
text: B.1 Useful lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 B.2 Main proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 B.3 Discussion on the irreducible error . . . . . . . . . . . . . . . . . . . . . . . . 21 B.4 Proof of Lemma 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 B.5 Proof of Lemma 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 B.6 Proof of Lemma 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

heading: C Proof of Theorem 2 26
text: C.1 Useful lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 C.2 Main proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 C.3 Proof of Lemma 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

heading: A ADDITIONAL EXPERIMENTAL SETUP AND STUDIES
text: We provide detailed experimental setup and additional evaluations in this section. The link to the GitHub repository is available at https://github.com/MortezaRamezani/llcg/.

heading: A.1 HARDWARE SPECIFICATION AND ENVIRONMENT
text: For all our experiments we use a single server equipped with 4 NVIDIA QUADRO RTX 8000 GPUs with driver version 460.80, two Intel Xeon 6230 CPU and 768GB of main memory using Ubuntu 18.04 running kernel 5.4.0. For fair comparison between all methods, we developed a unified framework for performing the experiments, using Pytorch 1.7 compiled for CUDA 11.2. We also used Pytorch Geometric 1.7.0 and Pytorch Sparse 0.6.8 for various GNN and sparse operations. Please refer to README.md for detailed instruction on how to install and run the code, and also dataset specific configurations.

heading: A.2 DATASET AND MODELS DETAILS
text: In all experiments, we used ADAM optimizer, and the learning rate, alongside the local machine mini-batch size (BL) and server mini-batch size (BS) are included in config files under scripts/configs folder. We used three datasets Flickr, Reddit, Yelp from Zeng et al. (2020) and OGB-Proteins, OGB-Arxiv, OGB-Products from Hu et al. (2020). For each dataset, we choose the aggregation method that has been shown to achieve the highest accuracy, and reported the results in Table 1. In addition, we also evaluate two frequently used architectures, GAT and APPNP to show the flexibility of our proposed method to other models. Note that some of OGB datasets can use more complicated models Sun & Wu (2021; 2020), however for the sake of simplicity, we only used popular GNN operators which show competitive results in most cases. Next we formally introduce various layers and operators used in our experiments, followed by the base architecture for each dataset in Table 2.
• GCN (G): Originally introduced by Kipf et. al Kipf & Welling (2017), the representation for node vi is calculated using:
h (`) i = σ
( ∑
j∈N (vi)
αi,jh (`−1) j W
(`) ) , (6)
where αi,j = 1√ deg(vi) deg (vj) for symmetric normalized Laplacian and αi,j = 1deg(vi) for row normalized Laplacian.
• SAGE (S): First introduce in Hamilton et al. (2017), SAGE learn different weights for the node and its neighbors and the representation is calculated by:
h (`) i = σ ( h (`−1) i W (`) 1 +
1 |N (vi)| ∑
j∈N (vi)
h (`−1) j W (`) 2
) (7)
Note that different operator such as concatenation can be used instead of addition in the above formula. However in our experiments we used addition for all SAGE layers.
• Linear (L): Applies a linear transformation to the node features. In other word, the graph structured is completely ignored in this layer and the representation is computed as follows.
h (`) i = h (`−1) i W (`) (8)
• BatchNorm (B): Applies batch normalization according to Ioffe & Szegedy (2015) using the following where γ and β are learnable parameters and is added for numerical stability.
h (`) i =
h (`−1) i − E(h (`−1) i )√
Var(h (`−1) i ) +
∗ γ + β (9)
• GAT: The Graph Attention layer proposed by (Velickovic et al., 2018) as follows. h (`) i = σ ( ∑
j∈N (vi)
αi,jh (`−1) j W
(`) )
(10)
where αi,j is the attention between node i and j and calculated by:
αi,j = exp(LeakyRelu(a[h(`−1)i W (`)||h(`−1)j W(`)]))∑ k∈N (vi) exp(LeakyRelu(a[h (`−1) i W (`)||h(`−1)k W(`)])) (11)
where a is the learnable weight vector and || indicates concatenation.
• APPNP: This network is proposed by (Klicpera et al., 2019), which is a combination of graph-agnostic prediction and label propagation:
h (`) i = βh (0) i + (1− β)
∑
j∈N (vi)
αi,jh (`−1) j (12)
where αi,j is same as defined in GCN, β is the teleport probability and h(0) is the output of first linear layer.
Effect of sampling at correction. Recall that LLCG requires full-neighbors for the server correction step for the convergence analysis, however, we find that server correction with neighbor sampling also works well in practice. As shown in Figure 7 and Figure 8, although server correction with neighbor sampling can introduce some randomness at the beginning of the training phase, the final accuracy of training is very close to the server correction with full-neighbors.

heading: A.3 MINIBATCH SELECTION FOR CORRECTION STEP
text: Recall that in a server correction step, the mini-batch is selected by sampling from the entire graph, uniformly at random to estimate an unbiased stochastic gradient of the global loss function. However, one might suggest to include more cut-edges, which are missing from the local machines, in the correction minibatch instead of uniform random sampling, to improve the performance of the global model. To this end, we conduct an experiments on two datasets (Reddit and Arxiv), where we compare the random minibatch (default setting in LLCG ) and minibatch with higher number of cut-edges, by selecting the nodes on the ends of cut-edges and building the minibatch from there. As shown in Figure 9, including more cut-edges nodes in the mini-batch does not make significant improvement when comparing to selecting mini-batchs using uniform sampling. This potentially happens due to the fact that sampling more cut-edges make the gradient of server correction step biased, while we need an unbiased gradient of full graph for correction steps at server.

heading: A.4 EFFECTIVENESS OF CORRECTION
text: To illustrate the effectiveness of correction steps and relate it to the benefits of global graph structure, we conduct the same experiments on Yelp and OGB-Products with 700K and 2.4M nodes, respectively and report the results in Figure 10.
Yelp dataset. In case of Yelp in Figure 10 (a), as we can see the PSGD-PA and GGS are performing quite similarly. To further investigate this, we compared the validation accuracy of Yelp using GCN against MLP (i.e., without utilizing the graph structure in training), where the GCN layers replaced by Linear layer (i.e LLL instead of GGG architecture) and plot the validation accuracy per iterations in Figure 10 (b). As it can be deducted from this figure, for this dataset MLP is performing as good as GCN, which means this dataset does not depend on the global structure of the graph, further explaining why there is no performance gap between PSGD-PA and GGS. In other words, no server correction is necessary and we can use S = 0 in this case for LLCG.
OGB-Products. Similarly, in case of OGB-Products in Figure 10 (c), we cannot see any noticeable accuracy drop due to distributed training. However, unlike Yelp, this is mainly due to very small κ which is caused by two main factors: (1) very small ratio of training nodes for this dataset, only 8% of the nodes are used for training, and (2) very small number of cut-edges (less than 7%) after applying METIS on this dataset. It is also worth mentioning that due to shallow structure of the model required for this dataset, we barely need multi-hop neighbors and consequently we cannot see any noticeable performance drop.

heading: A.5 EFFECTIVENESS IN LARGE-SCALE SETTINGS
text: In this section, we further investigate the effect of large number of local machines and large graphs on the performance of proposed algorithm, under the simulated environment with up to 16 local machines on OGB-Products and OGB-MAG240M Hu et al. (2021) datasets for node classification. Besides,
we also compare against subgraph approximation distributed GNN training algorithm Angerd et al. (2020) using 10% extra storage overhead and fully synchronous distributed GNN training.3
Node classification tasks. For the OGB-Products dataset, we train a 3-layer GraphSAGE model with learning rate 0.003 with 50 rounds of communications , which has the same hyper-parameter configuration and model architecture as OGB-Products’s leaderboard implementation . For the OGB-MAG240M dataset, we train a 2-layer skip-connected GraphSAGE model with learning rate 0.001 with 400 rounds of communications, which has the same hyper-parameter selection and model architecture as OGB-MAG240M’s leaderboard implementation. As shown in Figure 11, we can observe that
• PSGD-PA suffers from performance degeneration and has a large gap to the full sync training accuracy. Subgraph approximation might alleviate the issue to some extent, but requires significant storage overhead. Our proposal LLCG can bridge the gap between PSGD-PA and full sync training.
• The pure computation time4 is negligible when ignoring the communication time but just consider the computation time. Besides, the server correction step does not introduce significant computation overhead.
It is worth mentioning that we do not observe the performance gap between PSGD-PA and fully sync distributed training in Figure 10 (c) because of different number of layers (2-layers) and limited number of local machine (8 local machines) are used in the previous experiments. However, as the number of layers and number of local machines increase, the effect of ignoring cut-edges and data heterogeneity is becoming more serious, which is an interesting observation and worth exploring as a potential future direction.

heading: B PROOF OF THEOREM 1
text: In this section, we first introduce the useful lemmas in Section B.1, then process our proof of theorem in Section B.2. In particular, we show that solely averaging the local machine models and ignoring the global graph structure will suffer from an irreducible residual error.

heading: B.1 USEFUL LEMMAS
text: The following lemma gives the upper bound for the norm of the stochastic gradient on each local machine as the norm of the expectation of local gradient and the stochastic gradient variance scaled by the number of local machines P .
Lemma 1 Let ∇̃Llocalp (θp; ξp) be stochastic gradients such that E [∥∥∥∇̃Llocalp (θp; ξp)− E[∇̃Llocalp (θp; ξp)] ∥∥∥ 2 ] ≤ σ2var, (13)
Then, we have
E [∥∥∥ 1 P P∑
p=1
∇̃Llocalp (θp; ξp) ∥∥∥ 2] ≤ σ 2 var
P + E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θp, ξtp) ] ∥∥∥ 2]
(14)
The following lemma provides an upper bound on the difference of each local gradient ∇Llocalp (θp) compared to 1P ∑P q=1∇Llocalq (θq), which is a function of:
• The deviation of each local model to the virtual averaged model, i.e., 1P ∑P i=1 ‖θp − θ̄‖2.
• The difference of gradient computed with and without having access to the full graph information, i.e., κ2.
3Suppose there are Np nodes on the p-th local machine, then subgraph of size 10%×Np is sampled and stored on the pth local machine. Notice that 10% storage overhead is the maximum overhead recommended by Angerd et al. (2020), which is expected to bring the best accuracy performance.
4We record the time using Python’s time.time() function.
Lemma 2 For any P machines θp, p ∈ {1, . . . , P}, if we define θ̄ = 1P ∑P i=1 θp and let κ > 0 such that ‖∇Llocalp (θ)−∇L(θ)‖2 ≤ κ2, we have 1
P
P∑
p=1
∥∥∥∇Llocalp (θp)− 1
P
P∑
q=1
∇Llocalq (θq) ∥∥∥ 2 ≤ 8L 2
P
P∑
i=1
‖θp − θ̄‖2 + 8κ2. (15)
The following lemma provides an upper bound on the deviation of each local model to virtual averaged model, which plays an important role in the previous lemma.
Lemma 3 For all P machines with parameters θp, p ∈ {1, . . . , P}, if we define θ̄ = 1P ∑P i=1 θp and let κ > 0 such that ‖Llocalp (θ)−∇L(θ)‖2 ≤ κ2, we have T−1∑
t=0
1
P
P∑
p=1
E[‖θ̄t − θtp‖2] ≤ 4η2Kσ2var 1− 16L2η2K2T + 4η2K2σ2bias + 16η 2K2κ2 1− 16L2η2K2 T. (16)

heading: B.2 MAIN PROOF OF THEOREM 1
text: Equipped with the above lemmas, we are now ready to present the proof of Theorem 1. From the smoothness assumption, we have
E[L(θ̄t+1)] ≤ E[L(θ̄t)] + E[〈∇L(θ̄t), θ̄t+1 − θ̄t〉] + L 2 E[‖θ̄t+1 − θ̄t‖2]. (17)
Recall that θ̄t is defined as θ̄t = 1P ∑P i=1 θ t p for any t. Then, according to the update rule
θ̄t+1 = θ̄t − η P
P∑
p=1
∇̃Llocalp (θtp; ξtp), (18)
we have
E[〈∇L(θ̄t), θ̄t+1 − θ̄t〉] = −ηE [〈 ∇L(θ̄t), 1
P
P∑
p=1
∇̃Llocalp (θtp; ξtp) 〉]
= −η 〈 ∇L(θ̄t), 1
P
P∑
p=1
E [ ∇̃Llocalp (θtp; ξtp) ] 〉 .
(19)
We can upper bound the right hand side of Eq. 19 by
− η 〈 ∇L(θ̄t), 1
P
P∑
p=1
E[∇̃Llocalp (θtp; ξtp)] 〉
= (a) −η 2
( ‖∇L(θ̄t)‖2 + ∥∥∥ 1 P P∑
p=1
E[∇̃Llocalp (θtp; ξtp)] ∥∥∥ 2 − ∥∥∥∇L(θ̄t)− 1
P
P∑
p=1
E[∇̃Llocalp (θtp; ξtp)] ∥∥∥ 2)
≤ (b) −η 2
( ‖∇L(θ̄t)‖2 + ∥∥∥ 1 P P∑
p=1
E[∇̃Llocalp (θtp; ξtp)] ∥∥∥ 2 − 2κ2 − 4σ2bias − 4L2
P
P∑
p=1
‖θ̄t − θtp‖2 ) ,
(20) where (a) is due to 2〈x,y〉 = ‖x‖2 + ‖y‖2 − ‖x− y‖2 and (b) is due to
∥∥∥∇L(θ̄t)− 1 P P∑
p=1
E[∇̃Llocalp (θtp; ξtp)] ∥∥∥ 2
= ∥∥∥∥∥ 1 P P∑
p=1
( ∇L(θ̄t)−∇Llocalp (θ̄t) +∇Llocalp (θ̄t)− E[∇̃Llocalp (θtp; ξtp)] )∥∥∥∥∥ 2
≤ 2 P
P∑
p=1
‖∇L(θ̄t)−∇Llocalp (θ̄t)‖2 + 4
P
P∑
p=1
‖∇Llocalp (θ̄t)−∇Llocalp (θtp)‖2
+ 4
P
P∑
p=1
∥∥∥E[∇̃Llocalp (θ̄t; ξtp)]−∇Llocalp (θtp) ∥∥∥ 2
≤ (c)
2κ2 + 4σ2bias + 4L2
P
P∑
p=1
‖θ̄t − θtp‖2,
(21)
where (c) follows from smoothness assumption and the definition of κ2, σ2bias in Theorem 1.
Combining Eq. 19 and Eq. 20 gives us E[〈∇L(θ̄t), θ̄t+1 − θ̄t〉]
≤ −η 2 E[‖∇L(θ̄t)‖2]− η 2 E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ η(κ2 + 2σ2bias) + 2ηL2
P
P∑
p=1
E[‖θ̄t − θtp‖2].
(22)
According to the update rule
θ̄t+1 = θ̄t − η P
P∑
p=1
∇̃Llocalp (θtp; ξtp), (23)
we have
E[‖θ̄t+1 − θ̄t‖2] = η2E [∥∥∥ 1 P P∑
p=1
∇̃Llocalp (θtp; ξtp) ∥∥∥ 2] . (24)
Substituting Eq. 22 and Eq. 24 to Eq. 17, we have
E[L(θ̄t+1)] ≤ E[L(θ̄t)]− η 2 E[‖∇L(θ̄t)‖2]− η 2 E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ η(κ2 + 2σ2bias) + 2ηL2
P
P∑
p=1
E[‖θ̄t − θtp‖2]
+ η2L 2 E [∥∥∥ 1 P P∑
p=1
∇̃Llocalp (θtp; ξtp) ∥∥∥ 2] .
(25)
Dividing both sides by η2 and rearranging the terms yields
E[‖∇L(θ̄t)‖2] ≤ 2 η
( E[L(θ̄t)]− E[L(θ̄t+1)] ) − E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ 2(κ2 + 2σ2bias) + 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2] + ηLE [∥∥∥ 1 P P∑
p=1
∇̃Llocalp (θtp; ξtp) ∥∥∥ 2]
≤ (a)
2
η
( E[L(θ̄t)]− E[L(θ̄t+1)] ) − E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξpt ) ] ∥∥∥ 2]
+ 2(κ2 + 2σ2bias) + 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2]
+ ηL ( E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ σ2var P
)
= 2
η
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + (ηL− 1)E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ 2(κ2 + 2σ2bias) + 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2] + ηLσ2var P ,
(26) where (a) is due to Lemma 1.
Summing over t ∈ {0, . . . , T − 1} and dividing both side by T , we get 1
T
T−1∑
t=0
E[‖∇L(θ̄t)‖2]
≤ 2 ηT
( L(θ̄0)− L(θ?) ) + ηL− 1 T T−1∑
t=0
E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ 2(κ2 + 2σ2bias) + 4L2
T
T−1∑
t=0
1
P
P∑
p=1
E[‖θ̄t − θtp‖2] + ηLσ2var P
≤ (a)
2
ηT
( L(θ̄0)− L(θ?) ) + ηL− 1 T T−1∑
t=0
E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θtp, ξtp) ] ∥∥∥ 2]
+ ηLσ2var P
+ 2(κ2 + 2σ2bias) + L 2 16η 2Kσ2var (1− 16η2L2K2) + L 2 16η 2K2σ2bias + 64η 2K2κ2 (1− 16η2L2K2)
≤ (b)
2
ηT
( L(θ̄0)− L(θ?) ) + 2(κ2 + 2σ2bias) +
ηLσ2var P
+ L2 16η2Kσ2var (1− 16η2L2K2) + L 2 16η 2K2σ2bias + 64η 2K2κ2 (1− 16η2L2K2)
≤ (c)
2
ηT
( L(θ̄0)− L(θ?) ) + 2(κ2 + 2σ2bias) +
ηLσ2var P
+ 32L2η2Kσ2var + 128L 2η2K2κ2 + 32L2η2K2σ2bias,
(27)
where (a) is due to Lemma 3, (b) is due to 0 < η < 1L , and (c) is due to K ≤ √ 2
8Lη is selected to satisfy 1− 16η2L2K2 ≥ 12 . If we choose η = √ P√ T and K ≤ √ 2T 1/4 8LP 3/4 = O(T 1/4 P 3/4 ), then for any T ≥ L2P we have
1
T
T−1∑
t=0
E[‖∇L(θ̄t)‖2] = O (
1√ PT
) +O(κ2 + σ2bias), (28)
with R = O(P 3/4T 3/4) rounds of communication.

heading: B.3 DISCUSSION ON THE IRREDUCIBLE ERROR
text: To better understand why the error in the upper bound might be irreducible (i.e., the second term in RHS of Theorem 1 is independent of T), let’s first recall the key sources that cause gradient diversity∥∥∥∇L(θ̄t)− 1P ∑P p=1 E[∇̃Llocalp (θtp; ξtp)] ∥∥∥ 2 (in Eq. 21), in our setting:
• The κ2 term, which is the upper bound of 1P ∑P p=1 ‖∇L(θ̄t)−∇Llocalp (θ̄t)‖2 and is caused
by ignoring the cut-edges for local gradient computation;
• The σ2bias term, which is the upper bound of 1 P ∑P p=1 ∥∥∥E[∇̃Llocalp (θ̄t; ξtp)]−∇Llocalp (θtp) ∥∥∥ 2
and is caused by using neighbor sampling;
• The model divergence term 1P ∑P p=1 ‖θ̄t − θtp‖2 which is caused by the difference between
the model parameters and the virtual average model due to periodic averaging. Notice that the model divergence term also exists in distributed learning regardless of whether the dataset is graph or not and caused by infrequent synchronization.
Now, let’s have a closer look at the model divergence term. As shown in Lemma 3 the model divergence term 1P ∑P p=1 ‖θ̄t − θtp‖2 is further caused by three factors:
• The σ2var term (in Eq. 64), which is the upper bound of mini-batch sampling variance 1 P ∑P p=1 E [∥∥∥ ∑t−1 τ=t0 ( E[∇̃Llocalp (θτp ; ξτp )]−∇Llocalp (θτp ) )∥∥∥ 2]
• The κ2 term, which is the upper bound of 1P ∑P p=1 ‖∇L(θ̄t)−∇Llocalp (θ̄t)‖2 and is caused
by ignoring the cut-edges for local gradient computation;
• The σ2bias term, which is the upper bound of 1 P ∑P p=1 ∥∥∥E[∇̃Llocalp (θ̄t; ξtp)]−∇Llocalp (θtp) ∥∥∥ 2
and is caused by using neighbor sampling.
Fortunately, the model divergence term can be controlled by the number of local gradient update steps and learning rate, which is reducing with respect to the number of total training steps T and the number of local machine P , and it leads to the first term in our upper bound in Theorem 1, i.e., O(σ 2 var+σ 2 bias+κ 2
√ PT
) = O( 1√ PT ). However, unfortunately, κ2 and σ2bias in gradient diversity are isolated from the model diversity part, therefore are irreducible and results in O(σ2bias + κ2) in the second term in our upper bound in Theorem 1. Please refer to Figure 12 for an illustration.
Intuitively, these theoretical results make sense. During training, we are minimizing the loss without cut-edges 1P ∑P p=1 Llocal(θ), which is a different objective to the loss defined on the full-graph by taking the cut-edges into consideration L(θ). Therefore, solely by adding the number of training iterations, we cannot guarantee a small gradient of L(θ) by minimizing 1P ∑P p=1 Llocal(θ).

heading: B.4 PROOF OF LEMMA 1
text: By the definition of Llocalp (θp; ξp), we have
E [∥∥∥ 1 P P∑
p=1
∇̃Llocalp (θp; ξp) ∥∥∥ 2]
= E [∥∥∥ 1 P P∑
p=1
( ∇̃Llocalp (θp; ξp)− E [ ∇̃Llocalp (θp, ξtp) ]) ∥∥∥ 2 ] + E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θp, ξtp) ] ∥∥∥ 2]
≤ (a)
1
P 2
P∑
p=1
E [∥∥∥∇̃Llocalp (θp; ξp)− E [ ∇̃Llocalp (θp; ξp) ] ∥∥∥ 2 ] + E   ∥∥∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θp, ξtp) ]∥∥∥∥∥ 2  
≤ (b) σ2var P + E [∥∥∥ 1 P P∑
p=1
E [ ∇̃Llocalp (θp, ξtp) ] ∥∥∥ 2]
(29) where (a) is due to fact that each ∇̃Llocalp (θp; ξp)−E[∇̃Llocalp (θp; ξp)] is independent random vectors with zero mean and (b) is from Assumption 1.

heading: B.5 PROOF OF LEMMA 2
text: By the definition of local gradient∇Llocalp (θ), we have 1
P
P∑
p=1
∥∥∥∇Llocalp (θp)− 1
P
P∑
q=1
∇Llocalq (θq) ∥∥∥ 2
= 1
P
P∑
p=1
∥∥∥∇Llocalp (θp)−∇Llocalp (θ̄) +∇Llocalp (θ̄)−∇L(θ̄)
+∇L(θ̄)− 1 P
P∑
q=1
∇L(θq) + 1
P
P∑
q=1
∇L(θq)− 1
P
P∑
q=1
∇Llocalq (θq) ∥∥∥ 2
≤ (a)
1
P
P∑
p=1
4‖∇Llocalp (θp)−∇Llocalp (θ̄)‖2 + 1
P
P∑
p=1
4‖∇Llocalp (θ̄)−∇L(θ̄)‖2
+ 1
P
P∑
p=1
4 ∥∥∥∇L(θ̄)− 1
P
P∑
q=1
∇L(θq) ∥∥∥ 2 + 1
P
P∑
p=1
4 ∥∥∥ 1 P P∑
q=1
∇L(θq)− 1
P
P∑
q=1
∇Llocalq (θq) ∥∥∥ 2
≤ (b)
8L
P
P∑
p=1
‖θp − θ̄‖2 + 8κ2,
(30) where (a) is due to ‖∑ni=1 xi‖2 ≤ ∑n i=1 n‖xi‖ and (b) is due to the definition of κ.

heading: B.6 PROOF OF LEMMA 3
text: When (t mod K) = 0 we have θ̄t = θtp. When (t mod K) 6= 0 and t ≥ 1, let t0 < t be the largest iteration index that t0 mod K = 0. For any τ ∈ {t0 + 1, . . . , t}, we have
θτp − θτ−1p = −η∇̃Llocalp (θτ−1p ; ξτ−1p ). (31)
Summing over τ ∈ {t0 + 1, . . . , t}, we have
θtp = θ t0 p − η
t−1∑
τ=t0
∇̃Llocalp (θτp ; ξτp )
= (a) θ̄t0 − η
t−1∑
τ=t0
∇̃Llocalp (θτp ; ξτp ), (32)
where (a) is due to θ̄t = θtp when (t mod K) = 0.
Similarly, we have
θ̄t = θ̄t0 − η t−1∑
τ=t0
1
P
P∑
p=1
∇̃Llocalp (θτp ; ξτp ). (33)
Combining Eq. 32 and Eq. 33, we have 1
P
P∑
p=1
E[‖θ̄t − θtp‖2]
= η2
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇̃Llocalp (θτp ; ξτp )− 1
P
P∑
q=1
∇̃Llocalq (θτq ; ξτq ) )∥∥∥ 2]
≤ (a)
2η2 1
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
(( ∇̃Llocalp (θτp ; ξτp )−∇Llocalp (θτp ) ) − 1 P P∑
q=1
( ∇̃Llocalq (θτq ; ξτq )−∇Llocalq (θτq ) ))∥∥∥ 2]
︸ ︷︷ ︸ (A)
+ 2η2 1
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇Llocalp (θτp )− 1
P
P∑
q=1
∇Llocalq (θτq ) )∥∥∥ 2]
︸ ︷︷ ︸ (B)
,
(34) where (a) is due to addition and subtraction of ∇̃Llocalp (θτp ; ξτp ) and∇Llocalp (θτp ) and using ‖x+y‖2 ≤ ‖x‖2 + ‖y‖2. We can upper bound (A) in Eq. 63 by
(A) ≤ (a)
1
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇̃Llocalp (θτp ; ξτp )−∇Llocalp (θτp ) )∥∥∥ 2]
≤ 2 P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇̃Llocalp (θτp ; ξτp )− E[∇̃Llocalp (θτp ; ξτp )] )∥∥∥ 2]
+ 2
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( E[∇̃Llocalp (θτp ; ξτp )]−∇Llocalp (θτp ) )∥∥∥ 2]
≤ 2Kσ2var + 2K2σ2bias,
(35)
where (a) is due to 1
n
n∑
i=1
∥∥∥xi − 1
n
n∑
j=1
xj ∥∥∥ 2 = 1
n
n∑
i=1
‖xi‖2 − ∥∥∥ 1 n n∑
i=1
xi ∥∥∥ 2
≤ 1 n
n∑
i=1
‖xi‖2. (36)
We can further bound (B) in Eq. 63 by
(B) ≤ 1 P
P∑
p=1
(t− t0) t−1∑
τ=t0
E [∥∥∥∇Llocalp (θτp )− 1
P
P∑
q=1
∇Llocalq (θτq ) ∥∥∥ 2]
≤ K t−1∑
τ=t0
1
P
P∑
p=1
E [∥∥∥∇Llocalp (θτp )− 1
P
P∑
q=1
∇Llocalq (θτq ) ∥∥∥ 2]
≤ (a) K
t−1∑
τ=t0
(8L2 P P∑
i=1
‖θτp − θ̄τ‖2 + 8κ2 )
≤ 8L 2K
P
t−1∑
τ=t0
P∑
i=1
‖θτp − θ̄τ‖2 + 8K2κ2,
(37)
where (a) is due to Lemma 2 and the definition of κ.
By plugging Eq. 64 and Eq. 66 into Eq. 63, we have 1
P
P∑
p=1
E[‖θ̄t−θtp‖2] ≤ 4η2Kσ2var +4η2K2σ2bias +16η2L2K t−1∑
τ=t0
1
P
P∑
i=1
‖θτp − θ̄τ‖2 +16η2K2κ2.
(38)
By summing over t ∈ {0, . . . , T}, we have T−1∑
t=0
1
P
P∑
p=1 E[‖θ̄t − θtp‖2] ︸ ︷︷ ︸
(A)
≤ 4η2Kσ2varT + 4η2K2σ2biasT
+ 16η2L2K2 T−1∑
t=0
1
P
P∑
i=1
‖θtp − θ̄t‖2
︸ ︷︷ ︸ (A)
+16η2K2κ2T.
(39)
Collecting common terms and dividing both sides by (1− 16η2L2K2) gives (A) ≤ 4η 2Kσ2T
(1− 16η2L2K2) + 4η2K2σ2biasT + 16η 2K2κ2T (1− 16η2L2K2) . (40)

heading: C PROOF OF THEOREM 2
text: In the following, we first introduce the useful lemmas in Section C.1, then process our proof of theorem in Section C.2. In particular, we show that this residual error (which is cause by ignoring the cut-edges) can eliminated by running server correction steps after the parameter averaging on the server.

heading: C.1 USEFUL LEMMA
text: The following lemma provides an upper bound on the deviation of each local model to virtual averaged model, which is important to upper bound on the difference of each local gradient ∇Llocalp (θp) compared to 1P ∑P p=1∇Lp(θ).
Lemma 4 For all P machines with θp, p ∈ {1, . . . , P}, if we define θ̄ = 1P ∑P p=1 θp and let κ > 0 such that ‖Llocalp (θ)−∇L(θ)‖2 ≤ κ2, we have ∑
t∈Tlocal(r)
1
P
P∑
p=1
E[‖θ̄t − θtp‖2] ≤ 4η2(σ2bias + σ 2 var)K 2ρ2r + 16η2κ2K2ρ2r
1− 16η2L2K2ρ2r . (41)

heading: C.2 MAIN PROOF OF THEOREM 2
text: Equipped with the above lemma and results from Appendix B.1 , we are now ready to present the proof of Theorem 2.
Let Tlocal(r) and Tserver(r) as the iteration indicies that a local machine and server run on, which is defined as
Tlocal(r) =   k +   r−1∑
j=1
Kρj  + S(r − 1) : k = 1, . . . ,Kρr   
Tserver(r) =   s+   r∑
j=1
Kρj  + S(r − 1) : s = 1, . . . , S    ,
(42)
and let define Tlocal = Tlocal(1) ∪ . . . ∪ Tlocal(R) and Tglobal = Tglobal(1) ∪ . . . ∪ Tglobal(R). By the smoothness assumption, we have
E[L(θ̄t+1)] ≤ E[L(θ̄t)] + E[〈∇L(θ̄t), θ̄t+1 − θ̄t〉] + L 2 E[‖θ̄t+1 − θ̄t‖2]. (43)
Let first consider t ∈ Tserver(r), with the following update θ̄t+1 = θ̄t − γ∇̃L(θ̄t; ξt)
= θ̄t − γ P
P∑
p=1
∇̃Lfullp (θ̄t; ξtp). (44)
Therefore, we have
E[〈∇L(θ̄t), θ̄t+1 − θ̄t〉] = −γE [〈 ∇L(θ̄t), 1
P
P∑
p=1
∇̃Lfullp (θ̄t; ξtp) 〉]
= (a) −γE[‖∇L(θ̄t)‖2],
(45)
where the equality is due to 1P ∑P p=1 E[∇̃Lfullp (θ̄t; ξtp)] = ∇L(θ̄t) since all neighbors are used for the server correction steps.
Besides, by taking the norm on the both side of Eq. 44, we have the following equality
E[‖θ̄t+1 − θ̄t‖2] = γ2E [∥∥∥ 1 P P∑
p=1
∇̃Lfullp (θ̄t; ξtp) ∥∥∥ 2]
≤ (a) γ2E [∥∥∥ 1 P P∑
p=1
∇Lfullp (θ̄t) ∥∥∥ 2] + γ2 σ2global P
= γ2E[‖∇L(θ̄t)‖2] + γ2 σ2global P .
(46)
where inequality (b) hold due to Lemma 1. Notice that the result in Lemma 1 holds for both the “local” and the “global” setting since 1P ∑P p=1 E[∇̃Lfullp (θ̄t; ξtp)] = ∇L(θ̄t).
Substituting Eq. 45, 46 into Eq. 43, we know that for t ∈ Tserver(r), we have
E[L(θ̄t+1)] ≤ E[L(θ̄t)] + (γ2L 2 − γ ) E[‖∇L(θ̄t)‖2] + γ 2L 2 σ2global P
= E[L(θ̄t)] + γ 2
( γL− 1 ) E[‖∇L(θ̄t)‖2] + γ 2L
2 σ2global P − γ 2 E[‖∇L(θ̄t)‖2].
(47)
Dividing both sides by γ2 and rearranging terms yields
E[‖∇L(θ̄t)‖2] ≤ 2 γ
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + (γL− 1)E[‖∇L(θ̄t)‖2] + γL σ2global P . (48)
Then, let consider the local update steps where t ∈ Tlocal(r). According to Eq. 25 in proof of Theorem 1, we have
E[‖∇L(θ̄t)‖2] ≤ 2 η
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + (ηL− 1)E [∥∥∥ 1 P P∑
p=1
∇Llocalp (θtp) ∥∥∥ 2]
+ 2(κ2 + 2σ2bias) + 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2] + ηLσ2var P .
(49)
Let T = (∑R r=1Kρ r )
+ SR and summing over t ∈ {1, . . . , T}, combining Eq. 48 and Eq. 49 we have
T∑
t=1
E[‖∇L(θ̄t)‖2]
=
R∑
r=1
∑
t∈Tlocal(r)
E[‖∇L(θ̄t)‖2] + R∑
r=1
∑
t∈Tglobal(r)
E[‖∇L(θ̄t)‖2]
≤ R∑
r=1
∑
t∈Tglobal(r)
[ 2 γ ( E[L(θ̄t)]− E[L(θ̄t+1)] ) + (γL− 1)E[‖∇L(θ̄t)‖2] + γL σ2global P ]
+
R∑
r=1
∑
t∈Tlocal(r)
[2 η ( E[L(θ̄t)]− E[L(θ̄t+1)] ) + (ηL− 1)E [∥∥∥ 1 P P∑
p=1
∇Llocalp (θtp) ∥∥∥ 2]
+ 2(κ2 + σ2bias) + 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2] + ηLσ2var P
] .
(50)
Rearrange the above equation, we have T∑
t=1
E[‖∇L(θ̄t)‖2]
≤ R∑
r=1
  ∑
t∈Tglobal(r)
2
γ
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + ∑
t∈Tlocal(r)
2
η
( E[L(θ̄t)]− E[L(θ̄t+1)] )  
+
R∑
r=1
  ∑
t∈Tglobal(r)
(γL− 1)E[‖∇L(θ̄t)‖2] + ∑
t∈Tlocal(r)
(ηL− 1)E [∥∥∥ 1 P P∑
p=1
∇Llocalp (θtp) ∥∥∥ 2
 
+
R∑
r=1
  ∑
t∈Tglobal(r)
γL σ2global P + ∑
t∈Tlocal(r)
ηLσ2var P
 + ( R∑
r=1
Kρr ) (κ2 + 2σ2var)
+
R∑
r=1
∑
t∈Tlocal(r)
( 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2] ) .
(51)
Let define
Grglobal = min t∈Tglobal(r) E[‖∇L(θ̄t)‖2], Grlocal = min t∈Tlocal(r) E [∥∥∥ 1 P P∑
p=1
∇Llocalp (θtp) ∥∥∥ 2 . (52)
Our goal is to select the size of Tlocal(r) and Tglobal(r) such that the following inequality holds (Kρr) (κ2 + 2σ2bias) ≤ (1− γL) ∑
t∈Tglobal(r)
E[‖∇L(θ̄t)‖2] + (1− ηL) ∑
t∈Tlocal(r)
E [∥∥∥ 1 P P∑
p=1
∇Llocalp (θtp) ∥∥∥ 2
≤ (1− γL)SGrglobal + (1− ηL)KρrGrlocal (53)
After rearranging it, we will have
S ≥ Kρ r Grglobal(1− γL) ( κ2 + 2σ2bias − (1− ηL)Grlocal ) . (54)
Suppose Eq. 54 holds, we have T∑
t=1
E[‖∇L(θ̄t)‖2]
≤ R∑
r=1
  ∑
t∈Tglobal(r)
2
γ
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + ∑
t∈Tlocal(r)
2
η
( E[L(θ̄t)]− E[L(θ̄t+1)] )  
+
R∑
r=1
  ∑
t∈Tglobal(r)
γL σ2global P + ∑
t∈Tlocal(r)
ηLσ2var P
 + R∑
r=1
∑
t∈Tlocal(r)
( 4L2
P
P∑
p=1
E[‖θ̄t − θtp‖2] )
≤ (a)
R∑
r=1
  ∑
t∈Tglobal(r)
2
γ
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + ∑
t∈Tlocal(r)
2
η
( E[L(θ̄t)]− E[L(θ̄t+1)] )  
+
R∑
r=1
  ∑
t∈Tglobal(r)
γL σ2global P + ∑
t∈Tlocal(r)
ηLσ2var P
 
+ 4L2 R∑
r=1
( 4η2(σ2bias + σ 2 var)K 2ρ2r
1− 16η2L2K2ρ2r + 16η2κ2K2ρ2r 1− 16η2L2K2ρ2r ) ,
(55) where (a) is due to Lemma 4.
By selecting Kρr such that 1− 16η2L2K2ρ2r ≥ 12 , we have Kρr ≤ √ 2 8Lη and
1
T
T∑
t=1
E[‖∇L(θ̄t)‖2]
= 1
T
R∑
r=1
  ∑
t∈Tglobal(r)
2
γ
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + ∑
t∈Tglobal(r)
γL σ2global P
 
+ 1
T
R∑
r=1
  ∑
t∈Tlocal(r)
2
η
( E[L(θ̄t)]− E[L(θ̄t+1)] ) + ∑
t∈Tlocal(r)
ηLσ2var P
 
+ 4η2L2
( 8(σ2var + σ 2 bias) + 32κ 2 )
T
R∑
r=1
K2ρ2r.
(56)
Notice that the condition Kρr ≤ √ 2
8Lη implies R∑
r=1
K2ρ2r = K2 1− ρ2r 1− ρ2 ≤
R
32L2η2 . (57)
To this end, by selecting η = γ = √ P√ T , ∑R r=1K 2ρ2r ≤ RT 1/2 32L2P 3/2
, and σ2 = max{σ2var, σ2global}, we have
1
T
T∑
t=1
E[‖∇L(θ̄t)‖2] ≤ 2√ PT
( L(θ̄0)−L(θ̄?) ) + Lσ2√ PT + 4L2√ PT R∑
r=1
(8(σ2bias+σ 2 var)+32κ 2), (58)
and
S ≥ (
κ2 + 2σ2bias 1− L( √ P/T ) −Grlocal
) Kρr
Grglobal . (59)

heading: C.3 PROOF OF LEMMA 4
text: Recall that when t ∈ Tglobal(r) we have θ̄t = θtp. For any t ∈ Tlocal(r), let t0 + 1 denote as the first indices of Tlocal(r) and |Tlocal(r)| = Kρr. Then, for any τ ∈ {t0 + 1, . . . , t}, we have θτp − θτ−1p = −η∇̃Llocalp (θτ−1p ; ξτ−1p ) (60) Summing over τ ∈ {t0 + 1, . . . , t}, we have
θtp = θ t0 p − η
t−1∑
τ=t0
∇̃Llocalp (θτp ; ξτp )
= (a) θ̄t0 − η
t−1∑
τ=t0
∇̃Llocalp (θτp ; ξτp ), (61)
where (a) is due to θ̄t = θtp when t ∈ Tglobal(r). Similarly, we have
θ̄t = θ̄t0 − η t−1∑
τ=t0
1
P
P∑
p=1
∇̃Llocalp (θτp ; ξτp ). (62)
By combining Eq. 61 and Eq. 62, we have 1
P
P∑
p=1
E[‖θ̄t − θtp‖2]
= η2
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇̃Llocalp (θτp ; ξτp )− 1
P
P∑
q=1
∇̃Llocalq (θτq ; ξτq ) )∥∥∥ 2]
≤ (a)
2η2 1
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
(( ∇̃Llocalp (θτp ; ξτp )−∇Llocalp (θτp ) ) − 1 P P∑
q=1
( ∇̃Llocalq (θτq ; ξτq )−∇Llocalq (θτq ) ))∥∥∥ 2]
︸ ︷︷ ︸ (A)
+ 2η2 1
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇Llocalp (θτp )− 1
P
P∑
q=1
∇Llocalq (θτq ) )∥∥∥ 2]
︸ ︷︷ ︸ (B)
,
(63) where (a) follows by adding and substracting ∇̃Llocalp (θτp ; ξτp ),∇Llocalp (θτp ) and using ‖x + y‖2 ≤ ‖x‖2 + ‖y‖2. We can upper bound (A) in Eq. 63 by
(A) ≤ (a)
1
P
P∑
p=1
E [∥∥∥ t−1∑
τ=t0
( ∇̃Llocalp (θτp ; ξτp )−∇Llocalp (θτp ) )∥∥∥ 2]
≤ 2 P
P∑
p=1
t−1∑
τ=t0
E [∥∥∥∇̃Llocalp (θτp ; ξτp )− E[∇̃Llocalp (θτp ; ξτp )] ∥∥∥ 2]
+ 2
P
P∑
p=1
t−1∑
τ=t0
E [∥∥∥E[∇̃Llocalp (θτp ; ξτp )]−∇Llocalp (θτp ) ∥∥∥ 2]
≤ 2Kρr(σ2var + σ2bias),
(64)
where (a) is due to 1
n
n∑
i=1
∥∥∥xi − 1
n
n∑
j=1
xj ∥∥∥ 2 = 1
n
n∑
i=1
‖xi‖2 − ∥∥∥ 1 n n∑
i=1
xi ∥∥∥ 2
≤ 1 n
n∑
i=1
‖xi‖2. (65)
We can further bound (B) in Eq. 63 by
(B) ≤ 1 P
P∑
p=1
(t− t0) t−1∑
τ=t0
E [∥∥∥∇Llocalp (θτp )− 1
P
P∑
q=1
∇Llocalq (θτq ) ∥∥∥ 2]
≤ Kρr t−1∑
τ=t0
1
P
P∑
p=1
E [∥∥∥∇Llocalp (θτp )− 1
P
P∑
q=1
∇Llocalq (θτq ) ∥∥∥ 2]
≤ (a)
Kρr t−1∑
τ=t0
(8L2 P P∑
p=1
‖θτp − θ̄τ‖2 + 8κ2 )
≤ 8L 2Kρr
P
t−1∑
τ=t0
P∑
p=1
‖θτp − θ̄τ‖2 + 8K2ρ2rκ2,
(66)
where (a) is due to Lemma 2 and the definition of κ.
By plugging Eq. 64 and Eq. 66 into Eq. 63, we have 1
P
P∑
p=1
E[‖θ̄t−θtp‖2] ≤ 4η2Kρr(σ2var+σ2bias)+16η2L2Kρr t−1∑
τ=t0
1
P
P∑
p=1
‖θτp−θ̄τ‖2+16η2K2ρ2rκ2.
(67)
Let us define Tlocal = Tlocal(1) ∪ . . . ∪ Tlocal(R). By summing over t ∈ Tlocal, we have R∑
r=1
∑
t∈Tlocal(r)
1
P
P∑
p=1 E[‖θ̄t − θtp‖2] ︸ ︷︷ ︸
(A)
≤ R∑
r=1
( 4η2(σ2bias + σ 2 var)K 2ρ2r + 16η2κ2K2ρ2r )
+
R∑
r=1
16η2L2K2ρ2r ∑
t∈Tlocal(r)
1
P
P∑
p=1
‖θtp − θ̄t‖2
︸ ︷︷ ︸ (A)
.
(68)
Rearranging the terms gives R∑
r=1
(1−16η2L2K2ρ2r) ∑
t∈Tlocal(r)
1
P
P∑
p=1 E[‖θ̄t − θtp‖2] ︸ ︷︷ ︸
(A)
≤ R∑
r=1
( 4η2((σ2bias + σ 2 var))K 2ρ2r + 16η2κ2K2ρ2r ) .
(69)
Therefore, we conclude that ∑
t∈Tlocal(r)
1
P
P∑
p=1
E[‖θ̄t − θtp‖2] ≤ 4η2((σ2bias + σ 2 var))K 2ρ2r + 16η2κ2K2ρ2r
1− 16η2L2K2ρ2r . (70)
which completes the proof.

Review summaries:
The paper considers the problem of distributed training for graph learning tasks, under a setting where data privacy is significant for each individual machine and communication to/from a central parameter server is expensive. To preserve privacy each machine has only access to a distinct partition of the overall graph. The central server has access to the full graph. In the LLCG algorithm that the paper proposes, each machine trains on its local graph partition for some time before sending the parameters to the server. The server averages the received parameter, but additionally also does its own training using the full graph available to it. Theoretically the authors show that the proposed method avoids an error gap in the gradient norm that would exist if server correction is not performed. Experimental results show the proposed scheme performing similarly to GGS albeit with much lower communication costs.  

Training GNNs is challenging due to high communication costs or large memory overheads. This paper proposes a communication-efficient distributed GNN training technique named Learn Locally, Correct Globally (LLCG) to periodically model averaging on the server using locally trained models.  It also applies global server corrections to refine the locally learned models and solve the irreducible performance degradation caused by ignoring node dependency. This paper provides the convergence analysis and shows the proposed method can address the residual error. The experimental results show significant improvement compared to existing methods.

This paper deals with the problem of distributed training of GNNs. Existing methods are either communication-intensive (sampling) or do not achieve good performance (averaging).  The authors propose a novel method, dubbed "LLCG: Learn Locally Correct Globally". Essentially, this method captures the idea of transmitting only local averages but adds a centralized step on the server to account for global structural information lost in the subgraph partition.  The authors further provide theoretical convergence guarantees. They both show that just averaging leads to an insurmountable residual error that explains the poor performance of averaging methods, as well as prove that this residual error disappears when adding the global correction step.

This paper proposes a distributed training technique for GNN. This technique includes local computations done in parallel by several machines and a correction phase done by a centralized server. A theoretical analysis is given for this technique, showing that the server correction phase reduces some irreducible error that happens due to splitting the graph and doing a local computation on each subgraph. Several experiments are made on real datasets which show the merits of this technique over previous techniques in terms of performance, communication steps, and size.