MULTI-LEVEL LOCAL SGD: DISTRIBUTED SGD FOR HETEROGENEOUS HIERARCHICAL NETWORKS
Year: 2021

Abstract:
We propose Multi-Level Local SGD, a distributed stochastic gradient method for learning a smooth, non-convex objective in a multi-level communication network with heterogeneous workers. Our network model consists of a set of disjoint subnetworks, with a single hub and multiple workers; further, workers may have different operating rates. The hubs exchange information with one another via a connected, but not necessarily complete, communication network. In our algorithm, sub-networks execute a distributed SGD algorithm, using a hub-and-spoke paradigm, and the hubs periodically average their models with neighboring hubs. We first provide a unified mathematical framework that describes the Multi-Level Local SGD algorithm. We then present a theoretical analysis of the algorithm; our analysis shows the dependence of the convergence error on the worker node heterogeneity, hub network topology, and the number of local, sub-network, and global iterations. We illustrate the effectiveness of our algorithm in a multi-level network with slow workers via simulation-based experiments.

heading: 1 INTRODUCTION
text: Stochastic Gradient Descent (SGD) is a key algorithm in modern Machine Learning and optimization (Amari, 1993). To support distributed data as well as reduce training time, Zinkevich et al. (2010) introduced a distributed form of SGD. Traditionally, distributed SGD is run within a huband-spoke network model: a central parameter server (hub) coordinates with worker nodes. At each iteration, the hub sends a model to the workers. The workers each train on their local data, taking a gradient step, then return their locally trained model to the hub to be averaged. Distributed SGD can be an efficient training mechanism when message latency is low between the hub and workers, allowing gradient updates to be transmitted quickly at each iteration. However, as noted in Moritz et al. (2016), message transmission latency is often high in distributed settings, which causes a large increase in overall training time. A practical way to reduce this communication overhead is to allow the workers to take multiple local gradient steps before communicating their local models to the hub. This form of distributed SGD is referred to as Local SGD (Lin et al., 2018; Stich, 2019). There is a large body of work that analyzes the convergence of Local SGD and the benefits of multiple local training rounds (McMahan et al., 2017; Wang & Joshi, 2018; Li et al., 2019).
Local SGD is not applicable to all scenarios. Workers may be heterogeneous in terms of their computing capabilities, and thus the time required for local training is not uniform. For this reason, it can be either costly or impossible for workers to train in a fully synchronous manner, as stragglers may hold up global computation. However, the vast majority of previous work uses a synchronous model, where all clients train for the same number of rounds before sending updates to the hub (Dean et al., 2012; Ho et al., 2013; Cipar et al., 2013). Further, most works assume a hub-and-spoke model, but this does not capture many real world settings. For example, devices in an ad-hoc network may not all be able to communicate to a central hub in a single hop due to network or communication range limitations. In such settings, a multi-level communication network model may be beneficial. In flying ad-hoc networks (FANETs), a network architecture has been proposed to improve scalability by partitioning the UAVs into mission areas (Bekmezci et al., 2013). Here, clusters of UAVs have their own clusterheads, or hubs, and these hubs communicate through an upper level network, e.g., via satellite. Multi-level networks have also been utilized in Fog and Edge computing, a paradigm de-
∗T. Castiglia, A. Das, and S. Patterson are with the Department of Computer Science, Rensselaer Polytechnic Institute, 110 8th St, Troy, NY 12180, castit@rpi.edu, dasa2@rpi.edu, sep@cs.rpi.edu.
signed to improve data aggregation and analysis in wireless sensor networks, autonomous vehicles, power systems, and more (Bonomi et al., 2012; Laboratory, 2017; Satyanarayanan, 2017).
Motivated by these observations, we propose Multi-Level Local SGD (MLL-SGD), a distributed learning algorithm for heterogeneous multi-level networks. Specifically, we consider a two-level network structure. The lower level consists of a disjoint set of hub-and-spoke sub-networks, each with a single hub server and a set of workers. The upper level network consists of a connected, but not necessarily complete, hub network by which the hubs communicate. For example, in a Fog Computing application, the sub-network workers may be edge devices connected to their local data center, and the data centers act as hubs communicating over a decentralized network. Each subnetwork runs one or more Local SGD rounds, in which its workers train for a local training period, followed by model averaging at the sub-network’s hub. Periodically, the hubs average their models with neighbors in the hub network. We model heterogeneous workers using a stochastic approach; each worker executes a local training iteration in each time step with a probability proportional to its computational resources. Thus, different workers may take different numbers of gradient steps within each local training period. Note since MLL-SGD averages every local training period, regardless of how many gradient steps each worker takes, slow workers do not slow algorithm execution.
We prove the convergence of MLL-SGD for smooth and potentially non-convex loss functions. We assume data is distributed in an IID manner to all workers. Further, we analyze the relationship between the convergence error and algorithm parameters and find that, for a fixed step size, the error is quadratic in the number of local training iterations and the number of sub-network training iterations, and linear in the average worker operating rate. Our algorithm and analysis are general enough to encompass several variations of SGD as special cases, including classical SGD (Amari, 1993), SGD with weighted workers (McMahan et al., 2017), and Decentralized Local SGD with an arbitrary hub communication network (Wang & Joshi, 2018). Our work provides novel analysis of a distributed learning algorithm in a multi-level network model with heterogeneous workers.
The specific contributions of this paper are as follows. 1) We formalize the multi-level network model with heterogeneous workers, and we define the MLL-SGD algorithm for training models in such a network. 2) We provide theoretical analysis of the convergence guarantees of MLLSGD with heterogeneous workers. 3) We present an experimental evaluation that highlights our theoretical convergence guarantees. The experiments show that in multi-level networks, MLL-SGD achieves a marked improvement in convergence rate over algorithms that do not exploit the network hierarchy. Further, when workers have heterogeneous operating rates, MLL-SGD converges more quickly than algorithms that require all workers to execute the same number of training steps in each local training period.
The rest of the paper is structured as follows. In Section 2, we discuss related work. Section 3 introduces the system model and problem formulation. We describe MLL-SGD in Section 4, and we present our main theoretical results in Section 5. Proofs of these results are deferred to the appendix. We provide experimental results in Section 6. Finally, we conclude in Section 7.

heading: 2 RELATED WORK
text: Distributed SGD is a well studied subject in Machine Learning. Zinkevich et al. (2010) introduced parallel SGD in a hub-and-spoke model. Variations on Local SGD in the hub-and-spoke model have been studied in several works (Moritz et al., 2016; Zhang et al., 2016; McMahan et al., 2017). Many works have provided convergence bounds of SGD within this model (Wang et al., 2019b; Li et al., 2019). There is also a large body of work on decentralized approaches for optimization using gradient based methods, dual averaging, and deep learning (Tsitsiklis et al., 1986; Jin et al., 2016; Wang et al., 2019a). These previous works, however, do not address a multi-level network structure.
In practice, workers may be heterogeneous in nature, which means that they may execute training iterations at different rates. Lian et al. (2017) addressed this heterogeneity by defining a gossipbased asynchronous SGD algorithm. In Stich (2019), workers are modeled to take gradient steps at an arbitrary subset of all iterations. However, neither of these works address a multi-level network model. Grouping-SGD (Jiang et al., 2019) considers a scenario where workers can be clustered into groups, for example, based on their operating rates. Workers within a group train in a synchronous manner, while the training across different groups may be asynchronous. The system model differs
significantly from that in MLL-SGD in that as the model parameters are partitioned vertically across multiple hubs, and workers communicate with every hub.
Several recent works analyze Hierarchical Local SGD (HL-SGD), an algorithm for training a model in a hierarchical network. Different from MLL-SGD, HL-SGD assumes the hub network topology is a hub-and-spoke and also that workers are homogeneous. Zhou & Cong (2019) and Liu et al. (2020) analyze the convergence error of HL-SGD, while Abad et al. (2020) analyzes convergence time. Unlike HL-SGD, MLL-SGD accounts for an arbitrary hub communication graph, and MLL-SGD algorithm execution does not slow down in the presence of heterogeneous worker operating rates.
Several other works seek to encapsulate many variations of SGD under a single framework. Koloskova et al. (2020) created a generalized model that considers a gossip-based decentralized SGD algorithm where the communication network is time-varying. However, this work does not account for a multi-level network model nor worker heterogeneity. Wang et al. introduced the Cooperative SGD framework (Wang & Joshi, 2018), a model that includes communication reduction through local SGD steps and decentralized mixing between homogeneous workers. Cooperative SGD also allows for auxiliary variables. These auxiliary variables can be used to model SGD in a multi-level network, but only when sub-network averaging is immediately followed by hubs averaging with their neighbors in the hub network. Our model is more general; it considers heterogeneous workers and it allows for an arbitrary number of averaging rounds within each sub-network between averaging rounds across sub-networks, which is more practical in multi-level networks where inter-hub communication is slow or costly.

heading: 3 SYSTEM MODEL AND PROBLEM FORMULATION
text: In this section, we introduce our system model, the objective function that we seek to minimize, and the assumptions we make about the function.
We consider a set of D sub-networks D = {1, . . . , D}. Each sub-network d ∈ D has a single hub and a set of workersM(d), with |M(d) | = N (d). Workers inM(d) only communicate with their own hub and not with any other workers or hubs. We define the set of all workers in the system as M = ⋃D d=1M
(d). Let |M | = N . Each worker i holds a set S(i) of local training data. Let S = ⋃N i=1 S
(i). The set of all D hubs is denoted C. The hubs communicate with one another via an undirected, connected communication graph G = (C, E). Let Nd = {j | ed,j ∈ E} denote the set of neighbors of the hub in sub-network d in the hub graph G.
Let the model parameters be denoted by x ∈ Rn. Our goal is to find an x that minimizes the following objective function over the training set:
F (x) = 1 | S | ∑ s∈S f(x; s) (1)
where f(·) is the loss function. The workers collaboratively minimize this loss function, in part by executing local iterations of SGD over their training sets. For each executed local iteration, a worker samples a mini-batch of data uniformly at random from its local data. Let ξ be a randomly sampled mini-batch of data and let g(x; ξ) = 1|ξ| ∑ s∈ξ∇f(x; s) be the mini-batch gradient. For simplicity, we use g(x) instead of g(x; ξ) from here on. Assumption 1. The objective function and the mini-batch gradients satisfy the following:
1a The objective function F : Rn → R is continuously differentiable, and the gradient is Lipschitz with constant L > 0, i.e., ‖∇F (x)−∇F (y)‖2 ≤ L‖x− y‖2 for all x,y ∈ Rn.
1b The function F is lower bounded, i.e., F (x) ≥ Finf > −∞ for all x ∈ Rn.
1c The mini-batch gradients are unbiased, i.e., Eξ|x[g(x)] = ∇F (x) for all x ∈ Rn.
1d There exist scalars β ≥ 0 and σ ≥ 0 such that Eξ|x‖g(x)−∇F (x)‖22 ≤ β||∇F (x)||22+σ2 for all x ∈ Rn.
Assumption 1a requires that the gradients do not change too rapidly, and Assumption 1b requires that our objective function is lower bounded by some Finf . Assumptions 1c and 1d assume that
Algorithm 1 Multi-Level Local SGD
1: Initialize: y(d)1 for hubs d = 1, . . . , D 2: for k = 1, . . . ,K do 3: parallel for d ∈ D do 4: parallel for i ∈M(d) do 5: x(i)k ← y (d) k . Workers receive updated model from hub 6: for j = k, . . . , k + τ − 1 do 7: x(i)k+1 ← x (i) k − ηg (i) k . Local iteration (probabilistic) 8: end for 9: end parallel for
10: z(d) ← ∑ i∈M(d) v (i)x (i) k+1 . Hub d computes average of its workers’ models 11: if k mod q · τ = 0 then 12: y(d)k+1 ← ∑ j∈N (d) Hj,dz
(j) . Hub d averages its model with neighboring hubs 13: else 14: y(d)k+1 ← z(d) 15: end if 16: end parallel for 17: end for
the local data at each worker can be used as an unbiased estimate for the full dataset with the same bounded variance. These assumptions are common in convergence analysis of SGD algorithms (e.g., Bottou et al. (2018)).

