Year: 2021

Abstract:
In distributed learning, local SGD (also known as federated averaging) and its simple baseline minibatch SGD are widely studied optimization methods. Most existing analyses of these methods assume independent and unbiased gradient estimates obtained via with-replacement sampling. In contrast, we study shufflingbased variants: minibatch and local Random Reshuffling, which draw stochastic gradients without replacement and are thus closer to practice. For smooth functions satisfying the Polyak-Łojasiewicz condition, we obtain convergence bounds (in the large epoch regime) which show that these shuffling-based variants converge faster than their with-replacement counterparts. Moreover, we prove matching lower bounds showing that our convergence analysis is tight. Finally, we propose an algorithmic modification called synchronized shuffling that leads to convergence rates faster than our lower bounds in near-homogeneous settings.

heading: 1 INTRODUCTION
text: Distributed learning within the framework of federated learning (Konečnỳ et al., 2016; McMahan et al., 2017) has witnessed increasing interest recently. A key property of this framework is that models are trained locally using only private data on devices/machines distributed across a network, while parameter updates are aggregated and synchronized at a server.1 Communication is often the key bottleneck for federated learning, which drives the search for algorithms that can train fast while requiring less communication—see Li et al. (2020a); Kairouz et al. (2021) for recent surveys.
A basic algorithm for federated learning is local stochastic gradient descent (SGD), also known as federated averaging. The goal is to minimize the global objective that is an average of the local objectives. In local SGD, we have M machines and a server. After each round of communication, each of the M machines locally runs B steps of SGD on its local objective. Every B iterations, the server aggregates the updated local iterates from the machines, averages them, and then synchronizes the machines with the average. Convergence analysis of local SGD and its variants has drawn great interest recently (Dieuleveut & Patel, 2019; Haddadpour et al., 2019; Haddadpour & Mahdavi, 2019; Stich, 2019; Yu et al., 2019; Li et al., 2020b;c; Koloskova et al., 2020; Khaled et al., 2020; Spiridonoff et al., 2020; Karimireddy et al., 2020; Stich & Karimireddy, 2020; Qu et al., 2020).
Of the many, the biggest motivation for our paper comes from the line of work by Woodworth et al. (2020a;b; 2021). In (Woodworth et al., 2020a;b), minibatch SGD is studied as a simple yet powerful baseline for this intermittent communication setting. Instead of locally updating the iterates B times, minibatch SGD aggregates B gradients (evaluated at the last synced iterate) from each of the M machines, forms a minibatch of size MB, and then updates the shared iterate. Given the same M and B, local SGD and minibatch SGD have the same number of gradient computations per round of communication, so it is worthwhile to understand which converges faster. Woodworth et al. (2020a;b) point out that many existing analyses on local SGD show inferior convergence rate compared to minibatch SGD. Through their new upper and lower bounds, they identify regimes where local SGD can be faster than minibatch SGD.
While the theory of local and minibatch SGD has seen recent progress, there is still a gap between what is analyzed versus what is actually used. Most theoretical results assume independent and
1A distinctive feature of federated learning is that not all devices necessarily participate in the updates; however, we focus on the full participation setting in this paper.
unbiased gradient estimates obtained via with-replacement sampling of stochastic gradients (i.e., choosing training data indices uniformly at random). In contrast, most practitioners use withoutreplacement sampling, where they shuffle indices randomly and access them sequentially.
Convergence analysis of without-replacement methods is challenging because gradients sampled within an epoch lack independence. As a result, the standard theory based on independent gradient estimates does not apply to shuffling-based methods. While shuffling-based methods are believed to be faster in practice (Bottou, 2009), broad theoretical understanding of such methods remains elusive, except for noteworthy recent progress mainly focusing on the analysis of SGD (Gürbüzbalaban et al., 2019; Haochen & Sra, 2019; Nagaraj et al., 2019; Nguyen et al., 2020; Safran & Shamir, 2020; 2021; Rajput et al., 2020; 2021; Ahn et al., 2020; Mishchenko et al., 2020; 2021; Tran et al., 2021). These results indicate that in the large-epoch regime (where the number of epochs is greater than some threshold), without-replacement SGD converges faster than with-replacement SGD.

heading: 1.1 OUR CONTRIBUTIONS
text: We analyze convergence rates of without-replacement versions of local and minibatch SGD, where local component functions are reshuffled at every epoch. We call the respective algorithms local RR (Algorithm 1) and minibatch RR (Algorithm 2), and their with-replacement counterparts local SGD and minibatch SGD. Our key contributions are as follows:
• In Section 3, we present convergence bounds on minibatch and local RR for L-smooth functions satisfying the µ-Polyak-Łojasiewicz condition (Theorems 1 & 2). Our theorems give highprobability bounds, a departure from the common in-expectation bounds in the literature. We show that minibatch and local RR converge faster than minibatch and local SGD when the number of epochs is sufficiently large. We also identify a regime where local RR converges as fast as minibatch RR: when synchronization happens frequently enough and local objectives are not too heterogeneous. See also Appendix A for a detailed comparison with existing upper bounds.
• In Section 4, we prove that the upper bounds obtained in Section 3 are tight, in all factors except L and µ. We present Theorems 3 & 4 and Proposition 5 which show lower bounds that match the upper bound up to a factor of L2/µ2. Our lower bound on local RR indicates that if the synchronization interval B is too large, then local RR has no gain from parallel computation.
• In Section 5, we propose a simple modification called synchronized shuffling that allows us to bypass the lower bounds in Section 4, at the cost of a slight increase in communication. By having the server broadcast random permutations to local machines, we show that in near-homogeneous settings, the modified algorithms converge faster than the lower bounds (Theorems 6 & 7).
• In Appendix C, we present numerical experiments that corroborate our theoretical findings.

heading: 2 PROBLEM SETUP
text: Notation. For a natural number a ∈ N, let [a] := {1, 2, . . . , a}. Let Sa be the set of all permutations of [a]. Since our indices start from 1, we redefine the modulo operation between a ∈ Z and b ∈ N as a mod b := a− ba−1b cb, to make a mod b ∈ [b].
Optimization task. Consider M machines, each with its objective Fm(x) := 1N ∑N i=1 f m i (x), for m ∈ [M ]. The m-th machine has access only to the gradients of its own N local components fm1 (x), . . . , f m N (x). In this setting, we wish to minimize the global objective function which is an
average of the local objectives: F (x) := 1M ∑M m=1 F m(x) = 1MN ∑M m=1 ∑N i=1 f m i (x).
Further, we assume that each individual component function fmi is L-smooth, so that
fmi (y) ≤ fmi (x) + 〈∇fmi (x),y − x〉+ L2 ‖y − x‖ 2 , for all x,y ∈ Rd, (1)
and that the global objective F satisfies the µ-Polyak-Łojasiewicz (PŁ) condition.2
1 2 ‖∇F (x)‖ 2 ≥ µ(F (x)− F ∗) for all x ∈ Rd, where µ > 0. (2)
Algorithms. Under the above setting, we analyze local RR (Algorithm 1) and minibatch RR (Algorithm 2) and characterize their worst-case convergence rates.3 The algorithms are run overK epochs,
2PŁ functions can be thought as a nonconvex generalization of strongly convex functions. 3In Algorithms 1 and 2, consider SYNCSHUF as FALSE for now. We will discuss SYNCSHUF in Section 5.
Algorithm 1 Local RR (with and without SYNCSHUF) Input: Initialization y0, step-size η, # machines M , # components N , # epochs K, sync interval B.
1: Initialize xm1,0 := y0 for all m ∈ [M ]. 2: for k ∈ [K] do 3: if SYNCSHUF = TRUE then . Local RR with SYNCSHUF 4: Sample σ ∼ Unif(SN ), π ∼ Unif(SM ). 5: Set σmk (i) := σ((i+ N M π(m)) mod N) for all m ∈ [M ], i ∈ [N ]. 6: else . Local RR 7: Sample σmk ∼ Unif(SN ) independently and locally, for all m ∈ [M ]. 8: end if 9: for i ∈ [N ] do
10: for m ∈ [M ] do locally 11: Update xmk,i := x m k,i−1 − η∇fmσm
k (i)(x
m k,i−1).
12: end for 13: if B divides i then 14: Aggregate and average yk, i
B := 1 M
∑M m=1 x m k,i.
15: Synchronize xmk,i := yk, i B , for all m ∈ [M ]. 16: end if 17: end for 18: xmk+1,0 := yk,N
B , for all m ∈ [M ].
19: end for 20: return the last iterate yK,N
B .
Algorithm 2 Minibatch RR (with and without SYNCSHUF) Input: Initialization x0, step-size η, # machines M , # components N , # epochs K, sync interval B.
1: Initialize x1,0 := x0. 2: for k ∈ [K] do 3: if SYNCSHUF = TRUE then . Minibatch RR with SYNCSHUF 4: Sample σ ∼ Unif(SN ), π ∼ Unif(SM ). 5: Set σmk (i) := σ((i+ N M π(m)) mod N) for all m ∈ [M ], i ∈ [N ]. 6: else . Minibatch RR 7: Sample σmk ∼ Unif(SN ) independently and locally, for all m ∈ [M ]. 8: end if 9: for i ∈ [N
B ] do 10: Update xk,i := xk,i−1 − ηM ∑M
m=1
1 B ∑iB j=(i−1)B+1 ∇fmσm k
(j)(xk,i−1)︸ ︷︷ ︸ averaging done locally . 11: end for 12: xk+1,0 := xk,N
B .
13: end for 14: return the last iterate xK,N
B .
i.e., K passes over the entire component functions. At the beginning of epoch k, each machine m shuffles its local component functions {fmi }Ni=1 using a random permutation σmk ∼ Unif(SN ). In local RR, each machine makes B local RR updates to its iterate by sequentially accessing its shuffled component functions, before the server aggregates iterates from all the machines and then synchronizes the machines with the average iterate. In minibatch RR, instead of making B local updates, each machine collects B gradients evaluated at the last iterate, and the server aggregates them to make an update using theseMB gradients. Since these two algorithms use the same amount of communication and local gradients, minibatch RR is a simple yet powerful baseline for local RR.
Below, we collect our assumptions on the algorithm parameters used throughout the paper.
Assumption 1 (Algorithm parameters). We assume M ≥ 1, N ≥ 2, and K ≥ 1. Also, assume that B divides N . We restrict 1 ≤ B ≤ N2 for minibatch RR because B = N makes the algorithm equal to GD. We also assume 2 ≤ B ≤ N for local RR because B = 1 makes the two algorithms the same. We choose a constant step-size scheme, i.e., η > 0 is kept constant over all updates.
We next state assumptions on intra- and inter-machine deviations used in this paper.4
4Assumptions 2, 3 & 4 require that they hold for the whole Rd. We discuss ways to avoid it in Appendix D.7.
Assumption 2 (Intra-machine deviation). There exists ν ≥ 0 such that for allm ∈ [M ] and i ∈ [N ],
‖∇fmi (x)−∇Fm(x)‖ ≤ ν, for all x ∈ Rd.
Assumption 2 requires that the difference between the gradient of each local component function fmi (x) and its corresponding local objective function F
m(x) is uniformly bounded. It models the variance of local components fmi within each machine. While the uniform boundedness requirement may look strong, we use this assumption to prove high-probability upper bounds, which are stronger than the common in-expectation bounds. See Appendix A for comparisons with other assumptions, and also Appendix D.7 for ways to avoid uniform boundedness over the entire Rd.
The next two assumptions capture the deviation across different machines, i.e., the degree of heterogeneity, in two different levels of granularity: objective-wise and component-wise. Assumption 3 (Objective-wise inter-machine deviation). There exist τ ≥ 0 and ρ ≥ 1 such that
1 M ∑M m=1 ‖∇Fm(x)‖ ≤ τ + ρ ‖∇F (x)‖ , for all x ∈ Rd.
Assumption 3 models the heterogeneity by bounding the mean of ‖∇Fm‖ by a constant plus a multiplicative factor times ‖∇F‖. The assumption includes the homogeneous case (i.e., F 1 = · · · = FM = F ) by τ = 0 and ρ = 1. Assumption 3 is weaker than many other heterogeneity assumptions in the literature (e.g., Karimireddy et al. (2020)); see Appendix A for detailed comparisons.
Assumption 3 measures heterogeneity by only considering the local objectives Fm, not the local components fmi . We consider a more fine-grained notion of heterogeneity in Assumption 4:
Assumption 4 (Component-wise inter-machine deviation). For all i ∈ [N ], let f̄i := 1M ∑M m=1 f m i . There exist λ ≥ 0 such that for all m ∈ [M ] and i ∈ [N ],∥∥∇fmi (x)−∇f̄i(x)∥∥ ≤ λ, for all x ∈ Rd. Assumption 4 states that the gradients of the i-th components of local machines are “close” to each other. The assumption subsumes the component-wise homogeneous setting, i.e., f1i = f 2 i = · · · = fMi , by λ = 0. In distributed learning, this choice corresponds to the setting where each machine has the same training dataset. Assumption 4 with λ > 0 is also relevant to the case where each device has a slightly perturbed (e.g., by data augmentation techniques) version of a certain dataset. It is straightforward to check that Assumption 4 implies Assumption 3 with τ = λ and ρ = 1.
We conclude this section by defining the function classes we study in this paper. Definition 1 (Function classes). We consider two classes of global objective functions F , also taking into account their local objectives Fm and local components fmi . We assume throughout that f m i are differentiable and F is bounded from below.
Fobj(L, µ, ν, τ, ρ) := { F | F is µ-PŁ; fmi are L-smooth; F , Fm, fmi satisfy Assumptions 2 & 3 } ,
Fcmp(L, µ, ν, λ) := { F | F is µ-PŁ; fmi are L-smooth; F , Fm, fmi satisfy Assumptions 2 & 4 } .
Notice that Fobj(L, µ, ν, τ, ρ) ⊃ Fcmp(L, µ, ν, τ) for any ρ ≥ 1. We only make the PŁ assumption on the global objective F , not on the local objectives Fm nor on the local components fmi . Using L and µ, we define the condition number κ := L/µ ≥ 1.

heading: 3 CONVERGENCE ANALYSIS OF MINIBATCH AND LOCAL RR

heading: 3.1 UPPER BOUND FOR MINIBATCH RR
text: We first begin with the convergence result for minibatch RR on Fobj(L, µ, ν, τ, ρ), which exhibits a faster large-epoch rate compared to the single-machine setting. For upper bounds, we use Õ(·) to hide universal constants and logarithmic factors of 1δ , M , N , K, and B. Theorem 1 (Upper bound for minibatch RR). Suppose that minibatch RR has parameters satisfying Assumption 1. For any F ∈ Fobj(L, µ, ν, τ, ρ), consider running the algorithm using step-size η = B log(MNK
2) µNK for epochs K ≥ 6κ log(MNK 2). Then, with probability at least 1− δ,
F (xK,NB )− F ∗ ≤ F (x0)− F
∗
MNK2 + Õ
( L2
µ3 ν2 MNK2
) . (3)
Proof. The proof is in Appendix D.2. The key challenge in the convergence analysis of our shuffling-based method stems from the indices sampled within an epoch being dependent on each other. For example, if fm1 is accessed already, then the index i = 1 will not be used in later iterations of the epoch; this dependence significantly complicates the analysis. Our approach starts with realizing that for any permutation σ, ∑N i=1 f m σ(i) = NF
m. We decompose gradients ∇fmσmk (j)(xk,i−1) (see Line 10 of Algorithm 2) into ∇fmσmk (j)(xk,0) plus noise, then aggregate all updates over an epoch to get “one big step of GD plus noise”: xk+1,0 = xk,0 − ηN∇F (xk,0) + η2rk. We bound the noise rk using Lemma 8 (Appendix D.6), which is our extension of the Hoeffding-Serfling inequality to the mean ofM independent without-replacement sums of vectors; the lemma might be of independent interest too. Lemma 8 shows that averaging accumulated gradients over M machines reduces variance by M , which leads to the reduction by a factor of M in the bound (3).
Theorem 1 shows that for large enough epochs K & κ, minibatch RR converges at a rate of Õ( L 2ν2
µ3MNK2 ), with high probability. Compared to the large-epoch rate Õ( L2ν2
µ3NK2 ) of singlemachine RR (e.g., Ahn et al. (2020)), we see an additional factor M in the denominator, which highlights the advantage of multiple machines. If we compare against the with-replacement counterpart, it is known that for strongly convex and smooth F , the optimal convergence rate of minibatch SGD is Θ( ν 2
µMNK ), 5 which is worse than our bound (3) if K & κ2. Also notable is that the con-
vergence rate does not depend on the heterogeneity constants (i.e., τ and ρ from Assumption 3) of the local objective functions. This observation that minibatch RR is “immune” to heterogeneity is consistent with minibatch SGD in the with-replacement setting (Woodworth et al., 2020b).
Epoch vs communication complexity. One might wonder why (3) does not have the batch size B. In (3), we wrote convergence rates in terms of epochs K, which captures the gradient computation complexity because the same number of gradients are evaluated in a single epoch regardless ofB. If we are interested in communication complexity instead, we can write (3) in terms of the number of communication rounds R := NKB and get a rate of Õ( L2ν2N µ3MB2R2 ). From these, we can also discuss the overall cost of the algorithm. If the cost of a communication round is cc, and the cost of local gradient computations over an epoch is ce, then the total cost to obtain an -accurate solution is
Cminibatch( ) = Õ ( ccν √ N
B √ M + ceν√ MN
) , (4)
omitting L and µ for simplicity. The total cost shows that there is essentially no harm increasing the batch sizeB in minibatch RR, as we can get more accurate estimates of true gradients asB becomes larger. In the next subsection, we will see that this is not the case in local RR.
What about K . κ? We remark that all upper bounds in this paper hold only for the “large-epoch” regime, where K & κ. Such requirements are common in the literature of without-replacement SGD (Haochen & Sra, 2019; Nagaraj et al., 2019; Rajput et al., 2020; Ahn et al., 2020), and there is a recent result (Safran & Shamir, 2021) suggesting that faster convergence of without-replacement SGD may not be possible in the K . κ regime. We defer a more detailed discussion on this regime to Section 4, after Theorem 3.

heading: 3.2 UPPER BOUND FOR LOCAL RR
text: Next, we are interested in how fast local RR can converge, what is the optimal batch size B, and whether local RR can be as fast as minibatch RR.
Theorem 2 (Upper bound for local RR). Suppose that local RR has parameters satisfying Assumption 1. For any F ∈ Fobj(L, µ, ν, τ, ρ), consider running the algorithm using step-size η = log(MNK
2) µNK for epochs K ≥ 7ρκ log(MNK 2). Then, with probability at least 1− δ,
F (yK,NB )− F ∗ ≤ F (y0)− F
∗
MNK2 + Õ
( L2
µ3
( ν2
MNK2 +
ν2B
N2K2 +
τ2B2
N2K2
)) . (5)
5The optimal rate for (with-replacement) SGD after R iterations is Θ( ν 2
µR ) (see e.g., Rakhlin et al. (2012)).
With-replacement minibatching reduces the variance ν2 to ν 2
MB , andR = NK B . However, achieving the optimal
rate for last iterates typically requires carefully designed step-size schemes (Jain et al., 2019).
Proof. The proof is in Appendix D.3. We take the same “big GD step plus noise” approach as in Theorem 1; however, due to local updates, bounding the noise is much more involved. In the proof, we obtain the epoch update yk+1,0 = yk,0− ηN∇F (xk,0) + η2rk,1 + η2rk,2− η3rk,3, where rk,1 and rk,3 contain errors introduced by local updates. Noise from local updates accumulates over B iterations, which cannot be remedied by averaging over M machines. They result in two additional terms in the rate (5), one from intra-machine variance and the other from heterogeneity.

heading: 3.2.1 DISCUSSION OF THEOREM 2
text: Let us compare our high-probability bound (5) with existing in-expectation bounds. For strongly convex F , the corresponding last-iterate bound of local SGD is Õ( Lν 2
µ2MNK + L2ν2B µ3N2K2 + L2τ2B2 µ3N2K2 ) 6
(Khaled et al., 2020; Spiridonoff et al., 2020; Qu et al., 2020). Notice that (5) is better than this with-replacement bound when K & κ. For average iterates, there are known bounds Õ( ν 2
µMNK + Lν2B
µ2N2K2 + Lτ2B2 µ2N2K2 ) 6 (Koloskova et al., 2020; Woodworth et al., 2020b) which are smaller than the last-iterate bound by a factor of κ. It is unclear if averaging iterates could improve our rate, because most such analyses exploit Jensen’s inequality, which we cannot use for nonconvex F .
Dependence on τ and ρ. Out of the two heterogeneity constants τ and ρ (Assumption 3), ρ does not appear in (5), and it only affects the epoch requirement K & ρκ. Consider the case τ = 0 and ρ > 1, which is heterogeneous but in the “interpolation regime,” because ∇Fm(x) = 0 whenever ∇F (x) = 0. In such a case, the rate (5) is equal to the homogeneous case. Using B = Θ(N) is no better than single-machine. A close look at Theorem 2 reveals a rather surprising fact. Even in the homogeneous case (τ = 0), if we choose B = Θ(N), then local RR converges at the rate of Õ( 1NK2 ): the same rate as the single-machine RR! In Section 4, we show that this observation is not due to a suboptimal analysis; the rate Õ( 1NK2 ) is tight for B = Θ(N). Trade-off in the choice of B. As done for Theorem 1, we can compute from (5) that the total cost of local RR for -accuracy is (omitting L and µ for simplicity)
Clocal( ) = Õ ( cc ( ν √ N
B √ M + ν√ B + τ√
) + ce ( ν√ MN + ν √ B N √ + τB N √ )) . (6)
Note that for local RR, there exists a trade-off between communication and epoch complexity in the choice of B. If B is too small, this reduces the number of epochs required but increases communication costs. On the other hand, if B is too large, this reduces communication rounds but errors that accumulate in local updates get severer, resulting in the need for more epochs. Hence, the optimal choice of B must balance the two complexity measures. The existence of this trade-off is indeed different from minibatch RR where larger B always reduces the total cost Cminibatch( ).
When can local RR match minibatch RR? Comparing the convergence rates (3) and (5), we can identify some regimes in which local RR converges as fast as minibatch RR. In a nutshell, if machines are not too heterogeneous and communication happens frequently, then local RR can have the same upper bound as minibatch RR. For example, if B is chosen to be a constant, M . N , and τ . ν √ N/M , then the Õ( L 2ν2
µ3MNK2 ) term in (5) becomes the dominating factor and hence matches (3). Another example of such a regime is when B . NM and τ . ν √ M/N . Note that this comparison assumes that the same values of B are chosen for both algorithms. Also, such “frequent communication” regimes are favorable if the communication cost cc is small.
Can local RR ever beat minibatch RR? The upper bounds (3) and (5) indicate that local RR is always no better than minibatch RR, at least for the function class Fobj(L, µ, ν, τ, ρ). This is in fact consistent with Woodworth et al. (2020a;b), because the authors identify a regime where local SGD performs better than minibatch SGD for convex objective functions, but fail to do so for strongly convex functions. However, as was also pointed out in Woodworth et al. (2020a), there is a simple extreme scenario in which local RR can be faster: when ν ≈ τ ≈ 0 and ρ ≈ 1. In this case, we have fmi ≈ F for all m and i, so local RR corresponds to NK steps of GD, whereas minibatch RR corresponds to NKB steps of GD. Clearly, local RR will converge faster, exploiting the advantage of more updates. Finding out other such regimes is an important future direction.
6Due to differences in assumptions, many existing rates cannot be compared directly. These rates are the ones we consider “comparable” to our bound. See Appendix A for more detailed comparisons.

heading: 4 MATCHING LOWER BOUNDS
text: In Section 3, we presented large-epoch upper bounds (i.e., for K & κ) for constant step-size minibatch and local RR. In this section, we prove matching lower bounds to show that the upper bounds are tight, in all factors except L and µ. We use Ω(·) to hide universal constants in lower bounds.

heading: 4.1 LOWER BOUND FOR MINIBATCH RR
text: Theorem 3 (Lower bound for minibatch RR). Suppose that minibatch RR has parameters satisfying Assumption 1. Additionally, assume that N is a multiple of 2. Then, there exist large enough constants c1, c2 > 0 such that the following holds: For L and µ satisfying κ = Lµ ≥ c1, there exists a function F ∈ Fcmp(L, µ, ν, 0) such that for any constant step-size η,
E [ F (xK,NK )− F ∗ ] = Ω ( ν2 µMNK ) if K < c2κ, Ω ( ν2
µMNK2
) if K ≥ c2κ.
(7)
Proof. We prove Theorem 3 in Appendix E. The proof is an extension of Rajput et al. (2020); Safran & Shamir (2020; 2021) to minibatch RR. We will sketch some key intuitions after Theorem 4.
First notice that the function F is fromFcmp(L, µ, ν, 0), where all the machines are component-wise homogeneous. As seen in Definition 1, Fcmp(L, µ, ν, 0) ⊂ Fobj(L, µ, ν, τ, ρ) for any τ ≥ 0 and ρ ≥ 1, so Theorem 3 provides a lower bound for Fcmp(·) and Fobj(·), with arbitrary heterogeneity constants. We assume that N is even because we construct functions g1 and g2 such that fmi := g1 if i ≤ N2 , and f m i := g2 if i > N 2 . One can remove this assumption by using a zero function when N is odd (see e.g., Safran & Shamir (2020)). It is rather unsatisfactory that our theorem requires large enough constants c1 and c2; we believe a tighter analysis can relax this restriction.
Theorem 3 proves lower bounds for two different regimes: K & κ and K . κ. In the large-epoch regime (K & κ), we can observe that the lower bound Ω( ν 2
µMNK2 ) matches the upper bound (3) in Theorem 1, modulo a factor of κ2. Tightening the κ2 gap between upper and lower bounds is left for future work. In the small-epoch regime (K . κ), we observe that the lower bound Ω( ν 2
µMNK )
exactly matches the convergence rate of (with-replacement) minibatch SGD; hence, the lower bound implies that minibatch RR has no hope for faster convergence than minibatch SGD, at least in the constant step-size and small-epoch regime. This observation is in line with Safran & Shamir (2021).
Upper bounds for K . κ? Even for single-machine RR (M = 1), proving an upper bound that matches the small-epoch lower bound Ω( ν 2
µNK ) still remains a challenge. Nagaraj et al. (2019, The-
orem 2) prove an upper bound for non-quadratic strongly convex functions that matches Ω( ν 2
µNK ) if NK & κ2; however, they use suffix averaging, so it is not directly comparable to Theorem 3 which considers last iterates. Safran & Shamir (2021) prove upper bounds for quadratic strongly convex functions, but assume that their Hessian matrices commute. For noncommutative cases, proving a small-epoch upper bound seems to require some form of matrix AM-GM inequalities, whose availability is an open problem (Recht & Ré, 2012; Lai & Lim, 2020; De Sa, 2020; Yun et al., 2021). Remark 1 (Strong convexity in construction). We note that all lower bounds in this paper are constructed with strongly convex functions, a stronger assumption than PŁ functions (2). Thus, our lower bounds are also applicable to strong convexity counterparts of Fobj(·) and Fcmp(·).

