Tighter Convergence Bounds for Shuffled SGD via Primal-Dual Perspective

Authors: xufeng  cai
yin  cheuk
  lin
jelena  diakonikolas
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 18:50:32

Abstract: Stochastic gradient descent (SGD) is perhaps the most prevalent optimization method in modern machine learning. Contrary to the empirical practice of sampling from the datasets \emph{without replacement} and with (possible) reshuffling at each epoch, the theoretical counterpart of SGD usually relies on the assumption of \emph{sampling with replacement}. It is only very recently that SGD using sampling without replacement -- shuffled SGD -- has been analyzed with matching upper and lower bounds. However, we observe that those bounds are too pessimistic to explain often superior empirical performance of data permutations (sampling without replacement) over vanilla counterparts (sampling with replacement) on machine learning problems. Through fine-grained analysis in the lens of primal-dual cyclic coordinate methods and the introduction of novel smoothness parameters, we present several results for shuffled SGD on smooth and non-smooth convex losses, where our novel analysis framework provides tighter convergence bounds over all popular shuffling schemes (IG, SO, and RR). Notably, our new bounds predict faster convergence than existing bounds in the literature -- by up to a factor of $O(\sqrt{n})$, mirroring benefits from tighter convergence bounds using component smoothness parameters in randomized coordinate methods. Lastly, we numerically demonstrate on common machine learning datasets that our bounds are indeed much tighter, thus offering a bridge between theory and practice.