heading: 4 ALGORITHM
text: We now present our Multi-Level Local SGD (MLL-SGD) algorithm. The pseudocode is shown in Algorithm 1. Each sub-network trains in parallel and, periodically, the hubs average their models with neighboring hubs. The steps corresponding to Local SGD are shown in lines 5-10. Each hub and worker stores a copy of the model. For worker i ∈M(d), we denote its copy of the local model by x(i). We denote the model at hub d by y(d). The hub first sends its model to its workers, and the workers update their local models to match their hub’s model. Workers then execute multiple local training iterations, shown in line 7, to refine their local models independently. To represent the different rates of computation at each worker, we use a probabilistic approach. We assume that, in expectation, a worker i execute τ (i) local iterations for every τ time steps (τ (i) ≤ τ ). We thus define the N-vector p where each entry pi = τ (i)
τ is the probability with which worker i executes a local gradient step in each iteration k. Worker i updates its local model at iteration k as follows:
x (i) k+1 = x (i) k − ηg (i) k (2)
where η is the step size and g(i)k is a random variable such that
g (i) k =
{ g(x
(i) k ) w/ probability pi
0 w/ probability 1− pi. (3)
After τ time steps, the hub updates its model based on the models of its workers (line 10). For each worker i, we assign a positive weight w(i). Let v(i) be the weight for worker i normalized within its sub-network: v(i) = w
(i)∑ j∈M(d(i)) w (j) , where d(i) denotes the sub-network of worker i.
Each hub’s updates its model to be a weighted average over the workers’ models in its sub-network: y(d) = ∑ i∈M(d) v
(i)x(i). Weights may be assigned for different reasons. If all worker gradients are treated equally, then w(i) = 1 and v(i) = 1
N(d(i)) . We may also weight a worker’s gradient
proportional to its local dataset size, in which case w(i) = | S(i) | and v(i) = | S (i) |∑
r∈M(d(i)) | S (r) | . The
latter approach is used in Federated Averaging (McMahan et al., 2017).
After q iterations of Local SGD in each sub-network (q ·τ time steps), the hubs average their models with their neighbors in the hub communication network (line 12). The weight assigned to each hub’s
model is defined by a D ×D matrix H so that: y(d) = ∑
j∈N (d) Hj,dy
(j). (4)
Define the total weight in the network to be wtot = ∑ i∈M w
(i). Let b be a D-vector with each component d given by bd = ( ∑ i∈M(d) w
(i))/wtot. We assume H meets the following requirements. Assumption 2. The matrix H satisfies the following:
2a If (i, j) ∈ E, then Hi,j > 0. Otherwise, Hi,j = 0.
2b H is column stochastic, i.e., ∑D i=1 Hi,j = 1.
2c For all i, j ∈ D, we have biHi,j = bjHj,i.
Assumption 2 implies that H has one as a simple eigenvalue, with corresponding right eigenvector b and left eigenvector 1D. Further, all of its other eigenvalues have magnitude strictly less than 1 (since G is connected) (Rotaru & Nägeli, 2004). By defining H in this way, we ensure that the contributions from the workers’ gradients in each hub are incorporated in proportion to the workers’ weights. This weighted averaging approach allows us to naturally extend Federated Averaging to the multi-level network model.

heading: 5 ANALYSIS
text: We note that hubs are essentially stateless in MLL-SGD, as the hub models are copied to all workers after each sub-network or hub averaging. Thus, our analysis focuses on how worker models evolve. We first present an equivalent formulation of the MLL-SGD algorithm in terms of the evolution of the worker models. We then present our main result on the convergence of MLL-SGD.
The system behavior can be summarized by the following update rule for worker models:
Xk+1 = (Xk −ηGk) Tk (5)
where n×N matrix Xk = [x(1)k , . . . ,x (N) k ], n×N matrix Gk = [g (1) k , . . . , g (N) k ], andN×N matrix Tk is a time-varying operator that captures the three stages in MLL-SGD: local iterations, hub-andspoke averaging within each sub-network, and averaging across the hub network. We define Tk as follows:
Tk =  Z if k mod qτ = 0 V if k mod τ = 0 and k mod qτ 6= 0 I otherwise.
(6)
For local iterations, Tk = I, as there are no interactions between workers or hubs. For sub-network averaging, V is an N × N block diagonal matrix, with each block V(d) corresponding to a single sub-network d. The matrix V(d) is an N (d) ×N (d) matrix where each entry is V(d)i,j = v(i). Finally, we define an N ×N matrix Z that captures the sub-network averaging and hub network averaging in one operation that involves all workers. The components of Z are given by
Zi,j = Hd(i),d(j)v(i). (7)
Let a be an N -vector with each component ai = w (i)
wtot representing the weight of worker i, nor-
malized over all worker weights. We observe that Z and V satisfy the following: each have a right eigenvector of a and left eigenvector of 1TN with eigenvalue 1 and all other eigenvalues have magnitude strictly less than 1. The proof of these properties can be found in the appendix. These properties are necessary (but not sufficient) to ensure that the worker models converge to a consensus model, where each worker’s updates have been incorporated according to the worker’s weight.
As is common, we study an averaged model over all workers in the system (Yuan et al., 2016; Wang & Joshi, 2018). Specifically, we define a weighted average model:
uk = Xk a. (8)
We identify the recurrence relation of uk. If we multiply a on both sides of (5):
Xk+1 a = (Xk −ηGk) Tk a (9) uk+1 = uk − ηGk a (10)
uk+1 = uk − η N∑ i=1 aig (i) k (11)
where (10) follows from a being a right eigenvector of V and Z with eigenvalue 1. We note that uk is updated via a stochastic gradient descent step using a weighted average of several mini-batch gradients. Since F (·) may be non-convex, SGD may converge to a local minimum or saddle point. Thus, we study the gradients of uk as k increases.
We next provide the main theoretical result of the paper. Theorem 1. Under Assumptions 1 and 2, if η satisfies the following for all i ∈M:
(4pi − p2i − 2) ≥ ηL ( aipi(β + 1)− aip2i + p2i ) + 8L2η2q2τ2Γ (12)
where Γ = ζ1−ζ2 + 2 1−ζ + ζ (1−ζ)2 and ζ = max{|λ2(H)|, |λN (H)|}, then the expected square norm of the average model gradient, averaged over K iterations, is bounded as follows:
E
[ 1
K K∑ k=1 ‖∇F (uk)‖22
] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi
+ 4L2η2σ2q3τ3 ( 1
qτ − 1 K
)( ζ2
1− ζ2 +
2ζ
1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P (13)
K→∞−−−−→ σ2ηL N∑ i=1 a2ipi + 4L 2η2σ2q2τ2 ( ζ2 1− ζ2 + 2ζ 1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P (14)
where P = ∑N i=1 aipi.
The proof of Theorem 1 is provided in the appendix. The first term in (13) is the same as in centralized SGD (Bottou et al., 2018). As K → ∞, this term goes to zero. The second term is similar to centralized SGD as well. If the stochastic gradients have high variance, then the convergence error will be larger. This term is also related to the convergence error in distributed SGD (Bottou et al., 2018), which is equivalent to MLL-SGD when there is one sub-network, q = τ = 1, ai = 1/N , and pi = 1 for all i. MLL-SGD has a dependence on the probabilities of gradient steps and worker weights, replacing the 1N in the equivalent term in distributed SGD.
The third and fourth terms in (13) are additive errors that depend on the topology of the hub network. The value of ζ is given by the second largest eigenvalue of H, by magnitude, which is an indication of the sparsity of the hub network. When worker weights are uniform, a fully connected hub graph G will have ζ = 0, while a sparse G will typically have ζ close to 1. It is interesting to note that ζ only depends on H, and not Z or V, meaning the convergence error does not depend on how worker weights are distributed within sub-networks.
We also note the third and fourth terms depend on P, the weighted average probability of the workers. The convergence error increases as the average worker operating rate increases. This relation is expected as more local iterations will increase convergence error (Wang & Joshi, 2018). It is interesting to note that the convergence error does not depend on the distribution of p, meaning that a skewed and uniform distribution with the same average probability would have the same convergence error. We observe that the condition on η in (12) cannot always be satisfied given certain probabilities. Specifically, when there exists a pi ≤ 2 − √ 2 ≈ 0.59, then the left-hand side will be non-positive, and the inequality can no longer be satisfied. Although this may be a conservative bound, intuitively, when pi’s are below this threshold, the algorithm may not make sufficient progress in each time step to guarantee convergence.
The third and fourth terms also grow with q and τ , the number of local iterations per hub network averaging and sub-network averaging steps, respectively. The longer workers train locally without reconciling their models, the more their models will diverge, leading to larger convergence error. We can see that τ plays a slightly larger role in convergence error than q. For a given q · τ , meaning a given number of time steps between hub averaging steps, a larger τ leads to higher convergence error than a larger q would. Thus, there is a slight penalty to performing more local iterations between sub-network averaging steps. We explore this more in Section 6.
We note that when setting ai = 1/N and pi = 1 for all workers i, and setting q = 1, MLL-SGD reduces to Cooperative SGD. However, the bound in Theorem 1 differs from that of Cooperative SGD. Specifically, Theorem 1 has error terms dependent on τ2 as opposed to τ in Cooperative SGD. This discrepancy is due to accommodating all possible values of pi. More details can be found in Appendix C.4.
In the following corollary, we analyze the convergence rate of Algorithm 1 when η = 1 L √ K .
Corollary 1. Let η = 1 L √ K
and let q2τ2 ≤ √ K. If qτ < K, then
E
[ 1
K K∑ k=1 ‖∇F (uk)‖22
] ≤ O ( L√ K ) (F (x1)− Finf ]) +O ( σ2√ K ) (15)
Under the conditions given in Corollary 1, MLL-SGD achieves the same asymptotic convergence rate as Local SGD and HL-SGD.

heading: 6 EXPERIMENTS
text: In this section, we show the performance of MLL-SGD compared to algorithms that do not account for hierarchy and heterogeneous worker rates. We also explore the impact of the different algorithm parameters that show up in Theorem 1.
We use the EMNIST (Cohen et al., 2017) and CIFAR-10 (Krizhevsky et al., 2009) datasets. For all experiments, we provide results for training a simple Convolutional Neural Network (CNN) on EMNIST and training ResNet-18 on CIFAR-10. The CNN has two convolutional layers and two fully connected layers. We train the CNN with a step size of 0.01. For ResNet, we use a standard approach of changing the step size from 0.1 to 0.01 to 0.001 over the course of training (He et al., 2016). We conduct experiments using Pytorch 1.4.0 and Python 3.
We compare MLL-SGD with Distributed SGD, Local SGD, and HL-SGD. Distributed SGD is equivalent to MLL-SGD when there is one hub, q = τ = 1, and ai = 1/N and pi = 1 for all i, which means Distributed SGD averages all worker models at every iteration. Thus, we use Distributed
SGD as a baseline for convergence error and accuracy in some experiments. Local SGD is equivalent to MLL-SGD when ai = 1/N and pi = 1 for all i, when the hub network is fully connected, and q = 1. HL-SGD extends Local SGD to allow q > 1. For all experiments, we let τ = 32 for Local SGD. We let qτ = 32 for all HL-SGD and MLL-SGD variations to be comparable with Local SGD. In all experiments, we measure training loss and test accuracy of the averaged model uk every 32 iterations.
We first explore the effect of different values of τ and q in MLL-SGD. We configure a multi-level network with a fully connected hub network and with 10 hubs, each with 10 workers. We use two configurations for MLL-SGD, one with τ = 8 and q = 4, and one with τ = 4 and q = 8. Distributed SGD and Local SGD treat the hubs as pass-throughs, and average all workers every iteration and every τ iterations respectively. Workers are split into five groups of 20 workers each. Each group is assigned a percentage of the full dataset: 5%, 10%, 20%, 25%, and 40%. Workers within a group partition the data evenly. The workers weights are assigned based on dataset sizes. In Figures 1a and 1c we plot the training loss, and in Figures 1b and 1d we plot the test accuracy for the CNN and ResNet, respectively. We observe that as q increases, while keeping qτ = 32, MLL-SGD improves and approaches the Distributed SGD baseline. Thus, increasing the number of sub-network training rounds improves the convergence behavior of MLL-SGD. The benefit is more pronounced in training ResNet on CIFAR.
We next investigate how the number and sizes of the sub-networks impacts the convergence of MLLSGD. From a pool of 100 workers, we distribute them across 5, 10, and 20 sub-networks. The hub network is a path graph, which yields the largest ζ while keeping the network connected. This hub network topology the worst-case scenario in terms of the convergence bound. Note that as the number of hubs increases, the larger ζ becomes. We let ai = 1/N and pi = 1 for all workers i. We set q = 4 and τ = 8. We also include results using Local SGD with 1 hub and 100 workers. The results of this experiment are shown in Figures 2 and 3. In the case of the CNN, the difference in training loss is minimal among the MLL-SGD variations. In the case of ResNet we can see that as the number of hubs increase, the convergence rate decreases. This is in line with Theorem 1 since an increased number hubs corresponds with an increased ζ. Interestingly, despite the low hub network connectivity, MLL-SGD outperforms Local SGD. This shows that MLL-SGD still benefits from a hierarchy even when hub connectivity is sparse.
Next, we explore the impact of different distributions of worker operating rates. According to Theorem 1, the average probability across workers plays a role in the error bound. To see if this holds in practice, we compare four different MLL-SGD setups, all of which includes a complete hub network, 10 hubs, each with 10 workers, ai = 1/N , and an average probability amongst workers of 0.55: (i) all workers with a pi = 0.55 (Fixed); (ii) workers in each sub-network with probability ranging from 0.1 to 1 at steps of 0.1 (Uniform Distribution); (iii) 90 workers with pi = 0.5 and 10 workers with pi = 1 (Skewed 1); (iv) 90 workers with pi = 0.6 and 10 workers with pi = 0.1 (Skewed 2). We include a case where all workers have pi = 1 as a baseline (Prob=1). In Figures 4
and 5 we can see that in all cases except the baseline, the convergence rate is similar in both models. This is in line with our theoretical results, since all cases have the same average worker probability.
Finally, we compare the convergence time of MLL-SGD against algorithms that wait for slower workers: Local SGD and HL-SGD. We simulate real-time with time slots. In every time slot, each worker will take a gradient step with a probability pi. Note when pi = 1 for a worker i, the number of gradient steps taken will match the number of time slots T . Otherwise, the number of gradient steps taken will be T · pi in expectation. MLL-SGD will wait τ time slots before averaging worker models in a sub-network, regardless of the number of gradient steps taken, while Local SGD and HL-SGD will wait for all workers to take τ gradient steps. This approach allows us to compare the progress of each algorithm over time. In this experiment, we set pi = 0.9 for 90% of workers and pi = 0.6 for 10% of the workers. As in the previous experiments, we use a multi-level network with a fully connected hub network and with 10 hubs, each with 10 workers. We study MLL-SGD with two parameter settings, τ = 32, q = 1 and τ = 8, q = 4. We also include results for Local-SGD and HL-SGD. By comparing MLL-SGD with τ = 32, q = 1 with Local-SGD, we can evaluate the impact of using a local training period based on time rather than a number of worker iterations. By comparing MLL-SGD with τ = 8, q = 4 with HL-SGD, we can evaluate this impact in a multi-level network.
In Figures 6a and 6c, we plot the training loss, and in Figures 6b and 6d, we plot the test accuracy for the CNN and ResNet, respectively. We can see that MLL-SGD with q = 1 converges more quickly, in both loss and accuracy, than Local SGD, and that MLL-SGD with q = 4 converges more quickly than HL-SGD. These trends hold in both the CNN and ResNet models. The results show that in this experimental setup, waiting for slow workers is detrimental to the overall convergence time.

heading: 7 CONCLUSION
text: We have introduced MLL-SGD, a variation of Distributed SGD in a multi-level network model. Our algorithm incorporates the heterogeneity of worker devices using a stochastic approach. We provide theoretical analysis of the algorithm’s convergence, and we show how the convergence error depends on the average worker rate, the hub network topology, and the number of local, subnetwork averaging, and hub averaging steps. Finally, we provide experimental results that illustrate the effectiveness of MLL-SGD over Local SGD and HL-SGD. In future work, we plan to analyze the effects of non-IID data on convergence error.

heading: ACKNOWLEDGMENTS
text: This work is supported by the Rensselaer-IBM AI Research Collaboration (http://airc.rpi.edu), part of the IBM AI Horizons Network (http://ibm.biz/AIHorizons), and by the National Science Foundation under grants CNS 1553340 and CNS 1816307.

heading: A CODE REPOSITORY
text: The code used in our experiments can be found at: https://github.com/rpi-nsl/MLL-SGD. This code simulates a multi-level network with heterogeneous workers, and trains a model using MLL-SGD.

heading: B ADDITIONAL EXPERIMENTS
text: Our experiments in Section 6 explore how changing MLL-SGD parameters affect training on a nonconvex function. In this section, show the results of the same experiments on a convex loss function. We train a logistic regression model on the MNIST dataset (Bottou et al., 1994). We train a binary classification model with half the classes being 0 and the other half being 1 and use a step size of 0.2. We run all experiments for 32,000 iterations.
We rerun our first experiment from Figure 1 with logistic regression trained on MNIST. Figures 7a and 7b show the training loss and test accuracy, respectively. As with the non-convex functions, we can see that MLL-SGD with larger q approaches the Distributed SGD baseline.
We rerun our second experiment comparing different hub and worker distributions with logistic regression trained on MNIST. Figure 8 shows the training loss. The three variations of MLL-SGD do not show much difference in terms of convergence rate, indicating that ζ has little effect in this case. However, they still outperform Local SGD due to q being larger.
We rerun our third experiment comparing different worker operating rates distributions with logistic regression trained on MNIST. Figure 9 shows the training loss. As with the non-convex functions, all MLL-SGD variations with the same average probability have similar convergence rate.
We rerun our first experiment from Figure 6 with logistic regression trained on MNIST. Figures 10a and 10b show the training loss and test accuracy, respectively. We can see an improvement in convergence rate of MLL-SGD over both Local SGD and HL-SGD.

heading: C PROOF OF THEOREM 1
text: For our proof we adopt a similar approach to that in Wang & Joshi (2018). This section is structured as follows. We first define some notation and make some observations in Section C.1. Our supporting lemmas are stated in Section C.2. We close with the full proof of Theorem 1 in Section C.3.

heading: C.1 PRELIMINARIES
text: For simplicity of notation, we let ‖ · ‖ denote the l2 vector norm. Let the weighted Frobenius norm of an N ×M matrix X with an N -vector a be defined as follows:
‖X‖2Fa = ∣∣∣Tr((diag(a))1/2 XXT (diag(a))1/2)∣∣∣ = N∑
i=1 M∑ j=1 ai|xi,j |2. (16)
The matrix operator norm for a square matrix Q is defined as:
‖Q‖op = √ λmax(QT Q). (17)
We define the set of Bernoulli random variables Θ = {θ1k, . . . , θNk }, where
θik = { 1 with probability pi 0 with probability (1− pi).
Let Ξk = {ξ(1)k , ..., ξ (N) k } be the set of mini-batches used by the N workers at time step k. Without loss of generality, we assign a mini-batch to each worker, even if it does not execute a gradient step in that iteration. An equivalent definition of g(i)k is then
g (i) k = θ i kg(ξ (i) k ). (18)
For simplicity of notation, let Ek be equivalent to EΘk,Ξk|Xk .
We note that Assumption 1c implies:
Ek[g(i)k ] = piEk[g(x (i) k )] (19)
= pi∇F (x(i)k ). (20)
Further, when i 6= j:
Ek[(g(i)k ) Tg (j) k ] = pipjEk[(g(x (i) k ))] TEk[g(x(j)k )] (21)
= pipj∇F (x(i)k ) T∇F (x(j)k ). (22)
We also note that Assumption 1d implies: Ek ∥∥∥g(i)k −∇F (x(i)k )∥∥∥2 = Ek [∥∥∥g(i)k ∥∥∥2 + ∥∥∥∇F (x(i)k )∥∥∥2 − 2(g(i)k )T∇F (x(i)k )] (23)
= Ek ∥∥∥g(i)k ∥∥∥2 + ∥∥∥∇F (x(i)k )∥∥∥2 − 2Ek(g(i)k )T∇F (x(i)k ) (24)
= piEk ∥∥∥g(x(i)k )∥∥∥2 + ∥∥∥∇F (x(i)k )∥∥∥2 − 2piEkg(x(i)k )T∇F (x(i)k ) (25)
= piEk ∥∥∥g(x(i)k )∥∥∥2 + pi ∥∥∥∇F (x(i)k )∥∥∥2 − 2piEkg(x(i)k ) T∇F (x(i)k ) + (1− pi) ∥∥∥∇F (x(i)k )∥∥∥2 (26)
= piEk ∥∥∥g(x(i)k )−∇F (x(i)k )∥∥∥2 + (1− pi)∥∥∥∇F (x(i)k )∥∥∥2 (27)
≤ piβ ∥∥∥∇F (x(i)k )∥∥∥2 + piσ2 + (1− pi)∥∥∥∇F (x(i)k )∥∥∥2 (28)
= (pi(β − 1) + 1) ∥∥∥∇F (x(i)k )∥∥∥2 + piσ2. (29)
Finally, we define the weighted average stochastic gradient and the weighted average batch gradient as:
Gk = N∑ i=1 aig (i) k , Hk = N∑ i=1 ai∇F (x(i)k ).

heading: C.2 LEMMAS AND PROPOSITIONS
text: Next, we state our supporting lemmas and propositions. Proposition 1. The matrices Z and V satisfy the following properties:
1. Z and V each have a right eigenvector of a with eigenvalue 1.
2. Z and V each have a left eigenvector of 1TN with eigenvalue 1.
3. All other eigenvalues of Z and V have magnitude strictly less than 1.
Proof. Assumption 2 indicates that H is a Generalized Diffusion Matrix as defined in Rotaru & Nägeli (2004).
Recall Assumption 2:
Assumption 2. The matrix H satisfies the following:
2a If (i, j) ∈ E, then Hi,j > 0. Otherwise, Hi,j = 0.
2b H is column stochastic, i.e., ∑D i=1 Hi,j = 1.
2c For all i, j ∈ D, we have biHi,j = bjHj,i.
If we show this implies that Z and V are Generalized Diffusion Matrices with the same properties to those in Assumption 2 with vector a, then the properties in the proposition are satisfied.
Since H and b are non-negative, then Z is also non-negative. It is also clear that Z is column stochastic by construction. It is left to prove that:
Zi,j aj = Zj,i ai. (30)
Applying the definition of Z to the left side, we have:
Zi,j aj = Hd(i),d(j)v(i)aj (31)
Since we know that H is a Generalized Diffusion Matrix with vector b, we know that:
Hi,jbj = Hj,ibi (32)
Hi,j = Hj,i bi bj . (33)
Plugging this in for Hd(i),d(j), we have:
Zi,j aj = Hd(j),d(i) bd(i)
bd(j) v(i)aj (34)
= Hd(j),d(i)
∑ r∈M(d(i)) w (r)
wtot
wtot∑ r∈M(d(j)) w (r) w(i)∑ r∈M(d(i)) w (r) w(j) wtot (35)
= Hd(j),d(i) w(i)
wtot w(j)∑ r∈M(d(j)) w (r) (36)
= Hd(j),d(i)v (j)ai (37)
= Zj,i ai. (38)
Therefore, Z is a Generalized Diffusion Matrix.
We can show that V is also a Generalized Diffusion Matrix with the vector a. V is constructed to be non-negative and column stochastic. It is left to prove that
Vi,j aj = Vj,i ai. (39)
When i, j are outside a block V(d), then Vi,j = Vj,i = 0, so the equation is trivially satisfied. When within a block, in terms of w, we have:
Vi,j aj = Vj,i ai (40) w(i)∑
r∈M(d(i)) w (r)
w(j) wtot = w(j)∑ r∈M(d(j)) w (r) w(i) wtot . (41)
Noting that we are within a block, therefore d(i) = d(j), we can see that both sides are equal:
w(i)w(j) = w(j)w(i). (42)
Therefore, V is a Generalized Diffusion Matrix.
Proposition 2. Given a diffusion matrix H with the properties in Assumption 2, if Z constructed as follows,
Zi,j = Hd(i),d(j)v(i) (43) then the largest eigenvalues of Z are the eigenvalues of H, and zero otherwise.
Proof. In order to prove the relationship of the eigenvalues of Z and H, we prove the following two points separately:
1. The rank of Z is the same as H.
2. All non-zero eigenvalues of H are eigenvalues of Z with the same multiplicity.
For the rank of Z, we take a look at how each column is constructed. Consider column j of Z:
Zj = [H1,d(j)v(1), . . . ,H1,d(j)v(N (1)),H2,d(j)v (N(1)+1), . . . ,HD,d(j)v (N)]T . (44)
For two columns i and j where d(i) = d(j), these columns are identical. Therefore, the rank of Z will be, at most, the number of hubs, D. Further, we can see that the elements of a column j in Z are simply scaled elements of column d(j) in H. So any linearly dependent columns in H will also be linearly dependent in Z. Therefore, the rank of the two matrices are the same.
For the second point, we show there is a bijective mapping from eigenpairs of H to eigenpairs of Z. Let (λ,y) be an eigenpair of H (with λ 6= 0), i.e.
Hy = λy. (45)
Define the N -vector x with components xi = v(i)yd(i). We will show that Zx = λx. Looking at the i-th entry of the vector Zx, we have
(Zx)i = N∑ j=1 Zi,j xj . (46)
Applying the definition of Z and x, we obtain
(Zx)i = N∑ j=1 1 v(i) Hd(i),d(j)v (j)yd(j) (47)
= 1
v(i) D∑ l=1 Hd(i),lyl ∑ k∈M(l) v(k) (48)
= 1
v(i) D∑ l=1 Hd(i),lyl. (49)
Note that the m-th entry of the vector Hy equals ∑D l=1 Hm,lyl = λym. Applying this equality, we obtain
(Zx)i = 1
v(i) λyd(i) (50)
= λxi. (51)
Therefore, for any eigenpair (λ,y) of H, we can find an eigenpair (λ,x) of Z. It is left to prove that this mapping is a bijection.
Suppose eigenvalue λ of H has multiplicity k > 1. We consider any two of the k eigenpairs (λ, c) and (λ,d). Let the corresponding eigenpairs of Z be (λ, e) and (λ,f). We know that e 6= f because c and d are unique, and there must exist an index i such that v(i)cd(i) 6= v(i)dd(i). Therefore, the mapping of eigenpairs of H to eigenpairs of Z is a bijection.
Proposition 3. Given definition of Z and V in Proposition 1, it is the case that
Z V = V Z = Z . (52)
Proof. First, we prove that V Z = Z. Note that the i-th row of V contains either vi or zero. Looking at an arbitrary entry i, j of V Z we have:
(V Z)i,j = v(i) ∑
r∈Md(i)
Zr,j (53)
(V Z)i,j = v(i)Hd(i),d(i) (54) (V Z)i,j = Zi,j . (55)
Next we prove that Z V = Z. Note that for any row i in Z, Zi,j = Zi,k when d(j) = d(k).
(Z V)i,j = Zi,j N∑ r=1 Vr,j . (56)
Since V is column stochastic:
(Z V)i,j = Zi,j . (57)
Proposition 4. Let A = a 1T . Given our definition of Tk in (6),
Tk A = A Tk = A (58)
for all k.
Proof. We prove each of the three cases of Tk: I, V, and Z. Clearly, I A = A I = A. It is left to prove V A = A V = A and Z A = A Z = A.
We can see that Z A = A since a is a right eigenvector of Z with eigenvalue 1: Z A = Za 1T = a 1T = A. Similarly, we can see that A Z = a 1T Z = a 1T = A as 1T is a left eigenvector of Z. The same holds for V.
Lemma 1. Under Assumptions 1c and 1d, the variance of the weighted average stochastic gradient is bounded as follows:
Ek[‖Gk −Hk‖2] ≤ N∑ i=1 a2i [ (pi(β − 1) + 1) ∥∥∥∇F (x(i)k )∥∥∥2 + piσ2]
+ N∑ l=1 N∑ j 6=l alaj(1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k ). (59)

heading: Proof.
text: Ek[‖Gk −Hk‖2] = Ek ∥∥∥∥∥ N∑ i=1 ai(g (i) k −∇F (xk)) ∥∥∥∥∥ 2  (60)
= Ek  N∑ i=1 a2i ‖g (i) k −∇F (xk)‖ 2 + N∑ l=1 N∑ j 6=l alaj 〈 g (j) k −∇F (x (j) k ), g (l) k −∇F (x (l) k ) 〉 (61)
= N∑ i=1 a2iEk‖g (i) k −∇F (xk)‖ 2 + N∑ l=1 N∑ j 6=l alajEk [〈 g (j) k −∇F (x (j) k ), g (l) k −∇F (x (l) k ) 〉] .
(62)
Looking at the cross-terms in (62):
Ek [〈 g (j) k −∇F (x (j) k ), g (l) k −∇F (x (l) k ) 〉]
= Ek [ (g (j) k ) Tg (l) k ] − Ek [ (g (j) k ) T∇F (x(l)k ) ]
− Ek [ ∇F (x(j)k ) Tg (l) k ] +∇F (x(j)k ) T∇F (x(j)k ) (63)
= pjpl∇F (x(j)k ) T∇F (x(l)k )− pj∇F (x (j) k ) T∇F (x(l)k )
− pl∇F (x(j)k ) T∇F (x(l)k ) +∇F (x (j) k ) T∇F (x(j)k ) (64)
= (1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k ). (65)
Plugging (65) into (62) we have:
Ek[‖Gk −Hk‖2] = N∑ i=1 a2iEk‖g (i) k −∇F (xk)‖ 2
+ N∑ l=1 N∑ j 6=l alaj(1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k ) (66)
≤ N∑ i=1 a2i [ (pi(β − 1) + 1) ∥∥∥∇F (x(i)k )∥∥∥2 + piσ2]
+ N∑ l=1 N∑ j 6=l alaj(1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k ) (67)
where (67) follows from Assumption 1d and (29).
Lemma 2. Under Assumptions 1c and 1d, the squared norm of the stochastic gradients is bounded by:
Ek [ ‖Gk‖2 ] ≤ N∑ i=1 [ a2i ( pi(β + 1)− p2i ) + aip 2 i ] ∥∥∥∇F (x(i)k )∥∥∥2 + N∑ i=1 a2ipiσ 2 (68)

heading: Proof.
text: Ek [ ‖Gk‖2 ] = Ek [ ‖Gk − Ek[Gk]‖2 ] + ‖Ek[Gk]‖2 (69)
= Ek [ ‖Gk −
N∑ i=1 aipi∇F (x(i)k )‖ 2
] + ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2
(70)
= Ek ∥∥∥∥∥Gk − N∑ i=1 ai∇F (x(i)k ) + N∑ i=1 ai(1− pi)∇F (x(i)k ) ∥∥∥∥∥ 2 + ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2 . (71)
Applying the definition of Gk to (71) we get: Ek [ ‖Gk‖2 ] = Ek
∥∥∥∥∥ N∑ i=1 [ aig (i) k − ai∇F (x (i) k ) + ai(1− pi)∇F (x (i) k ) ]∥∥∥∥∥ 2 + ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2 (72)
= Ek [ N∑ i=1 ∥∥∥aig(i)k − ai∇F (x(i)k ) + ai(1− pi)∇F (x(i)k )∥∥∥2 ] + ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2
+ Ek [ N∑ j=1 N∑ l=1,l 6=j alaj 〈 (g (j) k −∇F (x (j) k )) + (1− pj)∇F (x (j) k ),
(g (l) k −∇F (x (l) k )) + (1− pl)∇F (x (l) k ) 〉] . (73)
Let the cross-terms in (73) be
CR = Ek [ N∑ j=1 N∑ l=1,l 6=j alaj 〈 (g (j) k −∇F (x (j) k )) + (1− pj)∇F (x (j) k ),
(g (l) k −∇F (x (l) k )) + (1− pl)∇F (x (l) k ) 〉] . (74)
We can simplify CR as follows:
CR = N∑ j=1 N∑ l=1,l 6=j alajEk [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k ))
+ (g (j) k −∇F (x (j) k )) T (1− pl)∇F (x(l)k ) + (1− pj)∇F (x (j) k ) T (g (l) k −∇F (x (l) k ))
+ (1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k )
T ]
(75)
= N∑ j=1 N∑ l=1,l 6=j alaj [ Ek [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] + (Ek[g(j)k ]−∇F (x (j) k )) T (1− pl)∇F (x(l)k )
+ (1− pj)∇F (x(j)k ) T (Ek[g(l)k ]−∇F (x (l) k ))
+ (1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k )
T ] . (76)
Applying Assumption 1c to (76), we get:
CR = N∑ j=1 N∑ l=1,l 6=j alaj [ Ek [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] + (pj − 1)(1− pl)∇F (x(j)k ) T∇F (x(l)k ) + (pl − 1)(1− pj)∇F (x (j) k ) T∇F (x(l)k )
+ (1− pj)(1− pl)∇F (x(j)k ) T∇F (x(l)k )
] (77)
= N∑ j=1 N∑ l=1,l 6=j alaj [ Ek [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] − (1− pl)(1− pj)∇F (x(j)k ) T∇F (x(l)k ) ] . (78)
Applying (78) to (73) we have:
Ek [ ‖Gk‖2 ] = Ek [ N∑ i=1 ∥∥∥aig(i)k − ai∇F (x(i)k ) + ai(1− pi)∇F (x(i)k )∥∥∥2 ]
+ N∑ j=1 N∑ l=1,l 6=j alaj [ Ek [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] − (1− pl)(1− pj)∇F (x(j)k ) T∇F (x(l)k ) ]
+ ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2 . (79)
Expanding the first term in (79) we have:
Ek [ ‖Gk‖2 ] = Ek [ N∑ i=1 ∥∥∥aig(i)k − ai∇F (x(i)k )∥∥∥2 ] + N∑ i=1 a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2
+ N∑ i=1 Ek [〈 aig (i) k − ai∇F (x (i) k ),ai(1− pi)∇F (x (i) k ) 〉]
︸ ︷︷ ︸ CR1
+ N∑ i=1 Ek [〈 ai(1− pi)∇F (x(i)k ),aig (i) k − ai∇F (x (i) k ) 〉]
︸ ︷︷ ︸ CR2
+ N∑ j=1 N∑ l=1,l 6=j alaj [ Ek [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] − (1− pl)(1− pj)∇F (x(j)k ) T∇F (x(l)k ) ]
+ ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2 . (80)
We simplify CR1:
CR1 = [〈
aipiEk(g(x(i)k )) T − ai∇F (x(i)k ),ai(1− pi)∇F (x (i) k ) 〉]
(81) = [〈 ai(pi − 1)∇F (x(i)k ),ai(1− pi)∇F (x (i) k ) 〉] (82)
= −a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2 . (83)
Similarly, for CR2:
CR2 = −a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2 . (84)
Plugging (83) and (84) back into (80):
Ek [ ‖Gk‖2 ] = Ek [ N∑ i=1 ∥∥∥aig(i)k − ai∇F (x(i)k )∥∥∥2 ] − N∑ i=1 a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2
+ N∑ j=1 N∑ l=1,l 6=j alaj [ Ek [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] − (1− pl)(1− pj)∇F (x(j)k ) T∇F (x(l)k ) ]
+ ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2 . (85)
We can simplify by observing that:
Ek [ ‖Gk −Hk‖2 ] = Ek [ N∑ i=1 ∥∥∥aig(i)k − ai∇F (x(i)k )∥∥∥2 ]
+ N∑ j=1 N∑ l=1,l 6=j alajEk [ (g (j) k −∇F (x (j) k )) T (g (l) k −∇F (x (l) k )) ] (86)
which gives us:
Ek [ ‖Gk‖2 ] = Ek [ ‖Gk −Hk‖2 ] − N∑ i=1 a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2 + ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2
− N∑ j=1 N∑ l=1,l 6=j alaj(1− pl)(1− pj)∇F (x(j)k ) T∇F (x(l)k ) (87)
Applying Lemma 1 to (87):
Ek [ ‖Gk‖2 ] ≤ N∑ i=1 a2i [ (pi(β − 1) + 1) ∥∥∥∇F (x(i)k )∥∥∥2 + piσ2]
− N∑ i=1 a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2 + ∥∥∥∥∥ N∑ i=1 aipi∇F (x(i)k ) ∥∥∥∥∥ 2 (88)
≤ N∑ i=1 a2i [ (pi(β − 1) + 1) ∥∥∥∇F (x(i)k )∥∥∥2 + piσ2]
− N∑ i=1 a2i (1− pi)2 ∥∥∥∇F (x(i)k )∥∥∥2 + N∑ i=1 aip 2 i ∥∥∥∇F (x(i)k )∥∥∥2 (89) =
N∑ i=1 [ a2i ( pi(β + 1)− p2i ) + aip 2 i ] ∥∥∥∇F (x(i)k )∥∥∥2 + N∑ i=1 a2ipiσ 2 (90)
where equation (89) follows from Jensen’s inequality.
Lemma 3. Under Assumption 1c, the expected inner product of the batch gradient and the weighted average stochastic gradient is equal to:
Ek[〈∇F (uk),Gk〉] = 1
2 ‖∇F (uk)‖2 + N∑ i=1 ai 2 ∥∥∥pi∇F (x(i)k )∥∥∥2 −
N∑ i=1 ai 2 ∥∥∥∇F (uk)− pi∇F (x(i)k )∥∥∥2 (91)

heading: Proof.
text: Ek[〈∇F (uk),Gk〉] = Ek [〈 ∇F (uk),
N∑ i=1 aig (i) k
〉] (92)
= 〈 ∇F (uk),
N∑ i=1 piai∇F (x(i)k )
〉 (93)
= N∑ i=1 ai 〈 ∇F (uk),pi∇F (x(i)k ) 〉 (94)
= N∑ i=1 ai 2 [ ‖∇F (uk)‖2 + ‖pi∇F (x(i)k )‖ 2 − ‖∇F (uk)− pi∇F (x(i)k )‖ 2 ] (95)
= 1
2 ‖∇F (uk)‖2 + N∑ i=1 ai 2 ∥∥∥pi∇F (x(i)k )∥∥∥2 −
N∑ i=1 ai 2 ∥∥∥∇F (uk)− pi∇F (x(i)k )∥∥∥2 (96) where (93) follows from (20), and (95) follows from the fact that, for arbitrary vectors y and z, 2yTz = ||y||2 + ||z||2 − ||y − z||2.
Lemma 4. Under Assumption 1, following the update rule given in (5), if all model parameters are initialized at the same x1, the expected weighted average gradient is bounded as follows:
E
[ 1
K K∑ k=1
‖∇F (uk)‖2 ] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi + 2L2 K K∑ k=1 E‖Xk(I−A)‖2Fa
− 1 K K∑ k=1 N∑ i=1 ai ( (4pi − p2i − 2)− ηL ( aipi(β + 1)− aip2i + p2i )) E ∥∥∥∇F (x(i)k )∥∥∥2 (97)
where A = a 1T .
Proof. According to Assumption 1a, Ek[F (uk+1)]− F (uk) ≤ Ek [ 〈∇F (uk)),uk+1 − uk〉+ L
2 ‖uk+1 − uk‖22
] (98)
= −ηEk[〈∇F (uk),Gk〉] + η2L 2 Ek [ ‖Gk‖22 ] . (99)
Plugging in Lemmas 2 and 3, we get:
Ek[F (uk+1)]− F (uk)
≤ −η
[ 1
2 ‖∇F (uk)‖2 + N∑ i=1 ai 2 p2i ∥∥∥∇F (x(i)k )∥∥∥2 − N∑ i=1 ai 2 ‖∇F (uk)− pi∇F (x(i)k )‖ 2
]
+ η2L
2 N∑ i=1 ( a2ipi(β + 1)− a2ip2i + aip2i ) ∥∥∥∇F (x(i)k )∥∥∥2 + σ2η2L2 N∑ i=1 a2ipi (100)
= −η 2 ‖∇F (uk)‖2 + η 2 N∑ i=1 ai‖∇F (uk)− pi∇F (x(i)k )‖ 2 + σ2η2L 2 N∑ i=1 a2ipi
− η 2 N∑ i=1 ai ( p2i − ηL ( aipi(β + 1)− aip2i + p2i )) ∥∥∥∇F (x(i)k )∥∥∥2 . (101)
After some rearranging, we obtain:
‖∇F (uk)‖2 ≤ 2 (F (uk)− Ek[F (uk+1)])
η + σ2ηL N∑ i=1 a2ipi + N∑ i=1 ai‖∇F (uk)− pi∇F (x(i)k )‖ 2
− N∑ i=1 ai ( p2i − ηL ( aipi(β + 1)− aip2i + p2i )) ∥∥∥∇F (x(i)k )∥∥∥2 . (102)
Taking the total expectation over all iterations:
E
[ 1
K K∑ k=1
‖∇F (uk)‖2 ] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi
+ 1
K K∑ k=1 N∑ i=1 aiE‖∇F (uk)− pi∇F (x(i)k )‖ 2
− 1 K K∑ k=1 N∑ i=1 ai ( p2i − ηL ( aipi(β + 1)− aip2i + p2i )) E ∥∥∥∇F (x(i)k )∥∥∥2 .
(103)
The third term in (103) can be bounded as:
N∑ i=1 aiE‖∇F (uk)− pi∇F (x(i)k )‖ 2
= N∑ i=1 aiE‖∇F (uk)−∇F (x(i)k ) + (1− pi)∇F (x (i) k )‖ 2 (104)
≤ N∑ i=1 [ 2aiE‖∇F (uk)−∇F (x(i)k )‖ 2 + 2ai(1− pi)2E‖∇F (x(i)k )‖ 2 ] (105)
≤ N∑ i=1 2aiL 2E‖uk − x(i)k ‖ 2 + N∑ i=1 2ai(1− pi)2E‖∇F (x(i)k )‖ 2 (106)
where (105) follows from the fact that ‖y + z‖2 ≤ 2‖y‖2 + 2‖z‖2, and (106) follows from (105) by Assumption 1a.
Recalling the definition of the weighted Frobenius norm and the definition of u, we can simplify the first term in (106):
N∑ i=1 2aiL 2E‖uk − x(i)k ‖ 2 = 2L2E‖uk1T − Xk ‖2Fa (107)
= 2L2E ∥∥Xk a1T − Xk∥∥2Fa (108) = 2L2E‖Xk(I−A)‖2Fa . (109)
Plugging (106) and (109) back into (103), we obtain:
E
[ 1
K K∑ k=1
‖∇F (uk)‖2 ] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi
+ 2L2
K K∑ k=1 E‖Xk(I−A)‖2Fa + 1 K K∑ k=1 N∑ i=1 2ai(1− pi)2E‖∇F (x(i)k )‖ 2
− 1 K K∑ k=1 N∑ i=1 ai ( p2i − ηL ( aipi(β + 1)− aip2i + p2i )) E ∥∥∥∇F (x(i)k )∥∥∥2 (110)
= 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi + 2L2 K K∑ k=1 E‖Xk(I−A)‖2Fa
− 1 K K∑ k=1 N∑ i=1 ai ( (4pi − p2i − 2)− ηL ( aipi(β + 1)− aip2i + p2i )) E ∥∥∥∇F (x(i)k )∥∥∥2 . (111)
Lemma 5. Given the properties of Z and V given in Propositions 1 and 2, it is the case that:∥∥Zj −A∥∥ op = ζj , ‖V−A‖op = 1, ‖I−A‖op = 1 (112)
where A = a 1T and ζ = max{|λ2(H)|, |λ(H)|}.
Proof. According to the definition of the matrix operator norm,∥∥Zj −A∥∥ op = √ λmax((Zj −A)T (Zj −A)) (113)
= √ λmax(Z2j −A Zj −Zj A + A) (114)
= √ λmax(Z2j −A) (115)
where (114) follows from Aj = A, and (115) follows from A Z = Z A = A.
We can simplify (115) further:
= √ λmax(Z2j −A2j) (116)
= √ λmax(Z−A)2j (117)
where (117) follows from the commutability of Z and A.
Based on Proposition 2, the non-zero eigenvalues of Z are the same as H. As shown in Lemma 6 of Rotaru & Nägeli (2004), for a matrix Z with the properties in Proposition 1, the spectral norm of Z−A is equal to ζ. Therefore: ∥∥Zj −A∥∥ op = √ ζ2j (118)
= ζj . (119)
Similarly for V:
‖V−A‖op = √ λmax((V−A)T (V−A)) (120)
= √ λmax(V−A V−V A + A) (121)
= √ λmax(V−A) (122)
(123)
where (121) follows from Aj = A and Vj = V, and (122) follows from A V = V A = A.
Note that the eigenvalues of each block V(d) are N (d) − 1 zeros and a one. The set of eigenvalues of V will include D ones. If D > 1, then based on Lemma 6 of Rotaru & Nägeli (2004) and Proposition 1, the spectral norm of V−A is 1, so
‖V−A‖op = √ 1 (124)
= 1. (125)
Since the eigenvalues of I are all 1, and I is commutable with A, we can similarly say: ‖I−A‖op = 1. (126)
Lemma 6. Given two matrices C ∈ RN×M and D ∈ RM×N , and an N -vector a,∣∣∣Tr((diag(a))1/2 C D (diag(a))1/2)∣∣∣ ≤ ‖C ‖Fa‖D ‖Fa . (127) Proof. We define the i-th row of C as cTi and the i-th column of D as di. We can rewrite the trace as:
Tr((diag(a)) 1/2 C D (diag(a))1/2) = N∑ i=1 M∑ j=1 ai Ci,j Dj,i (128)
= N∑ i=1 aic T i di. (129)
Placing a squared norm around (129), we can apply the Cauchy-Schwartz inequality:∣∣∣∣∣ N∑ i=1 aic T i di ∣∣∣∣∣ 2 ≤ ( N∑ i=1 ai‖cTi ‖2 )( N∑ i=1 ai‖di‖2 )
(130)
=  N∑ i=1 M∑ j=1 ai C2i,j  N∑ i=1 M∑ j=1 ai D2i,j  (131) = ‖C ‖2Fa‖D ‖ 2 Fa . (132)
Lemma 7. Given two matrices C ∈ RM×N and D ∈ RN×N , and an N -vector a, then ‖C D ‖Fa ≤ ‖C ‖Fa‖D ‖op. (133)
Proof. We define the i-th row of C as cTi and the set I = {i ∈ [1,M ] : ‖cTi ‖ 6= 0}. We can rewrite the squared Frobenius norm as:
‖C D ‖2Fa = M∑ i=1 ‖cTi D (diag(a)) 1/2 ‖2 (134)
= M∑ i∈I ‖cTi D (diag(a)) 1/2 ‖2 (135)
= M∑ i∈I ‖cTi (diag(a)) 1/2 ‖2 ‖c T i D (diag(a)) 1/2 ‖2 ‖cTi (diag(a)) 1/2 ‖2
(136)
≤ M∑ i∈I ‖cTi (diag(a)) 1/2 ‖2‖D ‖2op (137)
= ‖C ‖2Fa‖D ‖ 2 op. (138)

heading: C.3 PROOF OF THEOREM 1
text: We recall Theorem 1.
Theorem 1. Under Assumptions 1 and 2, if η satisfies the following for all i ∈M: (4pi − p2i − 2) ≥ ηL ( aipi(β + 1)− aip2i + p2i ) + 8L2η2q2τ2Γ (12)
where Γ = ζ1−ζ2 + 2 1−ζ + ζ (1−ζ)2 and ζ = max{|λ2(H)|, |λN (H)|}, then the expected square norm of the average model gradient, averaged over K iterations, is bounded as follows:
E
[ 1
K K∑ k=1 ‖∇F (uk)‖22
] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi
+ 4L2η2σ2q3τ3 ( 1
qτ − 1 K
)( ζ2
1− ζ2 +
2ζ
1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P (13)
K→∞−−−−→ σ2ηL N∑ i=1 a2ipi + 4L 2η2σ2q2τ2 ( ζ2 1− ζ2 + 2ζ 1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P (14)
where P = ∑N i=1 aipi.
We now give the proof of Theorem 1 using Lemmas 2-7.
Proof. Using our intermediate result from Lemma 4, we decompose Xk(I−A) using our recursive definition of Xk:
Xk(I−A) = (Xk−1−ηGk−1) Tk−1(I−A) (139) = Xk−1(I−A) Tk−1−ηGk−1(Tk−1−A) (140) = [(Xk−2−ηGk−2) Tk−2(I−A)] Tk−1−ηGk−1(Tk−1−A) (141) = [Xk−2(I−A) Tk−2−ηGk−2(Tk−2−A)] Tk−1−ηGk−1(Tk−1−A) (142) = Xk−2(I−A) Tk−2 Tk−1−ηGk−2(Tk−2 Tk−1−A)− ηGk−1(Tk−1−A). (143)
where (140) follows from the commutability of Tk and A by Proposition 4.
Continuing this, we end up with:
Xk(I−A) = X1(I−A) k−1∏ l=1 Tl−η k−1∑ s=1 Gs ( k−1∏ l=s Tl−A ) . (144)
Since all workers initialize their models to the same vector, X1(I−A) ∏k−1 l=1 Tk = 0, and thus we have:
E ‖Xk(I−A)‖2Fa = η 2E ∥∥∥∥∥ k−1∑ s=1 Gs ( k−1∏ l=s Tl−A )∥∥∥∥∥ 2
Fa
. (145)
Let k = jqτ + lτ + f , where j is the number of hub network averaging rounds, l is the number of sub-network averaging rounds since the last hub network averaging round, and f is the number of local iterations since the last sub-network averaging round. Define:
Φs,k−1 = k−1∏ l=s Tl .
Noting that Vj = V, and V Z = Z V = Z by Proposition 3, Φs,k−1 can be expressed as:
Φs,k−1 =  I jqτ + lτ < s < jqτ + lτ + f V jqτ < s ≤ jqτ + lτ Z (j − 1)qτ < s ≤ jqτ Z2 (j − 2)qτ < s ≤ (j − 1)qτ ... Zj 1 ≤ s ≤ qτ.
(146)
For r < j, let
Yr = (r+1)qτ∑ s=rqτ+1 Gs , Qr = (r+1)qτ∑ s=rqτ+1 ∇F (Xs)
We also let Yj1 = ∑jqτ+lτ s=jqτ+1 Gs, Yj2 = ∑jqτ+lτ+f s=jqτ+lτ+1 Gs, Qj1 = ∑jqτ+lτ s=jqτ+1∇F (Xs), and Qj2 =∑jqτ+lτ+f
s=jqτ+lτ+1∇F (Xs). With this in mind, we can split the sum in (145) into batches for each hub network averaging period:
qτ∑ s=1 Gs (Φs,k−1 − A) = Y0(Zj −A) (147)
2qτ∑ s=qτ+1 Gs (Φs,k−1 − A) = Y1(Zj−1−A) (148)
...
jqτ∑ s=(j−1)qτ+1 Gs (Φs,k−1 − A) = Yj−1(Z−A) (149)
jqτ+lτ+f∑ s=jqτ+1 Gs (Φs,k−1 − A) = Yj1(V−A) + Yj2(I−A). (150)
Summing this all together, we get: k−1∑ s=1 Gs (Φs,k−1 − A) = j−1∑ r=0 Yr(Zj−r −A) + Yj1(V−A) + Yj2(I−A). (151)
Plugging (151) into (145):
E ‖Xk(I−A)‖2Fa = η 2E ∥∥∥∥∥ j−1∑ r=0 Yr(Zj−r −A) + Yj1(V−A) + Yj2(I−A) ∥∥∥∥∥ 2
Fa
(152)
= η2E ∥∥∥∥∥ j−1∑ r=0 (Yr −Qr)(Zj−r −A) + (Yj1 −Qj1)(V−A)
+ (Yj2 −Qj2)(I−A) + j−1∑ r=0 Qr(Z j−r −A) + Qj1(V−A) + Qj2(I−A) ∥∥∥∥∥ 2
Fa
(153)
≤ 2η2E ∥∥∥∥∥ j−1∑ r=0 (Yr −Qr)(Zj−r −A) + (Yj1 −Qj1)(V−A) + (Yj2 −Qj2)(I−A)) ∥∥∥∥∥ 2
Fa︸ ︷︷ ︸ T1
+ 2η2E ∥∥∥∥∥ j−1∑ r=0 Qr(Z j−r −A) + Qj1(V−A) + Qj2(I−A) ∥∥∥∥∥ 2
Fa︸ ︷︷ ︸ T2
(154)
where (154) follows from the fact that ‖y + z‖2 ≤ 2‖y‖2 + 2‖z‖2. We first put a bound on T1:
T1 = 2η 2E ∥∥∥∥∥ j−1∑ r=0 (Yr −Qr)(Zj−r −A) + (Yj1 −Qj1)(V−A) + (Yj2 −Qj2)(I−A) ∥∥∥∥∥ 2
Fa
(155)
= 2η2 ( j−1∑ r=0 E ∥∥(Yr −Qr)(Zj−r −A)∥∥2Fa + E∥∥(Yj1 −Qj1)(V−A)∥∥2Fa
+ E ∥∥(Yj2 −Qj2)(I−A)∥∥2Fa
)
+ 2η2 j−1∑ n=0 j−1∑ l=0,l 6=n E ∣∣∣∣∣∣∣Tr ( (diag(a)) 1/2 (Zj−n−A)(Yn−Qn)T (Yl−Ql)(Zj−l−A) (diag(a)) 1/2 )
︸ ︷︷ ︸ TR ∣∣∣∣∣∣∣︸ ︷︷ ︸ TR0
+ 4η2 j−1∑ l=0 E ∣∣∣Tr((diag(a))1/2 (V−A)(Yj1 −Qj1)T (Yl−Ql)(Zj−l−A) (diag(a))1/2)∣∣∣︸ ︷︷ ︸
TR1 + 4η2 j−1∑ l=0 E ∣∣∣Tr((diag(a))1/2 (I−A)(Yj2 −Qj2)T (Yl−Ql)(Zj−l−A) (diag(a))1/2)∣∣∣︸ ︷︷ ︸
TR2 + 4η2E ∣∣∣Tr((diag(a))1/2 (V−A)(Yj1 −Qj1)T (Yj2 −Qj2)(I−A) (diag(a))1/2)∣∣∣︸ ︷︷ ︸
TR3
. (156)
TR can be bounded as:
TR ≤ ∥∥(Zj−n−A)(Yn−Qn)T∥∥Fa ∥∥(Yl−Ql)(Zj−l−A)∥∥Fa (157)
≤ ∥∥(Zj−n−A)∥∥
op ‖Yn−Qn‖Fa ‖Yl−Ql‖Fa ∥∥(Zj−l−A)∥∥ op
(158)
≤ ζ2j−n−l ‖Yn−Qn‖Fa ‖Yl−Ql‖Fa (159) ≤ 1 2 ζ2j−n−l [ ‖Yn−Qn‖ 2 Fa + ‖Yl−Ql‖ 2 Fa ] (160)
where (157) follows from Lemma 6, (158) follows from Lemma 7, and (159) follows from Lemma 5. We can similarly bound TR1 and TR3:
TR1 ≤ 2η2 j−1∑ l=0 ζj−l [ E ∥∥Yj1 −Qj1∥∥2Fa + E ‖Yl−Ql‖2Fa] (161)
TR2 ≤ 2η2 j−1∑ l=0 ζj−l [ E ∥∥Yj2 −Qj2∥∥2Fa + E ‖Yl−Ql‖2Fa] (162)
TR3 ≤ 2η2 [ E ∥∥Yj1 −Qj1∥∥2Fa + E∥∥Yj2 −Qj2∥∥Fa] . (163)
Summing TR0 through TR3, we get:
3∑ t=0 TRt ≤ η2 j−1∑ n=0 j−1∑ l=0,l 6=n ζ2j−n−l [ E ‖Yn−Qn‖ 2 Fa + E ‖Yl−Ql‖ 2 Fa ]
+ 2η2 j−1∑ l=0 ζj−l [ E ∥∥Yj1 −Qj1∥∥2Fa + E ‖Yl−Ql‖2Fa]
+ 2η2 j−1∑ l=0 ζj−l [ E ∥∥Yj2 −Qj2∥∥2Fa + E ‖Yl−Ql‖2Fa]
+ 2η2E ∥∥Yj1 −Qj1∥∥2Fa + 2η2E∥∥Yj2 −Qj2∥∥2Fa (164)
≤ 2η2 j−1∑ n=0 j−1∑ l=0,l 6=n ζ2j−n−lE ‖Yn−Qn‖ 2 Fa
+ 2η2 j−1∑ l=0 ζj−lE ‖Yl−Ql‖ 2 Fa + 2η2 j−1∑ l=0 ζj−lE ‖Yl−Ql‖ 2 Fa
+ 2η2 j∑ l=0 ζj−lE ∥∥Yj1 −Qj1∥∥2Fa + 2η2 j∑ l=0 ζj−lE ∥∥Yj2 −Qj2∥∥2Fa (165)
= 2η2 j−1∑ n=0 ζj−nE ‖Yn−Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l + 4η2 j−1∑ l=0 ζj−lE ‖Yl−Ql‖ 2 Fa
+ 2η2E ∥∥Yj1 −Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Yj2 −Qj2∥∥2Fa j∑
l=0
ζj−l (166)
where (165) follows from the symmetry of the n and l indices.
Plugging (166) back into (156):
T1 ≤ 2η2 j−1∑ r=0 E ‖(Yr −Qr)‖ 2 Fa ∥∥(Zj−r −A)∥∥2 op + 2η2E ∥∥(Yj1 −Qj1)∥∥2Fa ‖V−A‖2op
+ 2η2E ∥∥(Yj2 −Qj2)∥∥2Fa ‖I−A‖2op + 2η2 j−1∑
n=0
ζj−nE ‖Yn−Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l
+ 2η2E ∥∥Yj1 −Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Yj2 −Qj2∥∥2Fa j∑
l=0
ζj−l
+ 4η2 j−1∑ l=0 ζj−lE ‖Yl−Ql‖ 2 Fa
(167)
≤ 2η2 j−1∑ r=0 E ‖(Yr −Qr)‖ 2 Fa ζ2(j−r) + 2η2E ∥∥(Yj1 −Qj1)∥∥2Fa + 2η2E
∥∥(Yj2 −Qj2)∥∥2Fa + 2η2 j−1∑ n=0 ζj−nE ‖Yn−Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l
+ 2η2E ∥∥Yj1 −Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Yj2 −Qj2∥∥2Fa j∑
l=0
ζj−l
+ 4η2 j−1∑ l=0 ζj−lE ‖Yl−Ql‖ 2 Fa
(168)
where (167) follows from Lemma 7, and (168) follows from Lemma 5.
We further bound T1:
T1 ≤ 2η2 j−1∑ r=0 E ‖(Yr −Qr)‖ 2 Fa ζ2(j−r) + 2η2E ∥∥(Yj1 −Qj1)∥∥2Fa + 2η2E
∥∥(Yj2 −Qj2)∥∥2Fa + 2η2 j−1∑ n=0 ζj−nE ‖Yn−Qn‖ 2 Fa ζ 1− ζ
+ 2η2E ∥∥Yj1 −Qj1∥∥2Fa 11− ζ + 2η2E∥∥Yj2 −Qj2∥∥2Fa 11− ζ
+ 4η2 j−1∑ l=0 ζj−lE ‖Yl−Ql‖ 2 Fa
(169)
= 2η2 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) E ‖(Yr −Qr)‖ 2 Fa
+ 2η2 (
2− ζ 1− ζ ) E ∥∥Yj1 −Qj1∥∥2Fa + 2η2 ( 2− ζ 1− ζ ) E ∥∥Yj2 −Qj2∥∥2Fa (170)
where (169) follows from the summation formulae of a power series: j∑ l=0 ζj−l ≤ j∑ l=−∞ ζj−l ≤ 1 1− ζ , j−1∑ l=0 ζj−l ≤ j−1∑ l=−∞ ζj−l ≤ ζ 1− ζ . (171)
Taking a closer look at E ‖(Yr −Qr)‖ 2 Fa for 0 ≤ r < j:
E ‖(Yr −Qr)‖ 2 Fa = E ∥∥∥∥∥∥ (r+1)qτ∑ s=rqτ+1 (Gs−∇F (Xs)) ∥∥∥∥∥∥ 2
Fa
(172)
= N∑ i=1 aiE ∥∥∥∥∥∥ (r+1)qτ∑ s=rqτ+1 (gis −∇F (x(i)s )) ∥∥∥∥∥∥ 2
(173)
≤ N∑ i=1 aiqτ (r+1)qτ∑ s=rqτ+1 E ∥∥∥(gis −∇F (x(i)s ))∥∥∥2 (174)
≤ qτ  N∑ i=1 ai (r+1)qτ∑ s=rqτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2 + q2τ2σ2 N∑ i=1 aipi
(175)
= qτ  N∑ i=1 ai (r+1)qτ∑ s=rqτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2 + q2τ2σ2 P . (176) where (175) follows from Assumption 1d and (29).
Similarly, for r = j1 and r = j2:
E ∥∥(Yj1 −Qj1)∥∥2Fa ≤ lτ  N∑ i=1 ai jqτ+lτ∑ s=jqτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2 + l2τ2σ2 P (177) E ∥∥(Yj2 −Qj2)∥∥Fa ≤ (f − 1)  N∑ i=1 ai jqτ+lτ+f−1∑ s=jqτ+lτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2
 + (f − 1)2σ2 P . (178)
Plugging (176), (177), and (178) into (170), we can bound T1 as follows:
T1 ≤ 2η2σ2 ((
q2τ2 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ )) + ( 2− ζ 1− ζ )( l2τ2 + (f − 1)2 )) P
+ 2η2qτ j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) (r+1)qτ∑ s=rqτ+1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2
+ 2η2 (
2− ζ 1− ζ
) lτ  N∑ i=1 ai jqτ+lτ∑ s=jqτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2  + 2η2 ( 2− ζ 1− ζ ) (f − 1)  N∑ i=1 ai jqτ+lτ+f−1∑ s=jqτ+lτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2
 . (179) Referring back to Lemma 4, our goal is to sum T1 over k = 1, . . . ,K iterations. First, we sum over the j-th sub-network update period up to the j-th hub network averaging, for l = 0, . . . , q − 1 and f = 1, . . . , τ :