heading: 4.2 LOWER BOUNDS FOR LOCAL RR
text: In this subsection, we present lower bounds for local RR. We prove two bounds that correspond to homogeneous and heterogeneous cases. By combining the two bounds, we get a lower bound that matches our upper bound (5) in Theorem 2 up to a factor of κ2. Theorem 4 (Lower bound for local RR: homogeneous case). Suppose that local RR has parameters satisfying Assumption 1. Additionally, assume that B is a multiple of 4. Then, there exist large enough constants c3, c4 > 0 such that the following holds: For L and µ satisfying κ = Lµ ≥ c3, there exists a function F ∈ Fcmp(L, µ, ν, 0) such that for any constant step-size η,
E [ F (yK,NK )− F ∗ ] = Ω ( ν2 µMNK ) if K < max { c4κ, MB N } , Ω ( ν2
µMNK2 + ν2B µN2K2
) if K ≥ max { c4κ, MB N } .
(8)
Proof. The proof is in Appendix G. For the large-epoch lower bounds in Theorems 3 and 4, we use “skewed” quadratics fmi (x) = (L1x≤0 + µ1x>0) x2 2 + ziνx, where zi = +1 if i ≤ N 2 and zi = −1 otherwise. For x ≈ 0, the imbalance results in a “drift” towards positive x, whose strength is approximately proportional to the absolute value of partial sums of random permutations over N2 +1’s and N2 −1’s. By averaging the sums over M machines (minibatch RR), their absolute values shrink by 1√
M ; in contrast, if each machine makes local updates (local RR), the magnitude of the
drift cannot be reduced with M , because we average after local iterates already have taken B “big” steps. The proof uses techniques from Rajput et al. (2020).
Proposition 5 (Lower bound for local RR: heterogeneous case). Suppose that local RR has parameters satisfying Assumption 1. Additionally, assume that B is a multiple of 2 and κ = Lµ ≥ 2. Then, there exists a function F ∈ Fobj(L, µ, 0, τ, 1) such that for any constant step-size η,
E [ F (yK,NK )− F ∗ ] = Ω ( τ2B2
µN2K2
) . (9)
Proof. We note that Proposition 5 is almost identical to Theorem II of Karimireddy et al. (2020); however, we provide a proof specific to our algorithm in Appendix I.
Theorem 4 constructs a component-wise homogeneous function from Fcmp(L, µ, ν, 0) and Proposition 5 constructs a heterogeneous function from Fobj(L, µ, 0, τ, 1). Since Fcmp(L, µ, ν, 0) ∪ Fobj(L, µ, 0, τ, 1) ⊂ Fobj(L, µ, ν, τ, ρ) for any ρ ≥ 1, combining (8) and (9) for the K ≥ max{c4κ, MBN } case gives a lower bound Ω(max{ ν2 µMNK2 + ν2B µN2K2 , τ2B2
µN2K2 }) that matches the large-epoch upper bound (5) in Theorem 2, up to a factor of κ2. When κN & MB, c4κ becomes the dominating term in the max, in which case the threshold in (8) is Θ(κ). Tightening the κ2 gap as well as removing additional requirements such as κ ≥ c3 and κN &MB are left for future work. Using B = Θ(N) does not help, indeed. In Section 3.2.1, we observed that if B = Θ(N), then even in the homogeneous case (τ = 0), local RR converges at the rate of Õ( 1NK2 ). This is the same rate as single-machine RR, meaning the efforts byM−1 machines become meaningless. Our lower bound (8) shows that Õ( 1NK2 ) is in fact the best we can hope for (treating L and µ as constants). In order to make the best use of M machines, B should be smaller than Θ(N), as suggested in Section 3.2.1. In an existing work, Mishchenko et al. (2021) consider local RR with B = N as a special case of a proximal algorithm. In Theorem 8 of Mishchenko et al. (2021), the authors claim “the convergence bound improves with the number of devices involved” because the bound has a factor of M in the denominator. However, at least under our assumption, this is not the case; if we apply our Assumption 3 to upper-bound their σ∗, the term “Nσ2∗” in the numerator grows linearly with M . Hence, our bounds do not contradict Mishchenko et al. (2021); see Appendix A for details. Remark 2 (Small-epoch bound is likely loose). We note that while we focused on deriving a matching large-epoch lower bound, we did not try hard to tighten the small-epoch lower bound. Our small-epoch lower bound in (8) misses a term (such as ν
2B µN2K2 ) that corresponds to the error from
local updates. We leave investigations on small-epoch lower and upper bounds for future work.

heading: 5 SYNCHRONIZED SHUFFLING: HOW TO BYPASS LOWER BOUNDS
text: Recall from the total complexity of minibatch RR (4) that the total cost shrinks with a factor of 1√ M .
Using M machines, we are only getting a √ M -factor speedup. Ideally, we hope to see a linear speedup, i.e., cost inverse proportional to M . Hence, Theorem 1 falls short of achieving this goal, and our lower bound in Theorem 3 confirms that linear speedup is indeed impossible.
In this section, we show that the desired linear speedup is possible, at least in some special cases. We consider the component-wise near-homogeneous case (i.e., Assumption 4 with small λ) and discuss how a simple modification to minibatch and local RR can let us “break” the lower bounds and achieve linear speedup. This comes at a cost of broadcasting permutations: at the beginning of the k-th epoch, the server samples σ ∼ Unif(SN ) and π ∼ Unif(SM ), and broadcasts them to the machines. Then, local machines choose their permutations σmk to be shifted versions of σ, 7
i.e., σmk (i) := σ (( i+ NM π(m) ) mod N ) . We call this trick synchronized shuffling, denoted as
7We assume for simplicity that M divides N .
SYNCSHUF. Please revisit Algorithms 1 and 2 for the precise descriptions of the modified algorithms local RR with SYNCSHUF and minibatch RR with SYNCSHUF, respectively.
The intuition why this should help is simple. In the proof of RR, we aggregate the component gradients over an epoch (i.e., N iterations) to write it as a full gradient plus noise. If we are in the component-wise homogeneous setting and permutations are synchronized, then instead of aggregating N component gradients on a single machine, we can aggregate NM component gradients on M machines to get a full gradient. This allows us to reduce the “noise” from without-replacement sampling. We emphasize here that we do not necessarily set B = NM to get a full gradient every time; our analysis works for arbitrary B and M , as long as both divide N . See Appendix B for a detailed illustration of SYNCSHUF; also, see Appendix C for experiments showing its effectiveness.
The idea of synchronized shuffling is similar to approaches in distributed learning that shuffle and partition datasets and distribute them to local machines (see e.g., Lee et al. (2017); Meng et al. (2017)). In contrast, we do not communicate data, but communicate how to permute datasets stored in local machines. Meng et al. (2017, Theorem 3.3) provide an analysis for a distributed method similar to minibatch RR, but fail to show convergence to global minima in strongly convex cases. We also note that an independent concurrent result (Szlendak et al., 2021) uses the same idea as SYNCSHUF to build compressors for communication-efficient distributed optimization.

heading: 5.1 UPPER BOUNDS FOR MINIBATCH AND LOCAL RR WITH SYNCSHUF
text: With SYNCSHUF, we can show that the M ’s appearing in the convergence rates ((3) and (5)) in Theorems 1 and 2 can be replaced withM2, for a more stringent function classFcmp(·) that requires bounded component-wise inter-machine deviation (Assumption 4). Theorem 6 (Upper bound for minibatch RR with SYNCSHUF). Suppose that minibatch RR with SYNCSHUF has parameters satisfying Assumption 1. Additionally assume that M divides N . For any F ∈ Fcmp(L, µ, ν, λ), consider running the algorithm using step-size η = B log(M 2NK2) µNK for epochs K ≥ 6κ log(M2NK2). Then, with probability at least 1− δ,
F (xK,NB )− F ∗ ≤ F (x0)− F
∗
M2NK2 + Õ
( L2
µ3
( ν2
M2NK2 +
λ2
MK2
)) . (10)
The proof of Theorem 6 is presented in Appendix D.4. One can check that if the componentwise deviation constant λ satisfies λ . ν√
MN (i.e., near-homogeneous), then the rate (10) becomes
Õ( 1M2NK2 ). It is then easy to confirm that M machines reduce total costs by 1 M —a linear speedup.
A similar speedup can be shown for local RR. In Appendix D.5, we prove that Theorem 7 (Upper bound for local RR with SYNCSHUF). Suppose that local RR with SYNCSHUF has parameters satisfying Assumption 1. Additionally assume that M divides N . For any F ∈ Fcmp(L, µ, ν, λ), consider running the algorithm with step-size η = log(M 2NK2) µNK for epochs K ≥ 7κ log(M2NK2). Then, with probability at least 1− δ,
F (yK,NB )− F ∗ ≤ F (y0)− F
∗
M2NK2 + Õ
( L2
µ3
( ν2
M2NK2 +
ν2B
N2K2 +
λ2B2
N2K2 +
λ2
MK2
)) . (11)
We can similarly check that if B . NM2 and λ . ν√ MN , i.e., frequent communication and nearhomogeneity, then the Õ( 1M2NK2 ) term dominates in (11), and hence gives a linear speedup that matches the best rate of minibatch RR with SYNCSHUF (10). Nevertheless, we note again that for local RR, such a small B is favorable only when the communication cost cc is small (recall (6)).

heading: 6 CONCLUSION
text: We studied convergence bounds for local RR and minibatch RR, which are the practical withoutreplacement versions of local and minibatch SGD studied in the theory literature. For smooth functions satisfying the Polyak-Łojasiewicz condition, we showed large-epoch convergence bounds for minibatch and local RR that are faster than their with-replacement counterparts. We also proved matching lower bounds showing that our convergence analysis is tight. We also proposed a simple modification called synchronized shuffling that leads to convergence rates faster than our lower bounds in near-homogeneous settings. Immediate future research directions include extension to small-epoch regimes, as well as to general convex and nonconvex functions.

heading: ETHICS STATEMENT
text: This paper develops theoretical guarantees for popular distributed stochastic optimization algorithms. Therefore, the authors do not see any particular concerns related to its ethical aspects or future societal consequences.

heading: REPRODUCIBILITY STATEMENT
text: This paper is a theoretical work, without any experimental results. Definitions and assumptions are provided in Section 2. Our theoretical contributions as well as some additionally required assumptions are clearly stated in Sections 3, 4, and 5. Complete proofs of all the theorems are provided in the appendix.

heading: CONTENTS

heading: 1 Introduction 1
text: 1.1 Our contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

heading: 2 Problem setup 2

heading: 3 Convergence analysis of minibatch and local RR 4
text: 3.1 Upper bound for minibatch RR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.2 Upper bound for local RR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

heading: 4 Matching lower bounds 7
text: 4.1 Lower bound for minibatch RR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
4.2 Lower bounds for local RR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

heading: 5 Synchronized shuffling: how to bypass lower bounds 8
text: 5.1 Upper bounds for minibatch and local RR with SYNCSHUF . . . . . . . . . . . . 9

heading: 6 Conclusion 9

heading: A Comparisons with assumptions and rates in existing results 14

heading: B More detailed illustration of synchronized shuffling 18

heading: C Experimental results 19
text: C.1 SYNCSHUF improves convergence of minibatch/local RR . . . . . . . . . . . . . 20
C.2 Local RR becomes closer to single-machine RR as B → N . . . . . . . . . . . . 21
C.3 With- vs. without-replacement sampling . . . . . . . . . . . . . . . . . . . . . . . 22

heading: D Proofs of upper bounds 23
text: D.1 Proof outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
D.2 Proof of upper bound for minibatch RR (Theorem 1) . . . . . . . . . . . . . . . . 24
D.3 Proof of upper bound for local RR (Theorem 2) . . . . . . . . . . . . . . . . . . . 29
D.4 Proof of upper bound for minibatch RR with SYNCSHUF (Theorem 6) . . . . . . . 35
D.5 Proof of upper bound for local RR with SYNCSHUF (Theorem 7) . . . . . . . . . 39
D.6 A generalized vector-valued Hoeffding-Serfling inequality . . . . . . . . . . . . . 42
D.7 How can we avoid uniform bounds over Rd in our assumptions? . . . . . . . . . . 44

heading: E Proof of lower bound for minibatch RR (Theorem 3) 46
text: E.1 Lower bound for η ≤ BµNK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
E.2 Lower bound for η ≥ BµNK and η ≤ B 513LN . . . . . . . . . . . . . . . . . . . . . 47
E.3 Lower bound for η ≥ BµNK and η ≥ B 513LN . . . . . . . . . . . . . . . . . . . . . 50

heading: F Proofs of helper lemmas for Appendix E 52
text: F.1 Proof of Lemma 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
F.2 Proof of Lemma 12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
F.3 Proof of Lemma 13 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
F.4 Proof of Lemma 14 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
F.5 Proof of Lemma 15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

heading: G Proof of lower bound for local RR: homogeneous case (Theorem 4) 59
text: G.1 Lower bound for η ≤ 1µNK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
G.2 Lower bound for η ≥ 1µNK and η ≤ 1 1025LN . . . . . . . . . . . . . . . . . . . . . 61
G.3 Lower bound for η ≥ 1µNK and η ≥ 1 1025LN . . . . . . . . . . . . . . . . . . . . . 64

heading: H Proofs of helper lemmas for Appendix G 66
text: H.1 Proof of Lemma 16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
H.2 Proof of Lemma 17 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
H.3 Proof of Lemma 18 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
H.4 Proof of Lemma 19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

heading: I Proof of lower bound for local RR: heterogeneous case (Proposition 5) 73
text: I.1 Lower bound for 18µNK ≤ η ≤ 1 8µB and 1 8µB ≤ η ≤ 1 µ . . . . . . . . . . . . . . . 73
I.2 Lower bound for η ≤ 18µNK and η ≥ 1 µ . . . . . . . . . . . . . . . . . . . . . . . 75

heading: A COMPARISONS WITH ASSUMPTIONS AND RATES IN EXISTING RESULTS
text: In this section, we compare our assumptions and convergence bounds against other existing results mentioned in the main text. Most existing results that study independent and unbiased gradient estimates state their assumptions in terms of the expectation over the randomness in the estimate; for such assumptions, we adapt them to our finite sum setting in order to make for easier comparison.
Heterogeneity assumptions. We start by discussing our definition of objective-wise heterogeneity (Assumption 3), namely that there exist τ ≥ 0 and ρ ≥ 1 such that
1
M M∑ m=1 ‖∇Fm(x)‖ ≤ τ + ρ ‖∇F (x)‖ , for all x ∈ Rd. (12)
Perhaps the most relevant to this assumption is the (G,B)-BGD assumption that appears in Karimireddy et al. (2020): For all x ∈ Rd,
1
M M∑ m=1 ‖∇Fm(x)‖2 ≤ G2 +B2 ‖∇F (x)‖2 . (13)
Note that thanks to Jensen’s inequality and a2 + b2 ≤ (a+ b)2 for a, b ≥ 0, (13) implies( 1
M M∑ m=1 ‖∇Fm(x)‖ )2 ≤ 1 M M∑ m=1 ‖∇Fm(x)‖2 ≤ G2 +B2 ‖∇F (x)‖2 ≤ (G+B ‖∇F (x)‖)2,
and hence (12) with τ = G and ρ = B. Therefore, our Assumption 3 is weaker than the (G,B)BGD assumption. Several papers (Haddadpour & Mahdavi, 2019; Li et al., 2020b) use the same
assumption (13), with G = 0; therefore, Assumption 3 also subsumes the heterogeneity assumption posed in these papers. Note that G = 0 implies that, the minima for F are also the minima for Fm, for every m, and hence G = 0 results in a weak form of heterogeneity.
Some papers (Yu et al., 2019; Li et al., 2020c; Qu et al., 2020) assume bounded local gradients: for all m ∈ [M ] and x ∈ Rd,
1
N N∑ i=1 ‖∇fmi (x)‖ 2 ≤ G2, (14)
and this in fact implies the (G, 0)-BGD assumption (13). To see why, from Jensen’s inequality
‖∇Fm(x)‖2 := ∥∥∥∥∥ 1N N∑ i=1 ∇fmi (x) ∥∥∥∥∥ 2 ≤ 1 N N∑ i=1 ‖∇fmi (x)‖ 2 ≤ G2.
Therefore, (14) is a stronger assumption for objective-wise heterogeneity than Assumption 3.
In Theorem 3 of Woodworth et al. (2020b), the authors use the following assumption on heterogeneity: for all x ∈ Rd,
1
M M∑ m=1 ‖∇Fm(x)−∇F (x)‖2 ≤ ζ̄2. (15)
Noting that 1M ∑M m=1 ‖∇Fm(x)−∇F (x)‖ 2 = 1M ∑M m=1 ‖∇Fm(x)‖
2−‖∇F (x)‖2, we can see that (15) implies (13) and hence Assumption 2 (12), with τ = ζ̄ and ρ = 1.
Indeed, there are also some results that make weaker heterogeneity assumptions than ours (12), by requiring bounded deviation only at the global optimum x∗. Given the global optimum x∗ of F , Koloskova et al. (2020) define
ζ̄2∗ := 1
M M∑ m=1 ‖∇Fm(x∗)‖2 , (16)
and use this constant in their bounds. Khaled et al. (2020) also define a similar quantity that can capture heterogeneity, but does not provide a result on strongly convex functions in the heterogeneous setting. While assuming bounded ζ̄2∗ (16) is weaker than Assumption 3 in the sense that only a bound at x∗ is required, we note that these assumptions cannot be applied easily in nonconvex settings; in fact, for nonconvex (but not necessarily PŁ) functions, Koloskova et al. (2020) also use (13) as their heterogeneity assumption.
Intra-machine variance assumptions. We next consider our notion of intra-machine deviation (Assumption 2), namely that there exists ν ≥ 0 such that for all m ∈ [M ] and i ∈ [N ],
‖∇fmi (x)−∇Fm(x)‖ ≤ ν, for all x ∈ Rd. (17)
In many existing results using independent and unbiased gradient estimates (Karimireddy et al., 2020; Yu et al., 2019; Li et al., 2020c; Qu et al., 2020; Woodworth et al., 2020a;b; 2021), the bounded variance assumption is adopted: For all m ∈ [M ] and x ∈ Rd,
1
N N∑ i=1 ‖∇fmi (x)−∇Fm(x)‖ 2 ≤ σ2. (18)
We note that the bounded local gradients assumption (14) also implies (18), by
1
N N∑ i=1 ‖∇fmi (x)−∇Fm(x)‖ 2 ≤ 1 N N∑ i=1 ‖∇fmi (x)−∇Fm(x)‖ 2 + ‖∇Fm(x)‖2
= 1
N N∑ i=1 ‖∇fmi (x)‖ 2 ≤ G2. (19)
Some other papers (Haddadpour et al., 2019; Haddadpour & Mahdavi, 2019; Spiridonoff et al., 2020) consider a generalized version of (18), namely that for all m ∈ [M ] and x ∈ Rd,
1
N N∑ i=1 ‖∇fmi (x)−∇Fm(x)‖ 2 ≤ c ‖∇Fm(x)‖2 + σ2, (20)
for c, σ ≥ 0. There are other papers (Koloskova et al., 2020; Khaled et al., 2020) that use intramachine variance at the global minimum x∗ in their bounds. Koloskova et al. (2020) define
σ̄2∗ := 1
MN M∑ m=1 N∑ i=1 ‖∇fmi (x∗)−∇Fm(x∗)‖ 2 , (21)
for the global minimum x∗, and use it in their strong convexity and convexity bounds. Khaled et al. (2020) define
σ2opt := 1
MN M∑ m=1 N∑ i=1 ‖∇fmi (x∗)‖ 2 , (22)
which is used in their bound on strongly convex functions in homogeneous cases. Note that for homogeneous cases,∇Fm(x∗) = ∇F (x∗) = 0, so σ̄2∗ = σ2opt. We note that in contrast to our discussion on inter-machine deviation (Assumption 3), the intramachine variance assumptions in the existing literature are weaker than our Assumption 2. However, we utilize our stronger assumption to prove our high-probability upper bounds, which is a departure from in-expectation bounds in the literature.
Existing upper bounds on local SGD. In the discussion after Theorem 2, we mentioned some recent upper bounds on (with-replacement) local SGD. We make more detailed comparisons here. For the reader’s convenience, we restate our theorem on local RR below. Theorem 2 (Upper bound for local RR). Suppose that local RR has parameters satisfying Assumption 1. For any F ∈ Fobj(L, µ, ν, τ, ρ), consider running the algorithm using step-size η = log(MNK
2) µNK for epochs K ≥ 7ρκ log(MNK 2). Then, with probability at least 1− δ,
F (yK,NB )− F ∗ ≤ F (y0)− F
∗
MNK2 + Õ
( L2
µ3
( ν2
MNK2 +
ν2B
N2K2 +
τ2B2
N2K2
)) . (5)
We start with last-iterate bounds in homogeneous cases. Theorem 3 and Corollary 3 of Khaled et al. (2020) consider local SGD on µ-strongly convex and L-smooth F . Khaled et al. (2020) allow variable synchronization intervals, where the maximum interval is upper-bounded by B. In this setting, Corollary 3 of Khaled et al. (2020) shows that local SGD after T total local update steps yield
E [ ‖x̄T − x∗‖2 ] = Õ ( ‖x̄0 − x∗‖2
T 2 + σ2opt µ2MT + Lσ2opt(B − 1) µ3T 2
) , (23)
where x̄i is the average over all M machines’ i-th local iterates, and σopt is from (22). Noting that T corresponds to NK in local SGD and σopt corresponds to ν in Assumption 2, (23) is comparable to an upper bound8
E [F (x̄T )− F ∗] = Õ ( Lν2
µ2MNK +
L2ν2B
µ3N2K2
) . (24)
Here, we do not compare F (y0)−F∗MNK2 and L‖x̄0−x∗‖2 T 2 because in Theorem 2, the term F (y0)−F∗ MNK2 can be made arbitrarily “small” (e.g., F (y0)−F∗ (MNK)l for any l ∈ N) by changing the log factor in η.
A similar homogeneous, strongly convex, and smooth setting is considered in Spiridonoff et al. (2020), under a intra-machine variance assumption defined in (20). Theorem 1 of Spiridonoff et al. (2020) proves general theorem statement for arbitrary synchronization intervals. If we specialize to constant interval B, Corollary 1 of Spiridonoff et al. (2020) gives
E [F (x̄T )− F ∗] = Õ ( β2(F (x̄0)− F ∗)
T 2 +
Lσ2 µ2MT + L2σ2(B − 1) µ3T 2
) , (25)
8Note that an additional factor L is due to conversion from squared distance to function value.
where β ≥ 2κ2 is a constant defined to choose algorithm parameters such as the step-size. Noting again that T corresponds to NK and σ in (20) is comparable to ν in Assumption 2, (25) also translates to (24).
The next last-iterate bound we compare against is Qu et al. (2020). Theorem 1 of Qu et al. (2020) uses bounded intra-machine variance assumption (18) and bounded gradient assumption (14). Specializing Theorem 1 of Qu et al. (2020) to full device participation and uniform weight (p1 = · · · = pN = 1 N ) case, their bound reads
E [F (x̄T )− F ∗] = Õ ( Lσ2
µ2MT + L2G2B2 µ3T 2
) . (26)
Recalling that T corresponds to NK, σ to ν in Assumption 2, and G to the heterogeneity constant τ in Assumption 3 and also ν (due to (19)), (26) translates to
E [F (x̄T )− F ∗] = Õ ( Lν2
µ2MNK + L2(ν2 + τ2)B2 µ3N2K2
) . (27)
Comparing the local SGD last-iterate bounds (24) and (27) against our local RR bound (5) in Theorem 2, we can see that the last iterate of local RR satisfies a smaller upper bound as soon as K ≥ κ. Admittedly, this is not a fully rigorous comparison given the differences in assumptions and types of bounds; nevertheless, we believe that the comparison at least provides some degree of evidence for faster convergence of local RR than local SGD. It is also interesting to see that the “error from local updates” terms match in with- and without-replacement bounds.
Next, we review existing average-iterate bounds mentioned in the main text, which are better than the last-iterate bounds. Koloskova et al. (2020) present a unifying framework for analyzing distributed optimization algorithms over networks, which can specialize to local SGD. For µ-strongly convex and L-smooth F , Theorem 2 of Koloskova et al. (2020) shows that
E [F (x̂)− F ∗] = Õ ( LB ‖x0 − x∗‖2 exp ( − µT LB ) + σ̄2∗ µMT + Lσ̄2∗B µ2T 2 + Lζ̄2∗B 2 µ2T 2 ) , (28)
where x̂ is some weighted average of iterates and ζ̄2∗ and σ̄ 2 ∗ are defined in (16) and (21), respectively. Noting that T corresponds to NK, ζ̄∗ to τ in Assumption 3, and σ̄∗ to ν in Assumption 2 (although our assumptions are stronger), we can see that the bound (28) for large enough T can be translated into
E [F (x̂)− F ∗] = Õ ( ν2
µMNK +
Lν2B
µ2N2K2 +
Lτ2B2
µ2N2K2
) . (29)
Notice that (29) is smaller than the last-iterate bounds (24) and (27) by a factor of κ. Similarly, Theorem 3 of Woodworth et al. (2020b) proves that for µ-strongly convex and L-smooth F ,
E [F (x̃)− F ∗] = Õ
( L2 ‖x0 − x∗‖2
LT + µT 2 + σ̄2∗ µMT + Lσ2B µ2T 2 + Lζ̄2B2 µ2T 2
) , (30)
for some weighted average of iterates x̃. Here, σ̄2∗, σ 2, and ζ̄2 are as defined in (21), (18), and (15), respectively. Substituting ν to its comparable constants σ̄∗ and σ, and τ to ζ̄, we can similarly check that (30) can be converted to (29).
Comparison to Mishchenko et al. (2021). In Mishchenko et al. (2021), the authors study a proximal algorithm referred to as Proximal Random Reshuffling (ProxRR), and obtain a distributed optimization algorithm called FedRR as a special case. If we set R ≡ 0 in FedRR, the algorithm then is equal to local RR with B = N , i.e., the one that synchronizes only after one entire epoch.
Assuming that all component functions fmi (x) are µ-strongly convex 9 and L-smooth, and objectivewise homogeneity F 1 = · · · = Fm = F , the authors obtain Theorem 8 (Mishchenko et al., 2021), which states that
E [ ‖yK − x∗‖2 ] ≤ (1− ηµ)NK ‖y0 − x∗‖2 +
η2LMNσ2opt µM , (31)
9This is in fact quite strong compared to this paper, because we only assume F to be PŁ, not Fm nor fmi .
where y0 and yK are the initialization and last iterate of the algorithm, and σ2opt was defined above in (22). The termMNσ2opt in (31) corresponds to the term “Nσ 2 ∗” as per the notation in Mishchenko et al. (2021). If we apply our Assumption 3 to bound σ2opt, we get σ 2 opt ≤ ν2, which reduces the last term in (31) to η 2Lν2N µ . If we substitute η = log(MNK2) µNK to the bound (31), we get
E [ ‖yK − x∗‖2 ] ≤ ‖y0 − x ∗‖2
MNK2 + Õ
( Lν2
µ3NK2
) , (32)
which translates to the same convergence rate on F (yK) − F ∗ as single-machine RR. For the heterogeneous setting, applying Lemma 3 of Mishchenko et al. (2021) to Theorem 2, we can obtain
E [ ‖yK − x∗‖2 ] ≤ (1− ηµ)NK ‖y0 − x∗‖2
+ 2η2L
µM M∑ m=1
( N2 ‖∇Fm(x∗)‖2 + 1
4 N∑ i=1 ‖∇fmi (x∗)−∇Fm(x∗)‖ 2
)
= (1− ηµ)NK ‖y0 − x∗‖2 + 2η2LN2
µM M∑ m=1 ‖∇Fm(x∗)‖2 + η 2Lσ̄2∗N 2µ , (33)
where σ̄2∗ was defined in (21). Note that in Lemma 3 (Mishchenko et al., 2021), the function “Fm” in the authors’ notation is equal to NFm in our notation. Recall that the (G,B)-BGD assumption (13) is “comparable” to Assumption 3 (12). If we apply (13) to the bound (33), we get 1 M ∑M m=1 ‖∇Fm(x∗)‖
2 ≤ G2. Similarly, if we apply Assumption 3, we get σ̄2∗ ≤ ν2. Substituting these upper bounds and η = log(MNK
2) µNK to (33) gives
E [ ‖yK − x∗‖2 ] ≤ ‖y0 − x ∗‖2
MNK2 + Õ
( LG2
µ3K2 +
Lν2
µ3NK2
) , (34)
and after translating this into a bound on function value, we get an upper bound Õ( L 2ν2 µ3NK2 + L2G2 µ3K2 ) which in fact matches our upper bound (5) in Theorem 2 when we set B = N .
The two upper bounds obtained for homogeneous (32) and heterogeneous (34) settings indicate that, at least under assumptions on intra- and inter-machine deviation such as ours, the claimed advantage that “the convergence bound improves with the number of devices involved” (Mishchenko et al., 2021) is not achievable. As our lower bound shows, one needs to choose B smaller than N in order to get the most out of parallelism. That being said, since Mishchenko et al. (2021) is free of uniform intra- and inter-machine deviation assumptions, there may still exist certain scenarios where multiple machines can speed up performance even with B = N .

