Year: 2021

Abstract:
Training example order in SGD has long been known to affect convergence rate. Recent results show that accelerated rates are possible in a variety of cases for permutation-based sample orders, in which each example from the training set is used once before any example is reused. In this paper, we develop a broad condition on the sequence of examples used by SGD that is sufficient to prove tight convergence rates in both strongly convex and non-convex settings. We show that our approach suffices to recover, and in some cases improve upon, previous state-of-the-art analyses for four known example-selection schemes: (1) shuffle once, (2) random reshuffling, (3) random reshuffling with data echoing, and (4) Markov Chain Gradient Descent. Motivated by our theory, we propose two new example-selection approaches. First, using quasi-Monte-Carlo methods, we achieve unprecedented accelerated convergence rates for learning with data augmentation. Second, we greedily choose a fixed scan-order to minimize the metric used in our condition and show that we can obtain more accurate solutions from the same number of epochs of SGD. We conclude by empirically demonstrating the utility of our approach for both convex linear-model and deep learning tasks.

heading: 1 INTRODUCTION
text: To minimize a differentiable function f : Rd → R, stochastic gradient descent (SGD) iteratively updates a parameter vector w ∈ Rd starting at some w0 by running
wt+1 = wt − αt∇f(wt;xt), (1)
where αt is the step size at iteration t, xt is a data example (often a minibatch of data) chosen by some process—typically by subsampling a training dataset—for SGD to use at iteration t, and ∇f(wt;xt) is an example gradient, which we hope will be a good approximation for the gradient of the objective ∇f(wt). In the standard setup, the x’s are drawn from a dataset D of size n, and f(w) = 1n ∑ x∈D f(w;x), which is often referred to as Empirical Risk Minimization (ERM). The order in which the example sequence x0, x1, . . . are chosen is known to affect the convergence of SGD. For instance, compare so-called “random reshuffling” to with-replacement sampling: optimizing a strongly convex objective with T total iterations, random reshuffling samples the xt from D without replacement and achieves an accelerated convergence rate of 1/T 2, while with-replacement sampling yields a convergence rate of 1/T (Bottou, 2012; Recht & Ré, 2012; Gürbüzbalaban et al., 2021). Similar accelerated rates have been shown in other settings and for other example orders, such as shuffling the dataset once (Nguyen et al., 2020; Ahn et al., 2020; Mishchenko et al., 2020). However, these analyses have mostly focused only on specific example-selection schemes, and study only the case of ERM-type finite-sum objectives. This does not help us understand how new example orders (e.g., the data echoing method of Choi et al. (2019)) affect the convergence of SGD.
This paper develops a general condition on the example gradients themselves that is sufficient to provide a convergence rate for SGD. Intuitively, our main result is: the convergence rate of SGD depends on how fast the averages of consecutive example gradients∇f(w;xt) converge to the full objective gradient ∇f(w). With-replacement sampling causes the average of m consecutive gradient examples starting at any timestep τ , 1m ∑τ+m−1 t=τ
∇f(w;xt), to converge to∇f(w) at a rate of O(1/m) in terms of the norm squared: we show SGD with any example sequence that fulfills this condition will converge at the same asymptotic rate as with-replacement sampling. Alternatively, if that average
converges at the faster Õ(1/m2) rate typical of Quasi-Monte-Carlo (QMC) (Caflisch, 1998), then we show SGD enjoys the accelerated rate that random reshuffling gets. Our contributions are as follows:
• We propose a new condition on the example gradients—average gradient error—and provide convergence analysis using this general condition for both non-convex and strongly-convex problems. We justify the validity of this condition on synthetic experiments (Section 3). • We show that many commonly used example orderings—shuffle once, random reshuffling, data echoing, Markov Chain Gradient Descent—can be analyzed as special cases under our theoretical results, which match or improve upon their existing rates in the literature (Section 4). • We propose two new algorithms: (1) QMC-based data augmentation that transforms examples via a low-discrepancy sequence, improving generalization; and (2) a greedy algorithm that sorts the examples before each epoch based on our new average gradient error metric (Section 5). • Empirically, we evaluate our two algorithms on several image classification benchmarks including MNIST, CIFAR10/100 and ImageNet. We show with QMC-based data augmentation, a higher validation accuracy can be achieved without hyperparameter tuning—this suggests that QMC may be a good default driver to use with data augmentation for deep learning in general. Meanwhile, the greedy algorithm converges faster both in terms of iteration and wall-clock time (Section 6).

heading: 2 RELATED WORK
text: Example ordering in stochastic optimization. Traditional example ordering in SGD is carried out in a with-replacement fashion, which is used to ensure unbiased estimation of the full gradient (Robbins & Monro, 1951; Bach & Moulines, 2011; Zhang, 2004; Bottou et al., 2018; Drori & Shamir, 2020). Significant attention has been paid to importance sampling with respect to various measures, such as Lipschitz constants (Schmidt et al., 2017; Needell et al., 2014), example gradient norms and bounds (Zhao & Zhang, 2015; Alain et al., 2015; Papa et al., 2015; Lee et al., 2019), individual losses (Kawaguchi & Lu, 2020; Loshchilov & Hutter, 2015), and data heterogeneity (Lu et al., 2021). Without-replacement sampling, however, is more common in practice and empirically allows faster convergence (Bottou, 2012). Among the most popular without-replacement approaches are shuffle once (SO) (Bertsekas, 2011; Gürbüzbalaban et al., 2019) and random reshuffling (RR) (Ying et al., 2017). In theory, Recht & Ré (2012) undertook the first investigation on convergence of RR via the noncommutative arithmetic-geometric mean conjecture, to which subsequent works provide counter examples (Yun et al., 2021; De Sa, 2020). HaoChen & Sra (2019) performed an epoch-wise acceleration analysis on RR while Gürbüzbalaban et al. (2021) considered its convergence over infinite epochs. Safran & Shamir (2020) analyzed the lower bounds for both SO and RR methods, and their results are further polished by Mishchenko et al. (2020) via the Bregman divergence bound. Aside from manual ordering, another line of research focuses on perturbed example ordering from data echoing (Choi et al., 2019; Agarwal et al., 2020).
Quasi-Monte Carlo. Quasi-Monte Carlo (QMC) is a variant of Monte Carlo (MC) methods that uses a low-discrepancy sequence instead of a pseudorandom sequence. QMC has been successfully applied in a wide variety of domains including computer graphics (Keller, 1995), finance (Joy et al., 1996), and computational biology (Cieslak et al., 2008). In machine learning, using QMC in place of MC can significantly improve many techniques including variational inference (Buchholz et al., 2018; Liu & Owen, 2021), feature mapping (Yang et al., 2014; Avron et al., 2016), normalizing flows (Wenzel et al., 2018), deep learning based PDE (Chen et al., 2019), and time series analysis (Philipson et al., 2020). For stochastic optimization, early works like Homem-de Mello (2008) and Pennanen (2005) established the asymptotic convergence of a QMC sequence in terms of the training set size, while Jank (2005) proposed replacing MC with QMC in computing the E-step of the EM algorithm. Similar to our motivation, Buchholz et al. (2018) analyzed the convergence of SGD when samples are drawn using QMC. Their approach differs significantly from ours in that their method draws an independent unbiased length-b QMC sequence for each minibatch, while our examples come from contiguous subsequences of a length-T QMC sequence that is used across all iterations.

heading: 3 EXAMPLE-GRADIENT AVERAGES AND SGD CONVERGENCE
text: In this section, we describe our setup, define the average gradient error condition we are proposing, and state our main result. Our objective is to minimize a continuously differentiable function
f : Rd → R using examples xt from some set X retrieved at each iteration t. We make the usual assumption that both the loss gradient and the example gradients are L-Lipschitz continuous. Assumption 1 (L-Smoothness). For some L <∞, for any u, v ∈ Rd and any example x ∈ X ,
‖∇f(u;x)−∇f(v;x)‖ ≤ L · ‖u− v‖ and ‖∇f(u)−∇f(v)‖ ≤ L · ‖u− v‖.
We propose a new condition on the average gradient error. Informally, this condition bounds how averages of consecutive example gradients approximate the objective gradient. Formally, Assumption 2. In the context of Equation (1), we say that the example sequence x0, x1, . . . is (γ,C,Φ)-concentrating for γ ∈ [1, 2], C > 0, Φ ≥ 0, if for any timestep τ ≥ 0 and any m > 0,∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f(wτ ;xt)−∇f(wτ ) ∥∥∥∥∥ 2 ≤ 1 mγ ( C2 + Φ2‖∇f(wτ )‖2 ) , (2)
where wτ is the weight parameter vector arrived at after τ SGD update steps in Equation (1).
The constants C and Φ here may depend on the total number of iterations T : specifically, they may absorb logarithmic factors like log(τ +m)2s typical in a QMC error bound, where s is the dimension of the sample space. In addition, when the {xt}’s are random, we will show that this assumption holds with high probability. Furthermore, to show this assumption holds, it suffices to show that the inequality holds for any w: our requirement that it holds only for the specific wτ arrived at by SGD is weaker. We can develop intuition about Assumption 2 by considering familiar cases:
• For general with-replacement sampling, the sum in (2) is a sum of independent random variables; so, a concentration argument would yield Assumption 2 with γ = 1 with high probability. • For general without-replacement sampling from a dataset of size n, any whole-epoch subsequences from the sum in (2) will cancel to zero, so we expect that sum to have magnitudeO(n) independent of m, and thus Assumption 2 should hold with γ = 2 (Propositions 1 to 4). • For low-discrepancy sampling, we would expect Assumption 2 to hold with γ = 2, along with a multiplicative log(τ +m)2s term: this is the classic error rate for QMC (Proposition 6).
In Section 4, we will make this intuition rigorous with high probability under the assumption that the error of the example gradients is bounded. This intuition illustrates both how Assumption 2 can cover previously analyzed settings and how it can generalize to cases not previously studied (e.g., QMC data augmentation). Also observe that Assumption 2 is easily adapted to minibatch SGD: if it holds for the single-example case, it should also hold with modified constants for minibatches of size b (consisting of averages of b consecutive example gradients) by substituting m 7→ mb.1 Since this is straightforward, for simplicity of presentation we focus our theoretical analysis on the batch-size-1 case. Our analysis is based on the following key lemma, which bounds the evolution of SGD over an “analysis phase” of m steps. Lemma 1. Suppose our setup satisfies Assumptions 1 and 2 and that we use a constant step size α. For all timesteps τ ≥ 0, let m > 0 be some integer such that 3αm1−γ/2Φ ≤ αm ≤ 16L . Then the objective at timestep τ +m is bounded by
f(wτ+m) ≤ f (wτ )− 14αm‖∇f(wτ )‖ 2 + 2αm1−γC2.
By applying this lemma inductively and choosing the step size appropriately, we can derive convergence rates for SGD in a variety of settings. Note that while here for simplicity we state results for a fixed step size, Lemma 1 can also be used to get essentially the same rates for diminishing step size schemes, which we analyze in Appendices A.2 and A.3. In what follows, we let f∗ be the global minimum of f , let ∆ := f(w0)− f∗, and use Õ to hide logarithmic terms in the problem parameters such as C, Φ, L, ∆, and , while treating γ as a constant. Theorem 1 (Non-convex case). Suppose that our setup satisfies Assumptions 1 and 2, and let > 0 be any target error. Using SGD (1) with a constant step size α = 16L ⌈ (4C/ + 3Φ)2/γ ⌉−1 , the number of steps T needed to achieve mint=0,··· ,T−1‖∇f(wt)‖2 ≤ 2 is at most
T = ⌈ 48L∆ 2 ⌉ · ⌈( 4C + 3Φ )2/γ⌉ = Õ ( C 2/γL∆ 2+2/γ + Φ 2/γL∆ 2 + C 2/γ 2/γ + Φ2/γ ) .
1The only technical subtlety is that the iterates wτ of SGD would vary with minibatch sizes.
Observe that when γ = 1, this gives the standard −4 rate expected for non-convex SGD, and for γ = 2, this gives us the accelerated −3 rate of shuffling methods. We can get an even faster rate if f satisfies the µ-Polyak-Łojasiewicz (PL) condition ‖∇f(w)‖2 ≥ 2µ(f(w)− f∗), which generalizes strong convexity and has been applied before in the study of sample orders (Mishchenko et al., 2020; Ahn et al., 2020): in particular the following theorem holds for µ-strongly convex functions. Theorem 2. Suppose that f satisfies the µ-PL condition and our setup satisfies Assumptions 1 and 2. Let > 0 be any target error, and κ = L/µ be the condition number of the problem. Using SGD (1) with a constant step size α = 16L ⌈ (8C2/(µ 2) + 9Φ2)1/γ
⌉−1, the number of steps T needed to guarantee f(wT )− f∗ ≤ 2 is at most
T = ⌈ 12κ log ( 2∆ 2 )⌉ · ⌈( 8C2 µ 2 + 9Φ 2 )1/γ⌉ = Õ ( C2/γκ µ1/γ 2/γ + κΦ2/γ + κ ) .
Observe that when γ = 1, this recovers the ordinary T = κ −2 rate we usually get for strongly convex SGD, and when γ = 2 we get a faster rate of κ −1. These theorems together show that our Assumption 2 is sufficient to show the convergence of SGD, and the convergence-rate parameter γ of the assumption translates to affect the convergence rate of SGD. This validates our intuition that faster convergence of averages of consecutive example gradients to the full gradient∇f(w) leads to faster convergence of SGD.
Synthetic experiments. We quickly validate these results on a synthetic 10-dimensional strongly-convex problem. The first setting we consider is the expected risk minimization of R(w) = E [(xTw − y)2], where x ∼ N (0, Id) and y |x ∼ N (xTw∗, 1) for some optimal value w∗. We run SGD in an online fashion—at iteration t we draw samples (xt, yt) from the underlying distribution to compute∇R(wt;xt) used in the update. We compare drawing these samples independently at random against drawing using a QMC sequence (Sobol in [0, 1]d) via an inverse transform, using for both cases the same diminishing step size scheme selected to minimize the expected risk for the random-sampling case—the optimal step size scheme for vanilla SGD. The online plot of Fig. 1(a) shows that the convergence rate is strictly superior with QMC, which achieves a O(1/t2) rate compared to the O(1/t) rate of random sampling, which is what Theorem 2 predicts.
We also evaluate the offline setting, where we draw n independent examples from the same distribution to form a training set, and minimize the empirical risk Rn(w) = 1n ∑n i=1(x T iw − yi)2. This corresponds to a least squares problem with optimal solution w∗n. We run SGD in an epoch-wise manner for K epochs, in which we compare sampling from the training set uniformly with replacement (IID Uniform), random reshuffling (RR), shuffle once (SO), and sampling using one QMC sequence in [0, 1] of length T = nK followed by a mapping to example indices (Sobol). In the offline plot of Fig. 1(a), we see that the low-discrepancy methods all yield an accelerated rate compared to sampling with replacement, which again validates our theory. Additional details on the experiment setup, hyperparameters, and diminishing step size derivation can be found in Appendix A.1.3.
With the same synthetic setup, we also verified our bound in Assumption 2 by measuring the average gradient error over a sequence of examples. The fixed point wτ is arbitrarily set to the origin. For the online setting (left of Fig. 1(b)), we use the same set of examples as in the SGD experiments above starting at t = 0. Similarly, in the offline setting, we go through the examples epoch-wise. As we can see, the sample orderings given by QMC, RR and SO (offline only) indeed give us an accelerated rate of decrease in the average gradient error as we increase m (γ ≈ 2), justifying our main assumption.