q−1∑ l=0 τ∑ f=1 T1 ≤ 2η2σ2q3τ3 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) P
+ 2η2σ2 (
2− ζ 1− ζ
)( τ3 q(q − 1)(2q + 1)
6 + q τ(τ − 1)(2τ + 1) 6
) P
+ 2η2q2τ2 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) (r+1)qτ∑ s=rqτ+1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2
+ η2 (
2− ζ 1− ζ
) q(q − 1)τ2 N∑ i=1 ai j(qτ+1)∑ s=jqτ+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2
+ η2 (
2− ζ 1− ζ
) q2τ(τ − 1) N∑ i=1 ai j(qτ+1)+τ−1∑ s=j(qτ+1)+1 (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2 . (180)
Let:
Γr =
( ζ2(j−r) + 2ζj−r + ζj−r+1
1− ζ
) . (181)
Note that Γj = 3−2ζ1−ζ > 2−ζ 1−ζ . Using this inequality, we can bound the sum of the last three terms of (180) to get 2q2τ2 ∑j r=0 Γr:
q−1∑ l=0 τ∑ f=1 T1 ≤ 2η2σ2q3τ3 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) P
+ 2η2σ2 (
2− ζ 1− ζ
)( τ3 q(q − 1)(2q + 1)
6 + q τ(τ − 1)(2τ + 1) 6
) P
+ 2η2q2τ2 j∑ r=0 Γr (r+1)qτ∑ s=rqτ+1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2 . (182)
Summing (182) over the hub network averaging periods j = 0, . . . ,K/(qτ)− 1, we obtain:
K/(qτ)−1∑ j=0 q−1∑ l=0 τ∑ f=1 T1 ≤ 2η2σ2q3τ3 K/(qτ)−1∑ j=0 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) P
+ 2η2σ2K ( 2− ζ 1− ζ )( τ2 (q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6 ) P + 2η2q2τ2 K/(qτ)−1∑
j=0 j∑ r=0 Γr (r+1)qτ∑ s=rqτ+1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2 (183)
= 2η2σ2q3τ3 K/(qτ)−2∑
r=0 K/(qτ)−1∑ j=r+1 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) P
+ 2η2σ2K ( 2− ζ 1− ζ )( τ2 (q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6 ) P + 2η2q2τ2 K/(qτ)−1∑
r=0 K/(qτ)−1∑ j=r Γj  (r+1)qτ∑ s=rqτ+1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)s )∥∥∥2  . (184)
Applying the following summation formula to sum over Γj , we obtain
K/(qτ)−1∑ j=r ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) ≤ ∞∑ j=r ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) (185)
≤ 1 1− ζ2 + 2 1− ζ +
ζ
(1− ζ)2 . (186)
We let Γ = 11−ζ2 + 2 1−ζ + ζ (1−ζ)2 . We can also apply this following summation formula to the first term in (184):
K/(qτ)−1∑ j=r+1 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) ≤ ∞∑ j=r+1 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) (187)
≤ ζ 2
1− ζ2 +
2ζ
1− ζ +
1
(1− ζ)2 . (188)
Applying the summation formula in (188), plugging Γ in, and indexing the iterations in terms of k, we bound (184) as:
K∑ k=1 T1 ≤ 2η2σ2q3τ3 ( K qτ − 1 )( ζ2 1− ζ2 + 2ζ 1− ζ +
1
(1− ζ)2
) P
+ 2η2σ2K ( 2− ζ 1− ζ )( τ2 (q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6 ) P
+ 2η2q2τ2Γ K∑ k=1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)k )∥∥∥2 . (189)
Now we bound T2:
T2 = 2η 2E ∥∥∥∥∥ j−1∑ r=0 Qr(Z j−r −A) + Qj1(V−A) + Qj2(I−A) ∥∥∥∥∥ 2
Fa
(190)
= 2η2 j−1∑ r=0 E ∥∥Qr(Zj−r −A)∥∥2Fa + 2η2E∥∥Qj1(V−A)∥∥2Fa + 2η2E∥∥Qj2(I−A)∥∥2Fa
+ 2η2 j−1∑ n=0 j−1∑ l=0,l 6=n E Tr((diag(a))1/2 (Zj−n−A) QTn Ql(Zj−l−A) (diag(a))1/2)︸ ︷︷ ︸ TR′  ︸ ︷︷ ︸
TR′0
+ 4η2 j−1∑ l=0 E [ Tr ( (diag(a)) 1/2 (V−A) QTj1 Ql(Z j−l−A) (diag(a))1/2 )] ︸ ︷︷ ︸
TR′1
+ 4η2 j−1∑ l=0 E [ Tr ( (diag(a)) 1/2 (I−A) QTj2 Ql(Z j−l−A) (diag(a))1/2 )] ︸ ︷︷ ︸
TR′2 + 4η2E [ Tr ( (diag(a)) 1/2 (V−A) QTj1 Qj2(I−A) (diag(a)) 1/2 )]
︸ ︷︷ ︸ TR′3 . (191)
TR′ can be bounded as:
TR′ ≤ ∥∥∥(Zj−n−A) QTn∥∥∥
Fa ∥∥Ql(Zj−l−A)∥∥Fa (192) ≤ ∥∥(Zj−n−A)∥∥
op ‖Qn‖Fa ‖Ql‖Fa ∥∥(Zj−l−A)∥∥ op
(193)
≤ 1 2 ζ2j−n−l
[ ‖Qn‖ 2 Fa + ‖Ql‖ 2 Fa ] (194)
where (192) follows from Lemma 6. We can similarly bound TR′1 through TR ′ 3:
TR′1 ≤ 2η2 j−1∑ l=0 ζj−l [ E ∥∥Qj1∥∥2Fa + E ‖Ql‖2Fa] (195)
TR′2 ≤ 2η2 j−1∑ l=0 ζj−l [ E ∥∥Qj2∥∥2Fa + E ‖Ql‖2Fa] (196)
TR′3 ≤ 2η2E ∥∥Qj1∥∥2Fa + 2η2E ∥∥Qj2∥∥2Fa . (197)
Summing TR′0 through TR ′ 3, we get:
3∑ t=0 TR′t ≤ η2 j−1∑ n=0 j−1∑ l=0,l 6=n ζ2j−n−lE [ E ‖Qn‖ 2 Fa + E ‖Ql‖ 2 Fa ]
+ 2η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa + 2η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa
+ 2η2 j∑ l=0 ζj−lE ∥∥Qj1∥∥2Fa + 2η2 j∑ l=0 ζj−lE ∥∥Qj2∥∥2Fa (198)
≤ 2η2 j−1∑ n=0 ζj−nE ‖Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l + 4η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa
+ 2η2E ∥∥Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Qj2∥∥2Fa j∑
l=0
ζj−l (199)
where (199) follows from the symmetry of the indices n and l.
Plugging (199) back into (191):
T2 ≤ 2η2 j−1∑ r=0 E ∥∥Qr(Zj−r −A)∥∥2Fa + 2η2E∥∥Qj1(V−A)∥∥2Fa + 2η2E∥∥Qj2(I−A)∥∥2Fa + 2η2
j−1∑ n=0 ζj−nE ‖Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l + 4η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa
+ 2η2E ∥∥Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Qj2∥∥2Fa j∑
l=0
ζj−l (200)
≤ 2η2 j−1∑ r=0 E ‖Qr‖ 2 Fa ∥∥(Zj−r −A)∥∥2 op + 2η2E ∥∥Qj1∥∥2Fa ‖V−A‖2op + 2η2E∥∥Qj2∥∥2Fa ‖I−A‖2op
+ 2η2 j−1∑ n=0 ζj−nE ‖Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l + 4η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa
+ 2η2E ∥∥Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Qj2∥∥2Fa j∑
l=0
ζj−l (201)
≤ 2η2 j−1∑ r=0 ζj−rE ‖Qr‖ 2 Fa + 2η2E ∥∥Qj1∥∥2Fa + 2η2E∥∥Qj2∥∥2Fa
+ 2η2 j−1∑ n=0 ζj−nE ‖Qn‖ 2 Fa j−1∑ l=0,l 6=n ζj−l + 4η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa
+ 2η2E ∥∥Qj1∥∥2Fa j∑
l=0
ζj−l + 2η2E ∥∥Qj2∥∥2Fa j∑
l=0
ζj−l (202)
where (201) follows from Lemma 7, and (202) follows from Lemma 5.
We further bound T2:
T2 ≤ 2η2 j−1∑ r=0 ζj−rE ‖Qr‖ 2 Fa + 2η2E ∥∥Qj1∥∥2Fa + 2η2E∥∥Qj2∥∥2Fa
+ 2η2 j−1∑ n=0 ζj−nE ‖Qn‖ 2 Fa ζ 1− ζ + 4η2 j−1∑ l=0 ζj−lE ‖Ql‖ 2 Fa
+ 2η2E ∥∥Qj1∥∥2Fa 11− ζ + 2η2E∥∥Qj2∥∥2Fa 11− ζ (203)
≤ 2η2 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) E ‖Qr‖ 2 Fa
+ 2η2 (
2− ζ 1− ζ ) E ∥∥Qj1∥∥2Fa + 2η2 ( 2− ζ 1− ζ ) E ∥∥Qj2∥∥2Fa (204)
where (203) follows from the summation formulae of a power series in (171).
After applying the definition of Q to (203), we obtain:
T2 = 2η 2 j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) E ∥∥∥∥∥ qτ∑ s=1 ∇F (Xrqτ+s) ∥∥∥∥∥ 2
Fa + 2η2 (
2− ζ 1− ζ
) E ∥∥∥∥∥ lτ∑ s=1 ∇F (Xjqτ+s) ∥∥∥∥∥ 2
Fa
+ 2η2 (
2− ζ 1− ζ
) E ∥∥∥∥∥ f−1∑ s=1 ∇F (Xjqτ+lτ+s) ∥∥∥∥∥ 2
Fa
(205)
≤ 2η2qτ j−1∑ r=0 ( ζ2(j−r) + 2ζj−r + ζj−r+1 1− ζ ) qτ∑ s=1 E ‖∇F (Xrqτ+s)‖2Fa
+ 2η2lτ ( 2− ζ 1− ζ ) lτ∑ s=1 E ‖∇F (Xjqτ+s)‖2Fa
+ 2η2(f − 1) (
2− ζ 1− ζ ) f−1∑ s=1 E ‖∇F (Xjqτ+lτ+s)‖2Fa (206)
where (206) follows from (205) by Jensen’s inequality.
Summing over all iterates in the j-th sub-network update period, we obtain:
q−1∑ l=0 τ∑ f=1 T2 ≤ 2η2q2τ2 j−1∑ r=0
(( ζ2(j−r) + 2ζj−r + ζj−r+1
1− ζ ) qτ∑ s=1 E ‖∇F (Xrqτ+s)‖2Fa )
+ η2qτ(q − 1) (
2− ζ 1− ζ ) qτ∑ s=1 E ‖∇F (Xjqτ+s)‖2Fa
+ η2qτ(τ − 1) (
2− ζ 1− ζ ) τ−1∑ s=1 E ‖∇F (Xjqτ+qτ+s)‖2Fa (207)
≤ 2η2q2τ2 j∑ r=0 Γr qτ∑ s=1 E ‖∇F (Xrqτ+s)‖2Fa . (208)
Summing over all iterations and applying the summation bound in (186) to (208):
K/(qτ)−1∑ j=0 q−1∑ l=0 τ∑ f=1 T2 ≤ 2η2q2τ2Γ K∑ k=1 E ‖∇F (Xk)‖2Fa . (209)
Summing T1 and T2, we obtain
2L2
K K∑ k=1 E‖Xk(I−A)‖2Fa ≤ 2L2 K K∑ k=1 T1 + 2L2 K K∑ k=1 T2 (210)
≤ 4L2η2σ2q3τ3 ( 1
qτ − 1 K
)( ζ2
1− ζ2 +
2ζ
1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P
+ 8L2η2q2τ2Γ 1
K K∑ k=1 N∑ i=1 ai (pi(β − 1) + 1)E ∥∥∥∇F (x(i)k )∥∥∥2 . (211)
Plugging T1 and T2 back into Lemma 4, we arrive at
E
[ 1
K K∑ k=1
‖∇F (uk)‖2 ] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi + 2L2 K K∑ k=1 T1 + 2L2 K K∑ k=1 T2
− 1 K K∑ k=1 N∑ i=1 ai ( (4pi − p2i − 2)− ηL ( aipi(β + 1)− aip2i + p2i )) E ∥∥∥∇F (x(i)k )∥∥∥2 (212)
= 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi
+ 4L2η2σ2q3τ3 ( 1
qτ − 1 K
)( ζ2
1− ζ2 +
2ζ
1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P
− 1 K K∑ k=1 N∑ i=1 ai ( (4pi − p2i − 2)− ηL ( aipi(β + 1)− aip2i + p2i ) − 8L2η2q2τ2Γ ) E ∥∥∥∇F (x(i)k )∥∥∥2
(213)
If η satisfies the following for i = 1, . . . , N , (4pi − p2i − 2) ≥ ηL ( aipi(β + 1)− aip2i + p2i ) + 8L2η2q2τ2Γ (214) then we can simplify (213):
E
[ 1
K K∑ k=1
‖∇F (uk)‖2 ] ≤ 2 (F (x1)− Finf ])
ηK + σ2ηL N∑ i=1 a2ipi
+ 4L2η2σ2q3τ3 ( 1
qτ − 1 K
)( ζ2
1− ζ2 +
2ζ
1− ζ +
1
(1− ζ)2
) P
+ 4L2η2σ2 (
2− ζ 1− ζ
)( τ2
(q − 1)(2q + 1) 6 + (τ − 1)(2τ + 1) 6
) P . (215)

heading: C.4 COMPARISON TO COOPERATIVE SGD
text: We note that when setting ai = 1/N and pi = 1 for all workers i, and setting q = 1, MLL-SGD reduces to Cooperative SGD (Wang & Joshi, 2018). However, the bound in Theorem 1 differs when compared to the bound of Cooperative SGD. Specifically, Theorem 1 has error terms dependent on τ2 as opposed to τ .
This is due to the formulation of g(i)k . Namely:
Ek[g(i)k ] = piEk[g(x (i) k )] (216)
= pi∇F (x(i)k ). (217)
Because we cannot assume pi = 1, there are cross terms in the expressions in equations (156) and (173) that do not cancel out. Thus, we needed to use a more conservative analysis at these steps on the proof. This is the reason that plugging in a value of pi = 1 is not enough to recover the same bound as in Cooperative SGD. A similar discrepancy can be observed when comparing with Koloskova et al. (2020).

Review summaries:
This paper proposes a new variant of local SGD algorithm to make it be more realistic. In particular, (1) it allows workers to perform different number of local steps, depending on their computational resources; (2) workers are organized in a multi-level structure. Workers connected to one central hub can synchronize frequently and hubs are communicated in an infrequent and decentralized manor.

This paper extends (Wang & Joshi, 2018) and proposes MLL-SGD for training models in hierarchic networks, where the network consists multiple sub-networks, and each sub-network contains multiple workers. In the level of sub-networks, models can be averaged. In the level of workers, the local copies of models can be averaged within a sub-network; however, workers cannot communicate directly with those from a different sub-networks. In such setting, MLL-SGD is proved to enjoy certain convergence property.