heading: B MORE DETAILED ILLUSTRATION OF SYNCHRONIZED SHUFFLING
text: In this section, we provide a more detailed explanation on synchronized shuffling that we introduced in Section 5. For the illustration, let us consider the component-wise homogeneous case. Component-wise homogeneity means that all the machines have the same set of components: f1i = f 2 i = · · · = fMi =: fi for i ∈ [N ]. Hence, we have F = Fm for all m ∈ [M ] and our
goal is to minimize F = 1N ∑N i=1 fi.
In the proof of Theorem 1 (presented in Appendix D.2; see Appendix D.1 for a sketch), we add the component gradients over an epoch and then use the following key identity: for any permutation σ,
N∑ i=1 ∇fσ(i) = N∑ i=1 ∇fi = N∇F. (35)
We use this identity (35) to represent the per-epoch progress as “one big GD step plus noise.” For the rest of the proof we bound the “noise” term, and the key to bounding it is to upper bound the norm of summations of the following form, for i ∈ [N/B − 1] (see (41) and (43)):
1
M M∑ m=1 iB∑ j=1 ∇fσmk (j)(xk,0). (36)
To bound the norm, we decompose it into two terms using the triangle inequality∥∥∥∥∥∥ 1M M∑ m=1 iB∑ j=1 ∇fσmk (j)(xk,0) ∥∥∥∥∥∥≤ ∥∥∥∥∥∥ 1M M∑ m=1 iB∑ j=1 ∇fσmk (j)(xk,0)− iB∇F (xk,0) ∥∥∥∥∥∥+iB ‖∇F (xk,0)‖ , and we use concentration bounds (44) on the first term of the RHS, which gives a high-probability
upper bound on the RHS of the form Õ ( ν √
iB M
) + iB ‖∇F‖ (45).
The key to proving fast convergence is to make the upper bound above as small as possible. To make the bound (45) even smaller, we wish to be able to apply the identity (35) to the summation (36). However, under the standard way of choosing permutations σmk independently over machines, one cannot apply the identity because we do not sum over all j = 1, . . . , N . This limitation motivates our proposed technique synchronized shuffling, a manipulation on the choices of σmk that lets us prove even faster convergence.
Recall the definition of synchronized shuffling. At the beginning of the k-th epoch, the server samples σ ∼ Unif(SN ) and π ∼ Unif(SM ), and broadcasts them to the machines. Then, local machines choose their permutations σmk to be shifted versions of σ: σ m k (i) := σ (( i+ NM π(m) ) mod N ) .
Now setN = 6 andM = 3. Assume for simplicity that the permutation π ∈ S3 of machines satisfies π(m) = m for m ∈ [3].10 Suppose the server samples σ = (σ(1), σ(2), σ(3), σ(4), σ(5), σ(6)) and broadcasts it. Under synchronized shuffling, the local machines choose
σ1k = (σ(3), σ(4), σ(5), σ(6), σ(1), σ(2)),
σ2k = (σ(5), σ(6), σ(1), σ(2), σ(3), σ(4)),
σ3k = (σ(1), σ(2), σ(3), σ(4), σ(5), σ(6)).
One can see that each permutation is a shifted version of σ, with an offset that is a multiple of N M = 2.
Now consider adding ∇fσmk (i) over m = 1, 2, 3 and j = 1, 2 (in fact, any two consecutive j’s will do). By synchronized shuffling, we get a summation over all N = 6 component functions, which by (35) gives us the full gradient:
3∑ m=1 2∑ j=1 ∇fσmk (j) = 6∑ i=1 ∇fi = 6∇F.
The point here is that the permutation identity (35) can be applied to the summations (36) to further reduce their norm bounds. This is in contrast to sampling independent σmk ’s where one cannot apply (35). For this reason, synchronized shuffling significantly reduces the noise that comes from without-replacement sampling, thus resulting in faster convergence rates in (near-)homogeneous cases.