heading: 4 ANALYSIS OF EXISTING SCAN ORDERS
text: In this section, we illustrate the power of our approach by proving convergence rates for exampleselection methods proposed and analyzed in previous literature. For each method, we show Assumption 2 holds with high probability, and then by applying Theorems 1 and 2 to these results, we show how existing rates for these methods can be recovered and in some ways improved upon. To the best of our knowledge, these are the first high-probability results for shuffle once and random reshuffle for general non-convex optimization. When using a finite training set of examples, we let n denote its size, let the examples be indexed as x(0), x(1), . . . , x(n−1) (to avoid confusion with xt, the example used by SGD at step t), and let T = nK denote the total number of iterations after K epochs. We will also require the following standard assumption that bounds the error of a single example gradient.
Assumption 3. For all examples x ∈ X and points w ∈ Rd, there exists A,B ≥ 0 such that the gradient errors satisfy ‖∇f(w;x)−∇f(w)‖2 ≤ A2 +B2‖∇f(w)‖2.
Shuffle once (SO). In the shuffle-once variant of SGD, a single permutation σ of {0, . . . , n− 1} is chosen uniformly at random at the start, and the examples are used repeatedly in that order: explicitly, xt = x
(σ(t mod n)). One way to analyze shuffle-once is to prove Assumption 2 for permutation-based methods generally. Proposition 1. Let Assumptions 1 and 3 hold, and suppose that we are using a permutation-based method of sampling, that is, any method such that xkn, xkn+1, . . . , xkn+n−1 is a permutation of x(0), x(1), . . . , x(n−1) for all epochs k ≥ 0. Any such method satisfies Assumption 2 with γ = 2, C2 = n2A2 and Φ2 = n2B2.
This immediately lets us recover previous rates up to constant factors, but we can do better. In the case where we are learning over a bounded region, we can prove a stronger result for shuffle once. Proposition 2. Suppose that we are using the shuffle once variant of SGD to learn over a region B ∈ Rd of radius at most R, such that the iterates wt are guaranteed to remain within this region. Assume that for all w ∈ B and all examples x in the training set of size n, Assumption 1 (LSmoothness) and Assumption 3 hold. Then with probability at least 1− p, Assumption 2 holds with γ = 2, C2 = Õ(dA2(n+B2)), and Φ2 = Õ(ndB2).2
In comparison to previous results and to our rate implied by Proposition 1, this improves the dependence from n to √ nd, which is a significant improvement over the best rates for shuffle once available in the literature when the dimension is small relative the training set size. In particular, if we set B = 0 and consider small , our rate in the non-convex case becomes
T = Õ ( AL∆ √ nd
3
) , which implies 2 ≤ Õ ( (AL∆) 2/3(nd) 1/3
T 2/3
) = Õ (( d
nK2
)1/3) .
This rate matches that for shuffle once achieved in Nguyen et al. (2020, Corollary 1) in terms of , up to logarithmic factors. In the µ-PL (or strongly convex) case for B = 0, we again obtain a rate matching that of Nguyen et al. (2020), Ahn et al. (2020), and Mishchenko et al. (2020):
T = Õ ( κA √ nd µ ) , which implies 2 ≤ Õ ( κ2A2nd µT 2 ) = Õ ( κ2A2d µnK2 ) .
Random reshuffling (RR). Random reshuffling is similar to shuffle once, except that a new ordering is chosen at each epoch: sampling the dataset without replacement. Concretely, if σk denotes the permutation used by random reshuffling at the kth epoch, then xt = x(σbt/nc(t mod n)). Proposition 3. Suppose that we are using the random reshuffling variant of SGD. Assume that for all w ∈ Rd and all examples, Assumption 1 (L-Smoothness) and Assumption 3 hold. For some p ∈ (0, 1), set the constant step size to satisfy α ≤ ( max { 1460BnL · log ( 4e2T/p ) , 2nL })−1 . Then with
probability at least 1− p, Assumption 2 holds with γ = 2, C2 = Õ(nA2) and Φ2 = Õ(nB2).
Setting B = 0, our rate in the non-convex case now becomes T = Õ ( AL∆ √ n
3
) , which implies 2 ≤ Õ ( (AL∆) 2/3(n) 1/3
T 2/3
) = Õ (( 1
nK2
)1/3) .
2The probability p is only present in log terms in these expressions, so it does not appear in the Õ.
Here we match the best rate obtained by Mishchenko et al. (2020, Corollary 3) in the small setting. In the µ-PL (or strongly-convex) case, we get
T = Õ ( κA √ n µ ) , which implies 2 ≤ Õ ( κ2A2n µT 2 ) = Õ ( κ2A2 µnK2 ) .
As in shuffle once, here our rate for random reshuffling matches that obtained by Nguyen et al. (2020) and Ahn et al. (2020), as well as Mishchenko et al. (2020) albeit with a slightly worse dependency on κ. It’s worth noting that for simple quadratics, RR and SO are only faster than with-replacement SGD when K & 1/µ (Safran & Shamir, 2021).
Random reshuffling with data echoing. Data echoing is a technique that can be easily implemented in a machine learning training pipeline to increase throughput and improve performance. It was first introduced and tested empirically by Choi et al. (2019) and analyzed by Agarwal et al. (2020). The idea is to perform multiple SGD updates on each example xi (or minibatch) before proceeding to the next. By “echoing” examples we allow more time for upstream data loading and preprocessing, and consequently decrease downstream GPU idle time for gradient computation. For simplicity, we also use a fixed number of echos c as in Agarwal et al. (2020). Concretely, a c-echoed version of a sample order x̂ is given by xt = x̂bt/cc. Data echoing can essentially be applied to any example-ordering scheme, and here we provide one analysis under random reshuffling. The justification for Assumption 2 under random reshuffling with data echoing follows essentially without modification from the c = 1 version in the previous subsection.
Proposition 4. Suppose that we are using the random reshuffling variant of SGD, where each example is echoed c times using the same step size in RR. Under the same assumptions as in Proposition 3, with probability at least 1− p, Assumption 2 holds with γ = 2, C2 = Õ(cnA2) and Φ2 = Õ(cnB2).
It immediately follows that data echoing should get the same convergence rates we showed in Section 4 with A2 and B2 multiplied by c. Although Agarwal et al. (2020) also provided an analysis for data echoing, they require that the examples are sampled independently, rather than the random reshuffling setting that is more commonly-used: as a result, their analysis did not achieve the accelerated −3 rate that shuffled methods enjoy. Another advantage of our analysis is that the proof follows exactly from that of vanilla random reshuffling, from which the constant c simply propagates.
Markov chain gradient descent (MCGD). To illustrate the versatility of Assumption 2, we show how it can be satisfied by a problem where the objective is not a finite sum and where γ 6= 2. Consider f(w) = Eξ∼Ξ[f(w; ξ)] with some underlying distribution Ξ. Running SGD then requires that at each iteration t, we draw ∇f(wt; ξt) where ξt ∼ Ξ; however, sampling from Ξ can be intractable. The method of Markov Chain Gradient Descent addresses this problem by sampling the ξt from the trajectory of a single Markov chain with stationary distribution Ξ (Sun et al., 2018). The intuition is that although at early iterations the ξ’s have not converged to their true distribution, the iterates visited by SGD are also far from the optimum, thus larger approximation error in the early ξ’s is rather harmless. As we continue iterating, the Markov chain will mix as SGD converges. We show that the convergence of MCGD can be bounded in terms of the mixing time of that Markov chain.
Proposition 5. Suppose that we use samples xt from a Markov chain with mixing time tmix. Assume that for all w ∈ R and all examples xt, Assumption 3 holds. Then with probability at least 1 − p, Assumption 2 holds with γ = 1, C2 = Õ(A2t2mix), and Φ2 = Õ(B2t2mix).
It follows that for non-convex optimization in the B = 0 case, our convergence rate is given by T = Õ ( A2L∆t2mix
4
) , which implies 2 ≤ Õ ( Atmix √ L∆
T 1/2
) .
Sun et al. (2018, Theorem 2) use a diminishing step size O(1/tq) to obtain T = O( −2/1−q), where q ∈ (1/2, 1). In contrast, our rate is faster and holds with high probability instead of in expectation.

heading: 5 NEW EXAMPLE-SELECTION METHODS FOR FASTER CONVERGENCE
text: The analysis in Section 4 focused on recovering the convergence rates for SGD with known exampleordering algorithms. In this section, we propose two new example-selection approaches that allow faster convergence: QMC-based data augmentation and greedily minimizing the metric in (2).
Algorithm 1 Example-Ordered SGD via Greedily Minimizing Average Gradient Error Input: step size α, number of iterations T , random projection matrix Π, buffer for gradients estima-
tion: gi ← 0,∀i ∈ {0, · · · , n− 1}, g ← 0, initial weights w0, initial permutation σ0. 1: for t = 0, · · · , T/n− 1 do 2: Initialize : gi ← 0,∀i ∈ {0, · · · , n− 1}; g ← 0; I ← ∅. 3: for i = 0, · · · , n− 1 do 4: Update the model parameters: wtn+i+1 ← wtn+i − α∇f(wtn+i;x(σt(i))). 5: Update the buffers: gσt(i) ← Π∇f(wtn+i;x(σt(i))); g ← g + gσt(i). 6: end for 7: for i = 0, · · · , n− 1 do
8: σt+1(i)← arg min i∈{0,··· ,n−1}\I [∥∥∥∥ ∑ j∈I∪{i} (gj − g/n) ∥∥∥∥2]; I ← I ∪ {σt+1(i)}.
9: end for 10: end for 11: return wT
QMC-based data augmentation. In many scenarios where only limited examples are given, we want to augment the dataset for better generalization. More formally, given a transform function A that takes example x and a random variable ζ uniformly distributed in [0, 1]s as input, where s denotes the dimension of augmentation space, the augmented objective for a dataset D of size n is
f(w) = 1n ∑ x∈D Eζ∼U[0,1]sf(w;A(x, ζ)) = 1n ∑ x∈D ∫ Rs f(w;A(x, ζ)) dζ. (3)
The rationale is that by performing some reasonable random transformation on a given example, we assume the output would be another example that is identically distributed, and the expected value models an infinitely-large training set consisting of such transformed examples. For example, in an image classification task, we could set s = 1 and haveA(x, ζ) output the image x rotated by an angle of 20◦(2ζ − 1), modeling that a slight rotation of an image should preserve its label. QMC can approximate this expectation with a low-discrepancy sequence of ζt drawn from the s-dimensional unit cube [0, 1]s. Examples of such sequences include the Halton and Sobol sequences (Drmota & Tichy, 2006). QMC is especially favorable for data augmentation because s is usually small in most data augmentation methods. We propose to use QMC for data augmentation together with random reshuffling. Concretely, if ζ0, ζ1, . . . is our low-discrepancy sequence and σk denotes the permutation used by random reshuffling in the kth epoch, then we propose to use the example xt = A(x (σbt/nc(t mod n)) , ζbt/nc+σbt/nc(t mod n)). That is, when we sample example i in epoch k, we use the (k + i)th element of the low-discrepancy sequence. (This is not the only reasonable way of combining QMC and RR: it is just one way we found to work well.) In theory we would expect the approximation error here to decay at the rate Õ(1/m2), instead of the O(1/m) of the random sampling that is standard for data augmentation. To prove this rigorously, we make two additional assumptions which are commonly used in analyzing QMC sequences (Aistleitner & Dick, 2014). Assumption 4 (Bounded gradient variation). There exists a constant V > 0 such that for any fixed w ∈ Rd and x ∈ X , the example gradients as a function of ζ under QMC data augmentation, F (ζ) = ∇f(w;A(x, ζ)), has Hardy-Krause variation (Tezuka, 2000) at most V , that is VHK(F ) ≤ V . Assumption 5. The QMC sequence {ζt}t≥0 has low star-discrepancy (Owen, 2003): for all m > 0,
sup a∈[0,1]s ∣∣∣∣ 1mm−1∑ t=0 1{ζt ∈ [0, a)} − s∏ j=1 aj ∣∣∣∣ ≤ CQMC · log(m)sm , where [0, a) = {x ∈ [0, 1]s | 0 ≤ xj < aj , j = 1, . . . , d} for some constant CQMC.
Proposition 6. Suppose that we are using the random reshuffling variant of SGD with QMC data augmentation as described. Assume that for all w ∈ Rd and all examples, Assumptions 1 , 3, 4 and 5 hold for Equation 3. For some p ∈ (0, 1), set the step size to be a constant such that α ≤ ( max{1460BnL · log(4e2T/p ) , 2nL})−1. Then with probability at least 1− p, Assumption 2
holds with γ = 2, C2 = Õ(n2V 2C2QMC log(T )2s + nA2) and Φ2 = Õ(nB2).
Comparing it with Proposition 3, this QMC variant enjoys the same Õ(1/T 2) rate we get for vanilla random reshuffling: to our knowledge, this is the first accelerated rate for learning with data augmen-
tation. In addition to the theoretical advance, we demonstrate in Section 6 that our QMC variant can achieve better validation performance in practice in multiple applications.
Better example ordering via greedy selection. Taking a closer look at Assumption 2, the magnitude of the left-hand side plays a crucial role in the convergence. In the ERM setting, this motivates us to select a permuted example order that minimizes this expression, following which SGD would converge with minimized average gradient error. However, naively constructing such a sequence is tedious as iterating over all the τ,m > 0 and all permutations can be computationally intensive. In light of this, we apply several approximation techniques into the construction and formulate it into Algorithm 1. The first technique is to use stale gradients, i.e., using the gradients computed at each epoch to estimate the sequence used in the next epoch (line 5 in Algorithm 1). The intuition is that based on the smoothness of loss function (Assumption 1), we would expect the stale gradients to approximate the current gradients with tolerant approximation error as long as the step size is reasonably small. The second technique is random projection (line 9): we search a lower-dimensional space, which allows faster construction and reduces memory use (see similar strategies of using smoothness and dimension reduction techniques in (Caflisch, 1998)). Because this is a selection-by-permutation method, under Assumption 3 our greedy selection method will trivially satisfy Assumption 2 with γ = 2, C = nA, and Φ = nB (see Proposition 1).
6 EXPERIMENTS Table 1: Top1 validation accuracy (%) of ResNet18 on ImageNet. The original one is the standard benchmark provided by PyTorch.
Original Uniform (tuned) QMC (untuned)
69.76 70.19 70.48 In this section we evaluate our new algorithms on several deep learning benchmarks. First, we compare QMC-based data augmentation against IID-uniform augmentation on CIFAR10/100 and ImageNet datasets. Second, we compare greedy ordering (Algorithm 1) with RR and SO, and show how to further accelerate it with randomly projected sorting. Other details on the experimental setup can be found in Appendix A.1.1. Each experiment is repeated 10 times with consistent seeds among the algorithms.
QMC-based data augmentation. We start by training ResNet20 on CIFAR10 and CIFAR100, where discrete and continuous data augmentations are applied, respectively. Specifically, CIFAR10 uses random crop and random horizontal flip while CIFAR100 uses an additional random rotation of 15 degrees. To apply the QMC-based data augmentation, we first generate a Sobol sequence of appropriate dimension using the qmcpy package (Choi et al., 2020+), and then replace the pseudorandom sequence used in the original random augmentation pipelines with that Sobol sequence. We run the baseline IID-uniform method with finetuned hyperparameters (weight decay 1e-4), which reproduces the result from He et al. (2016) with an error rate 8.4%. Then we run QMC-base augmentation with the same hyperparameter (untuned) and finetuned counterparts (we perform grid search on weight decay in the range of {1e-4, 2e-4, 3e-4, 4e-4}). From Figure 2 we observe the QMCbased augmentation consistently outperforms the baseline methods, even without hyperparameter tuning. Comparing Figure 2 with He et al. (2016), we observe the QMC-based augmentation allows ResNet20 to reach comparable validation accuracy as ResNet44 while requiring only 40% as many parameters (0.27M vs 0.66M). We run a t-test on these results (in Appendix A.1.2) to show the validation accuracy from the two augmentation methods are different statistically significantly (p-value
p=7e-7 on CIFAR10 and p=0.036 on CIFAR100). Importantly, this improvement comes essentially for free, as generating low discrepancy sequences in low dimension has very little overhead.
We also evaluate our method on fine-tuning ResNet18 on ImageNet. We apply different augmentation methods on a pre-trained model and train for 5 additional epochs with step size 1e-4. We report their Top1 accuracies in Table 6, which shows that our QMC method improves the validation accuracy by 0.3% compared to the same number of epochs of fine-tuning using random sampling.
Better example ordering via greedy selection. In this section we evaluate Algorithm 1 on two benchmarks: Logistic Regression on MNIST and ResNet20 on CIFAR10. As discussed, sorting the stale gradients naively could incur substantial overhead on memory and computation. To mitigate this, we adopt two methods: random projection and QR decomposition. The former is mainly to reduce storage: we obtain a gradient computed at some time in an epoch and project it into a lower-dimensional space before storing it. Classic ways of projection include Gaussian projection or random sparsification: we adopt the latter as it does not require storing the projection matrix, which minimizes the storage cost. After we obtain all the stale gradients, we concatenate them into a matrix and perform QR decomposition before sorting, which allows us to sort in a low-dimensional space while preserving the order of gradients since the inner products between any two tensors will remain the same Gander (1980). We set the target dimension to be of 10% size of the original space. In the spirit of evaluating the applicability of Algorithm 1, we do not perform hyperparameter tuning in this section but reuse the ones tuned in the literature on Random Reshuffling.
We plot the results in Figure 3. In Figure 3(a) we observe greedy algorithms can consistently converge faster than RR and SO epoch-wise with optimizing the same loss function. When QR is used without projection, the algorithm is able to reach higher validation accuracy but converges slower with respect to the wall-clock time. On the other hand, when we use random sparsification additionally, the algorithm converges faster with respect to both epoch and wall-clock time without compromising the validation accuracy. For CIFAR10, we observe the greedy method can converge faster when optimizing the same loss as other baselines, and achieves higher validation accuracy when fine-tuned.

