summaries: # === Update ===  I have decided to maintain my score. This is very strong submission in my opinion. A brief overview of my thoughts is the following:  **Topic**: This paper comes in the middle of a wave of interest in permutation SGD and similar variants. It is highly topical and a good fit for the conference.  **Theoretical Results**: It is easy focus on FlipFlop and its limited analysis (i.e. quadratics only), but there are many other highly interesting results in this submission. Theorems 1 and 2 show that one-dimensional functions are "easy" for permutation-based optimization, while quadratics in higher dimensions are hard. Yes, this result requires a fixed step-size, but it is also the first analysis of permutation SGD conducted on the level of specific permutations. Theorem 3 refines these results to show that the improvement in Theorem 1 is specific to quadratics. The proof technique for these results is elegant (continuity of the composition of updates + IVT) and, as a I noted in my review, leads to useful intermediate results.  **Experiments** The experiments are limited to synthetic data, but this is a theoretical paper with ~37 pages of appendices already. I think that an intensive evaluation of FlipFlop is beyond the scope of this submission. Moreover, since the analysis considers quadratics only, I think that the nice results on logistic regression in Appendix H are quite encouraging.  # ======  This submission analyzes the effects of permutation choice on the convergence of the random-reshuffling, shuffle-once, and incremental variants of stochastic gradient descent (SGD). In particular, the goal of this work is to determine settings whether specific permutations can give faster convergence than random reshuffling and in what settings.  The authors show that, in the case of one-dimensional finite-sum functions with smooth Hessians, there exists a sequence of permutations for which SGD converges at a linear (i.e. exponential) rate.  However, they also prove that this phenomenon is specific to one-dimensional functions using a new, dimension-dependent lower bound for any permutation-based SGD method. This lower bound is tight with the known convergence rate of random reshuffling, implying this method is optimal in the general setting. Finally, the authors restrict themselves to finite-sums of convex quadratics and develop a simple heuristic permutation "schedule" giving provably faster convergence.  This approach, which is called FlipFlop, alternates between fresh permutations and the previous permutation, but in _reverse_ order. Small experiments confirm that FlipFlop accelerates convergence of permutation-based SGD variants.  This paper is motivated by the observed phenomenon that, in stochastic gradient descent (SGD), without-replacement sampling (random permutation) gives faster convergence than with-replacement sampling. The paper studies whether random permutations are optimal among permutation-based SGD, by considering different deterministic, random, or hybrid ways of generating permutations of input points.  Focusing on optimizing convex functions at a constant step size, the paper shows that: 1. there exist optimal permutations which converges exponentially faster than random permutations for 1-dimensional functions. 2. such improvement is not possible in higher dimensions or for strongly convex objectives, where random permutations are optimal. 3. by reversing the permutation every other epoch (flipflopping), convergence on quadratic functions improves for three permutation-based methods: Incremental Gradient Descent (deterministic), Random Reshuffle (random), and Single Shuffle (hybrid). This work investigates the optimality of random permutation as a scan ordering for SGD. They found that for general strongly convex functions with Lipschitz Hessian, random permutations are optimal in high dimension but not optimal in 1-dimension. For general convex quadratics, random permutations are also not optimal. Finally, the authors introduced a new technique termed FlipFlop that works by reversing the permutation of the previous epoch at every even epoch, and at the odd epochs the algorithm just follows its original permutation, whether it be random or cyclic. FlipFlop has been proven to improve the convergence of random reshuffling, single shuffle, and incremental gradient descent on quadratic functions. Experiments on a 100-dimensional quadratic objective and 1-dimensional logistic regression demonstrate that FlipFlop indeed converges faster than random reshuffling. This paper studies theoretical properties of permutation-based fixed-step size SGD for finite sum optimization. The main theorems state that  * For 1-d Hessian-smooth functions, there exists permutations such that the convergence rate is exponential in the number of iterations $K$ * For higher dimension, there exists strongly-convex function such that the convergence rate is at best $O(1/K^3)$ * If some component in the strongly-convex objective function is nonconvex, the convergence rate has a lower bound $O(1/K^2)$  The paper then proposes the algorithm FlipFlop, where adjacent epochs use permutations of reverse order. The algorithm is proved to improve upon random permutation when all the component functions are quadratic. The results are corroborated with simulations.
abstractText: A recent line of ground-breaking results for permutation-based SGD has corroborated a widely observed phenomenon: random permutations offer faster convergence than with-replacement sampling. However, is random optimal? We show that this depends heavily on what functions we are optimizing, and the convergence gap between optimal and random permutations can vary from exponential to nonexistent. We first show that for 1-dimensional strongly convex functions, with smooth second derivatives, there exist permutations that offer exponentially faster convergence compared to random. However, for general strongly convex functions, random permutations are optimal. Finally, we show that for quadratic, strongly-convex functions, there are easy-to-construct permutations that lead to accelerated convergence compared to random. Our results suggest that a general convergence characterization of optimal permutations cannot capture the nuances of individual function classes, and can mistakenly indicate that one cannot do much better than random.
references.authors: Kwangjun Ahn Chulhee Yun Suvrit Sra
references.title: Sgd with shuffling: optimal rates without component convexity and large epoch requirements, 2020
references.authors: Dimitri P Bertsekas
references.title: Incremental gradient, subgradient, and proximal methods for convex optimization: A survey
references.venue: Optimization for Machine Learning,
references.authors: Léon Bottou
references.title: Curiously fast convergence of some stochastic gradient descent algorithms. Unpublished open problem offered to the attendance of the SLDS
references.venue: URL http: //leon.bottou.org/papers/bottou-slds-open-problem-2009
references.authors: Sébastien Bubeck
references.title: Convex optimization: Algorithms and complexity
references.venue: Foundations and Trends® in Machine Learning,
references.authors: Christopher M De Sa
references.title: Random reshuffling is not always better
references.venue: Advances in Neural Information Processing Systems,
references.authors: M Gürbüzbalaban Asuman Ozdaglar Pablo A Parrilo
references.title: Convergence rate of incremental gradient and incremental newton methods
references.venue: SIAM Journal on Optimization,
references.authors: Mert Gürbüzbalaban Asu Ozdaglar Pablo A Parrilo
references.title: Why random reshuffling beats stochastic gradient descent
references.venue: Mathematical Programming,
references.authors: Jeff Haochen Suvrit Sra
references.title: Random shuffling beats sgd after finite epochs
references.venue: In International Conference on Machine Learning,
references.authors: Zehua Lai Lek-Heng Lim
references.title: Recht-ré noncommutative arithmetic-geometric mean conjecture is false
references.venue: In International Conference on Machine Learning,
references.authors: Konstantin Mishchenko Ahmed Khaled Peter Richtárik
references.title: Random reshuffling: Simple analysis with vast improvements
references.venue: ArXiv, abs/2006.05988,
references.authors: Dheeraj Nagaraj Prateek Jain Praneeth Netrapalli
references.title: Sgd without replacement: Sharper rates for general smooth convex functions
references.venue: In International Conference on Machine Learning,
references.authors: Angelia Nedić Dimitri Bertsekas
references.title: Convergence rate of incremental subgradient algorithms. In Stochastic optimization: algorithms and applications, pp. 223–264
references.authors: Yurii Nesterov
references.title: Introductory Lectures on Convex Optimization - A Basic Course, volume
references.venue: Applied Optimization. Springer,
references.authors: Lam M Nguyen Quoc Tran-Dinh Dzung T Phan Phuong Ha Nguyen Marten van Dijk
references.title: A unified convergence analysis for shuffling-type gradient methods
references.venue: arXiv preprint arXiv:2002.08246,
references.authors: Shashank Rajput Anant Gupta Dimitris Papailiopoulos
references.title: Closing the convergence gap of sgd without replacement
references.venue: In International Conference on Machine Learning,
references.authors: Benjamin Recht Christopher Ré
references.title: Beneath the valley of the noncommutative arithmetic-geometric mean inequality: conjectures, case-studies, and consequences
references.venue: arXiv preprint arXiv:1202.4184,
references.authors: Benjamin Recht Christopher Ré
references.title: Parallel stochastic gradient algorithms for large-scale matrix completion
references.venue: Mathematical Programming Computation,
references.authors: Itay Safran Ohad Shamir
references.title: How good is sgd with random shuffling
references.authors: Itay Safran Ohad Shamir
references.title: Random shuffling beats sgd only after many epochs on ill-conditioned problems
references.venue: arXiv preprint arXiv:2106.06880,
references.authors: M. Schneider
references.title: Probability inequalities for kernel embeddings in sampling without replacement
references.venue: In AISTATS,
references.authors: Ohad Shamir
references.title: Without-replacement sampling for stochastic gradient methods
references.venue: In Proceedings of the 30th International Conference on Neural Information Processing Systems,
references.authors: Trang H Tran Lam M Nguyen Quoc Tran-Dinh
references.title: Shuffling gradient-based methods with momentum
references.venue: arXiv preprint arXiv:2011.11884,
references.authors: Bicheng Ying Kun Yuan Ali H Sayed
references.title: Variance-reduced stochastic learning under random reshuffling
references.venue: IEEE Transactions on Signal Processing,
references.authors: Rajput
references.title: non-convexity, variable step sizes, as well as techniques like variance reduction, momentum, etc. A.1 LOWER BOUNDS FOR VARIABLE STEP SIZES All the existing lower bounds (to the best of our knowledge) work in the constant step size regime (Safran
references.authors: Ahn
references.title: αnμ. Remark 2. To avoid confusion, we remark here that, the term ‘Sk
sections.text: 1 INTRODUCTION
  Finite sum optimization seeks to solve the following:
  min x F (x) :=
  1
  n n∑ i=1 fi(x). (1)
  Stochastic gradient descent (SGD) approximately solves finite sum problems, by iteratively updating the optimization variables according to the following rule:
  xt+1 := xt − α∇fσt(xt), (2)
  where α is the step size and σt ∈ [n] = {1, 2, . . . , n} is the index of the function sampled at iteration t. There exist various ways of sampling σt, with the most common being with- and without-replacement sampling. In the former, σt is uniformly chosen at random from [n], and for the latter, σt represents the t-th element of a random permutation of [n]. We henceforth refer to these two SGD variants as vanilla and permutation-based, respectively.
  Although permutation-based SGD has been widely observed to perform better in practice (Bottou, 2009; Recht & Ré, 2012; 2013), the vanilla version has attracted the vast majority of theoretical analysis. This is because of the fact that at each iteration, in expectation the update is a scaled version of the true gradient, allowing for simple performance analyses of the algorithm, e.g., see (Bubeck et al., 2015).
  Permutation-based SGD has resisted a tight analysis for a long time. However, a recent line of breakthrough results provides the first tight convergence guarantees for several classes of convex functions F (Nagaraj et al., 2019; Safran & Shamir, 2019; Rajput et al., 2020; Mishchenko et al., 2020; Ahn et al., 2020; Nguyen et al., 2020). These recent studies mainly focus on two variants of permutation-based SGD where (1) a new random permutation is sampled at each epoch (also known as RANDOM RESHUFFLE) (Nagaraj et al., 2019; Safran & Shamir, 2019; Rajput et al., 2020), and (2) a random permutation is sampled once and is reused throughout all SGD epochs (SINGLE SHUFFLE) (Safran & Shamir, 2019; Mishchenko et al., 2020; Ahn et al., 2020).
  Perhaps interestingly, RANDOM RESHUFFLE and SINGLE SHUFFLE exhibit different convergence rates and a performance gap that varies across different function classes. In particular, when run
  for K epochs, the convergence rate for strongly convex functions is Õ(1/nK2) for both RANDOM RESHUFFLE and SINGLE SHUFFLE (Nagaraj et al., 2019; Ahn et al., 2020; Mishchenko et al., 2020). However, when run specifically on strongly convex quadratics, RANDOM RESHUFFLE experiences an acceleration of rates, whereas SINGLE SHUFFLE does not (Safran & Shamir, 2019; Rajput et al., 2020; Ahn et al., 2020; Mishchenko et al., 2020). All the above rates have been coupled by matching lower bounds, at least up to constants and sometimes log factors (Safran & Shamir, 2019; Rajput et al., 2020).
  Algorithm 1 Permutation-based SGD variants Input: Initialization x10, step size α, epochs K 1: σ = a random permutation of [n] 2: for k = 1, . . . ,K do 3: if IGD then 4: σk = (1, 2, . . . , n) 5: else if SINGLE SHUFFLE then 6: σk = σ 7: else if RANDOM RESHUFFLE then 8: σk = a random permutation of [n] 9: end if
  10: if FLIPFLOP and k is even then 11: σk = reverse of σk−1 12: end if
  13: for i = 1, . . . , n do 14: xki := x k i−1 − α∇fσki (x k i−1)  Epoch k15: end for 16: xk+10 := x k n 17: end for
  From the above we observe that reshuffling at the beginning of every epoch may not always help. But then there are cases where RANDOM RESHUFFLE is faster than SINGLE SHUFFLE, implying that certain ways of generating permutations are more suited for certain subfamilies of functions.
  The goal of our paper is to take a first step into exploring the relationship between convergence rates and the particular choice of permutations. We are particularly interested in understanding if random permutations are as good as optimal, or if SGD can experience faster rates with carefully crafted permutations. As we see in the following, the answer the above is not straightforward, and depends heavily on the function class at hand.
  Our Contributions: We define as permutation-based SGD to be any variant of the iterates in (2), where a permutation of the n functions, at the start of each epoch, can be generated deterministically, randomly, or with a combination of the two. For example, SINGLE SHUFFLE, RANDOM RESHUFFLE, and incremental gradient descent (IGD), are all permutation-based SGD variants (see Algorithm 1).
  We first want to understand—even in the absence of computational constraints in picking the optimal permutations—what is the fastest rate one can get for permutation-based SGD? In other words, are there permutations that are better than random in the eyes of SGD?
  Perhaps surprisingly, we show that there exist permutations that may offer up to exponentially faster conver-
  gence than random permutations, but for a limited set of functions. Specifically, we show this for 1-dimensional functions (Theorem 1). However, such exponential improvement is no longer possible in higher dimensions (Theorem 2), or for general strongly convex objectives (Theorem 3), where random is optimal. The above highlight that an analysis of how permutations affect convergence rates
  needs to be nuanced enough to account for the structure of functions at hand. Otherwise, in lieu of further assumptions, random permutations may just appear to be as good as optimal.
  In this work, we further identify a subfamily of convex functions, where there exist easy-togenerate permutations that lead accelerated convergence. We specifically introduce a new technique, FLIPFLOP, which can be used in conjunction with existing permutation-based methods, e.g., RANDOM RESHUFFLE, SINGLE SHUFFLE, or INCREMENTAL GRADIENT DESCENT, to provably improve their convergence on quadratic functions (Theorems 4, 5, and 6). The way that FLIPFLOP works is rather simple: every even epoch uses the flipped (or reversed) version of the previous epoch’s permutation. The intuition behind why FLIPFLOP leads to faster convergence is as follows. Towards the end of an epoch, the contribution of earlier gradients gets attenuated. To counter this, we flip the permutation for the next epoch so that every function’s contribution is diluted (approximately) equally over the course of two consecutive epochs. FLIPFLOP demonstrates that finding better permutations for specific classes of functions might be computationally easy. We summarize FLIPFLOP’s convergence rates in Table 1 and report the results of numerical verification in Section 6.2.
  Note that in this work, we focus on the dependence of the error on the number of iterations, and in particular, the number of epochs. However, we acknowledge that its dependence on other parameters like the condition number is also very important. We leave such analysis for future work.
  Notation: We use lowercase for scalars (a), lower boldface for vectors (a), and upper boldface for matrices (A).
  2 RELATED WORK
  Gürbüzbalaban et al. (2019a;b) provided the first theoretical results establishing that RANDOM RESHUFFLE and INCREMENTAL GRADIENT DESCENT (and hence SINGLE SHUFFLE) were indeed faster than vanilla SGD, as they offered an asymptotic rate ofO ( 1/K2 ) for strongly convex functions,
  which beats the convergence rate ofO (1/nK) for vanilla SGD whenK = Ω(n). Shamir (2016) used techniques from online learning and transductive learning theory to prove an optimal convergence rate of Õ(1/n) for the first epoch of RANDOM RESHUFFLE (and hence SINGLE SHUFFLE). Later, Haochen & Sra (2019) also established a non-asymptotic convergence rate of Õ ( 1
  n2K2 + 1 K3
  ) , when
  the objective function is quadratic, or has smooth Hessian.
  Nagaraj et al. (2019) used a very interesting iterate coupling based approach to give a new upper bound on the error rate of RANDOM RESHUFFLE, thus proving for the first time that for general strongly convex smooth functions, it converges faster than vanilla SGD in all regimes of n and K. This was followed by (Safran & Shamir, 2019), where the authors were able to establish the first lower bounds, in terms of both n and K, for RANDOM RESHUFFLE. However, there was a gap in these upper and lower bounds. The gap in the convergence rates was closed by Rajput et al. (2020), who showed that the upper bound given by Nagaraj et al. (2019) and the lower bound given by Safran & Shamir (2019) were both tight up to logarithmic terms. For SINGLE SHUFFLE, Mishchenko et al. (2020) and Ahn et al. (2020) showed an upper bound of Õ (
  1 nK2
  ) , which matched the lower bound given earlier by (Safran & Shamir, 2019), up to logarithmic
  terms. Ahn et al. (2020) and Mishchenko et al. (2020) also proved tight upper bounds for RANDOM RESHUFFLE, with a simpler analysis and using more relaxed assumptions than (Nagaraj et al., 2019) and (Rajput et al., 2020). In particular, the results by Ahn et al. (2020) work under the PŁ condition and do not require individual component convexity.
  INCREMENTAL GRADIENT DESCENT on strongly convex functions has also been studied well in literature (Nedić & Bertsekas, 2001; Bertsekas, 2011; Gürbüzbalaban et al., 2019a). More recently, Nguyen et al. (2020) provide a unified analysis for all permutation-based algorithms. The dependence of their convergence rates on the number of epochs K is also optimal for INCREMENTAL GRADIENT DESCENT, SINGLE SHUFFLE and RANDOM RESHUFFLE.
  There has also been some recent work on the analysis of RANDOM RESHUFFLE on non-strongly convex functions and non-convex functions. Specifically, Nguyen et al. (2020) and Mishchenko et al. (2020) show that even there, RANDOM RESHUFFLE outperforms SGD under certain conditions. Mishchenko et al. (2020) show that RANDOM RESHUFFLE and SINGLE SHUFFLE beat vanilla SGD
  on non-strongly convex functions after Ω(n) epochs, and that RANDOM RESHUFFLE is faster than vanilla SGD on non-convex objectives if the desired error is O(1/ √ n).
  Speeding up convergence by combining without replacement sampling with other techniques like variance reduction (Shamir, 2016; Ying et al., 2020) and momentum (Tran et al., 2020) has also received some attention. In this work, we solely focus on the power of “good permutations” to achieve fast convergence.
  3 PRELIMINARIES
  We will use combinations of the following assumptions: Assumption 1 (Component convexity). fi(x)’s are convex. Assumption 2 (Component smoothness). fi(x)’s are L-smooth, i.e.,
  ∀x,y : ‖∇fi(x)−∇fi(y)‖ ≤ L‖x− y‖.
  Note that Assumption 2 immediately implies that F also has L-Lipschitz gradients:
  ∀x,y : ‖∇F (x)−∇F (y)‖ ≤ L‖x− y‖. Assumption 3 (Objective strong convexity). F is µ-strongly convex, i.e.,
  ∀x,y : F (y) ≥ F (x) + 〈∇F (x),y − x〉+ 1 2 µ‖y − x‖2.
  Note that Assumption 3 implies
  ∀x,y : 〈∇F (x)−∇F (y),x− y〉 ≥ µ‖y − x‖2. (3)
  We denote the condition number by κ, which is defined as κ = Lµ . It can be seen easily that κ ≥ 1 always. Let x∗ denote the minimizer of Eq. (1), that is, x∗ = arg minx F (x).
  We will study permutation-based algorithms in the constant step size regime, that is, the step size is chosen at the beginning of the algorithm, and then remains fixed throughout. We denote the iterate after the i-th iteration of the k-th epoch by xki . Hence, the initialization point is x 1 0. Similarly, the permutation of (1, . . . , n) used in the k-th epoch is denoted by σk, and its i-th ordered element is denoted by σki . Note that if the ambient space is 1-dimensional, then we represent the iterates and the minimizer using non-bold characters, i.e. xki and x ∗, to remain consistent with the notation.
  In the following, due to lack of space, we only provide sketches of the full proofs, when possible. The full proofs of the lemmas and theorems are provided in the Appendix.
  4 EXPONENTIAL CONVERGENCE IN 1-DIMENSION
  In this section, we show that there exist permutations for Hessian-smooth 1-dimensional functions that lead to exponentially faster convergence compared to random. Assumption 4 (Component Hessian-smoothness). fi(x)’s have LH -smooth second derivatives, that is,
  ∀x, y : |∇2fi(x)−∇2fi(y)| ≤ LH |x− y|.
  We also define the following instance dependent constants: G∗ := maxi ‖∇fi(x∗)‖, D = max { ‖x10 − x∗‖, G ∗
  2L
  } , and G = G∗ + 2DL.
  Theorem 1. Let Assumptions 1,2,3 and 4 hold. LetD andG be as defined above. If α = µ8n(L2+LHG) , then there exists a sequence of permutations σ1, σ2, . . . , σK such that using those permutations from any initialization point x10 gives the error
  |xKn − x∗| ≤ (D + 4nαG)e−CK ,
  where C = µ 2
  16(L2+LHG) .
  An important thing to notice in the theorem statement is that the sequence of permutations σ1, σ2, . . . , σK only depends on the function, and not on the initialization point x10. This implies that for any such function, there exists a sequence of permutations σ1, σ2, . . . , σK , which gives exponentially fast convergence, unconditionally of the initialization. Note that the convergence rate is slower than Gradient Descent, for which the constant ‘C’ would be larger. However, here we are purely interested in the convergence rates of the best permutations and their (asymptotic) dependence on K.
  Proof sketch The core idea is to establish that there exists an initialization point x10 (close to the minimizer x∗), and a sequence of permutations such that that starting from x10 and using that sequence of permutation leads us exactly to the minimizer. Once we have proved this, we show that if two parallel runs of the optimization process are initialized from two different iterates, and they are coupled so that they use the exact same permutations, then they approach each other at an exponential rate. Thus, if we use the same permutation from any initialization point, it will converge to the minimizer at an exponential rate. See Figure 1 for a graphical
  depiction of this sketch. We note that the figure is not the result of an actual run, but only serves to explain the proof sketch.
  5 LOWER BOUNDS FOR PERMUTATION-BASED SGD
  The result in the previous section leads us to wonder if exponentially fast convergence can be also achieved in higher dimensions. Unfortunately, for higher dimensions, there exist strongly convex quadratic functions for which there does not exist any sequence of permutations that lead to an exponential convergence rate. This is formalized in the following theorem. Theorem 2. For any n ≥ 4 (n must be even), there exists a 2n + 1-dimensional strongly convex function F which is the mean of n convex quadratic functions, such that for every permutation-based algorithm with any step size,
  ‖xK−1n − x∗‖2 = Ω ( 1
  n3K2
  ) .
  This theorem shows that we cannot hope to develop constant step size algorithms that exhibit exponentially fast convergence in multiple dimensions.
  Proof sketch In this sketch, we explain a simpler version of the theorem, which works for n = 2 and 2-Dimensions. Consider F (x, y) = 12f1(x, y) + 1 2f2(x, y) such that
  f1(x, y) = x2
  2 − x+ y, andf2(x, y) =
  y2
  2 − y + x.
  Hence, F = 14 (x 2 + y2), and has minimizer at the origin. Each epoch has two possible permutations, σ = (1, 2) or σ = (2, 1). Working out the details manually, it can be seen that regardless of the permutation, xk+10 + y k+1 0 > x k 0 + y k 0 , that is, the sum of the co-ordinates keeps increasing. This can be used to get a bound on the error term ‖[xK yK ]>‖2. Next, we show that even in 1-Dimension, individual function convexity might be necessary to obtain faster rates than RANDOM RESHUFFLE.
  Theorem 3. There exists a 1-Dimensional strongly convex function F which is the mean of two quadratic functions f1 and f2, such that one of the functions is non-convex. Then, every permutationbased algorithm with constant step size α ≤ 1L gives an error of at least
  ‖xK−1n − x∗‖2 = Ω ( 1
  K2
  ) .
  Proof sketch The idea behind the sketch is to have one of the two component functions as strongly concave. This gives it the advantage that the farther away from its maximum the iterate is, the more it pushes the iterate away. Hence, it essentially results in increasing the deviation in each epoch. This leads to a slow convergence rate.
  In the setting where the individual fi may be non-convex, Nguyen et al. (2020) and Ahn et al. (2020) show that SINGLE SHUFFLE, RANDOM RESHUFFLE, and INCREMENTAL GRADIENT DESCENT achieve the error rate of Õ( 1K2 ), for a fixed n. In particular, their results only need that the component functions be smooth and hence their results apply to the function F from Theorem 3. The theorem above essentially shows that when n = 2, this is the best possible error rate, for any permutationbased algorithm - deterministic or random. Hence, at least for n = 2, the three algorithms are optimal when the component functions can possibly be non-convex. However, note that here we are only considering the dependence of the convergence rate on K. It is possible that these are not optimal, if we further take into account the dependence of the convergence rate on the combination of both n and K. Indeed, if we consider the dependence on n as well, INCREMENTAL GRADIENT DESCENT has a convergence rate of Ω(1/K2) (Safran & Shamir, 2019), whereas the other two have a convergence rate of Õ(1/nK2) (Ahn et al., 2020).
  6 FLIPPING PERMUTATIONS FOR FASTER CONVERGENCE IN QUADRATICS
  In this section, we introduce a new algorithm FLIPFLOP, that can improve the convergence rate of SINGLE SHUFFLE,RANDOM RESHUFFLE, and INCREMENTAL GRADIENT DESCENT on strongly convex quadratic functions.
  The following theorem gives the convergence rate of FLIPFLOP WITH SINGLE SHUFFLE: Assumption 5. fi(x)’s are quadratic. Theorem 4. If Assumptions 1, 2, 3, and 5 hold, then running FLIPFLOP WITH SINGLE SHUFFLE for K epochs, where K ≥ 80κ3/2 log(nK) max { 1, √ κ n } is an even integer, with step size α =
  10 log(nK) µnK gives the error E [∥∥xK−1n − x∗∥∥2] = Õ( 1n2K2 + 1nK4 ) . (4)
  For comparison, Safran & Shamir (2019) give the following lower bound on the convergence rate of vanilla SINGLE SHUFFLE: E [∥∥xK−1n − x∗∥∥2] = Ω( 1nK2 ) , (5)
  Note that both the terms in Eq. (4) are smaller than the term in Eq. (5). In particular, when n K2 and n is fixed as we vary K, the RHS of Eq. (5) decays as Õ( 1K2 ), whereas the RHS of Eq. (4) decays as Õ( 1K4 ). Otherwise, when K
  2 n and K is fixed as we vary n, the RHS of Eq. (5) decays as Õ( 1n ), whereas the RHS of Eq. (4) decays as Õ( 1 n2 ). Hence, in both the cases, FLIPFLOP WITH SINGLE SHUFFLE outperforms SINGLE SHUFFLE.
  The next theorem shows that FLIPFLOP improves the convergence rate of RANDOM RESHUFFLE: Theorem 5. If Assumptions 1, 2, 3, and 5 hold, then running FLIPFLOP WITH RANDOM RESHUFFLE for K epochs, where K ≥ 55κ log(nK) max { 1, √
  n κ
  } is an even integer, with step size α =
  10 log(nK) µnK gives the error E [∥∥xK−1n − x∗∥∥2] = Õ( 1n2K2 + 1nK5 ) .
  For comparison, Safran & Shamir (2019) give the following lower bound on the convergence rate of vanilla RANDOM RESHUFFLE:
  E [∥∥xK−1n − x∗∥∥2] = Ω( 1n2K2 + 1nK3 ) .
  Hence, we see that in the regime when n K, which happens when the number of components in the finite sum of F is much larger than the number of epochs, FLIPFLOP WITH RANDOM RESHUFFLE is much faster than vanilla RANDOM RESHUFFLE.
  Note that the theorems above do not contradict Theorem 2, because for a fixed n, both the theorems above give a convergence rate of Õ(1/K2). We also note here that the theorems above need the number of epochs to be much larger than κ, in which range Gradient Descent performs better than with- or without- replacement SGD, and hence GD should be preferred over SGD in that case. However, we think that this requirement on epochs is a limitation of our analysis, rather than that of the algorithms.
  Finally, the next theorem shows that FLIPFLOP improves the convergence rate of INCREMENTAL GRADIENT DESCENT. Theorem 6. If Assumptions 1, 2, 3, and 5 hold, then running FLIPFLOP WITH INCREMENTAL GD for K epochs, where K ≥ 36κ log(nK) is an even integer, with step size α = 6 lognKµnK gives the error
  E [∥∥xK−1n − x∗∥∥2] = Õ( 1n2K2 + 1K3 ) .
  For comparison, Safran & Shamir (2019) give the following lower bound on the convergence rate of vanilla INCREMENTAL GRADIENT DESCENT:
  E [∥∥xK−1n − x∗∥∥2] = Ω( 1K2 ) ,
  In the next subsection, we give a short sketch of the proof of these theorems.
  6.1 PROOF SKETCH
  In the proof sketch, we consider scalar quadratic functions. The same intuition carries over to multi-dimensional quadratics, but requires a more involved analysis. Let fi(x) := aix 2
  2 + bix+ c. Assume that F (x) := 1n ∑n i=1 fi(x) has minimizer at 0. This assumption is valid because it can be achieved by a simple translation of the origin (see (Safran & Shamir, 2019) for a more detailed explanation). This implies that ∑n i=1 bi = 0.
  For the sake of this sketch, assume x10 = 0, that is, we are starting at the minimizer itself. Further, without loss of generality, assume that σ = (1, 2, . . . , n). Then, for the last iteration of the first epoch,
  x1n = x 1 n−1 − αf ′n(x1n−1)
  = x1n−1 − α(anx1n−1 + bn) = (1− αan)x1n−1 − αbn.
  Applying this to all iterations of the first epoch, we get
  x1n = n∏ i=1 (1− αai)x10 − α n∑ i=1 bi n∏ j=i+1 (1− αaj). (6)
  Substituting x10 = 0, we get
  x1n = −α n∑ i=1 bi n∏ j=i+1 (1− αaj). (7)
  Note that the sum above is not weighted uniformly: b1 is multiplied by ∏n j=2(1− αaj), whereas bn is multiplied by 1. Because (1− αaj) < 1, we see that b1’s weight is much smaller than bn. If the weights were all 1, then we would get x0n = −α ∑n i=1 bi = 0, i.e., we would not move away from the minimizer. Since we want to stay close to the minimizer, we want the weights of all the bi to be roughly equal.
  The idea behind FLIPFLOP is to add something like −α ∑n i=1 bi ∏i j=1(1− αaj) in the next epoch, to counteract the bias in Eq. (7). To achieve this, we simply take the permutation that the algorithm used in the previous epoch and flip it for the next epoch. Roughly speaking, in the next epoch b1 will get multiplied by 1 whereas bn will get multiplied by ∏n−1 j=1 (1− αaj). Thus over two epochs, both get scaled approximately the same.
  To see more concretely, we look at the first order approximation of Eq. (7):
  x1n = −α n∑ i=1 bi n∏ j=i+1 (1− αaj) ≈ −α n∑ i=1 bi(1− n∑ j=i+1 αaj) = α 2 n∑ i=1 bi n∑ j=i+1 aj , (8)
  where in the last step above we used the fact that ∑n i=1 bi = 0. Now, let us see what happens in the second epoch if we use FLIPFLOP. Doing a recursion analogous to how we got Eq. (6), but reversing the order of functions, we get:
  x2n = n∏ i=1 (1− αan−i+1)x20 − α n∑ i=1 bn−i+1 n∏ j=i+1 (1− αan−j+1).
  Carefully doing a change of variables in the equation above, we get:
  x2n = n∏ i=1 (1− αai)x20 − α n∑ i=1 bi i−1∏ j=1 (1− αaj). (9)
  Note that the product in the second term in the equation above is almost complementary to the product in Eq. (7). This is because we flipped the order in the second epoch. This will play an important part in cancelling out much of the bias in Eq. (7). Continuing on from Eq. (9), we again do a linear approximation similar to before and substitute Eq. (8) (and use the fact that x20 = x 1 n):
  x2n ≈
  ( 1− α
  n∑ i=1 ai ) x20 + α 2 n∑ i=1 bi i−1∑ j=1 aj
  ≈ ( 1− α
  n∑ i=1 ai )α2 n∑ i=1 bi n∑ j=i+1 aj + α2 n∑ i=1 bi i−1∑ j=1 aj .
  We assume that α is small and hence we only focus on the quadratic terms:
  x2n = α 2 n∑ i=1 bi ∑ j 6=i aj +O(α 3)
  = α2  n∑ i=1 bi n∑ j=1 aj − α2( n∑ i=1 biai ) +O(α3).
  Now, since ∑n i=1 bi = 0, we get
  x2n ≈ −α2 ( n∑ i=1 biai ) +O(α3). (10)
  Now, comparing the coefficients of the α2 terms in Eq. (8) and Eq. (10), we see that the former has O(n2) terms whereas the latter has only n terms. This correction of error is exactly what eventually manifests into the faster convergence rate of FLIPFLOP.
  The main reason that the analysis for multidimensional quadratics is not as simple as the 1-dimensional analysis done above, is because unlike scalar multiplication, matrix multiplication is not commutative, and the AM-GM inequality is not true in higher dimensions (Lai & Lim, 2020; De Sa, 2020). One way to bypass this inequality is by using the following inequality for small enough α:∥∥∥∥∥ n∏ i=1 (I − αAi) n∏ i=1 (I − αAn−i+1)
  ∥∥∥∥∥ ≤ 1− αnµ. Ahn et al. (2020) proved a stochastic version of this (see Lemma 6 in their paper). We prove a deterministic version in Lemma 3 (in the Appendix).
  6.2 NUMERICAL VERIFICATION
  We verify the theorems numerically by running RANDOM RESHUFFLE, SINGLE SHUFFLE and their FLIPFLOP versions on the task of mean computation. We randomly sample n = 800 points from a 100-dimensional sphere. Let the points be xi for i = 1, . . . , n. Then, their mean is the solution to the following quadratic problem : arg minx F (x) = 1n ∑n i=1 ‖x− xi‖2. We solve this problem by using the given algorithms. The results are reported in Figure 2. The results are plotted in a log–log graph, so that we get to see the dependence of error on the power of K.
  Note that since the points are sampled randomly, INCREMENTAL GRADIENT DESCENT essentially becomes SINGLE SHUFFLE. Hence, to verify Theorem 6, we need ‘hard’ instances of INCREMENTAL GRADIENT DESCENT, and in particular we use the ones used in Theorem 3 in (Safran & Shamir, 2019). These results are also reported in a log–log graph in Figure 2.
  We also tried FLIPFLOP in the training of deep neural networks, but unfortunately we did not see a big speedup there. We ran experiments on logistic regression for 1-Dimensional artificial data, the results for which are in Appendix H.
  7 CONCLUSION AND FUTURE WORK
  In this paper, we explore the theoretical limits of permutation-based SGD for solving finite sum optimization problems. We focus on the power of good, carefully designed permutations and whether they can lead to a much better convergence rate than random. We prove that for 1-dimensional, strongly convex functions, indeed good sequences of permutations exist, which lead to a convergence rate which is exponentially faster than random permutations. We also show that unfortunately, this is not true for higher dimensions, and that for general strongly convex functions, random permutations might be optimal.
  However, we think that for some subfamilies of strongly convex functions, good permutations might exist and may be easy to generate. Towards that end, we introduce a very simple technique, FLIPFLOP, to generate permutations that lead to faster convergence on strongly convex quadratics. This is a black box technique, that is, it does not look at the optimization problem to come up with the permutations; and can be implemented easily. This serves as an example that for other classes of functions, there can exist other techniques for coming up with good permutations. Finally, note that we only consider constant step sizes in this work for both upper and lower bounds. Exploring regimes in which the step size changes, e.g., diminishing step sizes, is a very interesting open problem, which we leave for future work. We think that the upper and lower bounds in this paper give some important insights and can help in the development of better algorithms or heuristics. We strongly believe that under nice distributional assumptions on the component functions, there can exist good heuristics to generate good permutations, and this should also be investigated in future work.
  ETHICS STATEMENT
  This paper explores the theoretical limits of convergence rates of variants of popular SGD type algorithms. Thus, the authors do not think there are any ethical concerns regarding the paper.
  REPRODUCIBILITY STATEMENT
  The proof for all the theoretical claims are provided in the Appendix. The assumptions and notations are specified in the Preliminaries section (Section 3). The code for all the experiments is provided in the supplementary materials zip file.
  REFERENCES Kwangjun Ahn, Chulhee Yun, and Suvrit Sra. Sgd with shuffling: optimal rates without component
  convexity and large epoch requirements, 2020.
  Dimitri P Bertsekas. Incremental gradient, subgradient, and proximal methods for convex optimization: A survey. Optimization for Machine Learning, 2010(1-38):3, 2011.
  Léon Bottou. Curiously fast convergence of some stochastic gradient descent algorithms. Unpublished open problem offered to the attendance of the SLDS 2009 conference, 2009. URL http: //leon.bottou.org/papers/bottou-slds-open-problem-2009.
  Sébastien Bubeck et al. Convex optimization: Algorithms and complexity. Foundations and Trends® in Machine Learning, 8(3-4):231–357, 2015.
  Christopher M De Sa. Random reshuffling is not always better. Advances in Neural Information Processing Systems, 33, 2020.
  M Gürbüzbalaban, Asuman Ozdaglar, and Pablo A Parrilo. Convergence rate of incremental gradient and incremental newton methods. SIAM Journal on Optimization, 29(4):2542–2565, 2019a.
  Mert Gürbüzbalaban, Asu Ozdaglar, and Pablo A Parrilo. Why random reshuffling beats stochastic gradient descent. Mathematical Programming, pp. 1–36, 2019b.
  Jeff Haochen and Suvrit Sra. Random shuffling beats sgd after finite epochs. In International Conference on Machine Learning, pp. 2624–2633, 2019.
  Zehua Lai and Lek-Heng Lim. Recht-ré noncommutative arithmetic-geometric mean conjecture is false. In International Conference on Machine Learning, pp. 5608–5617. PMLR, 2020.
  Konstantin Mishchenko, Ahmed Khaled, and Peter Richtárik. Random reshuffling: Simple analysis with vast improvements. ArXiv, abs/2006.05988, 2020.
  Dheeraj Nagaraj, Prateek Jain, and Praneeth Netrapalli. Sgd without replacement: Sharper rates for general smooth convex functions. In International Conference on Machine Learning, pp. 4703–4711, 2019.
  Angelia Nedić and Dimitri Bertsekas. Convergence rate of incremental subgradient algorithms. In Stochastic optimization: algorithms and applications, pp. 223–264. Springer, 2001.
  Yurii Nesterov. Introductory Lectures on Convex Optimization - A Basic Course, volume 87 of Applied Optimization. Springer, 2004. ISBN 978-1-4613-4691-3. doi: 10.1007/978-1-4419-8853-9. URL https://doi.org/10.1007/978-1-4419-8853-9.
  Lam M Nguyen, Quoc Tran-Dinh, Dzung T Phan, Phuong Ha Nguyen, and Marten van Dijk. A unified convergence analysis for shuffling-type gradient methods. arXiv preprint arXiv:2002.08246, 2020.
  Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science & Business Media, 2006.
  Shashank Rajput, Anant Gupta, and Dimitris Papailiopoulos. Closing the convergence gap of sgd without replacement. In International Conference on Machine Learning, pp. 7964–7973. PMLR, 2020.
  Benjamin Recht and Christopher Ré. Beneath the valley of the noncommutative arithmetic-geometric mean inequality: conjectures, case-studies, and consequences. arXiv preprint arXiv:1202.4184, 2012.
  Benjamin Recht and Christopher Ré. Parallel stochastic gradient algorithms for large-scale matrix completion. Mathematical Programming Computation, 5(2):201–226, 2013.
  Itay Safran and Ohad Shamir. How good is sgd with random shuffling?, 2019.
  Itay Safran and Ohad Shamir. Random shuffling beats sgd only after many epochs on ill-conditioned problems. arXiv preprint arXiv:2106.06880, 2021.
  M. Schneider. Probability inequalities for kernel embeddings in sampling without replacement. In AISTATS, 2016.
  Ohad Shamir. Without-replacement sampling for stochastic gradient methods. In Proceedings of the 30th International Conference on Neural Information Processing Systems, pp. 46–54, 2016.
  Trang H Tran, Lam M Nguyen, and Quoc Tran-Dinh. Shuffling gradient-based methods with momentum. arXiv preprint arXiv:2011.11884, 2020.
  Bicheng Ying, Kun Yuan, and Ali H Sayed. Variance-reduced stochastic learning under random reshuffling. IEEE Transactions on Signal Processing, 68:1390–1408, 2020.
  A DISCUSSION AND FUTURE WORK
  Being the first paper (to the best of our knowledge) to theoretically analyze the optimality of random permutations, we limited our scope to a specific, but common theoretical setting - strongly convex functions with constant step size. We think that future work can generalize the results of this paper to settings of non-convexity, variable step sizes, as well as techniques like variance reduction, momentum, etc.
  A.1 LOWER BOUNDS FOR VARIABLE STEP SIZES
  All the existing lower bounds (to the best of our knowledge) work in the constant step size regime (Safran & Shamir (2019); Rajput et al. (2020); Safran & Shamir (2021)). Thus, generalizing the lower bounds to variable step size algorithms would be a very important direction for future research.
  However, the case when step sizes are not constant can be tricky to prove lower bounds, since the step size could potentially depend on the permutation, and the current iterate. A more reasonable setting to prove lower bounds could be the case when the step sizes follow a schedule over epochs, similar to what happens in practice.
  A.2 FASTER PERMUTATIONS FOR NON-QUADRATIC OBJECTIVES
  The analysis of FlipFlop leverages the fact that the Hessians of quadratic functions are constant. We think that the analysis of FlipFlop might be extended to strongly convex functions or even PŁ functions (which are non-convex in general), under some assumptions on the Lipschitz continuity of the Hessians, similar to how Haochen & Sra (2019) extended their analysis of quadratic functions to more general classes. A key take-away from FlipFlop is that we had to understand how permutation based SGD works specifically for quadratic functions, that is, we did a white-box analysis. In general, we feel that depending on the specific class of non-convex functions (say deep neural networks), practitioners would have to think about permutation-based SGD in a white-box fashion, to come up with better heuristics for shuffling.
  In a concurrent work, (https://openreview.net/pdf?id=7gWSJrP3opB), it is shown that by greedily sorting stale gradients, a permutation order can be found which converges faster for some deep learning tasks. Hence, there do exist better permutations than random, even for deep learning tasks.
  A.3 FLIPFLOP ON RANDOM COORDINATE DESCENT
  A shuffling scheme similar to FlipFlop has been used in random coordinate descent for faster practical convergence (see page 231 in Nocedal & Wright (2006)). This should be further investigated empirically and theoretically in future work. Even though the current analysis of FlipFlop does not directly go through for random coordinate descent, we think the analysis can be adapted to work.
  B PROOF OF THEOREM 1
  In certain places in the proof, we would need that α ≤ 14nL . To see how this is satisfied, note that we have assumed that α ≤ µ8n(L2+LHG) in the theorem statement. Using the inequality µ ≤ L in α ≤ µ4n(L2+LHG) gives that α ≤ 1 4n(L+(LHG/µ)) ≤ 14nL .
  In this proof, we assume that the minimizer of F is at 0 to make the analysis simpler. This assumption can be satisfied by simply translating the origin to the minimizer (Safran & Shamir, 2019).
  There are three main components in the proof:
  1. Say that an epoch starts off at the minimizer. We prove that there exists at least one pair of permutations such that if we do two separate runs of the epoch, the first using the first permutation and the second using the second, then at the end of that epoch, the iterates corresponding to the two runs end up on opposite sides of the minimizer.
  2. There exists a sequence of permutations and a point in the neighborhood of the minimizer, such that intializing at that point and using these sequence of permutations, we converge exactly to the minimizer.
  3. Starting from any other point, we can couple the iterates with the iterates which were shown in the previous component, to get that these two sequences of iterates come close to each other exponentially fast.
  We prove the first and second components in the Subsections B.1 and B.2; and conclude the proof in Subsection B.3 where we also prove the third component.
  B.1 PERMUTATIONS IN ONE EPOCH
  In this subsection, we prove that if x0, x1, . . . , xn are the iterates in an epoch such that x0 = 0, then there exists a permutation of functions such that xn ≥ 0. By the same logic, we show that there exists a permutation of functions such that xn ≤ 0. These will give us control over movement of iterates across epochs.
  Order the gradients at the minimizer, ∇fi(0), in decreasing order. WLOG assume that it is just ∇f1,∇f2, . . . ,∇fn. We claim that this permutation leads to xn ≥ 0. We will need the following intermediate result. Let yi, yi−1 be such that yi = yi−1 − α∇fi(yi−1). Assume α ≤ 1/L and yi−1 ≥ xi−1. Then,
  yi − xi = yi−1 − xi−1 − α(∇fi(yi−1)−∇fi(xi−1)) ≥ yi−1 − xi−1 − αL(yi−1 − xi−1) = (1− αL)(yi−1 − xi−1) ≥ 0, (11)
  that is, yi−1 ≥ xi−1 =⇒ yi ≥ xi. Because 0 is the minimizer, we know that ∑n i=1∇fi(0) = 0. Also, recall that x0 = 0. There can be two cases:
  1. ∀i ∈ [1, n] : xi < 0. This cannot be true because if ∀i : 1 ≤ i ≤ n− 1 : xi < 0, then
  xn = n∑ i=1 −α∇fi(xi−1) ≥ −α n∑ i=1 ∇fi(0) ≥ 0,
  where we used the fact that∇fi(x) ≤ ∇fi(y) if x ≤ y and fi is convex. 2. Thus, ∃i ∈ [1, n] : xi ≥ 0. Now, consider the sequence yi, yi+1, . . . , yn such that yi = 0
  and for j ≥ i+ 1, yj = yj−1 − α∇fj(yj−1). Then because α ≤ 1/L and xi ≥ yi = 0, we get that xj ≥ yj for all j ≥ i (Using Ineq. (11)).
  Hence, there is an i ≥ 1 such that xi ≥ yi = 0, and further xj ≥ yj for all j ≥ i. Next, we repeat the process above for yi, . . . , yn. That is, there can be two cases:
  1. ∀j ∈ [i+ 1, n] : yj < 0. This cannot be true because if ∀j : i+ 1 ≤ j ≤ n− 1 : yj < 0, then
  yn = n∑ j=i −α∇fj(yj) ≥ −α n∑ j=i ∇fj(0) ≥ 0.
  2. Thus, ∃j ∈ [i + 1, n] : yj ≥ 0. Now, consider the sequence zj , zj+1, . . . , zn such that zj = 0 and for k ≥ j + 1, zk = zk−1 − α∇fk(zk−1). Then because α ≤ 1/L and yj ≥ zj = 0, we get that yk ≥ zk for all k ≥ j (Using Ineq. (11)).
  Hence, there exists an integer j > i > 0 such that yj ≥ 0. We have already proved that xj ≥ yj . Thus, we have that xj ≥ 0. We can continue repeating this process (apply the same two cases above for zj , zj+1, . . . , zn, and so on), to get that xn ≥ 0. We define p to be this non-negative value of xn. Note that the following lemma gives us that the gradients are bounded by G
  Lemma 1. Define G∗ := maxi ‖∇fi(x∗)‖ and D = max { ‖x10 − x∗‖, G ∗
  2L
  } . If Assumptions 2 and
  3 hold, and α < 18κnL , then for any permutation-based algorithm (deterministic or random), we have
  ∀i, j, k :‖xki − x∗‖ ≤ 2D, and ‖∇fj(xki )‖ ≤ G∗ + 2DL.
  Because the gradients are bounded by G, we get that p ≤ nαG. Similarly, we can show that the reverse permutation leads to xn ≤ 0. We define q to be this nonpositive value of xn. Because we have assumed that the gradients are bounded by G, we get that q ≥ −nαG.
  B.2 EXACT CONVERGENCE TO THE MINIMIZER
  In this section, we show that there exists a point such that if we initialize there and follow a specific permutation sequence, we land exactly at the minimizer.
  In particular, we will show the following: There exists a point in [4q, 4p] such that if we initialize there and follow a specific permutation sequence, then we land exactly at the minimizer.
  We show this recursively. We will prove that there exists a point mK ∈ [4q, 4p] such that if the last epoch begins there, that is xK0 = m
  K , then we land at the minimizer at the end of the last epoch, that is, xKn = 0. Then, we will show that there exists a point m
  K−1 ∈ [4q, 4p] such that if the xK−10 = m K−1, then xK0 = x K−1 n = m
  K . Repeating this for K − 2, . . . , 1, we get that there exists a point m0 ∈ [4q, 4p] such that if we initialize the first epoch there, that is x10 = m0, then there is a permutation sequence such that ultimately xKn = 0.
  We prove that any point mj ∈ [4q, 4p] can be reached at the end of an epoch by beginning the epoch at some point mj−1 ∈ [4q, 4p], that is if xj−10 = mj−1, then xj−1n = mj .
  • Case: mj ∈ [p, 4p]. In this case, we show that mj−1 ∈ [0, 4p]. We have proved in the previous subsection that there exists a permutation σ such that if xj−10 = 0 then x j−1 n = p.
  Next, we have the following helpful lemma that we will also use later. Lemma 2. Let x0, x1, . . . , xn be a sequence of iterates in an epoch and y0, y1, . . . , yn be another sequence of iterates in an epoch such that both use the same permutation of functions. If α ≤ µ2n(L2+LHG) , then
  (1− nαL)|y0 − x0| ≤ (1− Lα)n|y0 − x0| ≤ |yn − xn| ≤ (
  1− 1 2 nµα
  ) |y0 − x0|.
  If we set x0 = 4p, y0 = 0 in Lemma 2 and we follow the permutation σ, then we get that xn − yn ∈ (x0 − y0) [ 1− αnL, 1− αnµ
  2 ] =⇒ xn − p ∈ (4p− 0) [ 1− αnL, 1− αnµ
  2 ] (Since using σ and y0 = 0 results in yn = p.)
  =⇒ xn ≥ 4p,
  where we used the fact that α ≤ 14nL is the last step. Thus, if xj−10 = 4p and we follow the permutation σ, then x j−1 n ≥ 4p. Next, note that
  xj−11 = x j−1 0 − α∇fσ(0)(x j−1 0 ). xj−12 = x j−1 1 − α∇fσ(1)(x j−1 1 )
  ...
  xj−1n = x j−1 n−1 − α∇fσ(n−1)(x j−1 n−1)
  Looking above, we see that xj−11 is a continuous function of x j−1 0 ; x j−1 2 is a continuous function of xj−11 ; and so on. Thus, using the fact that composition of continuous functions is continuous, we get that xj−1n is also a continuous function of x j−1 0 . We have shown that if xj−10 = 0, then x j−1 n = p and if x j−1 0 = 4p, then x j−1 n ≥ 4p. Thus, using the fact that that xj−1n is a continuous function of x j−1 0 , we get that for any point m
  j ∈ [p, 4p], there is at least one point mj−1 ∈ [0, 4p], such that xj−10 = mj−1 leads to xj−1n = mj . • Case: mj ∈ [4q, q]. We can apply the same logic as above to show that there is at least one point mj−1 ∈ [4q, 0], such that xj−10 = mj−1 =⇒ xj−1n = mj . • Case: mj ∈ [q, p]. WLOG assume that |q| < |p|. Let σq be the permutation such that if xj−10 = 0 and the epoch uses this permutation, then we end up at x j−1 n = q.
  If we set x0 = 4p, y0 = 0 in Lemma 2 and we follow the permutation σq , then we get that
  xn − yn ∈ (x0 − y0) [ 1− αnL, 1− αnµ
  2 ] =⇒ xn − q ∈ (4p− 0) [ 1− αnL, 1− αnµ
  2 ] (Since using σq and y0 = 0 results in yn = q.)
  =⇒ xn ≥ q + 4p(1− αnL) ≥ q + 3p ≥ 2p,
  where we used the fact that α ≤ 14nL is the last step. Thus, if xj−10 = 4p and we follow the permutation σq , then x j−1 n ≥ 2p. Thus, using similar argument of continuity as the first case, we know that there is a point mj−1 ∈ [0, 4p], such that xj−10 = mj−1 leads to xj−1n = mj when we use the permutation σq .
  B.3 SAME SEQUENCE PERMUTATIONS GET CLOSER
  In the previous subsection, we have shown that there exists a point m0 ∈ [4q, 4p] and a sequence of permutations σ1, σ2, . . . , σK such that if x10 = m
  0 and epoch j uses permutation σj , then xKn = 0. In this subsection, we will show that if x10 is initialized at any other point such that |x10| ≤ D, then using the same permutations σ1, σ2, . . . , σK gives us that |xKn | ≤ (D + nαG)e−K . For this we will repeatedly apply Lemma 2 on all the K epochs. Assume that x10 = ν 0 with |ν0| ≤ D.
  Let yji be the sequence of iterates such that y 1 0 = m 0 and uses the permutation sequence σ1, σ2, . . . , σK . Then, we know that yKn = 0. Let x j i be the sequence of iterates such that x 1 0 = ν 0 and uses the same permutation sequence σ1, σ2, . . . , σK .
  Then, using Lemma 2 gives us that |y1n − x1n| ≤ |ν0 −m0|(1− µαn 2 ). Thus, we get that |y 2 0 − x20| ≤ |ν0 −m0|(1 − µαn2 ). Again applying Lemma 2 gives us that |y 2 n − x2n| ≤ |ν0 −m0|(1 − µαn 2 )
  2. Therefore, after applying it K times, we get
  |yKn − xKn | ≤ |ν0 −m0| (
  1− µαn 2
  )K .
  Note that |x10| = |ν0| ≤ D, and |y10 | = |m0|, with m0 ∈ [4q, 4p]. We showed earlier in Subsection B.1 that |p| ≤ nαG and |q| ≤ nαG. Therefore,
  |yKn − xKn | ≤ |D + 4nαG| (
  1− µαn 2
  )K .
  Further, we know that yKn = 0. Thus, |xKn | ≤ |D + 4nαG| (
  1− µαn 2 )K ≤ |D + 4nαG|e− µαn 2 K .
  Substituting the value of α completes the proof. Next we prove the lemmas used in this proof.
  B.4 PROOF OF LEMMA 1
  We restate the lemma below. Lemma. Define G∗ := maxi ‖∇fi(x∗)‖ and D = max { ‖x10 − x∗‖, G ∗
  2L
  } . If Assumptions 2 and 3
  hold, and α < 18κnL , then for any permutation-based algorithm (deterministic or random), we have
  ∀i, j, k :‖xki − x∗‖ ≤ 2D, and ‖∇fj(xki )‖ ≤ G∗ + 2DL.
  This lemma says that for any permutation-based algorithm, the domain of iterates and the norm of the gradient stays bounded during the optimization. This means that we can assume bounds on norm of iterates and gradients, which is not true in general for unconstrained SGD. This makes analyzing such algorithms much easier, and hence this lemma can be of independent interest for proving future results for permutation-based algorithms. Remark 1. This lemma does not hold in general for vanilla SGD where sampling is done with replacement. Consider the example with two functions f1(x) = x2 − x, and f2(x) = x; and F (x) = f1(x) + f2(x). This satisfies Assumptions 2 and 3, but one may choose f2 consecutively for arbitrary many iterations, which will lead the iterates a proportional distance away from the minimizer. This kind of situation can never happen for permutation-based SGD because we see every function exactly once in every epoch and hence no particular function can attract the iterates too much towards its minimizer, and by the end of the epochs most of the noise gets cancelled out.
  Note that the requirement α < 18κnL is stricter than the usual requirement α = O ( 1 nL ) , but we believe the lemma should also hold in that regime. For the current paper, however, this stronger requirement suffices.
  Proof. To prove this lemma, we show two facts: If for some epoch k, ‖xk0 − x∗‖ ≤ D, then a) ∀i : ‖xki − x∗‖ ≤ 2D and b) ‖x k+1 0 − x∗‖ = ‖xkn − x∗‖ ≤ D. To see how a) and b) are sufficient to prove the lemma, assume that they are true. Then, since the first epoch begins inside the bounded region ‖x10−x∗‖ ≤ D, we see using b) that every subsequent epoch begins inside the same bounded region, that is ‖xk0 − x∗‖ ≤ D as well. Hence using a) we get that during these epochs, the iterates satisfy ‖xki − x∗‖ ≤ 2D, which is the first part of the lemma. Further, this bound together with the gradient Lipshitzness directly gives the upper bound G∗ + 2DL on the gradients. Thus, all we need to do to prove this lemma is to prove a) and b), which we do next.
  We will prove a) and b) for D = max{‖x10 − x∗‖, 2κnαG ∗ 1−4κnαL}. Once we do this, using α < 1 8κnL will give us the exact statement of the lemma.
  Let |xk0 − x∗| ≤ D for some epoch k. Then, we try to find the minimum number of iterations i needed so that |xki − x∗| ≥ 2D. Within this region, the gradient is bounded by G∗ + 2DL. Thus, the minimum number of iterations needed are 2D−Dα(G∗+2DL) . However,
  2D −D α(G∗ + 2DL) = 1 α(G ∗
  D + 2L)
  ≥ 1 α(G∗ 1−4κnαL2κnαG∗ + 2L)
  (Using the fact that D ≥ 2κnαG ∗
  1−4κnαL .)
  = 1
  α( 1−4κnαL2κnα + 2L)
  = 2κn
  ≥ 2n. Thus, the minimum number iterations needed to go outside the bound |xki − x∗| ≥ 2D is more than the length of the epoch. This implies that within the epoch, ‖xki − x∗‖ ≤ 2D, which proves a). We prove b) next:
  ‖xkn − x∗‖ = ∥∥∥∥∥ ( xk0 − α n∑ i=0 ∇fσki (x k i ) ) − x∗ ∥∥∥∥∥
  = ∥∥∥∥∥xk0 − x∗ − α n∑ i=0 ∇fσki (x k 0) + α n∑ i=0 (∇fσki (x k 0)−∇fσki (x k i )) ∥∥∥∥∥ Note that ∑n i=0∇fσki (x k 0) is just the sum of all component gradients at x
  k 0 , that is∑n
  i=0∇fσki (x k 0) = n∇F (xk0). Using this, we get
  ‖xkn − x∗‖ = ∥∥∥∥∥xk0 − x∗ − nα∇F (xk0) + α n∑ i=0 (∇fσki (x k 0)−∇fσki (x k i )) ∥∥∥∥∥ ≤ ∥∥xk0 − x∗ − nα∇F (xk0)∥∥+ α n∑
  i=0 ∥∥∥∇fσki (xk0)−∇fσki (xki )∥∥∥ (Triangle inequality.)
  ≤ ∥∥xk0 − x∗ − nα∇F (xk0)∥∥+ αL n∑
  i=0 ∥∥xk0 − xki ∥∥ , (12) where we used gradient Lipschitzness (Assumption 2) in the last step.
  To bound the first term above, we use the standard analysis of gradient descent on smooth, strongly convex functions as follows∥∥xk0 − x∗ − nα∇F (xk0)∥∥2 = ∥∥xk0 − x∗∥∥2 − 2nα〈xk0 − x∗,∇F (xk0)〉+ n2α2‖∇F (xk0)‖2
  ≤ ∥∥xk0 − x∗∥∥2 − 2nαµ‖xk0 − x∗‖2 + n2α2‖∇F (xk0)‖2
  (Using Ineq. (3)) = (1− nαµ) ∥∥xk0 − x∗∥∥2 + nα(nα‖∇F (xk0)‖2 − µ‖xk0 − x∗‖2)
  ≤ (1− nαµ) ∥∥xk0 − x∗∥∥2 + nα(nαL2‖xk0 − x∗‖2 − µ‖xk0 − x∗‖2)
  (Using gradient Lipschitzness) = (1− nαµ) ∥∥xk0 − x∗∥∥2 + nα(nαL2 − µ)‖xk0 − x∗‖2
  ≤ (1− nαµ) ∥∥xk0 − x∗∥∥2 ,
  where in the last step we used that α ≤ µnL2 since α ≤ 1 8κnL . Substituting this inequality in Ineq. (12), we get
  ‖xkn − x∗‖ ≤ √ 1− nαµ ∥∥xk0 − x∗∥∥+ αL n∑
  i=0 ∥∥xk0 − xki ∥∥ ≤ (
  1− 1 2 nαµ )∥∥xk0 − x∗∥∥+ αL n∑ i=0 ∥∥xk0 − xki ∥∥ . We have already proven a) that says that the iterates xki satisfy ‖xki − x∗‖ ≤ 2D. Using gradient Lipschitzness, this implies that the gradient norms stay bounded byG∗+2DL. Hence,
  ∥∥xk0 − xki ∥∥ ≤ αi(G∗ + 2DL). Using this,
  ‖xkn − x∗‖ ≤ (
  1− 1 2 nαµ )∥∥xk0 − x∗∥∥+ αL n∑ i=0 αi(G∗ + 2DL)
  ≤ (
  1− 1 2 nαµ
  ) D + αL n∑ i=0 αi(G∗ + 2DL)
  ≤ (
  1− 1 2 nαµ
  ) D + n2α2L(G∗ + 2DL)
  ≤ D,
  where we used the fact that D ≥ 2κnαG ∗
  1−4κnαL in the last step.
  B.5 PROOF OF LEMMA 2
  Without loss of generality, let σ = (1, 2, 3, . . . , n). This is only done for ease of notation. The analysis goes through for any other permutation σ too.
  First we show the lower bound. WLOG assume y0 > x0. Because α < 1/L, we have that ∀i, yi > xi by induction (see the equations below). Then,
  yi − xi = yi−1 − xi−1 − α(∇fi(yi−1)−∇fi(xi−1)) ≥ yi−1 − xi−1 − αL(yi−1 − xi−1) = (1− αL)(yi−1 − xi−1)
  ...
  = (1− αL)i(y0 − x0) ≥ (1− iαL)(y0 − x0). (13)
  Next we show the upper bound yn − xn = y0 − x0 − α n∑ i=1 (∇fi(yi−1)−∇fi(xi−1))
  = y0 − x0 − α n∑ i=1 (∇fi(y0)−∇fi(x0)) + α n∑ i=1 (∇fi(y0)−∇fi(x0)−∇fi(yi−1) +∇fi(xi−1))
  = y0 − x0 − nα(∇F (y0)−∇F (x0)) + α n∑ i=1 (∇fi(y0)−∇fi(x0)−∇fi(yi−1) +∇fi(xi−1))
  ≤ (1− nαµ)(y0 − x0) + α n∑ i=1 (∇fi(y0)−∇fi(x0)−∇fi(yi−1) +∇fi(xi−1)).
  (Using strong convexity)
  We use the fact that the function is twice differentiable: yn − xn = (1− nαµ)(y0 − x0) + α n∑ i=1 (∫ y0 x0 ∇2fi(t)dt− ∫ yi−1 xi−1 ∇2fi(t)dt ) = (1− nαµ)(y0 − x0)
  + α n∑ i=1 (∫ y0 x0 ∇2fi(t)dt− ∫ xi−1+(y0−x0) xi−1 ∇2fi(t)dt− ∫ yi−1 xi−1+(y0−x0) ∇2fi(t)dt ) = (1− nαµ)(y0 − x0)
  + α n∑ i=1 (∫ y0 x0 (∇2fi(t)−∇2fi(xi−1 − x0 + t))dt− ∫ yi−1 xi−1+(y0−x0) ∇2fi(t)dt ) .
  In the above, we used the convention that ∫ b a f(x)dx is the same as − ∫ a b f(x)dx if a > b.
  Now, we can use the Hessian Lipschitzness to bound the term as follows
  yn − xn ≤ (1− nαµ)(y0 − x0) + α n∑ i=1 (∫ y0 x0 LH |xi−1 − x0|dt− ∫ yi−1 xi−1+(y0−x0) ∇2fi(t)dt )
  ≤ (1− nαµ)(y0 − x0) + α n∑ i=1 (∫ y0 x0 LHGαndt− ∫ yi−1 xi−1+(y0−x0) ∇2fi(t)dt )
  = (1− nαµ)(y0 − x0) + LHGα2n2(x0 − y0)− α n∑ i=1 ∫ yi−1 xi−1+(y0−x0) ∇2fi(t)dt
  ≤ (1− nαµ)(y0 − x0) + LHGα2n2(x0 − y0) + α n∑ i=1 L((y0 − x0)− (yi − xi))
  ≤ (1− nαµ)(y0 − x0) + LHGα2n2(x0 − y0) + α n∑ i=1 L(iαL)(y0 − x0)
  (Using Ineq. (13).)
  ≤ (1− nαµ)(y0 − x0) + LHGα2n2(x0 − y0) + α2n2L2(y0 − x0).
  Thus, if we have α ≤ µ2n(L2+LHG) , then
  yn − xn ≤ (
  1− µnα 2
  ) (y0 − x0).
  C PROOF OF THEOREM 2
  To prove this theorem, we consider three step-size ranges and do a case by case analysis for each of them. We construct functions for each range such that the convergence of any permutation-based algorithm is “slow” for the functions on their corresponding step-size regime. The final lower bound is the minimum among the lower bounds obtained for the three regimes.
  In this proof, we will use different notation from the rest of the paper because we work with scalars in the proof, and hence the superscript will denote the scalar power, and not the epoch number. We will use xk,i to denote the i-th iterate of the the k-th epoch.
  We will construct three functions F1(x), F2(y), and F3(z), each the means of n component functions, such that
  • Any permutation-based algorithm on F1(x) with α ∈ [ 12(n−1)KL , 3 nL ] and initialization
  x1,0 = 0 results in
  ‖xK,n‖2 = Ω ( 1
  n3K2
  ) .
  F1 will be an n-dimensional function, that is x ∈ Rn. This function will have minimizer at 0 ∈ Rn. NOTE: This is the ‘key’ step-size range, and the proof sketch explained in the main paper corresponds to this function’s construction.
  • Any permutation-based algorithm on F2(y) with α ∈ [ 3nL , 1 L ] and initialization y1,0 = 0
  results in
  ‖yK,n‖2 = Ω ( 1
  n2 ) F2 will be an n-dimensional function, that is y ∈ Rn. This function will have minimizer at 0 ∈ Rn. The construction for this is also inspired by the construction for F1, but this is constructed significantly differently due to the different step-size range.
  • Any permutation-based algorithm on F3(z) with α /∈ [ 12(n−1)KL , 1 L ] and initialization
  z1,0 = 1 results in
  |zK,n|2 = Ω (1)
  F3 will be an 1-dimensional function, that is z ∈ R. This function will have minimizer at 0.
  Then, the 2n+ 1-dimensional function F ([x>,y>, z]>) = F1(x) + F2(y) + F3(z) will show bad convergence in any step-size regime. This function will have minimizer at 0 ∈ R2n+1. Furthermore,
  n− 1 n LI ∇2F1,∇2F2,∇2F3,∇2F 2LI,
  that is, F1, F2, F3, and F are all n−1n L-strongly convex and 2L-smooth.
  In the subsequent subsections, we prove the lower bounds for F1, F2, and F3 separately.
  NOTE: In Appendix C.4, we have discussed how the lower bound is actually invariant to initialization.
  C.1 LOWER BOUND FOR F1 , α ∈ [ 12(n−1)KL , 3 nL ]
  We will work in n-dimensions (n is even) and represent a vector in the space as z = [x1, y1, . . . , xn/2, yn/2]. These xi and yi are not related to the vectors used by F2 and F3 later, we only use xi and yi to make the proof for this part easier to understand.
  We start off by defining the n component functions: For i ∈ [n/2], define
  fi(z) := L
  2 x2i − xi + yi + ∑ j 6=i ( L 2 x2j + L 2 y2j ) , and
  gi(z) := L
  2 y2i − yi + xi + ∑ j 6=i ( L 2 x2j + L 2 y2j ) .
  Thus, F (z) := 1n (∑n/2 i=1 fi + ∑n/2 i=1 gi ) = ( n−1 n ) L 2 ‖z‖ 2. This function has minimizer at z∗ = 0.
  Let zk,j denote z at the j-th iteration in k-th epoch. We initialize at z1,0 = 0. For the majority of this proof, we will work inside a given epoch, so we will skip the subscript denoting the epoch. We denote the j-th iterate within the epoch as zj and the coordinates as zj = [xj,1, yj,1, . . . , xj,n/2, yj,n/2]
  In the current epoch, let σ be the permutation of {f1, . . . , fn/2, g1, . . . , gn/2} used.
  For any i ∈ [n/2], let p and q be indices such that σp = fi and σq = gi. Let us consider the case that p < q (the case when p > q will be analogous). Then, it can be seen that
  xn,i = (1− αL)n−1x0,i + (1− αL)n−p−1α− (1− αL)n−qα, and yn,i = (1− αL)n−1y0,i − (1− αL)n−pα+ (1− αL)n−qα. (14)
  Hence,
  xn,i + yn,i = (1− αL)n−1(x0,i + y0,i) + 2(1− αL)n−p−1α2L ≥ (1− αL)n−1(x0,i + y0,i) + 2(1− αL)n−1α2L.
  In the other case, when p > q, we will get the same inequality. Let xK,n,i and yK,n,i denote the value of xn,i and yn,i at the K-th epoch. Then, recalling that z was initialized to 0, we use the inequality above to get
  xK,n,i + yK,n,i ≥ (1− αL)(n−1)K · 0 + 2(1− αL)n−1α2L 1− (1− αL)(n−1)K
  1− (1− αL)n−1
  = 2(1− αL)n−1α2L1− (1− αL) (n−1)K
  1− (1− αL)n−1 (15)
  Since this inequality is valid for all i, we get that
  ‖zK,n‖2 = n/2∑ i=1 (x2K,n,i + y 2 K,n,i)
  ≥ n/2∑ i=1 1 2 (xK,n,i + yK,n,i) 2 ≥ n/2∑ i=1 1 2 ( 2(1− αL)n−1α2L1− (1− αL) (n−1)K 1− (1− αL)n−1 )2
  = n ( (1− αL)n−1α2L1− (1− αL) (n−1)K
  1− (1− αL)n−1
  )2 . (16)
  Note that if α ≤ 3nL and n ≥ 4, then (1− αL) n−1 ≥ 18 . Using this in (16), we get
  ‖zK,n‖2 ≥ n ( (1− αL)n−1α2L1− (1− αL) (n−1)K
  1− (1− αL)n−1
  )2
  ≥ n 128
  α4L2 ( 1− (1− αL)(n−1)K
  1− (1− αL)n−1 )2 = n
  128L2
  ( (αL)2
  1− (1− αL)n−1
  )2 ( 1− (1− αL)(n−1)K )2
  We consider two cases:
  1. Case A, αL ≤ 12(n+2) : It can be verified that (αL)2
  1−(1−αL)n−1 is an increasing function of α when αL ≤ 12(n+2) . Noting that we are working in the range when α ≥ 1 2(n−1)KL , then
  ‖zK,n‖2 ≥ n
  128L2
   ( 1 2(n−1)K )2 1− ( 1− 12(n−1)K )n−1  2 ( 1− (1− αL)(n−1)K )2
  ≥ n 128L2
   ( 1 2(n−1)K )2 1
  2(n−1)K (n− 1)
   2 ( 1− (1− αL)(n−1)K )2
  = n 512(n− 1)4K2L2 ( 1− (1− αL)(n−1)K )2
  ≥ n 512(n− 1)4K2L2
  ( 1− e−αL(n−1)K )2 ≥ n 512(n− 1)4K2L2 ( 1− e−1/2 )2
  = Ω
  ( 1
  n3K2L2
  ) .
  2. Case B, αL > 12(n+2) : In this case,
  ‖zK,n‖2 ≥ n
  128L2
   ( 1 2(n+2) )2 1  2 ( 1− (1− αL)(n−1)K )2
  ≥ n 128L2
   ( 1 2(n+2) )2 1  2 ( 1− e−αL(n−1)K )2
  ≥ n 128L2
   ( 1 2(n+2) )2 1  2 ( 1− e− 3(n−1)K 2(n+2)
  )2 = Ω ( 1
  n3L2
  ) .
  C.2 LOWER BOUND FOR F2 , α ∈ [ 3nL , 1 L ]
  We will work in n-dimensions and represent a vector in the space as y = [y1, . . . , yn].
  We start off by defining the n component functions: For i ∈ [n], define
  fi(y) := −yi + ∑ j 6=i
  ( Ly2j
  2 + yj n− 1
  )
  Thus, F (y) := 1n ∑n i=1 fi = ( n−1 n ) L 2 ‖y‖ 2. This function has minimizer at y∗ = 0.
  Let yk,j denote y at the j-th iteration in k-th epoch. We initialize at y1,0 = 0. For the majority of this proof, we will work inside a given epoch. We denote the j-th iterate within the epoch as yj and the coordinates as yj = [yj,1, . . . , yj,n]
  In the current epoch, let σ be the permutation of {f1, . . . , fn} used. Let i be the index such that σn = i, that is, i is the last element of the permutation σ. Then at the end of the epoch,
  yn,i = (1− αL)n−1y0,i + α− α
  n− 1 n−2∑ j=0 (1− αL)j
  = (1− αL)n−1y0,i + α− (1− (1− αL)n−1)
  (n− 1)L . (17)
  For some j ∈ [n], let s be the integer such that σs = j, that is j is the s-th element in the permutation σ. Then for any j and any epoch,
  yn,j = (1− αL)n−1y0,j + α(1− αL)n−s − α
  n− 1 n−2∑ j=0 (1− αL)j .
  Then,
  yn,j ≥ (1− αL)n−1y0,j − α
  n− 1 n−2∑ j=0 (1− αL)j
  = (1− αL)n−1y0,j − (1− (1− αL)n−1)
  (n− 1)L .
  Note that the above is independent of σ, and hence applicable for all epochs. Applying it recursively and noting that we initialized y1,0 = 0, we get
  yn,j ≥ − (1− (1− αL)n−1)
  (n− 1)L
  K∑ k=1 (1− αL)n
  ≥ − 1 (n− 1)L .
  Note that y0,i is just yn,i from the previous epoch. Hence we can substitute the inequality above for y0,i in (17). Thus,
  yn,i ≥ − (1− αL)n−1 (n− 1)L + α− (1− (1− αL) n−1) (n− 1)L
  = α− 1 (n− 1)L
  ≥ 3 nL − 1 (n− 1)L
  = Ω
  ( 1
  nL
  ) .
  This gives us that ‖yk,n‖2 = Ω ( 1 n2L2 ) for any k.
  C.3 LOWER BOUND FOR F3 , α /∈ [ 12(n−1)KL , 1 L ]
  Consider the function F (z) = 1n ∑n i=1 fi(z), where fi(z) = Lz
  2 for all i. Note that the gradient of any function at z is just −2Lz. Hence, regardless of the permutation, if we start the shuffling based gradient descent method at z1,0 = 1, we get that
  zK,n = (1− 2αL)nKz1,0 = (1− 2αL)nK .
  In the case when α ≤ 12(n−1)KL , we see that
  zK,n ≥ (
  1− 2 1 2(n− 1)KL L )nK ≥ (
  1− 1 (n− 1)K )nK = Ω(1),
  for n,K ≥ 2. Finally, in the case when α ≥ 1L , we see that
  |zK,n| = |1− 2αL|nK ≥ ∣∣∣∣1− 2 1LL ∣∣∣∣nK ≥ 1nK
  = Ω(1).
  C.4 DISCUSSION ABOUT INITIALIZATION
  The lower bound partitions the step size in three ranges -
  • In the step size ranges α ∈ [ 12(n−1)KL , 3 nL ] and α ∈ [ 3 nL , 1 L ], the initializations are done at
  the minimizer and it is shown that any permutation-based algorithm will still move away from the minimizer. The choice of initializing at the minimizer was solely for the convenience of analysis and calculations, and the proof should work for any other initialization as well. Furthermore, the effect of initializing at any arbitrary (non-zero) point will decay exponentially fast with epochs anyway. To see how, note that every epoch can be treated as n steps of full gradient descent and some noise, and hence the full gradient descent part will essentially keep decreasing the effect of initialization exponentially, and what we would be left with is the noise in each epoch. Thus, it was more convenient for us to just assume initialization at the minimizer and only focus on the noise in each epoch.
  • The step size range α /∈ [ 12(n−1)KL , 1 L ] can be divided into two parts, α ∈ [0, 1 2(n−1)KL ]
  and α ∈ [ 1L ,∞). For the range α ∈ [0, 12(n−1)KL ], we essentially show that the step size is too small to make any meaningful progress towards the minimizer. Hence, instead of initializing at 1, initializing at any other arbitrary (non-zero) point will also give the same slow convergence rate. For the range α ∈ [ 1L ,∞), we show that the optimization algorithm will in fact diverge since the step size is too large. Hence, even here, any other arbitrary (non-zero) point of initialize will also give divergence.
  D PROOF OF THEOREM 3
  Define f1(x) := L2 x 2 − x and f2(x, y) := −L4 x 2 + x. Thus, F (x, y) = L8 x 2. This function has minimizer at x∗ = 0. For this proof, we will use the convention that xi,j is the iterate after the i-th iteration of the j-th epoch. Further, the number in the superscript will be the scalar power. For example x2i,j = xi,j · xi,j .
  Initialize at x0,0 = 1L . Then at epoch k, there are two possible permutations: σ = (1, 2) and σ = (2, 1). In the first case , σ = (1, 2), we get that after the first iteration of the epoch,
  x1,k = x0,k − α∇f1(x0,k, y0,k) = (1− αL)x0,k + α,
  Continuing on, in the second iteration, we get
  x2,k = x1,k − α∇f2(x1,k, y1,k)
  = ( 1 + 1
  2 αL
  ) x1,k − α
  = ( 1 + 1
  2 αL
  ) ((1− αL)x0,k + α)− α
  = ( 1 + 1
  2 αL
  ) (1− αL)x0,k + 1
  2 α2L.
  Note that x0,k+1 = x2,k. Thus, x0,k+1 = ( 1 + 12αL ) (1− αL)x0,k + 12α 2L.
  Similarly, for the other possible permutation, σ = (2, 1), we get x0,k+1 = ( 1 + 12αL ) (1−αL)x0,k+ α2L. Thus, regardless of what permutation we use, we get that
  x0,k+1 ≥ ( 1 + 1
  2 αL
  ) (1− αL)x0,k + 1
  2 α2L ≥ (1− αL)x0,k +
  1 2 α2L.
  Hence, recalling that we initialized at x0,0 = 1L , we get
  xn,K = x0,K+1
  ≥ (1− αL)K 1 L + 1 2 α2L K−1∑ i=0 (1− αL)i
  = 1 L (1− αL)K + 1 2 α2L
  1− (1− αL)K
  1− (1− αL)
  = 1 L (1− αL)K + 1 2 α ( 1− (1− αL)K ) . (18)
  Now, if α ≥ 1KL , then 1 2 α ( 1− (1− αL)K ) ≥ α 8 = Ω
  ( 1
  KL
  ) ,
  and hence, x2n,K = Ω( 1 K2L2 ). Otherwise, if α ≤ 1 KL , then continuing on from (18),
  xn,K ≥ 1 L (1− αL)K + 1 2 α ( 1− (1− αL)K ) ≥ 1 L (1− αL)K
  ≥ 1 L e−αLK ≥ 1 L e−1
  = Ω(1/L),
  and hence in this case, x2n,K = Ω( 1 L2 ).
  E PROOF OF THEOREM 4 Let F (x) := 1n ∑n i=1 fi(x) be such that its minimizer it at the origin. This can be assumed without loss of generality because we can shift the coordinates appropriately, similar to (Safran & Shamir, 2019). Since the fi are convex quadratics, we can write them as fi(x) = 12x
  >Aix − b>i x + ci, where Ai are symmetric, positive-semidefinite matrices. We can omit the constants ci because they do not affect the minimizer or the gradients. Because we assume that the minimizer of F (x) is at the origin, we get that
  n∑ i=1 bi = 0. (19)
  Let σ = (σ1, σ2, . . . , σn) be the random permutation of (1, 2, . . . , n) generated at the beginning of the algorithm. Then for k ∈ (1, 2, . . . ,K/2), epoch 2k − 1 sees the n functions in the following sequence: (
  1 2 x>Aσ1x− b>σ1x, 1 2 x>Aσ2x− b>σ2x, . . . , 1 2 x>Aσnx− b>σnx
  ) ,
  whereas epoch 2k sees the n functions in the reverse sequence:( 1
  2 x>Aσnx− b>σnx,
  1 2 x>Aσn−1x− b>σn−1x, . . . , 1 2 x>Aσ1x− b>σ1x
  ) .
  We define Si := αAσi and ti = αbσi for convenience of notation. We start off by computing the progress made during an even indexed epoch. Since the even epochs use the reverse permutation, we get
  x2k+10 = x 2k n = x2kn−1 − α ( Aσ1x 2k n−1 − bσ1 ) (fσ1 is used at the last iteration of even epochs.)
  = (I − αAσ1)x2kn−1 + αbσ1 = (I − S1)x2kn−1 + t1.
  We recursively apply the same procedure as above to the whole epoch to get the following
  x2k+10 = (I − S1)x2kn−1 + t1 = (I − S1) ( (I − S2)x2kn−2 + t2 ) + t1
  = (I − S1)(I − S2)x2kn−2 + (I − S1)t2 + t1 ...
  = ( n∏ i=1 (I − Si) ) x2k0 + n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i, (20) where the product of matrices {Mi} is defined as ∏m i=lMi = MlMl+1 . . .Mm if m ≥ l and 1 otherwise. Similar to Eq. (20), we can compute the progress made during an odd indexed epoch. Recall that the only difference will be that the odd indexed epochs see the permutations in the order (σ1, σ2, . . . , σn) instead of (σn, σn−1, . . . , σ1). After doing the computation, we get the following equation:
  x2k0 = ( n∏ i=1 (I − Sn−i+1) ) x2k−10 + n∑ i=1 n−i∏ j=1 (I − Sn−j+1)  ti. Combining the results above, we can get the total progress made after the pair of epoch 2k − 1 and 2k:
  x2k+10 = ( n∏ i=1 (I − Si) ) x2k0 + n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i =
  ( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ) x2k−10 +  n∏ j=1 (I − Sj)  n∑ i=1 n−i∏ j=1 (I − Sn+1−j)  ti +
  n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i. (21) In the sum above, the first term will have an exponential decay, hence we need to control the next two terms. We denote the sum of the terms as z (see the definition below) and we will control its norm
  later in this proof.
  z :=  n∏ j=1 (I − Sj)  n∑ i=1 n−i∏ j=1 (I − Sn+1−j)  ti + n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i =
  n∑ i=1  n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  ti + n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i. To see where the iterates end up after K epochs, we simply set 2k = K in Eq. 21 and then keep applying the equation recursively to preceding epochs. Then, we get
  xKn = x K+1 0 = ( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ) xK−10 + z
  = (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))2 xK−30 + ( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ) z + z
  ...
  = (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10 + K 2 −1∑ k=0 (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))k z.
  Taking squared norms and expectations on both sides, we get
  E[‖xKn ‖2] = E ∥∥∥∥∥∥ (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10
  + K 2 −1∑ k=0 (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))k z ∥∥∥∥∥∥ 2 
  ≤ 2E  ∥∥∥∥∥∥ (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10 ∥∥∥∥∥∥ 2 
  + 2E  ∥∥∥∥∥∥ K 2 −1∑ k=0 (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))k z ∥∥∥∥∥∥ 2 
  (Since (a+ b)2 ≤ 2a2 + 2b2)
  ≤ 2E  ∥∥∥∥∥∥ (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10 ∥∥∥∥∥∥ 2 
  + 2E  ‖z‖ K2 −1∑
  k=0 ∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1) ∥∥∥∥∥ k 2  .
  We assumed that the functions fi have L-Lipschitz gradients (Assumption 2). This translates to Ai having maximum eigenvalue less than L. Hence, if α ≤ 1/L, we get that I − αAi is positive semi-definite with maximum eigenvalue bounded by 1. Hence, ‖I − Si‖ ≤ 1. Using this and the fact that for matrices M1 and M2, ‖M1M2‖ ≤ ‖M1‖‖M2‖, we get that E [∥∥xKn ∥∥2] ≤ 2E  ∥∥∥∥∥∥ (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10 ∥∥∥∥∥∥ 2 
  + 2E  ‖z‖ K2 −1∑
  k=0 ( n∏ i=1 ‖I − Si‖ n∏ i=1 ‖I − Sn−i+1‖ )k2 
  ≤ 2E  ∥∥∥∥∥∥ (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10 ∥∥∥∥∥∥ 2 + 2E  ‖z‖ K2 −1∑ k=0 1 2 
  = 2E  ∥∥∥∥∥∥ (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))K/2 x10 ∥∥∥∥∥∥ 2 + K2 2 E [ ‖z‖2 ] .
  We handle the two terms above separately. For the first term, we have the following bound: Lemma 3. If α ≤ 18κL min { 2, √ κ n } , then∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1)
  ∥∥∥∥∥ ≤ 1− αnµ Note that K ≥ 80κ3/2 log(nK) max { 1, √ κ n } =⇒ α ≤ 18κL min { 2, √ κ n } .
  We also have the following lemma that bounds the expected squared norm of z.
  Lemma 4. If α ≤ 1L , then E [ ‖z‖2 ] ≤ 2n2α4L2(G∗)2 + 170n5α6L4G2 log n,
  where G∗ = maxi ‖bi‖, and the expectation is taken over the randomness of z.
  Using these lemmas, we get that
  E [∥∥xKn ∥∥2] ≤ 2 (1− nαµ)K/2 ∥∥x10∥∥2 +K2n2α4L2G2 + 85K2n5α6L4G2 log n
  ≤ 2e− 12nαµK ∥∥x10∥∥2 +K2n2α4L2G2 + 85K2n5α6L4G2 log n.
  Substituting α = 10 lognKµnK gives us the result.
  E.1 PROOF OF LEMMA 3
  We define (S̃1, . . . , S̃n) := (S1, . . . ,Sn) and (S̃n+1, . . . , S̃2n) := (Sn, . . . ,S1). Then,∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1) ∥∥∥∥∥ = ∥∥∥∥∥ 2n∏ i=1 (I − S̃i) ∥∥∥∥∥ =
  ∥∥∥∥∥∥I − 2n∑ i=1 S̃i + ∑ i<j S̃iS̃j − . . . ∥∥∥∥∥∥ ≤
  ∥∥∥∥∥∥I − 2n∑ i=1 S̃i + ∑ i<j S̃iS̃j ∥∥∥∥∥∥+ ∥∥∥∥∥∥ ∑ i<j<k S̃iS̃jS̃k ∥∥∥∥∥∥+ . . . ≤
  ∥∥∥∥∥∥I − 2n∑ i=1 S̃i + ∑ i<j S̃iS̃j ∥∥∥∥∥∥+ ∑ i<j<k ∥∥∥S̃i∥∥∥∥∥∥S̃j∥∥∥∥∥∥S̃k∥∥∥+ . . . .
  Note that ∑ i<j S̃iS̃j = 2 ∑ i 6=j SiSj + n∑ i=1 S2i
  = 2 ( n∑ i=1 Si )( n∑ i=1 Si ) − n∑ i=1 S2i .
  Substituting this and noting that ∑2n i=1 S̃i = 2 ∑n i=1 Si, we get∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1) ∥∥∥∥∥ ≤ ∥∥∥∥∥I − 2 n∑ i=1 Si + 2 ( n∑ i=1 Si )( n∑ i=1 Si ) − n∑ i=1 S2i
  ∥∥∥∥∥ + ∑ i<j<k
  ∥∥∥S̃i∥∥∥ ∥∥∥S̃j∥∥∥ ∥∥∥S̃k∥∥∥+ . . . ≤
  ∥∥∥∥∥I − 2 n∑ i=1 Si + 2 ( n∑ i=1 Si )( n∑ i=1 Si )∥∥∥∥∥+ ∥∥∥∥∥ n∑ i=1 S2i ∥∥∥∥∥ + ∑ i<j<k
  ∥∥∥S̃i∥∥∥ ∥∥∥S̃j∥∥∥ ∥∥∥S̃k∥∥∥+ . . . . Let us denote T := ∑n i=1 Si. Then, we know by Assumptions 2 and 3 that T has eigenvalues in [nαµ, nαL]. Then, as long as α ≤ 14nL , we get that∥∥∥∥∥I − 2 n∑ i=1 Si + 2 ( n∑ i=1 Si )( n∑ i=1 Si
  )∥∥∥∥∥ = ∥∥I − 2T + 2T 2∥∥ ≤ 1− 3
  2 nαµ.
  Substituting this, we get∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1) ∥∥∥∥∥ ≤ ( 1− 3 2 nαµ ) + ∥∥∥∥∥ n∑ i=1 S2i ∥∥∥∥∥+ ∑ i<j<k ∥∥∥S̃i∥∥∥∥∥∥S̃j∥∥∥∥∥∥S̃k∥∥∥+ . . . . By Assumption 2, we know that ‖Ai‖ ≤ L. Hence, ‖S̃i‖ ≤ αL. Hence,∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1) ∥∥∥∥∥ ≤ ( 1− 3 2 nαµ ) + nα2L2 + (( 2n 3 ) α3L3 + ( 2n 4 ) α4L4 + . . . )
  ≤ 1− 3 2 nαµ+ nα2L2 + 2n∑ i=3 (2nαL)i
  ≤ 1− 3 2 nαµ+ nα2L2 +
  8n3α3L3
  1− 2nαL
  ≤ 1− 3 2 nαµ+ nα2L2 + 16n3α3L3. (Since α ≤ 1/4nL.)
  Finally, as long as α ≤ 14κL and α ≤ 1 8nL √ κ ,∥∥∥∥∥ n∏ i=1 (I − Si) n∏ i=1 (I − Sn−i+1) ∥∥∥∥∥ ≤ 1− nαµ. E.2 PROOF OF LEMMA 4
  We start off by computing the first order expansion of z. We have the following lemma for this:
  Lemma 5.
  z = n∑ i=1 Siti+ 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i ) ,
  where (S̃1, . . . , S̃n) := (S1, . . . ,Sn) and (S̃n+1, . . . , S̃2n) := (Sn, . . . ,S1).
  The proof of this lemma is quite algebraic and hence has been pushed to the end, in Appendix E.4. The strategy is to bound ‖Siti‖, ‖ ∑2n−j i=1 ti‖, and ‖ ∑n−j i=1 tn+1−i‖. Hence, we apply Lemma 5 and use triangle inequality:
  E [ ‖z‖2 ] = E ∥∥∥∥∥∥ n∑ i=1 Siti + 2n−1∑ j=n+1 j−1∑ k=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti )
  + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )∥∥∥∥∥∥ 2 
  ≤ E  n∑ i=1 ‖Si‖‖ti‖+ 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 ‖I − S̃p‖ ) ‖S̃l‖‖S̃j‖ ( 2n−j∑ i=1 ‖ti‖ )
  + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 ‖I − Sp‖ ) ‖Sl‖‖Sj‖ ( n−j∑ i=1 ‖tn+1−i‖ )2  .
  Now, we recall that ‖Si‖ ≤ αL and ‖ti‖ ≤ αG. Because α ≤ 1/L, we also get that ‖I − Si‖ ≤ 1. Using these,
  E[‖z‖2] ≤ E  n∑ i=1 α2LG+ 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 1 ) αLαL ∥∥∥∥∥ 2n−j∑ i=1 ti ∥∥∥∥∥ +
  n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 1 ) αLαL ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ 2 
  ≤ E  nα2LG+ 2nα2L2 2n−1∑
  j=n+1 ∥∥∥∥∥ 2n−j∑ i=1 ti ∥∥∥∥∥+ nα2L2 n−1∑ j=1 ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ 2 
  = n2α4L2E  G+ 2L 2n−1∑
  j=n+1 ∥∥∥∥∥ 2n−j∑ i=1 ti ∥∥∥∥∥+ L n−1∑ j=1 ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ 2 
  ≤ 2n2α4L2 G2 + L2E  2 2n−1∑
  j=n+1 ∥∥∥∥∥ 2n−j∑ i=1 ti ∥∥∥∥∥+ n−1∑ j=1 ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ 2   .
  (Since (a+ b)2 ≤ 2a2 + 2b2)
  Using Hoeffding-Serfling inequality for bounded random vectors (Schneider, 2016, Theorem 2), we get the following lemma Lemma 6. For all j, l ∈ [1, n] we have that
  E ∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ 2  ≤ 18jα2(G∗)2 log(n)
  E [∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ ∥∥∥∥∥ l∑ i=1 ti ∥∥∥∥∥ ] ≤ 18 √ jlα2(G∗)2 log(n),
  where G∗ = maxi ‖bi‖, and the expectation is taken over the randomness of ti.
  Writing out the expansion of ( 2 ∑2n−1 j=n+1 ∥∥∥∑2n−ji=1 ti∥∥∥+∑n−1j=1 ∥∥∥∑n−ji=1 tn+1−i∥∥∥)2 and using the lemma above on the individual terms, we get
  E [ ‖z‖2 ] ≤ 2n2α4L2G2 + 2n2α4L4(90α2n3G2 log n).
  E.3 PROOF OF LEMMA 6
  This proof is similar to the one in (Ahn et al., 2020). Define G∗ := maxi ‖bi‖. We use the following theorem adapted to our setting Theorem 7. [(Schneider, 2016, Theorem 2)] With probability at least 1− δn ,∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ ≤ αG∗ √ 8j ( 1− j − 1 n ) log 2n δ .
  Then taking a union bound over j = 1, . . . , n, we get that with probability at least 1− δ,
  ∀j ∈ [1, n] : ∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ ≤ αG∗ √ 8j ( 1− j − 1 n ) log 2n δ ≤ αG∗ √ 8j log 2n δ .
  Then, for the complementary event (which happens with probability δ), we use the fact that ‖ti‖ = ‖αbσi‖ ≤ αG∗ to get the following:
  ∀j ∈ [1, n] : ∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ ≤ j∑ i=1 ‖ti‖ ≤ αG∗j.
  Now, choose δ = 1/n. Then, we get that
  E ∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ 2  ≤ (1− 1 n ) 8jα2G2 log(2n2) + 1 n (αG∗j)2
  ≤ 18jα2(G∗)2 log n.
  Similarly, we can also get
  E [∥∥∥∥∥ j∑ i=1 ti ∥∥∥∥∥ ∥∥∥∥∥ l∑ i=1 ti ∥∥∥∥∥ ] ≤ 18 √ jlα2(G∗)2 log(n).
  E.4 PROOF OF LEMMA 5
  As in the lemma’s statement, define (S̃1, . . . , S̃2n) as (S̃1, . . . , S̃n) = (S1, . . . ,Sn) and (S̃n+1, . . . , S̃2n) = (Sn, . . . ,S1). As a reminder, the term z is defined as follows:
  z := n∑ i=1  n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  ti + n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i. (22) First, we analyze the first term in z. Towards that end, we start by expanding(∏n
  j=1(I − Sj) )(∏n−i j=1(I − Sn+1−j) )
  : n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  = 2n−i∏ j=1 (I − S̃j)  (23)
  = 2n−i−1∏ j=1 (I − S̃j) (I − S̃2n−i)
  = 2n−i−1∏ j=1 (I − S̃j) − 2n−i−1∏ j=1 (I − S̃j)  S̃2n−i. Similarly, we expand the term (∏2n−i−1 j=1 (I − S̃j) ) and then recursively keep doing it to get the following: n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  = 2n−i−1∏ j=1 (I − S̃j) − 2n−i−1∏ j=1 (I − S̃j)
   S̃2n−i =
  2n−i−2∏ j=1 (I − S̃j) − 2n−i−2∏ j=1 (I − S̃j)  S̃2n−i−1 − 2n−i−1∏ j=1 (I − S̃j)  S̃2n−i = ...
  = (I − S̃1)− (I − S̃1)S̃2 − . . .
  − 2n−i−2∏ j=1 (I − S̃j)  S̃2n−i−1 − 2n−i−1∏ j=1 (I − S̃j)  S̃2n−i = I −
  2n−i∑ j=1 ( j−1∏ l=1 (I − S̃l) ) S̃j .
  Note that the term ∏j−1 l=1 (I − S̃l) above is similar to the RHS of Eq. (23). Hence, we repeat the process again on this term to get the following n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  = I − 2n−i∑ j=1 ( I − j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃l ) S̃j
  = I − 2n−i∑ j=1 S̃j + 2n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j . (24)
  Using this in the first term ∑n i=1 (∏n j=1(I − Sj) )(∏n−i j=1(I − Sn+1−j) ) ti(in Eq. (22)), we get n∑ i=1  n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  ti = n∑ i=1 I − 2n−i∑ j=1 S̃j + 2n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j
   ti =
  n∑ i=1 ti − n∑ i=1 2n−i∑ j=1 S̃jti + n∑ i=1 2n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃jti.
  Now, we use the fact that ∑n i=1 bi = 0 (Eq. (19)) to get that ∑n i=1 ti = 0. Then, n∑ i=1  n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  ti = − n∑ i=1 2n−i∑ j=1 S̃jti + n∑ i=1 2n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃jti
  For convenience, we define M̃j := ∑j−1 l=1 (∏l−1 p=1(I − S̃p) ) S̃lS̃j . Then,
  n∑ i=1  n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  ti = − n∑ i=1 2n−i∑ j=1 S̃jti + n∑ i=1 2n−i∑ j=1 M̃jti
  = − n∑ i=1 2n−i∑ j=1 S̃jti + n∑ i=1 n∑ j=1 M̃jti + n∑ i=1 2n−i∑ j=n+1 M̃jti = − n∑ i=1 2n−i∑ j=1 S̃jti + n∑ j=1 M̃j n∑ i=1 ti + n∑ i=1 2n−i∑ j=n+1 M̃jti = − n∑ i=1 2n−i∑ j=1 S̃jti + n∑ i=1 2n−i∑ j=n+1 M̃jti. (Since ∑n i=1 ti = 0)
  Note that ∑n i=1 ∑2n−i j=n+1 M̃jti = ∑2n−1 j=n+1 M̃j (∑2n−j i=1 ti ) . Hence,
  n∑ i=1  n∏ j=1 (I − Sj) n−i∏ j=1 (I − Sn+1−j)  ti = − n∑ i=1 2n−i∑ j=1 S̃jti + 2n−1∑ j=n+1 M̃j ( 2n−j∑ i=1 ti )
  = − n∑ i=1 n∑ j=1 S̃jti − n∑ i=1 2n−i∑ j=n+1 S̃jti + 2n−1∑ j=n+1 M̃j ( 2n−j∑ i=1 ti ) = − n∑ j=1 S̃j n∑ i=1 ti − n∑ i=1 2n−i∑ j=n+1 S̃jti + 2n−1∑ j=n+1 M̃j ( 2n−j∑ i=1 ti ) = − n∑ i=1 2n−i∑ j=n+1 S̃jti + 2n−1∑ j=n+1 M̃j ( 2n−j∑ i=1 ti ) (Since ∑n i=1 ti = 0) = − n∑ i=1 2n−i∑ j=n+1 S̃jti + 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) = − n∑ i=1 n∑ j=i+1 Sjti + 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) . (25)
  Next we analyze the second term in z. For this, we start by expanding ∏n−i j=1(I − Sj) in a similar way as Eq. (24)
  n−i∏ j=1 (I − Sj) = I − n−i∑ j=1 Sj + n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj (26)
  Using this, we get n∑ i=1 n−i∏ j=1 (I − Sj) tn+1−i = n∑ i=1 I − n−i∑ j=1 Sj + n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj  tn+1−i =
  n∑ i=1 tn+1−i − n∑ i=1 n−i∑ j=1 Sjtn+1−i + n∑ i=1 n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSjtn+1−i
  = − n∑ i=1 n−i∑ j=1 Sjtn+1−i + n∑ i=1 n−i∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSjtn+1−i,
  where we used the fact that ∑n i=1 ti = 0 in the last equality. For convenience, we define Mj :=∑j−1
  l=1 (∏l−1 p=1(I − Sp) ) SlSj . Then,
  n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i = − n∑ i=1 n−i∑ j=1 Sjtn+1−i + n∑ i=1 n−i∑ j=1 Mjtn+1−i.
  Since ∑n i=1 ∑n−i j=1 Mjtn+1−i = ∑n−1 j=1 Mj (∑n−j i=1 tn+1−i ) , we get
  n∑ i=1 n−i∏ j=1 (I − Sj)  tn+1−i = − n∑ i=1 n−i∑ j=1 Sjtn+1−i + n−1∑ j=1 Mj ( n−j∑ i=1 tn+1−i )
  = − n∑ l=1 l−1∑ j=1 Sjtl + n−1∑ j=1 Mj ( n−j∑ i=1 tn+1−i ) = − n∑ i=1 i−1∑ j=1 Sjti + n−1∑ j=1 Mj ( n−j∑ i=1 tn+1−i ) = − n∑ i=1 i−1∑ j=1 Sjti + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i ) .
  (27)
  Finally, substituting Eq. (25) and (27) in the definition of z (Eq. (22)), we get z = − n∑ i=1 n∑ j=i+1 Sjti + 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti )
  − n∑ i=1 i−1∑ j=1 Sjti + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )
  = − n∑ i=1 i−1∑ j=1 Sj + n∑ j=i+1 Sj  ti +
  2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )
  = − n∑ i=1  n∑ j=1 Sj  ti + n∑ i=1 Siti
  + 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )
  = −  n∑ j=1 Sj  n∑ i=1 ti + n∑ i=1 Siti
  + 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )
  = n∑ i=1 Siti + 2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) + n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i ) ,
  where we used the fact that ∑n i=1 ti = 0 in the last equality.
  F PROOF OF THEOREM 5
  The proof is similar to that of Theorem 4, except for that here we leverage the independence of random permutations in every other epoch. The setup is also the same as Theorem 4, but we explain it again here nevertheless, for the completeness of this proof.
  Let F (x) := 1n ∑n i=1 fi(x) be such that its minimizer it at the origin. This can be assumed without loss of generality because we can shift the coordinates appropriately, similar to (Safran & Shamir, 2019). Since the fi are convex quadratics, we can write them as fi(x) = 12x
  >Aix − b>i x + ci, where Ai are symmetric, positive-semidefinite matrices. We can omit the constants ci because they do not affect the minimizer or the gradients. Because we assume that the minimizer of F (x) is at the origin, we get that
  n∑ i=1 bi = 0. (28)
  Let σk = (σk1 , σ k 2 , . . . , σ k n) be the random permutation of (1, 2, . . . , n) sampled in epoch 2k − 1. Then epoch 2k − 1 sees the n functions in the reverse sequence:( 1
  2 x>Aσk1x− b > σk1 x,
  1 2 x>Aσk2x− b > σk2 x, . . . , 1 2 x>Aσknx− b > σkn x
  ) ,
  whereas epoch 2k sees the n functions in the reverse sequence:( 1
  2 x>Aσknx− b > σkn
  x, 1
  2 x>Aσkn−1x− b > σkn−1
  x, . . . , 1
  2 x>Aσk1x− b > σk1 x
  ) .
  We define Ski := αAσki and t k i = αbσki for convenience of notation. We start off by computing the progress made duting an even indexed epoch. Since the even epochs use the reverse permutation of σk, we get
  x2k+10 = x 2k n = x2kn−1 − α ( Aσk1x 2k n−1 − bσk1 ) (fσk1 used at the last iteration of epoch 2k.)
  = ( I − αAσk1 ) x2kn−1 + αbσk1
  = (I − Sk1 )x2kn−1 + tk1 .
  We recursively apply the same procedure as above to the whole epoch to get the following
  x2k+10 = (I − Sk1 )x2kn−1 + tk1 = (I − Sk1 ) ( (I − S2)x2kn−2 + tk2 ) + tk1
  = (I − Sk1 )(I − Sk2 )x2kn−2 + (I − S1)tk2 + tk1 ...
  = ( n∏ i=1 (I − Ski ) ) x2k0 + n∑ i=1 n−i∏ j=1 (I − Skj )  tkn+1−i, (29) where the product of matrices {Mi} is defined as ∏m i=lMi = MlMl+1 . . .Mm if m ≥ l and 1 otherwise. Similar to Eq. (29), we can compute the progress made during an odd indexed epoch. Recall that the only difference will be that the odd indexed epochs see the permutations in the order (σk1 , σ k 2 , . . . , σ k n) instead of (σ k n, σ k n−1, . . . , σ k 1 ). After doing the computation, we get the following equation:
  x2k0 = ( n∏ i=1 (I − Skn−i+1) ) x2k−10 + n∑ i=1 n−i∏ j=1 (I − Skn−j+1)  tki . Combining the results above, we can get the total progress made after the pair of epoch 2k − 1 and 2k:
  x2k+10 = ( n∏ i=1 (I − Ski ) ) x2k0 + n∑ i=1 n−i∏ j=1 (I − Skj )  tkn+1−i
  = ( n∏ i=1 (I − Ski ) )( n∏ i=1 (I − Skn−i+1) ) x2k−10 +  n∏ j=1 (I − Skj )  n∑ i=1 n−i∏ j=1 (I − Skn+1−j)  tki +
  n∑ i=1 n−i∏ j=1 (I − Skj )  tkn+1−i. (30) In the sum above, the first term will have an exponential decay, hence we need to control the next two terms. Similar to Theorem 4, we denote the sum of the terms as zk:
  zk :=  n∏ j=1 (I − Skj )  n∑ i=1 n−i∏ j=1 (I − Skn+1−j)  tki + n∑ i=1 n−i∏ j=1 (I − Skj )  tkn+1−i =
  n∑ i=1  n∏ j=1 (I − Skj ) n−i∏ j=1 (I − Skn+1−j)  tki + n∑ i=1 n−i∏ j=1 (I − Skj )  tkn+1−i. To see where the iterates end up after K epochs, we simply set 2k = K in Eq. 30 and then keep applying the equation recursively to preceding epochs. Then, we get
  xKn = x K+1 0 = ( n∏ i=1 (I − S K 2 i ) )( n∏ i=1 (I − S K 2 n−i+1) ) xK−10 + z K 2
  = (( n∏ i=1 (I − S K 2 i ) )( n∏ i=1 (I − S K 2 n−i+1) ))(( n∏ i=1 (I − S K 2 −1 i ) )( n∏ i=1 (I − S K 2 −1 n−i+1) )) xK−30
  + ( n∏ i=1 (I − S K 2 i ) )( n∏ i=1 (I − S K 2 n−i+1) ) z K 2 −1 + z K 2
  ...
  = K2 −1∏ k=0 ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )x10 +
  K 2 −1∑ k=0 ( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k.
  Taking squared norms and expectations on both sides, we get
  E[‖xKn ‖2] = E ∥∥∥∥∥∥ K2 −1∏
  k=0 ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )x10 +
  K 2 −1∑ k=0 ( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k ∥∥∥∥∥∥ 2 
  ≤ 2E  ∥∥∥∥∥∥ K2 −1∏
  k=0 ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )x10 ∥∥∥∥∥∥ 2 
  + 2E  ∥∥∥∥∥∥ K 2 −1∑ k=0 ( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k ∥∥∥∥∥∥ 2  , (31)
  where we used the fact that (a+ b)2 ≤ 2a2 + 2b2. Next, we expand the second term above to get
  E[‖xKn ‖2] ≤ 2E  ∥∥∥∥∥∥ K2 −1∏
  k=0 ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )x10 ∥∥∥∥∥∥ 2 
  + 2 K 2 −1∑ k=0 E ∥∥∥∥∥ ( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k ∥∥∥∥∥ 2  + 4 ∑
  0≤k′<k≤K2 −1
  E [〈( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k,
  k′−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′〉  . (32)
  We handle each of the three terms separately. Let Nk :=(∏n i=1(I − S K 2 −k i ) )(∏n i=1(I − S K 2 −k n−i+1) ) . The the first term can be written as:
  E  ∥∥∥∥∥∥ K2 −1∏
  k=0 ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )x10 ∥∥∥∥∥∥ 2  = E  ∥∥∥∥∥∥ K2 −1∏ k=0 Nk x10 ∥∥∥∥∥∥ 2 
  = E (x10)> K2 −1∏
  k=0
  Nk >K2 −1∏ k=0 Nk x10 
  = E (x10)> K2 −1∏
  k=1
  Nk >N>0 N0 K2 −1∏
  k=1
  Nk x10 
  = E (x10)> K2 −1∏
  k=1
  Nk > E [N>0 N0] K2 −1∏
  k=1
  Nk x10  ,
  where, in the last line, we used the fact that the permutations in every epoch are independent of the permutations in other epochs. Next, we have the following lemma that bounds the spectral norm of E [ N>k Nk ] for any k ∈ [0, K2 − 1]:
  Lemma 7. For any 0 ≤ α ≤ 316nL min{1, √
  n κ},∥∥E[N>k Nk]∥∥ ≤ 1− αnµ
  Note that K ≥ 55κ log(nK) max { 1, √
  n κ } =⇒ α ≤ 316nL min{1, n κ}, and thus this lemma.
  Using it, we get
  E  ∥∥∥∥∥∥ K2 −1∏
  k=0 ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )x10 ∥∥∥∥∥∥ 2 
  = E (x10)> K2 −1∏
  k=1
  Nk > E [N>0 N0] K2 −1∏
  k=1
  Nk x10 
  ≤ (1− nαµ)E (x10)> K2 −1∏
  k=1
  Nk >K2 −1∏ k=1 Nk x10 
  = (1− nαµ)E (x10)> K2 −1∏
  k=2
  Nk > E [N>1 N1] K2 −1∏
  k=2
  Nk x10 
  ...
  ≤ (1− nαµ)K/2‖x10‖2. (33)
  Next we look at the second term in Ineq. (32):
  E ∥∥∥∥∥ ( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k ∥∥∥∥∥ 2 
  ≤ E [( k−1∏ l=0 ( n∏ i=1 ‖I − S K 2 −l i ‖ 2 )( n∏ i=1 ‖I − S K 2 −l n−i+1‖ 2 )) ‖zK2 −k‖2 ]
  ≤ E [∥∥∥zK2 −k∥∥∥2] ,
  where in the last step we used that ‖I − Ski ‖ ≤ 1 for all i, k. To see why this is true, recall that Ski = αAσki . Further by Assumption 2, ‖Aσki ‖ ≤ L and hence as long as α ≤ 1/L, we have that ‖I − Ski ‖ ≤ 1.
  Next, note that for any k we can apply Lemma 4 on E[‖zK2 −k‖2]. Hence, we get the following bound on the second term:
  2 K 2 −1∑ k=0 E ∥∥∥∥∥ ( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k ∥∥∥∥∥ 2 
  ≤ 2K 2 (2n2α4L2(G∗)2 + 170n5α6L4(G∗)2 log n). (34)
  Finally, we focus on the third term in Ineq. (32). We have the following lemma that gives an upper bound for it: Lemma 8. Let α ≤ 12nL and n > 6. Then for 0 ≤ k ′ < k ≤ K2 − 1,
  E [〈( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k,k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′〉 
  ≤ 1000n2α4L2(G∗)2 + 2000n6α7L5(G∗)2 log n.
  Using Lemma 8, and inequalities (33) and (34) in Ineq. (32), we get E[‖xKn ‖2] ≤ 2(1− nαµ)K/2‖x10‖2 + 2n2Kα4L2(G∗)2 + 170n5Kα6L4(G∗)2 log n
  + 1000n2K2α4L2(G∗)2 + 2000n6K2α7L5(G∗)2 log n
  ≤ 2(1− nαµ)K/2‖x10‖2 + 1002n2Kα4L2(G∗)2 + 2170n6K2α7L4(G∗)2 log n. Substituting α = 10 lognKµnK gives us the desired result.
  F.1 PROOF OF LEMMA 7 Recall that we have defined Nk := (∏n i=1(I − S K 2 −k i ) )(∏n i=1(I − S K 2 −k n−i+1) ) . Hence,
  N>k = (( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) ))>
  = ( n∏ i=1 (I − S K 2 −k n−i+1) )>( n∏ i=1 (I − S K 2 −k i ) )>
  = ( n∏ i=1 (I − S K 2 −k i ) > )( n∏ i=1 (I − S K 2 −k n−i+1) > )
  = ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) ) = Nk,
  where we used the fact that S K 2 −k i are symmetric. Hence Nk is symmetric. Then,
  ‖E[N>k Nk]‖ = max x:‖x‖=1 E[x>N>k Nkx]
  = max x:‖x‖=1
  E[x>NkNkx].
  Next, note that ‖Nk‖ ≤ 1 as long as α ≤ 1nL . This combined with the fact that Nk is symmetric gives us that x>NkNkx ≤ x>Nkx. Hence, we get
  ‖E[N>k Nk]‖ ≤ max x:‖x‖=1 E[x>Nkx]
  = ‖E[Nk]‖
  = ∥∥∥∥∥E [( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k n−i+1) )]∥∥∥∥∥ = ∥∥∥∥∥∥E ( n∏
  i=1
  (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k i ) )>∥∥∥∥∥∥ . To complete the proof, we apply the following lemma from Ahn et al. (2020): Lemma 9 (Lemma 6, Ahn et al. (2020)). For any 0 ≤ α ≤ 316nL min{1,
  n κ} and k ∈ [K],∥∥∥∥∥∥E ( n∏ i=1 (I − S K 2 −k i ) )( n∏ i=1 (I − S K 2 −k i ) )>∥∥∥∥∥∥ ≤ 1− αnµ. Remark 2. To avoid confusion, we remark here that, the term ‘Sk’ in the paper Ahn et al. (2020) is the same as the term ‘ ∏n i=1(I − Ski )’ in our paper, and hence the original lemma statement in their paper looks different from what is written above.
  F.2 PROOF OF LEMMA 8
  We begin by decomposing the product into product of independent terms, similar to proof of Lemma 8 in Ahn et al. (2020). However, after that we diverge from their proof since we use FLIPFLOP specific analysis.
  E [〈( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k,k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′〉 
  = E (zK2 −k)>(k+1∏ l=k′ ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ))>  k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )>
  k′−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′  .
  Since k′ < k, we get that ( z K 2 −k )> , (∏k+1
  l=k′ (∏n i=1(I − S K 2 −l i ) )(∏n i=1(I − S K 2 −l n−i+1) ))> and k′∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )>k′−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ are independent. Hence, we can write the expectation as product of expectations:
  E [〈( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k,k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′〉 
  = E [( z K 2 −k )>] E (k+1∏ l=k′ ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ))> E
    k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′  .
  Applying the Cauchy-Schwarz inequality on the decomposition above, we get
  E [〈( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k,k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′〉 
  ≤ ∥∥∥E [zK2 −k]∥∥∥ ∥∥∥∥∥E [ k+1∏ l=k′ ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )]∥∥∥∥∥∥∥∥∥∥∥∥E   k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1)
  )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ ∥∥∥∥∥∥
  ≤ ∥∥∥E [zK2 −k]∥∥∥E[k+1∏
  l=k′ ( n∏ i=1 ‖I − S K 2 −l i ‖ )( n∏ i=1 ‖I − S K 2 −l n−i+1‖ )] ∥∥∥∥∥∥∥E   k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1)
  )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ ∥∥∥∥∥∥
  ≤ ∥∥∥E [zK2 −k]∥∥∥∥∥∥∥∥∥∥E   k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1)
  )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ ∥∥∥∥∥∥ ,
  where in the last step we used that ‖I − Ski ‖ ≤ 1 for all i, k. To see why this is true, recall that Ski = αAσki . Further by Assumption 2, ‖Aσki ‖ ≤ L and hence as long as α ≤ 1/L, we have ‖I − Ski ‖ ≤ 1. For the two terms in the product above, we have the following lemma:
  Lemma 10. If n > 6 and α ≤ 1nL , then ‖E[zK2 −k]‖ ≤ 28nα2LG∗ + 9α5L4n4G∗ √ 2n log n.
  Lemma 11. If n > 6 and α ≤ 12nL , then∥∥∥∥∥∥∥E   k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ ∥∥∥∥∥∥
  ≤ 32nα2LG∗ + 24α5L4n4G∗ √ 2n log n.
  Finally, using these lemma we get
  E [〈( k−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )) z K 2 −k,k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′〉 
  ≤ ( 28nα2LG∗ + 9α5L4n4G∗ √ 2n log n ) · ( 32nα2LG∗ + 24α5L4n4G∗ √ 2n log n )
  ≤ 896n2α4L2(G∗)2 + 960α7L5n5(G∗)2 √ 2n log n+ 432α10L8n9(G∗)2 log n
  ≤ 1000n2α4L2(G∗)2 + 2000α7L5n6(G∗)2 log n.
  F.3 PROOF OF LEMMA 10
  Since we are dealing with just a single epoch, we will skip the superscript. Using Lemma 5, we get
  ‖E[z]‖ = ∥∥∥∥∥∥E [ n∑ i=1 Siti ] + E  2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti ) + E
  n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )∥∥∥∥∥∥ ≤
  n∑ i=1 E [‖Si‖‖ti‖] + 2n−1∑ j=n+1 j−1∑ l=1 ∥∥∥∥∥E [( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti )]∥∥∥∥∥
  + n−1∑ j=1 j−1∑ l=1 ∥∥∥∥∥E [( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )∥∥∥∥∥ ] . (35)
  Define G∗ := maxi ‖bi‖. Then ‖Si‖‖ti‖ = ‖αAσki ‖‖αbσki ‖ ≤ α 2LG∗ and hence
  n∑ i=1 E [‖Si‖‖ti‖] ≤ nα2LG∗. (36)
  Next we bound the other two terms. Using Eq. (26), we get that for any l < j,
  E [( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )]
  = E [( I −
  l−1∑ p=1 Sp + l−1∑ p=1 p−1∑ q=1 ( q−1∏ r=1 (I − Sq) ) SrSp ) SlSj ( n−j∑ i=1 tn+1−i )]
  = n−j∑ i=1 E [SlSjtn+1−i]− E [ l−1∑ p=1 SpSlSj ( n−j∑ i=1 tn+1−i )]
  + E [ l−1∑ p=1 p−1∑ q=1 ( q−1∏ r=1 (I − Sq) ) SrSpSlSj ( n−j∑ i=1 tn+1−i )]
  = ∑ i>j,l E [SlSjti]− ∑ p<l,j<i E [SpSlSjti] + l−1∑ p=1 p−1∑ q=1 E [( q−1∏ r=1 (I − Sq) ) SrSpSlSj ( n−j∑ i=1 tn+1−i )] = ∑ i>j,l E [SlSjE[ti|Sl,Sj ]]− ∑ p<l,j<i E [SpSlSjE[ti|Sl,Sj ,Sp]]
  + l−1∑ p=1 p−1∑ q=1 E [( q−1∏ r=1 (I − Sq) ) SrSpSlSj ( n−j∑ i=1 tn+1−i )] .
  Since ∑n i=1 ti = 0, and we use uniform random permutations, E[ti|Sl,Sj ] = ∑ ti 6=tl, ti 6=tj ti n−2 = −tj−tl n−2 . Similarly, E[ti|Sl,Sj ,Sp] =
  −tj−tl−tp n−3 . Hence,∥∥∥∥∥E [( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )]∥∥∥∥∥ ≤ ∑ i>j,l ‖E [SlSjE[ti|Sl,Sj ]]‖+ ∑ p<l,j<i ‖E [SpSlSjE[ti|Sl,Sj ,Sp]]‖
  + l−1∑ p=1 p−1∑ q=1 E [( q−1∏ r=1 ‖I − Sq‖ ) ‖Sr‖‖Sp‖‖Sl‖‖Sj‖ ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ ]
  ≤ ∑ i>j,l E [ ‖Sl‖‖Sj‖ ‖tl + tj‖ n− 2 ] + ∑ p<l,j<i E [ ‖Sp‖‖Sl‖‖Sj‖ ‖tl + tj + tp‖ n− 3 ]
  + l−1∑ p=1 p−1∑ q=1 E [( q−1∏ r=1 ‖I − Sq‖ ) ‖Sr‖‖Sp‖‖Sl‖‖Sj‖ ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ ]
  ≤ 2n n− 2 α3L2G∗ + 3n2 2(n− 3) α4L3G∗ + α4L4 l−1∑ p=1 p−1∑ q=1 E [∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ ]
  (Since ‖Si‖ ≤ αL, ‖ti‖ ≤ αG∗) ≤ 4α3L2G∗ + 3nα4L3G∗ + α5L4n2G∗ √ 18n log n, (37)
  where we used Lemma 6 and the assumption that n > 6 in the last step.
  The third term in Ineq. (35) can be handled similarly. For any l, j:
  E [( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti )]
  = E [( I −
  l−1∑ p=1 S̃p + l−1∑ p=1 p−1∑ q=1 ( q−1∏ r=1 (I − S̃q) ) S̃rS̃p ) S̃lS̃j ( 2n−j∑ i=1 ti )]
  = 2n−j∑ i=1 E [ S̃lS̃jti ] − l−1∑ p=1 2n−j∑ i=1 E [ S̃pS̃lS̃jti ] + E [ l−1∑ p=1 p−1∑ q=1 ( q−1∏ r=1 (I − S̃q) ) S̃rS̃pS̃lS̃j ( 2n−j∑ i=1 ti )] .
  Now, it is easy to see that i 6= j and i 6= 2n − j + 1. Then, if i = l or i = 2n − l + 1 we use for that case the fact that ‖E[S̃lS̃jti]‖ ≤ α3L2G∗. For all other i, we can again use that E[ti|S̃l, S̃j ] = −tl−tjn−2 if l ≤ n or E[ti|S̃l, S̃j ] =
  −t2n−l+1−tj n−2 otherwise. Similarly, we can bound∥∥∥E [S̃pS̃lS̃jti]∥∥∥. Proceeding in a way similar to how we derived Ineq. (37), we get∥∥∥∥∥E [( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti
  )]∥∥∥∥∥ ≤
  2n−j∑ i=1 ∥∥∥E [S̃lS̃jti]∥∥∥+ l−1∑ p=1 2n−j∑ i=1 ∥∥∥E [S̃pS̃lS̃jti]∥∥∥ +
  l−1∑ p=1 p−1∑ q=1 ∥∥∥∥∥E [( q−1∏ r=1 (I − S̃q) ) S̃rS̃pS̃lS̃j ( 2n−j∑ i=1 ti )]∥∥∥∥∥ ≤ α3L2G∗ + 2n
  n− 2 α3L2G∗ + α4L3G∗ + nα4L3G∗ +
  3n2
  2(n− 3) α4L3G∗
  + l−1∑ p=1 p−1∑ q=1 ∥∥∥∥∥E [( q−1∏ r=1 (I − S̃q) ) S̃rS̃pS̃lS̃j ( 2n−j∑ i=1 ti )]∥∥∥∥∥ ≤ 5α3L2G∗ + 5nα4L3G∗
  + l−1∑ p=1 p−1∑ q=1 ∥∥∥∥∥E [( q−1∏ r=1 ‖I − S̃q‖ ) ‖S̃r‖‖S̃p‖‖S̃l‖‖S̃j‖ ( 2n−j∑ i=1 ti )]∥∥∥∥∥ ≤ 5α3L2G∗ + 5nα4L3G∗ + α5L4n2G∗ √ 18n log n. (38)
  Substituting Ineq. (36), (37) and (38) into (35), we get ‖E[z]‖ ≤ nα2LG∗ + 10n2α3L2G∗ + 10n3α4L3G∗ + 2α5L4n4G∗ √ 18n log n
  + 4n2α3L2G∗ + 3n3α4L3G∗ + α5L4n4G∗ √ 18n log n
  ≤ 28nα2LG∗ + 9α5L4n4G∗ √ 2n log n,
  where we used the assumption that α ≤ 1nL in the last step.
  F.4 PROOF OF LEMMA 11
  Define the matrix M as
  M := k′−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )>k′−1∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) .
  Since M is independent of (∏n
  i=1(I − S K 2 −k
  ′ i ) )(∏n i=1(I − S K 2 −k ′ n−i+1) ) and z K 2 −k ′ , using the
  tower rule, we get
  E   k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ 
  = E (( n∏ i=1 (I − S K 2 −k ′ i ) )( n∏ i=1 (I − S K 2 −k ′ n−i+1) ))> E[M ]z K 2 −k ′  . We will now drop the superscript K2 − k
  ′ for convenience. Hence, we need to control the following term:
  E (( n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) ))> E[M ]z  We define (S̃1, . . . , S̃2n) as (S̃1, . . . , S̃n) = (S1, . . . ,Sn) and (S̃n+1, . . . , S̃2n) = (Sn, . . . ,S1). Then, we use Eq. (24) to get((
  n∏ i=1 (I − Si) )( n∏ i=1 (I − Sn−i+1) )) = 2n∏ i=1 (I − S̃i)
  = I − 2n∑ j=1 S̃j + 2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j .
  Note that I − ∑2n j=1 S̃j is a constant matrix. Since S̃j = αAσj , we have that ‖S̃j‖ ≤ αL by
  Assumption 2. Hence, α ≤ 12nL then ‖I − ∑2n j=1 S̃j‖ ≤ 1. Further, α ≤ 1/L implies that ‖I − Si‖ ≤ 1, which implies ‖M‖ ≤ 1. Hence,∥∥∥∥∥∥E (( n∏
  i=1
  (I − Si) )( n∏ i=1 (I − Sn−i+1) ))> E[M ]z ∥∥∥∥∥∥ ≤
  ∥∥∥∥∥∥∥E  I − 2n∑ j=1 S̃j > E[M ]z  ∥∥∥∥∥∥∥
  + ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]z  ∥∥∥∥∥∥∥
  = ∥∥∥∥∥∥∥ I − 2n∑ j=1 S̃j > E[M ]E [z] ∥∥∥∥∥∥∥+ ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]z  ∥∥∥∥∥∥∥
  ≤ ‖E [z]‖+ ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]z  ∥∥∥∥∥∥∥ .
  We can apply Lemma 10 to bound ‖E [z]‖. So, we focus on the other term. Using Lemma 5,∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]z  ∥∥∥∥∥∥∥
  ≤ ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]( n∑ i=1 Siti ) ∥∥∥∥∥∥∥
  + ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]  2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti )  ∥∥∥∥∥∥∥
  + ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ] n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )  ∥∥∥∥∥∥∥
  = ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]( n∑ i=1 Siti )∥∥∥∥∥∥∥ +
  ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]  2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j ( 2n−j∑ i=1 ti )  ∥∥∥∥∥∥∥
  + ∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ] n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − Sp) ) SlSj ( n−j∑ i=1 tn+1−i )  ∥∥∥∥∥∥∥
  ≤ E  2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 ‖I − S̃p‖ ) ‖S̃l‖‖S̃j‖ E[‖M‖]∥∥∥∥∥ n∑ i=1 Siti ∥∥∥∥∥ + E
   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 ‖I − S̃p‖ ) ‖S̃l‖‖S̃j‖ E[‖M‖]  2n−1∑ j=n+1 j−1∑ l=1 ( l−1∏ p=1 ‖I − S̃p‖ ) ‖S̃l‖‖S̃j‖ ∥∥∥∥∥ 2n−j∑ i=1 ti ∥∥∥∥∥ 
  + E  2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 ‖I − S̃p‖ ) ‖S̃l‖‖S̃j‖ E[‖M‖] n−1∑ j=1 j−1∑ l=1 ( l−1∏ p=1 ‖I − Sp‖ ) ‖Sl‖‖Sj‖ ∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥  .
  Now, we use that ‖M‖ ≤ 1, ‖I − S̃i‖ ≤ 1, ‖S̃i‖ ≤ αL and ‖t̃i‖ ≤ αG∗:∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]z  ∥∥∥∥∥∥∥
  ≤ 4n3α4L3G∗ + 4n4α4L4E [∥∥∥∥∥ 2n−j∑ i=1 ti ∥∥∥∥∥ ] + n4α4L4E [∥∥∥∥∥ n−j∑ i=1 tn+1−i ∥∥∥∥∥ ] .
  Using Lemma 6, we get∥∥∥∥∥∥∥E   2n∑ j=1 j−1∑ l=1 ( l−1∏ p=1 (I − S̃p) ) S̃lS̃j > E[M ]z  ∥∥∥∥∥∥∥ ≤ 4n3α4L3G∗ + 15n4α5L4G∗ √ 2n log n.
  Putting everything together,∥∥∥∥∥∥∥E   k′∏ l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) )> k′−1∏
  l=0 ( n∏ i=1 (I − S K 2 −l i ) )( n∏ i=1 (I − S K 2 −l n−i+1) ) zK2 −k′ ∥∥∥∥∥∥
  ≤ (28nα2LG∗ + 9α5L4n4G∗ √ 2n log n) + (4n3α4L3G∗ + 15n4α5L4G∗ √ 2n log n)
  ≤ 32nα2LG∗ + 24α5L4n4G∗ √ 2n log n.
  G PROOF OF THEOREM 6
  Proof. We start off by defining the error term
  rk = ( n∑ i=1 ∇fi ( x2k−1i−1 ) − n∑ i=1 ∇fi ( x2k−10 )) + ( n∑ i=1 ∇fn−i+1 ( x2ki−1 ) − n∑ i=1 ∇fn−i+1 ( x2k−10 )) ,
  where k ∈ [K/2]. This captures the difference between true gradients that the algorithms observes, and the gradients that a full step of gradient descent would have seen.
  For two consecutive epochs of FLIPFLOP WITH INCREMENTAL GD, we have the following inequality
  ‖x2kn − x∗‖2 = ‖x2k−10 − x∗‖2 − 2α 〈 x2k−10 − x∗, n∑ i=1 ∇fi ( x2k−1i−1 ) + n∑ i=1 ∇fn−i+1 ( x2ki−1 )〉
  + α2 ∥∥∥∥∥ n∑ i=1 ∇fi ( x2k−1i−1 ) + n∑ i=1 ∇fn−i+1 ( x2ki−1 )∥∥∥∥∥ 2
  = ‖x2k−10 − x∗‖2 − 2α 〈 x2k−10 − x∗, 2n∇F (x 2k−1 0 ) 〉 − 2α 〈 x2k−10 − x∗, rk 〉 + α2
  ∥∥2n∇F (x2k−10 ) + rk∥∥2 ≤ ‖x2k−10 − x∗‖2 − 4nα [ Lµ
  L+ µ ‖x2k−10 − x∗‖2 +
  1
  L+ µ ∥∥∇F (x2k−10 )∥∥2] − 2α 〈 x2k−10 − x∗, rk 〉 + α2
  ∥∥2n∇F (x2k−10 ) + rk∥∥2 ≤ ‖x2k−10 − x∗‖2 − 4nα [ Lµ
  L+ µ ‖x2k−10 − x∗‖2 +
  1
  L+ µ ∥∥∇F (x2k−10 )∥∥2] − 2α 〈 x2k−10 − x∗, rk 〉 + 8α2n2
  ∥∥∇F (x2k−10 )∥∥2 + 2α2 ∥∥rk∥∥2 = ( 1− 4nα Lµ
  L+ µ
  ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 8α2n2 )∥∥∇F (x2k−10 )∥∥2 − 2α 〈 x2k−10 − x∗, rk 〉 + 2α2
  ∥∥rk∥∥2 , (39) where the first inequality is due to Theorem 2.1.11 in Nesterov (2004) and the second one is simply (a+ b)2 ≤ 2a2 + 2b2. What remains to be done is to bound the two terms with rk dependence. Firstly, we give a bound on the norm of rk:
  ‖rk‖ = ∥∥∥∥∥ ( n∑ i=1 ∇fi ( x2k−1i−1 ) − n∑ i=1 ∇fi ( x2k−10 )) + ( n∑ i=1 ∇fn−i+1 ( x2ki−1 ) − n∑ i=1 ∇fn−i+1 ( x2k−10 ))∥∥∥∥∥ ≤
  n∑ i=1 ∥∥∇fi (x2k−1i−1 )−∇fi (x2k−10 )∥∥+ n∑ i=1 ∥∥∇fn−i+1 (x2ki−1)−∇fn−i+1 (x2k−10 )∥∥ . Next, we will use the smoothness assumption and bounded gradients property (Lemma 1).
  ‖rk‖ ≤ L n∑ i=1 ∥∥x2k−1i−1 − x2k−10 ∥∥+ L n∑ i=1 ∥∥x2ki−1 − x2k−10 ∥∥ ≤ LGα
  n∑ i=1 i+ LGα n∑ i=1 (n+ i)
  = n(2n− 1)αGL.
  Hence, ∥∥rk∥∥2 ≤ 4n4α2G2L2. (40) For the rk term, we need a more careful bound. Since the Hessian is constant for quadratic functions, we use Hi to denote the Hessian matrix of function fi(·). We start off by using the definition of rk:
  rk = ( n∑ i=1 ∇fi ( x2k−1i−1 ) − n∑ i=1 ∇fi ( x2k−10 )) + ( n∑ i=1 ∇fn−i+1 ( x2ki−1 ) − n∑ i=1 ∇fn−i+1 ( x2k−10 ))
  = n∑ i=1 ( ∇fi ( x2k−1i−1 ) −∇fi ( x2k−10 )) + n∑ i=1 ( ∇fn−i+1 ( x2ki−1 ) −∇fn−i+1 ( x2k−10 )) =
  n∑ i=1 Hi ( x2k−1i−1 − x 2k−1 0 ) + n∑ i=1 Hn−i+1 ( x2ki−1 − x2k−10 ) ,
  where we used the fact that for a quadratic function f with Hessian H , we have that∇f(x)−∇f(y) = H(x − y). After that, we express x2k−1i−1 − x 2k−1 0 and x 2k i−1 − x 2k−1 0 as sum of gradient descent steps:
  rk = n∑ i=1 Hi i−1∑ j=1 −α∇fj(x2k−1j−1 ) + n∑ i=1 Hn−i+1  n∑ j=1 −α∇fj(x2k−1j−1 ) + i−1∑ j=1 −α∇fn−j+1(x2kj−1)  = −α
  n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−1j−1 ) − α n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2kj−1)  − α
  n∑ i=1 Hn−i+1  n∑ j=1 ∇fj(x2k−1j−1 )  = −α
  n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−10 ) − α n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2k0 )  − α
  n∑ i=1 Hn−i+1  n∑ j=1 ∇fj(x2k−10 )  − α
  n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−1j−1 )− fj(x 2k−1 0 )  − α
  n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2kj−1)−∇fn−j+1(x2k0 )  − α
  n∑ i=1 Hn−i+1  n∑ j=1 ∇fj(x2k−1j−1 )−∇fj(x 2k−1 0 )  = −2α
  n∑ i=1 Hi  n∑ j=1 ∇fj(x2k−10 ) + α n∑ i=1 Hi∇fi(x2k−10 )
  − α n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−1j−1 )− fj(x 2k−1 0 )  − α
  n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2kj−1)−∇fn−j+1(x2k0 ) 
  − α n∑ i=1 Hi  n∑ j=1 ∇fj(x2k−1j−1 )−∇fj(x 2k−1 0 )  Next, we use the fact that ∑n j=1∇fj(x) = n∇F (x). We will also again use the fact that for a quadratic function f with Hessian H , we have that∇f(x)−∇f(y) = H(x− y):
  rk = −2α n∑ i=1 Hi(n∇F (x2k−10 )) + α n∑ i=1 Hi(∇fi(x2k−10 )−∇fi(x∗)) + α n∑ i=1 Hi∇fi(x∗)
  − α n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−1j−1 )− fj(x 2k−1 0 )  − α
  n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2kj−1)−∇fn−j+1(x2k0 )  − α
  n∑ i=1 Hi  n∑ j=1 ∇fj(x2k−1j−1 )−∇fj(x 2k−1 0 )  = −2α
  ( n∑ i=1 Hi )2 (x2k−10 − x∗) + α n∑ i=1 H2i (x 2k−1 0 − x∗) + α n∑ i=1 Hi∇fi(x∗)
  − α n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−1j−1 )− fj(x 2k−1 0 )  − α
  n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2kj−1)−∇fn−j+1(x2k0 )  − α
  n∑ i=1 Hi  n∑ j=1 ∇fj(x2k−1j−1 )−∇fj(x 2k−1 0 )  = ak + bk,
  where the random variables ak, bk as
  ak := −2α ( n∑ i=1 Hi )2 (x2k−10 − x∗) + α n∑ i=1 H2i (x 2k−1 0 − x∗) + α n∑ i=1 Hi∇fi(x∗), and bk := −α n∑ i=1 Hi i−1∑ j=1 ∇fj(x2k−1j−1 )− fj(x 2k−1 0 )
   − α
  n∑ i=1 Hn−i+1 i−1∑ j=1 ∇fn−j+1(x2kj−1)−∇fn−j+1(x2k0 )  − α
  n∑ i=1 Hi  n∑ j=1 ∇fj(x2k−1j−1 )−∇fj(x 2k−1 0 )  . Again, using smoothness assumption and bounded gradients property (Lemma 1) we get,
  ‖bk‖ ≤ 3α2L2Gn3. (41) Next, we decompose the inner product of x2k−10 − x∗ and E [ rk ] in Eq. (39):
  −2α 〈 x2k−10 − x∗, rk 〉 = −2α 〈 x2k−10 − x∗,ak + bk 〉
  = −2α 〈 x2k−10 − x∗,ak 〉 − 2α 〈 x2k−10 − x∗, bk 〉 (42)
  For the first term in (42),
  −2α 〈 x2k−10 − x∗,ak 〉 = 4α2 〈 x2k−10 − x∗, ( n∑ i=1 Hi )2 (x2k−10 − x∗) 〉
  − 2α2 〈 x2k−10 − x∗, n∑ i=1 H2i (x 2k−1 0 − x∗) 〉
  − 2α2 〈 x2k−10 − x∗, n∑ i=1 Hi∇fi(x∗) 〉
  ≤ 4α2 〈 x2k−10 − x∗, ( n∑ i=1 Hi )2 (x2k−10 − x∗) 〉
  − 2α2 〈 x2k−10 − x∗, n∑ i=1 Hi∇fi(x∗) 〉
  = 4α2n2‖∇F (x2k−10 )‖2 − 2α2 〈 x2k−10 − x∗, n∑ i=1 Hi∇fi(x∗) 〉 ≤ 4α2n2‖∇F (x2k−10 )‖2 + 2α2nLG‖x 2k−1 0 − x∗‖. (43)
  For the second term in (42), we use Cauchy-Schwarz and Ineq. (41) −2α 〈 x2k−10 − x∗, bk 〉 ≤ 6α3L2Gn3‖x2k−10 − x∗‖. (44)
  Substituting (43) and (44) back to (42), we get −2α 〈 x2k−10 − x∗, rk 〉 ≤ 4α2n2‖∇F (x2k−10 )‖2 + 2α2nLG‖x 2k−1 0 − x∗‖
  + 6α3L2Gn3‖x2k−10 − x∗‖. (45)
  Substituting (40) (45) back to (39), we finally get a recursion bound for one epoch: ‖x2kn − x∗‖2 ≤ (
  1− 4nα Lµ L+ µ
  ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 8α2n2 )∥∥∇F (x2k−10 )∥∥2 − 2α 〈 x2k−10 − x∗, rk 〉 + 2α2
  ∥∥rk∥∥2 ≤ (
  1− 4nα Lµ L+ µ
  ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 8α2n2 )∥∥∇F (x2k−10 )∥∥2 4α2n2‖∇F (x2k−10 )‖2 + 2α2nLG‖x 2k−1 0 − x∗‖+ 6α3L2Gn3‖x 2k−1 0 − x∗‖
  + 8α4n4G2L2
  = ( 1− 4nα Lµ
  L+ µ
  ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 12α2n2 )∥∥∇F (x2k−10 )∥∥2 + 2α2nLG(1 + 3αLn2)‖x2k−10 − x∗‖+ 8α4n4G2L2.
  Next, we use the fact that 2ab ≤ λa2 + b2/λ (for any λ > 0) on the term 2α2nLG(1 + 3αLn2)‖x2k−10 − x∗‖ to get that
  2α2nLG(1 + 3αLn2)‖x2k−10 − x∗‖ ≤ ( α2nLG(1 + 3αLn2) )2 /(nαµ) + nαµ‖x2k−10 − x∗‖2
  = µ−1α3nL2G2(1 + 3αLn2)2 + nα‖x2k−10 − x∗‖2.
  Substituting this back we get, ‖x2kn − x∗‖2 ≤ (
  1− 4nα Lµ L+ µ
  ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 12α2n2
  )∥∥∇F (x2k−10 )∥∥2
  + µ−1α3nL2G2(1 + 3αLn2)2 + nαµ‖x2k−10 − x∗‖2 + 8α4n4G2L2 ≤ (1− 2nαµ+ nαµ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 12α2n2 )∥∥∇F (x2k−10 )∥∥2 + µ−1α3nL2G2(1 + 3αLn2)2 + 8α4n4G2L2 (Since µ ≤ L)
  = (1− nαµ) ‖x2k−10 − x∗‖2 − ( 4nα 1
  L+ µ − 12α2n2 )∥∥∇F (x2k−10 )∥∥2 + µ−1α3nL2G2(1 + 3αLn2)2 + 8α4n4G2L2
  Now, substituting the values of α and the bound on K, we get that 4nα 1L+µ − 12α 2n2 ≥ 0 and hence, ‖x2kn − x∗‖2 ≤ (1− nαµ) ‖x2k−10 − x∗‖2 + µ−1α3nL2G2(1 + 3αLn2)2 + 8α4n4G2L2 Now, iterating this for K/2 epoch pairs, we get
  ‖xKn − x∗‖2 ≤ (1− nαµ) K/2 ‖x10 − x∗‖2 +
  K
  2 µ−1α3nL2G2(1 + 3αLn2)2 + 4Kα4n4G2L2
  ≤ e−nαµK/2‖x10 − x∗‖2 + K
  2 µ−1α3nL2G2(1 + 3αLn2)2 + 4Kα4n4G2L2
  ≤ e−nαµK/2‖x10 − x∗‖2 +Kµ−1α3nL2G2 + 9Kµ−1α5n5L4G2 + 4Kα4n4G2L2 (Since (1 + a)2 ≤ 2 + 2a2)
  Substituting α = 6 lognKµnK gives us the desired result.
  H ADDITIONAL EXPERIMENTS
  Although our theoretical guarantees for FLIPFLOP only hold for quadratic objectives, we conjecture that FLIPFLOP might be able to improve the convergence performance for other classes of functions, whose Hessians are smooth near their minimizers. To see this, we also ran some experiments on 1-dimensional logistic regression. As we can see in Figure 3, the convergence rates are very similar to those on quadratic functions. The data was synthetically generated such that the objective function becomes strongly convex and well conditioned near the minimizer. Note that logistic loss is not strongly convex on linearly separable data. Therefore, to make the loss strongly convex, we ensured that the data was not linearly separable. Essentially, the dataset was the following: all the datapoints were z = ±1, and their labels were y = 1z>0 with probability 3/4 and y = 1z<0 with probability 1/4. Framing this as an optimization problem, we have
  min x F (x) := E [−y log(h(xz))− (1− y) log(1− h(xz))] ,
  where h(xz) = 1/(1 + e−xz). Note that x = − log 3 is the minimizer of this function, which is helpful because we can use it to compute the exact error. Similar to the experiment on quadratic functions, n was set to 800 and step size was set in the same regime as in Theorems 4, 5, and 6.