heading: C EXPERIMENTAL RESULTS
text: In this section, we present some simple numerical experiments that support our theoretical analysis. We evaluate the performance of the algorithms considered in this paper on the “hard instance” constructed in our lower bounds (Theorems 3 and 4).
Our hard instance F ∈ Fcmp(L, µ, ν, 0) is a function in the component-wise homogeneous setting, where all the machines have the same set of local component functions: f1i = f 2 i = · · · = fMi =: fi.
In the proofs of Theorems 3 and 4, we construct the global objective F (x) = 1N ∑N i=1 fi(x) as the following:
fi(x) := (L1x≤0 + µ1x>0) x2
2 + ziνx, zi = { +1 if 1 ≤ i ≤ N/2, −1 if N/2 < i ≤ N.
10In fact, permuting the machines by π is not required in the component-wise homogeneous setting (i.e., when λ = 0 in Assumption 4).
With this set of component functions, the global objective F (x) = (L1x≤0+µ1x>0)x 2
2 is µ-strongly convex and L-smooth with a unique global minimizer at x = 0.
We compare the performance of the algorithms on this problem instance, with L = 100, µ = 1, ν = 1, N = 768, and M = 16, while varying the choice of B ∈ {1, 4, 16, 64, 256} and K ∈ {1, 3, 5, 7, 10, 30, 50, 70, 100, 300, 500, 700, 1000}. For each value of B and K, we run the algorithms for K epochs (KN/B communication rounds for with-replacement algorithms) starting at x0 = 0 and return the values of F evaluated at the last iterates. Note that the algorithms are not deterministic, because the sampling/shuffling schemes are random. In order to account for randomness, for each combination of (algorithm, B,K) we execute 20 independent runs of the algorithm and plot the mean, first quartile, and third quartile of the final objective values.
In the subsequent subsections, we compare the following seven algorithms with constant step-sizes.
• Minibatch RR, η = B log(MNK 2)
µNK ;
• Local RR, η = log(MNK 2)
µNK ;
• Minibatch RR with SYNCSHUF, η = B log(MNK 2)
µNK ;
• Local RR with SYNCSHUF, η = log(MNK 2)
µNK ;
• Single-machine RR with minibatch size B, η = log(NK 2)
µNK ;
• With-replacement minibatch SGD, η = B log(MNK 2)
µNK ;
• With-replacement local SGD, η = log(MNK 2)
µNK .

heading: C.1 SYNCSHUF IMPROVES CONVERGENCE OF MINIBATCH/LOCAL RR
text: In Figure 1, we compare minibatch RR and local RR, with and without SYNCSHUF. Each plot in Figure 1 shows how the methods’ performance changes with K, for a fixed value B. Each point on the curve is the mean of the final objective function values over 20 independent runs of the corresponding algorithm with the specific B and K, and its error bar indicates the first and third quartiles.
Recall from our Theorems 1, 2, 6, and 7 that the four methods satisfy the following convergence bounds, in homogeneous settings (i.e., τ = λ = 0):
• Minibatch RR: Õ ( L2
µ3 ν2 MNK2
) ,
• Local RR: Õ ( L2
µ3
( ν2
MNK2 + ν2B N2K2
)) ,
• Minibatch RR with SYNCSHUF: Õ ( L2
µ3 ν2 M2NK2
) ,
• Local RR with SYNCSHUF: Õ ( L2
µ3
( ν2
M2NK2 + ν2B N2K2
)) .
In fact, if B = 1, local RR is identical to minibatch RR. Figure 1(a) confirms that this is indeed true, and also that the versions with SYNCSHUF outperforms the ones without SYNCSHUF. This corroborates the additional M factor speedup in our bounds. In Figure 1(b) and 1(c), we can see that as B increases, the performance of local RR with SYNCSHUF degrades and becomes closer to local RR without SYNCSHUF. This shows that the Õ ( L2
µ3 ν2B N2K2
) term starts to dominate. Also,
as we increase B further, in Figure 1(c) and 1(d) we see that local RR (without SYNCSHUF) also starts to degrade and its gap between minibatch RR becomes larger. Again, this means that the Õ ( L2
µ3 ν2B N2K2
) term becomes the dominant factor in the local RR bound. The performance of mini-
batch RR, with and without SYNCSHUF, stays relatively independent ofB. One thing to note is that for large values of B, the small-epoch behavior of minibatch RR looks rather unstable. The choice of step-size η = B log(MNK
2) µNK seems to cause overshooting when B is large and K is small. We
note that this does not contradict our convergence analysis because our theorems only characterize the large-epoch behavior (K above certain thresholds) of the algorithms. Perhaps in the small-epoch regime, our choice of η is not necessarily optimal and a smaller η is needed to prevent overshooting.

heading: C.2 LOCAL RR BECOMES CLOSER TO SINGLE-MACHINE RR AS B → N
text: Our next set of plots presented in Figure 2 provides a comparison of minibatch RR, local RR, and single-machine RR (i.e., minibatch RR with M = 1). In Theorems 2 and 4, we showed that when B = Θ(N), then the convergence of local RR becomes just as fast as the single-machine RR.
Indeed, we can observe from Figure 2 that this is really the case. As B increases, the curve of local RR moves closer and closer to that of single-machine RR, especially in the large-epoch regime.

heading: C.3 WITH- VS. WITHOUT-REPLACEMENT SAMPLING
text: Lastly, in Figure 3 we compare the with-replacement and without-replacement versions of minibatch/local SGD. In all plots, we can see that the without-replacement versions outperform withreplacement ones, at least for our problem instance. It is also intriguing to note that the two versions perform very similarly in the small-epoch regime (forK up to∼ 10), but without-replacement starts
to outperform for larger K’s. This observation supports our theoretical prediction from Theorem 3 that in the small-epoch regime, minibatch RR can at best perform as fast as minibatch SGD.

heading: D PROOFS OF UPPER BOUNDS
text: In this section, we provide proofs of our upper bounds stated in Sections 3 and 5. We start by describing a high-level proof outline that we use for all the proofs presented in this section (Appendix D.1). In the subsequent subsections, we prove Theorems 1, 2, 6, and 7, in the order they appeared in the main text. The next subsection (Appendix D.6) states and proves a key lemma that gives concen-
tration bounds for the mean of multiple without-replacement sums of vectors. This general-purpose lemma can be of independent interest and it can prove useful in various other settings. Lastly, in Appendix D.7 we discuss how we can modify the theorem statements to remove the requirement in Assumptions 2, 3, and 4 that they must hold for the entire Rd. Notation. Throughout this section, we use the product notation ∏
in a slightly unconventional manner. For indices i ≤ j and square matrices Ai,Ai+1, . . . ,Aj−1,Aj , we use ∏i l=jAl to denote
the matrix product AjAj−1 · · ·Ai+1Ai. If i > j, then ∏i l=jAl = I .

heading: D.1 PROOF OUTLINE
text: The proofs of upper bounds follow a common structure, consisting of the following three steps:
1. writing one epoch as one step of GD plus noise; 2. getting a high-probability upper bound on the noise term using concentration inequalities; 3. obtaining the convergence rate using the bounds on the noise term.
We first unroll the update equations over an epoch, and write an epoch of the algorithms as one step of gradient descent plus noise:11
xk+1,0 = xk,0 − ηN∇F (xk,0) + η2rk.
Substituting the above to the definition of L-smoothness of F and arranging terms, we obtain
F (xk+1,0)− F (xk,0)
≤ 〈∇F (xk,0),xk+1,0 − xk,0〉+ L
2 ‖xk+1,0 − xk,0‖2
≤ − ηN ‖∇F (xk,0)‖2 + η2 ‖∇F (xk,0)‖ ‖rk‖+ η2L
2 ‖N∇F (xk,0) + ηrk‖2
≤ (−ηN + η2LN2) ‖∇F (xk,0)‖2 + η2 ‖∇F (xk,0)‖ ‖rk‖+ η4L ‖rk‖2 . (37)
The next step is to get high-probability upper bounds on the term ‖rk‖. This is done by applying our concentration inequality lemma (Lemma 8) to partial without-replacement sums of component gradients. As a result, we will get upper bounds on ‖rk‖ and ‖rk‖2, for k = 1, . . . ,K, which hold with probability at least 1− δ. In the last part, we substitute the high-probability bounds to (37) and invoke the definition of PŁ functions to get a per-epoch progress bound. We then unroll the per-epoch inequality for all epochs k = 1, . . . ,K. Arranging the terms in the resulting inequality gives our desired convergence bound that holds with probability at least 1− δ.

heading: D.2 PROOF OF UPPER BOUND FOR MINIBATCH RR (THEOREM 1)
text: One epoch as one step of GD plus noise. To simplify the notation throughout the proof, we will prove the same convergence rate for a rescaled update rule and step-size:
xk,i := xk,i−1 − η
M M∑ m=1 iB∑ j=(i−1)B+1 ∇fmσmk (j)(xk,i−1) (38)
for i ∈ [N/B] and k ∈ [K], and η = log(MNK 2)
µNK . Note that the gradient term is scaled up by B and the step-size is scaled down by B. We will prove the convergence rate for this equivalent algorithm.
We start the proof by unrolling the update equations over an epoch and expressing the progress as
xk+1,0 = xk,0 − ηNF (xk,0) + η2rk,
i.e., one step of full gradient descent plus some noise.
11In case of local RR, we replace xk,0 with yk,0.
To this end, we decompose the gradient∇fmσmk (j)(xk,i−1) into the signal∇f m σmk (j) (xk,0) and noise:
∇fmσmk (j)(xk,i−1) = ∇f m σmk (j) (xk,0) +∇fmσmk (j)(xk,i−1)−∇f m σmk (j) (xk,0) =∇fmσmk (j)(xk,0) + [∫ 1
0
∇2fmσmk (j)(xk,0 + t(xk,i−1 − xk,0))dt ] (xk,i−1 − xk,0),
where ∇2f(x) denotes the Hessian of f at x, whenever it exists. We remark that the integral exists, due to the following reason. Since we assumed that each fmσmk (j) is differentiable and smooth, its gradient ∇fmσmk (j) is Lipschitz continuous, and hence absolutely continuous. This means that ∇fmσmk (j) is differentiable almost everywhere (i.e., ∇ 2fmσmk (j) (x) exists a.e.), and the fundamental theorem of calculus for Lebesgue integral holds; hence the integral exists.
To simplify notation, we define the following for all i ∈ [N/B]:
gi := 1
M M∑ m=1 iB∑ j=(i−1)B+1 ∇fmσmk (j)(xk,0),
Hi := 1
M M∑ m=1 iB∑ j=(i−1)B+1 ∫ 1 0 ∇2fmσmk (j)(xk,0 + t(xk,i−1 − xk,0))dt,
so that we can write (38) as xk,i = xk,i−1 − ηgi − ηHi(xk,i−1 − xk,0). (39)
From L-smoothness of fmi ’s, it is straightforward to check that ‖Hi‖ ≤ LB. Unrolling (39) for i = 1, . . . , N/B, it turns out that we can write
xk+1,0 = xk,0 − η N/B∑ i=1  i+1∏ j=N/B (I − ηHj)  gi. Due to summation by parts, the following identity holds:
N/B∑ i=1 aibi = aN/B N/B∑ j=1 bj − N/B−1∑ i=1 (ai+1 − ai) i∑ j=1 bj .
We apply this to the last term, by substituting ai = ∏i+1 j=N/B(I − ηHj) and bi = gi:
η N/B∑ i=1  i+1∏ j=N/B (I − ηHj)  gi = η
N/B∑ j=1 gj − η N/B−1∑ i=1  i+2∏ t=N/B (I − ηHt)− i+1∏ t=N/B (I − ηHt)  i∑ j=1 gj
= ηN∇F (xk,0)− η2 N/B−1∑ i=1  i+2∏ t=N/B (I − ηHt) Hi+1 i∑ j=1
gj︸ ︷︷ ︸ =:rk .
With the “noise” rk defined as above, we can write xk+1,0 = xk,0 − ηN∇F (xk,0) + η2rk, as desired. Next, it follows from L-smoothness of F that
F (xk+1,0)− F (xk,0)
≤ 〈∇F (xk,0),xk+1,0 − xk,0〉+ L
2 ‖xk+1,0 − xk,0‖2
≤ − ηN ‖∇F (xk,0)‖2 + η2 ‖∇F (xk,0)‖ ‖rk‖+ η2L
2 ‖N∇F (xk,0) + ηrk‖2
≤ (−ηN + η2LN2) ‖∇F (xk,0)‖2 + η2 ‖∇F (xk,0)‖ ‖rk‖+ η4L ‖rk‖2 , (40)
where the last inequality used ‖a + b‖2 ≤ 2 ‖a‖2 + 2 ‖b‖2.
Bounding noise term using concentration. It is left to bound ‖rk‖. We have
‖rk‖ = ∥∥∥∥∥∥ N/B−1∑ i=1  i+2∏ t=N/B (I − ηHt) Hi+1 i∑ j=1 gj ∥∥∥∥∥∥ ≤ N/B−1∑ i=1 ∥∥∥∥∥∥  i+2∏ t=N/B (I − ηHt) Hi+1 i∑ j=1 gj
∥∥∥∥∥∥ ≤ LB(1 + ηLB)N/B
N/B−1∑ i=1 ∥∥∥∥∥∥ i∑ j=1 gj ∥∥∥∥∥∥ , (41) where the last step used ‖Hi‖ ≤ LB for i ∈ [N/B]. Recall from the theorem statement that K ≥ 6κ log(MNK2) and η = log(MNK
2) µNK . This means that
(1 + ηLB)N/B = ( 1 + κB log(MNK2)
NK
)N/B ≤ ( 1 + B
6N
)N/B ≤ e1/6. (42)
Now, we use Lemma 8 to bound the norm of
i∑ j=1 gj = 1 M M∑ m=1 iB∑ j=1 ∇fmσmk (j)(xk,0). (43)
Note that for any epoch k, the permutations σ1k, . . . , σ M k are independent of the first iterate xk,0 of the epoch, and hence independent of all ∇fmi (xk,0). Therefore, we can apply Lemma 8 to the partial sum (43), with vmi ← ∇fmi (xk,0), n ← iB, and δ ← BδNK . By Lemma 8, with probability at least 1− BδNK , we have∥∥∥∥∥∥ 1iBM M∑ m=1 iB∑ j=1 ∇fmσmk (j)(xk,0)−∇F (xk,0) ∥∥∥∥∥∥ ≤ ν √ 8 log 2NKBδ iBM . (44)
Using this concentration bound, with probability at least 1− BδNK we have∥∥∥∥∥∥ i∑
j=1
gj ∥∥∥∥∥∥ ≤ iBν √ 8 log 2NKBδ iBM + iB ‖∇F (xk,0)‖ = ν √ 8iB log 2NKBδ M + iB ‖∇F (xk,0)‖ . (45)
We can now substitute (42) and (45) to (41) to get
‖rk‖ ≤ e1/6LB N/B−1∑ i=1
ν √
8iB log 2NKBδ M + iB ‖∇F (xk,0)‖  ≤ e 1/6 √ 8LνB3/2
M1/2
∫ N/B 1 √ tdt √ log 2NK Bδ + e1/6LB2 ‖∇F (xk,0)‖ N/B−1∑ i=1 i
≤ 5Lν(N 3/2 −B3/2)
2M1/2
√ log 2NK
Bδ + LN(N −B) ‖∇F (xk,0)‖ , (46)
which holds with probability at least 1 − δK , due to the union bound over i = 1, . . . , N/B − 1. The bound (46) holds for all k ∈ [K] with probability 1 − δ if we apply the union bound over k = 1, . . . ,K. Next, by (a+ b)2 ≤ 2a2 + 2b2, we have
‖rk‖2 ≤ 25L2ν2(N3/2 −B3/2)2
2M log
2NK
Bδ + 2L2N2(N −B)2 ‖∇F (xk,0)‖2 , (47)
which also holds for all k ∈ [K] with probability at least 1− δ.
Getting a high-probability convergence rate. Given our high-probability bounds (46) and (47), we can substitute them to (40) and get
F (xk+1,0)− F (xk,0) ≤ (−ηN + η2LN2) ‖∇F (xk,0)‖2 + η2 ‖∇F (xk,0)‖ ‖rk‖+ η4L ‖rk‖2
≤ ( −ηN + η2LN2 + η2LN(N −B) + 2η4L3N2(N −B)2 ) ‖∇F (xk,0)‖2
+ 5η2Lν(N3/2 −B3/2)
2M1/2
√ log 2NK
Bδ ‖∇F (xk,0)‖
+ 25η4L3ν2(N3/2 −B3/2)2
2M log
2NK
Bδ . (48)
The second term in the RHS of (48) can be bounded using ab ≤ a 2 2 + b2 2 :
5η2Lν(N3/2 −B3/2) 2M1/2
√ log 2NK
Bδ ‖∇F (xk,0)‖
=
( η1/2N1/2
2 ‖∇F (xk,0)‖
)( 5η3/2Lν(N3/2 −B3/2)
M1/2N1/2
√ log 2NK
Bδ
)
≤ ηN 8 ‖∇F (xk,0)‖2 +
25η3L2ν2(N3/2 −B3/2)2
2MN log
2NK
Bδ .
Putting this inequality to (48) and noting N −B ≤ N gives F (xk+1,0)− F (xk,0) ≤ ( −7
8 ηN + 2η2LN2 + 2η4L3N4
) ‖∇F (xk,0)‖2
+ 25η3L2(1 + ηLN)ν2(N3/2 −B3/2)2
2MN log
2NK
Bδ . (49)
Recall from K ≥ 6κ log(MNK2) and η = log(MNK 2)
µNK that ηLN ≤ 1 6 . Since the inequality
− 78z + 2z 2 + 2z4 ≤ − 12z holds on z ∈ [0, 1 6 ], we have
−7 8 ηN + 2η2LN2 + 2η4L3N4 ≤ −1 2 ηN.
Applying this bound to (49) results in
F (xk+1,0)− F (xk,0) ≤ − ηN
2 ‖∇F (xk,0)‖2 +
15η3L2ν2(N3/2 −B3/2)2
MN log
2NK
Bδ .
We now recall that F is µ-PŁ, so ‖∇F (xk,0)‖2 ≥ 2µ(F (xk,0)− F ∗):
F (xk+1,0)− F ∗ ≤ (1− ηµN)(F (xk,0)− F ∗) + 15η3L2ν2(N3/2 −B3/2)2
MN log
2NK
Bδ . (50)
Recall that (50) holds for all k ∈ [K], with probability 1− δ. Therefore, by unrolling the inequality,
F (xK,NB )− F ∗ ≤ (1− ηµN)K(F (x0)− F ∗)
+ 15η3L2ν2(N3/2 −B3/2)2
MN log
2NK
Bδ K−1∑ k=0 (1− ηµN)k
≤ (1− ηµN)K(F (x0)− F ∗) + 15η2L2ν2(N3/2 −B3/2)2
µMN2 log
2NK
Bδ . (51)
Lastly, substituting η = log(MNK 2)
µNK gives
F (xK,NB )− F ∗ ≤ F (x0)− F
∗
MNK2 +
15L2ν2(N3/2 −B3/2)2 log 2NKBδ log 2(MNK2)
µ3MN4K2
= F (x0)− F ∗
MNK2 + Õ
( L2ν2
µ3 1 MNK2
) . (52)
Getting an in-expectation bound from the high-probability bound. We conclude this subsection by briefly describing how we can obtain an in-expectation bound from the high-probability bound we just proved. Recall that the bound we proved above holds under the event E that all the concentration bounds used throughout the proof hold. The key to proving an in-expectation bound is to obtain an upper bound under its complement Ec, i.e., conditioned on the event that at least one of our concentration bounds does not hold. We do so by repeating the same proof without ever using the Hoeffding-Serfling bounds (Lemma 8). Of course, this leads to a much looser bound, but we can choose δ to be small enough so that the desired bound Õ ( L2ν2
µ3 1 MNK2
) holds in expectation.
For the version without concentration bounds, the proof proceeds in the same way until it starts diverge at (44). Instead of applying concentration inequalities, we loosely bound the quantity as the following: ∥∥∥∥∥∥ 1iBM M∑ m=1 iB∑ j=1 ∇fmσmk (j)(xk,0)−∇F (xk,0)
∥∥∥∥∥∥ =
∥∥∥∥∥∥ 1iBM M∑ m=1 iB∑ j=1 ∇fmσmk (j)(xk,0)− 1 M M∑ m=1 Fm(xk,0) ∥∥∥∥∥∥ ≤ 1 M M∑ m=1 ∥∥∥∥∥∥ 1iB iB∑ j=1 ∇fmσmk (j)(xk,0)− F m(xk,0)
∥∥∥∥∥∥ ≤ ν. With this bound, the RHS of the upper bound (45) on ‖ ∑i j=1 gj‖ becomes iBν + iB ‖∇F (xk,0)‖. This results in the bounds on ‖rk‖ and ‖rk‖2 (corresponding to (46) and (47)) that read
‖rk‖ ≤ LνN(N −B) + LN(N −B) ‖∇F (xk,0)‖ , ‖rk‖ ≤ 2L2ν2N2(N −B)2 + 2L2N2(N −B)2 ‖∇F (xk,0)‖2 .
The rest is substituting the bounds above to (40), and going through the same steps to obtain the final bound. The resulting bound that corresponds to (51) is
F (xK,NB )− F ∗ ≤ (1− ηµN)K(F (x0)− F ∗) +
7η2L2ν2(N −B)2
3µ ,
which, by substituting η = log(MNK 2)
µNK , yields
F (xK,NB )− F ∗ ≤ F (x0)− F
∗
MNK2 + Õ
( L2ν2
µ3 1 K2
) . (53)
To finish the proof of in-expectation bound, choose δ = 1MN . Recall that the probabilistic event E occurs when all our concentration bounds hold. Conditioned on E, which occurs with probability at least 1 − 1MN , the tighter bound (52) holds, with log 1 δ replaced by log(MN). The complement event Ec occurs with probability at most 1MN , under which the looser bound (53) is true. Thus, in expectation,
E [ F (xK,NB )− F ∗ ] = P(E)E [ F (xK,NB )− F ∗ | E ] + P(Ec)E [ F (xK,NB )− F ∗ | Ec ]
≤ E [ F (xK,NB )− F ∗ | E ] + 1 MN E [ F (xK,NB )− F ∗ | Ec ] ≤ 3(F (x0)− F ∗)
2MNK2 + Õ
( L2ν2
µ3 1 MNK2
) .
For the remaining high-probability upper bounds proved in the paper, we can similarly follow this process to obtain matching (up to log factors) in-expectation upper bounds.

heading: D.3 PROOF OF UPPER BOUND FOR LOCAL RR (THEOREM 2)
text: One epoch as one step of GD plus noise. The update rule of local RR can be written as the following. For k ∈ [K], i ∈ [N ], and m ∈ [M ],
xmk,i :=
{ xmk,i−1 − η∇fmσmk (i)(x m k,i−1) if B does not divide i,
1 M ∑M m=1(x m k,i−1 − η∇fmσmk (i)(x m k,i−1)) =: yk,i/B if B divides i.
(54)
Recall that xmk,0’s are the initial points of an epoch and they all the same regardless of the machine m. We define yk,0 := x1k,0. For any i in the range of (l − 1)B + 1 ≤ i ≤ lB for some l ∈ [M ], we will use yk,l−1 as the “pivot” and decompose the gradients into the ones evaluated at yk,l−1 plus noise terms.
∇fmσmk (i)(x m k,i−1) = ∇fmσmk (i)(yk,l−1) +∇f m σmk (i) (xmk,i−1)−∇fmσmk (i)(yk,l−1) =∇fmσmk (i)(yk,l−1) + [∫ 1
0
∇2fmσmk (i)(yk,l−1 + t(x m k,i−1 − yk,l−1))dt ] ︸ ︷︷ ︸
=:Hmi
(xmk,i−1 − yk,l−1),
where the integral Hmi exists due to the reason discussed in Appendix D.2. Also note from Lsmoothness of fmi ’s that ‖Hmi ‖ ≤ L. Using the decomposition, one can unroll the updates (54) and write yk,l in terms of yk,l−1 in the following way:
yk,l = yk,l−1 − η
M M∑ m=1 lB∑ i=(l−1)B+1  i+1∏ j=lB (I − ηHmj ) ∇fmσmk (i)(yk,l−1), (55) for l = 1, . . . , N/B. Next, we again decompose the gradient ∇fmσmk (i)(yk,l−1), this time using yk,0 as the pivot:
∇fmσmk (i)(yk,l−1) = ∇f m σmk (i) (yk,0) +∇fmσmk (i)(yk,l−1)−∇f m σmk (i) (yk,0) =∇fmσmk (i)(yk,0) + [∫ 1
0
∇2fmσmk (i)(yk,0 + t(yk,l−1 − yk,0))dt ]
︸ ︷︷ ︸ =:H̃mi
(yk,l−1 − yk,0).
This decomposition allows us to rewrite (55) in the following form
yk,l = yk,l−1 − ηtl − ηSl(yk,l−1 − yk,0), (56)
where
tl := 1
M M∑ m=1 lB∑ i=(l−1)B+1  i+1∏ j=lB (I − ηHmj ) ∇fmσmk (i)(yk,0), (57) Sl := 1
M M∑ m=1 lB∑ i=(l−1)B+1  i+1∏ j=lB (I − ηHmj )  H̃mi . (58) Unrolling (56) for l = 1, . . . , N/B then gives the progress over an epoch:
yk+1,0 = yk,0 − η N/B∑ l=1  l+1∏ j=N/B (I − ηSj)  tl. (59) As done in the proof of Theorem 1 (Appendix D.2), we will express (59) as one step of GD on F plus some noise. Of course, the noise terms here will be more complicated to handle than they were in Theorem 1. Due to summation by parts, the following identity holds:
N/B∑ l=1 albl = aN/B N/B∑ j=1 bj − N/B−1∑ l=1 (al+1 − al) l∑ j=1 bj .
We apply this to the last term of (59), by substituting ai = ∏l+1 j=N/B(I − ηSj) and bl = tl:
η N/B∑ l=1  l+1∏ j=N/B (I − ηSj)  tl = η N/B∑ l=1 tl − η2 N/B−1∑ l=1  l+2∏ j=N/B (I − ηSj) Sl+1 l∑ j=1 tj . (60)
We also apply the summation by parts to the inner summation of tl’s (57):
tl := 1
M M∑ m=1 lB∑ i=(l−1)B+1  i+1∏ j=lB (I − ηHmj ) ∇fmσmk (i)(yk,0) = 1
M M∑ m=1 lB∑ i=(l−1)B+1 ∇fmσmk (i)(yk,0)
− η M M∑ m=1 lB−1∑ i=(l−1)B+1 ( i+2∏ t=lB (I − ηHmt ) ) Hmi+1 i∑ j=(l−1)B+1 ∇fmσmk (j)(yk,0) (61)
Substituting (61) to (60) gives
yk+1,0 = yk,0 − ηN∇F (yk,0) + η2rk,1 + η2rk,2 − η3rk,3, where rk,1, rk,2, and rk,3 are noise terms defined as
rk,1 := 1
M N/B∑ l=1 M∑ m=1 lB−1∑ i=(l−1)B+1 ( i+2∏ t=lB (I − ηHmt ) ) Hmi+1 i∑ t=(l−1)B+1 ∇fmσmk (t)(yk,0),
rk,2 := 1
M N/B−1∑ l=1  l+2∏ j=N/B (I − ηSj) Sl+1 M∑ m=1 lB∑ t=1 ∇fmσmk (t)(yk,0),
rk,3 := 1
M N/B−1∑ l=1  l+2∏ j=N/B (I − ηSj) Sl+1× l∑
j=1 M∑ m=1 jB−1∑ i=(j−1)B+1  i+2∏ t=jB (I − ηHmt ) Hmi+1 i∑ t=(j−1)B+1 ∇fmσmk (t)(yk,0).
Defining rk := rk,1 + rk,2 − ηrk,3, it follows from L-smoothness of F that F (yk+1,0)− F (yk,0)
≤ 〈∇F (yk,0),yk+1,0 − yk,0〉+ L
2 ‖yk+1,0 − yk,0‖2
≤ − ηN ‖∇F (yk,0)‖2 + η2 ‖∇F (yk,0)‖ ‖rk‖+ η2L
2 ‖N∇F (yk,0) + ηrk‖2
≤ (−ηN + η2LN2) ‖∇F (yk,0)‖2 + η2 ‖∇F (yk,0)‖ ‖rk‖+ η4L ‖rk‖2 . (62)
Bounding noise terms using concentration. We next bound ‖rk‖ by bounding each ‖rk,1‖, ‖rk,2‖, and ‖rk,3‖. From this point on, we write gmt := ∇fmσmk (t)(yk,0) to simplify notation.
‖rk,1‖ = ∥∥∥∥∥∥ 1M N/B∑ l=1 M∑ m=1 lB−1∑ i=(l−1)B+1 ( i+2∏ t=lB (I − ηHmt ) ) Hmi+1 i∑ t=(l−1)B+1 gmt ∥∥∥∥∥∥ ≤ 1 M N/B∑ l=1 M∑ m=1 lB−1∑ i=(l−1)B+1 ∥∥∥∥∥∥ ( i+2∏ t=lB (I − ηHmt ) ) Hmi+1 i∑ t=(l−1)B+1 gmt
∥∥∥∥∥∥ ≤ L(1 + ηL) B
M
N/B∑ l=1 M∑ m=1 lB−1∑ i=(l−1)B+1 ∥∥∥∥∥∥ i∑ t=(l−1)B+1 gmt ∥∥∥∥∥∥ , (63)
where we used ‖Hmi ‖ ≤ L. Recall from the theorem statement that K ≥ 7ρκ log(MNK2) and η = log(MNK
2) µNK . This means that
(1 + ηL)B = ( 1 + κ log(MNK2)
NK
)B ≤ ( 1 + 1
7ρN
)B ≤ ( 1 + 1
7B
)B ≤ e1/7. (64)
Also note from the definition of Sl (58) that ‖Sl‖ ≤ LB(1 + ηL)B ≤ e1/7LB, which we use to get similar bounds for the next two terms rk,2 and rk,3.
‖rk,2‖ = ∥∥∥∥∥∥ 1M N/B−1∑ l=1  l+2∏ j=N/B (I − ηSj) Sl+1 M∑ m=1 lB∑ t=1 gmt ∥∥∥∥∥∥ ≤ e 1/7LB(1 + e1/7ηLB)N/B
M
N/B−1∑ l=1 ∥∥∥∥∥ M∑ m=1 lB∑ t=1 gmt ∥∥∥∥∥ , (65) and we can bound
(1 + e1/7ηLB)N/B = ( 1 + e1/7κB log(MNK2)
NK )N/B ≤ ( 1 + e1/7B
7ρN
)N/B ≤ exp ( e1/7
7
) . (66)
We similarly bound the norm of the last noise term rk,3:
‖rk,3‖ ≤ e1/7LB(1 + e1/7ηLB)N/B
M ×
N/B−1∑ l=1 ∥∥∥∥∥∥ l∑ j=1 M∑ m=1 jB−1∑ i=(j−1)B+1  i+2∏ t=jB (I − ηHmt ) Hmi+1 i∑ t=(j−1)B+1 gmt ∥∥∥∥∥∥ ≤ e 1/7L2B(1 + e1/7ηLB)N/B(1 + ηL)B
M
N/B−1∑ l=1 l∑ j=1 M∑ m=1 jB−1∑ i=(j−1)B+1 ∥∥∥∥∥∥ i∑ t=(j−1)B+1 gmt ∥∥∥∥∥∥ ≤ 11L 2B
7M N/B−1∑ l=1 l∑ j=1 M∑ m=1 jB−1∑ i=(j−1)B+1 ∥∥∥∥∥∥ i∑ t=(j−1)B+1 gmt ∥∥∥∥∥∥ , (67) where the last inequality used (64), (66), and e2/7 exp(e1/7/7) ≤ 11/7. Given the bounds (63), (65), and (67), we now use Lemma 8 to get high-probability bounds for the partial sums of gmt that appear in the bounds. For any i satisfying (j−1)B+1 ≤ i ≤ jB−1, where j ∈ [N/B], and for any m ∈ [M ], the following bound holds with probability at least 1− δ2MNK :∥∥∥∥∥∥ 1i− (j − 1)B i∑ t=(l−1)B+1 gmt − 1 N N∑ t=1 gmt
∥∥∥∥∥∥ = ∥∥∥∥∥∥ 1i− (j − 1)B i∑
t=(l−1)B+1
∇fmσmk (t)(yk,0)−∇F m(yk,0) ∥∥∥∥∥∥ ≤ ν √ 8 log 4MNKδ i− (j − 1)B .
From this, with probability at least 1− δ2MNK we have∥∥∥∥∥∥ i∑
t=(j−1)B+1
gmt ∥∥∥∥∥∥ ≤ ν √ 8(i− (j − 1)B) log 4MNKδ + (i− (j − 1)B) ‖∇F m(yk,0)‖ . (68)
Similarly, for l ∈ [N/B − 1], the following bound holds with probability at least 1− Bδ2NK :∥∥∥∥∥ 1lBM M∑ m=1 lB∑ t=1 gmt − 1 MN M∑ m=1 N∑ t=1 gmt ∥∥∥∥∥
= ∥∥∥∥∥ 1lBM M∑ m=1 lB∑ t=1 ∇fmσmk (t)(yk,0)−∇F (yk,0) ∥∥∥∥∥ ≤ ν √ 8 log 4NKBδ lBM ,
which gives us ∥∥∥∥∥ M∑ m=1 lB∑ t=1 gmt ∥∥∥∥∥ ≤ ν √ 8lBM log 4NKBδ + lBM ‖∇F (yk,0)‖ (69)
By applying the union bound, with probability at least 1− δK , the bound (68) holds for all m ∈ [M ] and i ∈ ⋃N/B j=1 [(j − 1)B + 1 : jB − 1], and the bound (69) holds for all l ∈ [N/B − 1].
We now substitute the bounds (68) and (69) to (63), (65), and (67) to get upper bounds for ‖rk,1‖, ‖rk,2‖, and ‖rk,3‖, respectively. First,
‖rk,1‖
≤ e 1/7L
M N/B∑ l=1 M∑ m=1 lB−1∑ i=(l−1)B+1 ( ν √ 8(i− (l − 1)B) log 4MNKδ + (i− (l − 1)B) ‖∇F m(yk,0)‖ )
= e1/7 √ 8LνN
B
( B−1∑ i=1 √ i )√ log 4MNK δ + e1/7LN B ( B−1∑ i=1 i )( 1 M M∑ m=1 ‖∇Fm(yk,0)‖ )
≤ 24LνN(B 3/2 − 1)
11B
√ log 4MNK
δ + 3LN(B − 1) 5 (τ + ρ ‖∇F (yk,0)‖) , (70)
where the last inequality used ∑B−1 i=1 √ i ≤ ∫ B 1 √ zdz = 23 (B
3/2 − 1) and Assumption 3. For the next noise term, we have e1/7(1 + e1/7ηLB)N/B ≤ 7/5, so
‖rk,2‖ ≤ 7LB
5M N/B−1∑ l=1 ( ν √ 8lBM log 4NKBδ + lBM ‖∇F (yk,0)‖ )
= 7 √ 8LνB3/2
5M1/2 N/B−1∑ l=1 √ l √log 4NK Bδ + 7LB2 5 N/B−1∑ l=1 l  ‖∇F (yk,0)‖ ≤ 8Lν(N
3/2 −B3/2) 3M1/2
√ log 4NK
Bδ + 7LN(N −B) 10 ‖∇F (yk,0)‖ . (71)
Lastly,
‖rk,3‖ ≤ 11L2B
7M N/B−1∑ l=1 l∑ j=1 M∑ m=1 jB−1∑ i=(j−1)B+1 ( ν √ 8(i− (j − 1)B) log 4MNKδ
+ (i− (j − 1)B) ‖∇Fm(yk,0)‖ )
= 11 √ 8L2νB
7
N/B−1∑ l=1 l (B−1∑ i=1 √ i )√ log 4MNK δ
+ 11L2B
7 N/B−1∑ l=1 l (B−1∑ i=1 i )( 1 M M∑ m=1 ‖∇Fm(yk,0)‖ )
≤ 3L 2νN(N −B)(B3/2 − 1)
2B
√ log 4MNK
δ
+ 2L2N(N −B)(B − 1)
5 (τ + ρ ‖∇F (yk,0)‖). (72)
Recalling the definition rk := rk,1 + rk,2− ηrk,3, we get an upper bound for ‖rk‖ from (70), (71), and (72):
‖rk‖ ≤ ‖rk,1‖+ ‖rk,2‖+ η ‖rk,3‖ ≤ Lν √ log 4MNK
δ
( 24N(B3/2 − 1)
11B + 8(N3/2 −B3/2) 3M1/2 + 3ηLN(N −B)(B3/2 − 1) 2B ) + Lτ ( 3N(B − 1)
5 + 2ηLN(N −B)(B − 1) 5 ) + L ‖∇F (yk,0)‖ ( 3ρN(B − 1)
5 + 7N(N −B) 10 + 2ηLρN(N −B)(B − 1) 5
) . (73)
Recall again that we have K ≥ 7ρκ log(MNK2) and η = log(MNK 2)
µNK , so ηLN ≤ 1/7. Using this and N −B ≤ N , we can further simplify (73).
‖rk‖ ≤ Lν √ log 4MNK
δ
( 12N(B3/2 − 1)
5B + 8(N3/2 −B3/2) 3M1/2 ) ︸ ︷︷ ︸
=:Φ
+ 2LτN(B − 1)
3
+ L ‖∇F (yk,0)‖ (
2ρN(B − 1) 3 + 7N(N −B) 10 ) ≤ LνΦ √ log 4MNK
δ + 2LτN(B − 1) 3 + 7LρN2 5 ‖∇F (yk,0)‖ , (74)
which holds with probability at least 1− δK . The bound (74) holds for all k ∈ [K] with probability 1 − δ if we apply the union bound over k = 1, . . . ,K. Next, by (a + b + c)2 ≤ 3a2 + 3b2 + 3c2, we have
‖rk‖2 ≤ 3L2ν2Φ2 log 4MNK
δ +
4L2τ2N2(B − 1)2
3 +
147L2ρ2N4
25 ‖∇F (yk,0)‖2 , (75)
which also holds for all k ∈ [K] with probability at least 1− δ.
Getting a high-probability convergence rate. Given our high-probability bounds (74) and (75), we can substitute them to (62) and get
F (yk+1,0)− F (yk,0) ≤ (−ηN + η2LN2) ‖∇F (yk,0)‖2 + η2 ‖∇F (yk,0)‖ ‖rk‖+ η4L ‖rk‖2
≤ ( −ηN + η2LN2 + 7η 2LρN2
5 +
147η4L3ρ2N4
25
) ‖∇F (yk,0)‖2
+ η2LνΦ √ log 4MNK
δ ‖∇F (yk,0)‖+ 3η4L3ν2Φ2 log
4MNK
δ
+ 2η2LτN(B − 1)
3 ‖∇F (yk,0)‖+
4η4L3τ2N2(B − 1)2
3 . (76)
The following terms in the RHS of (76) can be bounded using ab ≤ a 2 2 + b2 2 :
η2LνΦ √ log 4MNK
δ ‖∇F (yk,0)‖
= ( η1/2N1/2√
8 ‖∇F (yk,0)‖
)(√ 8η3/2LνΦ
N1/2
√ log 4MNK
δ
)
≤ ηN 16 ‖∇F (yk,0)‖2 +
4η3L2ν2Φ2
N log
4MNK
δ , (77)
2η2LτN(B − 1) 3 ‖∇F (yk,0)‖
= ( η1/2N1/2√
8 ‖∇F (yk,0)‖
)( 2 √
8η3/2LτN1/2(B − 1) 3
)
≤ ηN 16 ‖∇F (yk,0)‖2 +
16η3L2τ2N(B − 1)2
9 . (78)
Substituting (77) and (78) to (76) results in
F (yk+1,0)− F (yk,0) ≤ ( −7
8 ηN + η2LN2 +
7η2LρN2
5 +
147η4L3ρ2N4
25
) ‖∇F (yk,0)‖2
+ η3L2ν2(4 + 3ηLN)Φ2
N log
4MNK
δ +
4η3L2τ2(4 + 3ηLN)N(B − 1)2
9 . (79)
Again, we have K ≥ 7ρκ log(MNK2) and η = log(MNK 2)
µNK , so ηLρN ≤ 1 7 . Since the inequality
− 78z + 12 5 z 2 + 14725 z 4 ≤ − 12z holds on z ∈ [0, 1 7 ], we have
− 7 8 ηN + η2LN2 +
7η2LρN2
5 +
147η4L3ρ2N4
25
≤ − 7 8 ηN +
12η2LρN2
5 +
147η4L3ρ3N4
25 ≤ −1 2 ηN.
Substituting this inequality to (79), together with 4 + 3ηLN ≤ 317 < 9 2 , yields
F (yk+1,0)− F (yk,0) ≤ − ηN
2 ‖∇F (yk,0)‖2 +
9η3L2ν2Φ2
2N log
4MNK
δ + 2η3L2τ2N(B − 1)2.
We now recall that F is µ-PŁ, so ‖∇F (yk,0)‖2 ≥ 2µ(F (yk,0)− F ∗): F (yk+1,0)− F ∗ ≤ (1− ηµN)(F (yk,0)− F ∗)
+ 9η3L2ν2Φ2
2N log
4MNK
δ + 2η3L2τ2N(B − 1)2. (80)
Recall that (80) holds for all k ∈ [K], with probability 1− δ. Therefore, by unrolling the inequality,
F (yK,NB )− F ∗ ≤ (1− ηµN)K(F (y0)− F ∗)
+
( 9η3L2ν2Φ2
2N log
4MNK
δ + 2η3L2τ2N(B − 1)2 )K−1∑ k=0 (1− ηµN)k
≤ (1− ηµN)K(F (y0)− F ∗)
+ 9η2L2ν2Φ2
2µN2 log
4MNK
δ +
2η2L2τ2(B − 1)2
µ . (81)
Recall that Φ := 12N(B 3/2−1) 5B + 8(N3/2−B3/2) 3M1/2 , hence
Φ2 ≤ 288N 2(B3/2 − 1)2
25B2 +
128(N3/2 −B3/2)2
9M .
Substituting this inequality and also η = log(MNK 2)
µNK gives
F (yK,NB )− F ∗
≤ F (y0)− F ∗
MNK2 +
2L2τ2(B − 1)2
µ3N2K2 log2(MNK2)
+ 9L2ν2
2µ3N4K2 log
4MNK
δ log2(MNK2)
( 288N2(B3/2 − 1)2
25B2 +
128(N3/2 −B3/2)2
9M ) = F (y0)− F ∗
MNK2 + Õ
( L2τ2
µ3 B2 N2K2
) + Õ ( L2ν2
µ3 B N2K2
) + Õ ( L2ν2
µ3 1 MNK2
) ,
with probability at least 1− δ. This finishes the proof.

heading: D.4 PROOF OF UPPER BOUND FOR MINIBATCH RR WITH SYNCSHUF (THEOREM 6)
text: The first part (“One epoch as one step of GD plus noise”) of the proof is identical to that of Theorem 1. We start from the second part.
Bounding noise term using concentration. It is left to bound ‖rk‖. As seen in (41), we have
‖rk‖ ≤ LB(1 + ηLB)N/B N/B−1∑ i=1 ∥∥∥∥∥∥ i∑ j=1 gj ∥∥∥∥∥∥ . (82) Recall from the theorem statement thatK ≥ 6κ log(MNK2) and η = log(MNK
2) µNK . This means that
(1 + ηLB)N/B = ( 1 + κB log(MNK2)
NK
)N/B ≤ ( 1 + B
6N
)N/B ≤ e1/6. (83)
Next, we bound the norm of i∑
j=1
gj = 1
M M∑ m=1 iB∑ j=1 ∇fmσmk (j)(xk,0), (84)
exploiting our modification SYNCSHUF as well as Lemma 8. For each ∇fmσmk (j)(xk,0), we first add and subtract its corresponding ∇f̄σmk (j)(xk,0), where f̄i := 1 M ∑M m=1 f m i as defined in Assump-
tion 4. This way, (84) can be decomposed into two sums ∑i j=1 gj = 1 M (pi + qi), where
pi := M∑ m=1 iB∑ j=1 ∇f̄σmk (j)(xk,0),
qi := M∑ m=1 iB∑ j=1 ∇fmσmk (j)(xk,0)−∇f̄σmk (j)(xk,0).
Using this decomposition, we will derive high-probability bounds for ‖pi‖ and ‖qi‖. To simplify expressions to follow, we decompose iBM (i.e., the total number of component gradients that are summed up) into a multiple of N and the remainder. Let
α(i) :=
⌊ iBM
N
⌋ , β(i) := iBM −Nα(i),
so that iBM is decomposed into Nα(i) and the remainder 0 ≤ β(i) < N . Using this new notation, we can write pi as
pi = M∑ m=1 Nα(i) M∑ j=1 ∇f̄σmk (j)(xk,0) + M∑ m=1 iB∑ j=
Nα(i) M +1
∇f̄σmk (j)(xk,0). (85)
Here, recall that with SYNCSHUF, we defined σmk (j) := σ((j+ N M π(m)) mod N). With this choice of “shifted” permutations, one can notice that {σmk (j)} M, NM m=1,j=1 = [N ], meaning that adding f̄σmk (j) for m ∈ [M ] and j ∈ [N/M ] results in the sum of all N f̄i’s. In fact, this happens if we sum over m ∈ [M ] and any N/M consecutive j’s. From this observation and F = 1N ∑N i=1 f̄i, (85) can be written as
pi = Nα(i)∇F (xk,0) + M∑ m=1 iB∑ j=
Nα(i) M +1
∇f̄σ((j+NmM ) mod N)(xk,0). (86)
Assume for now that β(i) > 0, i.e., Nα(i) < iBM . The summation in the second term of RHS in (86) is a without-replacement sum (note that the indices j + NmM do not overlap) of β(i) terms.
Hence, it is equal in distribution to ∑β(i) j=1∇f̄σ(j)(xk,0). Also, from Assumption 2, it can be easily
checked that for any i ∈ [N ] ∥∥∇f̄i(xk,0)−∇F (xk,0)∥∥ ≤ ν. These observations mean that we can apply Lemma 8 to get a concentration bound, with M ← 1, v1i ← ∇f̄i(xk,0), n← β(i), and δ ← Bδ2NK . By Lemma 8, with probability at least 1− Bδ 2NK (over
the randomness in σ), we have∥∥∥∥∥∥∥ 1 β(i) M∑ m=1 iB∑ j=
Nα(i) M +1
∇f̄σ((j+NmM ) mod N)(xk,0)−∇F (xk,0) ∥∥∥∥∥∥∥ ≤ ν √ 8 log 4NKBδ β(i) . (87)
Combining (86) and (87), we get the following upper bound on ‖pi‖, which holds with probability at least 1− Bδ2NK .
‖pi‖ ≤ ‖iBM∇F (xk,0)‖+ ∥∥∥∥∥∥∥ M∑ m=1 iB∑ j=
Nα(i) M +1
∇f̄σ((j+NmM ) mod N)(xk,0)− β(i)∇F (xk,0) ∥∥∥∥∥∥∥ ≤ iBM ‖∇F (xk,0)‖+ ν √ 8β(i) log 4NK
Bδ
≤ iBM ‖∇F (xk,0)‖+ ν √ N √ 8 log 4NK
Bδ , (88)
where the last inequality used β(i) < N . Also recall that we assumed β(i) > 0 in order to use Lemma 8 and derive (88). However, note that even with β(i) = 0, the bound (88) trivially holds.
We next bound ‖qi‖. This time, we will apply Lemma 8 to the permutation π over the local machines. To do this, we will condition on a fixed instantiation of the permutation σ and derive a high-probability bound that holds with conditional probability at least 1 − Bδ2NK . The conditional probability is at least 1− Bδ2NK irrespective of the choice of σ, so we can conclude that the (unconditional) probability that our bound holds is also at least 1− Bδ2NK . Without loss of generality, choose the instantiation σ(l) = l for all l ∈ [N ]. With this σ, we have σmk (j) := (j + N M π(m)) mod N , so the vector qi reads
qi = M∑ m=1 iB∑ j=1 ∇fm (j+ NM π(m)) mod N (xk,0)−∇f̄(j+ NM π(m)) mod N (xk,0). (89)
Let us consider rewriting this summation as the sum over l ∈ [N ], where l appears in the subscript of the component functions. One can check that
l = ( j + N
M π(m)
) mod N ⇔ N
M | (l − j) and π(m) = (l − j)M N mod M,
where a | b denotes “a divides b.” From this, we can rewrite (89) as
qi = N∑ l=1 ∑ j∈[iB] N M |(l−j) ( ∇fπ −1((l−j)MN mod M) l (xk,0)−∇f̄l(xk,0) ) ︸ ︷︷ ︸
=:qi,l
. (90)
By the same reasoning above and below (86), we can see from (90) that for a given index l ∈ [N ], the cardinality of the set Jl := {j ∈ [iB] : NM | (l − j)} is either α(i) or α(i) + 1. From this, we notice that each qi,l is a without-replacement sum of α(i) or α(i) + 1 terms. For now, suppose α(i) > 0. For each qi,l, we can apply Lemma 8 to it, and show that with probability (conditioned on the instantiation σ(l) = l) at least 1− Bδ2N2K , we have
‖qi,l‖ ≤ λ √ 8(α(i)) log 4N 2K Bδ if |Jl| = α(i), λ √
8(α(i) + 1) log 4N 2K Bδ if |Jl| = α(i) + 1.
Note that cases in the RHS are all bounded from above by λ √
16α(i) log 4N 2K Bδ . Applying union
bound on all l ∈ [N ], we get that with probability at least 1− Bδ2NK , we have ‖qi‖ ≤ λN √ 16α(i) log 4N2K
Bδ ≤ λN
√ 16iBM
N log
4N2K
Bδ = 4λ
√ iBMN log 4N2K
Bδ . (91) Now consider the case α(i) = 0. Recall from the definition α(i) := ⌊ iBM N ⌋ that α(i) = 0 implies iBM < N . In this case, qi,l = 0 for N − iBM indices l satisfying |Jl| = 0, and ‖qi,l‖ ≤ λ for the remaining iBM l’s satisfying |Jl| = 1. Summing up, ‖qi‖ is bounded from above by λiBM , which is in fact less than the upper bound in (91). Therefore, the bound (91) holds even for α(i) = 0.
Recall that our goal was to find a bound on the norm of ∑i j=1 gj = 1 M (pi + qi). From the highprobability bounds obtained in (88) and (91), with probability at least 1− BδNK ,∥∥∥∥∥∥ i∑
j=1
gj ∥∥∥∥∥∥ ≤ iB ‖∇F (xk,0)‖+ ν √ N M √ 8 log 4NK Bδ + 4λ √ iBN M √ log 4N2K Bδ . (92)
We can now substitute (83) and (92) to (82) to get
‖rk‖ ≤ e1/6LB N/B−1∑ i=1
( iB ‖∇F (xk,0)‖+ ν √ N
M
√ 8 log 4NK
Bδ + 4λ
√ iBN
M
√ log 4N2K
Bδ
)
≤ e1/6LB2 ‖∇F (xk,0)‖ N/B−1∑ i=1 i+ e1/6 √ 8LνN1/2(N −B) M √ log 4NK Bδ
+ 4e1/6LλB3/2N1/2
M1/2
∫ N/B 1 √ tdt √ log 4N2K Bδ
≤ LN(N −B) ‖∇F (xk,0)‖+ 7LνN1/2(N −B)
2M
√ log 4NK
Bδ
+ 7LλN1/2(N3/2 −B3/2)
2M1/2
√ log 4N2K
Bδ , (93)
which holds with probability at least 1 − δK , due to the union bound over i = 1, . . . , N/B − 1. The bound (93) holds for all k ∈ [K] with probability 1 − δ if we apply the union bound over k = 1, . . . ,K. Next, by (a+ b+ c)2 ≤ 3a2 + 3b2 + 3c2, we have
‖rk‖2 ≤ 3L2N2(N −B)2 ‖∇F (xk,0)‖2 + 147L2ν2N(N −B)2
4M2 log
4NK
Bδ
+ 147L2λ2N(N3/2 −B3/2)2
4M log
4N2K
Bδ , (94)
which also holds for all k ∈ [K] with probability at least 1− δ.
Getting a high-probability convergence rate. Given our high-probability bounds (93) and (94), we can substitute them to (40) and get
F (xk+1,0)− F (xk,0) ≤ (−ηN + η2LN2) ‖∇F (xk,0)‖2 + η2 ‖∇F (xk,0)‖ ‖rk‖+ η4L ‖rk‖2
≤ ( −ηN + η2LN2 + η2LN(N −B) + 3η4L3N2(N −B)2 ) ‖∇F (xk,0)‖2
+ 7η2LνN1/2(N −B)
2M
√ log 4NK
Bδ ‖∇F (xk,0)‖+
147η4L3ν2N(N −B)2
4M2 log
4NK
Bδ
+ 7η2LλN1/2(N3/2 −B3/2)
2M1/2
√ log 4N2K
Bδ ‖∇F (xk,0)‖
+ 147η4L3λ2N(N3/2 −B3/2)2
4M log
4N2K
Bδ . (95)
Two terms in the RHS of (95) can be bounded using ab ≤ a 2 2 + b2 2 :
7η2LνN1/2(N −B) 2M
√ log 4NK
Bδ ‖∇F (xk,0)‖
=
( η1/2N1/2
2 √ 2 ‖∇F (xk,0)‖
)( 7 √
2η3/2Lν(N −B) M
√ log 4NK
Bδ
)
≤ ηN 16 ‖∇F (xk,0)‖2 +
49η3L2ν2(N −B)2
M2 log
4NK
Bδ , (96)
7η2LλN1/2(N3/2 −B3/2) 2M1/2
√ log 4N2K
Bδ ‖∇F (xk,0)‖
=
( η1/2N1/2
2 √ 2 ‖∇F (xk,0)‖
)( 7 √
2η3/2Lλ(N3/2 −B3/2) M1/2
√ log 4N2K
Bδ
)
≤ ηN 16 ‖∇F (xk,0)‖2 +
49η3L2λ2(N3/2 −B3/2)2
M log
4N2K
Bδ , (97)
Putting inequalities (96) and (97) to (95) and noting N −B ≤ N gives F (xk+1,0)− F (xk,0) ≤ ( −7
8 ηN + 2η2LN2 + 3η4L3N4
) ‖∇F (xk,0)‖2
+ 49η3L2(4 + 3ηLN)ν2(N −B)2
4M2 log
4NK
Bδ
+ 49η3L2(4 + 3ηLN)λ2(N3/2 −B3/2)2
4M log
4N2K
Bδ . (98)
Recall from K ≥ 6κ log(M2NK2) and η = log(M 2NK2)
µNK that ηLN ≤ 1 6 . Since the inequality
− 78z + 2z 2 + 3z4 ≤ − 12z holds on z ∈ [0, 1 6 ], we have
−7 8 ηN + 2η2LN2 + 3η4L3N4 ≤ −1 2 ηN.
Applying this bound to (98) results in
F (xk+1,0)− F (xk,0) ≤ − ηN
2 ‖∇F (xk,0)‖2 +
56η3L2ν2(N −B)2
M2 log
4NK
Bδ
+ 56η3L2λ2(N3/2 −B3/2)2
M log
4N2K
Bδ
We now recall that F is µ-PŁ, so ‖∇F (xk,0)‖2 ≥ 2µ(F (xk,0)− F ∗):
F (xk+1,0)− F ∗ ≤ (1− ηµN)(F (xk,0)− F ∗) + 56η3L2ν2(N −B)2
M2 log
4NK
Bδ
+ 56η3L2λ2(N3/2 −B3/2)2
M log
4N2K
Bδ (99)
Recall that (99) holds for all k ∈ [K], with probability 1− δ. Therefore, by unrolling the inequality, and using ∑K−1 k=0 (1− ηµN)k ≤ 1 ηµN , we get
F (xK,NB )− F ∗ ≤ (1− ηµN)K(F (x0)− F ∗) +
56η2L2ν2(N −B)2
µM2N log
4NK
Bδ
+ 56η2L2λ2(N3/2 −B3/2)2
µMN log
4N2K
Bδ (100)
Lastly, substituting η = log(M 2NK2)
µNK to (100) gives
F (xK,NB )− F ∗ ≤ F (x0)− F
∗
M2NK2 +
56L2ν2(N −B)2 log 4NKBδ log 2(M2NK2)
µ3M2N3K2
+ 56L2λ2(N3/2 −B3/2)2 log 4N 2K Bδ log 2(M2NK2)
µ3MN3K2
= F (x0)− F ∗
M2NK2 + Õ
( L2
µ3
( ν2
M2NK2 +
λ2
MK2
)) .

heading: D.5 PROOF OF UPPER BOUND FOR LOCAL RR WITH SYNCSHUF (THEOREM 7)
text: The first part (“One epoch as one step of GD plus noise”) of the proof is identical to that of Theorem 2. The first part defines our “noise” rk as the sum of three terms rk := rk,1 + rk,2− ηrk,3. We start from the second part.
Bounding noise terms using concentration. We next bound ‖rk‖ by bounding each ‖rk,1‖, ‖rk,2‖, and ‖rk,3‖. We have already seen from (63), (64), (65), (66), and (67) in Appendix D.3 that
‖rk,1‖ ≤ e1/7L
M N/B∑ l=1 M∑ m=1 lB−1∑ i=(l−1)B+1 ∥∥∥∥∥∥ i∑ t=(l−1)B+1 ∇fmσmk (t)(yk,0) ∥∥∥∥∥∥ , (101) ‖rk,2‖ ≤ 7LB
5M N/B−1∑ l=1 ∥∥∥∥∥ M∑ m=1 lB∑ t=1 ∇fmσmk (t)(yk,0) ∥∥∥∥∥ , (102) ‖rk,3‖ ≤ 11L2B
7M N/B−1∑ l=1 l∑ j=1 M∑ m=1 jB−1∑ i=(j−1)B+1 ∥∥∥∥∥∥ i∑ t=(j−1)B+1 ∇fmσmk (t)(yk,0) ∥∥∥∥∥∥ . (103) As in the previous subsections, the key is to bound the norm of the partial sums of ∇fmσmk (t)(yk,0) using Lemma 8. For the summations appearing in (101) and (103), we apply Lemma 8 in the same way as (68). For any i satisfying (j − 1)B + 1 ≤ i ≤ jB − 1, where j ∈ [N/B], and for any m ∈ [M ], the following bound holds with probability at least 1− δ3MNK :∥∥∥∥∥∥ i∑ t=(j−1)B+1 ∇fmσmk (t)(yk,0) ∥∥∥∥∥∥ ≤ ν √ 8(i− (j − 1)B) log 6MNKδ
+ (i− (j − 1)B) ‖∇Fm(yk,0)‖ . (104) For the summation that appear in (102), we use the techniques from Appendix D.4. For each l ∈ [N/B − 1], we can similarly decompose ∑M m=1 ∑lB t=1∇fmσmk (t)(yk,0) into pl + ql, where
pl := M∑ m=1 lB∑ t=1 ∇f̄σmk (t)(yk,0),
ql := M∑ m=1 lB∑ t=1 ∇fmσmk (t)(yk,0)−∇f̄σmk (t)(yk,0).
As done in Appendix D.4, we can follow the same steps and show a high-probability bound, which is a slightly different version of (88): with probability at least 1− Bδ3NK ,
‖pl‖ ≤ lBM ‖∇F (yk,0)‖+ ν √ N √ 8 log 6NK
Bδ . (105)
Similarly, for ‖ql‖ we can show a slight modification of (91): ‖ql‖ ≤ 4λ √ lBMN log 6N2K
Bδ , (106)
which holds with probability at least 1− Bδ3NK . Combining (105) and (106), with probability at least 1− 2Bδ3NK , we have∥∥∥∥∥ M∑ m=1 lB∑ t=1 ∇fmσmk (t)(yk,0) ∥∥∥∥∥ ≤ lBM ‖∇F (yk,0)‖+ ν√N √ 8 log 6NK Bδ
+ 4λ √ lBMN log 6N2K
Bδ . (107)
By applying the union bound, with probability at least 1− δK , the bound (104) holds for allm ∈ [M ] and i ∈ ⋃N/B j=1 [(j − 1)B + 1 : jB − 1], and the bound (107) holds for all l ∈ [N/B − 1].
We now substitute the bounds (104) and (107) to (101), (102), and (103) to get upper bounds for ‖rk,1‖, ‖rk,2‖, and ‖rk,3‖, respectively. For ‖rk,1‖ and ‖rk,3‖, we can apply the same calculations as in (70) and (72), modulo the fact that Assumption 3 is now implied by Assumption 4, with constants τ = λ and ρ = 1. We obtain
‖rk,1‖ ≤ 24LνN(B3/2 − 1)
11B
√ log 6MNK
δ + 3LN(B − 1) 5 (λ+ ‖∇F (yk,0)‖) , (108)
‖rk,3‖ ≤ 3L2νN(N −B)(B3/2 − 1)
2B
√ log 6MNK
δ
+ 2L2N(N −B)(B − 1)
5 (λ+ ‖∇F (yk,0)‖). (109)
For ‖rk,2‖, we have
‖rk,2‖ ≤ 7LB
5M N/B−1∑ l=1
( lBM ‖∇F (yk,0)‖+ ν √ N √ 8 log 6NK
Bδ + 4λ
√ lBMN log 6N2K
Bδ
)
= 7LB2
5 N/B−1∑ l=1 l  ‖∇F (yk,0)‖+ 7√8LνN1/2(N −B) 5M √ log 6NK Bδ
+ 28LλB3/2N1/2
5M1/2
N/B−1∑ l=1 √ l √log 6N2K Bδ
≤ 7LN(N −B) 10 ‖∇F (yk,0)‖+ 4LνN1/2(N −B) M
√ log 6NK
Bδ
+ 15LλN1/2(N3/2 −B3/2)
4M1/2
√ log 6N2K
Bδ . (110)
Recalling the definition rk := rk,1 + rk,2 − ηrk,3, we get an upper bound for ‖rk‖ from (108), (109), and (110):
‖rk‖ ≤ ‖rk,1‖+ ‖rk,2‖+ η ‖rk,3‖ ≤ Lν √ log 6MNK
δ
( 24N(B3/2 − 1)
11B + 4N1/2(N −B) M + 3ηLN(N −B)(B3/2 − 1) 2B ) + Lλ ( 3N(B − 1)
5 + 15N1/2(N3/2 −B3/2) 4M1/2
√ log 6N2K
Bδ + 2ηLN(N −B)(B − 1) 5
)
+ L ‖∇F (yk,0)‖ (
3N(B − 1) 5 + 7N(N −B) 10 + 2ηLN(N −B)(B − 1) 5
) . (111)
Recall again that we have K ≥ 7κ log(M2NK2) and η = log(M 2NK2)
µNK , so ηLN ≤ 1/7. Using this and N −B ≤ N , we can further simplify (111).
‖rk‖ ≤ Lν √ log 6MNK
δ
( 12N(B3/2 − 1)
5B + 4N1/2(N −B) M ) ︸ ︷︷ ︸
=:Φν
+ Lλ √ log 6N2K
Bδ
( 2N(B − 1)
3 + 15N1/2(N3/2 −B3/2) 4M1/2 ) ︸ ︷︷ ︸
=:Φλ
+ L ‖∇F (yk,0)‖ (
2N(B − 1) 3 + 7N(N −B) 10 ) ≤ LνΦν √ log 6MNK
δ + LλΦλ
√ log 6N2K
Bδ +
7LN2
5 ‖∇F (yk,0)‖ , (112)
which holds with probability at least 1− δK . The bound (112) holds for all k ∈ [K] with probability 1 − δ if we apply the union bound over k = 1, . . . ,K. Next, by (a + b + c)2 ≤ 3a2 + 3b2 + 3c2, we have
‖rk‖2 ≤ 3L2ν2Φ2ν log 6MNK
δ + 3L2λ2Φ2λ log
6N2K
Bδ +
147L2N4
25 ‖∇F (yk,0)‖2 , (113)
which also holds for all k ∈ [K] with probability at least 1− δ.
Getting a high-probability convergence rate. Given our high-probability bounds (112) and (113), we can substitute them to (62) and get
F (yk+1,0)− F (yk,0) ≤ (−ηN + η2LN2) ‖∇F (yk,0)‖2 + η2 ‖∇F (yk,0)‖ ‖rk‖+ η4L ‖rk‖2
≤ ( −ηN + η2LN2 + 7η 2LN2
5 +
147η4L3N4
25
) ‖∇F (yk,0)‖2
+ η2LνΦν √ log 6MNK
δ ‖∇F (yk,0)‖+ 3η4L3ν2Φ2ν log
6MNK
δ
+ η2LλΦλ √ log 6N2K
Bδ ‖∇F (yk,0)‖+ 3η4L3λ2Φ2λ log
6N2K
Bδ . (114)
The following terms in the RHS of (114) can be bounded using ab ≤ a 2 2 + b2 2 :
η2LνΦν √ log 6MNK
δ ‖∇F (yk,0)‖
= ( η1/2N1/2√
8 ‖∇F (yk,0)‖ )(√ 8η3/2LνΦν N1/2 √ log 6MNK δ )
≤ ηN 16 ‖∇F (yk,0)‖2 + 4η3L2ν2Φ2ν N log 6MNK δ , (115)
η2LλΦλ √ log 6N2K
Bδ ‖∇F (yk,0)‖
= ( η1/2N1/2√
8 ‖∇F (yk,0)‖ )(√ 8η3/2LλΦλ N1/2 √ log 6N2K Bδ )
≤ ηN 16 ‖∇F (yk,0)‖2 + 4η3L2λ2Φ2λ N log 6N2K Bδ . (116)
Substituting (115) and (116) to (114) results in
F (yk+1,0)− F (yk,0) ≤ ( −7
8 ηN +
12η2LN2
5 +
147η4L3ρ2N4
25
) ‖∇F (yk,0)‖2
+ η3L2ν2(4 + 3ηLN)Φ2ν
N log
6MNK δ + η3L2λ2(4 + 3ηLN)Φ2λ N log 6N2K Bδ . (117)
Again, we have K ≥ 7κ log(M2NK2) and η = log(M 2NK2)
µNK , so ηLN ≤ 1 7 . Since the inequality
− 78z + 12 5 z 2 + 14725 z 4 ≤ − 12z holds on z ∈ [0, 1 7 ], we have
−7 8 ηN +
12η2LρN2
5 +
147η4L3ρ3N4
25 ≤ −1 2 ηN.
Substituting this inequality to (117), together with 4 + 3ηLN ≤ 317 < 9 2 , yields
F (yk+1,0)− F (yk,0) ≤ − ηN
2 ‖∇F (yk,0)‖2 + 9η3L2ν2Φ2ν 2N log 6MNK δ
+ 9η3L2λ2Φ2λ
2N log
6N2K
Bδ .
We now recall that F is µ-PŁ, so ‖∇F (yk,0)‖2 ≥ 2µ(F (yk,0)− F ∗):
F (yk+1,0)− F ∗ ≤ (1− ηµN)(F (yk,0)− F ∗) + 9η3L2ν2Φ2ν
2N log
6MNK
δ
+ 9η3L2λ2Φ2λ
2N log
6N2K
Bδ . (118)
Recall that (118) holds for all k ∈ [K], with probability 1−δ. Therefore, by unrolling the inequality, and using ∑K−1 k=0 (1− ηµN)k ≤ 1 ηµN , we get
F (yK,NB )− F ∗ ≤ (1− ηµN)K(F (y0)− F ∗) + 9η2L2ν2Φ2ν 2µN2 log 6MNK δ
+ 9η2L2λ2Φ2λ
2µN2 log
6N2K
Bδ . (119)
Recall that Φν := 12N(B3/2−1) 5B + 4N1/2(N−B) M and Φλ := 2N(B−1) 3 + 15N1/2(N3/2−B3/2) 4M1/2 , hence
Φ2ν ≤ 288N2(B3/2 − 1)2
25B2 +
32N(N −B)2
M2 ,
Φ2λ ≤ 8N2(B − 1)2
9 +
225N(N3/2 −B3/2)2
2M .
Substituting these inequalities and also η = log(M 2NK2)
µNK gives
F (yK,NB )− F ∗
≤ F (y0)− F ∗
M2NK2 +
9L2ν2 log 6MNKδ log 2(M2NK2)
2µ3N4K2
( 288N2(B3/2 − 1)2
25B2 +
32N(N −B)2
M2 ) + 9L2λ2 log 6N 2K Bδ log 2(M2NK2)
2µ3N4K2
( 8N2(B − 1)2
9 +
225N(N3/2 −B3/2)2
2M ) = F (y0)− F ∗
M2NK2 + Õ
( L2
µ3
( ν2B
N2K2 +
ν2
M2NK2 +
λ2B2
N2K2 +
λ2
MK2
)) ,
with probability at least 1− δ. This finishes the proof.

heading: D.6 A GENERALIZED VECTOR-VALUED HOEFFDING-SERFLING INEQUALITY
text: We extend the vector-valued Hoeffding-Serfling inequality proved in Schneider (2016) to account for the mean of multiple independent without-replacement sums.
Lemma 8. Suppose there are MN vectors {vmi } M,N m=1,i=1 ∈ Rd that satisfy ‖vmi − v̄m‖ ≤ ν for m ∈ [M ], where v̄m := 1N ∑N i=1 v m i . Consider M independently and uniformly sampled permutations σ1, . . . , σM ∼ Unif(SN ). For any n ≤ N −1, with probability at least 1−δ, we have∥∥∥∥∥ 1Mn M∑ m=1 n∑ i=1 vmσm(i) − 1 M M∑ m=1 v̄m ∥∥∥∥∥ ≤ ν √ 8(1− n−1N ) log 2 δ Mn . (120)
Proof. The proof is an extension of Theorem 2 of Schneider (2016) which proves the M = 1 case for vectors in smooth separable Banach spaces. We prove our extended concentration inequality for Rd, but we note that the proof technique can be applied directly to general smooth separable Banach spaces, as done in Schneider (2016). Below, we state a special case of Theorem 3 of Pinelis (1992) and Theorem 3.5 of Pinelis (1994), because this Rd case serves our purpose.
Lemma 9 (Pinelis (1992; 1994)). Suppose that a sequence of random variables {xj}j≥0 is a martingale taking values in Rd, and ∑∞ j=1 ess sup ‖xj − xj−1‖
2 ≤ c2 for some c > 0. Then, for λ > 0,
P ( sup{‖xj‖ : j ≥ 0} ≥ λ ) ≤ 2 exp ( − λ 2
2c2
) .
The proof of Lemma 8 proceeds by defining a sequence of random variables {xj}, showing that it is a martingale, and applying Lemma 9 to prove our concentration bound. For m ∈ [M ], define index functions km : N ∪ {0} → [0 : n] in the following way:
km(j) := max{0,min{n, j − (m− 1)n}} =  0 if j ≤ (m− 1)n, j − (m− 1)n if (m− 1)n+ 1 ≤ j ≤ mn, n if j ≥ mn+ 1.
Using these index functions, we introduce the following sequence of random variables {xj}:
xj := M∑ m=1
1
N − km(j) km(j)∑ i=1 (vmσm(i) − v̄ m),
and show that this is a martingale, i.e.,
E[xj | x1, . . . ,xj−1] = xj−1 (121) for all j ≥ 1. Notice first that by definition of km’s we have xMn = xMn+1 = xMn+2 = . . . , so (121) is trivially satisfied for all j > Mn. Next, for any j satisfying (l − 1)n + 1 ≤ j ≤ ln where l ∈ [M ], we have
xj = l−1∑ m=1 1 N − n n∑ i=1 (vmσm(i) − v̄ m) +
1
N − kl(j) kl(j)∑ i=1 (vlσl(i) − v̄ l)
= xj−1 +
( 1
N − kl(j) − 1 N − kl(j) + 1 ) kl(j)−1∑ i=1 (vlσl(i) − v̄ l)
+ 1
N − kl(j) (vlσl(kl(j)) − v̄ l) (122)
Now note that for any k ∈ [N − 1], we have
E[vlσl(k) − v̄ l | σl(1), . . . , σl(k − 1)] =
1
N − k + 1 N∑ i=k (vlσl(i) − v̄ l)
= − 1 N − k + 1 k−1∑ i=1 (vlσl(i) − v̄ l),
where the last equality used ∑N i=1(v l σl(i) − v̄l) = 0. Using applying this fact to (122) and noting
1 N−kl(j) − 1 N−kl(j)+1 = 1 (N−kl(j))(N−kl(j)+1) ,
E[xj | x1, . . . ,xj−1] = xj−1 + 1
(N − kl(j))(N − kl(j) + 1) kl(j)−1∑ i=1 (vlσl(i) − v̄ l)
+ 1 N − kl(j) E [ vlσl(kl(j)) − v̄ l | x1, . . . ,xj−1 ]
= xj−1,
hence proving that {xj}j≥0 is a martingale. We now apply Lemma 9 to our {xj}. For j such that (l − 1)n+ 1 ≤ j ≤ ln, notice from (122) that
(N − kl(j))(xj − xj−1) = 1
(N − kl(j) + 1) kl(j)−1∑ i=1 (vlσl(i) − v̄ l) + (vlσl(kl(j)) − v̄l)
= − 1 (N − kl(j) + 1)  N∑ i=kl(j) (vlσl(i) − v̄ l) + (vlσl(kl(j)) − v̄l), which leads to
(N − kl(j)) ‖xj − xj−1‖ ≤ νmin{kl(j)− 1, N − kl(j) + 1}
N − kl(j) + 1 + ν ≤ 2ν,
by the triangle inequality. From this, we get the bound c2 in the statement of Lemma 9:
∞∑ j=1 ess sup ‖xj − xj−1‖2 = Mn∑ j=1 ess sup ‖xj − xj−1‖2 ≤M n∑ k=1 4ν2 (N − k)2
= 4ν2M
(N − n)2 + 4ν2M N−1∑ k=N−n+1 1 k2 ≤ 4ν 2M (N − n)2 + 4ν2M(n− 1) (N − n)N = 4ν2Mn (N − n)2 ( 1− n− 1 N ) ,
where the second inequality used the inequality that ∑b k=a+1 1 k2 ≤ b−a a(b+1) (Serfling, 1974, Lemma 2.1). Now, applying Lemma 9 to {xj} with c2 = 4ν 2Mn (N−n)2 ( 1− n−1N ) gives
P ( ‖xMn‖ ≥ λ ) ≤ P ( sup{‖xj‖ : j ≥ 0} ≥ λ ) ≤ 2 exp ( − λ 2(N − n)2
8ν2Mn ( 1− n−1N )) . (123) Recall from the definition of {xj} that
xMn = 1
N − n M∑ m=1 n∑ i=1 (vmσm(i) − v̄ m).
Substituting λ = Mn N−n to (123) gives
P (∥∥∥∥∥ 1Mn M∑ m=1 n∑ i=1 (vmσm(i) − v̄ m) ∥∥∥∥∥ ≥ ) ≤ 2 exp ( − Mn 2 8ν2 ( 1− n−1N )) , which finishes the proof.
D.7 HOW CAN WE AVOID UNIFORM BOUNDS OVER Rd IN OUR ASSUMPTIONS?
In Section 2, we introduced Assumptions 2, 3, and 4 on the intra- and inter-machine deviation. The assumptions required that inequalities such as ‖∇fmi (x)−∇Fm(x)‖ ≤ ν hold for all x ∈ Rd. In this subsection, we discuss more on this strong requirement “entire Rd.”
In fact, the entire-Rd requirement is posed in our assumptions to simplify the exposition of the main results, and is not strictly necessary. One can easily check from our proofs that the assumptions are only applied to the beginning iterates xk,0 (for minibatch RR) or yk,0 (for local RR) of epochs. Hence, if these iterates lie in a bounded set, then the constants ν, τ , ρ, and λ may become much smaller, depending on problem instances. Actually, if we explicitly assume that the iterates lie in a compact set S,12 then Assumptions 2–4 are even guaranteed to hold for some constants; e.g., for Assumption 2, we can choose
ν := max m∈[M ],i∈[N ],x∈S
‖∇fmi (x)−∇Fm(x)‖ ,
since the maximum always exists.
However, assuming that the iterates lie in a specific set S can be problematic because the distance that the iterates travel depend on the objective functions. One cannot know a priori if all iterates will stay in a fixed set S; hence, explicitly assuming bounded iterates should be avoided.
12This bounded iterates assumption is indeed used in some existing results such as Haochen & Sra (2019); Nagaraj et al. (2019); Rajput et al. (2020); Ahn et al. (2020).
Then, a natural question is whether we can prove bounded iterates under some reasonable conditions, instead of assuming it. We point out that this can be done by applying the technique developed in Ahn et al. (2020, Theorem 1) to our upper bound theorems. Using the technique, a modified version of our Theorem 1 can be written as follows:
Theorem 10 (Best-iterate version of Theorem 1). Suppose that minibatch RR has parameters satisfying Assumption 1. Assume that all local component functions fmi are L-smooth, the global objective function F is µ-PŁ, and the set of global minima of F is nonempty and compact. Consider running the algorithm using step-size η = B log(MNK
2) µNK and initialization x1,0 := x0, for epochs
K ≥ 6κ log(MNK2). Then, with probability at least 1− δ,
min k∈[K+1]
F (xk,0)− F ∗ ≤ F (x0)− F ∗
MNK2 + Õ
( L2
µ3 ν2 MNK2
) ,
where the constant ν <∞ is defined as
ν := sup x:F (x)≤F (x0) max i∈[N ] max m∈[M ]
‖∇fmi (x)−∇Fm(x)‖ . (124)
In the theorem, we used xK+1,0 to denote the last iterate of the algorithm xK,NB . Theorem 10 differs from Theorem 1 in three aspects: 1) it considers the best-iterate, not the last-iterate; 2) it additionally assumes that the set of global minima of F is nonempty and compact, which always holds if F is strongly convex; and 3) it does not rely on Assumption 2, but instead “proves” it for the F (x0)-sublevel set of F (124). Note that the constant ν (124) can be much smaller than the uniform bound required to make Assumption 2 hold for the entire Rd. For Theorems 2, 6, and 7, we can also apply similar techniques to prove best-iterate bounds with smaller intra- and inter-machine deviation constants ν, τ , ρ, and λ; we omit the precise statements.
We conclude this subsection with the proof of Theorem 10.
Proof. The proof follows that of Ahn et al. (2020, Theorem 1).
Existence of ν. We first show the existence of ν < ∞ (124). The global objective function F is µ-PŁ. If we denote the set of global minima of F as X∗, the set X∗ is nonempty and compact by assumption. Then, by Karimi et al. (2016, Theorem 2) µ-PŁ functions satisfy quadratic growth, i.e., denoting by x∗ the closest global minimum in X∗ to the point x,
F (x)− F ∗ ≥ 2µ ‖x− x∗‖2 .
Define the sublevel set S := {x | F (x) ≤ F (x0)}. Due to the quadratic growth property, we have F (x0)− F ∗ ≥ F (x)− F ∗ ≥ 2µ ‖x− x∗‖2 for all x ∈ S. This implies that
S := {x ∈ Rd | F (x) ≤ F (x0)} ⊂ { x ∈ Rd | ‖x− x∗‖2 ≤ F (x0)− F ∗
2µ
} .
Since we assumed that X∗ is compact, S is also bounded, and hence compact. Now, for anym ∈ [M ] and i ∈ [N ], ‖∇fmi (x)−∇Fm(x)‖ is a continuous function on a compact set S, so there must exist a constant νmi <∞ such that ‖∇fmi (x)−∇Fm(x)‖ ≤ νmi for all x ∈ S. Taking the maximum of νmi over all m and i gives ν.
Proving the best-iterate bound. With the constant ν (124), if all the iterates {xk,0}k∈[K+1] stay within the sublevel set S := {x | F (x) ≤ F (x0)}, one can consider Assumption 2 to be true with constant ν. From this observation, we consider two cases:
1. All the iterates {xk,0}k∈[K+1] stay in the sublevel set S.
2. There exists an iterate xk,0 /∈ S.
In fact, the first case can be proven by exactly the same steps as Theorem 1, described in Appendix D.2.
For the second case, suppose that there exists an iterate xk,0 that escapes the sublevel set S. Let k′ ∈ {2, . . . ,K + 1} be the first such k. Then, since xk′−1,0 is still in S, it follows from (50) in Appendix D.2 that we have
F (xk′,0)−F ∗ ≤ (1−ηµN)(F (xk′−1,0)−F ∗) + 15η3L2ν2(N3/2 −B3/2)2
MN log
2NK
Bδ . (125)
However, the fact that xk′,0 /∈ S and xk′−1,0 ∈ S implies
F (xk′,0) > F (x0) ≥ F (xk′−1,0). (126)
Combining the two bounds (125) and (126), we get
0 < −ηµN(F (xk′−1,0)− F ∗) + 15η3L2ν2(N3/2 −B3/2)2
MN log
2NK
Bδ ,
which implies
min k∈[K+1]
F (xk,0)− F ∗ ≤ F (xk′−1,0)− F ∗ < 15η2L2ν2(N3/2 −B3/2)2
µMN2 log
2NK
Bδ .
Substituting η = log(MNK 2)
µNK 13 gives the desired bound and finishes the proof.