heading: 7 CONCLUSION
text: We present a unified analysis on example orderings used in SGD, which generalizes several widelyused orderings in the literature. We propose a greedy algorithm that allows faster convergence via constructing a better example order with approximate sorting techniques, as well as QMC-based augmentation that achieves higher validation accuracy on multiple benchmarks. One potential future direction is designing example orderings to more efficiently minimize the average gradient errors.

heading: A APPENDIX

heading: TABLE OF CONTENTS
text: A.1 Experiment details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.1.1 Additional details for Section 6 . . . . . . . . . . . . . . . . . . . . . . . 14
A.1.2 t-test for data augmentation results . . . . . . . . . . . . . . . . . . . . . . 14
A.1.3 Small experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.1.4 Optimal step size derivation . . . . . . . . . . . . . . . . . . . . . . . . . 16
A.2 Convergence analysis: diminishing step size . . . . . . . . . . . . . . . . . . . . . 18
A.2.1 Non-convex case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
A.2.2 Strongly-convex case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
A.3 Proof for Lemma 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
A.4 Convergence analysis: constant step size . . . . . . . . . . . . . . . . . . . . . . . 28
A.4.1 Non-convex case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
A.4.2 Strongly-convex case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
A.5 Justifications for Assumption 2 under various example orderings . . . . . . . . . . 30
A.5.1 Arbitrary permutation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
A.5.2 Shuffle once . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
A.5.3 Random reshuffling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
A.5.4 Random reshuffling with data echoing . . . . . . . . . . . . . . . . . . . . 37
A.5.5 Markov chain gradient descent . . . . . . . . . . . . . . . . . . . . . . . . 38
A.5.6 QMC-based data augmentation with random reshuffling . . . . . . . . . . 40
A.6 Miscellaneous lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

heading: A.1 EXPERIMENT DETAILS

heading: A.1.1 ADDITIONAL DETAILS FOR SECTION 6
text: In Section 6, all the training scripts are implemented via PyTorch1.6 and run on a single machine configured with an 2.6GHz 4-core Intel (R) Xeon(R) CPU, 16GB memory and NVIDIA GeForce GTX 1080Ti with CUDA 10.1.
In the example ordering comparison, we use the same seed among different algorithms in the same run so as to guarantee every algorithm works with the same loss function.
In the ImageNet experiment, the standard step size schedule for ImageNet training is starting at 0.1 and decaying by 10 every 30 epochs. The pre-trained model is the trained model at epoch 90. Naturally, our learning rate should be 1e-4 by the same schedule. Other hyperparameters are adopted by the open source implementation: https://github.com/pytorch/examples/tree/ master/imagenet.
In the data augmentation comparison, obviously other augmentation techniques can be used. The strategies we used are taken from open source implementation https://github.com/ akamaster/pytorch_resnet_cifar10 and https://github.com/weiaicunzai/ pytorch-cifar100 that can reproduce the validation accuracy in He et al. (2016), so that our comparison can be consistent with the correct benchmarks. Below we include the convergence plot for the experiment on ImageNet with ResNet18, where each algorithm is repeated three times with seeds uniformly selected from [0, 1000].
A.1.2 t-TEST FOR DATA AUGMENTATION RESULTS
We now perform t-test on the validation accuracy on two data augmentation results, to show the improvement from QMC is statistically significant. We include validation accuracy at epoch 200 for CIFAR10 and best accuracy for CIFAR100 in Table A.1.2. We compute the p-value between IID-Uniform and tuned QMC, and found the p-value on CIFAR10 and CIFAR100 to be 7e-7 and 0.036, respectively. Since they are both smaller than 0.05, we reject the null hypothesis which concludes they are of the same mean.