Main paper: 
 introduction: 
 Originally proposed in [38] , SGD has been broadly studied in the machine learning literature due to its effectiveness in large-scale settings, where full gradient computations are often computationally prohibitive. When applied to unconstrained finite-sum problems f i (x), (P) SGD performs the update x t = x t-1 -η∇f it (x t-1 ) for i t ∈ [n] ([n] := {1, . . . , n}), in each iteration t. Traditional theoretical analysis for SGD builds upon the assumption of sampling i t ∈ [n] with replacement according to a fixed distribution p = (p 1 , . . . , p n ) ⊤ over [n], which leads to Figure 1 : An illustration of the convergence behaviour of shuffled SGD for logistic regression problems on LIBSVM datasets luke, leu and a9a, where we use step sizes from existing bounds and our work. Due to randomness, we average over 20 runs for each plot and include a ribbon around each line to show its variance. However, as suggested by the concentration of L (see Section 4.1 and Appendix E), the variance across multiple runs is negligible, hence the ribbons are not observable.E it [∇f it (x t-1 )/(np it )] = ∇f (x t-1 ), and thus much of the (deterministic) gradient descent-style analysis can be transferred to this setting. By contrast, no such connection between the component and the full gradient can be established for shuffled SGD -which employs sampling without replacement -making its analysis much more challenging. As a result, despite its fundamental nature, there were no non-asymptotic convergence results for shuffled SGD until a very recent line of work [2, 12, 20, 21, 30, 31, 34, 35, 42] . All existing results consider general finite sum problems, with the same regularity condition constant (Lipschitz constant of f i or its gradient) assumed for all the component functions. As a result, the obtained convergence bounds are typically no better than for (full) gradient descent, and are only better than the bounds for SGD with replacement sampling if the algorithm is run for many full passes over the data [30, 34] .Furthermore, there is a large gap between the empirical performance of shuffled SGD and the predicted convergence rates from prior work [20, 30] . One cause for this discrepancy are overly pessimistic bounds on the step size in prior work, which are of order 1/(nL max ), where L max is the maximum smoothness constant over components f i in (P). In practice, the step sizes are tuned to achieve better convergence bounds than predicted by the current theory. We illustrate how restrictions on the step size affect convergence of shuffled SGD (with random permutations in each epoch) in Fig. 1 , where we plot the resulting optimality gap over full data passes when shuffled SGD is applied to logistic regression problems on standard datasets. To compare the effect of the step size η from prior work and our work, we choose take η = 1/( √ 2nL max ) based on [30] , and η = 1/(n L L) from our work, where L, L are our novel fine-grained, data-dependent smoothness parameters defined in Section 3 for smooth convex finite-sum problems with linear predictors. As can be observed from Fig. 1 , larger step sizes resulting from our theory lead to faster convergence of shuffled SGD and, as a result, our convergence bounds better predict the performance of shuffled SGD.Building on these insights, we introduce a fined-grained theoretical analysis to transparently show how the structure of the data and the possibly different Lipschitz constants of the component functions or their gradients affect the performance of shuffled SGD, thus providing a better explanation of the heuristic success of shuffled SGD in modern machine learning.SGD (with replacement) has been extensively studied in many settings (see e.g., [1, 9, 10, 38] for convex optimization). Compared to SGD, shuffled SGD usually exhibits faster convergence in practice [8, 37] , and is easier and more efficient to implement [5] . For each epoch k, shuffled SGDstyle algorithms perform incremental gradient updates based on the sample ordering (permutation of the data points) denoted by π (k) . There are three main choices of data permutations: (i) π (k) ≡ π for some fixed permutation of [n] for all epochs, where shuffled SGD reduces to the incremental gradient (IG) method; (ii) π (k) ≡ π where π is randomly chosen only once, at the beginning of the first epoch, referred to as the shuffle-once (SO) scheme; (iii) π (k) randomly generated at the beginning of each epoch, referred to as random reshuffling (RR).For general smooth convex settings, the convergence of shuffled SGD has been established only recently. For the number of epochs K sufficiently large, [31] proved a convergence rate Table 1 : Comparison of our results with state of the art, in terms of individual gradient oracle complexity required to output x out with E[f (x out ) -f (x * )] ≤ ϵ, where ϵ > 0 is the target error and x * is the optimal solution. Here, σ 2 * = 1 n n i=1 ∥∇f i (x * )∥ 2 2 , D = ∥x 0 -x * ∥ 2 , and generalized linear model refers to objectives of the form f (x) = 1 n n i=1 ℓ i (a ⊤ i x) as defined in Section 3. Parameters Lg , Lg are defined in Section 2 and satisfy Lg ≤ 1 n n i=1 L i and Lg ≤ L max . Parameters L, L, and Ḡ are defined in Section 3, and are discussed in the text of this section.PAPER COMPLEXITY ASSUMPTIONS STEP SIZE NGUYEN ET AL. [34] CHA ET AL. [12] (RR)O nLmaxD 2 ϵ + √ nLmaxσ * D 2 ϵ 3/2 fi: Lmax-SMOOTH, CONVEX O 1 nLmaxMISHCHENKO ET AL. [30] SHAMIR [42] (RR/SO)O B2 G 2 max ϵ 2 ) (K = 1, n = Ω(1/ϵ 2 ) ℓi: Gmax-LIPSCHITZ, CONVEX B-BOUNDED ITERATES, ∥ai∥ ≤ 1 GENERALIZED LINEAR MODEL O 1 √ n [Ours, Theorem 3] (RR/SO) O n ḠD 2 ϵ 2 ℓi: Gi-LIPSCHITZ, CONVEX GENERALIZED LINEAR MODEL O 1 n √ ḠK O(1/ √ nK)for RR, which leads to the complexity matching SGD. This result was later improved to O(1/(n 1/3 K 2/3 )) by [12, 30, 34] for K sufficiently large and with bounded variance assumed at the minimizer, while the same rate holds for SO [30] . These results were complemented by matching lower bounds in [12] , under sufficiently small step sizes as utilized in prior work. The results in [30, 34] require restricted O(1/(nL)) step sizes and reduce to O(1/K) for small K, acquiring the same iteration complexity as full-gradient methods. Unlike in strongly convex settings, we are not aware of any follow-up work with improvements under small K for smooth convex settings.The major difficulty in analyzing shuffled SGD comes from characterizing the difference between the intermediate iterate and the iterate after one full data pass, for which current analysis (see e.g., [30] in smooth convex settings) uses the global smoothness constant with a triangle inequality. Such a bound may be too pessimistic and fail capturing the nuances of intermediate progress of shuffled SGD, which leads to a small step size and large K restrictions. To provide a more fine-grained analysis that narrows the theory-practice gap for shuffled SGD, we notice that such a proof difficulty is reminiscent of the analysis of cyclic block coordinate methods relating the partial gradients to the full one. This natural connection was further emphasized in studies of cyclic methods with random permutations [24, 47] ; however, these results were limited to convex quadratics. More generally, it is possible to interpret shuffled SGD as a primal-dual method performing cyclic updates on the dual side (see (PD) in Section 2.1 and (PL-PD) in Section 3). We note here that prior work on dual coordinate methods [41] provided theoretical guarantees only for the algorithms that choose the dual coordinate to optimize uniformly at random, while the cyclic variant (related to shuffled SGD) had only been studied numerically up until this work. Further discussion of related work appears in Appendix A.In this work, we study the convergence rates of shuffled SGD in various settings through a unified primal-dual perspective, making intriguing connections to cyclic coordinate methods. This analysis framework is novel and allows us to leverage cyclic bias accumulation techniques on the dual side to obtain fine-grained convergence bounds. The obtained bounds mirror the improvements in randomized coordinate methods, which come from different coordinate smoothness parameters. While coordinate methods are no better than full-gradient methods in the worst case, on typical problem instances, they are much faster and the improvements come precisely from a more finegrained view of smoothness. We see a similar phenomenon in our analysis, which highlights the usefulness of the fine-grained smoothness characterizations introduced in our work.We provide improved bounds for all three popular data permutation strategies RR, SO and IG, in smooth convex settings. When the problem objective narrows to empirical risk minimization with linear predictors, we are able to exploit the data-dependent structure and uncouple the linear and nonlinear parts of the objective function, allowing us to provide tighter data-dependent bounds, up to a factor of O( √ n). Moreover, we show that our techniques extend to non-smooth convex settings, providing improved bounds over existing work.We summarize our results and compare them to the state of the art in Table 1 . As is standard, all complexity results in Table 1 are expressed in terms of individual (component) gradient evaluations. They represent the number of gradient evaluations required to construct a solution with (expected) optimality gap ϵ, given a target error ϵ > 0.Extensions to mini-batching and IG. When presenting our results for general finite-sum problems (in Section 2), we consider simple updates without mini-batching for ease of presentation and to avoid introducing excessive notation. However, we emphasize that all our results can be extended to shuffled SGD with mini-batching. Our results are also the first to provide convergence bounds that demonstrate benefits of mini-batching in shuffled SGD. For completeness and generality, the proofs in the appendix are carried out for mini-batch settings with arbitrary batch sizes b ∈ {1, . . . , n}. Thus, all the results stated in Section 2 can be recovered by setting b = 1. Moreover, our framework can provide similar fine-grained convergence bounds for IG. However, as IG is not as commonly used in practice compared to RR and SO and due to space constraints, we only present our results for RR and SO in the main body and include the results for IG in the appendix.We consider a real d-dimensional Euclidean space (R d , ∥ • ∥) where d is finite and ∥ • ∥ is the ℓ 2 -norm. For a vector x, we let x j denote its j-th coordinate. For any positive integer m, we use [m] to denote the set {1, 2, . . . , m}. Given a matrix A, ∥A∥ := sup x∈R d ,∥x∥≤1 ∥Ax∥ denotes its operator norm. For a positive definite matrix Λ, ∥ • ∥ Λ denotes the Mahalanobis norm, ∥x∥ Λ := ⟨Λx, x⟩. We use I to denote the identity matrix, and diag(v) to denote the diagonal matrix with vector v on the main diagonal. For any j ∈ [n], we define I j↑ as the matrix obtained from the identity matrix I by setting the first j diagonal elements to zero, and let I j be the matrix with only the j-th diagonal element nonzero and equal to 1. To handle the cases with random data permutations, we use the following definitions corresponding to the data permutation π = {π 1 , π 2 , . . . , π n } of [n]: A π := a π1 , a π2 , . . . , a πn ⊤ permuting the rows based on π given a matrix A = [a 1 , a 2 , . . . , a n ] ⊤ , and v π := v π1 , v π2 , . . . , v πn ⊤ permuting the coordinates/subvectors based on π given a vectorv = (v 1 , v 2 , . . . , v n ) ⊤ . Algorithm 1 Shuffled SGD (Primal-Dual View, General Convex Smooth) 1: Input: Initial point x 0 ∈ R d , step size {η k } > 0, number of epochs K > 0 2: for k = 1 to K do 3:Generate some permutation π k of [n] (either deterministic or random) 4:x k-1,1 = x k-1for i = 1 to n in the ordering of π k do 6:y i k = arg max y i ∈R d y i , x k-1,i -f * i (y i ) 7: x k-1,i+1 = arg min x∈R d y i k , x + 1 2η k ∥x -x k-1,i ∥ 2 = x k-1,i -η k ∇f (x k-1,i ) 8:end for 9:x k = x k-1,n+1 , y k = y 1 k , y 2 k , . . . , y n k ⊤ 10: end for 11: Return: xK = K k=1 η k x k / K k=1 η k 2.1 Primal-dual view of shuffled SGD Problem (P) can be reformulated into a primal-dual form using the standard Fenchel conjugacy argument (see, e.g., [13, 14] ),EQUATIONwhere we slightly abuse the notation to denote y = (y 1 , . . . ,y n ) ⊤ ∈ R nd and f * i is the convex conjugate of f i defined by f * i (y) = sup x∈R d ⟨y, x⟩ -f i (x) . We let y x = (y 1 x , . . . , y n x ) ⊤ ∈ R nd be the conjugate pair of x ∈ R d , i.e., y i x = arg max y∈R d {⟨y, x⟩-f * i (y)}, and we denote y * = y x * . Given a primal-dual pair (x, y), the primal-dual gap of (PD) is defined by Gap(x, y) = max (u,v) {L(x, v) -L(u, y)}. In particular, we consider the pair (x, y * ) for x ∈ R d , and bound Gap v (x, y * ) := L(x, v) -L(x * , y * ) for an arbitrary but fixed v. To finally obtain the function value gap f (x) -f (x * ) for (P), we only need to choose v = arg max w L(x, w) = y x .Using this primal-dual formulation and standard convex conjugacy arguments, we can equivalently write the standard shuffled SGD algorithm in a primal-dual form as summarized in Algorithm 1.Improved bounds with new smoothness constants. To prove a convergence bound for shuffled SGD in this general setting, we first construct an upper estimate of Gap v (x k , y * ) for some fixed v to be set later, as summarized in the following lemma. Lemma 1. Under Assumption 1, for any k ∈ [K], the iterates {y i k } n i=1 and {x k-1,i } n+1 i=1 generated by Algorithm 1 satisfyEQUATIONwhereE k := η k Gap v (x k , y * ) + 1 2n ∥x * -x k ∥ 2 2 -1 2n ∥x * -x k-1 ∥ 2 2 , v k = v π (k), and y * ,k = y * ,π (k) are the (block-wise) permuted vectors based on the permutation π k at the k-th epoch.We note that the first termT 1 := η k n n i=1 y i k , x k -x k-1,i+1from Lemma 1 can be aggregated into the terms capturing the primal progress within one epoch and cancelled by the last term in Eq. ( 1). The precise bound on T 1 and its proof are provided in Lemma 10 in Appendix C.1. The second termT 2 := η k n n i=1 v i k -y i k , x k -x k-1,irequires us to relate the intermediate iterate x k-1,i to the iterate x k after one full data pass, which corresponds to a partial sum of the component gradients, each at different iterates {x k-1,j } n j=i . In contrast to prior analyses (e.g., Mishchenko et al. [30] ) using the global smoothness and triangle inequality to bound this partial sum, we provide a tighter bound on T 2 that tracks the progress of the cyclic update on the dual side, in the aggregate.To simplify the notation in the following lemmas and to clearly compare our results, we introduce the following novel definitions of smoothness constants for shuffled SGD:EQUATIONwhere I (di) = di j=d(i-E[T 2 ] ≤ E η 3 k n Lg π (k) Lg π (k) ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k + η 3 k (n + 1) Lg 6 σ 2 * ,whereT 2 := η k n n i=1 v i k -y i k , x k -x k-1,i , v k = v π (k) and y * ,k = y * ,π (k) .With Lemmas 1 and 2 in tow, we are ready to present the main result of this section.Theorem 1. Under Assumptions 1 and 2, if η k ≤ 1 n 2 Lg π (k) Lg π (k) and H K = K k=1 η k , the output xK of Algorithm 1 with uniformly random (RR/SO) shuffling satisfies E[H K (f ( xK ) -f (x * ))] ≤ 1 2n ∥x 0 -x * ∥ 2 2 + K k=1 η 3 k (n + 1) Lg 6 σ 2 * .As a consequence, for any ϵ > 0, there exists a choice of a constant step size η k = η for whichE[f ( xK ) -f (x * )] ≤ ϵ after O n √ Lg Lg ∥x0-x * ∥ 2 2 ϵ + √ n Lg σ * ∥x0-x * ∥ 2 2 ϵ 3/2individual gradient queries. 
limitations: 
 Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: Yes, see detailed discussion in the introdcution. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
tighter bounds for convex finite-sum problems with linear predictors: 
 To study the effect of the structure of the data on the convergence of shuffled SGD, we sharpen the focus from general finite-sum problems to convex finite-sum with linear predictors:EQUATIONwherea i ∈ R d (i ∈ [n]) are data vectors and ℓ i : R → R are convex and either smooth or Lipschitz nonsmooth functions associated with the linear predictors ⟨a i , x⟩ for i ∈ [n]. In addition to their explicit dependence on the data, it is worth noting that problems of the form (PL) cover most of the standard convex ERM problems where shuffled SGD is commonly applied, such as support vector machines, least absolute deviation, least squares, and logistic regression.Problem (PL) admits an explicit primal-dual formulation using the standard Fenchel conjugacy argument (see, e.g., [13, 14] ),min x∈R d max y∈R n L(x, y) := 1 n ⟨Ax, y⟩ - 1 n n i=1 ℓ * i (y i ) = 1 n n i=1 a ⊤ i xy i -ℓ * i (y i ) , (PL-PD)where A = [a 1 , a 2 , . . . , a n ] ⊤ ∈ R n×d is the data matrix and ℓ * i : R → R is the convex conjugate of ℓ i . This observation allows us to again interpret without-replacement SGD updates as cyclic coordinate updates on the dual side. Note that due to the objective structure in (PL), the primal-dual formulation (PL-PD) can decouple the linear (a ⊤ i x) and the non-linear (ℓ i ) parts within individual loss functions f i . We redefine the conjugate pair of x ∈ R d to be y x = (y 1x , . . . , y n x ) ⊤ ∈ R n , with y i x = arg max y i ∈R {y i a ⊤ i x -ℓ * i (y i )}. In this section, we consider shuffled SGD with mini-batch estimators of size b and assume without loss of generality that n = bm for some positive integer m. The detailed primal-dual view of shuffled SGD adapted to (PL-PD) and mini-batch estimators is provided in Alg. 2 in Appendix C.Throughout this subsection, we make the following (standard) assumptions, corresponding to Assumptions 1 and 2 from Section 2. Assumption 3. Each ℓ i is convex andL i -smooth (i ∈ [n]), i.e., |ℓ ′ i (x) -ℓ ′ i (y)| ≤ L i |x -y| for any x, y ∈ R. There exists a minimizer x * ∈ arg min x∈R d f (x).We remark that Assumption 3 implies that both f and each component functionf i (x) = ℓ i (a ⊤ i x) are L max -smooth, where L max = max i∈[n] L i ∥a i ∥ 22 . Assumption 3 also implies that each convex conjugate ℓ * i is 1 Li -strongly convex [3] . In the following, we let Λ = diag(L 1 , L 2 , . . . , L n ), andΛ π = diag L π 1 , L π 2 , . . . , L π n , given a permutation π of [n].We further assume bounded variance at x * , same as prior work [30, 34, 45, 46] .Assumption 4. σ 2 * := 1 n n i=1 ∥∇f i (x * )∥ 2 = 1 n n i=1 (ℓ ′ i (a ⊤ i x * )) 2 ∥a i ∥ 2 2 is bounded.Improved bounds with new smoothness constants. Our convergence bounds depend on the smoothness parameters defined in Eq. ( 3) below. We provide a detailed discussion on how these parameters relate to traditional smoothness parameters both in the worst case and on typical datasets, in Section 4.1, with additional numerical results provided in Appendix E.EQUATIONwhere I (j) := bj i=b(j-1)+1 I i . In comparison to the smoothness constants defined in Eq. ( 2) for general finite-sum problems, we note that the constants in Eq. (3) applying to generalized linear models are tighter and more informative estimates, as the data matrix A and the smoothness constants from the nonlinear part Λ are separated in Eq. (3). Thus, the constants Lπ and Lπ directly depend on the data matrix, which explicitly demonstrates how the structure of the data affects the convergence of shuffled SGD. The following theorem states the convergence of Algorithm 2 with these new refined smoothness constants, while its proof is provided in Appendix C. E[H K (f ( xK ) -f (x * ))] ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .As a result, given ϵ > 0, there exists a constant step sizeη k = η such that E[f ( xK ) -f (x * )] ≤ ϵ after O n √ L L∥x0-x * ∥ 2 2 ϵ + (n-b)(n+b) n(n-1) √ n Lσ * ∥x0-x * ∥ 2 2 ϵ 3/2individual gradient queries.A few remarks are in order here. When b = n, we recover the standard guarantee of gradient descent, which serves as a sanity check as in this case the algorithm reduces to standard gradient descent.When ϵ = Ω((n-b)(n+b)σ 2 * n 2 (n-1) L ), the resulting complexity is O n √ L L∥x0-x * ∥ 2 2 ϵ. Observe that this case can happen when either ϵ is large (compared to, say, 1/n) or when σ * is small (it is, in fact, possible for σ * to be zero, which happens, for example, when the data rows are linearly independent). Unlike in bounds from previous work, we observe from our bounds the benefit of using shuffled SGD compared to full gradient descent, where the difference is by a factor that can be as large as √ n, as we have discussed in the introduction (see also Section 4) . When ϵ = O((n-b)(n+b)σ 2 * n 2 (n-1)L ), the second term in our complexity bound dominates. In this case, when b = 1, we recover the state of the art results from [12, 30, 34] , while for b > 1 our bound provides the Ω n(n-1)(n-b)(n+b) • LL -factor improvement, providing insights into benefits from the mini-batching strategy commonly used in practice.In non-smooth settings, we make the following standard assumption. Assumption 5. Each ℓ i is convex andG i -Lipschitz (i ∈ [n]), i.e., |ℓ i (x) -ℓ i (y)| ≤ G i |x -y| for any x, y ∈ R; thus |g i (x)| ≤ G i where g i (x) ∈ ∂ℓ i (x). There exists a minimizer x * ∈ arg min x∈R d f (x). If Assumption 5 holds, each ℓ i (a ⊤ i x) is also G max -Lipschitz with respect to x, where G max = max i∈[n] G i ∥a i ∥ 2 .To state our results, we defineEQUATIONπn , given a data permutation π of [n]. We now extend our analysis of Algorithm 1 to convex nonsmooth Lipschitz settings, where the conjugate functions ℓ * i (y i ) are only convex. Proceeding as in Lemma 1, we obtain a bound on the primal-dual gap similar to (1), but lose two retraction terms induced by smoothness. Instead of cancelling the corresponding error terms like in the smooth case, we rely on the boundedness of the subgradients to bound these terms under a sufficiently small step size, which is common in nonsmooth Lipschitz settings. Similar to Section 2, we introduce the following quantities to obtain a tighter guarantee with respect to the data matrix and Lipschitz constants:= 1 mn Γ 1/2 π m j=1 I b(j-1)↑ A π A ⊤ π I b(j-1)↑ Γ 1/2 π 2 , Gπ := 1 b Γ 1/2 π m j=1 I (j) A π A ⊤ π I (j) Γ 1/2 π 2 .We discuss the improvements in convergence from Ĝπ and Gπ in Section 4, while the convergence of Algorithm 2 is described in Theorem 3, with its proof deferred to Appendix D. E[H K (f ( xK ) -f (x * ))] ≤ 1 2n ∥x 0 -x * ∥ 2 2 + K k=1 2η 2 k n Ḡ,As a result, for any ϵ > 0, there exists a step sizeη k = η such that E[f ( xK ) -f (x * )] ≤ ϵ after O n Ḡ∥x0-x * ∥ 2 2 ϵ 2individual gradient queries.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
discussion of our new smoothness constants and numerical results: 
 To succinctly explain where our improvements come from, we now consider (PL) where ℓ i is 1-smooth and b = 1, ignoring the gains from the mini-batch estimators (for large K) and our softer guarantee that handles individual smoothness constants. For this specific case, L = L max = max 1≤i≤n ∥a i ∥ 2 , and thus our results for the smooth case and the RR and SO variants match state of the art in the second term, which dominates when there are many (K = Ω(L 2 max D 2 n σ 2 *)) epochs. When there areK = O( L 2 max D 2 n σ 2 *) epochs in the SO and RR variants or for all regimes of K in the IG variant, the difference between our and state of the art bounds comes from the constant L that replaces L max , and our improvement is by a factor L max / L. Note that O( nLmax ϵ ) from prior bounds, which is the dominating term in the small K regime, is even worse than the complexity of full gradient descent, as the full gradient Lipschitz constant of f in this case is1 n ∥AA ⊤ ∥ 2 ≤ L max .Given a worst-case permutation π, and denoting by A π the data matrix A with its rows permuted according to π, our constant L can be bounded above by L max using the following sequence of inequalities:EQUATIONwhere (i) holds by the triangle inequality, (ii) holds because the operator norm of the matrix I (j-1)↑ A π A ⊤ π I (j-1)↑ (equal to the operator norm of the bottom right (n -j + 1)× (n -j + 1) submatrix of A π A ⊤ π ) is always at most ∥A π A ⊤ π ∥ = ∥AA ⊤ ∥,for any permutation π, and (iii) holds by bounding above the operator norm of a symmetric matrix by its trace. Hence L is never larger than L max , but can generally be much smaller, due to the sequence of inequality relaxations in (4) . While each of these inequalities can be loose, we emphasize that (iii) is almost always loose, by a factor that can be as large as n.As a specific example where L is smaller than L max by a factor of n, consider the example of Gaussian data, where we draw n i.i.d. standard Gaussian vectors from N (0, I d ) and take d = n. By standard concentration results, with high probability, all columns/rows of A π in this case are near-orthogonal (see, e.g., [7, Chapter] ) and ∥a i ∥ 2 2 ≈ d = n for all i. As a result, the operator norm to trace inequality (iii) is loose by a factor d = n, with high probability. Note that in this example all individual smoothness parameters of components f i are essentially the same (w.h.p.) and equal ∥a i ∥ 2 2 , thus the improvement of our bound on the smoothness parameter does not come from averaging but from the structure of the data. This observation is important for contrasting the results from Section 2 and Section 3. In particular, focusing solely on the finite sum structure and ignoring the structure of the data matrix would provide no improvements in the resulting convergence bounds.As further evidence, we empirically evaluate L max / L on 15 large-scale machine learning datasets and demonstrate that on those datasets L max / L is of the order n α , for α ∈ [0.15, 0.96] (see Sec. 4.1 for more details), providing strong evidence of a tighter guarantee as a function of n.For the nonsmooth settings, by a similar sequence of inequalities, we can show that Ḡ ≤ G 2 max , which can be loose by a factor 1/n due to the operator norm to trace inequality. Thus, our bound is never worse than what would be obtained from the full subgradient method, but can match the bound of standard SGD, or even improve 1 upon it for at least some data matrices A.In this section, we provide empirical evidence to support our claim about usefulness of the new convergence bounds obtained in our work. In particular, we conduct numerical evaluations to compare L to the classical smoothness constant L on synthetic datasets and on popular machine learning benchmark datasets.For a more streamlined comparison and to focus on the dependence on the data matrix, we assume that the loss functions ℓ i all have the same smoothness constant, which leads toL max / L = (max 1≤i≤n {∥a i ∥ 2 })/ 1 n 2 ∥ n j=1 I (j-1)↑ A π A ⊤ π I (j-1)↑ ∥ 2 .Since the scale of the smoothness constant of the loss functions is irrelevant for the ratio L max / L in this case, for simplicity, we take it to equal one. Note that assuming different smoothness constants over component loss functions would only make our bound better compared to related work (see Eq. ( 3) and the discussion following it).We also compare L and L max on a number of benchmarking datasets from LIBSVM [15] , MNIST [17] , CIFAR10 [22] , and Broad Bioimage Benchmark Collection [28] . For each dataset, we generate a uniformly random permutation π for the data matrix A and compute Lπ . We repeat this procedure 1000 times for all datasets and display the average L max / Lπ in Table 2 , except for e2006train, CIFAR10, MNIST, and BBBC005 where we do 20 repetitions due to limitations of Finally, as a justification for using the empirical mean of Lπ over random permutations π in the results displayed in Table 2 , we observe in our evaluations that the values of L max / Lπ are fairly concentrated around their empirical mean values. Histogram plots showing the empirical distributions of L max / Lπ for each of the datasets are provided in Appendix E.We conclude with a few additional remarks. Our results indicate that the structure of the data is important for predicting behavior of popular machine learning methods such as variants of shuffled SGD considered in our work, and thus should be incorporated in their study: as demonstrated in the Gaussian data example, considering simple finite sum structure and ignoring the dependence on the data can lead to overly pessimistic bounds. Thus it would be interesting to provide a further theoretical study of shuffled SGD that incorporates distributional assumptions for the data. Additionally, as mentioned in the previous paragraph, we empirically observed that permutation-dependent parameter Lπ concentrates around its mean for permutations generated uniformly at random. Thus, it would be interesting to consider whether our theoretical results can be strengthened to depend on the mean value of Lπ (as opposed to maximum). We leave such considerations for future work.Outline. The supplementary material of the paper is organized as follows:• Section A provides a brief survey on shuffled SGD and its related work.• Section B presents the proofs related to the smooth convex setting from Section 2, where we only assume each component function f i to be convex and L i -smooth. • Section C presents the proofs related to the smooth convex setting with linear predictors from Section 3. • Section D presents the proofs related to the non-smooth convex setting with linear predictors from Section 3.• Section E presents the full details of the computational experiments performed in the paper.In this section, we continue with the discussion on the background of shuffled SGD from Section 1. We would like to briefly recall that shuffled SGD usually performs better in practice when compared to SGD, and is also easier and more efficient to implement. However, in terms of the theoretical analysis, sampling without replacement introduces the sampling bias at each iteration, making it difficult to approximate shuffled SGD by full gradient descent. Using empirical observations, shuffled SGD was conjectured to converge much faster than SGD with replacement, based on the noncommutative arithmetic-geometric mean inequality conjecture [36] , which was later proved to be false [16, 23] . As a consequence, whether or not shuffled SGD can be faster than SGD at least in some regimes remained open [8] until a breakthrough result in [20] , where it was shown that for the class of smooth strongly convex optimization problems, the convergence of the RR variant of shuffled SGD is essentially of the order-(1/K 2 ) for K full passes of the data (also called epochs), which is faster than order-(1/nK) convergence of SGD for sufficiently large K. This bound for the smooth strongly convex case was later improved under various regimes and additional assumptions [2, 21, 30, 31, 34, 42] , while the tightest of those bounds were matched by lower bounds in [12, 35, 39, 51] .Since our results are for the general (non-strongly) convex regimes, in this section we focus on the results that apply to those (convex, smooth or nonsmooth Lipschitz) regimes. For convex nonsmooth Lipschitz problems, we are only aware of the results in [42] . These results are only useful when the number of data passes K is small and the number of component functions n is large, as they contain an irreducible order-1 √ n error, and are not directly comparable to our results. For the IG variant of SGD without replacement (deterministic order), asymptotic convergence was established in [6, 29] , with further convergence results for both smooth and nonsmooth settings provided in [18, 25, 30, 32, 34, 50] . As IG does not benefit from randomization, it is known to have a worse convergence bound than RR under the Lipschitz Hessian assumption [18, 21] , which was also shown in more general settings [30] .In this paper, we viewed shuffled SGD as a primal-dual method where the updates are performed on the dual side in a cyclic manner, thus we can leverage techniques from general cyclic methods. However, in contrast to randomized methods (corresponding to standard SGD), cyclic methods are usually more challenging to analyze [33] , basic variants exhibit much worse worst-case complexity than even full gradient methods [4, 19, 26, 26, 40, 44, 48, 49] , with more refined results being established only recently [11, 27, 43] . While the inspiration for our work came from these recent results [11, 27, 43] , they are completely technically disjoint. First, all these results rely on nonstandard block Lipschitz assumptions, which are not present in our work. Second, all of them leverage proximal gradient-style cyclic updates to carry out the analysis, which is inapplicable in our case for the cyclic updates on the dual side, as otherwise the method would not correspond to (shuffled) SGD. Finally, [27, 43] utilize extrapolation steps, which would break the connection to shuffled SGD in our setting, while [11] relies on a gradient descent-type descent lemma, which is impossible to establish in our setting.In this section, we consider the general finite-sum setting where we assume that each component function f i is convex and smooth, and derive the refined analysis under this setting. Here we focus on the smooth convex problems as prior work did [30, 34] , since smoothness is essential to showing the advantage of shuffled SGD [31] over SGD, otherwise the rate of SGD is optimal. In particular, we study the general smooth convex finite-sum problem (P)min x∈R d f (x) := 1 n n i=1 f i (x) , (P)where each f i is convex and smooth. (P) is equivalent tomin x∈R d max y∈R nd L(x, y) := 1 n n i=1 y i , x -f * i (y i ) = 1 n ⟨Ex, y⟩ - 1 n n i=1 f * i (y i ) , (PD)where we slightly abuse the notation in this section and use y i ∈ R d to be the i-th d elements of the vector y such that y = (y 1 , . . . ,y n ) ⊤ ∈ R nd , E = [I d , . . . , I d n ] ⊤ ∈ R nd×d is the vertical concatenation of n identity matrices I d ∈ R d×d and f * i is the convex conjugate of f i defined by f i (x) = sup y i ∈R d y i , x -f * i (y i ).In the following, we consider the mini-batch estimator of batch size b, and let y (i) ∈ R bd denote the vector comprised of the i th bd elements of y. For simplicity and without loss of generality, we assume that n = bm for some positive integer m, so that y = (y (1) , . . . , y (m) ) ⊤ . Note that if choosing b = 1, our setting is the same as the ones in [30, 34] . Then we have the primal-dual view of shuffled SGD scheme for general smooth convex minimization as in Alg. 1, whereEQUATIONn } of [n] at the k-th epoch, we use the same notation of v k = (v π (k) 1 , . . . , v π (k) n ) ⊤ ∈ R nd , y * ,k = (y π (k) 1 * , . . . , y π (k) n * ) ⊤ ∈ R nd as in previous sections except now each v π (k) i , y π (k) i *are d-dimensional subvectors. Further, we denote the permuted smoothness constant matrices byΛ k = diag(L π (k) 1 , . . . , L π (k) 1 d , . . . , L π (k) n , . . . , L π (k) n d ) ∈R nd×nd , and we use I for I nd ∈ R nd×nd throughout this section.New smoothness constants and comparisons. We first recall the new smoothness constants for any permutation π of [n], defined in Eq. ( 2):Lg π := 1 mn Λ 1/2 π m i=1 I bd(i-1)↑ EE ⊤ I bd(i-1)↑ Λ 1/2 π 2 , Lg = max π Lg π , Lg π := 1 b Λ 1/2 π m i=1 I (di) EE ⊤ I (di) Λ 1/2 π 2 , Lg = max π Lg π ,whereI (di) = bdi j=bd(i-1)+1 I j .To compare Lg π and L := max i∈[n] L i , we make use of the Kronecker product with notation ⊗ defined byA ⊗ B =    A 11 B • • • A 1n B . . . . . . . . . A m1 B • • • A nn B   for two matrices A ∈ R m×n and B ∈ R p×q . The following lemma states a useful fact for the Kronecker product.Lemma 3. For square matrices A and B of sizes p and q and with eigenvalues λ i (i ∈ [p]) and µ j (j ∈ [q]) respectively, the eigenvalues ofA ⊗ B are λ i µ j for i ∈ [p], j ∈ [q].We now use the following chain of inequalities to compare Lg π and L for any permutation π of [n]:Lg π = 1 mn Λ 1/2 π m i=1 I bd(i-1)↑ EE ⊤ I bd(i-1)↑ Λ 1/2 π 2 ≤ 1 n Λ 1/2 EE ⊤ Λ 1/2 2 = 1 n (l π l ⊤ π ) ⊗ I d 2 (i) = 1 n n i=1 L i ≤ L,where we define l π = L π1 , L π2 , . . . , L πn ⊤ . For (i), we use Lemma 3 and notice that the eigenvalues of I d all equal 1, while the largest eigenvalue ofl π l ⊤ π = ∥l∥ 2 2 = n i=1 L i , so the operator norm of (l k l ⊤ k ) ⊗ I d is n i=1 L i . To compare Lgπ and L, we notice thatΛ 1/2 π m i=1 I (di) EE ⊤ I (di) Λ 1/2 π = m i=1 I (di) Λ 1/2 π EE ⊤ Λ 1/2 π I (di)is a block diagonal matrix whose operator norm is the maximum of the operator norms over its diagonal block submatrices, so we haveLg π = 1 b max i∈[m] I (di) Λ 1/2 π EE ⊤ Λ 1/2 π I (di) = 1 b max i∈[m] I (di) (l π l ⊤ π ) ⊗ I d I (di) (i) = max i∈[m] 1 b b j=1 L π b(i-1)+j ≤ L,where for (i) we use Lemma 3 for each submatrix (l(i) π l (i)⊤ π) ⊗ I d and l(i) π = (0, . . . , 0, L π b(i-1)+1 , . . . , L π bi , 0, . . . , 0) ⊤ . Similar to the case of generalized linear models, the inequality is tight when b = 1 but can be loose for other values of b.Before proceeding to the omitted proofs, we first state the following standard definitions and first-order characterization of strong convexity, for completeness. Definition 1. A function f : R d → R is said to be µ-strongly convex with parameter µ > 0, if for any x, y ∈ R d and any λ ∈ (0, 1):f (λx + (1 -λ)y) ≤ λf (x) + (1 -λ)f (y) - µ 2 λ(1 -λ)∥x -y∥ 2 2 .Lemma 4. Let f : R d → R be a continuous µ-strongly convex function with µ > 0. Then, for any x, y ∈ R d :f (y) ≥ f (x) + ⟨g x , y -x⟩ + µ 2 ∥x -y∥ 2 2 ,where g x ∈ ∂f (x), and ∂f (x) is the subdifferential of f at x.We also include the following lemma on the variance bound under without-replacement sampling, which is useful for our proof. i∈B ∇f i (x) -∇f (x) 2 2 = n -b b(n -1) E i ∥∇f i (x) -∇f (x)∥ 2 2 .Proof. We first expand the square on the left-hand side, as followsE B 1 b i∈B ∇f i (x) -∇f (x) 2 2 = 1 b 2 E B i,i ′ ∈B ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ = 1 b 2 E B i,i ′ ∈B,i̸ =i ′ ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ + 1 b E i ∥∇f i (x) -∇f (x)∥ 2 2 .Since the batch B is sampled uniformly and without replacement from [n], the probability that any pair(i, i ′ ) from [n] with i ̸ = i ′ is in B is b(b-1) n(n-1). By the linearity of expectation, we haveE B i,i ′ ∈B,i̸ =i ′ ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ = E B i,i ′ ∈[n],i̸ =i ′ 1 i,i ′ ∈B ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ = i,i ′ ∈[n],i̸ =i ′ E B 1 i,i ′ ∈B ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ = b(b -1) n(n -1) i,i ′ ∈[n],i̸ =i ′ ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ ,where 1 is the indicator function such that 1 i,i ′ ∈B = 1 if both i, i ′ ∈ B and is equal to zero otherwise.Hence, we obtainE B 1 b i∈B ∇f i (x) -∇f (x) 2 = b -1 bn(n -1) i,i ′ ∈[n],i̸ =i ′ ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ + 1 b E i ∥∇f i (x) -∇f (x)∥ 2 = b -1 bn(n -1) i,i ′ ∈[n] ⟨∇f i (x) -∇f (x), ∇f i ′ (x) -∇f (x)⟩ + n -b b(n -1) E i ∥∇f i (x) -∇f (x)∥ 2 (i) = n -b b(n -1) E i ∥∇ j f i (x) -∇ j f (x)∥ 2 ,where (i) is due to f = 1 n n i=1 f i having the finite sum structure. Now we provide the omitted proofs from Section 2.Lemma 6. Under Assumption 1, for any k ∈ [K], the iterates {yEQUATIONwhereE k := η k Gap v (x k , y * ) + b 2n ∥x * -x k ∥ 2 2 -b 2n ∥x * -x k-1 ∥ 2 2 .Proof. We first note that based on Line 6 of Alg. 1, we haveE ⊤ b y (i) , x k-1,i - b j=1 f * π (k) b(i-1)+j (y j ) = b j=1 y j , x k-1,i -f * π (k) b(i-1)+j (y j ) .Since the max problem defining y k is separable, we have for b(i -1) + 1 ≤ j ≤ bi and i ∈ [m]y j k = arg max y j ∈R d y j , x k-1,i -f * π (k) j (y j ) , which leads to x k-1,i ∈ ∂f * π (k) j (y j k ). Further, since each component function f * j is 1 Lj -strongly convex thus for b(i -1) + 1 ≤ j ≤ bi, we also have f * π (k) j (v j k ) ≥ f * π (k) j (y j k ) + x k-1,i , v j k -y j k + 1 2L π (k) j ∥v j k -y j k ∥ 2 ,which leads toL(x k , v) = 1 n m i=1 E ⊤ b v (i) k , x k-1,i - bi j=b(i-1)+1 f * π (k) j (v j k ) + 1 n m i=1 E ⊤ b v (i) k , x k -x k-1,i ≤ 1 n m i=1 E ⊤ b y (i) k , x k-1,i - bi j=b(i-1)+1 f * π (k) j (y j k ) + 1 n m i=1 E ⊤ b v (i) k , x k -x k-1,i - 1 2n ∥y k -v k ∥ 2 Λ -1 k .Using the same argument, asx * ∈ ∂f * i (y i * ) for i ∈ [n], we have f * π (k) i (y i k ) ≥ f * π (k) i (y i * ,k ) + x * , y i k -y i * ,k + 1 2L π (k) i ∥y i k -y i * ,k ∥ 2 .Thus,L(x * , y * ) = 1 n m i=1 E ⊤ b y (i) * ,k , x * - bi j=b(i-1)+1 f * π (k) j (y j * ,k ) ≥ 1 n m i=1 E ⊤ b y (i) k , x * - bi j=b(i-1)+1 f * π (k) j (y j k ) + 1 2n ∥y k -y * ,k ∥ 2 Λ -1 k = 1 n m i=1 E ⊤ b y (i) k , x * + b 2η k ∥x * -x k-1,i ∥ 2 - b 2η k ∥x * -x k-1,i ∥ 2 - bi j=b(i-1)+1 f * π (k) j (y j k ) + 1 2n ∥y k -y * ,k ∥ 2 Λ -1 k .Using the updating scheme of x k-1,i+1 and noticing that ϕi k (x) = E ⊤ b y (i) k , x + b 2η k ∥x-x k-1,i ∥ 2 is bη k -strongly convex and minimized at x k-1,i+1 , we haveE ⊤ b y (i) k , x * + b 2η k ∥x * -x k-1,i ∥ 2 ≥ E ⊤ b y (i) k , x k-1,i+1 + b 2η k ∥x k-1,i+1 -x k-1,i ∥ 2 + b 2η k ∥x k-1,i+1 -x * ∥ 2 ,which leads toL(x * , y * ) ≥ 1 n m i=1 E ⊤ b y (i) k , x k-1,i+1 + b 2η k ∥x k-1,i+1 -x k-1,i ∥ 2 - bi j=b(i-1)+1 f * π (k) j (y i k ) + b 2nη k m i=1 ∥x k-1,i+1 -x * ∥ 2 -∥x k-1,i -x * ∥ 2 + 1 2n ∥y k -y * ,k ∥ 2 Λ -1 k = 1 n m i=1 E ⊤ b y (i) k , x k-1,i+1 + b 2η k ∥x k-1,i+1 -x k-1,i ∥ 2 - bi j=b(i-1)+1 f * π (k) j (y j k ) + b 2nη k ∥x k -x * ∥ 2 -∥x k-1 -x * ∥ 2 + 1 2n ∥y k -y * ,k ∥ 2 Λ -1 k .Hence, combining the bounds on L(x k , v) and L(x * , y * ) and lettingE k := η k L(x k , v) -L(x * , y * ) + b 2n ∥x k -x * ∥ 2 - b 2n ∥x k-1 -x * ∥ 2 ,we obtainE k ≤ η k n m i=1 E ⊤ b y (i) k , x k-1,i -x k-1,i+1 + η k n m i=1 E ⊤ b v (i) k , x k -x k-1,i - η k 2n ∥y k -v k ∥ 2 Λ -1 k - η k 2n ∥y k -y * ,k ∥ 2 Λ -1 k - b 2n m i=1 ∥x k-1,i+1 -x k-1,i ∥ 2 = η k n m i=1 E ⊤ b y (i) k , x k -x k-1,i+1 + η k n m i=1 E ⊤ b v (i) k -y (i) k , x k -x k-1,i - η k 2n ∥y k -v k ∥ 2 Λ -1 k - η k 2n ∥y k -y * ,k ∥ 2 Λ -1 k - b 2n m i=1 ∥x k-1,i+1 -x k-1,i ∥ 2 ,thus completing the proof.We note that the first inner product termT 1 := η k n m i=1 E ⊤ b y (i) k , x k -x k-1,i+1in Eq. ( 5) can be cancelled by the last negative term -b 2n m i=1 ∥x k-1,i+1 -x k-1,i ∥ 2 therein, as precisely proved in Lemma 10 of Appendix C. In the following subsections, we continue our analysis and handle the remaining terms in Eq. ( 5) according to different shuffling and derive the final complexity.We introduce the following lemma to bound the second inner product term T 2 :=η k n m i=1 E ⊤ b v (i) k -y (i)k , x k -x k-1,i in Lemma 6 when there are random permutations.Lemma 7. Under Assumptions 1 and 2, for any k ∈ [K], the iterates {y(i) k } m i=1 and {x k-1,i } m+1 i=1generated by Algorithm 1 with uniformly random shuffling (RR/SO) satisfyE[T 2 ] ≤ E η 3 k n Lg π (k) Lg π (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k + η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * ,whereT 2 := η k n m i=1 E ⊤ b v (i) k -y (i) k , x k -x k-1,i . Proof. First note that x k -x k-1,i = m j=i (x k-1,j+1 -x k-1,j ) = -η k b m j=i E ⊤ b y (j) k = -η k b E ⊤ I bd(i-1)↑ y k , so we have η k n m i=1 E ⊤ b v (i) k -y (i) k , x k -x k-1,i = η k n m i=1 E ⊤ b v (i) k -y (i) k , m j=i (x k-1,j+1 -x k-1,j ) = - η 2 k bn m i=1 E ⊤ I (di) (v k -y k ), E ⊤ I bd(i-1)↑ y k = - η 2 k bn m i=1 E ⊤ I (di) (v k -y k ), E ⊤ I bd(i-1)↑ (y k -y * ,k ) I1 - η 2 k bn m i=1 E ⊤ I (di) (v k -y k ), E ⊤ I bd(i-1)↑ y * ,k I2 .For the term I 1 , we use Young's inequality with α > 0 to be set later and obtainEQUATIONFurther, notice thatEQUATIONwhere for the last inequality we use Cauchy-Schwarz inequality. Using the same argument, we can boundEQUATIONThus, combining ( 6)-( 8) and choosing α = 2η k Lg π (k) , we obtainI 1 ≤ η 3 k m Lg π (k) Lg π (k) b ∥y k -y * ,k ∥ 2 Λ -1 k + η k 4n ∥v k -y k ∥ 2 Λ -1 k .For the term I 2 , we again apply Young's inequality with β > 0 to be set later and obtainI 2 = - η 2 k bn m i=1 E ⊤ I (di) (v k -y k ), E ⊤ I bd(i-1)↑ y * ,k ≤ η 2 k β 2bn m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 + η 2 k 2bnβ m i=1 ∥E ⊤ I (di) (v k -y k )∥ 2 ≤ η 2 k β 2bn m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 + η 2 k Lg π (k) 2nβ ∥v k -y k ∥ 2 Λ -1 k .Choosing β = 2η k Lg and using the fact that Lg π (k) ≤ Lg , we haveI 2 ≤ η 3 k Lg bn m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 + η k 4n ∥v k -y k ∥ 2 Λ -1 k .Hence, combining the above two estimates with m = n/b, we haveT 2 ≤ η 3 k Lg bn m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 + η 3 k n Lg π (k) Lg π (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k .First, consider the RR scheme. Taking conditional expectation on both sides w.r.t. the randomness up to but not including k-th epoch, we haveE k [T 2 ] ≤ η 3 k Lg bn E k m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 + E k η 3 k n Lg π (k) Lg π (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k .For the first term, since the only randomness comes from the permutation π (k) , we can proceed as in the proof of Lemma 11 and obtainη 3 k Lg bn E k m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 (i) = η 3 k Lg bn m i=1 E π (k) ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 = η 3 k Lg bn m i=1 (n -b(i -1)) 2 E π (k) E ⊤ I bd(i-1)↑ y * ,k n -b(i -1) 2 (ii) ≤ η 3 k Lg bn m i=1 (n -b(i -1)) 2 b(i -1) (n -b(i -1))(n -1) σ 2 * = η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * ,where we use the linearity of expectation for (i), and (ii) is due to Lemma 5 and the definitionσ 2 * = 1 n n i=1 ∥∇f i (x * )∥ 2 = 1 n n i=1 ∥y i * ∥ 2 .Then taking expectation w.r.t. all randomness on both sides, we obtainE[T 2 ] ≤ E η 3 k n Lg π (k) Lg π (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k + η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * .Finally, we remark that the above argument for bounding the termη 3 k Lg bn E k m i=1 ∥E ⊤ I bd(i-1)↑ y * ,k ∥ 2 also applies to the SO scheme, in which case there is only one random permutation at the very beginning that induces the randomness.We state the final convergence rate and complexity in the following theorem and provide the proof for completeness. E[H K (f ( xK ) -f (x * ))] ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * .As a consequence, for any ϵ > 0, there exists a choice of a constant step size η k = η for whichE[f ( xK ) -f (x * )] ≤ ϵ after O n √ Lg Lg ∥x0-x * ∥ 2 2 ϵ + (n-b)(n+b) n(n-1) √ n Lg σ * ∥x0-x * ∥ 2 2 ϵ 3/2 gradient queries.Proof. Combining the bounds in Lemma 10 and 2 and plugging them into Eq. ( 5), we obtainE[E k ] ≤ E η 3 k n Lg π (k) Lg π (k) b 2 - η k 2n ∥y k -y * ,k ∥ 2 Λ -1 k + η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * .For the stepsizeη k such that η k ≤ b n 2 Lg π (k) Lg π (k), we haveη 3 k n Lg π (k) Lg π (k) b 2 -η k 2n ≤ 0, thus E[E k ] ≤ η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * .Using our definition of E k and telescoping from k = 1 to K, we haveE K k=1 η k Gap v (x k , y * ) ≤ b 2n ∥x * -x 0 ∥ 2 2 - b 2n E[∥x * -x K ∥ 2 2 ] + K k=1 η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * . Noticing that L(x, v) is convex in x for a fixed v, we have Gap v ( xK , y * ) ≤ K k=1 η k Gap v (x k , y * )/H K , where xK = K k=1 η k x k /H K and H K = K k=1 η k , which leads to E H K Gap v ( xK , y * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 η 3 k Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * .Further choosing v = y xK , we obtainEQUATIONTo analyze the individual gradient oracle complexity, we choose constant stepsizes η≤ b n √ 2 Lg Lg, then Eq. ( 9) will becomeE[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * .Without loss of generality, we assume that b ̸ = n, otherwise the method and its analysis reduce to (full) gradient descent. We consider the following two cases:• "Small K" case: if η = b n √ 2 Lg Lg ≤ 3b 3 (n-1)∥x0-x * ∥ 2 2 n(n-b)(n+b) Lg Kσ 2 * 1/3, we haveE[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * ≤ Lg Lg √ 2K ∥x 0 -x * ∥ 2 2 + 1 2 (n -b)(n + b) n 2 (n -1) 1/3 ( Lg ) 1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 3 1/3 K 2/3 . • "Large K" case: if η = 3b 3 (n-1)∥x0-x * ∥ 2 2 n(n-b)(n+b) Lg Kσ 2 * 1/3 ≤ b n √ 2 Lg Lg, we haveE[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 Lg (n -b)(n + b) 6b 2 (n -1) σ 2 * ≤ (n -b)(n + b) n 2 (n -1) 1/3 ( Lg ) 1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 3 1/3 K 2/3 .Combining these two cases by setting η = min , we obtainE[f ( xK ) -f (x * )] ≤ Lg Lg √ 2K ∥x 0 -x * ∥ 2 2 + (n -b)(n + b) n 2 (n -1) 1/3 ( Lg ) 1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 3 1/3 K 2/3 .Hence, to guarantee E[f ( xK ) -f (x * )] ≤ ϵ for ϵ > 0, the total number of individual gradient evaluations will benK ≥ max n 2 Lg Lg ∥x 0 -x * ∥ 2 2 ϵ , (n -b)(n + b) n -1 1/2 2 3/2 ( Lg ) 1/2 σ * ∥x 0 -x * ∥ 2 2 3 1/2 ϵ 3/2 ,as claimed.In this subsection, we provide the convergence results for incremental gradient descent which does not involve random permutations. We first prove the technical lemma below to bound the term 5) of Lemma 6.T 2 := η k n m i=1 E ⊤ b v (i) k -y (i) k , x k -x k-1,i in Eq. (Lemma 8. For any k ∈ [K], the iterates {y(i)k } m i=1 and {x k-1,i } m+1 i=1 generated by Algorithm 1 with fixed data ordering satisfyEQUATIONProof. Proceeding as in the proof of Lemma 7, we haveη k n m i=1 E ⊤ b v (i) -y (i) k , x k -x k-1,i = η k n m i=1 E ⊤ b v (i) -y (i) k , m j=i (x k-1,j+1 -x k-1,j ) = - η 2 k bn m i=1 E ⊤ I (di) (v -y k ), E ⊤ I bd(i-1)↑ y k = - η 2 k bn m i=1 E ⊤ I (di) (v -y k ), E ⊤ I bd(i-1)↑ (y k -y * ) I1 - η 2 k bn m i=1 E ⊤ I (di) (v -y k ), E ⊤ I bd(i-1)↑ y * I2 .For both terms I 1 and I 2 , we apply Young's inequality with α = 2η k Lg 0 and obtainEQUATIONandEQUATIONWe now show that the termη 3 k Lg 0 nb m i=1 ∥E ⊤ I bd(i-1)↑ y * ∥ 2 2 is no larger than either η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 or η 3 k (n-b) 2 b 2 Lg 0 σ 2 * . This is trivial when b = n as E ⊤ I 0↑ y * = n i=1 y i * = 0. When b < n, to show the former one, we have m i=1 ∥E ⊤ I bd(i-1)↑ y * ∥ 2 2 ≤ Λ 1/2 m i=1 I bd(i-1)↑ EE ⊤ I bd(i-1)↑ Λ 1/2 2 ∥y * ∥ 2 Λ -1 = mn Lg 0 ∥y * ∥ 2 Λ -1 = n 2 b Lg 0 ∥y * ∥ 2 Λ -1 .To prove the latter one, we notice that ∥y j * ∥ 2 2 ≤ b(m -1) m-1 i=1 bi j=1 ∥y j * ∥ 2 2 = b(m -1) m-1 i=1 bi j=b(i-1)+1 (m -i)∥y j * ∥ 2 2 ≤ b(m -1) 2 (m-1)b i=1 ∥y i * ∥ 2 2 .Further noticing that(m-1)b i=1 ∥y i * ∥ 2 2 ≤ n i=1 ∥y i * ∥ 2 = nσ 2 * , we have η 3 k Lg 0 nb m i=1 ∥E ⊤ I bd(i-1)↑ y * ∥ 2 2 ≤ η 3 k Lg 0 nb b(m -1) 2 nσ 2 * = η 3 k Lg 0 (n -b) 2 b 2 σ 2 * .The same bound also captures the case b = n and leads toEQUATIONHence, combining Eq. ( 11)-( 13), we obtainI 2 ≤ η 3 k n b 2 Lg 0 Lg 0 ∥y k -y * ∥ 2 Λ -1 + η k 2n ∥v -y k ∥ 2 Λ -1 + min η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 Lg 0 σ 2 * ,which finishes the proof.We are now ready to state our convergence results for IGD in the following theorem, with its proof provided for completeness. H K f ( xK ) -f (x * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 min η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 Lg 0 σ 2 * .As a consequence, for any ϵ > 0, there exists a choice of a constant step size η k = η such thatf ( xK ) -f (x * ) ≤ ϵ after O n √ Lg 0 Lg 0 ∥x0-x * ∥ 2 2 ϵ + min √ n Lg 0 Lg 0 ∥y * ∥ Λ -1 , (n-b) √ Lg 0 σ * ∥x0-x * ∥ 2 2 ϵ 3/2gradient queries.Proof. Combining the bounds in Lemma 10 and 8 and plugging them into Eq. ( 5) in Lemma 6 without random permutations, we haveE k ≤ η 3 k n Lg 0 Lg 0 b 2 - η k 2n ∥y k -y * ∥ 2 Λ -1 + min η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 Lg 0 σ 2 * . If η k ≤ b n √ 2 Lg 0 Lg 0 , we have η 3 k n Lg 0 Lg 0 b 2 -η k 2n ≤ 0, thus E k ≤ min η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 Lg 0 σ 2 * .Using the definition of E k and telescoping from k = 1 to K, we obtainK k=1 η k Gap v (x k , y * ) ≤ b 2n ∥x * -x 0 ∥ 2 2 - b 2n ∥x * -x K ∥ 2 2 + K k=1 min η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 Lg 0 σ 2 * . Noticing that L(x, v) is convex w.r.t. x, we have Gap v ( xK , y * ) ≤ K k=1 η k Gap v (x k , y * )/H K , where xK = K k=1 η k x k /H K and H K = K k=1 η k , so we obtain H K Gap v ( xK , y * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 min η 3 k n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 Lg 0 σ 2 * ,Further choosing v = y xK , we obtainEQUATIONTo analyze the individual gradient oracle complexity, we choose constant stepsizes η ≤ b n √ 2 Lg 0 Lg 0 and assume b < n without loss of generality, then Eq. ( 14) becomesf ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + min η 2 n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 , η 2 (n -b) 2 b 2 Lg 0 σ 2 * . When Lg 0 ∥y * ∥ 2 Λ -1 ≤ (n-b) 2 n σ 2 * , we set η = min b n √ 2 Lg 0 Lg 0 , b 3 ∥x0-x * ∥ 2 2 2n 2 Lg 0 Lg 0 K∥y * ∥ 2 Λ -1 1/3and consider the following two possible cases:• "Small K" case: if η = b n √ 2 Lg 0 Lg 0 ≤ b 3 ∥x0-x * ∥ 2 2 2n 2 Lg 0 Lg 0 K∥y * ∥ 2 Λ -1 1/3, we havef ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 ≤ Lg 0 Lg 0 √ 2K ∥x 0 -x * ∥ 2 2 + Lg 0 Lg 0 1/3 ∥y * ∥ 2/3 Λ -1 ∥x 0 -x * ∥ 4/3 2 2 2/3 n 1/3 K 2/3 . • "Large K" case: if η = b 3 ∥x0-x * ∥ 2 2 2n 2 Lg 0 Lg 0 K∥y * ∥ 2 Λ -1 1/3 ≤ b √ 2 Lg 0 Lg 0 , we have f ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 n b 2 Lg 0 Lg 0 ∥y * ∥ 2 Λ -1 ≤ 2 1/3 Lg 0 Lg 0 1/3 ∥y * ∥ 2/3 Λ -1 ∥x 0 -x * ∥ 4/3 2 n 1/3 K 2/3 .Combining these two cases, we havef ( xK ) -f (x * ) ≤ Lg 0 Lg 0 √ 2K ∥x 0 -x * ∥ 2 2 + 2 1/3 Lg 0 Lg 0 1/3 ∥y * ∥ 2/3 Λ -1 ∥x 0 -x * ∥ 4/3 2 n 1/3 K 2/3 .Hence, to guarantee E[f ( xK ) -f (x * )] ≤ ϵ for ϵ > 0, the total number of required individual gradient evaluations will beEQUATION)When (n-b) 2 n σ 2 * ≤ Lg 0 ∥y * ∥ 2 Λ -1 , we set η = min b n √ 2 Lg 0 Lg 0 , b 3 ∥x0-x * ∥ 2 2 2n(n-b) 2 Lg 0 Kσ 2 * 1/3and consider the two cases as below:• "Small K" case: if η = b n √ 2 Lg 0 Lg 0 ≤ b 3 ∥x0-x * ∥ 2 2 2n(n-b) 2 Lg 0 Kσ 2 * 1/3, we havef ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 (n -b) 2 b 2 Lg 0 σ 2 * ≤ Lg 0 Lg 0 √ 2K ∥x 0 -x * ∥ 2 2 + (n -b) 2/3 ( Lg ) 1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 2 2/3 n 2/3 K 2/3 . • "Large K" case: if η = b 3 ∥x0-x * ∥ 2 2 2n(n-b) 2 Lg 0 Kσ 2 * 1/3 ≤ b n √ 2 Lg 0 Lg 0 , we have f ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 (n -b) 2 b 2 Lg 0 σ 2 * ≤ 2 1/3 (n -b) 2/3 ( Lg ) 1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 n 2/3 K 2/3 .Combining these two cases, we obtainf ( xK ) -f (x * ) ≤ Lg 0 Lg 0 √ 2K ∥x 0 -x * ∥ 2 2 + 2 1/3 (n -b) 2/3 ( Lg ) 1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 n 2/3 K 2/3 .To guarantee E[f ( xK ) -f (x * )] ≤ ϵ for ϵ > 0, the total number of required individual gradient evaluations will beEQUATION)Combining Eq. ( 15) and Eq. ( 16), we finally havenK ≥ n 2 Lg 0 Lg 0 ∥x 0 -x * ∥ 2 2 ϵ + min 4n 1/2 Lg 0 Lg 0 1/2 ∥y * ∥ Λ -1 ∥x 0 -x * ∥ 2 2 ϵ 3/2 , 4(n -b)( Lg 0 ) 1/2 σ * ∥x 0 -x * ∥ 2 2 ϵ 3/2 ,thus finishing the proof.Before proceeding to the omitted proofs for the smooth convex settings in finite-sum with linear predictors, we first recall its primal-dual reformulation, then state the specialized version of a primaldual shuffled SGD algorithm in Algorithm 2. Recall that (PL) admits an explicit reformulation using convex conjugates of ℓ i :min x∈R d max y∈R n L(x, y) := 1 n ⟨Ax, y⟩ - 1 n n i=1 ℓ * i (y i ) = 1 n n i=1 a ⊤ i xy i -ℓ * i (y i ) (PL-PD)where y i x = arg max y i ∈R {y i a ⊤ i x -ℓ * i (y i )} (different from the general smooth convex finite-sum settings in Section 2 and Appendix B). Further, for notational convenience, we assume that the partition is ordered, in the sense that for 1 ≤ j < j ′ ≤ m, max i∈S j i < min i ′ ∈S j ′ i ′ . 2 We denote by y (j) the subvector of y ∈ R n indexed by the elements of S j , and by A (j) the submatrix obtained from A ∈ R n×d by selecting the rows indexed by S j .Based on the formulation (PL-PD), we view shuffled SGD as a primal-dual method with block coordinate updates on the dual side, as summarized in Algorithm 2, for completeness. To see the equivalence, in i-th inner iteration of k-th epoch, we first update the i-th block y(i) k ∈ R b of the dual vector y k-1 ∈ R n based on x k-1,ias in Line 6. Since the dual update has a decomposable structure, this maximization step corresponds to computing the (sub)gradients {ℓ ′π (k) j (a ⊤ π (k) j x k-1,i )} bi j=b(i-1)+1at x k-1,i for the batch of individual losses indexed by {π (k) j } bi j=b(i-1)+1 . Then in Line 7, we perform a minimization step using y (i) k to compute x k-1,i+1 on the primal side. Combining these two steps, we havex k-1,i+1 = x k-1,i -η k b bi j=b(i-1)+1 ℓ ′ π (k) j (a ⊤ π (k) j x k-1,i )a π (k) j, which is exactly the original primal shuffled SGD updating scheme. x k-1,1 = x k-1for i = 1 to m do 6:y (i) k = arg max y∈R b y ⊤ A (i) k x k-1,i - b j=1 ℓ * π (k) b(i-1)+j (y j ) 7: x k-1,i+1 = arg min x∈R d y (i)⊤ k A (i) k x + b 2η k ∥x -x k-1,i ∥ 2 8:end for 9: EQUATIONEQUATIONProof. By Line 6 in Alg. 2, we have y(i) k = arg max y∈R b y ⊤ A (i) k x k-1,i - b j=1 ℓ * π (k) b(i-1)+j (y j ) for i ∈ [m]. Notice that since y ⊤ A (i) k x k-1,i - b j=1 ℓ * π (k) b(i-1)+j (y j ) = b j=1 y j a ⊤ π (k) b(i-1)+j x k-1,i -ℓ * π (k) b(i-1)+j (y j )is separable, we havey j k = arg max y∈R {ya ⊤ π (k) j x k-1,i -ℓ * π (k) j (y)} for b(i -1) + 1 ≤ j ≤ bi, thus a ⊤ π (k) j x k-1,i ∈ ∂ℓ * π (k) j (y j k ). Since ℓ * i is 1Li -strongly convex by Assumption 3, then by Lemma 4 we obtain for b(i -1) + 1 ≤ j ≤ bi ℓ * π (k) j v j k ≥ ℓ * π (k) j (y j k ) + a ⊤ π (k) j x k-1,i v j k -y j k + 1 2L π (k) j v j k -y j k 2 ,which leads toEQUATION)Using the same argument for L(x * , y * ) as a ⊤ j x * ∈ ∂ℓ * j (y j * ) for j ∈ [n], we haveEQUATION)Adding and substracting the term b2nη k m i=1 ∥x * -x k-1,i ∥ 22 on the R.H.S. of Eq. ( 19), we obtainL(x * , y * ) ≥ 1 n m i=1 y (i)⊤ k A (i) k x * + b 2η k ∥x * -x k-1,i ∥ 2 2 - bi j=b(i-1)+1 ℓ * π (k) j (y j k ) - b 2nη k m i=1 ∥x * -x k-1,i ∥ 2 2 + 1 2n ∥y k -y * ,k ∥ 2 Λ -1 k .By Line 7 of Alg. 2, we havex k-1,i+1 = arg min x∈R d y (i)⊤ k A (i) k x + b 2η k ∥x -x k-1,i ∥ 2 2 . Fur- ther noticing that ϕ (i) k (x) := y (i)⊤ k A (i) k x + b 2η k ∥x -x k-1,i ∥ 2 2 is b η k -strongly convex w.r.t.x and ∇ϕ(i) k (x k-1,i+1 ) = 0, we have ϕ (i) k (x * ) ≥ ϕ (i) k (x k-1,i+1 ) + b 2η k ∥x * -x k-1,i+1 ∥ 2 2 ,which leads toEQUATIONwhere we telescope from i = 1 to m for the termm i=1 ∥x * -x k-1,i+1 ∥ 2 2 -∥x * -x k-1,i ∥ 2 2, and use the definitions that x k = x k-1,m+1 and x k-1 = x k-1,1 for (i).Combining the bounds from Eq. ( 18) and Eq. ( 20) and denotingE k := η k L(x k , v) -L(x * , y * ) + b 2n ∥x * -x k ∥ 2 2 - b 2n ∥x * -x k-1 ∥ 2 2 ,we obtainE k ≤ η k n m i=1 y (i)⊤ k A (i) k (x k-1,i -x k-1,i+1 ) + η k n m i=1 v (i)⊤ k A (i) k (x k -x k-1,i ) - η k 2n ∥y k -v k ∥ 2 Λ -1 k - η k 2n ∥y k -y * ,k ∥ 2 Λ -1 k - b 2n m i=1 ∥x k-1,i -x k-1,i+1 ∥ 2 2 = η k n m i=1 y (i)⊤ k A (i) k (x k -x k-1,i+1 ) + η k n m i=1 (v (i) k -y (i) k ) ⊤ A (i) k (x k -x k-1,i ) - η k 2n ∥y k -v k ∥ 2 Λ -1 k - η k 2n ∥y k -y * ,k ∥ 2 Λ -1 k - b 2n m i=1 ∥x k-1,i -x k-1,i+1 ∥ 2 2 ,thus completing the proof.Lemma 10. For any k ∈ [K], the iterates {y(i) k } m i=1 and {x k-1,i } m+1 i=1 in Algorithm 2 satisfy T 1 = b 2n m i=1 ∥x k-1,i -x k-1,i+1 ∥ 2 2 - b 2n ∥x k-1 -x k ∥ 2 2 .Proof. By Line 7 in Alg. 2, we haveA (i)⊤ k y (i) k = b η k (x k-1,i -x k-1,i+1 ). Further noticing that x k -x k-1,i+1 = - m j=i+1 (x k-1,j -x k-1,j+1), we obtainT 1 := η k n m i=1 y (i)⊤ k A (i) k (x k -x k-1,i+1 ) = - b n m-1 i=1 m j=i+1 ⟨x k-1,i -x k-1,i+1 , x k-1,j -x k-1,j+1 ⟩ = b 2n m i=1 ∥x k-1,i -x k-1,i+1 ∥ 2 - b 2n m i=1 (x k-1,i -x k-1,i+1 ) 2 ,thus completing the proof.Lemma 11. Under Assumption 4, for any k ∈ [K], the iterates {y(i) k } m i=1 and {x k-1,i } m+1 i=1generated by Algorithm 2 with uniformly random shuffling satisfyE[T 2 ] ≤ E η 3 k n Lπ (k) Lπ (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k + η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .Proof. By Line 7 in Alg. 2, we havex k-1,i -x k-1,i+1 = η k b A (i)⊤ k y (i)k . Using the definition of I j↑ for 0 ≤ j ≤ n -1 as in Section 1, we obtainx k -x k-1,i = - m j=i (x k-1,j -x k-1,j+1 ) = - η k b m j=i A (j)⊤ k y (j) k = - η k b A k I b(i-1)↑ y k .Also, we haveA (i)⊤ k (v (i) k -y (i) k ) = A k I (i) (v k -y k) by the definition of I (i) in Section 3. Combining these two observations, we haveEQUATIONwhere we make a decomposition w.r.t. y * ,k in (i). For the first term in Eq. ( 21), we use Young's inequality for α > 0 and haveEQUATION)Expanding the squares and rearranging the terms in Eq. ( 23), we haveEQUATIONwhere we use Cauchy-Schwarz inequality for (i). Using a similar argument, we also haveη 2 k 2bnα m i=1 ∥A ⊤ k I (i) (v k -y k )∥ 2 2 ≤ η 2 k 2bnα Λ 1/2 k m i=1 I (i) A k A ⊤ k I (i) Λ 1/2 k 2 ∥v k -y k ∥ 2 Λ -1 k .By the definitions of Lπ (k) and Lπ (k) , and choosing α = 2η k Lπ (k) in Eq. ( 23), we obtainEQUATIONFor the second term in Eq. ( 22), we apply Young's inequality with β > 0 and proceed as above:- η 2 k bn m i=1 A ⊤ k I b(i-1)↑ y * ,k , A ⊤ k I (i) (v k -y k ) ≤ η 2 k β 2bn m i=1 ∥A ⊤ k I b(i-1)↑ y * ,k ∥ 2 2 + η 2 k 2bnβ m i=1 ∥A ⊤ k I (i) (v k -y k )∥ 2 2 ≤ η 2 k β 2bn m i=1 ∥A ⊤ k I b(i-1)↑ y * ,k ∥ 2 2 + η 2 k 2nβ Lπ (k) ∥v k -y k ∥ 2 Λ -1 k .Noticing that Lπ (k) ≤ L, we choose β = 2η k L and obtainEQUATIONCombining Eq. ( 25) and Eq. ( 26), we haveEQUATIONWe first assume the RR scheme. Taking conditional expectation w.r.t. the randomness up to but not including k-th epoch, we haveE k [T 2 ] ≤ η 3 k L nb E k m i=1 ∥A ⊤ k I b(i-1)↑ y * ,k ∥ 2 2 + E k η 3 k n Lπ (k) Lπ (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k .For the first termη 3 k L nb E k m i=1 ∥A ⊤ k I b(i-1)↑ y * ,k ∥ 2 2, the only randomness is from the random permutation π (k) . In this case, each term A ⊤ k I b(i-1)↑ y * ,k can be considered as a sum of a batch sampled without replacement from {y j * a j } j∈[n] , while n j=1 y j * a j = 0 as x * is the minimizer, we then can use Lemma 5 and obtainEQUATION≤ η 3 k L nb m i=1 n -b(i -1) 2 b(i -1) n -b(i -1) (n -1) σ 2 * = η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * ,where (i) is due to the linearity of expectation, and we use our definitionσ 2 * = 1 n n j=1 (y j * ) 2 ∥a j ∥ 2 2 = E j ∥y j * a j ∥ 2 2for (ii). Taking expectation w.r.t. all the randomness on both sides and using the law of total expectation, obtainE[T 2 ] ≤ E η 3 k n Lπ (k) Lπ (k) b 2 ∥y k -y * ,k ∥ 2 Λ -1 k + η k 2n ∥v k -y k ∥ 2 Λ -1 k + η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .For the SO scheme, since there is only one random permutation generated at the very beginning, we can take expectation w.r.t. all the randomness on both sides of ( 27), and the randomness for the termη 3 k L nb E m i=1 ∥A ⊤ k I b(i-1)↑ y * ,k ∥ 2 2is only from the initial random permutation. So the above argument still applies to this case, and we complete the proof. E[H K (f ( xK ) -f (x * ))] ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .As a result, given ϵ > 0, there exists a constant step sizeη k = η such that E[f ( xK ) -f (x * )] ≤ ϵ after O n √ L L∥x0-x * ∥ 2 2 ϵ + (n-b)(n+b) n(n-1) √ n Lσ * ∥x0-x * ∥ 2 2 ϵ 3/2individual gradient queries.Proof. Combining the bounds in Lemma 10 and 11 and plugging them into Eq. ( 17), we obtainE[E k ] ≤ E η 3 k n Lπ (k) Lπ (k) b 2 - η k 2n ∥y k -y * ,k ∥ 2 Λ -1 k + η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .For the stepsizeη k such that η k ≤ b n 2 Lπ (k) Lπ (k), we haveη 3 k n Lπ (k) Lπ (k) b 2 -η k 2n ≤ 0, thus E[E k ] ≤ η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .Noticing thatE k = η k Gap v (x k , y * ) + b 2n ∥x * -x k ∥ 2 2 -b 2n ∥x * -x k-1 ∥ 2 2and telescoping from k = 1 to K, we haveE K k=1 η k Gap v (x k , y * ) ≤ b 2n ∥x * -x 0 ∥ 2 2 - b 2n E[∥x * -x K ∥ 2 2 ] + K k=1 η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * . Noticing that L(x, v) is convex w.r.t. x, we have Gap v ( xK , y * ) ≤ K k=1 η k Gap v (x k , y * )/H K , where xK = K k=1 η k x k /H K and H K = K k=1 η k , which leads to E H K Gap v ( xK , y * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 η 3 k L(n -b)(n + b) 6b 2 (n -1) σ 2 * .Further choosing v = y xK , we obtainEQUATIONTo analyze the individual gradient oracle complexity, we choose constant stepsizes η ≤ b n √ 2 L L , then Eq. ( 28) will becomeE[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 L(n -b)(n + b) 6b 2 (n -1) σ 2 * .Without loss of generality, we assume that b ̸ = n, otherwise the method and its analysis reduce to (full) gradient descent. We consider the following two cases:• "Small K" case: if η = b n √ 2 L L ≤ 3b 3 (n-1)∥x0-x * ∥ 2 2 n(n-b)(n+b) LKσ 2 * 1/3, we haveE[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 L(n -b)(n + b) 6b 2 (n -1) σ 2 * ≤ L L √ 2K ∥x 0 -x * ∥ 2 2 + 1 2 (n -b)(n + b) n 2 (n -1) 1/3 L1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 3 1/3 K 2/3 . • "Large K" case: if η = 3b 3 (n-1)∥x0-x * ∥ 2 2 n(n-b)(n+b) LKσ 2 * 1/3 ≤ b n √ 2 L L , we have E[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 L(n -b)(n + b) 6b 2 (n -1) σ 2 * ≤ (n -b)(n + b) n 2 (n -1) 1/3 L1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 3 1/3 K 2/3 .η = min b n √ 2 L L , 3b 3 (n-1)∥x0-x * ∥ 2 2 n(n-b)(n+b) LKσ 2 * 1/3, we obtainE[f ( xK ) -f (x * )] ≤ L L √ 2K ∥x 0 -x * ∥ 2 2 + (n -b)(n + b) n 2 (n -1) 1/3 L1/3 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 3 1/3 K 2/3 .Hence, to guarantee E[f ( xK ) -f (x * )] ≤ ϵ for ϵ > 0, the total number of individual gradient evaluations will benK ≥ max n 2 L L∥x 0 -x * ∥ 2 2 ϵ , (n -b)(n + b) n -1 1/2 2 3/2 L1/2 σ * ∥x 0 -x * ∥ 2 2 3 1/2 ϵ 3/2 ,as claimed.We now provide the proof for convergence of IGD in the smooth convex settings. We first prove the following technical lemma, which bounds the inner product termT 2 := η k n m i=1 v (i) - y (i) k ⊤ A (i) (x k -x k-1,i ) without random permutations involved. Lemma 12. For any k ∈ [K], the iterates {y (i) k } mi=1 and {x k-1,i } m+1 i=1 generated by Algorithm 2 with fixed data ordering satisfyEQUATIONProof. Proceeding as in Lemma 11, we haveEQUATIONEQUATIONFor both terms in Eq. ( 30) and Eq. ( 31), we use Young's inequality for α = 2η k L0 > 0 and proceed as in Eq. ( 24) to obtainEQUATIONandEQUATIONwhere again we used α = 2η k L0 . We then prove the term 33) is no larger than the minimum ofη 3 k L0 nb m i=1 ∥A ⊤ I b(i-1)↑ y * ∥ 2 2 in Eq. (η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 and η 3 k (n-b) 2 b 2 L0 σ 2 * .Note that when b = n, we have A ⊤ I (0)↑ y * = 0, so this term disappears. When b < n, the former one can be derived as in Eq.( 24), which givesm i=1 ∥A ⊤ I b(i-1)↑ y * ∥ 2 2 ≤ Λ 1/2 m i=1 I b(i-1)↑ AA ⊤ I b(i-1)↑ Λ 1/2 2 ∥y * ∥ 2 Λ -1 = mn L0 ∥y * ∥ 2 Λ -1 = n 2 b L0 ∥y * ∥ 2 Λ -1 .For the latter one, we notice that (m -i)∥y j * a j ∥ 2 2 ≤ b(m -1) 2 (m-1)b i=1 ∥y j * a j ∥ 2 2 .By the definition thatEQUATIONNote that the bound in Eq. ( 34) equals to zero when b = n, which recovers the case of full gradient descent, so we haveEQUATIONCombining Eq. ( 32)- (35) , we obtainT 2 ≤ η 3 k n b 2 L0 L0 ∥y k -y * ∥ 2 Λ -1 + η k 2n ∥v -y k ∥ 2 Λ -1 + min η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 L0 σ 2 * ,thus finishing the proof. H K f ( xK ) -f (x * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 min η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 L0 σ 2 * .As a consequence, given ϵ > 0, there exists a constant step size η k = η such that f ( xK ) -f (x * ) ≤ ϵ after the number of gradient queries bounded byO n √ L0 L0∥x0-x * ∥ 2 2 ϵ + min √ n L0 L0∥y * ∥ Λ -1 , (n-b) √ L0σ * ∥x0-x * ∥ 2 2 ϵ 3/2 .Proof. Proceeding as in Lemmas 9 and 10, but without random permutations, we haveEQUATION)Using the bound in Lemma 12 and applying Eq. ( 29) into Eq. ( 36), we obtainE k ≤ η 3 k n L0 L0 b 2 - η k 2n ∥y k -y * ∥ 2 Λ -1 + min η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 L0 σ 2 * . If η k ≤ b n √ 2 L0 L0, we haveη 3 k n L0 L0 b 2 -η k 2n ≤ 0, thus E k ≤ min η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 L0 σ 2 * .Noticing thatE k = η k Gap v (x k , y * ) + b 2n ∥x * -x k ∥ 2 2 -b 2n ∥x * -x k-1 ∥ 2 2and telescoping from k = 1 to K, we haveK k=1 η k Gap v (x k , y * ) ≤ b 2n ∥x * -x 0 ∥ 2 2 - b 2n ∥x * -x K ∥ 2 2 + K k=1 min η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 L0 σ 2 * .Noticing that L(x, v) is convex w.r. H K Gap v ( xK , y * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 min η 3 k n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 3 k (n -b) 2 b 2 L0 σ 2 * ,Further choosing v = y xK , we obtain EQUATIONf ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + min η 2 n b 2 L0 L0 ∥y * ∥ 2 Λ -1 , η 2 (n -b) 2 b 2 L0 σ 2 * . When L0 ∥y * ∥ 2 Λ -1 ≤ (n-b) 2 n σ 2 * , we set η = min b n √ 2 L0 L0 , b 3 ∥x0-x * ∥ 2 2 2n 2 L0 L0K∥y * ∥ 2 Λ -1 1/3and consider the following two possible cases:• "Small K" case: if η = b n √ 2 L0 L0 ≤ b 3 ∥x0-x * ∥ 2 2 2n 2 L0 L0K∥y * ∥ 2 Λ -1 1/3, we havef ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 n b 2 L0 L0 ∥y * ∥ 2 Λ -1 ≤ L0 L0 √ 2K ∥x 0 -x * ∥ 2 2 + L1/3 0 L1/3 0 ∥y * ∥ 2/3 Λ -1 ∥x 0 -x * ∥ 4/3 2 2 2/3 n 1/3 K 2/3 . • "Large K" case: if η = b 3 ∥x0-x * ∥ 2 2 2n 2 L0 L0K∥y * ∥ 2 Λ -1 1/3 ≤ b √ 2 L0 L0, we havef ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 n b 2 L0 L0 ∥y * ∥ 2 Λ -1 ≤ 2 1/3 L1/3 0 L1/3 0 ∥y * ∥ 2/3 Λ -1 ∥x 0 -x * ∥ 4/3 2 n 1/3 K 2/3 .Combining these two cases, we havef ( xK ) -f (x * ) ≤ L0 L0 √ 2K ∥x 0 -x * ∥ 2 2 + 2 1/3 L1/3 0 L1/3 0 ∥y * ∥ 2/3 Λ -1 ∥x 0 -x * ∥ 4/3 2 n 1/3 K 2/3 .Hence, to guarantee E[f ( xK ) -f (x * )] ≤ ϵ for ϵ > 0, the total number of individual gradient evaluations will beEQUATION)When (n-b) 2 n σ 2 * ≤ L0 ∥y * ∥ 2 Λ -1 , we set η = min b n √ 2 L0 L0 , b 3 ∥x0-x * ∥ 2 2 2n(n-b) 2 L0Kσ 2 * 1/3and consider the two cases as below:• "Small K" case: if η = b n √ 2 L0 L0 ≤ b 3 ∥x0-x * ∥ 2 2 2n(n-b) 2 L0Kσ 2 * 1/3, we havef ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 (n -b) 2 b 2 L0 σ 2 * ≤ L0 L0 √ 2K ∥x 0 -x * ∥ 2 2 + (n -b) 2/3 L1/3 0 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 2 2/3 n 2/3 K 2/3 . • "Large K" case: if η = b 3 ∥x0-x * ∥ 2 2 2n(n-b) 2 L0Kσ 2 * 1/3 ≤ b n √ 2 L0 L0, we havef ( xK ) -f (x * ) ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + η 2 (n -b) 2 b 2 L0 σ 2 * ≤ 2 1/3 (n -b) 2/3 L1/3 0 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 n 2/3 K 2/3 .Combining these two cases, we obtainf ( xK ) -f (x * ) ≤ L0 L0 √ 2K ∥x 0 -x * ∥ 2 2 + 2 1/3 (n -b) 2/3 L1/3 0 σ 2/3 * ∥x 0 -x * ∥ 4/3 2 n 2/3 K 2/3 .To guarantee E[f ( xK ) -f (x * )] ≤ ϵ for ϵ > 0, the total number of individual gradient evaluations will beEQUATION)Combining Eq. (38) and Eq. ( 39), we finally havenK ≥ n 2 L0 L0 ∥x 0 -x * ∥ 2 2 ϵ + min 4n 1/2 L1/2 0 L1/2 0 ∥y * ∥ Λ -1 ∥x 0 -x * ∥ 2 2 ϵ 3/2 , 4(n -b) L1/2 0 σ * ∥x 0 -x * ∥ 2 2 ϵ 3/2 ,thus finishing the proof.Before we prove Theorem 3 in convex Lipschitz settings, for completeness, we first recall the following standard first-order characterization of convexity. Lemma 13. Let f : R d → R be a continuous convex function. Then, for any x, y ∈ R d :f (y) ≥ f (x) + ⟨g x , y -x⟩ ,where g x ∈ ∂f (x), and ∂f (x) is the subdifferential of f at x.The following technical lemma provides a primal-dual gap bound in convex nonsmooth settings.Lemma 14. For any k ∈ [K], the iterates {y(i)k } m i=1 and {x k-1,i } m+1 i=1 generated by Algorithm 2 satisfyEQUATIONwhereE k := η k L(x k , v) -L(x * , y * ) + b 2n ∥x * -x k ∥ 2 2 -b 2n ∥x * -x k-1 ∥ 2 2 .Proof. By the same argument as in the proof for Lemma 9, we know that a ⊤π (k) j x k-1,i ∈ ∂ℓ * π (k) j(y j k ) for b(i -1) + 1 ≤ j ≤ bi, then by Lemma 13 we haveℓ * π (k) j v j k ≥ ℓ * π (k) j (y j k ) + a ⊤ π (k) j x k-1,i v j k -y j k ,which leads toEQUATIONUsing the same argument for L(x * , y * ) as a ⊤ j x * ∈ ∂ℓ * j (y j * ) for j ∈ [n], we haveEQUATIONAdding and substracting the term b2nη k m i=1 ∥x * -x k-1,i ∥ 22 on the R.H.S. of Eq. ( 42), we obtainL(x * , y * ) ≥ 1 n m i=1 y (i)⊤ k A (i) k x * + b 2η k ∥x * -x k-1,i ∥ 2 2 - bi j=b(i-1)+1 ℓ * π (k) j (y j k ) - b 2nη k m i=1 ∥x * -x k-1,i ∥ 2 2 .Denote ϕ(i) k (x) := y (i)⊤ k A (i) k x + b 2η k ∥x -x k-1,i ∥ 2 2 , which is b η k -strongly convex w.r.t. x. Noticing that x k-1,i+1 = arg min x∈R d y (i)⊤ k A (i) k x + b 2η k ∥x -x k-1,i ∥ 2 by Line 7 of Alg. 2, we have ∇ϕ (i) k (x k-1,i+1 ) = 0, which leads to ϕ (i) k (x * ) ≥ ϕ (i) k (x k-1,i+1 ) + b 2η k ∥x * -x k-1,i+1 ∥ 2 2 .Thus, we obtainEQUATIONwhere (i) is by telescopingm i=1 ∥x * -x k-1,i+1 ∥ 2 2 -∥x * -x k-1,i ∥ 2 2 and using x k = x k-1,m+1 and x k-1 = x k-1,1, which both hold by definition.Combining the bounds from Eq. ( 41) and Eq. (43) , and denotingE k := η k L(x k , v) -L(x * , y * ) + b 2n ∥x * -x k ∥ 2 2 - b 2n ∥x * -x k-1 ∥ 2 2 ,we finally obtainE k ≤ η k n m i=1 y (i)⊤ k A (i) k (x k-1,i -x k-1,i+1 ) + η k n m i=1 v (i)⊤ k A (i) k (x k -x k-1,i ) - b 2n m i=1 ∥x k-1,i -x k-1,i+1 ∥ 2 2 = η k n m i=1 y (i)⊤ k A (i) k (x k -x k-1,i+1 ) + η k n m i=1 (v (i) k -y (i) k ) ⊤ A (i) k (x k -x k-1,i ) - b 2n m i=1 ∥x k-1,i -x k-1,i+1 ∥ 2 2 ,thus completing the proof.Note that we can still use Lemma 10 to bound the first inner product term in Eq. ( 40), as we are studying the same algorithm. The following lemma provides a bound on the second inner product term T EQUATIONEQUATION)Proof. Proceeding as in Lemma 11, we haveT 2 := η k n m i=1 v (i) k -y (i) k ⊤ A (i) k (x k -x k-1,i ) = - η 2 k bn m i=1 A ⊤ k I b(i-1)↑ y k , A ⊤ k I (i) (v k -y k ) .Using Young's inequality for some α > 0 and proceeding as in Eq. ( 24), we obtainT 2 ≤ η 2 k α 2bn m i=1 ∥A ⊤ k I b(i-1)↑ y k ∥ 2 2 + η 2 k 2bnα m i=1 ∥A ⊤ k I (i) (v k -y k )∥ 2 2 ≤ η 2 k nα 2b 2 Ĝπ (k) ∥y k ∥ 2 Γ -1 k + η 2 k 2nα Gπ (k) ∥v k -y k ∥ 2 Γ -1 k ,where we use our definitions that Ĝπ (k) :=1 mn Γ 1/2 k m j=1 I b(j-1)↑ A k A ⊤ k I b(j-1)↑ Γ 1/2 k 2 and Gπ (k) := 1 b Γ 1/2 k m j=1 I (j) A k A ⊤ k I (j) Γ 1/2 k 2 . It remains to choose α = 2b n Gk Ĝkto finish the proof.We are now ready to prove Theorem 3 for the convergence of shuffled SGD in the convex nonsmooth Lipschitz settings. E[H K (f ( xK ) -f (x * ))] ≤ 1 2n ∥x 0 -x * ∥ 2 2 + K k=1 2η 2 k n Ḡ,As a result, for any ϵ > 0, there exists a step sizeη k = η such that E[f ( xK ) -f (x * )] ≤ ϵ after O n Ḡ∥x0-x * ∥ 2 2 ϵ 2individual gradient queries.Proof. To simplify the presentation of our analysis, we first assume ∥v∥ 2 Γ -1 ≤ n, which will be later verified by our choice of v = y xK and Assumption 5.Combining the bounds in Lemma 10 and 15 and plugging them into Eq. ( 40), we haveEQUATIONwhere we use Young's inequality for ∥v k -y k ∥ 2Γ -1 k and ∥v k ∥ Γ -1 k = ∥v∥ 2 Γ -1as v is a fixed vector for (i), and (ii) is due to ∥y k ∥ 2 Γ -1 k ≤ n by Assumption 5 and assuming that ∥v∥ 2 Γ -1 ≤ n. Proceeding as the proof for Theorem 2, we first assume the RR scheme and take conditional expectation w.r.t. the randomness up to but not including k-th epoch, then we obtainE k [E k ] ≤ 2η 2 k nE k Ĝπ (k) Gπ (k) b .Since the randomness only comes from the random permutation π (k) , we haveE k [E k ] ≤ 2η 2 k nE π [ Ĝπ Gπ ] b .For notational convenience, we denote Ḡ = E π [ Ĝπ Gπ ], and further take expectation w.r.t. all the randomness on both sides and use the law of total expectation to obtainEQUATIONFor the SO scheme, there is one random permutation π generated at the very beginning such that π (k) = π for all k ∈ [K]. So we can directly take expectation w.r.t. all the randomness on both sides of Eq. ( 45), with the randomness only from π, which leads to the same bound as Eq. ( 46)with E Ĝπ (k) Gπ (k) = E π Ĝπ Gπ . Note that for incremental gradient (IG) descent, we can let Ḡ = Ĝ0 G0 without randomness involved, where Ĝ0 = Ĝπ (0) and G0 = Gπ (0) w.r.t. the initial, fixed permutation π (0) of the data matrix A.Noticing thatE k = η k Gap v (x k , y * ) + b 2n ∥x * -x k ∥ 2 2 -b 2n ∥x * -x k-1 ∥ 2 2 and telescoping from k = 1 to K, we have E K k=1 η k Gap v (x k , y * ) ≤ b 2n ∥x * -x 0 ∥ 2 2 - b 2n E[∥x * -x K ∥ 2 2 ] + K k=1 2η 2 k n Ḡ b . Noticing that L(x, v) is convex wrt x, we have Gap v ( xK , y * ) ≤ K k=1 η k Gap v (x k , y * )/H K , where xK = K k=1 η k x k /H K and H K = K k=1 η k , so we obtain E H K Gap v ( xK , y * ) ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 2η 2 k n Ḡ b .Further choosing v = y xK , which also verifies ∥v∥ 2Γ -1 = ∥y xK ∥ 2 Γ -1 ≤ n by Assumption 5, we obtain E[H K (f ( xK ) -f (x * ))] ≤ b 2n ∥x 0 -x * ∥ 2 2 + K k=1 2η 2 k n Ḡ b .To analyze the individual gradient oracle complexity, we choose constant stepsize η. Then, the above bound becomesE[f ( xK ) -f (x * )] ≤ b 2nηK ∥x 0 -x * ∥ 2 2 + 2nη Ḡ b . Choosing η = b∥x0-x * ∥2 2n √ K Ḡ , we have E[f ( xK ) -f (x * )] ≤ 2 √ Ḡ∥x 0 -x * ∥ 2 √ K .Hence, given ϵ > 0, to ensure E[f ( xK ) -f (x * )] ≤ ϵ, the total number of individual gradient evaluations will benK ≥ 4n Ḡ∥x 0 -x * ∥ 2 2 ϵ 2 ,thus completing the proof.We now briefly discuss this result. The total number of individual gradient queries isO n Ḡ∥x0-x * ∥ 2 2 ϵ 2, which appears independent of the batch size, but this is actually not the case, as the parameter Ḡ = E π [ Ĝπ Gπ ] depends on the block partitioning, due to Eq. ( 4). When b = n, as a sanity check, we recover the standard guarantee of (full) subgradient descent, which is expected, as in this case shuffled SGD reduces to subgradient descent. When b = 1, however, the bound is worse than the corresponding bound for standard SGD, by a factor O(n Ḡ/G 2 ). By a similar sequence of inequalities as in Eq. ( 4), this factor is never worse than n, but it is typically much smaller, taking values as small as 1. We note that it is not known whether a better bound is possible for shuffled SGD in this setting, as the only seemingly tighter upper bound from [42] applies only for constant K, when n = Ω( 1 ϵ 2 ), and under an additional boundedness assumption for the algorithm iterates.We implement the computation of L and L max in Julia, a high-performance scientific computation programming language, and compute matrix operator norms using the default settings in the Julia Arpack Package. However, limited by computational memory and time constraint, our selection of datasets is focused on moderately large-scale datasets of n in the order of O(10 5 ). We also include comparisons of small datasets such as a1a and sonar.We first study the gap between Lπ and L max for different batch sizes b, as shown in Figure 2 . As in Section 4.1, we focus on their dependence on the data matrix, and assume that the loss functions ℓ i all have the same smoothness constant. In this case, the ratio L max / Lπ that characterizes the gap between Lπ and L max will become L max / Lπ = (max1≤i≤n {∥a i ∥ 2 2 })/ 1 b ∥ m j=1 I (j) A π A ⊤ π I (j) ∥ 2 .In particular, we run experiments on standard Gaussian data of size (n, d). We fix the dimension d = 500, and vary the number of samples with n = 100, 500, 1000, 2000. In Figure 2 , we plot the ratio L max / Lπ versus the batch size b for 100 different random permutations π, where the dotted lines represent the mean values and the filled regions indicate the standard deviation of permutations. We observe that the ratio L max / Lπ is concentrated around its empirical mean and exhibits b α (α ∈ [0.74, 0.87]) growth as the batch size b increases. In particular, if we choose b = √ n, the ratio can be O(n 0.4 ). In this subsection, we include histograms in Figure 3 to illustrate the spread of L max / Lπ with respect to random permutations, for completeness. We observe that in all the examples L max / Lπ is concentrated around its empirical mean. The following plots are normalized, with y-axis representing the empirical probability density. The x-axis represents L max / Lπ . Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: N/A Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We use public benchmarking datasets from LIBSVM [15] , MNIST [17] , CIFAR10 [22] , and Broad Bioimage Benchmark Collection [28] , and have properly cited and credited the asset's creators.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: N/A Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: N/A Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification: N/A Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.This is because it is possible for inequalities (i) and (ii) to be loose, in addition to (iii).This is without loss of generality, as it can be achieved by reordering the rows in the data matrix. 
 
 Appendix:  


Decision: The paper studies shuffling-based SGD on convex finite-sum optimization problems. Through a primal-dual analysis on the primal algorithm SGD, the authors develop new convergence theorems for random-reshuffling and shuffle-once SGD with alternative smoothness constants that can improve upon the existing bounds that contain the max smoothness constant. The authors then develop convergence bounds for smooth and non-smooth convex generalized linear model cases, demonstrating similar but more fine-grained bounds that can better reflect the data matrix.

The paper presents nontrivial and novel analysis inspired from the connection of finite-sum SGD with coordinate descent in the dual space. The bounds improve upon existing results in the literature and can offer tighter bounds especially when there is a large gap between the max smoothness constant and the "average" smoothness. The discussion provided in Section 4 is inspiring, in that it suggests that we should take the structure of data into account in the convergence analysis to get tighter results.

The newly proposed constants $\tilde L$, $\hat L$, etc…  are not readily interpretable, as the authors also admitted in the discussion. The authors provide some numerical experiments demonstrating that they can be smaller than $L_\max$, but I recommend the authors to supplement the paper with some theoretical analysis of the proposed constants in some specialized scenarios, e.g., when $L_i$'s decay as $i^{-\alpha}$. 

Overall, the reviews were mostly positive and the contributions look solid; hence I recommend acceptance of this paper.

Reviews: review: paper_summary: paper_summary: This paper focuses on SGD with shuffling for finite-sum minimization problems. While there exists tight convergence upper and lower bounds for SGD assuming a global smoothness constant $L_{\max}$, the proposed results give a more fine-grained analysis in terms of the component-wise smoothness constants $L_i$ to show possibly faster convergence rates using Fenchel dual arguments. By replacing the $L_{\max}$’s with smaller constants $\tilde{L}, \hat{L}$, the new results can improve the iteration complexity by up to a factor of $L_{\max}/\hat{L} = O(n)$ as demonstrated in the linear predictor example. This work also extends to non-smooth, Lipschitz objectives in the case of linear predictors, and provides empirical evidence supporting the claims.
strength_weakness: strength_weakness: - Focusing on the results, the only big difference with previous work would be the constant $L_{\max}$ being replaced with $\tilde{L}, \hat{L}$, which means that the improvements are directly linked with how small the ratios $\tilde{L}/L_{\max}, \hat{L}/L_{\max}$ could actually be. While Table 2 suggests that $L_{\max}$ could be quite large (for the linear predictor case), it is quite hard to see how ‘small’ the values $\tilde{L}, \hat{L}$ are in first glance due to the relatively complicated definitions. The paper explains the gap for Gaussian data, but I think this example is a bit too ‘specific’.
- Moreover, based on my understanding, the definitions of $\tilde{L}, \hat{L}$ for the general convex case and the linear predictor case seem to be a bit different. It is still unclear whether these fine-grained bounds could imply a significantly *faster* convergence rate, especially for the general convex case.
- Minor Stuff: I might have missed this, but it seems that the definition of $\boldsymbol{A}_{\pi}$ (which I think is $\\boldsymbol{A}\_{\\pi}  = [a\_{\pi_1} \cdots a\_{\pi_n}]^{\top}$) is missing in Section 3.1. Also, the proofs are all written for batch size $b$, might be confusing to readers: maybe it would have been better if the main body also contained statements including $b$.
- TYPO: Line 9 of Algorithm 1, $m$ → $n$
- TYPO: Line 473 mini-bath → mini-batch
questions: questions: - Is there a way to intuitively understand when or why $\tilde{L}/L_{\max}, \hat{L}/L_{\max}$ could be ‘much’ smaller than $1$ (with possible $n^{\alpha}$-ish dependencies), both for the  Or are there at least any examples apart from the linear predictor that demonstrate significantly small $\tilde{L}/L_{\max}, \hat{L}/L_{\max}$?
- Based on my knowledge, in coordinate descent algorithms it is common to use *different* step sizes for different coordinates, say $\eta_i = 1/L_i$, to improve dependencies on the $L$’s. Would it be possible to do something similar for this case if we take the primal-dual approach?
- Can you elaborate a bit more on why the inequality $(iii)$ of $(4)$ (in Section 4) is loose by a factor of $n$ in most cases? My understanding of the linear predictor case is that the norm-trace gap corresponds to the largest VS sum of singular values, and hence if the singular values are even, then the gap can be as loose as a factor of $n$, while if we have a single large $\lambda_{\max}$ then the inequality can be tight. I wonder if (i) this interpretation is correct, and (ii) whether something similar could hold for the $\tilde{L}^g, \hat{L}^g$’s in the general convex case.
- Minor question: Is it true Algorithm 1 is completely equivalent to the standard SGD with random reshuffling (without the $\boldsymbol{y}$’s), and the dual variables are just there for illustrative purposes?
limitations: limitations: The authors adequately addressed the limitations, and there seems to be no potential negative societal impact.
paper_summary: paper_summary: This paper considers the primal-dual aspect of SGD using sampling without replacement (shuffled SGD). The authors present results for shuffled SGD on smooth and non-smooth convex settings with tight convergence bounds for several shuffling schemes (IG, SO, and RR). In some specific settings, the convergence rate can be improved with a factor of $\sqrt{n}$. The authors perform experiments to demonstrate that their bound is tight and better than prior work.
strength_weakness: strength_weakness: The improvement in the convergence is not significant as it is only changing the Lipschitz constant. 
The authors should be more transparent about the setting where there is $\sqrt{n}$ improvement as it seems like this only applies to linear predictors (the claim in abstract, line 14-17 is not clear). The authors should include that context in your statements. Another weakness is that the bounds in Section 3 depend on the data.
questions: questions: Why do you call RR and SO uniformly random shuffling? 

How the results in section 3 reduces and compares to general convex finite-sum problems when $a_i$ are all-1 vectors i.e. loss function applies on $x$?

I have read the author responses.
limitations: limitations: NA
paper_summary: paper_summary: This paper studied the convergence of random shuffled SGD through the lens of dual coordinate descent. By leveraging the analysis in coordinate descent, the author(s) derived a rate that is $O(\sqrt{n})$ faster than existing rate.
strength_weakness: strength_weakness: - The primal-dual relationship between SGD and coordinate descent has a quite long history, e.g. [1]. The author(s) should mention these works in the related work section or section 2.

[1] Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization, Shai Shalev-Shwartz and Tong Zhang, 2013.
questions: questions: The convergence analysis of coordinate descent usually have to assume the objective is smooth, which means the primal objective is strongly convex. But the strong convexity assumption was not used in this paper, I wonder how does the author(s) circumvent this assumption?
limitations: limitations: Limitations are not well-addressed. I do not see any potential negative impact of this paper.
paper_summary: paper_summary: This paper aims at improving existing bounds on random reshuffling. While lower bounds are tight in the worst case, refined smoothness definitions allow to take a larger step-size in favorable cases, which in turns allows faster convergence. These smoothness definitions involve a supremum over permutation-dependent partial averaging of individual smoothness constants. 
 
The dependence on these refined smoothness constants is obtained through a primal-dual interpretation of SGD, though the convergence analysis does not consist in applying existing primal-dual analyses to this formulation. As a side note, results include mini-batching, which is nice to have. 

For generalized linear models, individual covariances can be averaged instead of individual smoothness constants, thus allowing to gain even more especially in high-dimensions. A non-smooth extension is given, performing similar averagings over individual Lipschitz constants. 

Typos:
Line 27: why does SGD rely on a standard basis vector?? This looks like a mix between SGD and coordinate descent. 
L45: we choose take 
Algo 1: shouldn't this be a $\nabla f_i$ on line 7?
strength_weakness: strength_weakness: - $\tilde{L}$, the $L_{max}$-like constant, is the one that appears in Theorems 1 and 2. In particular, the more "average smoothness" constant (from which most of the improvement seems to come) only affect higher-order terms by using larger step-sizes to be used, but does not change anything for small step-sizes. 
- While results are strong for generalized linear models (taking into account the actual spectrum of each Hessian), general results are quite looser (though they still improve on existing ones).
questions: questions: 1) In the end, everything seems to depend on primal quantities. What makes an actual difference in the analysis that could not have been done in the primal (or why was it more intuitive to do it in the primal)? Why can't the full analyses be done by replacing the $y^i_k with $\nabla f_i(x_{k,i})$?

In particular, what technical challenges do you expect moving to the non-convex setting? In the end all quantities can, but the primal-dual interpretation of SGD does not hold anymore since biconjugacy is used at some point if I'm not mistaken).

2) What are the connections with existing (primal-dual) CD analyses? I see that the primal-dual reformulation allows to study a setting that is close (although different since CD would not get any residual variance). It is said in Appendix A that results are technically disjoints, so then why is this inspiration helpful still? 

3) It seems that the discussion from Section 4 only applies to the generalized linear models case, because it does not seem possible to benefit from structure in the general case. Do you have any intuition for why such assumptions cannot be made using Hessians, and thus benefit from similar improvements? Additionally, how would Table 2 look with these bounds instead of the generalized linear models ones?

4) Do you actually expect to depend on the smoothness for the average permutation (instead of worst-case)? Is this complicated to analyze directly because of the dependence between iterates within a single epoch?
limitations: limitations: See weaknesses

Comments: title: Thank you
Thanks for your response
A small request for clarification
References for the rebuttal
References for the rebuttal
comment: We thank the reviewer for their feedback, indicating **no ethical concerns** in our paper. Regarding the summarized technical concerns from the initial reviews, we note that they have been addressed during the rebuttal and discussion periods, as acknowledged by the three reviewers who have responded. We further note that Reviewer EeyJ found their concerns satisfactorily addressed and increased their score accordingly.
We thank the reviewer for their feedback indicating **no ethical concerns**, and for pointing us to potential applications of our work in other domains, which we plan to explore.
Thank you for acknowledging our response and for the clarification regarding the Mahalanobis semi-norm extension.
Thank you for carefully considering our response and acknowledging that your questions have been addressed. We appreciate your increased support of our work.
I would like to thank the author(s) for the detailed response. My questions are addressed, I have raised my score slightly.
Thank you for carefully answering my questions. 

I believe the Mahalanobis semi-norm extension (which as you point out allows to bypass the second order differentiability requirement, though might be less tight for non-GLMs) would be nice to have, but I don't consider it compulsory.
Thank you for reading through our rebuttal. We are glad to learn you found our answers detailed and precise.

We noticed that you kept your score as 'borderline' and wanted to kindly ask if you could provide some additional insight into your reasoning. Understanding your perspective would be valuable for us to address any remaining concerns and improve the work further.

Thank you in advance,
Authors
Thank you for your detailed and precise answers.
### References

[1] Chambolle, A. and Pock, T. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120–145, 2011.

[2] Mishchenko, K., Khaled, A., and Richtarik, P. Random reshuffling: Simple analysis with vast improvements. In Proc. NeurIPS’20, 2020.

[3] Koloskova, A., Doikov, N., Stich, S. U., & Jaggi, M. On Convergence of Incremental Gradient for Non-Convex Smooth Functions. In Proc. ICML'24, 2024.

[4] Wang P, Zarifis N, Diakonikolas I, Diakonikolas J. Robustly learning a single neuron via sharpness. In Proc. ICML’23, 2023.

[5] Lee, C.-P. and Wright, S. J. Random permutations fix a worst case for cyclic coordinate descent. IMA Journal of Numerical Analysis, 39(3):1246–1275, 2019.

[6] Wright, S. and Lee, C.-P. Analyzing random permutations for cyclic coordinate descent. Mathematics of Computation, 89(325):2217–2248, 2020.

[7] Xufeng Cai, Chaobing Song, Stephen J Wright, and Jelena Diakonikolas. Cyclic block coordinate descent with variance reduction for composite nonconvex optimization. In Proc. ICML'23, 2023.

[8] Chaobing Song and Jelena Diakonikolas. Cyclic Coordinate Dual Averaging with Extrapolation for Generalized Variational Inequalities. SIAM Journal on Optimization, 2023.

[9] Cheuk Yin Lin, Chaobing Song, and Jelena Diakonikolas. Accelerated cyclic coordinate dual averaging with extrapolation for composite convex optimization. In Proc. ICML'23, 2023.
### References

[1] Jaeyoung Cha, Jaewook Lee, and Chulhee Yun. Tighter lower bounds for shuffling SGD: Random permutations and beyond. In Proc. ICML'23, 2023.

[2] Konstantin Mishchenko, Ahmed Khaled, and Peter Richtárik. Random reshuffling: Simple analysis with vast improvements. In Proc. NeurIPS’20, 2020.

[3] 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. The Journal of Machine Learning Research, 22(1):9397–9440, 2021.

[4] Xufeng Cai, Chaobing Song, Stephen J Wright, and Jelena Diakonikolas. Cyclic block coordinate descent with variance reduction for composite nonconvex optimization. In Proc. ICML'23, 2023.