heading: E PROOF OF LOWER BOUND FOR MINIBATCH RR (THEOREM 3)
text: For Theorem 3, we consider three step-size ranges and do case analysis for each of them. We construct functions for each corresponding step-size regime such that the convergence of minibatch RR is “slow” for the functions on their corresponding step-size regime. The final lower bound is the minimum among the lower bounds obtained for the three regimes. More concretely, we will construct three one-dimensional functions F1(x), F2(x), and F3(x) satisfying L-smoothness (1), µ-PŁ condition (2), and Assumption 2 such that14
• Minibatch RR on F1(x) with η ≤ BµNK and initialization x0 = ν µ results in
E[F1(xK,NB )] = Ω ( ν2
µ
) .
• Minibatch RR on F2(x) with η ≥ BµNK and η ≤ B 513LN and initialization x0 = 0 results in
E[F2(xK,NB )] = Ω (
ν2
µMNK2
) .
Note that the step-size range requires K ≥ 513κ, hence this lower bound occurs only in the “large-epoch” regime, i.e., K & κ.
• Minibatch RR on F3(x) with η ≥ BµNK and η ≥ B 513LN and initialization x0 = 0 results in
E[F3(xK,NB )] = Ω ( ν2
µMNK
) .
Then, the three dimensional function F ([x, y, z]>) = F1(x) + F2(y) + F3(z) will show bad convergence in any step-size regime. Furthermore,
µI min(∇2F1,∇2F2,∇2F3)I ∇2F max(∇2F1,∇2F2,∇2F3)I LI,
13Recall that this is different from η = B log(MNK 2)
µNK in the theorem statement, because for the proofs, we
consider an equivalent “rescaled” version of minibatch RR defined in the beginning of Appendix D.2. 14In fact, the functions constructed in this theorem are µ-strongly convex, which is stronger than µ-PL required in Definition 1.
that is, if F1, F2 and F3 are µ-strongly convex and L-smooth, then so is F . Moreover, since the component functions in each coordinate are designed to satisfy Assumption 2 with ν, the resulting three dimensional function F also satisfies Assumption 2 with √ 3ν.
Since the final lower bound is the minimum among the lower bounds obtained in the step-size ranges, the lower bound becomes Ω ( ν2
µMNK2
) ifK ≥ 513κ, and Ω ( ν2
µMNK
) ifK < 513κ (in which case
the second step-size range does not exist).
In the subsequent subsections, we prove the lower bounds for F1, F2, and F3 separately.
E.1 LOWER BOUND FOR η ≤ BµNK
Consider the case where every function at every machine is the same: for all i ∈ [N ] and m ∈ [M ], fmi (x) := µx2 2 . Hence, F1(x) = µx2 2 .
Let xk,0 and xk,NB denote the iterates where the k-th epoch starts and ends respectively. Then,
xk+1,0 = xk,NB = (1− ηµ)NB xk,0.
Initializing at x1,0 = νµ and unrolling this for K epochs, we get
xK,NB = (1− ηµ)NKB · ν µ ≥ ( 1− B NK )NK B · ν µ ≥ ν 4µ ,
since N ≥ 2, K ≥ 1, and B divides N . Hence, F1(xK,NB ) = Ω( ν2 µ ).
E.2 LOWER BOUND FOR η ≥ BµNK AND η ≤ B 513LN
For most part of this subsection, we consider iterates within a single epoch, and hence we will omit the subscripts denoting epochs. Let x0 denote the iterate at the beginning of the epoch, and xi denote the iterate after the i-th communication round in that epoch. In our construction, each machine will have the same set of component functions, that is, there will be no inter-machine deviation. We therefore omit the superscriptm from the local component functions fmi . The function we construct for the lower bound and its component functions are as follows:
F2(x) := 1
N  N2∑ i=1 f+1(x) + N∑ i=N2 +1 f−1(x)  , where f+1(x) := (L1x≤0 + µ1x>0) x2
2 + νx, and
f−1(x) := (L1x≤0 + µ1x>0) x2
2 − νx
Note that the function F2(x) = (L1x≤0 + µ1x>0)x 2
2 is µ-strongly convex and L-smooth with minimizer at 0, and also satisfies Assumption 2.
Let σm be a random permutation of N2 +1’s and N 2 −1’s. Then, machine m computes gradients on f−1 and f+1 in the order given by σm. Let σmj denote the j-th ordered element of σ m. Then,
∇fσmj (x) = (L1x≤0 + µ1x>0)x+ νσ m j .
Hence, the last iterate of an epoch, xN B , is given by
xN B − x0 = N B−1∑ i=0 − η MB M∑ m=1 (i+1)B∑ j=iB+1 ∇fσmj (xi)  =
N B−1∑ i=0 − η MB M∑ m=1 (i+1)B∑ j=iB+1 ((L1xi≤0 + µ1xi>0)xi + νσ m j ) 
= N B−1∑ i=0 − η MB M∑ m=1 (i+1)B∑ j=iB+1 (L1xi≤0 + µ1xi>0)xi  (Since∑Nj=1 σmj = 0) = −η
N B−1∑ i=0 (L1xi≤0 + µ1xi>0)xi. (127)
Thus, E[xN B − x0] = −η ∑N B−1 i=0 E[(L1xi≤0 + µ1xi>0)xi]. We want to prove that E[xNB ] keeps increasing over an epoch, that is E[xN B − x0] > 0 when x0 is close enough to the minimizer 0.
For this, we first consider the case where the first iterate x0 of the epoch satisfies x0 ≥ 0. The x0 < 0 case will be considered later. For the case x0 ≥ 0, we will show that whenever x0 is small, the expected amount of update made in the (i+ 1)-th iteration, E[(L1xi≤0 +µ1xi>0)xi], is negative in the first half of the epoch and not too big in the second half.
We use the following lemmas, proven in Appendices F.1 and F.2, respectively. Lemma 11. For x0 ≥ 0, 0 ≤ i ≤ b N2B c, η ≤ B 513LN , and L µ ≥ 7695,
E[(L1xi≤0 + µ1xi>0)xi] ≤ 6
7 Lx0 −
ηLν
1536
√ i
MB .
Lemma 12. For x0 ≥ 0, 0 ≤ i ≤ NB − 1, and η ≤ B 513LN , E[(L1xi≤0 + µ1xi>0)xi] ≤ µ ( 1 + 513iηL
512
) x0 + 513ηµν
512
√ i
MB .
The key intuition is that, for Lµ big enough, we can use the lemmas above in (127) to get
E[xN B − x0] = −η N B−1∑ i=0 (L1xi≤0 + µ1xi>0)xi
≈ Ω ( η N
B ηLν
√ N/B
MB
) ,
whenever |x0| is small. Multiplying the above by K (for K epochs) will give us the required lower bound (up to factors of Lµ ). We will make this approximate calculation precise in the rest of the proof.
Using the two lemmas above in (127), we get that E[xN
B − x0]
= − η N B−1∑ i=0 E[(L1xi≤0 + µ1xi>0)xi] = − η b N2B c∑ i=0 E[(L1xi≤0 + µ1xi>0)xi]− η N B−1∑
i=b N2B c+1
E[(L1xi≤0 + µ1xi>0)xi]
≥ − η b N2B c∑ i=0
( 6
7 Lx0 −
ηLν
1536
√ i
MB
)
− η N B−1∑
i=b N2B c+1
( µ ( 1 + 513iηL
512
) x0 + 513ηµν
512
√ i
MB
) . (128)
Since iηL ≤ ηLNB ≤ 1 513 , µ ≤ L 7695 , and N/B ≥ 2, the following bound holds: b N2B c∑ i=0 6 7 L+ N B−1∑
i=b N2B c+1
µ ( 1 + 513iηL
512
)
≤ (⌊ N
2B
⌋ + 1 ) 6L
7 +
( N B − ⌊ N 2B ⌋ − 1 ) L 7695 ( 1 + 1 512 ) ≤ 6LN
7B +
LN 7680B ≤ 7LN 8B . (129)
Also, note that b N2B c ≥ N 3B whenever N/B ≥ 2. We have
b N2B c∑ i=0 √ i ≥ ∫ b N2B c 0 √ tdt = 2 3 (⌊ N 2B ⌋)3/2 ≥ 2 3 ( N 3B )3/2 = 2N3/2 9 √ 3B3/2 , (130)
N B−1∑
i=b N2B c+1
√ i ≤
∫ N B
b N2B c+1
√ tdt ≤ 2
3
[( N
B
)3/2 − (⌊ N
2B
⌋ + 1 )3/2]
≤ 2 3
[( N
B
)3/2 − ( N
2B
)3/2] = (2 √ 2− 1)N3/2
3 √ 2B3/2 ≤ N
3/2
2B3/2 . (131)
Substituting the bounds (129), (130), and (131) into (128), and using µ ≤ L7695 ,
E[xN B − x0] ≥ −
7ηLN
8B x0 +
η2LνN3/2
6912 √ 3M1/2B2 − 513η
2µνN3/2
1024M1/2B2
≥ −7ηLN 8B x0 + η2LνN3/2 56000M1/2B2 . (132)
For the other case x0 < 0, we have the following Lemma, which we prove in Appendix F.3: Lemma 13. If η ≤ B513NL and an epoch starts at x0 < 0, then
E[xN B | x0 < 0] ≥
( 1− 7ηLN
8B
) x0.
Further, if the first epoch of the algorithm is initialized at 0, then for any starting iterate x0 of any following epoch, we have P(x0 ≥ 0) ≥ 1/2.
Using (132) and Lemma 13 we get
E[xN B ] = P(x0 ≥ 0)E[xN B | x0 ≥ 0] + P(x0 < 0)E[xN B | x0 < 0] ≥ P(x0 ≥ 0) ((
1− 7ηLN 8B
) x0 + η2LνN3/2
56000M1/2B2
) + P(x0 < 0) ( 1− 7ηLN
8B
) x0
≥ (
1− 7ηLN 8B
) x0 + η2LνN3/2
112000M1/2B2 .
Thus far, we have characterized the expected per-epoch update, starting from the initial iterate x0 and iterating until the last iterate xN
B of the epoch. Now recall that we run the algorithm for K epochs.
Using xk,i to denote the i-th iterate of the k-th epoch, we get a lower bound on the expectation of the last iterate xk,NB if we initialize at x1,0 = 0:
E[xK,NB ] ≥ ( 1− 7ηLN 8B )K x1,0 + η2LνN3/2 112000M1/2B2 K−1∑ k=0 ( 1− 7ηLN 8B )k
= η2LνN3/2
112000M1/2B2
1− ( 1− 7ηLN8B )K
7ηLN 8B
= ηνN1/2
98000M1/2B
( 1− ( 1− 7ηLN
8B
)K)
≥ ηνN 1/2
98000M1/2B
( 1− ( 1− 7L
8µK
)K) . (Since η ≥ BµNK )
Note that since Lµ ≥ 7695 and K ≥ 513L µ (which is implied by B µNK ≤ η ≤ B 513LN ),
1− (
1− 7L 8µK
)K ≥ 1− e− 7L 8µ ≥ 1− e−6733 ≈ 1.
Hence, we get from η ≥ BµNK that
E[xK,NB ] = Ω ( ηνN1/2
M1/2B
) = Ω ( ν
µM1/2N1/2K
) ,
and by Jensen’s inequality, we finally have
E[F (xK,NB )] ≥ 1 2 E[µx2 K,NB ] = Ω(µE[xK,NB ] 2) = Ω
( ν2
µMNK2
) .
E.3 LOWER BOUND FOR η ≥ BµNK AND η ≥ B 513LN
Similar to earlier parts of the proof, here as well, each machine will have the same component functions, that is, there will be no inter-machine deviation. The proof uses a similar construction as Safran & Shamir (2020; 2021):
F3(x) := 1
N  N2∑ i=1 f+1(x) + N∑ i=N2 +1 f−1(x)  , where f+1(x) := Lx2
2 + νx, and f−1(x) :=
Lx2
2 − νx.
Hence, F3(x) = Lx 2
2 , and has its minimizer at 0.
We first compute the expected “progress” over a given epoch. For simplicity, let us omit the subscript for epochs for now. Let x0 denote the iterate at the beginning of the epoch and xi denote the iterate after the i-th communication round in that epoch. For a given epoch, let σm be the permutation of N 2 +1’s and N 2 −1’s, sampled by machine m. Then,
xN B = xN B−1 − η MB M∑ m=1 N∑ j=(NB−1)B+1 (LxN B−1 + νσmj )
= (1− ηL)xN B−1 − ην MB M∑ m=1 N∑ j=(NB−1)B+1 σmj
= · · ·
= (1− ηL)NB x0 − ην
MB N B∑ i=1 (1− ηL)NB−i M∑ m=1 iB∑ j=(i−1)B+1 σmj .
For the rest of this proof, x2i refers to the square of the i-th iterate. Then,
E[x2N B
] = (1− ηL) 2NB x20 − 2ην(1− ηL)NB x0
MB E  NB∑ i=1 (1− ηL)NB−i M∑ m=1 iB∑ j=(i−1)B+1 σmj  + η2ν2
M2B2 E   NB∑ i=1 (1− ηL)NB−i M∑ m=1 iB∑ j=(i−1)B+1 σmj 2 
= (1− ηL) 2NB x20 + η2ν2
M2B2 E   NB∑ i=1 (1− ηL)NB−i M∑ m=1 iB∑ j=(i−1)B+1 σmj 2  , (133)
where we used the fact that E[σmj ] = 0. Further, because σm and σm ′
are independent and identically distributed for different m and m′, we get that
E   NB∑ i=1 (1− ηL)NB−i M∑ m=1 iB∑ j=(i−1)B+1 σmj 2 
=E   M∑ m=1 N B∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σmj 2 
= M∑ m=1 E   NB∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σmj 2 
+ ∑ m6=m′ E  NB∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σmj E  NB∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σm ′ j  =ME
  NB∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σ1j 2  ,