heading: A.1.3 SMALL EXPERIMENTS
text: For both experiments in Section 3, we generated w∗ from a standard Normal distribution. The minibatch size is 1, and the total number of iterations is 107, with n = 104 over 1000 epochs for the offline setting. For all variants, we use the theoretical step size optimized for SGD with replacement (corresponding to IID Uniform) on this particular problem, given by
αt+1 = αt(1− αt)
1− α2t (d+ 2) with α0 =
‖w∗‖2
‖w∗‖2(d+ 2) + d ,
which are derived in Appendix A.1.4 below. To obtain low-discrepancy samples from the Gaussian distribution, we use the inverse transform method. First we obtain the QMC sequences (Sobol in our
case) ζt ∈ [0, 1]s where s is the appropriate dimension, and then use that in the inverse CDF function of a Gaussian distribution to obtain the corresponding Gaussian sample.
In addition to using synthetic data, we also performed an offline version of the experiment Figure 1(b) on a real dataset, a6a from the LIBSVM repository (Chang & Lin, 2011). The dataset contains n=11220 examples with d=124 features (including bias), and all labels are binary. We use logistic regression with `2 regularization (with λ=1e-4) as the empirical risk:
Rn(w) = 1
n n∑ i=1 log(1 + exp(−yixTiw)) + λ 2 ‖w‖2.
To measur the left hand side of Assumption 2, we again fix wτ to be at the origin, and we take the number of epochs to be 1000. The results are in Figure 5 for one run with an arbitrarily-set seed, and we draw similar conclusions as observed with synthetic data (see Figure 1(b) in the main paper).

heading: A.1.4 OPTIMAL STEP SIZE DERIVATION
text: In this section, we derive the optimal step-size sequence used in our synthetic toy example from Section 3. Recall our setup x ∼ N (0, Id), y = xTw∗ + where w∗ is the true model parameters that we are trying to recover, and ∼ N (0, 1) is some intrinsic error. The goal is to minimize the expected risk
f(w) = 1 2 E [ (xTw − y)2 ] .
By using the tower rule of expectation conditioning on x, this objective can be written as
f(w) = 1
2 (‖w − w∗‖2 + 1).
The full gradient and the example gradient given a pair (x, y) drawn from the above distribution are
∇f(w) = w − w∗, ∇f(w;x, y) = (xTw − y)x = (xT(w − w∗)− )x.
Define the gradient error to be
∆g := ∇f(w;x, y)−∇f(w) = (xxT − Id)(w − w∗)− x. (4) Let us first analyze the expected gradient error conditioned on any randomness in w that might arise from the sampling history:
E [ ‖∆g‖2 ] = (w − w∗)TE [ (xxT − Id)2 ] − 2E[ xT(xxT − Id)](w − w∗) + E[ 2xTx].
The last term is simply E[ 2xTx] = d since is independent of x. The middle term cancels to 0 again by independence and zero-mean of . For the first term,
E [ (xxT − Id)2 ] = E[xxTxxT − 2xxT + I] = (d+ 2)Id − 2Id + Id = (d+ 1)Id.
Together, we have
E [ ‖∆g‖2 ] = (d+ 1)‖w − w∗‖2 + d. (5)
Now let us derive an suboptimality gap recursion for the SGD update step using the step size αt:
wt+1 − w∗ = wt − αt∇f(wt;x, y)− w∗
= (wt − w∗) + αt∇f(wt)− αt∇f(wt)− αt∇f(wt;x, y) = (1− αt)(wt − w∗)− αt∆gt ,
where ∆gt is the expression in Eq. (4) with wt as input. The expected suboptimality gap conditioned on all randomness up to timestep t is then
E [ ‖wt+1 − w∗‖2 ] = (1− αt)2E [ ‖wt − w∗‖2 ] − 2(1− αt)E[((wt − w∗)T∆gt ] + α2tE [ ‖∆gt ‖ 2 ] .
Observe that under iid uniform sampling, E[∆gt ] = 0, and using Eq. (5) gives us E [ ‖wt+1 − w∗‖2 ] = (1− αt)2‖wt − w∗‖2 + α2t (d+ 1)‖wt − w∗‖ 2 + α2td.
Taking an expectation over the entire history to remove conditioning and let ρt := E [ ‖wt − w∗‖2 ] ,
ρt+1 = (1− αt)2ρt + α2t (d+ 1)ρt + α2td = (1− 2αt + α2t (d+ 2))ρt + α2td. (6)
Note that the RHS is convex in αt, we can differentiate to minimize the expected suboptimality at every iteration ρt+1 = E [ ‖wt+1 − w∗‖2 ] :
0 = (−2 + 2αt(d+ 2))ρt + 2αtd⇒ αt = ρt
(d+ 2)ρt + d .
This also gives us an expression for ρt in terms of αt:
ρt = d
α−1t − (d+ 2) ,
which combined with Eq. (6) gives us
ρt+1 = (1− αt)ρt = (1− αt) d
α−1t − (d+ 2) .
Finally, the optimal step-size sequence can be implemented via the following recursion
1
αt+1 = d+ 2 +
d
ρt+1 = d+ 2 + α−1t − (d+ 2) 1− αt
1
α0 = d+ 2 +
d
ρ0 = d+ 2 +
d
‖w0 − w∗‖2 .
It is easy to verify that this is indeed a decreasing sequence. If we initialize at w0 = 0, then our step sizes are given by
αt+1 = αt(1− αt)
1− α2t (d+ 2) with α0 =
‖w∗‖2
‖w∗‖2(d+ 2) + d .

heading: A.2 CONVERGENCE ANALYSIS: DIMINISHING STEP SIZE
text: In the main text, our convergence results focused on the constant step size regime for simplicity of presentation and ease of comparison to prior works, as they often lack diminishing step size results for permutation-based SGD variants. Here we formally present the convergence rate under our main Assumption 2 using a diminishing step size sequence. In Theorems 3 and 4 that we present below, the rate is of the same order as what we obtained under a constant step size (Theorems 1 and 2) up to a factor of log(T ) that typically arises when using a diminishing step size.
Our analysis for the diminishing step size setting is accomplished by breaking the total number of iterations into phases. Suppose we want to run a total of T iterations of SGD updates. We will break the analysis into I number of phases. In each phase i = 0, . . . I − 1, we run K(i)m(i) number of iterations with constant step size α(i). The step size is decayed at the beginning of each phase. For this, we define the triplets (α(i),m(i),K(i)), where { m(i) } is an increasing sequence.
Recall that γ and Φ are constants in Assumption 2. In the non-convex case, the inner interval length is chosen to be
m(i) = d3(Φ + 1)2/γe2i,
where as for the strongly-convex setting it is
m(i) = d3(Φ + 1)2/γ · ei/γe.
For both function classes the diminishing step-size is set to
α(i) = 1
6Lm(i) . (7)

heading: A.2.1 DIMINISHING STEP SIZE: NON-CONVEX CASE
text: Theorem 3 (Diminishing step size: non-convex case). Suppose that our setup satisfies Assumptions 1 and 2, and let > 0 be any target error. After T iterations of SGD (Eq. (1)) with a diminishing step-size (Eq. (7)), we can obtain
min t=0,...,T−1
‖∇f(wt)‖2 ≤ 12 log2 ( 16
3(Φ + 1) T
) L∆ + C2(3(Φ + 1))−2
T γ/1+γ(12(Φ + 1)2)−γ/1+γ − 1 = O
( log(T )
T γ/1+γ
) .
This implies the number of gradient evaluations to achieve mint‖∇f(wt)‖2 ≤ 2 is
T = Õ( −2(1+γ) γ ).
Proof. To apply Lemma 1 to the iterations within a particular phase i, we need
3Φ ≤ mγ/2 and α ≤ 1 6Lm
where we have temporarily dropped the index (i) for convenience. Lemma 1 gives us the following bound on the objective value between any length-m number of iterations within that phase: for τ ∈ {0,m− 1, 2m− 1, . . . , (K − 1)m− 1},
f(wτ+m) ≤ f(wτ )− η
4 ‖∇f(wτ )‖2 +
2C2
η α2m2−γ
≤ f(wτ )− αm
4 ‖∇f(wτ )‖2 + 2C2αm1−γ ,
since the step size is constant throughout one phase. Summing over K such intervals from τ = 0 followed by a telescope,
αm
4 K−1∑ k=0 ‖∇f(wmk)‖2 ≤ f(w0)− f(wmK−1) + 2C2Kαm1−γ .
We now restore our phase index,
α(i)m(i)K(i)
4
1
K(i) K(i)−1∑ k=0 ‖∇f(wmk)‖2 ≤ f(w0)− f(wm(i)K(i)−1) + 2C2K(i)α(i)(m(i))1−γ
= f(w0)− f(wm(i)K(i)−1) + C2
3L K(i)(m(i))−γ ,
where we have chosen the largest possible α(i) = 1/6Lm(i) for all i. Summing over I phases and letting ∆ := f(w0)− f∗,
I−1∑ i=0 α(i)m(i)K(i) 4 1 K(i) K(i)−1∑ k=0 ‖∇f(wm(i)k)‖ 2 ≤ I∆ + C 2 3L I−1∑ i=0 K(i)(m(i))−γ
where the left hand side can be further lower bounded with I−1∑ i=0 α(i)m(i)K(i) 4 min t=0,...,T−1 ‖∇f(wt)‖2 ≤ I∆ + C2 3L I−1∑ i=0 K(i)(m(i))−γ ,
using the inclusion that { 0,m(i) − 1, . . . ,K(i)m(i) − 1 }I−1 i=0 ⊂ {0, . . . , T − 1}. Now let us choose
m(i) = d(3(Φ + 1))2/γe2i
K(i) = d2iγe.
Here the choice of m(i) guarantees that
(m(i)) γ/2 = ( d(3(Φ + 1))2/γe2i )γ/2 = d(3(Φ + 1))2/γeγ/2(2i)γ/2
≥ (3(Φ + 1))2iγ/2 ≥ 3Φ.
Next, observe that for the last term in the previous sum,
I−1∑ i=0 K(i)(m(i))−γ = I−1∑ i=0 d2iγe(d(3(Φ + 1))2/γe2i)−γ
≤ I−1∑ i=0 ( 2iγ + 1 ) ((3(Φ + 1)) 2/γ2i)−γ
≤ I−1∑ i=0 2(2iγ)((3(Φ + 1)) 2/γ2i)−γ
= I−1∑ i=0 2(3(Φ + 1))−22iγ2−iγ = 2I(3(Φ + 1))−2.
Combining with the above,
I−1∑ i=0 α(i)m(i)K(i) 4 min t=0,...,T−1 ‖∇f(wt)‖2 ≤ I∆ + 2C2 3L · I(3(Φ + 1))−2. (8)
Furthermore, using our choice of α(i)
I−1∑ i=0 α(i)m(i)K(i) 4 = 1 6L I−1∑ i=0 K(i).
Re-arranging Eq. (8),
min t=0,...,T−1
‖∇f(wt)‖2 ≤ 6IL∆ + 23C 2 · I(3(Φ + 1))−2∑I−1 i=0 K (i) . (9)
The total number of iterations is given by
T = I−1∑ i=0 K(i)m(i)
= I−1∑ i=0 d2iγe · d(3(Φ + 1))2/γe2i
≥ (3(Φ + 1))2/γ I−1∑ i=0 2(γ+1)i
= (3(Φ + 1)) 2/γ 2 (γ+1)I − 1 2γ+1 − 1 ≥ (3(Φ + 1))2/γ 2 (γ+1)I − 1 2γ+1 .
Solving for I and using lg to denote log2, we obtain I ≤ ( 1
γ + 1
) lg ( 2γ+1
(3(Φ + 1))2/γ T + 1
) ≤ ( 1
γ + 1
) lg ( 2γ+2
(3(Φ + 1))2/γ T
) .
Moreover, we can also upper bound T using similar arguments,
T = I−1∑ i=0 d2iγe · d(3(Φ + 1))2/γe2i
≤ I−1∑ i=0 2(2iγ) · 2(3(Φ + 1))2/γ2i
= 4(3(Φ + 1)) 2/γ I−1∑ i=0 2(γ+1)i
= 4(3(Φ + 1)) 2/γ 2
(γ+1)I
2γ+1 − 1 ,
which gives
I ≥ ( 1
γ + 1
) lg  2γ+1 − 14(3(Φ + 1))2/γ︸ ︷︷ ︸ =:Γ T  We are left to bound
I−1∑ i=0 K(i) = I−1∑ i=0 d2iγe ≥ I−1∑ i=0 2iγ = 2γI − 1 2γ − 1
≥ 2 γ/1+γ lg(ΓT ) − 1
2γ
= (ΓT )γ/1+γ − 1
2γ
=
( 2γ+1−1 4(3(Φ+1))2/γ T )γ/1+γ − 1
2γ
using our lower bound for I . Substituting this and the upper bound for I into Eq. (9),
min t=0,...,T−1
‖∇f(wt)‖2 ≤ 6IL∆ + 23C 2 · I(3(Φ + 1))−2∑I−1 i=0 K (i)
≤ ( 1
γ + 1
) lg ( 2γ+2
(3(Φ + 1))2/γ T
) 6L∆ + 23C
2(3(Φ + 1))−2( 2γ+1−1
4(3(Φ+1)) 2/γ
T )γ/1+γ −1
2γ ≤ 12 lg ( 16
3(Φ + 1) T
) L∆ + C2(3(Φ + 1))−2(
2γ+1−1 4(3(Φ+1))2/γ
T )γ/1+γ − 1
≤ 12 lg ( 16
3(Φ + 1) T
) L∆ + C2(3(Φ + 1))−2
T γ/1+γ(12(Φ + 1)2)−γ/1+γ − 1 = O ( log(T )
T γ/1+γ
) ,
which yields our convergence rate in the nonconvex setting.

heading: A.2.2 DIMINISHING STEP SIZE: STRONGLY-CONVEX CASE
text: Theorem 4 (Diminishing step size: strongly-convex case). Suppose f is µ-Polyak-Łojasiewicz (PL), and that our setup satisfies Assumptions 1 and 2 Let > 0 be any target error, and κ = L/µ be the condition number of the problem. After T iterations of SGD (Eq. (1)) with a diminishing step-size (Eq. (7)), we can obtain
f(wT )− f∗ ≤ (
e1/γ − 1 2d12κe(3(Φ + 1))2/γ
)−γ( ∆ + 4C2
µ
(3(Φ + 1))−2
e(1− ρ)d12κe
) = O ( 1
T γ
) .
This implies the number of gradient evaluations to achieve f(wT )− f∗ ≤ 2 is
T = O( −2/γ).
Proof. We will begin with the same analysis technique as used in the non-convex case. To apply Lemma 1 to the iterations within a particular phase i, we need
3Φ ≤ mγ/2 and α ≤ 1 6Lm
dropping the index (i) for convenience and will re-introduce it later when appropriate. Lemma 1 gives us the following bound on the objective value between any length-m number of iterations within that phase: for τ ∈ {0,m− 1, . . . , (K − 1)m− 1},
f(wτ+m) ≤ f(wτ )− η
4 ‖∇f(wτ )‖2 +
2C2
η α2m2−γ
≤ f(wτ )− αm
4 ‖∇f(wτ )‖2 + 2C2αm1−γ ,
since the step size is constant throughout one phase. Strong-convexity (or the Polyak-Łojasiewicz (PL) inequality) of f implies ‖∇f(w)‖2 ≥ 2µ(f(w) − f∗) for all w ∈ Rd. Using this while subtracting f∗ on both sides leads to
f(wτ+m)− f∗ ≤ (
1− αmµ 2
) (f(wτ )− f∗) + 2C2αm1−γ .
Applying this recursively K times gives
f(wmK)− f∗ ≤ (
1− αmµ 2
)K (f(w0)− f∗) + 2C2αm1−γ K−1∑ k=0 ( 1− αmµ 2 )k ≤ (
1− αmµ 2
)K (f(w0)− f∗) + 2C2αm1−γ ∞∑ k=0 ( 1− αmµ 2 )k = (
1− µ 12L
)K (f(w0)− f∗) + 2C2αm1−γ 2
αmµ
= ( 1− 1
12κ
)K (f(w0)− f∗) + 4C2m−γµ−1,
where we have used
α ≤ 1 6Lm =⇒ αmµ 2 ≤ µ 12L < 1 as µ ≤ L.
We now restore the phase index. Letting Ti := ∑i−1 i′=0m (i′)K(i ′) be the total number of iterations passed after i phases so that T = TI , and ∆t := f(wt)− f∗,
f(wTi+1)− f∗ ≤ (
1− 1 12κ
)K(i) ∆Ti + 4C2
µ (m(i))−γ .
Now let us choose
m(i) = d(3(Φ + 1))2/γ · ei/γe
K(i) = K = d12κe = d −1 log(1− 1/12κ) e ≥ 12.
Here the choice of m(i) guarantees that
(m(i)) γ/2 = ( d(3(Φ + 1))2/γ · ei/γe )γ/2 ≥ ( (3(Φ + 1)) 2/γ )γ/2 e i/2 ≥ 3Φ.
Using the constant K across all phases, our recursion can be simplified to f(wTi+1)− f∗ ≤ (
1− 1 12κ
)K ∆Ti + 4C2
µ (m(i))−γ .
Since this holds for any i = 0, . . . , I − 1, applying recursion over I phases yields,
f(wT )− f∗ ≤ (1− ρ)KI∆0 + 4C2
µ I−1∑ i=0 (1− ρ)Ki(m(I−i))−γ (10)
where we defined ρ := 1/12κ. Observe that for our choice of m(i) and K,
T = I−1∑ i=0 K(i)m(i) = K I−1∑ i=0 d(3(Φ + 1))2/γ · ei/γe
≤ K I−1∑ i=0 2(3(Φ + 1)) 2/γ · ei/γ
= 2K(3(Φ + 1)) 2/γ I−1∑ i=0 e i/γ
= 2K(3(Φ + 1))2/γ e1/γ − 1 (e I/γ − 1) ≤ 2K(3(Φ + 1)) 2/γ e1/γ − 1 e I/γ ,
which implies
I ≥ γ log (
e1/γ − 1 2K(3(Φ + 1))2/γ T
) ,
and so
(1− ρ)KI ≤ exp(−ρKI) ≤ exp ( −ρKγ · log ( e1/γ − 1
2K(3(Φ + 1))2/γ T
))
note that
ρKγ = γρd −1 log(1− ρ) e ≥ γ.
This gives us
(1− ρ)KI ≤ (
e1/γ − 1 2K(3(Φ + 1))2/γ
)−γ T−γ . (11)
Substituting this into Eq. (10),
f(wT )− f∗ ≤ ∆0 (
e1/γ − 1 2K(3(Φ + 1))2/γ
)−γ T−γ + 4C2
µ I−1∑ i=0 (1− ρ)Ki(m(I−i))−γ . (12)
For the last term,
I−1∑ i=0 (1− ρ)Ki(m(I−i))−γ = I−1∑ i=0 (1− ρ)Kid(3(Φ + 1))2/γ · eI−i/γe−γ
≤ (3(Φ + 1))−2 I−1∑ i=0 (1− ρ)Kiei−I
= (3(Φ + 1))−2e−I I−1∑ i=0 (e(1− ρ)K)i
= (3(Φ + 1))−2e−I (e(1− ρ)K)I
e(1− ρ)K
= (3(Φ + 1))−2 (1− ρ)KI
e(1− ρ)K .
Substituting these into Eq. (12), and using Eq. (11) to bound (1− ρ)KI gives us the convergence rate for strongly-convex functions using a diminishing step size:
f(wT )− f∗ ≤ ∆(1− ρ)KI + 4C2
µ (3(Φ + 1))−2
(1− ρ)KI
e(1− ρ)K ≤ (1− ρ)KI ( ∆ + 4C2
µ
) (3(Φ + 1))−2
e(1− ρ)K ≤ (
e1/γ − 1 2d12κe(3(Φ + 1))2/γ
)−γ( ∆ + 4C2
µ
(3(Φ + 1))−2
e(1− ρ)d12κe
) T−γ
= O ( 1
T γ
) .

heading: A.3 PROOF FOR LEMMA 1
text: We now prove the lemma that bounds the evolution of SGD over a length-m analysis phase. This lemma is the key to our convergence analyses in both the diminishing and constant step size regimes. Before we proceed, we first state and prove the following bound on the gradient error when scaled by a nonincreasing sequence. Lemma 2. Suppose Assumption 2 holds. If {ρt} is a deterministic, nonincreasing, and nonnegative sequence, then∥∥∥∥∥ τ+m−1∑ t=τ ρt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 ≤ ρ2τ ·m2−γ · ( C2 + Φ2‖∇f(wτ )‖2 ) .
Proof. Let βt = ρt − ρt+1 for t ∈ {τ, τ + 1, . . . , τ +m− 2}, and let βτ+m−1 = ρt+m−1. Observe that these are all nonnegative, and∥∥∥∥∥ τ+m−1∑ t=τ ρt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2
= ∥∥∥∥∥ τ+m−1∑ t=τ τ+m−1∑ k=t βk (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2
= ∥∥∥∥∥ τ+m−1∑ k=τ k∑ t=τ βk (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2
= ρ2τ ∥∥∥∥∥ τ+m−1∑ k=τ βk ρτ k∑ t=τ (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 .
Applying Jensen’s inequality using ∑τ+m−1 k=τ βk = ρτ ,∥∥∥∥∥ τ+m−1∑ t=τ ρt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2
≤ ρ2τ τ+m−1∑ k=τ βk ρτ ∥∥∥∥∥ k∑ t=τ (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2
= ρ2τ τ+m−1∑ k=τ βk ρτ (k − τ + 1)2 ∥∥∥∥∥ 1k − τ + 1 k∑ t=τ (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 .
Applying Assumption 2 on the squared norm,∥∥∥∥∥ τ+m−1∑ t=τ ρt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2
≤ ρ2τ τ+m−1∑ k=τ βk ρτ (k − τ + 1)2 · 1 (k − τ + 1)γ ( C2 + Φ2‖∇f(wτ )‖2 ) = ρ2τ
τ+m−1∑ k=τ βk ρτ · (k − τ + 1)2−γ · ( C2 + Φ2‖∇f(wτ )‖2 ) ≤ ρ2τ
τ+m−1∑ k=τ βk ρτ ·m2−γ · ( C2 + Φ2‖∇f(wτ )‖2 ) = ρ2τ ·m2−γ · ( C2 + Φ2‖∇f(wτ )‖2 ) ,
which is what we wanted to show.
We now recall the statement for Lemma 1. In the main paper, we state the lemma with constant step size. In the following proof, we prove it holds for non-decreasing step size, where constant step size will hold naturally as a special case. For this, we define η := ∑τ+m−1 t=τ αt. Lemma 1. Suppose our setup satisfies Assumptions 1 and 2 and that we use a constant step size α For all timesteps τ ≥ 0, let m > 0 be some integer such that 3ατm1−γ/2Φ ≤ η ≤ 16L . Then the objective at timestep τ +m is bounded by
f(wτ+m) ≤ f (wτ )− η
4 ‖∇f(wτ )‖2 +
2 η α2τm 2−γC2.
Proof. Let η and g denote
η := τ+m−1∑ t=τ αt and g := 1 η τ+m−1∑ t=τ αt∇f(wt;xt).
This means that wτ+m = wτ − ηg. From Assumption 1 (L-Smoothness),
f(wτ+m) = f (wτ − ηg)
≤ f (wτ )− η 〈∇f(wτ ), g〉+ η2L
2 ‖g‖2
= f (wτ )− η
2 ‖∇f(wτ )‖2 −
η 2 ‖g‖2 + η 2 ‖g −∇f(wτ )‖2 +
η2L
2 ‖g‖2
≤ f (wτ )− η
2 ‖∇f(wτ )‖2 +
η 2 ‖g −∇f(wτ )‖2,
where the last inequality follows because we assumed ηL ≤ 1/6 (< 1). Next, observe that η
2 ‖g −∇f(wτ )‖2
= 1
2η ‖ηg − η∇f(wτ )‖2
= 1
2η ∥∥∥∥∥ τ+m−1∑ t=τ αt (∇f(wt;xt)−∇f(wτ )) ∥∥∥∥∥ 2
≤ 1 η ∥∥∥∥∥ τ+m−1∑ t=τ αt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 + 1 η ∥∥∥∥∥ τ+m−1∑ t=τ αt (∇f(wt;xt)−∇f(wτ ;xt)) ∥∥∥∥∥ 2 .
To bound the second term,
1
η ∥∥∥∥∥ τ+m−1∑ t=τ αt (∇f(wt;xt)−∇f(wτ ;xt)) ∥∥∥∥∥ 2 = η ∥∥∥∥∥ τ+m−1∑ t=τ αt η (∇f(wt;xt)−∇f(wτ ;xt)) ∥∥∥∥∥ 2
≤ η τ+m−1∑ t=τ αt η ‖(∇f(wt;xt)−∇f(wτ ;xt))‖2
≤ τ+m−1∑ t=τ αtL 2‖wt − wτ‖2.
Combining the above,
η 2 ‖g −∇f(wτ )‖2 ≤ 1 η ∥∥∥∥∥ τ+m−1∑ t=τ αt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 + τ+m−1∑ t=τ αtL 2‖wt − wτ‖2. (13)
Applying Lemma 2 on the first term in Eq. (13) gives
1
η ∥∥∥∥∥ τ+m−1∑ t=τ αt (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 ≤ 1 η α2τm 2−γ ( C2 + Φ2‖∇f(wτ )‖2 ) .
For the second term in Eq. (13),
τ+m−1∑ t=τ αtL 2‖wt − wτ‖2 = τ+m−1∑ t=τ αtL 2 ∥∥∥∥∥ t−1∑ u=τ αu∇f(wu;xu) ∥∥∥∥∥ 2
≤ 3 τ+m−1∑ t=τ αtL 2 ∥∥∥∥∥ t−1∑ u=τ αu (∇f(wu;xu)−∇f(wτ ;xu)) ∥∥∥∥∥ 2
+ 3 τ+m−1∑ t=τ αtL 2 ∥∥∥∥∥ t−1∑ u=τ αu (∇f(wτ ;xu)−∇f(wτ )) ∥∥∥∥∥ 2
+ 3 τ+m−1∑ t=τ αtL 2 ∥∥∥∥∥ t−1∑ u=τ αu∇f(wτ ) ∥∥∥∥∥ 2
≤ 3 τ+m−1∑ t=τ αtL 4 ( t−1∑ u=τ αu ) t−1∑ u=τ αu‖wu − wτ‖2
+ 3 τ+m−1∑ t=τ αtL 2 · α2τm2−γ ( C2 + Φ2‖∇f(wτ )‖2 )
+ 3 τ+m−1∑ t=τ αtL 2 ( t−1∑ u=τ αu )2 ‖∇f(wτ )‖2,
where the second term in the last inequality follows from Lemma 2. Continuing,
τ+m−1∑ t=τ αtL 2‖wt − wτ‖2 ≤ 3 τ+m−1∑ t=τ αtL 4η t−1∑ u=τ αu‖wu − wτ‖2
+ 3 τ+m−1∑ t=τ αtL 2α2τm 2−γ ( C2 + Φ2‖∇f(wτ )‖2 ) + 3η3L2‖∇f(wτ )‖2
= 3η2L2 τ+m−1∑ t=τ αtL 2‖wt − wτ‖2
+ 3L2ηα2τm 2−γ ( C2 + Φ2‖∇f(wτ )‖2 ) + 3η3L2‖∇f(wτ )‖2.
This implies that
(1− 3η2L2) τ+m−1∑ t=τ αtL 2‖wt − wτ‖2 ≤ 3L2C2ηα2τm2−γ
+ ( 3η3L2 + 3L2ηα2τm 2−γΦ2 ) ‖∇f(wτ )‖2.
And because by our assumption, α2τm 2−γΦ2 ≤ η
2
9 , it follows that
(1− 3η2L2) τ+m−1∑ t=τ αtL 2‖wt − wτ‖2 ≤ 3L2C2ηα2τm2−γ + 10 3 η3L2‖∇f(wτ )‖2.
Since we also assumed η ≤ 16L , it follows that 3η 2L2 ≤ 112 , and so(
1− 1 12 ) τ+m−1∑ t=τ αtL 2‖wt − wτ‖2 ≤ 3L2C2ηα2τm2−γ + 5 54 η‖∇f(wτ )‖2,
which gives
τ+m−1∑ t=τ αtL 2‖wt − wτ‖2 ≤ 12 11 · 3L2C2ηα2τm2−γ + 12 11 · 5 54 η‖∇f(wτ )‖2
= 36
11 L2C2ηα2τm
2−γ + 10
99 η‖∇f(wτ )‖2.
So, putting this all together,
f(wτ+m) ≤ f (wτ )− η
2 ‖∇f(wτ )‖2 +
1 η α2τm
2−γ ( C2 + Φ2‖∇f(wτ )‖2 ) + 36
11 L2C2ηα2τm
2−γ + 10
99 η‖∇f(wτ )‖2
= f (wτ )− η
2 ‖∇f(wτ )‖2 +
1 η α2τm 2−γC2
+ 1
η α2τm 2−γΦ2‖∇f(wτ )‖2
+ 36
11 L2C2ηα2τm
2−γ + 10
99 η‖∇f(wτ )‖2
≤ f (wτ )− η
2 ‖∇f(wτ )‖2 +
1 η α2τm 2−γC2
+ 1
9 η‖∇f(wτ )‖2
+ 4L2C2ηα2τm 2−γ +
10 99 η‖∇f(wτ )‖2
= f (wτ )− η
2 ‖∇f(wτ )‖2 +
1 η α2τm 2−γC2
+ 4L2C2ηα2τm 2−γ +
21 99 η‖∇f(wτ )‖2
≤ f (wτ )− η
2 ‖∇f(wτ )‖2 +
1 η α2τm 2−γC2
+ 1
4η C2α2τm
2−γ + 1
4 η‖∇f(wτ )‖2
≤ f (wτ )− η
4 ‖∇f(wτ )‖2 +
2 η α2τm 2−γC2.
This is what we wanted to show.

heading: A.4 CONVERGENCE ANALYSIS: CONSTANT STEP SIZE

heading: A.4.1 CONSTANT STEP SIZE: NON-CONVEX CASE
text: Theorem 1 (Non-convex case). Suppose that our setup satisfies Assumptions 1 and 2, and let > 0 be any target error. Using SGD (1) with a constant step size α = 16L ⌈ (4C/ + 3Φ)2/γ ⌉−1 , the number of steps T needed to achieve mint=0,··· ,T−1‖∇f(wt)‖2 ≤ 2 is at most
T = ⌈ 48L∆ 2 ⌉ · ⌈( 4C + 3Φ )2/γ⌉ = Õ ( C 2/γL∆ 2+2/γ + Φ 2/γL∆ 2 + C 2/γ 2/γ + Φ2/γ ) .
Proof. To satisfy the requirements of Lemma 1 in the constant step size, constant m case, we need
3αm1− γ/2Φ ≤ mα ≤ 1
6L .
This breaks apart into first a constraint on only m:
3Φ ≤ mγ/2, (14) and then a constraint on α in terms of m:
α ≤ 1 6Lm .
In the nonconvex setting, we invoke Lemma 1 followed by summing up over K phases and telescoping,
αm
4 K−1∑ k=0 ‖∇f(wmk)‖2 ≤ f(w0)− f∗ + 2K αm α2m2−γC2
= f(w0)− f∗ + 2Kαm1−γC2, which gives a rate of
1
K K−1∑ k=0 ‖∇f(wmk)‖2 ≤ 4(f(w0)− f∗) αmK + 8m−γC2.
If we want to minimize the right side, observe that we will always want to set α as large as possible, i.e. set α = 1/6Lm. This gives
1
K K−1∑ k=0 ‖∇f(wmk)‖2 ≤ 24L(f(w0)− f∗) K + 8m−γC2.
Now, for this to all be less than 2, it suffices for
24L(f(w0)− f∗) K ≤ 2 2 and 8m−γC2 ≤ 2 2 .
The former occurs when
K ≥ 48L(f(w0)− f ∗)
2 ,
while the latter happens when
4C
≤ mγ/2. (15)
If we let ∆ = f(w0)− f∗, using Eqs. (14) and (15) and taking the minimum m and K required, we can bound the number of iterations as
T = mK ≥ ⌈ 48L∆
2
⌉ · ⌈( 4C + 3Φ )2/γ⌉ .
So,
T = O ( C2/γL∆
2+2/γ + L∆Φ2/γ 2 + C2/γ 2/γ + Φ 2/γ
) ,
where we suppose that γ acts as a constant.

heading: A.4.2 CONSTANT STEP SIZE: STRONGLY-CONVEX CASE
text: Theorem 2. Suppose that f satisfies the µ-PL condition and our setup satisfies Assumptions 1 and 2. Let > 0 be any target error, and κ = L/µ be the condition number of the problem. Using SGD (1) with a constant step size α = 16L ⌈ (8C2/(µ 2) + 9Φ2)1/γ ⌉−1, the number of steps T needed to guarantee f(wT )− f∗ ≤ 2 is at most
T = ⌈ 12κ log ( 2∆ 2 )⌉ · ⌈( 8C2 µ 2 + 9Φ 2 )1/γ⌉ = Õ ( C2/γκ µ1/γ 2/γ + κΦ2/γ + κ ) .
Proof. The strongly convex setting has the same constraints on m and α as the non-convex setting. Recall that µ-strong convexity or µ-PL of f implies that for all x, ‖∇f(x)‖2 ≥ 2µ(f(x) − f∗). Applying this to the result of Lemma 1 gives
f(wτ+m) ≤ f (wτ )− αmµ
2 (f(wτ )− f∗) + 2αm1−γC2,
which is equivalent to f(wτ+m)− f∗ ≤ (
1− αmµ 2
) (f(wτ )− f∗) + 2αm1−γC2.
Applying this recursively over K intervals of length m starting from τ = 0,
f(wmK)− f∗ ≤ (
1− αmµ 2
)K (f(w0)− f∗) + 2αm1−γC2 K−1∑ k=0 ( 1− αmµ 2 )k ≤ exp ( −αmµK
2
) (f(w0)− f∗) + 2αm1−γC2 ∞∑ k=0 ( 1− αmµ 2 )k = exp ( −αmµK
2
) (f(w0)− f∗) + 2αm1−γC2 ( 1− ( 1− αmµ
2 ))−1 = exp ( −αmµK
2
) (f(w0)− f∗) + 2αm1−γC2 · 2
αmµ
= exp ( −αmµK
2
) (f(w0)− f∗) + 4C2
µ m−γ .
As in the non-convex case, it best advantages us to set α as large as possible. Setting α = 1/6Lm gives f(wmK)− f∗ ≤ exp ( − µK
12L
) (f(w0)− f∗) + 4C2
µ m−γ .
For this to be less than 2, it suffices for each of these two terms to be less than 2/2. In the latter case, we would need
8C2 µ 2 ≤ mγ . (16)
The former case holds when
K ≥ 12L µ log
( 2∆
2
) = 12κ log ( 2∆
2
) ,
where κ = L/µ. Applying a ceiling, using the minimum requirement on K and m from Eq. (16) above and Eq. (14) required in Lemma 1, we can lower bound the total number of iterations as
T = mK ≥ ⌈ 12κ log ( 2∆
2
)⌉ · ⌈( 8C2
µ 2 + 9Φ2
)1/γ⌉ .
Ignoring logarithmic terms gives
T = Õ ( κ ( C2
µ 2
)1/γ + κΦ 2/γ ) .

heading: A.5 JUSTIFICATIONS FOR ASSUMPTION 2 UNDER VARIOUS EXAMPLE ORDERINGS

heading: A.5.1 PROOF FOR PROPOSITION 1 (ARBITRARY PERMUTATION)
text: Proposition 1. Let Assumptions 1 and 3 hold, and suppose that we are using a permutation-based method of sampling, that is, any method such that xkn, xkn+1, . . . , xkn+n−1 is a permutation of x(0), x(1), . . . , x(n−1) for all epochs k ≥ 0. Any such method satisfies Assumption 2 with γ = 2, C2 = n2A2 and Φ2 = n2B2.
Proof. This result follows trivially from Assumption 3. Consider an arbitrary sum of gradient errors going from τ to τ +m− 1:
τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ),
where σ(t) is the index into the training set given by the permutation σ used at time step t. Since the interval {τ, τ + 1, . . . , τ +m− 1} is arbitrary, this m can potentially be greater than n, the epoch size. We can split this interval up as follows. Let τ1 be the first epoch boundary in the interval, such that all t going from τ to τ1 − 1 are within the same epoch as wτ , or else τ1 = τ +m if there is no epoch boundary in the interval. Let τ2 be the last epoch boundary in the interval, such that all t going from τ2 to τ +m−1 are within a later epoch than wτ (it may be the case that τ1 = τ2). Then the sum
τ2−1∑ t=τ1 ∇f(wτ ;xσ(t))−∇f(τ)
must be zero, since this interval goes over full epochs. This leaves us with∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤
∥∥∥∥∥ τ1−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥+ ∥∥∥∥∥ τ+m−1∑ t=τ2 ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤ n ( A2 +B2‖∇f(wτ )‖2 )1/2 .
where we have used the fact that the intervals {τ, . . . , τ1} and {τ2, τ +m− 1} are of length at most bn/2c (if such interval is of length greater than bn/2c, then bounding its norm is equivalent to bounding the sums of the remaining terms, for which there are at most bn/2c of them). Therefore∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2 ≤ 1 m2 n2 ( A2 +B2‖∇f(wτ )‖2 ) ,
and so Assumption 2 is satisfied with γ = 2, C2 = n2A2, and Φ2 = n2B2.

heading: A.5.2 PROOF FOR PROPOSITION 2 (SHUFFLE ONCE)
text: Proposition 2. Suppose that we are using the shuffle once variant of SGD to learn over a region B ∈ Rd of radius at most R, such that the iterates wt are guaranteed to remain within this region. Assume that for all w ∈ B and all examples x in the training set of size n, Assumption 1 (LSmoothness) and Assumption 3 hold. Then with probability at least 1− p, Assumption 2 holds with γ = 2, C2 = Õ(dA2(n+B2)), and Φ2 = Õ(ndB2).3
Proof. We begin by invoking Lemma 6 (Permuted vector Hoeffding inequality) with
Xi,j = ∇f(w;xi,j)−∇f(w)( A2 +B2‖∇f(w)‖2 )1/2 3The probability p is only present in log terms in these expressions, so it does not appear in the Õ.
where we use xi,j to denote the j-th example after time step i, for all i ∈ {1, . . .m}, j ∈ {1, . . . , n}, with n ≥ m > 0. Note that here we are using a fixed w ∈ B that does not depend on the permutation σ used in shuffle once. Clearly, ∑n j=1Xi,j = 0 for all i due to the periodicity in shuffle once, and ‖Xi,j‖ ≤ 1 by our assumption. Therefore the requirements of Lemma 6 are satisfied, which gives the following high probability bound for some > 0 on any sequence drawn without replacement:
P ∥∥∥∥∥ m∑ t=1 ∇f(w;xσ(t))−∇f(w) ∥∥∥∥∥ 2 ≥ 2 ( A2 +B2‖∇f(w)‖2 ) ≤ 2e2 exp(− 2 32m ) .
Using the periodicity property of shuffle once, i.e. if σ is the drawn permutation at the beginning of training then σ(t+ n) = σ(t), which allows us to arbitrarily shift the starting time step in the above:
P ∥∥∥∥∥ τ+m−1∑ t=τ ∇f(w;xσ(t))−∇f(w) ∥∥∥∥∥ 2 ≥ 2 ( A2 +B2‖∇f(w)‖2 ) ≤ 2e2 exp(− 2 32m ) .
Consider an arbitrary interval of length m out of the permutation. As all of the n(n−1)/2 intervals are equivalent to an interval of size at most bn/2c, by the sum to 0 property, we only need to consider intervals with m ≤ n/2. This gives us
P ∥∥∥∥∥ τ+m−1∑ t=τ ∇f(w;xσ(t))−∇f(w) ∥∥∥∥∥ 2 ≥ 2 ( A2 +B2‖∇f(w)‖2 ) ≤ 2e2 exp(− 2 16n ) .
By a union bound across all intervals we have
P ∃τ,m ∈ Z, ∥∥∥∥∥ τ+m−1∑ t=τ ∇f(w;xσ(t))−∇f(w) ∥∥∥∥∥ 2 ≥ 2 ( A2 +B2‖∇f(w)‖2 ) ≤ e2n2 exp ( − 2
16n
) .
Unfortunately, we cannot take a union bound over Rn or B (or even just over {w0, . . . , wT }), since the set may contain iterates visited by shuffle once that depend on σ, and thus Lemma 6 may not hold. Instead, suppose that we cover the whole region B with balls of radius δ. We will be able to do this with (1 + 2R/δ)d balls by Lemma 7 ( -net lemma). So, if the centers of the balls form a setW , then
P ∃τ,m ∈ Z, ŵ ∈ W, ∥∥∥∥∥ τ+m−1∑ t=τ ∇f(ŵ;xσ(t))−∇f(ŵ) ∥∥∥∥∥ 2 ≥ 2 ( A2 +B2‖∇f(ŵ)‖2 ) ≤ e2n2 ( 1 + 2R
δ
)d exp ( − 2
16n
) .
Note that the centers we use to cover the region is completely independent of the running algorithm, therefore the union bound over Lemma 6 can be applied. Next, consider some wτ that is not necessarily the center of a ball. The function
τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ )
is Ln-Lipschitz continuous, because each of the components is L-Lipschitz, and we can sum up only at most n/2 of them. Let ŵ ∈ W be such that ‖wτ − ŵ‖ ≤ δ. Adding and subtracting ∇f(ŵ;σ(i)) and ∇f(ŵ), applying triangle inequality followed by Lipschitz continuity, we have∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤ 2Lnδ + ∥∥∥∥∥ τ+m−1∑ t=τ ∇f(ŵ;xσ(t))−∇f(ŵ)
∥∥∥∥∥ ≤ 2Lnδ + √ A2 +B2‖∇f(ŵ)‖2. (w.h.p.)
Adding and subtracting∇f(wτ ) from ∇f(ŵ) and bound using δ and Lipschitz continuity again,∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤ 2Lnδ + √ A2 +B2‖∇f(ŵ −∇f(wτ ) +∇f(wτ ))‖2
≤ 2Lnδ + √ A2 + 2B2L2δ2 + 2B2‖∇f(wτ )‖2,
and consequently∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2 ≤ 4L2n2δ2 + 2 2 ( A2 + 2B2L2δ2 + 2B2‖∇f(wτ )‖2 ) .
Note that this inequality will fail to hold with probability at most
e2n2 ( 1 + 2R
δ
)d exp ( − 2
16n
) .
If we want this to fail with probability less than some p, it suffices to set
2 = 16n ( log ( e2n2
p
) + d log ( 1 + 2R
δ
)) .
This gives∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2 ≤ 4L2n2δ2
+ 32n · ( log ( e2n2
p
) + d log ( 1 + 2R
δ
)) · ( A2 + 2B2L2δ2 + 2B2‖∇f(wτ )‖2 ) .
If we set δ such that L2nδ2 = A2, then we can bound this with∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2 ≤ 4nA2
+ 32n · ( log ( e2n2
p
) + d log ( 1 + 2RL √ n
A
)) · ( A2 + 2A2B2
n + 2B2‖∇f(wτ )‖2
) ,
which gives∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2
≤ 32n m2 · ( log ( e3n2 p ) + d log ( 1 + 2RL √ n A )) · ( A2 + 2A2B2 n + 2B2‖∇f(wτ )‖2 ) .
Hiding logarithmic terms involving n, p,R, L,A, the shuffle once setup satisfies the requirement of Assumption 2 that with probability at least 1− p,∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2 = Õ ( nd m2 ( A2 + 2A2B2 n + 2B2‖∇f(wτ )‖2
)) = Õ ( 1
m2
( ndA2 + 2dA2B2 + 2ndB2‖∇f(wτ )‖2 )) = Õ ( 1
m2
( dA2(n+ 2B2) + 2ndB2‖∇f(wτ )‖2 )) ,
and so the parameters are γ = 2, C2 = Õ(dA2(n+B2)), and Φ2 = Õ(ndB2).

heading: A.5.3 PROOF FOR PROPOSITION 3 (RANDOM RESHUFFLING)
text: We first introduce a lemma that bounds the distance between iterates obtained from SGD with random permutation and that obtained from deterministic gradient descent, over one epoch. Lemma 3. Consider a single epoch of SGD with constant step size α > 0, where the examples come from a random permutation over the training set. Let Assumption 1 (L-Smoothness) and the bounded gradient error Assumption 3 be satisfied. Without loss of generality, assume that the epoch starts at time t = 0 at w0. Let the sequence ut be defined by u0 = w0 and
ut+1 = ut − α∇f(ut),
while
wt+1 = wt − α∇f(wt;xt)
for some xt chosen from the permutation. With probability at least (1− δ) it will hold that for all T ∈ {1, . . . , n}, if we set
α ≤ min
{ 1
64e(e2 + 1)BnL log (
2e2n δ
) , 1 2nL
} ,
then
‖wT − uT ‖2 ≤ 128α2n(e2 + 1)2 · log ( 2e2n
δ
)( A2 +B2e2‖∇f(wT )‖2 ) .
Proof. Given our setup, observe that we can write the difference between these sequences as
wt+1 − ut+1 = wt − ut − α (∇f(wt;xt)−∇f(ut)) = wt − ut − α (∇f(wt;xt)−∇f(ut;xt))− α (∇f(ut;xt)−∇f(ut)) ,
such that summing this up and using w0 = u0 gives
wT − uT = −α T−1∑ t=0 (∇f(wt;xt)−∇f(ut;xt))− α T−1∑ t=0 (∇f(ut;xt)−∇f(ut)) ,
and so
‖wT − uT ‖ ≤ α T−1∑ t=0 ‖∇f(wt;xt)−∇f(ut;xt)‖+ α ∥∥∥∥∥ T−1∑ t=0 (∇f(ut;xt)−∇f(ut)) ∥∥∥∥∥. ≤ αL
T−1∑ t=0 ‖wt − ut‖+ α ∥∥∥∥∥ T−1∑ t=0 (∇f(ut;xt)−∇f(ut)) ∥∥∥∥∥. Recall the update rule of ut, with Assumption 1 (L-Smoothness), we obtain
‖∇f(ut+1)−∇f(ut)‖ ≤ αL‖∇f(ut)‖.
By the reverse triangle inequality, this implies
‖∇f(ut)‖ ≤(1− αL)−1‖∇f(ut+1)‖ ‖∇f(ut+1)‖ ≤(1 + αL)‖∇f(ut)‖,
which further implies for any T ∈ {1, · · · , n} and t ∈ {0, · · · , T − 1},
‖∇f(ut)‖ ≤(1− αL)−n‖∇f(uT )‖
= ( 1 + αL
1− αL
)n ‖∇f(uT )‖
≤ exp ( αLn
1− αL
) ‖∇f(uT )‖
≤ exp (2αLn) ‖∇f(uT )‖ ≤e‖∇f(uT )‖,
while for any t ∈ {T + 1, · · · , n}, ‖∇f(ut)‖ ≤(1 + αL)n‖∇f(uT )‖
≤e‖∇f(uT )‖, where we apply the condition that α ≤ 1/(2nL). For a given T , consider the following vector set
Xj = ∇f(uj ;xj)−∇f(uj)√ A2 +B2e2‖∇f(uT )‖ ,∃j ∈ {0, · · · , n},
it can be easily verified that they sum to zero.
Now we apply Lemma 6 (Permuted vector Hoeffding inequality) , for any γ ≥ 0,
P (∥∥∥∥∥ T−1∑ t=0 (∇f(ut;xt)−∇f(ut)) ∥∥∥∥∥ ≥ γ√A2 +B2e2‖∇f(uT )‖2 ) ≤ 2e2 exp ( − γ 2 32T ) ≤ 2e2 exp ( − γ 2
32n ) as T ≤ n. Now, this holds for just one T . By a union bound,
P ( ∃T ∈ {1, . . . , n}, ∥∥∥∥∥ T−1∑ t=0 (∇f(ut;xt)−∇f(ut)) ∥∥∥∥∥ ≥ γQT ) ≤ 2e2n exp ( − γ 2 32n ) .
where QT = √ A2 +B2e2‖∇f(uT )‖2. Now, if we set γ such that
2e2n exp ( − γ 2
32n
) = δ ⇒ γ2 = 32n log ( 2e2n
δ
) ,
then we get that
P ( ∃T ∈ {1, . . . , n}, ∥∥∥∥∥ T−1∑ t=0 (∇f(ut;xt)−∇f(ut)) ∥∥∥∥∥ ≥ γQT ) ≤ δ.
In this case, we will have that with probability at least (1− δ),
‖wT − uT ‖ ≤ αL T−1∑ t=0 ‖wt − ut‖+ αγQT .
If we let ρ0 = 0 and
ρT = αL T−1∑ t=0 ρt + αγQT ,
then ‖wT − uT ‖ ≤ ρT . Here, if T > 0, ρT+1 − ρT = αLρT + αγ(QT+1 −QT ),
on the other hand, obviously ρ1 = αγQ1,
ρT =αγ T−1∑ k=0 (1 + αL)T−k−1(Qk+1 −Qk)
=αγ ( T∑ k=1 (1 + αL)T−kQk − T−1∑ k=0 (1 + αL)T−k−1Qk )
≤αγ ( (1 + αL)
T−1∑ k=0 (1 + αL)T−k−1Qk − T−1∑ k=0 (1 + αL)T−k−1Qk +QT
)
≤αγ ( αLe1/2
T−1∑ k=0 Qk +QT ) ≤αγ(e2 + 1)QT
=αγ(e2 + 1) √ A2 +B2e2‖∇f(wT )‖2
≤αγ(e2 + 1)(A+Be‖∇f(wT )‖).
Put it back we obtain
‖wT − uT ‖ ≤αγ(e2 + 1) (A+Be‖∇f(uT )‖) ≤αγ(e2 + 1) (A+BeL‖wT − uT ‖+Be‖∇f(wT )‖) ,
which gives
(1− αγ(e2 + 1)BeL)‖wT − uT ‖ ≤ αγ(e2 + 1) (A+Be‖∇f(wT )‖) . If we require
α ≤ 1 64e(e2 + 1)BnL log ( 2e2n δ ) Squaring and substituting the value of γ gives
‖wT − uT ‖2 ≤ 128α2n(e2 + 1)2 · log ( 2e2n
δ
)( A2 +B2e2‖∇f(wT )‖2 ) with probability at least (1− δ).
We now provide justifications to Assumption 2 for the random reshuffling scheme. Proposition 3. Suppose that we are using the random reshuffling variant of SGD. Assume that for all w ∈ Rd and all examples, Assumption 1 (L-Smoothness) and Assumption 3 hold. For some p ∈ (0, 1), set the constant step size to satisfy α ≤ ( max { 1460BnL · log ( 4e2T/p ) , 2nL })−1 . Then with
probability at least 1− p, Assumption 2 holds with γ = 2, C2 = Õ(nA2) and Φ2 = Õ(nB2).
Proof. For some γ > 0 (but different from the γ of Lemma 3, consider the event that for some specific epoch k, for some wτ ∈ {w0, w1, . . . , wn(k−1)} ∪ {u0, u1, . . . , un}, where the ui are the u from Lemma 3 for epoch k, and for some τ and mk such that n(k − 1) ≤ τ < τ +mk ≤ nk,∥∥∥∥∥ τ+mk−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤ γ √ A2 +B2‖∇f(wτ )‖2.
Since all the xσ(t) are from the k-th epoch, but wτ is independent of any randomness in the k-th epoch (as it is either a point visited in a previous epoch, or a value from the u sequence which depends only on the position at the start of the k-th epoch and not on any k-th epoch randomness), it follows that we can apply Lemma 6 (Permuted vector Hoeffding inequality) on either this sum or, alternatively, the terms from epoch k but not in the sum (the terms left out) to get that
P (∥∥∥∥∥ τ+mk−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≥ γ √ A2 +B2‖∇f(wτ )‖2 )
≤ 2e2 exp ( − γ 2
32 min(mk, n−mk)
) .
As mk ≤ n, it follows that
P (∥∥∥∥∥ τ+mk−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≥ γ √ A2 +B2‖∇f(wτ )‖2 )
≤ 2e2 exp ( − γ 2
16n
) .
Now, by a union bound the probability that there exists some τ , wτ , and mk such that the average gradient error is large is bounded by
P ( ∃ ∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(w) ∥∥∥∥∥ ≥ γ √ A2 +B2‖∇f(wτ )‖2 )
≤ 2e2nT 2 exp ( − γ 2
16n
) ,
where T is the total number of iterations across all epochs (and we assume that we finish all epochs so n divides T ). This follows from the fact that there are at most T such τ that we could take on, at most T values of wτ that can be adopted for each, and at most n values mk can take on. If we set γ such that
2e2nT 2 exp ( − γ 2
16n
) = p
2 ⇒ γ2 = 16n log
( 4e2nT 2
p
) , (17)
then with probability at least (1− p/2), it will hold that for all epochs, wτ , τ , and mk,∥∥∥∥∥ τ+mk−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤ γ √ A2 +B2‖∇f(wτ )‖2. (18)
Additionally, set the δ in Lemma 3 to be pn/2T . By a union bound on the result of Lemma 3, across all T/n epochs, with probability at least (1− p/2), it must follow that
‖wT − uT ‖2 ≤ 128α2n(e2 + 1)2 · log ( 4e2T
p
)( A2 +B2e2‖∇f(wT )‖2 ) , (19)
for the corresponding ut sequence for all epochs. Therefore, both of these inequalities Eqs. (18) and (19) hold together with probability at least (1− p). Now, consider an arbitrary sum of gradient errors going from τ to τ +m− 1. Note since the interval here is arbitrary, this m can be different from mk and potentially greater than n. We can split this interval up as follows. Let τ1 be the first epoch boundary in the interval, such that all t going from τ to τ1 − 1 are within the same epoch as wτ , or else τ1 = τ +m if there is no epoch boundary in the interval. Let τ2 be the last epoch boundary in the interval, such that all t going from τ2 to τ +m− 1 are within a later epoch than wτ (it may be the case that τ1 = τ2). It follows that∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ )
∥∥∥∥∥ ≤
∥∥∥∥∥ τ1−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ +
∥∥∥∥∥ τ2−1∑ t=τ1 ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥+ ∥∥∥∥∥ τ+m−1∑ t=τ2 ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥. Observe that since the second of these sums must go over some number of full epochs, its value must be 0. Therefore,∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ )
∥∥∥∥∥ ≤
∥∥∥∥∥ τ1−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥+ ∥∥∥∥∥ τ+m−1∑ t=τ2 ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥. Observe that the term in the first sum must be nL/2-Lipschitz continuous in w, because it can be written as the sum of at most bn/2c terms each of which is L-Lipschitz (either as the actual terms of the sum, or else the terms left out of the sum). So, add and subtract∇f(uτ ;xσ(t)) and ∇f(uτ ),∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ )
∥∥∥∥∥ ≤ nL‖wτ − uτ‖
+ ∥∥∥∥∥ τ1−1∑ t=τ ∇f(uτ ;xσ(t))−∇f(uτ ) ∥∥∥∥∥+ ∥∥∥∥∥ τ+m−1∑ t=τ2 ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥.
Now applying Eq. (19) on the first term and Eq. (18) on the last two terms gives∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ ≤ nL · √ 128α2n(e2 + 1)2 · log ( 4e2T
p
)( A2 +B2e2‖∇f(wτ )‖2 ) + 2γ √ A2 +B2‖∇f(wτ )‖2.
Squaring both sides for simplicity, we can bound this with∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2
≤ 256α2n3L2(e2 + 1)2 · log ( 4e2T
p
)( A2 +B2e2‖∇f(wτ )‖2 ) + 8γ2(A2 +B2‖∇f(wτ )‖2),
where we have also used Eq. (17) for γ. If we apply our requirement that αLn ≤ 1/2, we get∥∥∥∥∥ τ+m−1∑ t=τ ∇f(wτ ;xσ(t))−∇f(wτ ) ∥∥∥∥∥ 2 ≤ 64n(e2 + 1)2 · log ( 4e2T p )( A2 +B2e2‖∇f(wτ )‖2 ) + 128n(A2 +B2‖∇f(wτ )‖2) · log ( 4e2nT 2
p
) .
It follows that for random reshuffling with probability 1− p, if we set
α ≤min  164(e2 + 1)BenL log ( 4e2Tp ) , 1 2nL  = [ max { 1460BnL log ( 4e2T
p
) , 2nL }]−1 ,
random reshuffling satisfies the requirement of Assumption 2 with γ = 2, C2 = Õ(nA2), and Φ2 = Õ(nB2), where Õ(·) hides logarithmic terms in n, T, p.

heading: A.5.4 PROOF FOR PROPOSITION 4 (RANDOM RESHUFFLING WITH DATA ECHOING)
text: The proof for Proposition 4 is nearly a repeat of the random reshuffling proof in Proposition 3. The trick is to re-define an epoch when data echoing is used. Instead of referring to an epoch as a random permutation of the n examples in the training set as in vanilla random reshuffling, here we define the cn samples as an epoch, where each example σ(i) is repeated c times. For instance, let c = 3, then the k-th epoch with permutation σk is given by the sequence
. . . , xσk(1), xσk(1), xσk(1), xσk(2), . . . , xσk(n−1), xσk(n), xσk(n), xσk(n)︸ ︷︷ ︸ examples used in epoch k , . . . .
As a result of this redefinition, Lemma 3 needs to be modified such that wherever n appears, we now need cn. We omit the proof for the following lemma as it is a trivial adaptation from that of Lemma 3.
Lemma 4. Consider a single epoch of SGD with constant step size α > 0, where the examples come from a random permutation over the training set, each echoed c times. Without loss of generality, assume that the epoch starts at time t = 0 at w0. Let the sequence ut be defined by u0 = w0 and
ut+1 = ut − α∇f(ut),
while
wt+1 = wt − α∇f(wt;xt)
for some xt chosen from the permutation. Under the same assumptions as Lemma 3, with probability at least (1− δ) it will hold that for all T ∈ {1, . . . , cn},
‖wT − uT ‖2 ≤ 128α2cn(e2 + 1)2 · log ( 2e2cn
δ
)( A2 +B2e2‖∇f(wT )‖2 ) .
The justifications for Assumption 2 for the random reshuffling with data echoing scheme can also be obtained from Appendix A.5.3 similarly. At appropriate places one should invoke Lemma 4 instead of Lemma 3, and replace n with cn. We omit the proof as it is again a trivial modification.

heading: A.5.5 PROOF FOR PROPOSITION 5 (MARKOV CHAIN GRADIENT DESCENT (MCGD))
text: To justify Assumption 2 for MCGD, we will first need the following lemma. Lemma 5. Let F be any vector-valued measurable function, and let x0, x1, . . . be a sequence of samples from a Markov chain with mixing time tmix and stationary distribution π starting from an arbitrary initial distribution. If the function is constrained such that ‖F (x)‖ ≤ 1 for all x, and if we also have EX∼π[F (X)] = 0, then for any δ ∈ (0, 1),
P (∥∥∥∥∥ m−1∑ i=0 F (xi) ∥∥∥∥∥ ≥ 5tmix √ 2m log ( 2e2 δ )) ≤ δ,
Proof. Consider the Doob martingale
Wk = E [ m−1∑ i=0 F (xi) | Fk ] ,
where Fk contains all randomness up to timestep k, i.e. x0, x1, . . . , xk, and so (as usual for a Doob martingale) the martingale property is trivially satisfied using repeated conditioning:
E [Wk+1 | Fk+1] = E [ E [ m−1∑ i=0 F (xi) | Fk+1 ] | Fk ] = E [ m−1∑ i=0 F (xi) | Fk ] = Wk.
Observe that the sum we want is Wm = E [∑m−1 i=0 F (xi) | Fm ] = ∑m−1 i=0 F (xi), and that this sum has increments
Wk+1 −Wk = E [ m−1∑ i=0 F (xi) | Fk+1 ] −E [ m−1∑ i=0 F (xi) | Fk ]
= E [ k+1∑ i=0 F (xi) | Fk+1 ] + E [ m−1∑ i=k+2 F (xi) | Fk+1 ]
−E [ k∑ i=0 F (xi) | Fk ] −E [ m−1∑ i=k+1 F (xi) | Fk ]
= k+1∑ i=0 F (xi)− k∑ i=0 F (xi) + m−1∑ i=k+2 E [F (xi) | Fk+1]− m−1∑ i=k+1 E [F (xi) | Fk]
= F (xk+1) + m−1∑ i=k+2 E [F (xi) | Fk+1]− m−1∑ i=k+1 E [F (xi) | Fk] .
Now, observe that since the mixing time of the Markov chain is tmix, for any i ≥ k, if µ denotes the distribution of xi conditioned on Fk, then using results from Levin & Peres (2017, Section 4.5)
‖µ− π‖TV ≤ 2 −b(i−k)/tmixc.
It follows from this and the fact that F is bounded that
‖E [F (xi) | Fk]‖ ≤ 2−b(i−k)/tmixc.
Therefore, ∥∥∥∥∥ m−1∑ i=k+2 E [F (xi) | Fk+1] ∥∥∥∥∥ ≤ m−1∑ i=k+2 ‖E [F (xi) | Fk+1]‖
≤ m−1∑ i=k+2 2−b(i−k−1)/tmixc
≤ ∞∑
i=k+2
2−b(i−k−1)/tmixc = 2tmix − 1,
and similarly for the last term. It follows that
‖Wk+1 −Wk‖ ≤ 4tmix.
Therefore, by the vector Azuma’s inequality of Hayes (2005, Theorem 1.8), for any a > 0, P (‖Wm −W0‖ ≥ 4tmixa) ≤ 2e2 exp ( − a 2
2m
) .
On the other hand, by the same reasoning as before,
‖W0‖ = ∥∥∥∥∥E [ m−1∑ i=0 F (xi) | F0 ]∥∥∥∥∥ ≤ ‖F (x0)‖+
∥∥∥∥∥E [ m−1∑ i=1 F (xi) | F0 ]∥∥∥∥∥ ≤ 1 + 2tmix − 1 = 2tmix.
So,
P (‖Wm‖ ≥ 2tmix + 4tmixa) ≤ 2e2 exp ( − a 2
2m
) .
Now, setting a such that
2e2 exp ( − a 2
2m
) = δ ⇒ a2 = 2m log ( 2e2
δ
) ,
and noting that this makes a ≥ 2, we get
P (∥∥∥∥∥ m−1∑ i=0 F (xi) ∥∥∥∥∥ ≥ 5tmix √ 2m log ( 2e2 δ )) ≤ δ,
which is what we wanted to show.
The justification for Assumption 2 then follows straightforwardly. Proposition 5. Suppose that we use samples xt from a Markov chain with mixing time tmix. Assume that for all w ∈ R and all examples xt, Assumption 3 holds. Then with probability at least 1 − p, Assumption 2 holds with γ = 1, C2 = Õ(A2t2mix), and Φ2 = Õ(B2t2mix).
Proof. Observe that we need Lemma 5 to hold for all subintervals of examples, of which there are only at most T 2. So, setting δ to be p/T 2 in Lemma 5, we can show that for all τ,m, the probability∥∥∥∥∥ τ+m−1∑ t=τ (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 ≥ 25 ( A2 +B2‖∇f(wτ )‖2 ) t2mix2m log ( 2e2T 2 p )
is at least 1− p. Equivalently,∥∥∥∥∥ 1m τ+m−1∑ t=τ (∇f(wτ ;xt)−∇f(wτ )) ∥∥∥∥∥ 2 ≤ 50 m t2mix ( A2 +B2‖∇f(wτ )‖2 ) · log ( 2e2T 2 p ) .
It follows that MCGD satisfies the requirements of Assumption 2 with γ = 1,
C2 = 50A2t2mix log
( 2e2T 2
p
) and Φ2 = 50B2t2mix log ( 2e2T 2
p
) .
The big-Õ expressions in the proposition statement immediately follow.

heading: A.5.6 JUSTIFICATION FOR QMC-BASED DATA AUGMENTATION WITH RANDOM RESHUFFLING
text: Before we begin with the proof, let us first present the introductory material necessary on quasi-Monte Carlo methods. Similar to Monte Carlo integration, QMC is also used for numerical integration but using low-discrepancy sequences instead of pseudorandom number sequences. Concretely, the problem is to approximate the integral of a function f over some s-dimensional hypercube [0, 1]s,
I(f) := ∫ [0,1]s f(x)dx
using the average of the function evaluated at a sequence of points x1, . . . , xm,
Im(f) := 1
m m∑ i=1 f(xi).
The approximation error rate is defined as = |I(f)− Im(f)|, and it is well-known that in the case of Monte Carlo integration where the xi’s are drawn uniformly at random from [0, 1]s, the error rate is 2 = O(1/m). If the sequence of xi’s has low star-discrepancy, which is defined as
D∗m := sup U ∣∣∣∣∣ 1m m∑ i=1 1{xi ∈ U} − Volume(U) ∣∣∣∣∣ where U = ∏s j=1[0, bi] for bi ∈ [0, 1), and the volume is measured using the s-dimensional Lebesgue measure. Intuitively, the smaller this quantity is the more evenly the sequence of points covers the space. Some popular low-discrepancy sequences include the Halton sequence, Sobol sequence, van der Corput sequence, etc., which are essentially deterministic sequences that are cleverly constructed to mimic random numbers but in fact have low star-discrepancy. For instance, the Halton sequence satisfies D∗m = O((logm) s /m).
To ensure fast convergence of Im to I as we increase m, in addition to using low discrepency sequences we also need f to be relatively well-behaved. For this, the Hardy-Krause variation VHK is often used, for which we refer the reader to Drmota & Tichy (2006, Definition 1.13) for its detailed characterization. Most importantly, if f has finite VHK on [0, 1]s, then the Koksma-Hlawka inequality guarantees that
≤ VHKD∗m. This implies that using quasi-Monte Carlo with, for instance, the Halton sequence, to integrate a function f with bounded Hardy-Krause variation, the error rate would be 2 = O((logm)2s/m2). In comparison to the Monte Carlo error rate, the QMC error rate can be much faster when the dimensionality s is relatively small. For an in-depth exposition of the related materials we recommend Owen (2003) and Drmota & Tichy (2006).
We are now ready to restate our proposition from Section 5, in which we describe our QMC data augmentation setup. Proposition 6. Suppose that we are using the random reshuffling variant of SGD with QMC data augmentation as described. Assume that for all w ∈ Rd and all examples, Assumptions 1 , 3, 4 and 5 hold for Equation 3. For some p ∈ (0, 1), set the step size to be a constant such that α ≤ ( max{1460BnL · log(4e2T/p ) , 2nL})−1. Then with probability at least 1− p, Assumption 2
holds with γ = 2, C2 = Õ(n2V 2C2QMC log(T )2s + nA2) and Φ2 = Õ(nB2).
Proof. Recall that the example used in the t-th iteration is being transformed as
xt = A(x(σbt/nc(t mod n)), ζbt/nc+σbt/nc(t mod n)).
We start from the average gradient error term, that is,∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f (wτ ;A (xt, ζt))−∇f(wτ ) ∥∥∥∥∥ 2
= ∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f (wτ ;A (xt, ζt))− 1 n n∑ i=1 E ζ∼[0,1]s ∇f(wτ ;A ( x(i); ζ ) ) ∥∥∥∥∥ 2
≤ 2 ∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f (wτ ;A (xt, ζt))− 1 m τ+m−1∑ t=τ E ζ∼[0,1]s ∇f(wτ ;A (xt; ζ)) ∥∥∥∥∥ 2
+ 2 ∥∥∥∥∥ 1m τ+m−1∑ t=τ E ζ∼[0,1]s ∇f(wτ ;A (xt; ζ))− 1 n n∑ i=1 E ζ∼[0,1]s ∇f(wτ ;A ( x(i); ζ ) ) ∥∥∥∥∥ 2 ,
where the first step follows from the definition of f in the data augmentation setup (Eq. (3)). In the second step, the first norm relates to the QMC approximation error while the second norm relates to the RR analysis. From the analysis in Section A.5.3 we know the second norm can be bounded asymptotically (and probabilistically) by Õ(nA2) + Õ(nB2)‖∇f(wτ )‖2/m2. Now we analyze the first term. Since we use a contiguous QMC subsequence on each individual example, for the period of [τ, τ +m− 1], we define τj and mj as the starting point and length of example x(j) being chosen during this period, such that
∑n j=1mj = m. With this notation, we can now rewrite the first norm as∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f (wτ ;A (xt, ζt))− E ζ∼[0,1]s ∇f(wτ ;A (xt; ζ)) ∥∥∥∥∥ 2
= ∥∥∥∥∥∥ 1m n∑ j=1 τj+mj−1∑ t=τj ∇f ( wτ ;A ( x(j), ζt )) − E ζ∼[0,1]s ∇f ( wτ ;A ( x(j); ζ ))∥∥∥∥∥∥ 2
≤ n m2 n∑ j=1 ∥∥∥∥∥∥ τj+mj−1∑ t=τj ∇f ( wτ ;A ( x(j), ζt )) − E ζ∼[0,1]s ∇f ( wτ ;A ( x(j); ζ ))∥∥∥∥∥∥ 2
≤ 2n(τj +mj) 2
m2
n∑ j=1 ∥∥∥∥∥ 1τj +mj τj+mj−1∑ t=0 ∇f ( wτ ;A ( x(j), ζt )) − E ζ∼[0,1]s ∇f ( wτ ;A ( x(j); ζ ))∥∥∥∥∥ 2
+ 2nτ2j m2 n∑ j=1 ∥∥∥∥∥ 1τj τj−1∑ t=0 ∇f ( wτ ;A ( x(j), ζt )) − E ζ∼[0,1]s ∇f ( wτ ;A ( x(j); ζ ))∥∥∥∥∥ 2 .
Now we can use the Koksma–Hlawka inequality (Aistleitner & Dick, 2014) to bound these norms, which gives us∥∥∥∥∥ 1m τ+m−1∑ t=τ ∇f (wτ ;A (xt, ζt))− E ζ∼[0,1]s ∇f(wτ ;A (xt; ζ)) ∥∥∥∥∥ 2
≤ O nC2QMCV 2 m2 n∑ j=1 · ( log(τj +mj) 2s + log(τj) 2s ) ≤ O(n2C2QMCV 2 log(T )2s m2 ) .
Putting it together, we have shown that random reshuffling with QMC data augmentation satisfies Assumption 2 with γ = 2, C2 = Õ ( n2C2QMCV 2 log(T )2s + nA2 ) , and Φ2 = Õ ( nB2 ) .

heading: A.6 MISCELLANEOUS LEMMAS
text: We now collect some technical lemmas used in our analyses.
The following lemma is a Hoeffding-type concentration bound on the sums of random permutations of vectors. This lemma is particularly useful for simplifying the logic of the proofs of our shuffling propositions above, as it frees us from having to use Doob-martingale-type arguments throughout.
Lemma 6 (Permuted vector Hoeffding inequality). Let n ≥ m > 0 be some integers, and let Xi,j for i ∈ {1, . . . ,m} and j ∈ {1, . . . , n} be vectors in Rd. We also require for all i, j, ‖Xi,j‖ ≤ 1, and that for all i,
n∑ j=1 Xi,j = 0.
Then, if σ : {1, . . . , n} → {1, . . . , n} is a random permutation and that the individual Xi,j’s do not depend on σ,
P (∥∥∥∥∥ m∑ i=1 Xi,σ(i) ∥∥∥∥∥ ≥ a ) ≤ 2e2 exp ( − a 2 32m ) .
Note that this lemma trivially also applies to random subsamples of vector sums, by just letting Xi,j = Yj for the desired vector sequence Yj .
Proof. This proof is adapted from the proof of Theorem 4.3 in Bercu et al. (2015). Their approach is more sophisticated and tends to a Bernstein-type inequality, but is specialized to only scalars. Consider sampling the elements of the permutation one at a time, and letFj be the filtration containing the random variables σ(1), σ(2), . . . , σ(j) but not σ(j + 1), . . . , σn. For k ≤ n, define the process
Wk = E [ m∑ i=1 Xi,σ(i) | Fk ] .
Observe that this must be a (vector) martingale process, as it is a Doob martingale. Explicitly, we can write it as
Wk = E [ k∑ i=1 Xi,σ(i) + 1 n− k m∑ i=k+1 ∑ j /∈{σ(1),...,σ(k)} Xi,j | Fk ]
= E  k∑ i=1 Xi,σ(i) − 1 n− k m∑ i=k+1 k∑ j=1 Xi,σ(j) | Fk  ,
where here we used the fact that theXi,j’s sum to 0 over j = [n]. Thus, the increments of this process will have
Wk −Wk−1 = E  k∑ i=1 Xi,σ(i) − 1 n− k m∑ i=k+1 k∑ j=1 Xi,σ(j) | Fk  −E
k−1∑ i=1 Xi,σ(i) − 1 n− k + 1 m∑ i=k k−1∑ j=1 Xi,σ(j) | Fk−1  = Xk,σ(k) − 1
n− k E  m∑ i=k+1 k∑ j=1 Xi,σ(j) | Fk  + 1
n− k + 1 E  m∑ i=k k−1∑ j=1 Xi,σ(j) | Fk−1  = Xk,σ(k) − 1
n− k m∑ i=k+1 Xi,σ(k)
− ( 1
n− k − 1 n− k + 1 ) m∑ i=k+1 k−1∑ j=1 Xi,σ(j) + 1 n− k + 1 k−1∑ j=1 Xk,σ(j),
where we have repeatedly applied the martingale property. Now to bound all of these, first by our assumption ∥∥Xk,σ(k)∥∥ ≤ 1. Also,
1
n− k ∥∥∥∥∥ m∑
i=k+1
Xi,σ(k) ∥∥∥∥∥ ≤ m− kn− k ≤ 1. Next, again using the sum to 0 assumption, we have
1
n− k + 1 ∥∥∥∥∥∥ k−1∑ j=1 Xk,σ(j) ∥∥∥∥∥∥ = 1n− k + 1 ∥∥∥∥∥∥− n∑ j=k Xk,σ(j) ∥∥∥∥∥∥ ≤ n− k + 1n− k + 1 = 1, and finally, by a combination of these bounds,(
1 n− k − 1 n− k + 1 )∥∥∥∥∥∥ m∑
i=k+1 k−1∑ j=1 Xi,σ(j) ∥∥∥∥∥∥ = 1(n− k)(n− k + 1) ∥∥∥∥∥∥ m∑ i=k+1 k−1∑ j=1 Xi,σ(j) ∥∥∥∥∥∥ ≤ 1
(n− k)(n− k + 1) · (m− k) · (n− k + 1)
≤ 1.
It follows that the increment satisfies
‖Wk −Wk−1‖ ≤ 4
with probability 1. We can now apply the vector Azuma’s inequality of Hayes (2005, Theorem 1.8). Applying this to Wm gives the result that, for any a,
P (‖Wm‖ ≥ a) ≤ 2e2 exp ( − a 2
32m
) .
This proves the lemma.
The next lemma we state is a standard result showing that we can cover a region of radius R with a bounded number of balls of radius .
Lemma 7 ( -net lemma). For any region D of radius R in d-dimensional space, and any > 0, there exists a subset S ⊆ D such that S is of size at most
|S| ≤ ( 1 + 2R )d ,
and for every point x ∈ D, there exists a x̂ ∈ S such that ‖x− x̂‖ ≤ .
Proof. Consider the following procedure. As long as there are points in D that are not within a distance of an existing point in S, choose one such point arbitrarily and add it to S. Observe that with this construction, any two points in S must be at a distance greater than from each other. This means that if we center a ball with radius /2 at each of the points in S, these balls will be disjoint.
The total volume of these balls will be |S| · V1 · ( /2)d, where V1 is the volume of the unit ball in d-dimensional space. However, the centers of all these balls must lie in D, and so the balls themselves must all lie within a slightly larger region of radius R+ /2. This region will have volume V1 · (R+ /2)d. This shows that our process must eventually stop, implying that S must exist, and gives us a bound on the size of S as
|S| ≤ (R+ /2) d
( /2) d
= ( 1 + 2R )d ,
and the proof is complete.

Review summaries:
1: This paper theoretically builds the following connection: if the averages of consecutive stochastic gradients converge faster to the full gradient, then the SGD with the corresponding sampling strategy will have a faster convergence rate.  2: Applying this theoretical finding on a few commonly used sampling strategies, the convergence rates established by prior works can be recovered.  3: Two algorithms are proposed based on this theory. 

The paper proposes a new condition (metric) on the sequence of samples seen by SGD (with or without replacement), which can be used to prove tight convergence rates for strongly convex and non-convex objectives. The metric quantifies how quickly the average gradient of a subsequence of length $m$ converges to the true gradient. This metric is computed for popular schemes like shuffle once, and random reshuffling, as well as for random reshuffling with echoing and Markov chain gradient descent. Then, the paper proposes to combine quasi-Monte Carlo based data augmentation with random reshuffling, proves the condition for it, and shows empirically that it performs well. The paper also proposes a technique which greedily tries to optimize the metric in the proposed condition.

The paper studies the dependency of SGD convergence on order of examples. For smooth loss functions, the paper bounds SGD convergence based on averages of consecutive example gradients. Specifically better the average of consecutive example gradients approximate the objective gradient, faster the SGD convergence is. The paper recovers (and in some cases improves) the convergence results of various example selection schemes.   Based on the framework. the paper also introduces couple of new example selection sequences: QMC based data augmentation and greedy algorithm that optimizes on  average gradient metric. The authors also propose using stale gradients and random projection for faster approximate greedy example selection.  The paper evaluates their theoretical guarantees on several deep learning benchmarks and demonstrate the faster convergence with respect to other example selection techniques such as random sampling, shuffle once and random reshuffling.