where the last equality used the fact that E[σmj ] = 0 for all m ∈ [M ] and i ∈ [N ], and that σm are identically distributed. Since we only consider the permutation σ1 (i.e., the one for machine 1) from now on, we henceforth omit the superscript. Substituting this to (133) gives
E[x2N B
] = (1− ηL) 2NB x20 + η2ν2
MB2 E   NB∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σj 2 
︸ ︷︷ ︸ =:Φ
. (134)
From (134), we have calculated the per-epoch expected update. Recall that we run the algorithm for K epochs. Using xk,i to denote the i-th iterate of the k-th epoch, we get a lower bound on the expectation of the last iterate xk,NB squared:
E[x2 k,NB
] = (1− ηL) 2NKB x21,0 + η2ν2
MB2 Φ K−1∑ k=0 (1− ηL) 2NkB ≥ η 2ν2 MB2 Φ, (135)
where the inequality used x1,0 = 0 and ∑K−1 k=0 (1 − ηL) 2Nk B ≥ (1 − ηL)0 = 1. Next, we analyze the expectation term, i.e., Φ, defined in (134).
Φ := E   NB∑ i=1 (1− ηL)NB−i iB∑ j=(i−1)B+1 σj 2 
= N∑ j=1 (1− ηL)2(NB−b j−1 B c−1)σ2j + ∑ j 6=j′ (1− ηL)NB−b j−1 B c−1(1− ηL)NB−b j′−1 B c−1E[σjσj′ ]
Noting that σ2j = 1 and E[σjσj′ ] = − 1N−1 , we get
Φ = B N B−1∑ j=0 (1− ηL)2j − 1 N − 1 ∑ j 6=j′ (1− ηL)NB−b j−1 B c−1(1− ηL)NB−b j′−1 B c−1
= B N B−1∑ j=0 (1− ηL)2j − 1 N − 1   N∑ j=1 (1− ηL)NB−b j−1 B c−1 2 − N∑ j=1 (1− ηL)2(NB−b j−1 B c−1) 
= B N B−1∑ j=0 (1− ηL)2j − 1 N − 1
B2 NB−1∑
j=0
(1− ηL)j 2 −B NB−1∑
j=0
(1− ηL)2j  = B2
N − 1 N B N B−1∑ j=0 (1− ηL)2j − NB−1∑ j=0 (1− ηL)j 2 
= B2(NB − 1) N − 1 (1 + 1N B − 1 ) N B−1∑ j=0 (1− ηL)2j − 1 N B − 1 NB−1∑ j=0 (1− ηL)j 2  . (136)
Note that the term in the parenthesis is exactly the right hand side of Equation (23) in Safran & Shamir (2020) modulo n and α replaced with NB and ηL, respectively. Hence, by Lemma 1 of Safran & Shamir (2020), we have(
1 + 1
N B − 1 ) N B−1∑ j=0 (1− ηL)2j − 1 N B − 1 NB−1∑ j=0 (1− ηL)j 2≥ c ·min{ 1 ηL , η2L2N3 B3 } ,
(137)
for some universal constant c > 0. Using the fact that η ≥ B513LN , it is easy to check that the RHS of (137) is lower-bounded by c ′
ηL , where c ′ > 0 is a universal constant. Combining (135), (136), and
(137) gives
E[x2 k,NB
] ≥ η 2ν2 MB2 · B2(NB − 1) N − 1 · c ′ ηL = c′ην2 LMB N −B N − 1 .
Since 2B divides N , we have B ≤ N/2. Since N ≥ 2, we have N−BN−1 ≥ N 2(N−1) ≥ 1 2 . Using this and the fact that η ≥ BµNK , we get
E[F3(xk,NB )] = L 2 E[x2 k,NB ] ≥ c
′ν2
4µMNK .

heading: F PROOFS OF HELPER LEMMAS FOR APPENDIX E

heading: F.1 PROOF OF LEMMA 11
text: First, if i = 0 then the lemma trivially holds, because x0 ≥ 0 gives
E[(L1x0≤0 + µ1x0>0)x0] = µx0 ≤ 6
7 Lx0.
The inequality holds because Lµ ≥ 7695.
For the rest of the proof, we consider the case 1 ≤ i ≤ N2B . By the law of total expectation we have
E[(L1xi≤0 + µ1xi>0)xi] = P  M∑ m=1 iB∑ j=1 σmj > 0 E (L1xi≤0 + µ1xi>0)xi ∣∣∣∣∣∣ M∑ m=1 iB∑ j=1 σmj > 0  + P
 M∑ m=1 iB∑ j=1 σmj ≤ 0 E (L1xi≤0 + µ1xi>0)xi ∣∣∣∣∣∣ M∑ m=1 iB∑ j=1 σmj ≤ 0  ≤ P
 M∑ m=1 iB∑ j=1 σmj > 0 LE xi ∣∣∣∣∣∣ M∑ m=1 iB∑ j=1 σmj > 0  + P
 M∑ m=1 iB∑ j=1 σmj ≤ 0 µE xi ∣∣∣∣∣∣ M∑ m=1 iB∑ j=1 σmj ≤ 0  , (138)
where the last inequality used the fact that (L1t≤0 + µ1t>0)t ≤ Lt and (L1t≤0 + µ1t>0)t ≤ µt for any t ∈ R.
Define E := ∑M m=1 ∑iB j=1 σ m j . We handle each of the two expectations in (138) separately. We first bound E [xi|E > 0].
E [xi|E > 0] = E x0 − i−1∑ j=0 η MB M∑ m=1 (j+1)B∑ k=jB+1 (νσmk + (L1xj≤0 + µ1xj>0)xj) ∣∣∣∣∣∣E > 0 
= E x0 − i−1∑ j=0 η MB M∑ m=1 (j+1)B∑ k=jB+1 (νσmk + (L1xj≤0 + µ1xj>0)(xj − x0)) ∣∣∣∣∣∣E > 0 
+ E − i−1∑ j=0 η MB M∑ m=1 (j+1)B∑ k=jB+1 (L1xj≤0 + µ1xj>0)x0 ∣∣∣∣∣∣E > 0 
= x0E 1− η i−1∑ j=0 (L1xj≤0 + µ1xj>0) ∣∣∣∣∣∣E > 0 − ην MB E [E|E > 0]
− η i−1∑ j=0 E [ (L1xj≤0 + µ1xj>0)(xj − x0) ∣∣E > 0]
≤ x0E 1− η i−1∑ j=0 (L1xj≤0 + µ1xj>0) ∣∣∣∣∣∣E > 0 − ην MB E [E|E > 0]
+ ηL i−1∑ j=0 E[|xj − x0| | E > 0]. (139)
Next, we use the following lemma to bound the conditional expectations that arise in (139). This lemma is proven in Appendix F.4 and it may be of independent interest to readers.
Lemma 14. For m ∈ [M ], let σm be a random permutation of N2 +1’s and N 2 −1’s. Then, for any i ≤ N2 and k ≤ B 2 , we have
1
64
(√ i
M + √ k
) ≤ E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  .
Furthermore, for any 0 ≤ i ≤ N and 0 ≤ k ≤ N satisfying i+ k ≤ N , we have
E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  ≤√ i M + √ k.
Lastly, for any 0 ≤ i ≤ N2 and 0 ≤ k ≤ B 2 satisfying i+ k ≥ 1, we have
P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj > 0  = P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj < 0  ≥ 1 6 .
Lemma 14 implies that 1ME [E|E > 0] ∈ [ 1 64 √ iB M , √ iB M ] and P(E > 0) = P(E < 0) ≥ 1/6. From this, we get
ην MB E [E|E > 0] ≥ ην 64
√ i
MB , (140)
E[|xj − x0| | E > 0] ≤ E[|xj − x0|] P(E > 0) ≤ 6E[|xj − x0|]. (141)
Also, since η ≤ BLN we have
1− iηµ ≥ 1− η i−1∑ j=0 (L1xj≤0 + µ1xj>0) ≥ 1− ηLN B ≥ 0,
which implies that
x0E 1− η i−1∑ j=0 (L1xj≤0 + µ1xj>0) ∣∣∣∣∣∣E > 0  ≤ (1− iηµ)x0. (142)
Substituting (140), (141), and (142) to (139), we obtain
E [xi|E > 0] ≤ (1− iηµ)x0 − ην
64
√ i
MB + 6ηL i−1∑ j=0 E[|xj − x0|]. (143)
Next, we have the following lemma that we can apply to E[|xj − x0|]. Proof of Lemma 15 can be found in Appendix F.5. Lemma 15. For x0 ≥ 0, 0 ≤ i ≤ NB − 1 and η ≤ B 513LN ,
E[|xi − x0|] ≤ 513
512 ην
√ i
MB +
513 512 iηLx0.
Applying this lemma to (143), we get
E [xi|E > 0] ≤ (1− iηµ)x0 − ην
64
√ i
MB +
1539η2Lν 256 √ MB i−1∑ j=0 √ j + 1539η2L2x0 256 i−1∑ j=0 j
≤ (1− iηµ)x0 − ην
64
√ i
MB +
513i3/2η2Lν
128 √ MB
+ 1539i2η2L2
512 x0
= ( 1− iηµ+ 1539i 2η2L2
512
) x0 − ( 1
64 − 513iηL 128
) ην √ i
MB ≤ (
1− iηµ+ 3iηL 512
) x0 − ην
128
√ i
MB . (144)
where we got the last inequality by using the fact that iηL ≤ ηLNB ≤ 1 513 , which follows from η ≤ B513LN . So far, we have obtained an upper bound for E [xi|E > 0]. Recall that there is another conditional expectation in (138) that we want to bound, namely E [xi|E ≤ 0]. We bound it below, using the tools developed so far. For i ≤ N2B ,
E [xi|E ≤ 0] = x0 + E [xi − x0 | E ≤ 0] ≤ x0 + E [|xi − x0| | E ≤ 0]
≤ x0 + E [|xi − x0|] P(E ≤ 0)
≤ x0 + 6E [|xi − x0|] (Using Lemma 14)
≤ x0 + 1539ην
256
√ i
MB + 1539iηLx0 256
(Using Lemma 15)
≤ ( 1 + 1539iηL
256
) x0 + 1539ην
256
√ i
MB (145)
Using (144) and (145) in (138), we get that for i ≤ N2B :
E[(L1xi≤0 + µ1xi>0)xi]
≤P (E > 0)LE [xi|E > 0] + P (E ≤ 0)µE [xi|E ≤ 0]
≤P (E > 0)L (( 1− iηµ+ 3iηL
512
) x0 − ην
128
√ i
MB
)
+ P (E ≤ 0)µ (( 1 + 1539iηL
256
) x0 + 1539ην
256
√ i
MB
) . (146)
From Lemma 14, note that 16 ≤ P (E > 0) ≤ 5 6 and 1 6 ≤ P (E ≤ 0) ≤ 5 6 . We use these inequalities, along with iηL ≤ ηLNB ≤ 1 513 and L µ ≥ 7695, to bound the terms appearing in (146).
P (E > 0)L (
1− iηµ+ 3iηL 512
) x0 + P (E ≤ 0)µ ( 1 + 1539iηL
256
) x0
≤ 5 6 L
( 1 + 1
87552
) x0 + 5
6 · L 7695
( 1 + 3
256
) x0 ≤ 6
7 Lx0. (147)
We also have
− P (E > 0) ηLν 128
√ i
MB + P (E ≤ 0) 1539ηµν 256
√ i
MB
≤ − ηLν 768
√ i
MB +
2565ηµν
512
√ i
MB
≤ − ηLν 1536
√ i
MB , (148)
where we used the assumption Lµ ≥ 7695. Substituting (147) and (148) to (146), we get
E[(L1xi≤0 + µ1xi>0)xi] ≤ 6
7 Lx0 −
ηLν
1536
√ i
MB ,
as desired.

heading: F.2 PROOF OF LEMMA 12
text: E[(L1xi≤0 + µ1xi>0)xi] ≤ µE[xi] = µx0 + µE[xi − x0] ≤ µx0 + µE[|xi − x0|]
≤ µx0 + 513
512 iηLµx0 +
513 512 ηµν
√ i
MB . (Using Lemma 15)

heading: F.3 PROOF OF LEMMA 13
text: We consider iterates within a single epoch, and hence we will omit the subscripts denoting epochs. In our construction, each machine has the same set of component functions, that is, there will be no inter-machine deviation. We therefore omit the superscript m from the local component functions. Consider the function
G2(x) := 1
N  N2∑ i=1 g+1(x) + N∑ i=N2 +1 g−1(x)  , where g+1(x) := Lx2
2 + νx, and g−1(x) :=
Lx2
2 − νx.
Hence, G2(x) = Lx 2
2 . We will prove the lemma by coupling iterates corresponding to F2 and G2. In particular, we will perform minibatch RR on F2 and G2 such that both start the given epoch at x0 and all the corresponding machines use the same random permutations. Let xi,F be the iterate after the i-th round of communication for F2 and xi,G be the iterate after the i-th round of communication for G2. We use mathematical induction to prove that xi,F ≥ xi,G for all i = 0, . . . , NB . After that, we will use this to prove our desired statement E[xN
B ,F | x0 < 0] ≥ (1 − 7ηLN8B )x0. Let σ m be a random permutation of N2 +1’s and N 2 −1’s.
Base case. x0,F ≥ x0,G since both start the epoch at the same point x0.
Inductive case. There can be three cases:
• Case 1: xi,F ≥ xi,G ≥ 0. Then,
xi+1,F − xi+1,G = xi,F − xi,G − η
MB M∑ m=1 (i+1)B∑ j=iB+1 (∇fσmj (xi,F )−∇gσmj (xi,G))
= xi,F − xi,G − η
MB M∑ m=1 (i+1)B∑ j=iB+1 ( µxi,F + νσ m j − Lxi,G − νσmj ) = xi,F − xi,G − η (µxi,F − Lxi,G) = xi,F (1− ηµ)− xi,G(1− ηL) ≥ 0.
• Case 2: 0 ≥ xi,F ≥ xi,G. Then,
xi+1,F − xi+1,G = xi,F − xi,G − η
MB M∑ m=1 (i+1)B∑ j=iB+1 (∇fσmj (xi,F )−∇gσmj (xi,G))
= xi,F − xi,G − η
MB M∑ m=1 (i+1)B∑ j=iB+1 ( Lxi,F + νσ m j − Lxi,G − νσmj ) = xi,F − xi,G − η (Lxi,F − Lxi,G) = xi,F (1− ηL)− xi,G(1− ηL) ≥ 0.
• Case 3: xi,F ≥ 0 ≥ xi,G. Then,
xi+1,F − xi+1,G = xi,F − xi,G − η
MB M∑ m=1 (i+1)B∑ j=iB+1 (∇fσmj (xi,F )−∇gσmj (xi,G))
= xi,F − xi,G − η
MB M∑ m=1 (i+1)B∑ j=iB+1 ( µxi,F + νσ m j − Lxi,G − νσmj ) = xi,F − xi,G − η (µxi,F − Lxi,G) = xi,F (1− ηµ)− xi,G(1− ηL).
Note that since η ≤ BLN , xi,F (1 − ηµ) ≥ 0 and xi,G(1 − ηL) ≤ 0, which proves that xi+1,F − xi+1,G ≥ 0.
Thus, we see that xi+1,F ≥ xi+1,G. Further, by linearity of expectation and gradient, it is easy to check that
E[xN B ,G ] = E[xN B−1,G − η∇G2(xN B−1,G )]
= (1− ηL)E[xN B−1,G ] = · · ·
= (1− ηL)NB x0.
Using the result that xN B ,F ≥ xN B ,G which we proved above, we get E[xN B ,F ] ≥ (1 − ηL)NB x0 for any initial iterate x0. Specifically for x0 < 0, this implies E[xN B ,F | x0 < 0] ≥ (1− ηL) N B x0.
Further, since ηL ≤ B513N , we have (1− ηL) N B ≤ 1− 7ηLN8B . This is because 1− 7zN 8B − (1− z) N B is nonnegative on the interval [ 0, 1− (7/8) 1 N/B−1 ] , and 1− (7/8)
1 N/B−1 ≥ B513N for all N B ≥ 2. To
see why, note that (1 − 1513(n−1) ) n−1 ≥ 78 for all n ≥ 2, and this gives 1 − (7/8) 1 n−1 ≥ 1513(n−1) , which then implies 1 − (7/8) 1
n−1 ≥ 1513n for all n ≥ 2. Therefore, for x0 < 0, we have E[xNB ,F | x0 < 0] ≥ (1− 7ηLN8B )x0. For the last statement of the lemma, note that by symmetry of the function G2, if we initialize the Algorithm 2 at 0, then for any starting iterate of an epoch we have P(x0,G ≥ 0) ≥ 1/2. This combined with the fact that xi,F ≥ xi,G gives us that P(x0,F ≥ 0) ≥ 1/2.

heading: F.4 PROOF OF LEMMA 14
text: For m = 1, . . . ,M , let σm be a random permutation of N2 +1’s and N 2 −1’s. Then, we first show that for any i ≤ N/2 and k ≤ B/2,
1
64
(√ i
M + √ k
) ≤ E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  .
To prove the lower bound, we will use Khintchine’s inequality along with Lemma 12 from Rajput et al. (2020). Let us define random variables am := | 1M ∑i j=1 σ m j |, xm := sign( ∑i j=1 σ m j ), bM :=
| ∑i+k j=i+1 σ M j |, and yM := sign( ∑i+k j=i+1 σ M j ). For xm, if the sum ∑i j=1 σ m j = 0 then xm is +1 with probability 0.5 and −1 with probability 0.5. Ties occurring in yM are also broken similarly. We can note that xm’s and yM are i.i.d. Rademacher random variables, which allows us to apply Khintchine’s inequality, Then, by Khintchine’s inequality,
E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  = E[∣∣∣∣∣ M∑ m=1 amxm + bMyM ∣∣∣∣∣ ]
≥ 1√ 2 E ( M∑ m=1 a2m + b 2 M )1/2 .

heading: By applying ‖z‖2 ≥ 1√d‖z‖1 for z ∈ R
text: d twice, we get
1√ 2 E ( M∑ m=1 a2m + b 2 M )1/2 ≥ 1 2 E ( M∑ m=1 a2m )1/2 + bM  ≥ 1 2 E [ 1√ M M∑ m=1 am + bM ] .
Next, noticing that am’s are i.i.d.,
1 2 E [ 1√ M M∑ m=1 am + bM ] = 1 2 E [√ MaM + bM ] = 1 2 E  1√ M ∣∣∣∣∣∣ i∑ j=1 σMj ∣∣∣∣∣∣+ ∣∣∣∣∣∣ i+k∑ j=i+1 σMj ∣∣∣∣∣∣ 
≥ 1 64
(√ i
M + √ k
) . (Lemma 12 from Rajput et al. (2020))
Note that Lemma 12 from Rajput et al. (2020) has the requirement that N ≥ 256. However, that requirement is for the entire lemma to hold, whereas we need only the first inequality in the lemma. For that, the requirement is simply N ≥ 8. Further, note that for N = 2, 4, and 6 it can be manually verified that the required inequalities in Lemma 12 of Rajput et al. (2020) hold. Hence, this lemma holds for all even N .
The upper bound comes from Jensen’s inequality:
E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  ≤ 1 M E ∣∣∣∣∣∣ M∑ m=1 i∑ j=1 σmj ∣∣∣∣∣∣ + E ∣∣∣∣∣∣ i+k∑ j=i+1 σMj ∣∣∣∣∣∣ 
≤ 1 M √√√√√√E   M∑ m=1 i∑ j=1 σmj 2 + √√√√√√E   i+k∑ j=i+1 σMj 2 
= 1
M √√√√√√ M∑ m=1 E   i∑ j=1 σmj 2 + √√√√√√E   i+k∑ j=i+1 σMj 2 
(Since ∑i j=1 σ m j for m = 1, 2 . . . are mean 0 and independent.)
= √√√√√√ 1M E   i∑ j=1 σMj 2 + √√√√√√E   i+k∑ j=i+1 σMj 2 
= √√√√√ 1 M i+∑ j 6=l E[σMj σMl ] +√k +∑ j 6=l E[σMj σMl ].
E[σMj σMl ] ≤ 0 because σMj and σMl are negatively correlated. Hence, we get
E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  ≤√ i M + √ k,
as desired.
Next, it is left to show that for 0 ≤ i ≤ N2 and 0 ≤ k ≤ B 2 satisfying i+ k ≥ 1, we have
P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj > 0  = P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj < 0  ≥ 1 6 .
By symmetry, proving the equality is straightforward, and hence it is sufficient prove that
P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj = 0  ≤ 2 3 . (149)
For this, it in fact suffices to show that
P  l∑ j=1 σ1j = 0  ≤ 2 3 for all 1 ≤ l ≤ N − 1, (150)
because (149) can be derived from (150). We first explain why (150) implies (149), and then show (150).
Suppose (150) is true. Then,
• Case 1: If i = 0, then (149) becomes P( ∑k j=1 σ M j = 0) ≤ 23 , which is true due to (150).
• Case 2: If M = 1, then (149) becomes P( ∑i+k j=1 σ 1 j = 0) ≤ 23 , which is true due to (150). • Case 3: If i ≥ 1 and M ≥ 2, then we can consider two events that partition the probability space:
1. E1 := { ∑M m=2 ∑i j=1 σ m j +M ∑i+k j=i+1 σ M j = 0}. Conditioned on this event E1,
P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj = 0 | E1  = P  i∑ j=1 σ1j = 0  ≤ 2 3
due to independence of machines and (150). 2. E2 := { ∑M m=2 ∑i j=1 σ m j +M ∑i+k j=i+1 σ M j 6= 0}. Conditioned on this event E2, let
c := ∑M m=2 ∑i j=1 σ m j +M ∑i+k j=i+1 σ M j . Then,
P  M∑ m=1 i∑ j=1 σmj +M i+1∑ j=i+1 σMj = 0 | E2  = P  i∑ j=1 σ1j = −c  .
However, by symmetry, P (∑i
j=1 σ 1 j = −c
) = P (∑i j=1 σ 1 j = c ) ≤ 12 .
From these two events, we conclude that (149) must hold.
It is now left to prove (150). It is clear that P( ∑i j=1 σ 1 j = 0) = 0 for all odd i, so we assume that
i is even. Also note that P( ∑i j=1 σ 1 j = 0) = P( ∑N j=i+1 σ 1 j = 0) = P( ∑N−i j=1 σ 1 j = 0), since σ 1 j ’s sum to zero. Therefore, for the rest of the proof, we can focus on even i’s in the range 2 ≤ i ≤ N2 . Note that P( ∑i j=1 σ m j = 0) is just the probability of having i 2 +1’s and i 2 −1’s in the first i spots in a random shuffling of N2 +1’s and N 2 −1’s. This is equivalent to choosing i 2 indices (for +1) out of the first i, and then choosing N−i2 indices out of the remaining N − i. Thus,
P  i∑ j=1 σmj = 0  = ( ii/2) · (N−iN−i2 )( N N/2
) . Note that the term above is a decreasing function of i for i ≤ N/2. Hence, putting i = 2 to the RHS we get
P  i∑ j=1 σmj = 0  ≤ (21) · (N−2N−22 )( N N/2 ) = N 2(N − 1) ≤ 2 3 ,
where the inequality holds for N ≥ 4.

heading: F.5 PROOF OF LEMMA 15
text: E[|xi − x0|] = E ∣∣∣∣∣∣ ηMB i−1∑ j=0 M∑ m=1 (j+1)B∑ k=jB+1 νσmk + (L1xj<0 + µ1xj≥0)xj ∣∣∣∣∣∣ 
≤ ην B
√ iB
M + ηE ∣∣∣∣∣∣ i−1∑ j=0 (L1xj<0 + µ1xj≥0)xj ∣∣∣∣∣∣  (By Lemma 14)
≤ ην √ i
MB + ηL i−1∑ j=0 E[|xj |]
≤ ην √ i
MB + iηLx0 + ηL i−1∑ j=0 E[|xj − x0|].
Let h(i) := ην √
i MB + iηLx0 + ηL ∑i−1 j=0 h(j), starting with h(0) = 0. Then using induction, it
can be seen that E[|xi − x0|] ≤ h(i). Further, since h(i) is an increasing function of i, we get h(i) ≤ ην √ i MB + iηLx0
1− iηL .
Since i ≤ NB and η ≤ B 513LN , we get that 1 1−iηL ≤ 513 512 , so E[|xi−x0|] ≤ 513 512ην
√ i
MB + 513 512 iηLx0.

heading: G PROOF OF LOWER BOUND FOR LOCAL RR: HOMOGENEOUS CASE (THEOREM 4)
text: Recall that Theorem 4 gives the bound for local RR in the homogeneous setting, where all machines have the same local objectives. Similar to Theorem 3, we consider three step-size ranges and do case analysis for each of them. We construct functions for each corresponding step-size regime such that the convergence of local RR is “slow” for the functions on their corresponding step-size regime. The final lower bound is the minimum among the lower bounds obtained for the three regimes. More
concretely, we will construct three one-dimensional functions F1(x), F2(x), and F3(x) satisfying L-smoothness (1), µ-PŁ condition (2), and Assumption 2 such that15
• Local RR on F1(x) with η ≤ 1µNK and initialization y0 = ν µ results in
E[F1(yK,NB )] = Ω ( ν2
µ
) .
• Local RR on F2(x) with η ≥ 1µNK and η ≤ 1 1025LN and initialization y0 = 0 results in
E[F2(yK,NB )] = Ω (
ν2
µMNK2 +
ν2B
µN2K2
) .
Note that the step-size range requires K ≥ 1025κ, hence this lower bound occurs only in the “large-epoch” regime, i.e., K & κ.
• Local RR on F3(x) with η ≥ 1µNK and η ≥ 1 1025LN and initialization y0 = 0 results in
E[F3(yK,NB )] = Ω ( ν2
µMNK
) .
Then, the three dimensional function F ([x, y, z]>) = F1(x) + F2(y) + F3(z) will show bad convergence in any step-size regime. Furthermore,
µI min(∇2F1,∇2F2,∇2F3)I ∇2F max(∇2F1,∇2F2,∇2F3)I LI,
that is, if F1, F2 and F3 are µ-strongly convex and L-smooth, then so is F . Moreover, since the component functions in each coordinate are designed to satisfy Assumption 2 with ν, the resulting three dimensional function F also satisfies Assumption 2 with √ 3ν.
Since the final lower bound is the minimum among the lower bounds obtained in the step-size ranges, the lower bound becomes Ω ( ν2
µMNK2 + ν2B µN2K2 ) if K ≥ 1025κ and K ≥ MBN (this inequality is
required to make sure ν 2B µN2K2 ≤ ν2 µMNK ) and Ω ( ν2 µMNK ) otherwise.
In the subsequent subsections, we prove the lower bounds for F1, F2, and F3 separately.
G.1 LOWER BOUND FOR η ≤ 1µNK
Consider the case where every function at every machine is the same: for all i ∈ [N ] and m ∈ [M ], fmi (x) := µx2 2 . Hence, F1(x) = µx2 2 .
Since all fmi ’s are the same, the local updates in all the machines are identical. Hence, for this subsection we omit the superscript for local machines. Let xk,0 and xk,N denote the local iterates at the beginning and end of the k-th epoch. Then,
xk,N = (1− ηµ)Nxk,0.
Initializing at x1,0 = νµ and repeating this for K epochs, we get that after K epochs, the last iterate yK,NB = xK,N satisfies
yK,NB = xK,N = (1− ηµ)NK ·
ν µ ≥ ( 1− 1 NK )NK · ν µ ≥ ν 4µ ,
since N ≥ 2, and K ≥ 1. Hence, F1(yK,NB ) = Ω( ν2 µ ).
15Again, the functions constructed in this theorem are µ-strongly convex, which is stronger than µ-PL required in Definition 1. Also, our functions satisfy Assumption 3 with τ = 0, ρ = 1.
G.2 LOWER BOUND FOR η ≥ 1µNK AND η ≤ 1 1025LN
For most part of this subsection, we consider iterates within a single epoch, and hence we will omit the subscripts denoting epochs. Let x0 denote the iterate at the beginning of the epoch (which is the same across all the machines), and xmi denote the i-th local iterate for machine m. After every B local iterates, the server aggregates the local iterates xmiB , computes their average yi := 1 M ∑M m=1 x m iB , and synchronizes all the machines x m iB := yi.
Let y0 denote the iterate at the beginning of the epoch (which is the same across all the machines xm0 = y0), and x m i denote the i-th local iterate at machine m. In our construction, each machine will have the same set of component functions, that is, there will be no inter-machine deviation. We therefore omit the superscriptm from the local component functions fmi . The function we construct for the lower bound and its component functions are as follows:
F2(x) := 1
N  N2∑ i=1 f+1(x) + N∑ i=N2 +1 f−1(x)  , where f+1(x) := (L1x≤0 + µ1x>0) x2
2 + νx, and
f−1(x) := (L1x≤0 + µ1x>0) x2
2 − νx
Note that the function F2(x) = (L1x≤0 + µ1x>0)x 2
2 is µ-strongly convex and L-smooth with minimizer at 0, and also satisfies Assumption 2.
Let σm be a random permutation of N2 +1’s and N 2 −1’s. Then, machine m computes gradients on f−1 and f+1 in the order given by σm. Let σmj denote the j-th ordered element of σ m. Then,
∇fσmj (x) = (L1x≤0 + µ1x>0)x+ νσ m j .
Hence, the last iterate of an epoch, yN B , is given by
yN B − y0 = N B−1∑ i=0 − η M M∑ m=1 B−1∑ j=0 ∇fσmiB+j+1(x m iB+j)  = − η
M N B−1∑ i=0 M∑ m=1 B−1∑ j=0 ((L1xmiB+j≤0 + µ1xmiB+j>0)x m iB+j + νσ m iB+j+1)
= − η M N B−1∑ i=0 M∑ m=1 B−1∑ j=0 (L1xmiB+j≤0 + µ1xmiB+j>0)x m iB+j ,
where, in the last line, we used the fact that ∑N j=1 σ m j = 0.
Recall that in the construction, each machine has the same component functions. Hence,
E[yN B − y0] = −
η
M E NB−1∑ i=0 M∑ m=1 B−1∑ j=0 (L1xmiB+j≤0 + µ1xmiB+j>0)x m iB+j  = −η
N B−1∑ i=0 B−1∑ j=0 E [ (L1x1iB+j≤0 + µ1x1iB+j>0)x 1 iB+j ] , (151)
where the last equality holds because the iterates xmiB+j are identically distributed across different m ∈ [M ]. Hence, we need to bound E[(L1x1iB+j≤0 +µ1x1iB+j>0)x 1 iB+j ]. As we did for Theorem 3, we want to prove that E[yN B ] keeps increasing over an epoch, that is E[yN B − y0] > 0 when y0 is close enough to the minimizer 0.
For this, we first consider the case where the first iterate y0 of the epoch satisfies y0 ≥ 0. The y0 < 0 case will be considered later. For the case y0 ≥ 0, we will show that whenever y0 is small, the expected amount of update made in the (iB+j+1)-th iteration, E[(L1x1iB+j≤0+µ1x1iB+j>0)x 1 iB+j ], is negative if i ≤ b N2B c and B 4 ≤ j ≤ B 2 , and not too big otherwise.
We use the following lemmas, proven in Appendices H.1 and H.2, respectively. Lemma 16. For y0 ≥ 0, 0 ≤ i ≤ b N2B c, B 4 ≤ j ≤ B 2 , η ≤ 1 1025LN , and L µ ≥ 15375 2 ,
E[(L1x1iB+j≤0 + µ1x1iB+j>0)x 1 iB+j ] ≤
6 7 Ly0 − ηLν 1536
(√ iB
M + √ j ) .
Lemma 17. For y0 ≥ 0, 0 ≤ i ≤ NB − 1, 0 ≤ j ≤ B − 1, and η ≤ 1 1025LN ,
E[(L1x1iB+j≤0 + µ1x1iB+j>0)x 1 iB+j ] ≤ µ
( 1+ 1025(iB + j)ηL
1024
) y0+ 1025ηµν
1024
(√ iB
M + √ j ) .
Next, we apply the two lemmas above in (151). For now, we consider the case N/B ≥ 2. The case B = N will be handled separately at the end of this subsection. For simplicity of notation, define EiB+j := E[(L1x1iB+j≤0 + µ1x1iB+j>0)x 1 iB+j ]. We will divide the summation in (151) into four groups; for one of them we can apply Lemma 16, and for the other three we apply Lemma 17.
E[yN B − y0] = −η N B−1∑ i=0 B−1∑ j=0 EiB+j
= − η b N2B c∑ i=0 B 4 −1∑ j=0 EiB+j + b N2B c∑ i=0 B 2∑
j=B4
EiB+j + b N2B c∑ i=0 B−1∑ j=B2 +1 EiB+j + N B−1∑ i=b N2B c+1 B−1∑ j=0 EiB+j  ≥ − η
b N2B c∑ i=0 B 4 −1∑ j=0
( µ ( 1 + 1025(iB + j)ηL
1024
) y0 + 1025ηµν
1024
(√ iB
M + √ j
))
− η b N2B c∑ i=0 B 2∑
j=B4
( 6
7 Ly0 −
ηLν
1536
(√ iB
M + √ j
))
− η b N2B c∑ i=0 B−1∑ j=B2 +1
( µ ( 1 + 1025(iB + j)ηL
1024
) y0 + 1025ηµν
1024
(√ iB
M + √ j
))
− η N B−1∑
i=b N2B c+1
B−1∑ j=0
( µ ( 1 + 1025(iB + j)ηL
1024
) y0 + 1025ηµν
1024
(√ iB
M + √ j
))
≥ − η b N2B c∑ i=0 B 2∑
j=B4
( 6
7 Ly0 −
ηLν
1536
(√ iB
M + √ j
))
− η N B−1∑ i=0 B−1∑ j=0
( µ ( 1 + 1025(iB + j)ηL
1024
) y0 + 1025ηµν
1024
(√ iB
M + √ j )) , (152)
where the last inequality is true because the RHS of the inequality in Lemma 17 is nonnegative. First consider the terms in (152) that involve y0. Since (iB + j)ηL ≤ ηLN ≤ 11025 , µ ≤ 2L 15375 , N/B ≥ 2, and B ≥ 4, we have the following loose bound:
b N2B c∑ i=0 B 2∑
j=B4
6 7 L+ N B−1∑ i=0 B−1∑ j=0 µ ( 1 + 1025(iB + j)ηL 1024 )
≤ (⌊ N
2B
⌋ + 1 )( B
4 + 1
) 6
7 L+
N B ·B · 2L 15375
( 1 + 1
1024 ) ≤ N B · B 2 · 6 7 L+ LN 7680 ≤ 7LN 8 . (153)
We next bound the terms in (152) that involve summation of square roots. From N/B ≥ 2, we have b N2B c ≥ N 3B and b N 2B c+ 1 ≥ N 2B , so
b N2B c∑ i=0 B 2∑
j=B4
(√ iB
M + √ j ) = ( B 4 + 1 )√ B M b N2B c∑ i=0 √ i+ (⌊ N 2B ⌋ + 1 ) B2∑ j=B4 √ j
≥ B 3/2
4M1/2 ∫ b N2B c 0 √ tdt+ N 2B ∫ B 2
B 4 −1
√ tdt
≥ B 3/2
6M1/2
( N
3B
)3/2 + N
3B
[( B
2
)3/2 − ( B
4
)3/2]
= N3/2
18 √ 3M1/2 +
(2 √ 2− 1)NB1/2
24 . (154)
For the other sum, we have N B−1∑ i=0 B−1∑ j=0 (√ iB M + √ j ) = B3/2 M1/2 N B−1∑ i=0 √ i+ N B B−1∑ j=0 √ j
≤ B 3/2
M1/2
∫ N B
0
√ tdt+ N
B ∫ B 0 √ tdt
= 2B3/2
3M1/2
( N
B
)3/2 + 2N
3B B3/2
= 2N3/2
3M1/2 +
2NB1/2
3 . (155)
Substituting the bounds (153), (154), and (155) into (152), and using µ ≤ L40000 ,
E[yN B − y0] ≥ −
7ηLN
8 y0 +
η2Lν
1536
( N3/2
18 √ 3M1/2 +
(2 √ 2− 1)NB1/2
24
)
− 1025η 2µν
1024
( 2N3/2
3M1/2 +
2NB1/2
3 ) ≥ −7ηLN
8 y0 +
η2Lν
240000
( N3/2
M1/2 +NB1/2
) . (156)
For the other case y0 < 0, we have the following lemma, a local RR counterpart of Lemma 13. In Appendix H.3, we prove the following: Lemma 18. If η ≤ 11025LN and an epoch starts at y0 < 0, then
E [ yN B | y0 < 0 ] ≥ ( 1− 7ηLN 8 ) y0.
Further, if the first epoch of the algorithm is initialized at 0, then for any starting iterate y0 of any following epoch, we have P(y0 ≥ 0) ≥ 1/2.
Using (156) and Lemma 18, we get
E[yN B ]
=P(y0 ≥ 0)E[yN B | y0 ≥ 0] + P(y0 < 0)E[yN B | y0 < 0]
≥P(y0 ≥ 0) ((
1− 7ηLN 8
) y0 + η2Lν
240000
( N3/2
M1/2 +NB1/2
)) + P(y0 < 0) ( 1− 7ηLN
8
) y0
≥ (
1− 7ηLN 8
) y0 + η2Lν
480000
( N3/2
M1/2 +NB1/2
) .
Thus far, we have characterized the expected per-epoch update, starting from the initial iterate y0 and iterating until the last iterate yN
B of the epoch. Now recall that we run the algorithm for K
epochs. Using yk,i to denote the i-th aggregated iterate of the k-th epoch, we get a lower bound on the expectation of the last iterate yk,NB if we initialize at y1,0 = 0:
E[yK,NB ] ≥ ( 1− 7ηLN 8 )K y1,0 + η2Lν 480000 ( N3/2 M1/2 +NB1/2 )K−1∑ k=0 ( 1− 7ηLN 8 )k
= η2Lν
480000
( N3/2
M1/2 +NB1/2
) 1− (1− 7ηLN8 )K 7ηLN
8
= ην
420000
( N1/2
M1/2 +B1/2
)( 1− ( 1− 7ηLN
8
)K)
≥ ην 420000
( N1/2
M1/2 +B1/2
)( 1− ( 1− 7L
8µK
)K) . (Since η ≥ 1µNK )
Note that since Lµ ≥ 40000 and K ≥ 1025L µ (which is implied by 1 µNK ≤ η ≤ 1 1025LN ),
1− (
1− 7L 8µK
)K ≥ 1− e− 7L 8µ ≥ 1− e−35000 ≈ 1.
Hence, we get from η ≥ 1µNK that E[yK,NB ] = Ω ( ηνN1/2 M1/2 + ηνB1/2 ) = Ω ( ν µM1/2N1/2K + νB1/2 µNK ) ,
and by Jensen’s inequality, we finally have
E[F (yK,NB )] ≥ 1 2 E[µy2 K,NB ] = Ω(µE[yK,NB ] 2) = Ω
( ν2
µMNK2 +
ν2B
µN2K2
) .
Recall that, from the paragraph below Lemmas 16 and 17 to this point, we have assumed N/B ≥ 2. We handle the case B = N now. In this case, notice that all the √ iB M terms that appear in (152) disappear, because we always have i = 0. Therefore, the proof goes through in the same why,
modulo the fact that we do not have the terms that originate from the √
iB M terms in (152). Therefore,
we can show
E[F (yK,NB )] ≥ 1 2 E[µy2 K,NB ] = Ω(µE[yK,NB ] 2) = Ω
( ν2
µNK2
) .
G.3 LOWER BOUND FOR η ≥ 1µNK AND η ≥ 1 1025LN
Similar to earlier parts of the proof, here as well, each machine will have the same component functions, that is, there will be no inter-machine deviation. The proof uses a similar construction as Safran & Shamir (2020; 2021):
F3(x) := 1
N  N2∑ i=1 f+1(x) + N∑ i=N2 +1 f−1(x)  , where f+1(x) := Lx2
2 + νx, and f−1(x) :=
Lx2
2 − νx.
Hence, F3(x) = Lx 2
2 , and has its minimizer at 0.
We first compute the expected “progress” over a given epoch. For simplicity, let us omit the subscript for epochs for now. Let y0 denote the iterate at the beginning of the epoch (which is the same across all the machines xm0 = y0) and x m i denote the i-th local iterate for machine m. After every B local iterates, the server aggregates the local iterates xmiB , computes their average yi := 1 M ∑M m=1 x m iB , and synchronizes all the machines xmiB := yi.
For the epoch, let σm be the permutation of N2 +1’s and N 2 −1’s sampled by machine m. Upon receiving the aggregated iterate xmiB = yi, each machine m performs B local updates. Unrolling the local update rules, the iterate after the B updates (and before synchronization) can be written as follows:
xm(i+1)B = x m (i+1)B−1 − η∇fσm(i+1)B (x m (i+1)B−1)
= xm(i+1)B−1 − η(Lx m (i+1)B−1 + νσ m (i+1)B)
= (1− ηL)xm(i+1)B−1 − ηνσ m (i+1)B
= · · · = (1− ηL)BxmiB − ην (i+1)B∑ j=iB+1 (1− ηL)(i+1)B−jσmj = (1− ηL)Byi − ην (i+1)B∑ j=iB+1 (1− ηL)(i+1)B−jσmj ,
After synchronization, we get
yi+1 := 1
M M∑ m=1 xm(i+1)B = (1− ηL) Byi − ην M M∑ m=1 (i+1)B∑ j=iB+1 (1− ηL)(i+1)B−jσmj .
Unrolling the equation above from yN B (the final iterate of the epoch, after synchronization) to y0 (the starting iterate), we get
yN B = (1− ηL)ByN B−1 − ην M M∑ m=1 N∑ j=N−B+1 (1− ηL)N−jσmj
= · · · = (1− ηL)Ny0 − N B∑ i=1 (1− ηL)N−iB ην M M∑ m=1 iB∑ j=(i−1)B+1 (1− ηL)iB−jσmj
= (1− ηL)Ny0 − ην
M M∑ m=1 N∑ j=1 (1− ηL)N−jσmj .
Then, by squaring both sides and taking expectations,
E[y2N B
] = (1− ηL)2Ny20 − 2ην(1− ηL)Nx0
M E [ M∑ m=1 N∑ i=1 (1− ηL)N−iσmi ]
+ η2ν2
M2 E ( M∑ m=1 N∑ i=1 (1− ηL)N−iσmi )2 = (1− ηL)2Ny20 + η2ν2
M2 E ( M∑ m=1 N∑ i=1 (1− ηL)N−iσmi )2 , (157)
where we used the fact that E[σmi ] = 0. Further, because σm and σm ′
are independent‘ and identically distributed for different m and m′, we get that
E (ην M M∑ m=1 N∑ i=1 (1− ηL)N−iσmi )2 =
M∑ m=1 E ( N∑ i=1 (1− ηL)N−iσmi )2+ ∑ m 6=m′ E [ N∑ i=1 (1− ηL)N−iσmi ] E [ N∑ i=1 (1− ηL)N−iσm ′ i ]
=ME ( N∑ i=1 (1− ηL)N−iσ1i )2 , where the last equality used the fact that E[σmj ] = 0 for all m ∈ [M ] and i ∈ [N ], and that σm are identically distributed. Since we only consider the permutation σ1 (i.e., the one for machine 1) from now on, we henceforth omit the superscript. Substituting this to (157) gives
E[y2N B
] = (1− ηL)2Ny20 + η2ν2
M E ( N∑ i=1 (1− ηL)N−iσi )2 ︸ ︷︷ ︸
:=Φ
. (158)
From (158), we have calculated the per-epoch expected update, because the final iterate yN B is also the initial iterate of the next epoch. Recall that we run the algorithm for K epochs. We now use yk,i to denote the iterate after the i-th communication round in the k-th epoch. Using (158), we get a lower bound on the expectation of the last iterate yK,NB squared:
E[y2 K,NB
] = (1− ηL)2NKy21,0 + η2ν2
M Φ K−1∑ k=0 (1− ηL)2Nk ≥ η 2ν2 M Φ. (159)
where the inequality used that we initialize at y1,0 = 0 and ∑K−1 k=0 (1− ηL)2Nk ≥ (1− ηL)0 = 1. Next, we bound the expectation term, i.e., Φ, defined in (158). Using Lemma 1 from Safran & Shamir (2020) with n and α replaced with N and ηL respectively, we have
Φ ≥ c ·min { 1
ηL , η2L2N3
} , (160)
for some universal constant c > 0. Using the fact that η ≥ 11025LN , it is easy to check that the RHS of (160) is lower-bounded by c ′
ηL , where c ′ > 0 is a universal constant. Combining (159) and (160)
gives
E[y2 K,NB
] ≥ η 2ν2 M · c ′ ηL = c′ην2 LM .
Also using the fact that η ≥ 1µNK , we get
E[F3(yK,NB )] = L 2 E[y2 K,NB ] ≥ c
′ν2
2µMNK .

heading: H PROOFS OF HELPER LEMMAS FOR APPENDIX G

heading: H.1 PROOF OF LEMMA 16
text: The proof of Lemma 16 is similar to its minibatch RR counterpart, Lemma 11. From the given 0 ≤ i ≤ NB − 1, 0 ≤ j ≤ B − 1, define k := iB + j, in order to simplify notation. We also define E := 1M (∑M m=1 ∑iB l=1 σ m l ) + ∑k l=iB+1 σ 1 l .
By the law of total expectation we have
E[(L1x1k≤0 + µ1x1k>0)x 1 k] = P (E > 0)E [ (L1x1k≤0 + µ1x1k>0)x 1 k ∣∣∣E > 0]
+ P (E ≤ 0)E [ (L1x1k≤0 + µ1x1k>0)x 1 k ∣∣∣E ≤ 0] ≤ P (E > 0)LE [ x1k ∣∣E > 0]+ P (E ≤ 0)µE [x1k∣∣E ≤ 0] , (161)
where the last inequality used the fact that (L1t≤0 + µ1t>0)t ≤ Lt and (L1t≤0 + µ1t>0)t ≤ µt for any t ∈ R. We handle each of the two expectations in (161) separately. We first bound E [ x1k ∣∣E > 0]. Recall from the definition of algorithm iterates that
x1k − y0 = x1iB+j − y0 = − η
M M∑ m=1 iB−1∑ l=0 ∇fσml+1(x m l )− η k−1∑ l=iB fσ1l+1(x 1 l ). (162)
Expanding (162) using the definition of∇fσml+1 ’s, we obtain E [ x1k ∣∣E > 0]
=E [ y0 − η
M M∑ m=1 iB−1∑ l=0 ( νσml+1 + (L1xml ≤0 + µ1xml >0)x m l ) −η
k−1∑ l=iB ( νσ1l+1 + (L1x1l≤0 + µ1x1l>0)x 1 l )∣∣∣∣∣E > 0 ]
=E [ y0 − η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)(x m l − y0)
−η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0)(x 1 l − y0) ∣∣∣∣∣E > 0 ] − ηνE [E|E > 0]
− E
[ η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)y0 + η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0)y0 ∣∣∣∣∣E > 0 ]
= y0E [ 1− η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)− η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0) ∣∣∣∣∣E > 0 ]
− E
[ η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)(x m l − y0)
+η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0)(x 1 l − y0) ∣∣∣∣∣E > 0 ] − ηνE [E|E > 0]
≤ y0E [ 1− η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)− η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0) ∣∣∣∣∣E > 0 ]
+ ηL(M − 1)
M
iB−1∑ l=0 E [ |x2l − y0| | E > 0 ] + ηL M iB−1∑ l=0 E [ |x1l − y0| | E > 0 ] + ηL
k−1∑ l=iB E [ |x1l − y0| | E > 0 ] − ηνE [E|E > 0] ,
where the last inequality used the fact that (L1t≤0 + µ1t>0)t ≤ Lt and (L1t≤0 + µ1t>0)t ≤ µt for any t ∈ R; and that for different m = 2, . . . ,M , the local iterates xml are identically distributed conditioned on E > 0. Next, we use the fact that for any nonnegative random variable v and event ∆, we have that E[v|∆] ≤ E[v]/P(∆). Hence,
E [ x1k ∣∣E > 0]
≤ y0E [ 1− η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)− η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0) ∣∣∣∣∣E > 0 ]
+ ηL(M − 1)
M
iB−1∑ l=0
E [ |x2l − y0| ] P(E > 0) + ηL M iB−1∑ l=0 E [ |x1l − y0| ] P(E > 0) + ηL k−1∑ l=iB E [ |x1l − y0| ] P(E > 0)
− ηνE [E|E > 0]
= y0E [ 1− η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)− η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0) ∣∣∣∣∣E > 0 ]
− ηνE [E|E > 0] + ηL k−1∑ l=0
E [ |x1l − y0| ] P(E > 0)
(163)
where the last equality used the fact that for different m ∈ [M ], the local iterates xml are identically distributed when they are not conditioned.
Next, we use Lemma 14 again to bound the conditional expectations that arise in (163). We restate the lemma for the reader’s convenience. Lemma 14. For m ∈ [M ], let σm be a random permutation of N2 +1’s and N 2 −1’s. Then, for any i ≤ N2 and k ≤ B 2 , we have
1
64
(√ i
M + √ k
) ≤ E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  .
Furthermore, for any 0 ≤ i ≤ N and 0 ≤ k ≤ N satisfying i+ k ≤ N , we have
E ∣∣∣∣∣∣  1 M M∑ m=1 i∑ j=1 σmj + i+k∑ j=i+1 σMj ∣∣∣∣∣∣  ≤√ i M + √ k.
Lastly, for any 0 ≤ i ≤ N2 and 0 ≤ k ≤ B 2 satisfying i+ k ≥ 1, we have
P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj > 0  = P  M∑ m=1 i∑ j=1 σmj +M i+k∑ j=i+1 σMj < 0  ≥ 1 6 .
Lemma 14 implies that E [E|E > 0] ∈ [
1 64 (√ iB M + √ j ) , √ iB M + √ j ]
and P(E > 0) = P(E < 0) ≥ 1/6. From this, we get
ηνE [E|E > 0] ≥ ην 64
(√ iB
M + √ j ) , (164)
E [ |x1l − y0| ] P(E > 0) ≤ 6E [ |x1l − y0| ] . (165)
Also, since η ≤ 1LN we have
1− kηµ ≥ 1− η M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)− η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0) ≥ 1− ηLN ≥ 0,
which implies that
y0E [ 1− η
M M∑ m=1 iB−1∑ l=0 (L1xml ≤0 + µ1xml >0)− η k−1∑ l=iB (L1x1l≤0 + µ1x1l>0) ∣∣∣∣∣E > 0 ] ≤(1− kηµ)y0.
(166)
Substituting (164), (165), and (166) to (163), we obatin
E [ x1k ∣∣E > 0] ≤ (1− kηµ)y0 − ην
64
(√ iB
M + √ j ) + 6ηL k−1∑ l=0 E [ |x1l − y0| ] . (167)
Next, we have the following lemma that we can apply to E[|x1l − y0|]. Proof of Lemma 19 can be found in Appendix H.4.
Lemma 19. For y0 ≥ 0, 0 ≤ i ≤ NB − 1, 0 ≤ j ≤ B − 1, and η ≤ 1 1025LN ,
E[|x1iB+j − y0|] ≤ 1025
1024 ην
(√ iB
M + √ j ) + 1025 1024 (iB + j)ηLy0.
Applying this lemma to (167) and arranging the bounds (recall that k := iB + j), we get
1024
1025 k−1∑ l=0 E [ |x1l − y0| ] ≤ ην
k−1∑ l=0 (√ Bbl/Bc M + √ l −Bbl/Bc ) + ηLy0 k−1∑ l=0 l
= ην iB−1∑ l=0 (√ Bbl/Bc M + √ l −Bbl/Bc )
+ ην iB+j−1∑ l=iB (√ Bbl/Bc M + √ l −Bbl/Bc ) + ηLy0 k−1∑ l=0 l
= ην ( B
i−1∑ l=0
√ lB
M + j
√ iB
M + i B−1∑ l=0 √ l + j−1∑ l=0 √ l
) + ηLy0
k−1∑ l=0 l. (168)
The terms in (168) can be bounded using ∑c−1 l=0 √ l ≤ ∫ c 0 √ tdt:
1024
1025 k−1∑ l=0 E [ |x1l − y0| ] ≤ ην ( 2i3/2B3/2 3M1/2 + i1/2jB1/2 M1/2 + 2iB3/2 3 + 2j3/2 3 ) + k2ηLy0 2
≤ ην
( (2iB + 3j)
3
√ iB
M +
4iBj1/2
3 +
2j3/2
3
) + k2ηLy0
2
= ην
( (2iB + 3j)
3
√ iB
M +
(4iB + 2j)
3
√ j ) + k2ηLy0
2
≤ 4kην 3
(√ iB
M + √ j ) + k2ηLy0 2 , (169)
where the second last inequality used B4 ≤ j (and hence B 1/2 ≤ 2j1/2), and the last inequality used 2iB + 3j ≤ 4(iB + j) = 4k and 4iB + 2j ≤ 4(iB + j) = 4k. Substituting (169) to (167), we get
E [ x1k ∣∣E > 0] (170)
≤ (1− kηµ)y0 − ην
64
(√ iB
M + √ j ) + 1025kη2Lν 128 (√ iB M + √ j ) + 3075k2η2L2y0 1024
= ( 1− kηµ+ 3075k 2η2L2
1024
) y0 − ( 1
64 − 1025kηL 128
) ην (√ iB M + √ j )
≤ (
1− kηµ+ 3kηL 1024
) y0 − ην
128
(√ iB
M + √ j ) . (171)
The last inequality here used kηL ≤ ηLN ≤ 11025 , which follows from η ≤ 1 1025LN . Thus far, we have obtained an upper bound for E [ x1k ∣∣E > 0]. Recall that there is another conditional expectation in (161) that we want to bound, namely E [ x1k ∣∣E ≤ 0]. We bound it below, using the tools developed so far. For i ≤ N2B and B4 ≤ j ≤ B2 ,
E [ x1k ∣∣E ≤ 0] = y0 + E [x1k − y0 | E ≤ 0]
≤ y0 + E [ |x1k − y0| | E ≤ 0 ] ≤ y0 + E [ |x1k − x0| ] P(E ≤ 0)
≤ y0 + 6E [ |x1k − x0| ] (Using Lemma 14)
≤ y0 + 3075ην
512
(√ iB
M + √ j ) + 3075kηLy0 512
(Using Lemma 19)
≤ ( 1 + 3075kηL
512
) y0 + 3075ην
512
(√ iB
M + √ j ) . (172)
Using (171) and (172) in (161), we get that for i ≤ N2B and B 4 ≤ j ≤ B 2 :
E[(L1x1k≤0 + µ1x1k>0)x 1 k] ≤P (E > 0)LE [ x1k ∣∣E > 0]+ P (E ≤ 0)µE [x1k∣∣E ≤ 0]
≤P (E > 0)L (( 1− kηµ+ 3kηL
1024
) y0 − ην
128
(√ iB
M + √ j
))
+ P (E ≤ 0)µ (( 1 + 3075kηL
512
) y0 + 3075ην
512
(√ iB
M + √ j )) . (173)
From Lemma 14, note that 16 ≤ P (E > 0) ≤ 5 6 and 1 6 ≤ P (E ≤ 0) ≤ 5 6 . We use these inequalities, along with kηL ≤ ηLN ≤ 11025 and L µ ≥ 15375 2 , to bound the terms appearing in (173).
P (E > 0)L (
1− kηµ+ 3kηL 1024
) y0 + P (E ≤ 0)µ ( 1 + 3075kηL
512
) y0
≤ 5 6 L
( 1 +
3
1049600
) y0 + 5
6 · 2L 15375
( 1 + 3
512
) y0 ≤ 6
7 Ly0. (174)
We also have
− P (E > 0) ηLν 128
(√ iB
M + √ j ) + P (E ≤ 0) 3075ηµν 512 (√ iB M + √ j )
≤ − ηLν 768
(√ iB
M + √ j ) + 5125ηµν 1024 (√ iB M + √ j )
≤ − ηLν 1536
(√ iB
M + √ j ) , (175)
where we used the assumption Lµ ≥ 15375 2 . Substituting (174) and (175) to (173), we get
E[(L1x1k≤0 + µ1x1k>0)x 1 k] ≤
6 7 Ly0 − ηLν 1536
(√ iB
M + √ j ) ,
which finishes the proof.

heading: H.2 PROOF OF LEMMA 17
text: E[(L1x1iB+j≤0 + µ1x1iB+j>0)x 1 iB+j ] ≤ µE[x1iB+j ]
= µy0 + µE[x1iB+j − y0] ≤ µy0 + µE[|x1iB+j − x0|]
≤ µy0 + 1025(iB + j)ηLµ
1024 y0 +
1025ηµν
1024
(√ iB
M + √ j ) ,
where the last inequality used Lemma 19.

heading: H.3 PROOF OF LEMMA 18
text: We consider iterates within a single epoch, and hence we omit the subscripts denoting epochs. In our construction, each machine has the same set of component functions, that is, there will be no inter-machine deviation. We therefore omit the superscript m from the local component functions. Consider the function
G2(x) := 1
N  N2∑ i=1 g+1(x) + N∑ i=N2 +1 g−1(x)  , where g+1(x) := Lx2
2 + νx, and g−1(x) :=
Lx2
2 − νx.
Hence, G2(x) = Lx 2
2 . We prove the lemma by coupling iterates corresponding to F2 and G2. In particular, we perform local RR on F2 and G2 such that both start the given epoch at y0 and all the corresponding machines use the same random permutations. Let xmiB+j,F and x m iB+j,G denote the iterates (for (iB + j)-th iteration at machine m) for F2 and G2 respectively. We use mathematical induction to prove that xmiB+j,F ≥ xmiB+j,G for all i = 0, . . . , NB −1 and j = 0, . . . , B and machines m = 1, . . . ,M . After that, we will use this to prove our desired statement E[yN
B ,F | y0 < 0] = E[ 1M ∑M m=1 x m N,F | y0 < 0] ≥ (1− 7ηLN 8B )y0.
Let σm be a random permutation of N2 +1’s and N 2 −1’s. First we consider i = 0 and 0 ≤ j ≤ B.
Base case. For the base case, we know that xm0,F ≥ xm0,G, since xm0,F = xm0,G = y0 for all m.
Inductive case. There can be three cases:
• Case 1: xmiB+j,F ≥ xmiB+j,G ≥ 0. Then,
xmiB+j+1,F − xmiB+j+1,G =xmiB+j,F − xmiB+j,G − η(∇fσmiB+j+1(xiB+j,F )−∇gσmiB+j+1(xiB+j,G))
=xmiB+j,F − xmiB+j,G − η ( µxmiB+j,F + νσ m iB+j+1 − LxmiB+j,G − νσmiB+j+1 ) =xmiB+j,F − xmiB+j,G − η ( µxmiB+j,F − LxmiB+j,G
) =xmiB+j,F (1− ηµ)− xmiB+j,G(1− ηL) ≥ 0.
• Case 2: 0 ≥ xmiB+j,F ≥ xmiB+j,G. Then,
xmiB+j+1,F − xmiB+j+1,G =xmiB+j,F − xmiB+j,G − η(∇fσmiB+j+1(xiB+j,F )−∇gσmiB+j+1(xiB+j,G))
=xmiB+j,F − xmi,G − η ( LxmiB+j,F + νσ m iB+j+1 − LxmiB+j,G − νσmiB+j+1 ) =xmiB+j,F − xmiB+j,G − η ( LxmiB+j,F − LxmiB+j,G
) =xmiB+j,F (1− ηL)− xmiB+j,G(1− ηL) ≥ 0.
• Case 3: xiB+j,F ≥ 0 ≥ xiB+j,G. Then,
xmiB+j+1,F − xmiB+j+1,G =xmiB+j,F − xmiB+j,G − η(∇fσmiB+j+1(xiB+j,F )−∇gσmiB+j+1(xiB+j,G))
=xmiB+j,F − xmiB+j,G − η ( µxmiB+j,F + νσ m iB+j+1 − LxmiB+j,G − νσmiB+j+1 ) =xmiB+j,F − xmiB+j,G − η ( µxmiB+j,F − LxmiB+j,G
) =xmiB+j,F (1− ηµ)− xmiB+j,G(1− ηL).
Note that since η ≤ 1LN , we get that x m iB+j,F (1 − ηµ) ≥ 0 and xmiB+j,G(1 − ηL) ≤ 0, which proves that xmiB+j+1,F − xmiB+j+1,G ≥ 0.
Thus, we see that xmiB+j+1,F ≥ xmiB+j+1,G for all the three cases, which proves by mathematical induction that xmiB+j,F ≥ xmiB+j,G for all 0 ≤ j ≤ B and i = 0. Note that this implies that, the aggregated averages y1,F := 1M ∑M m=1 x m B,F and y1,G := 1 M ∑M m=1 x m B,G satisfy y1,F ≥ y1,G. Hence, after synchronization is complete, we get that for i = 1 and j = 0, xmiB+j,F ≥ xmiB+j,G for all machines m. This proves the base case for i = 1. Now, we can repeat the Inductive cases for 1 ≤ j ≤ B and i = 1, and thereby prove that y2,F ≥ y2,G. Continuing on this process, we get that xmiB+j,F ≥ xmiB+j,G for all 0 ≤ j ≤ B and 0 ≤ i ≤ NB − 1, and consequently, yNB ,F ≥ yNB ,G. Further, by linearity of expectation and gradient, it is easy to check that for any machine m,
E[yN B ,G ] = (1− ηL)Ny0.
Using the result that yN B ,F ≥ yN B ,G which we proved above, we get E[yN B ,F ] ≥ (1 − ηL)Ny0 for any initial iterate y0. Specifically for y0 < 0, this implies E[yN B ,F | y0 < 0] ≥ (1− ηL)Ny0.
Further, since ηL ≤ 11025N , we have (1− ηL) N ≤ 1− 7ηLN8 . This is because 1− 7zN 8 − (1− z) N is nonnegative on the interval [ 0, 1− (7/8) 1 N−1 ] , and 1 − (7/8)
1 N−1 ≥ 11025N for all N ≥ 2. To
see why, note that (1− 11025(n−1) ) n−1 ≥ 78 for all n ≥ 2, and this gives 1− (7/8) 1 n−1 ≥ 11025(n−1) , which then implies 1− (7/8) 1
n−1 ≥ 11025n for all n ≥ 2. Therefore, for y0 < 0, we have E[yNB ,F | y0 < 0] ≥ (1− 7ηLN8 )y0. For the last statement of the lemma, note that by symmetry of the function G2, if we initialize Algorithm 1 at 0, then for any starting iterate of an epoch we have P(y0,G ≥ 0) ≥ 1/2. This combined with the fact that yi,F ≥ yi,G gives us that P(y0,F ≥ 0) ≥ 1/2.

heading: H.4 PROOF OF LEMMA 19
text: E[|x1iB+j − y0|]
=E [∣∣∣∣∣ ηM M∑ m=1 iB−1∑ l=0 νσml+1 + (L1xml <0 + µ1xml ≥0)x m l + η iB+j−1∑ l=iB νσ1l+1 + (L1x1l<0 + µ1x1l≥0)x 1 l ∣∣∣∣∣ ]
≤ ην
(√ iB
M + √ j ) + ηE [∣∣∣∣∣ iB+j−1∑ l=0 (L1x1l<0 + µ1x1l≥0)x 1 l ∣∣∣∣∣ ]
(By Lemma 14)
≤ ην
(√ iB
M + √ j ) + ηL iB+j−1∑ l=0 E[|x1l |]
≤ ην
(√ iB
M + √ j ) + (iB + j)ηLy0 + ηL iB+j−1∑ l=0 E[|x1l − y0|].
Now define
h(k) := ην (√ Bbk/Bc M + √ k −Bbk/Bc ) + kηLy0 + ηL k−1∑ l=0 h(l).
In terms of i and j, note that k corresponds to k = iB + j. Then using induction, it can be seen that E[|x1k − y0|] ≤ h(k). Further, since h(k) is an increasing function of k, we get
h(k) = ην (√ Bbk/Bc M + √ k −Bbk/Bc ) + kηLy0 + ηL k−1∑ l=0 h(l)
≤ ην (√ Bbk/Bc M + √ k −Bbk/Bc ) + kηLy0 + kηLh(k)
=⇒ h(k) ≤ ην
(√ Bbk/Bc M + √ k −Bbk/Bc ) + kηLy0
1− kηL .
Since k ≤ N and η ≤ 11025LN , we get that
E[|x1iB+j − y0|] ≤ 1025
1024 ην
(√ iB
M + √ j ) + 1025 1024 (iB + j)ηLy0,
as desired.

heading: I PROOF OF LOWER BOUND FOR LOCAL RR: HETEROGENEOUS CASE (PROPOSITION 5)
text: Recall that Proposition 5 gives the bound for local RR in the heterogeneous setting, where different machines have different local objectives. In this section, we construct examples where there is no intra-machine variation (i.e., fm1 = f m 2 = · · · = fmN for all m ∈ [M ]), but there is certain level of heterogeneity among different machines.
Similar to the other two lower bounds, we consider four step-size ranges and do case analysis for each of them. This time, we construct a single function F for these step-size regimes such that the convergence of local RR is “slow” for F . The final lower bound is the minimum among the lower bounds obtained for the four regimes. More concretely, we will construct a one-dimensional function F (x) satisfying L-smoothness (1), µ-PŁ condition (2), and Assumption 3 such that16
• Local RR on F (x) with η ≤ 18µNK and initialization y0 = τ µ results in
E[F (yK,NB )] = Ω ( τ2
µ
) .
• Local RR on F (x) with 18µNK ≤ η ≤ 1 8µB and initialization y0 = 0 results in
E[F (yK,NB )] = Ω ( τ2B2
µN2K2
) .
• Local RR on F (x) with 18µB ≤ η ≤ 1 µ and initialization y0 = 0 results in
E[F (yK,NB )] = Ω ( τ2
µ
) .
• Local RR on F (x) with η ≥ 1µ and initialization y0 = τ µ results in
E[F (yK,NB )] = Ω ( τ2
µ
) .
In the subsequent subsections, we prove the lower bounds for F for the four step-size intervals.
I.1 LOWER BOUND FOR 18µNK ≤ η ≤ 1 8µB AND 1 8µB ≤ η ≤ 1 µ
We first consider the two intervals in the middle, because they are more interesting cases. The global objective function F and its local objective functions are as follows.
F (x) := 1
M  M2∑ i=1 f1(x) + M∑ i=M2 +1 f2(x)  , where f1(x) := −τx, and f2(x) := µx2 + τx
In this construction, M/2 machines will have the function f1 as their N local component functions (and hence their local objective functions) and the other M/2 machines will have the function f2.
16Again, the functions constructed in this theorem are µ-strongly convex, which is stronger than µ-PL required in Definition 1. Also, our functions satisfy Assumption 2 with ν = 0.
Then, B local RR updates in each machine corresponds to B updates using either f1 or f2. If we start from xmiB = yi, the B local updates on machine m result in
xm(i+1)B = { xmiB + ητB if machine m has f1, (1− 2ηµ)BxmiB − ητ ∑B−1 j=0 (1− 2ηµ)j if machine m has f2.
Taking the average of the M machines, we get that
yi+1 = 1
M M 2 (yi + ητB) + M 2 (1− 2ηµ)Byi − ητ B−1∑ j=0 (1− 2ηµ)j 
= 1
2
( 1 + (1− 2ηµ)B ) yi + ητ
2 B − B−1∑ j=0 (1− 2ηµ)j  .
Since there are total NKB such communication rounds over K epochs, at the end of the run we have
yK,NB =
( 1
2
( 1 + (1− 2ηµ)B ))NKB y0
+ ητ
2 B − B−1∑ j=0 (1− 2ηµ)j  NKB −1∑ l=0 ( 1 2 ( 1 + (1− 2ηµ)B ))l
= ητ
2
( B − 1− (1− 2ηµ) B
2ηµ
) 1− ( 1 2 ( 1 + (1− 2ηµ)B ))NK B
1− 12 (1 + (1− 2ηµ)B) , (176)
where we used initialization y0 = 0. Having defined the function and calculated its last iterate (176), let us now handle the two step-size regimes separately.
We first consider 18µNK ≤ η ≤ 1 8µB . In this case, we exploit the fact that
1− 2ηµB + η2µ2B2 ≤ (1− 2ηµ)B ≤ 1− 2ηµB + 4η2µ2B2, (177)
when 0 ≤ η ≤ 18µB . To see why, consider substituting z := 2ηµ. Then h1(z) := 1−Bz +B 2z2 − (1 − z)B has h′′1(z) ≥ 0 on z ∈ [0, 1], h′1(0) = 0, and h1(0) = 0, implying that h1(z) ≥ 0 on z ∈ [0, 1]. On the other hand, let h2(z) := 1 − Bz + B 2z2 4 − (1 − z) B . If B = 2, then h2 ≡ 0. If B > 2, then it can be checked that h2(z) ≤ 0 for small enough interval [0, 14B ]. Using (177) on (176),
yK,NB = ητ
2
( B − 1− (1− 2ηµ) B
2ηµ
) 1− ( 1 2 ( 1 + (1− 2ηµ)B ))NK B
1− 12 (1 + (1− 2ηµ)B)
≥ ητ 2
( B − 1− (1− 2ηµB + η 2µ2B2)
2ηµ
) 1− ( 1 2 ( 1 + 1− 2ηµB + 4η2µ2B2 ))NK B
1− 12 (1 + 1− 2ηµB + η2µ2B2)
= ητ
2
( ηµB2
2
) 1− ( 1− ηµB + 2η2µ2B2 )NK B
ηµB − 12η2µ2B2
≥ η 2µτB2
4 ·
1− ( 1− ηµB + 14ηµB )NK B
ηµB
≥ ητB 4
( 1− ( 1− 3ηµB
4
)NK B ) ,
where we used ηµB ≤ 18 . Now, substituting η ≥ 1 8µNK to above, we obtain
yK,NB ≥ ητB
4
( 1− ( 1− 3ηµB
4
)NK B )
≥ τB 32µNK
( 1− ( 1− 3B
32NK
)NK B )
≥ (1− e −3/32)τB
32µNK .
Therefore, F (yK,NB ) = Ω ( τ2B2 µN2K2 ) .
Next, consider 18µB ≤ η ≤ 1 µ . We take a close look at the term that appears in (176):
B − B−1∑ j=0 (1− 2ηµ)j = B − 1− (1− 2ηµ) B 2ηµ .
For this term, we would like to find a lower bound which holds for all η ∈ [ 18µB , 1 µ ]. To this end, consider substituting z := 2ηµ. Then, the function
h3(z) := B − B−1∑ j=0 (1− z)j (178)
is increasing on z ∈ [ 14B , 1], and we have
h3 ( 1 4B ) = B − 4B ( 1− ( 1− 14B )B) ≤ h3(z), for all z ∈ [ 14B , 1]. Using B ≥ 2, h3( 14B ) can be lower-bounded as
h3 ( 1 4B ) ≥ B − 4B ( 1− ( 1− 18 )2) = B16 .
Next, for z ≥ 1, the derivative of h3(z) := B − ∑B−1 j=0 (1 − z)j = B − 1−(1−z)B z is h ′ 3(z) = 1−(1−z)B−1((B−1)z+1) z . Since B ≥ 2 is assumed to be even, it is easy to check that h ′ 3(z) ≥ 0 for z ≥ 1, which means that h3 keeps increasing on [1, 2]. Therefore, we conclude that B 16 ≤ h3 ( 1 4B ) ≤ h3(z), for all z ∈ [ 1 4B , 2 ] ,
and hence
B − B−1∑ j=0 (1− 2ηµ)j ≥ B 16 ,
for all 18µB ≤ η ≤ 1 µ . Using yK,NB from (176), we get
yK,NB = ητ
2 B − B−1∑ j=0 (1− 2ηµ)j  NKB −1∑ l=0 ( 1 2 ( 1 + (1− 2ηµ)B ))l
≥ ητ 2 B − B−1∑ j=0 (1− 2ηµ)j  ≥ ητB 32 ≥ τ 256µ .
Here, we used the fact that ∑NK
B −1 l=0 ( 1 2 ( 1 + (1− 2ηµ)B ))l ≥ ( 12 (1 + (1− 2ηµ)B))0 = 1. Hence, we obtain F (yK,NB ) = Ω( τ2 µ ), finishing the proof.
I.2 LOWER BOUND FOR η ≤ 18µNK AND η ≥ 1 µ
We now conclude with the “extreme” step-size regimes. We consider the same function F as in the previous subsection, but with a different initialization y0 = τµ .
For F , recall from (176) that
yK,NB =
( 1
2
( 1 + (1− 2ηµ)B ))NKB y0
+ ητ
2 B − B−1∑ j=0 (1− 2ηµ)j  NKB −1∑ l=0 ( 1 2 ( 1 + (1− 2ηµ)B ))l . (179)
This time, we want to lower-bound the second term on the RHS of (179) with zero and focus on the first term. To this end, we revisit our discussion on h3 (178). It is easy to check that h3 is in fact increasing on the entire [0,∞), and h3(0) = 0. This shows B − ∑B−1 j=0 (1 − 2ηµ)j ≥ 0 for any
η ≥ 0. Next, since B is even, 12 ( 1 + (1− 2ηµ)B ) ≥ 0 for any η ≥ 0. This gives
yK,NB ≥ ( 1
2
( 1 + (1− 2ηµ)B ))NKB y0 = ( 1
2
( 1 + (1− 2ηµ)B ))NKB τ µ . (180)
First, consider the interval 0 ≤ η ≤ 18µNK . Recall from (177) that for this η,
(1− 2ηµ)B ≥ 1− 2ηµB + η2µ2B2 ≥ 1− 2ηµB, (181)
so yK,NB ≥ ( 1
2
( 1 + (1− 2ηµ)B ))NKB τ µ ≥ (1− ηµB)NKB τ µ ≥ ( 1− B 8NK )NK B τ µ ≥ 7τ 8µ ,
since NKB ≥ 1. Hence, F (yK,NB ) = Ω( τ2 µ ).
Finally, if η ≥ 1µ , then we have (1− 2ηµ) B ≥ 1, so
yK,NB ≥ ( 1
2
( 1 + (1− 2ηµ)B ))NKB τ µ ≥ τ µ .
As a result, F (yK,NB ) = Ω( τ2 µ ).

Review summaries:
This paper proposes two variants of stochastic gradient algorithms without replacement. For smooth functions satisfying the PŁ condition, the authors showed that the proposed shuffling-based variants converge faster than their with-replacement counterparts (for the case with large number of epochs). Moreover, the authors showed that their provided convergence analysis is tight for the case with not large number of epochs. 

This paper studies two popular algorithms: local SGD and minibatch SGD using the data shuffling technique, which means sampling without replacement. Authors provide analysis under PL condition and show that in some cases these methods with shuffling outperform classical local SGD and minibatch SGD. Additionally, this paper provides lower bounds for minibatch RR and local RR in homogeneous and heterogeneous settings. At the end of the paper, a new technique of using synchronized permutations is proposed. In an almost homogeneous setting, this approach can improve rates.