Shuffling the Data, Extrapolating the Step: Sharper Bias In Constant Step-Size SGD
Authors: Konstantinos Emmanouilidis
Emmanouil-Vasileios Vlatakis-Gkaragkounis
Rene Vidal
Venue: ICLR 2026 Poster (2026)
Keywords: variational inequalities
multi-agent optimization
stochastic algorithms
Primary area: optimization
TL;DR: Cubic Acceleration in Finite-Sum VIPs

Abstract: From adversarial robustness to multi-agent learning, many machine learning tasks can be cast as finite-sum min–max optimization  or, more generally, as variational inequality problems (VIPs). Owing to their simplicity and scalability, stochastic gradient methods with constant step size are widely used, despite the fact that they converge only up to a constant term. Among the many heuristics adopted in practice, two classical techniques have recently attracted attention to mitigate this issue: *Random Reshuffling* of data and *Richardson–Romberg extrapolation* across iterates. Random Reshuffling sharpens the mean-squared error (MSE) of the estimated solution, while Richardson-Romberg extrapolation acts orthogonally, providing a second order reduction in its bias. In this work, we show that their composition is strictly better than both, not only maintaining the enhanced MSE guarantees but also yielding an even greater cubic refinement in the bias.
To the best of our knowledge, our work provides the first theoretical guarantees for such a synergy in structured non-monotone VIPs. Our analysis proceeds in two steps: (i)
we smooth the discrete noise induced by reshuffling
and leverage tools from continuous-state Markov chain theory to establish a novel law of large numbers and a central limit theorem for its iterates; and (ii) we employ spectral tensor techniques to prove that extrapolation 
debiases and sharpens the asymptotic behavior even under the biased gradient oracle induced by reshuffling. Finally, extensive experiments validate our theory, consistently demonstrating substantial speedups in practice.

Published as a conference paper at ICLR 2026
SHUFFLING THE DATA, STRETCHING THE STEP-SIZE: SHARPER BIAS IN CONSTANT STEP-SIZE SGD
Konstantinos Emmanouilidis' Emmmmmmiilllssssmmmmaiilaaaaaaa
Vlatakis-Gkaragkkuniis René Vidal"
Department of Computer Science
"University of Pennsylvania, # Univeesityy of Wisconsin-Madison
ABSTRACT
From adversarial robustness to multi-agent learning, many machine learning tasks can be cast as finite-sum min-max optimization or, more generally, as variational inequality problems (VIPs). Owing to their simplicity and scalability, stochastic gradient methods with constant step size are widely used, despite the fact that they converge only up to constant term. Among the many heuristics adopted in practice, two classical techniques have recently attracted attention to mitigate this issue: Random Reshuffling of data and SettaardoooRmmmbe extrapolation across iterates. Random Reshuffling sharpens the mean-squared error (MSE) of the esti- mated solution, while Bechardson-Rombeee
extrapolation acts orthogonally, pro- viding a second order reduction in its bias. In this work, we show that their com- position is strictly better than both, not only maintaining the enhanced MSE guar- antees but also yielding an even greater cubic refinement in the bias. To the best of our knowledge, our work provides the first theoretical guarantees for such a synergy in structured non-monotone VIPs. Our analysis proceeds in two steps: (i) we smooth the discrete noise induced by reshuffling and leverage tools from continuous-state Markov chain theory to establish a novel law of large numbers and a central limit theorem for its iterates; and (ii) we employ spectral tensor tech- niques prove that extrapolation debiases and sharpens the asymptotic behavior even under the biased gradient oracle induced by reshuffling. Finally, extensive experiments validate our theory, consistently demonstrating substantial speedups in practice.
INTRODUCTION
Variational inequality problems (VIPs) (Stampacchia, 1964) are a unifying framework that extends beyond classical loss minimization to encompass min-max optimization, complementarity prob- lems (Dantzig & Cottle, 1968; Facchinei & Pang, 2003), equilibrium computation in games, and general fixed-point formulations (Bauschke & Combettes, 2017). In recent years, VIPs have gained significant traction in ML and data science, especially due to their broad potential applicability in domains where minimizing single empirical loss is insufficient. Notable examples include gener- ative adversarial networks (Goodfellow et al., 2014; Arjovsky et al., 2017), multi-agent and robust reinforcement learning (Namkoong & Duchi, 2016; Wang al., 2021; Giannou et al., 2022), and auction theory (Syrgkanis et al., 2015).
In practice, many of these tasks reduce to finite-sum formulations, where the objective depends on large collection of data samples agents. In such settings, stochastic gradient methods have become the workhorse of lurge-scale learning (Bottou et al., 2018). By exploiting the finite-sum structure, stochastic gradient descent (SGD) and its variants replace expensive full-gradient computations with inexpensive updates on a few components, enabling scalability to massive datasets.
While the theoretical underpinnings of SGD date back to the seminal work of Robbins & Monro (1951) and Kiefer & Wolfowitz (1952), the initial results involved an asymptotic analysis under vanishing step size. Consequent works (Ljung, 1978; Benaïm, 2006; Rakhlin et al., 201 I; Raginsky et al., 2017; Azizian et al., 2024; Malick & Mertikopoulos, 2024) have extended the theoretical guarantees of SGD under different design choices, much of its practical success can be traced to a handful of seemingly "low=level" heuristics (Bottou, 20 2012b): step-size schedules (constant vs. decaying), data ordering (with without resampling), and iterate selection (average vs. last iterate).

---

Published as a conference paper at ICLR 2026
To facilitate analysis, the community has typically adopted a ceteris paribus perspective-isoleeee one design choice at a time while holding the rest fixed an approach that clarifies individual effects but obscures their interaction.
A particularly important case is the use a constant step size ", which is popular in practice since it simplifies tuning, quickly erases the dependence initialization, and yields fast early progress (Yu et al., 2021). However, its fundamental drawback is that convergence halts at a non- vanishing error. Even in strongly convex problems with unique solution x², the last iterate of SGD typically satisfies:
MSE(SGD) = lim sup E[| - x= = O()) and bias(SGD) = lim sup [| - x"1| = O()).
k-
k-000
That is, the iterates stabilize a distance from the optinum on the order of the step size.
To mitigate this limitation, practitioners often turn to debiasing heuristics. A prominent example is random reshuffling (RR1), or without-replacemmnt sampling, where each data point is visited exactly once per epoch. Unlike classical with-rrelaccment SGD, which may resample or skip points, RR enforces a random full pass that closely mirrors large-scale training in practice (Bottou et al., 2018). Despite the dependence it induces across samples, recent work has established faster convergence guarantees for RR, in both minimization (Ahn et al., 2020; Gürbüzbalaban et al., 2021; Cai et al., 2023) and VIPs (Mishchenko et al., 2020b; Emmanouilidis et al., 2024), along with sharper MSE bounds from (nx) to (()²). However, the question of whether the bias term improves has remained open, Indeed, recall that for any estimator the mean squared error decomposes as
MSE
E[I|a x= [l²]
||
E[a] x= + Var(i),
so that (ias($) < / MS) Under this trivial bound, SGD-RRJ guarantees improved MSE compared to vanilla SGD, but does not necessarily yield smaller bias.
Orthogonal to reshuffting, another classical idea from numerical analysis has recently re-emerged CittardsonRRobbrr (RR2) extrapolation. Its principle is simple yet powerful: Run the algorithm of your choice with two different step sizes and combine their outputs so that the leading bias term cancels. Concretely, whenever the bias admits an expansion of the form bias(7) = A7 + O(7a)) with & > 1, running the stochastic approximation at two step sizes gives:
in stochastic optimization:
x²6o - Ay + 0(7")) and 2 - x 2A + 0(7").
Extrapolating these iterates then yields:
Textr x= 2r% 2 - 2 2A9 +
Originally introduced for accelerating discretization schemes in stochastic differential equations (Hildebrand, 1987; Talay & Tubaro, 1990; Bally & Talay, 1996), RR2 has since been applied to optimization, improving constant-step methods from SGD (Durmus et al., 2016; Dieuleveut al., 2020: Mangold et al., 2024; Sheshukova et al., 2024) to Q-learning and two-timescale stochastic approximation methods (Huo et al., 2023; Kwon et al., 2024; Zhang & Xie, 2024; Allmeier & Gast, 2024). However, despite its conceptual simplicity and empirical success, the theoretical foundations of RR2 for stochastic VIPs remain nascent (Vlatakis- Saaakkii-GGraakoonnii et al., 2024).
In summary, the known bias rates of the heuristics remain limited when they are applied in isolation: For unconstrained strongly monotone VIPs, RR, is known to sharpen MSE bounds (from O()) to O()²)) but does not, in general, guarantee improved bias order; on the other hand RR2 alone attains 3/23/2) bias VlattakisGKaarakkounii et al., 2024)!. However, it remains unclear whether one can combine effectively the aforementioned heuristics both, ideally reaching ()²²) bias. This raises the natural question:
unified method that outperforms
What new phenomena arise when these heuristics
constant step sizes, random reshuffling, and Richardson extrapolation interact simultaneously?
Addressing this question is delicate, because reshufiling introduces biased stochastic oracle whose discrete, permutation-driven noise structure is not covered by existing analyses of extrapolation,
'The (()2) rate in Vlatakis- Vakakis-Gaaragkounis et al. (2024, Sec. 5, Thm. 6) is obtained via reduction to Dieuleveut et al. (2020, Sec. 3, Thm, 4), which requires additional noise assumptions met in our setting.
2

---

Published as a conference paper at ICLR 2026
which predominantly assume unbiased or continuously distributed perturbations (Dieuleveut et al., 2020; Sheshukova et al., 2024; Vlataki-GGaaagkkone et al., 2024).
Our contributions. In this work, unify the RR1 and RBg heuristics into principled algorithmic framework, showing that their composition yiclds a level of bias reduction unattainable by either heuristic alone. Our main result (Theorem 4.6) shows that for quasi-strongly monotone smooth VIPs, our combined method (SGD-RRy@RR, Algorithm 1) cancels all lower-order terms in the bias expansion, yielding an asymptotic bias of order ((7²).
Our contributions are as follows:
Convergence of SGD-RR;, under perturbations. We first analyze the RR, component, proving ponential convergence with a bias that is linear (quadratic) under weak (quasi-strong) monotonicity, robust even to preprocessing perturbations.
Epoch-level Markov chain viewpoint of SGD-RR). We prove that the dynamics of RRi admit an invariant stationary distribution and establish a lar-o--large-numbers (LLN) and central-limit- theorem (CLT) for the per-epoch iterates of SGD-RRJ.
Cubic Bias Refinement through the synergy of RRy We extrapolate the iterates produced from SGD-RR) and analyze the combination of the heuristics. Through a spectral analysis of the reshuffling-nndus induced biased oracle, we prove that the first ((73) bias guarantee for our algorithm in quasi strong monotone VIPs.
PROBLEM SETUP AND ASSUMPTIONS
Variational inequalities. Let's recall first the basic framework of finite-sum variational inequal- ities (VIs), which will underlie our analysis. Let X S Rd be nonempty closed convex set and F Rd Rd single-valued operator. The variational inequality problem VI(X, F) asks for point x² E X such that
(F(x"), x-)))
E
(VIP)
In our setting, we focus on the unconstrained finite-sum case, i.c. X x and F(x) I E. (x(,), respectively, where each F, : IRd + IRd typically represents the gradient of loss function for data point in some dataset D. To build intuition, we illustrate the framework through a few canonical examples below:
Example 2.1: Solving Non-linear equations. A solution x² to the (VIP) corresponds to a root of the equation F(x) = 0. Therefore, we can cast any non-linear equation as a specific instantiation of the Variational Inequality framework. Well-known examples include the Navier-Stoker equations in computational dynamics (Hao, 2021).
Example 2.2: Empirical Risk Minimization. For any -moohh loss function É RRd
solution x² to the (VIP) with F(x) = Vl(x) is a critical point (KKT solution) to the associated empirical risk minimization problem, consisting the comerstone of machine learning objectives. Example 2.3: Nash Equilibria & Saddle-point Problems. Consider N players, each having an action set in Rdd and convex cost function c; Rd R. A Nash Equilibrium (NE) is a joint-action profile x² (xr)) that satisfies
cxxx") co Vi, E RR.
(NE)
For convex cost functions c; Rd + R. a (NE) coincides with the solution a (VIP) with operator = =
(V, In the particular case of two players and (quasi) convex-concave objective
L: Rd x Rd + IR, the solution I = x5) to the (VIP) with F(x) (V (xx), -VC() is saddle point of satisfying
C(ri, x2) & C(zi x3)) < C(x1, Væ 22 E Rdd
Saddle-point problems and applications of (NE) are ubiquitous, from training Generative Ad- versarial Networks (GANs) to multi-agent reinforcement learning and auction/bandit problems (Daskalakis al., 2017; Zhang et al., 2021; Pfau Vinyals, 2016).
Blanket assumptions. We now state the standing assumptions for our analysis, beginning with the existence of a solution x" to (VIP).
3

---

Published as a conference paper at ICLR 2026
Assumption 2.1. The solution set X of (VIP) is nonempty and there exists x' E X" R E R such that || 222 S R.
The next assumption introduces the class of operators F of the associated (VIP) for which our stochastic gradient algorithms will be analyzed.
Assumption 2.2 (A-weak je-uasii strong monotonicity). The operator F is A-weak µ-quast strongly monotone, i.c. there exist 2 0, μ > 0 such that for some x² E 2 it holds that
(F(x), x x²) 2 u=| x x²1| X,
Va E Rd.
(1)
Assumption 2.2 for À coincides with the well-known
notions of quasi-strong monotonicity (Loizou et al., 2020), strong stability condition (Mertikopoulos & Zhou, 2019), and strongly coherent VIPs (Song et al., 2020) in the opti- mization literature. can be seen as a relaxation of the clas- sical notion of strong monotonicity/convexiii
which re- Vx, x' E Rdd. For À > Assumption 2.2 represents further relaxation, motivated by dissipative dynamical systems and weakly convex optimization (Raginsky al., 2017; Erdogdu al., 2018), and it encompasses non-monotone games as well as variety of problems in statistical learning theory (Tan Vershynin, 2023).
quires (((z)=F(z*)) x x-x) 2 -=|==~21|2 x'|1
simple example of function satisfying Assumption 2.2 is f v) (1² + 7 sin(()) + xy cos(y))) where the assumption holds with (p. X) = (1,25).
Figure
common assumption in the literature of smooth optimiza- tion that we will utilize is that the operators in the finite-sum structure of the (VIP) are Lipschitz continuous.
Assumption 2.3 (Lipschitz continuity). Each F; is Lipschii:
[[((()) - Fi < L |2] 22l1,
Ven, 22 E RR. iE [m].
with Lmax maXie[1] L..
Unlike standard analyses assuming unbiased oracles with bounded variance, e.g.. (Loizou et al., 2021; Hsieh al., 2019; Lin et al., 2020; Mishchenko al., 2020b), random reshuffling induces bias via inter-step dependence. Such conditions may fail even for simple quadraties. Instead, we work directly with Lipschitz continuity and impose only mild moment bound:
Assumption 2.4 (Bounded moments the solution). At some x= E X", the oracle values have finite second and fourth moments:
[1(2")) xo,
oc.
Assumption 2.4 is mild: it does not require global boundedness of gradients. but only that the oracle values (x admit finite second and fourth moments at the solution. Building on this, we extend the variance bound of Emmanoullidis et al. (2024, Prop. A.2, 16) to higher-order moments.
Proposition 2.5. Let Assumptions 2.1 hold. Then, for any Rd it holds that
a [R(x)) F(x) <
n
Σ
L x
x11 112 + 20
n
n
1
128
(ii) [F(x)) - F(x)11 <
L | x²1 + 128o
==1
Having introduced the necessary assumptions, we, next, provide the proposed algorithm that will achieve the higher order bias refinement in results.
4

---

Published as a conference paper at ICLR 2026
3 THE ALGORITHMIC FRAMEWORK
While there are many conceivable ways to interleave RR₂ and RR,, both intra- and inter-epoch, we adopt the most natural and practically motivated design. In modern pipelines, RR1 is the workhorse the low-level training stage, while RR₂ often employed as a black-box refinement higher level, allowing parallelization and modular integration.
Accordingly, we study stochastic gradient algorithms that sample via random reshuffling to generate stochastic oracles of gradients/operators. At the start of each epoch k >> 0, random permutation wis [n] is drawn, prescribing the order in which data points are processed. The algorithm then performs the classical SGD update:
PreP (ieerroossssschiOrreeeede.i
SSCD-RR-033,(mmer--op)
where StochOracle (xf denotes either the stochastic gradient (in minimization problems) or the operator value Fug. (x%) (in the general VI case), indexed by the w data point and PreProcess{] preprocessing routine implementing calibrated Gaussian smoothing the input. Then, the final iterate of each epoch becomes the starting point of the next, and the procedure repeats.
Algorithm RR2RR
Require: Initial point E Re:: step size y > 0: epochs 1; dataset size
SToCHORACLE(2: i) returns (((2) (minimization) or operator value (VI):
PREPROCESS(y: i) adds calibrated Gaussian smoothing on g (e.g.. U2 ~ N'0 ()) do
1: for k= 0,
epoch
2: for = 27 do
Parallel iterations with two step-sizes
4: Draw random permutation 4i of [n]
3:
1...., do
b inner loop (reshuffled pass)
5:
1+1 k. +- Riin R ax (x))
6:
end
7:
11)1 <- =
0
baseline next-start (used for analysis)
8:
end for
9:
Sx+1 2 2x 2F,[20]
> outer loop extrapolation epoch end)
OR
ix+t < (2 Emejy xb/ =m/,)//)/
Alternatiee: ( (extrapolation at epoch's 's averages)
11: end for
return Rt
D (optionally average ( across epochs)
On the necessity of smoothing. A key challenge with reshuffling is that, after one epoch, the cu- mulative gradient estimator is biased, unlike sampling with replacement, which is unbiased and analytically simpler. The induced noise is also discrete, tied to permutations. To handle this, we introduce a calibrated Gaussian perturbation that smooths the discrete reshuffling noise into well- behaved proxy while preserving variance, moments, and bias order. In practice, the perturbation has negligible effect across datasets; clarifying its precise dependence on dataset size is an interesting direction for future work. For completeness, the supplement also includes brief sketch showing how our results extend even without this step.
Finally, at the end of each epoch we apply RR2, yielding the extrapolated update:
2 x
N
(SGD- ( (0
In Section 4, we prove that this combination achieves provable O(>²) bias- to the best of our knowledge, the first such result.
4 THEORETICAL GUARANTEES
We present the theoretical guarantees for the synergy of both heuristics in our algorithm, establishing cubic refinement in the bias attained, We start by analyzing the RR, component (inner loop of
5

---

Published as a conference paper at ICLR 2026
Aspect
Emmanouilidis et al. (2024)
Our work
Baseline Algorithm
SEG
SGDA
Model Assumptions (Smoothness)
P--L, Lipschitz
P.-L, Lipschitz
Model Assumptions (Drift)
F je-Strongly monotone
F stroilyyrongly monotone RR, @ RR;
Main hcuristic
RRI only
Asymptotic Bias order
O( + O
O(y")
Asymptotic MSE order
O(7*)
Condition number
Worse than vanilla-SEG EG-structure RR
Same as vanilla-SGDA
Mechanism
Bias cancellation (RR) RRa)
Table I: Summary of key differences between Emmanouilidis et al. (2024) and our results.
Algorithm 1), which serves the basis for performing RR2 (outer loop of Algorithm 1) on the produced iterates and describe the interplay of both heuristics the bias of Algorithm
4.1
CONVERGENCE GUARANTEES FOR RR1
Our first result concerns the convergence of the RR, component for A-weak μe-quasi strongly mono- tone VIPs. More specifically, in order effectively utilize Markov Chain tools study the RR, variant under preprocessing perturbations (inner loop of Algorithm 1), attributing for simplicity of exposition of our results the name Perturbed SGD-RR1.
Theorem 4.1. Let Assumptions 2.1-2.3 hold. If ? S Ymar: then the iterates of Perturbed SGD-RR, satisfy
[11=++1- 1.1 $
7m (1-
k+1
12 x²112
8 Dima 8A
L?
o²
+
2
112
/
where o² 1 EIE [[F(2*)] and Ymaz
min
3nTmmux
146 1,22
12n mmaz
Remark 1. Theorem 4.1 establishes exponential convergence up to a bias of O(?20? + A)) where the SA term inherent (Yu et al., 2021). For fair comparison we focus on the quasi-strongly mono- tone case (X 0), which already generalizes strong convexity. Our rate recovers known results for strongly monotone operators (Das al., 2022; Emmanouilidis et al., 2024) and extendy them to weak monotonicity. In this regime, reshuffling attains a much smaller bias than the O(7r72) of with-replacement SGD (Loizou et al., 2020; Gower et al., 2019), converging
tighter neighbor- hood. This sharper bias also yields faster accuracy rates: with 7 1/(nK) with-replacement SGD reaches (1)(nnK) accuracy (Das et al., 2022; Mishchenko et al., 2020a), while reshufflin accelerates to ((1/nn²²)), further support for its empirical success.
In the sequel, we view the algorithmic trajectory through the prism of Markov chain theory. This perspective enables a finer dissection of the reshuffling bias and, mutatis mutandis, equips us with the machinery to construct consistent estimators for performance statistics. The Markovian frame- work arises naturally, the method progresses from to I&+1 in a state-dependent fashion. The connection between stochastic approximation and Markov processes-trace such Robbins & Monro (1951); Pflug (1986)-ha has fueled a rich literature for algorithms with un- biased oracles. Random reshuffling, however, generates systematically biased oracles, necessitating genuine departure from this canonical line of analysis.
back to early works
For readers accustomed only to classical finite-state Markov chains, the transition mechanism is usu- ally represented by directed graph with fixed transition probabilities. In our setting, the analogue is the transition kernel P(x, A) = Pr Tnext E A Fnow x] A E B(Rd), where B(Rd) denotes the Borel sets of Rdd As in the finite-state case, it is highly desirable that this kernel remain invariant over time this is the property of home-homogeneity""
Empirically, for sufficiently large datasets the effect of discrete noise in smooth problems is negligible, making the preprocessing step unnecessary. detailed study of this effect lies beyond the scope of this pa- per, whose focus is instead the first systematic treatment of the interaction between Random reshuffling and Richardson extrapolation.
2If time-homogencic fails, process can be still Markovian in the sense that the future depends only the present, but its statistical regularity vary with time, complicating both analysis and long-run guarantees.

---

Published as a conference paper at ICLR 2026
At the step level, reshuffling destroys homogeneity: the transition kernel varies with the permutation index, making the process non-stationary. Fortunately, this irregularity vanishes at the epoch scale:
" "After one reshuffled pass, the law of the next iterate depends only on the epoch 's starting point and the drawn permutation, but not its position within the permutation."
Thus, the sequence of epoch-level iterates (æk J2z forms a bona fide time-homeen hime-homggeneous Markov chain, forming the basis for the asymptotic analysis of the RB₂ extrapolation component
(o)
Lemma 4.2 (Epoch-level homogeneity and kerne),, Fix y > and E N. Then the Perturbed SGD-RR, can be described at each epoch k as: Draw wi. uniformly from Gn and sel
Ik+1 = ( + Uk: Ux ~ N0,, (0,
where (x, w) denotes the endpoint of one reshuffled pass started at x with permutation w (i.e., the map induced by inner updates with step size 7). Then (2x)x20 time-homogeneous Markov chain on Rd with transition kernel
P(r,
1 H n! WEB,
E Av. (I,, 2) dy,
Σ
A B(R"),
where 6(;m,)) E) is the d-variate Gaussian density with mean m and covariance Σ.
By verifying irreducibility, aperiodicity, and positive Harris recurrence (Meyn & Tweedie, 2012), establish unique invariant distribution Ty, geometric convergence in total variation to it, and concentration of scalar observables (admissible test functions) around x="
Theorem 4.3. Under Assumptions 2.1-2.3, run Perturbed (GD-RR/with S (Ix k22 admits unique stationary distribution Thy E P2(R)), and additionally:
Tmaxx
Then
(i) [[((a=)] - ES [[(x)]] & p)*
Ve: |<(x)] - (1 + ||=]]),
(ii) EEwn, [[(x)] ((=))) I LVVC
Ve: L Lipschitz functions,
for some OC, (0,1), O((SSESSD RRI)) and Tmax defined in Theorem 4.1
geometrically in total variation to Ty Item (ii bounds the gap between the expectation of a measurement under n.y and its value at the solution x". Intuitively, if the method converged exactly to x", these expectations would coincide.
Remark 2. Item (i) of Theorem 4.3 shows that Perturbed SGD-RRR,covveree
The result of Theorem 4.3 follows from Foster-Lyapunov drift condition combined with mi- norization argument, showing that the induced Markov chain satisfies the standard ergodicity criteria in the spirit of Yu et al. (2021); Vlatakis- Vlatakis-Gkaragken Gkaragkounis al. (2024). Beyond geometric ergodicity, one may also ask whether the chain admits asymptotic statistical estimation functionals of its trajectory. By invoking the Birkhoff Khinnf-KKinccin ergodic theorem for continuous-stad Markov chains, we establish both a Law of Large Numbers (LLN) and a Central Limit Theorem (CLT) for empirical averages of test functions evaluated along the epoch iterates.
Theorem 4.4 (LLN and CLT for Perturbed SGD-RR1). Suppose Assumptions 2.1-2.3 hold and run Perturbed SGD-RR, witth 7 Tinax. (cf. Theorem 4.1). Let é : Rd - IR be any test function such that [(())) S (1 1 + [|x]|(2) and E EZ-T, [e(x)] < co. Then for the epoch-level iterates, holds that:
T-1
{(x+)
2y Exra, [((x)]
1/2
(xx) - Lont, ((x)) 1( 2 0(0,02, (0),
1=0
(LLN)
(CLT)
where 02, (c) = limy + Fr,, [S7] and S7 = T1 (f(x+) - ExNZ, 1(xx)))²
'On the augmented space Rd Gn. the chain ((z,, we 1420 is also time-homogeneous with kernel (((z,w), A B) = Sa Φ(v; H(x,w) E.a)) dy nf The above formulation is convenient for verifying Lyapunov-Foster and minorization criteria, since the coupling with uniform perturbation remains independent.
[B]

---

Published as a conference paper at ICLR 2026
4.2 CONVERGENCE GUARANTEES FOR RR
Having established the role of RR, within stochastic algorithms, we now examine its interplay with RR2 and the effect of combining these heuristics on the bias term. The previous results hold for the full class of weakly quasi-strongly monotone problems with À > O. To sharpen our understanding, Assumption 2.2). which already covers a broad range of non-monotone regimes (Loizou et al., 2020). A key step in our analysis is to bound higher-order moments of the deviation between RR; iterates and the solution of (VIP), thereby showing that the bias of Perturbed SGD-RR1 is linear in the step size with quadratic corrections.
we focus on the quasi-strongly monotone case (A
Technically, our analysis relies on two delicate ingredients that go beyond straightforward general- izations. (i) A spectral study of the full-puss operator (Lemma F.2), which approximates the under- lying map F of the VIP. This connection between RR, and the multi-step extragradient literature may be of independent interest, but its proof requires nontrivial handling of spectral properties across reshuffled passes. (ii) A combinatorial lemma (Lemma E.2) that bounds fourth order moments of finite-sum subsets of vectors. While reminiscent of Mishchenko et al. (2020b, Lemma 1, Sec. 7), our result demands substantially more intricate manipulations to accommodate the dependencies introduced by sampling without replacement.
Lemma 4.5. Let = and Assumptions 2.1-2.4 hold. If S Tmmx (cf. Lemma E.4), then SGD-RR.) lim sup [maa] x² C(x))(())))))))) 0
bias(Perturri
Remark 3. For classical SGD, the bias takes the form bias (SGD) = C(x")> + ((225)) (Dieuleveut al., 2020). Hence, while RB₁ retains the same first-order term, it improves the higher-order contribution and simultaneously yields sharper mean-squared error guarantees.
Building on this fact, we construct refined trajectory via the debiasing scheme RR2. Our final result shows that the combined scheme attains exponentially fast a provable asymptotic O(73) bias:
Theorem 4.6. Under the assumptions of Lemma 4.5, the output of Algorithm satisfies
Last-iterate version (line 9):
||
E [[xx]
- x"1l S c(1 - p)* +7()))
Averaged-iterate version (line 10):
E
1 K
Fm
x
S c/p k + (07))
where o E 10, c < oo (ef. Theorem 4.3).
Remark 4. Although the last-iterate estimator is often preferred in theory, in practice trade-off emerges vs ergodic-average: full-epoch or tailed averaging (the Polyak-Ruppert scheme (Polyak & Juditsky, 1992)) achieves improved variance properties, asymptotically captured by Theorem 4.4.
5 EXPERIMENTS
In this section, we conduct a series experiments demonstrating the benefits from the synergy of the two heuristics empirically. More specifically, focus in the strongly monotone setting and compare the relative error and bias attained by 4 variants: the classical SGD(A) algorithm using uniform with-replacement sampling (denoted as SGDA in the plots), the one equipped with RRi, the one equipped with RR2 and the method utilizing both heuristics. For each experiment, we report the average of 5 trials/runs and plot the relative error log the algorithm.
with respect to the iterations of |ize=""
2g 717
Two-player Zero-Sum Games. In the strongly monotone case, we consider the two-player zero-sum game from Emmanouilidis et al. (2024) and Loizou al. (2021), consisting strongly convex strongly concave quadratic of the form
min max f(1 x2) - A B,s2
! 2 2
1 C, + + af cT2
ZIER ER
For the interested reader, additional details on the experimental setup and the procedure used to sample the matrices A,, B,, Cs are provided in Appendix G.
8

---

Published as a conference paper at ICLR 2026
On the Rate of Convergence. In the first set of experiments, we aim to validate empirically the result of Theorem 4.1 by running SGDA with RR₁ and using the step sizes described by theory. We conduct experiments for multiple conditions
4 with value k {1, 10) and μ 1. In Figure 2, observe that the algorithm with RR, converges linearly to a neighbourhood around the solution x² and the neighbourhood depends on the step size used, validating in this way the results of Theorem 4.1. We have run experiments also for stepsizes that are larger than the ones predicted in theory, observing similar behaviour of the optimization algorithm. Additionally, we have performed an ablation study in Wasserstein GANs (Emmanouilidis et al., 2024; Daskalakis et al., 2017), showing that the performance benefit of the proposed heuristic is universal in many other common optimization algorithms used VIs (Appendix G).
Strongly Monotone Probtem [x=1]
Strangly Monotone Problem (x 5)
Strongly Monotone Problem (k= 101
10"
101
S686-A81
S2DARKI
1---
50D4-m2
5004-192
SGBA (B824881)
10-
30
52E 853-
10-
(X00)
Epochs (x100)
Epochs (x10))
Figure 2: Comparison of different heuristics. The combination of RR2@RR) converges linearly to neighborhood of the solution, validating the established theoretical results (Theorem 4.1). Even when using the last iterates, the combination of RR2ORR, converges to smaller relative error in comparison to the other variants (classical SGDA, RR1, RR2). This validates that the bias of Algorithm is improved even when R iterates are used.
Efficient Statistics & Empirical Concentration. This set of experiments examines the central limit theorem (CLT) and aims to validate empirically the theoretical results established in Theorem 4.4. The value of the game, which is zero, is used as the test value for which we observe the averaged evaluations after T (100, 500, 1000} iterations respectively. In particular, we run the algorithm with the step size suggested by Theorem 4.4 and maintain for the total number of iterations the sum of the evaluations, normalized with VT. We run the experiment for T = 2000 trials/runs and plot the corresponding histograms. In Figure 3, we observe that the histograms tend to concentrate the value of the game as the number of iterations increase. Additionally, we examine the effect of the step size concentration of the observed distributions.
Concentration Around the Game Value
Effect of Stepsize on CLT
Density
Density
6502
1835
2000
-30
-5
20
-35
53 15 200
4200
3000 -200
-300
530
200 300
Figure 3: Validation of concentration around the mean and effect of the number of iterations and step sized selected. The average of the game values tend to concentrate more around the mean for larger number of iterations and smaller step sizes. The right plot indicates the effect of two different step sizes 7 E {0.1, 0.001). showing that for smaller step sizes the corresponding distribution attains higher concentration around the mean of the values.
6 CONCLUSION
In summary, our work establishes that the synergy of random reshufiling and extrapolation yields principled reduction of bias, culminating in accelerated convergence guarantees for structured non- monotone VIPs. By combining Markov chain techniques, spectral amalysis, and higher-moment bounds, we provide the first rigorous evidence that these heuristics can be synergistically integrated rather than studied in isolation. This perspective bridges long-standing gap between practice and

---

Published as a conference paper at ICLR 2026
theory, offering a systematic framework that extends naturally broad class of constant step-size stochastic methods. We view this as a foundation for new generation of analyses where practi- cal heuristics not only empirically verified but also theoretically grounded to deliver provable performance improvements in complex stochastic optimization landscapes.
ACKNOWLEDGMENTS
Konstantinos Emmanouilidis and René Vidal acknowledge support from the NSF-Simons grant 2031985 "Research Collaborations the Mathematical and Scientific Foundations of Deep Learn- ing". Emmanouil- Vlatakis-Gkarakkounii acknowledges support from the startup funding Wisconsin Madison university under the grant (233-AAN5596 GF000020971).
REFERENCES
Rafik Aguech, Eric Moulines, and Pierre Priouret. On perturbation approach for the analysis of stochastic tracking algorithms. SIAM Journal on Control and Optimization, 39(3):872-89) 2000.
Kwangjun Ahn, Chulhee Yun, and Suvrit Sra. Sgd with shuffling: optimal without component convexity and large epoch requirements. In NeurlPS, 2020.
Sebastian Allmeier and Nicolas Gast. Computing the bias of constant-step stochastic approximation with markovian noise. Advances in Neural Information Processing Systems, 37:137873-137999 2024.
Kimon Antonakopoulos, Panayotis Mertikopoulos, Georgios Piliouras, and Xiao Wang. Adagrad avoids saddle points. In International Conference on Machine Learning, pp. 731-771. PMLR, 2022.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In ICML, 2017.
Waiss Azizian, Franck lutzeler, Jérôme Malick, and Panayotis Mertikopoulos. What is the long- run distribution of stochastic gradient descent? large deviations analysis. arXiv preprint arXiv:2406.09241. 2024.
Vlad Bally and Denis Talay. The law of the euler scheme for stochastic differential equations: I. convergence rate of the distribution function. Probability theory and related fiels, 104(1):43-60. 1996.
HH Bauschke and PL Combettes. Convex analysis and monotone operator theory in hilbert spaces, 2011 springer. New York, 2017.
Michel Benaim. A dynamical system approach to stochastic approximations. SIAM Journal Control and Optimization, 34(2):437447, 1996.
Michel Benaïm. Dynamics of stochastic approximation algorithms. In Seminaire de probabilites XXXIII, pp. 1-68. Springer, 2006.
Michel Benaïm and Morris Hirsch. Dynamics of morse-smale um processes. Ergodic Theory and Dynamical Systems, 15(6):1005-1) 1995.
Yoshua Bengio. Practical Recommendations for Gradient Based Training of Deep Architectures. Neural Networks: Tricks of the Trude, pp. 437-478, 2012. ISSN 1611-3344 1-349. doi: 10. 1001/
778-3-642-35288-8.22
Dimitri P. Bertsekas. Incremental Gradient, Subgradient, and Proximal Methods for Convex Opti- mization: A Survey. In Suvrit Sra. Sebastan Nowozin, and Stephen J. Wright (eds.). Optimization for Machine Learning. chapter 4. The MIT Press, 2011. ISBN 9780262298773.
Dimitri P. Bertsekas and John N. Tsitsiklis. Gradient Convergence in Gradient methods with Error.. SIAM Journal on Optimization, 10(3):277642, January 2000a. doi: 10.1137/
s1052623497331063
10

---

Published as a conference paper at ICLR 2026
Dimitri P Bertsekas and John N Tsitsiklis. Gradient convergence in gradient methods with errors. SIAM Journal on Optimization, 10(3):627-64,, 2000b.
Léon Bottou. Curiously fast convergence of some stochastic gradient descent algorithms. Un- published open problem offered to the attendance of the SLDS 2009 conference, 2009. URL probl probllen em 2mm-009.
hittp://e /tt.://elon.../ .oon.oooouu..... org/ pappers rs
Léon Bottou. Stochastic gradient descent tricks. In Neural Networks, 2012a.
Léon Bottou. Stochastic gradient descent tricks. In Neural networks: tricks of the trade: second edition, pp. 421 421-436. Springer, 2012b.
Léon Bottou, Frank Curtis, and Jorge Nocedal. Optimization methods for large-scale machine learning. SIAM review, 60(2):223-311, 2018.
Odile Brandière and Marie Duflo. Les algorithmes stochastiques contournent-ils les pieges? In Annales de I'IHP Probabilités et statistiques, volume 32, pp. 395-427, 1996.
Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra. Sampling from a log-concave distribution with compact support with proximal langevin monte carlo. In Conference on learning theory, pp. 3123342. PMLR, 2017.
Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from log-concave distribution with projected langevin monte carlo. Discrete Computational Geometry, 59(4):757-783, 2018.
Alexander Bukharin et al. Robust multi-agent reinforcement learning via adversarial regularization: Theoretical foundation and stable algorithms. In Advances in Neural Information Processing Systems (NeurIPS), volume 36, pp. 68121-68133, 2023.
Xufeng Cai, Cheuk Yin Lin, and Jelena Diakonikolas. Empirical risk minimization with shuffled sgd: A primal-dual perspective and improved bounds. AXXi:2330.112998, 2023.
Jerry Chee and Panos Toulis. Convergence diagnostics for stochastic gradient descent with constant learning rate. In International Conference on Artificial Intelligence and Statistics, pp. 1485- PMLR, 2018.
George HG Chen and R Tyrrell Rockafellar: Convergence rates in SIAM Journal on Optimization, 7(2):421-444, 1997.
splitting.
Xiang Cheng, Niladri S Chatterji, Yasin Yadkor. Peter Bartlett, and Michael Jor- dan. Sharp convergence rates for langevin dynamics in the nonconvex setting. arXiv preprint arXiv:1805.0164 2018a.
Xiang Cheng, Niladri Chatterji, Peter I. Bartlett, and Michael Jordan. Underdamped langevin meme: non-asymptotic analysis. In Conference on learning theory, pp. 300-323, PMLR, 2018b.
Hanseul Cho and Chulhee Yun. SGDA with shuffling: faster convergence for nonconvex-pl minimax optimization. In ICLR, 2023.
Sayantan Choudhury, Eduard Gorbunov, and Nicolas Loizon. Single-call stochastic extragradient methods for structured non-monotone variational inequalitics: Improved analysis under weaker conditions. In NeurIPS, 2023.
Kai Lai Chung. On a stochastic approximation method. The Annals of Mathematical Statistics, pp- 463-483, 1954.
Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(3):651- 676. 2017.
Amak Dalalyan and Avetik Karagulyan. User-friende frienffre guarantees for the langevin monte carlo with inaccurate gradient. Stochastic Processes and their Applications, 129(12):5278-5)) 2019.

---

Published as a conference paper at ICLR 2026
Arnak S Dalalyan and Lionel Riou-Durand. On sampling from log-concave density using kinetic langevin diffusions. 2020.
George B Dantzig and RW Cottle. Complementary pivot theory of. mathematical programming. Mathematics of the decision sciences, part, 1:115-136, 1968.
Aniket Das, Bernhard Schölkopf, and Michael Muehlebach. Sampling without replacement leads to faster rates in finite-sum minimax optimization. In NeurlPS, 2022.
Constantinos Daskalakis, Andrew Ilyas, Vasilis Syrgkanis, and Haoyang Zeng. Training gans with optimism. arXiv preprint arXiv:1711.00144 2017.
Constantinos Daskalakis, Stratis Skoulakis, and Manolis Zampetakis. The complexity of constrained min-max optimization. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 1466-1478, 2021.
Aymeric Dieuleveut, Alain Durmus, and Francis Bach. Bridging the gap between constant step size stochastic gradient descent and markov chains. 2020.
Yoel Drori and Ohad Shamir. The Complexity of Finding Stationary Points with Stochastic Gradient Descent. arXiv preprint arXiv:1910.0188. 2019.
Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai Gradient descent finds global minima of deep neural networks. In International conference on machine learning, pp. 1675- 1685. PMLR, 2019.
Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted langevin algorithm. 2017.
Alain Durmus, Umut Simsekli, Eric Moulines, Roland Badcau, and Gaël Richard. Stochastic gra- dient richardson-rombler markov chain monte carlo. Advances neural information processing systems, 29, 2016.
Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu. Log-concave sampling: Metropolis- hastings algorithms are fast. Journal of Machine Learning Research, 20(183):1-42 2019.
Konstantinos Emmanouilidis, René Vidal, and Nicolas Loizou. Stochastic extragradient with ran- dom reshuffling: Improved convergence for variational inequalities. In International Conference on Artificial Intelligence and Statistics, pp. 6682-3690. PMLR, 2024.
Murat A Erdogdu and Rasa Hosseinzadeh. On the convergence of langevin monte carlo: The inter- play between tail growth and smoothness. In Conference on Learning Theory, pp. 1776-1822. PMLR, 2021.
Murat A Erdogdu, Lester Mackey, and Ohad Shamir. Global non-convex optimization with dis- cretized diffusions. Advances in Neural Information Processing Systems, 31, 2018.
Vaclav Fabian. On asymptotic normality in stochastic approximation. The Annals of Mathematical Statistics, pp. 1327-1332, 1968.
Francisco Facchinei and Jong- Shi Pang. Finite-dimensional variational inequalities and comple- mentarity problems. Springer, 2003.
Yixin Fang, Jinfeng Xu, and Lei Yang. Online bootstrap confidence intervals for the stochastic gradient descent estimator. Journal of Machine Learning Research, 19(78):1-21, 2018.
Jean-Claude Fort and Gilles Pages. Asymptotic behavior of markovian stochastic algorithm with constant step. SIAM journal on control and optimization, 37(5):1456-1485 1999.
Angeliki Giannou, Kyriakos Lotidis, Panayotis Mertikopoulos, and Emmanouil- Vasileios Vlatakis- Gkaragkounis. On the convergence policy gradient methods to nash equilibria in general stochastic games. Advances in Neural Information Processing Systems, 35:7128-7144 2022.
12

---

Published as a conference paper at ICLR 2026
Paul W Goldberg and Matthew Katzman. Ppad-complete pure approximate nash equilibria in lips- chitz games. In International Symposium on Algorithmic Game Theory, pp. 169-186. Springer, 2022.
I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio. Generative adversarial nets. In NeurIPS, 2014.
Eduard Gorbunov, Nicolas Loizou, and Gauthier Gidel. Extragradient method: O (1/k) last-iterate convergence for monotone variational inequalities and connections with cocoercivity. In A/STA7S, 2022a.
Eduard Gorbunov, Adrien Taylor, and Gauthier Gidel. Last-iterate convergence of optimistic gradi- method for monotone variational inequalities. In NeurIPS, 2022b.
Robert Mansel Gower, Nicolas Loizou, Xun Qian, Alibek Sailanbayev, Egor Shulgin, and Peter Richtarik. Sgd: General analysis and improved rates. In AISTATS, 2019.
Luigi Grippo. A class of unconstrained minimization methods for neural network training. Opti- mization Methods and Software, 4(2):135-16) 150, January 1994. doi: 10. 10010//055679900005588888
Mert Gürbüüzaalaban, Asuman Özdağlar, and Pablo A. Parrilo. Why random reshuffling beats
stochastic gradient descent. Mathematical Programming, Oct 2019. ISSN 1436-4646. doi:
10.10077s100107-019-044444
Mert Gürbüübalaaan, Asu Ozdaglar, and Pablo A Parrilo. Why random reshuffling beats stochastic gradient descent. Mathematical Programming, 186:49-84, 2021.
Andi Han et al. Nonconvex-noncooeee min-max optimization on riemannian manifolds. Transac- fions Machine Learning Research, 2023.
Wenrui Hao. A gradient descent method for solving system of nonlinear equations. Applied Mathematics Letters, L12:106739 2021.
Jeff Haochen and Suvrit Sra. Random Shuffling Beats SGD after Finite Epochs. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 2624-2633, Long Beach, California, USA, 09-15 Jun 2019. PMLR.
Francis Begnaud Hildebrand. Introduction numerical analysis. Courier Corporation, 1987.
Ya-Ping Hsich, Panayotis Mertikopoulos, and Volkan Cevher. The limits of min-max optimization algorithms: Convergence to spurious non-critical sets. In International Conference on Machine PMLR, 2021.
Learning, pp.
Ya-Ping Hsich, Mohammad Reza Karimi Jaghargh, Andreas Krause, and Panayotis Mertikopou- los. Riemannian stochastic optimization methods avoid strict saddle points. Advances in Neural Information Processing Systems, 022295800296011 2023.
Yu- fu-Guan Hsieh, Franck Iutzeler, Jérôme Malick, and Panayotis Mertikopoulos. On the convergence of single-call stochastic extra-gradient methods. In NeurIPS, 2019.
Dongyan Huo, Yudong Chen, and Qiaomin Xie. Bias and extrapolation markovian linear stochas- tic approximation with constant stepsizes. In Abstract Proceedings of the 2023 ACM SIGMET- RICS International Conference on Measurement and Modeling of Computer Systems, pp. 81-82, 2023.
UC IEOR. The linear complementarity problem, lemke algorithm, perturbation, and the complexity class ppad. 2011.
Chi Jin, Praneeth Netrapalli, and Michael I. Jordan. What is local optimality in nonconvex- nonconcave minimax optimization? In Proceedings of the 37th International Conference Machine Learning (ICML), PMLR, 2020.
Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the fokker- planck equation. SIAM journal on mathematical analysis, 29(1):1-17 1998.
13

---

Published as a conference paper at ICLR 2026
Jack Kiefer and Jacob Wolfowitz. Stochastic estimation of the maximum of regression function. The Annals of Mathematical Statistics, pp. 462-466, 1952.
Yuri Kifer. Random perturbations of dynamical systems. Nonlinear Problems in Future Particle Accelerators, 189, 1988.
Beomsu Kim and Junghoo Sco. Semi-implicit hybrid gradient methods with application to adver- sarial robustness. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS), PMLR, 2022.
Galina M Korpelevich. The extragradient method for finding saddle points and other problems. Matecon, 12:747-756, 1976.
Jeongyeol Kwon, Luke Dotson, Yudong Chen, and Qiaomin Xie. Two-timescale linear stochastic approximation: Constant stepsizes go long way. arXiv preprint arXiv:2410.13067 2024.
Zehua Lai and Lek-Heng Lim. Recht-Ré Noncommutative Arilhnetie-Gcometrie Mean Conjecture is False. arXiv preprint arXiv:2006.01510, 2020.
Xuechen Li, Yi Wu, Lester Mackey, and Murat A Erdogdu. Stochastic runge-kutta accelerates langevin monte carlo and beyond. Advances in neural information processing systems, 32, 2019.
Tianyi Lin, Zhengyuan Zhou, Panayotis Mertikopoulos, and Michael Jordan. Finite-time last-iterate convergence for multi-agent learning in games. In ICML, 2020.
Chaoyue Liu, Dmitriy Drusvyatskiy, Misha Belkin, Damek Davis, and Yian Ma. Aiming towards the minimizers: fast convergence of sgd for overparametrized problems. Advances in neural information processing systems, 36:60748-60767, 2023.
Lennart L.jung. Strong convergence of a stochastic approximation algorithm. The Annals of Statis- tics, 6(3):680-696, 1978.
Lennart Ljung. Analysis of recursive stochastic algorithms. IEEE transactions automatic control, 22(4):551-575, 2003.
Nicolas Loizou, Hugo Berard, Alexia Dolicoeur-Martinee Pascal Vincent, Simon Lacoste-Julien, and Ioannis Mitliagkas. Stochastic hamiltonian gradient methods for smooth games. In ICML, 2020.
Nicolas Loizou, Hugo Berard, Gauthier Gidel, Ioannis Mitliagkas, and Simon Lacoste-Julien. Stochastic gradient descent-ascent and consensus optimization for smooth games: Convergence analysis under expected co-coercivity. In NeurIPS, 2021.
Zhi-Quan Luo. On the Convergence of the LMS Algorithm with Adaptive Learning Rate for Linear Feedforward Networks. Neural Computation, 3(2):226-24) 22::22266245,10 June 1991. doi: 19.1162/ncc.19911. 3.2.226.
Siyuan Ma, Raef Bassily, and Mikhail Belkin. The power of interpolation: Understanding the effec- tiveness of sgd in modern per-parametrized learning. In International Conference Machine Learning. pp. 3325-332 3334. PMLR, 2018.
Jérôme Malick and Panayotis Mertikopoulos. The global convergence time of stochastic gradi- ent descent in non-convex landscapes: Sharp estimates via large deviations. arXiv preprint: 2503.16396. 2024.
Olvi L. Mangasarian and Mikhail V. Solodov. Serial and parallel backpropagation convergence via nonmonotone perturbed minimization. Optimization Methods and Software, 4(2):103-116; 1994. doi: 10:088/105567884400555811
Paul Mangold, Alain Durmus, Aymeric Diculeveut, Sergey Samsonov, and Eric Moulines. Refined analysis of federated averaging's bias and federated richardson-romberg extrapolation. arXiv preprint 0rXiv:24/2.01/88, 2024.
Panayotis Mertikopoulos and Zhengyuan Zhou. Learning in games with continuous action sets and unknown payoff functions. Mathematical Programming, 173(1):465-50) 2019.
14

---

Published as a conference paper at ICLR 2026
Panayotis Mertikopoulos, Nadav Hallak, Ali Kavis, and Volkan Cevher. On the almost conver- gence of stochastic gradient descent in non-convex problems. Advances in Neural Information Processing Systems, 33:1117-1128, 2020.
Panayotis Mertikopoulos, Ya-Ping Hsieh, and Volkan Cevher. A unified stochastic approximation framework for learning in games. Mathematical Programming, 203(1):559 -09, 2024.
Lars Mescheder, Scbastian Nowozin, and Andreas Geiger. The numerics of gans. In NeurIPS, 2017.
Sean Meyn and Richard L Tweedie. Markov chains and stochastic stability. Springer Science & Business Media, 2012.
Konstantin Mishchenko, Ahmed Khaled, and Peter Richtárik, Random reshuffling: Simple analysis with vast improvements. In NeurlPS, 2020a.
Konstantin Mishchenko, Dmitry Kovalev, Egor Shulgin, Peter Richtárik, and Yura Malitsky. Revis- iting stochastic extragradient. In AISTATS, 2020b.
Dheeraj Nagaraj. Prateek Jain, and Pranceth Netrapalli. SGD without Replacement: Sharper Rates for General Smooth Convex Functions. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceed- ings of Machine Learning Research, pp. 4703 4711, Long Beach, California, USA, 09-15 Jun 2019, PMLR.
Hongscok Namkoong and John C Duchi. Stochastic gradient methods for distributionally robust optimization with f-divergences. In NeurIPS, 2016.
Angelia Nedić and Dimitri P. Bertsekas. Incremental Subgradient Methods for Nondiffcrentiabl Optimization, SIAM Journal on Optimization, 12(1):109-138, January 2001. doi: 10.1137/ $1052623499362111.
Deanna Needell, Rachel Ward, and Nati Srebro. Stochastic Gradient Descent, Weighted Sampling, and the Randomized Kaczmarz algorithm. In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger (eds.), Advances in Neural Information Processing Systems 27, pp. 1017-1025. Curran Associates, Inc., 2014.
Phuong Ha Nguyen, Lam Nguyen, and Marten van Dijk. Tight Dimension Independent Lower Bound the Expected Convergence Rate for Diminishing Step Sizes in SGD. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 3660-3669. Curran Associates, Inc., 2019.
Quynh N Nguyen and Marco Mondelli. Global convergence of deep networks with one wide layer followed by pyramidal topology. Advances in Neural Information Processing Systems, 33:11961- 11972, 2020.
1961-
Christos H Papadimitriou, Emmanouil- -ananouil-Vaaiiliios Vlatakis- Vlatakis-Gkaragken Gkaragkounss, and Manolis Zampetakis. The computational complexity of multi-player concave games and kakutani fixed points. arXiv preprint arXiv:2207.07557. 2022.
Robin Pemantle. Nonconvergence to unstable points in urn models and stochastic approximations. The Annals of Probability, pp. 698-712, 1990.
David Pfau and Oriol Vinyals. Connecting generative adversarial networks and actor-critic methods. arXiv preprint arXiv: 110.0096,0 2016.
Georg Ch Pflug. Stochastic minimization with constant step-size: asymptotic laws. SIAM Journal on Control and Optimization, 24(4):655-666, 1986.
Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM journal on control and optimization, 30(4):838-88) 855, 1992.
Priouret and A Yu Veretenikov. A remark on the stability of the Ims tracking algorithm. Stochastic analysis and applications, 16(1):119-129 1998.
15

---

Published as a conference paper at ICLR 2026
Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis. In Conference on Learning Theory, pp- 1674-1703. PMLR. 2017.
Shashank Rajput, Anant Gupta, and Dimitris Papailiopoulos. Closing the convergence gap of SGD without replacement. arXiv preprint arXiv:2002. 10400, 2020.
Alexander Rakhlin, Ohad Shamir, and Karthik Sridharan. Making gradient descent optimal for strongly convex stochastic optimization. arXiv preprint arXiv: / 109.5667, 2011.
Alexander Rakhlin, Ohad Shamir, and Karthik Sridharan. Making Gradient Descent Optimal for Strongly Convex Stochastic Optimization. In Proceedings of the 29th International Coference on International Conference on Machine Learning, ICML' 12, pp. 1571-1578, Madison, WI, USA, 2012. Omnipress. ISBN 9781450312851
Benjamin Recht and Christopher Ré. Toward a noncommutative arithmetic-geomen mean in- equality: Conjectures, case-studies, and consequences. In S. Mannor, N. Srebro, and R. C. Williamson (eds.): Proceedings of the 25th Annual Conference on Learning Theory, volume 23, pp. 11.1-11.24, 2012. Edinburgh, Scotland.
Benjamin Recht and Christopher Ré. Parallel Stochastic Gradient Algorithms for Large-Scale Ma- trix Completion. Mathematical Programming Compulation, 5(2):201-226, April 2013. doi: 10.10077611332-013-000533
Lewis Fry Richardson. Ix. the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Philosophical Transactions of the Royal Society London. Series A, containing papers of mathematical or physical character, 2 2100459-470:307--555 1911.
Herbert Robbins and Sutton Monro. A stochastic approximation method. The annals of mathemati- cal statistics, pp. 400-407, 1951.
Werner Romberg. Vereinfachte numerische integration. Norske Vid. Selsk. Forh., 28:30-36, 1955.
David Ruppert. Efficient estimations from slowly convergent robbins-monro process. Technical report, Cornell University Operations Research and Industrial Engineering, 1988.
Jerome Sacks. Asymptotic distribution of stochastic approximation procedures. The Annals of Mathematical Statistics, 29(2):373-405, 1958.
Itay Safran and Ohad Shamir. How good is sgd with random shuffling? In COLT, 2020.
Mark Schmidt and Nicolas Le Roux. Fast convergence of stochastic gradient descent under strong growth condition. arXiv preprint arXiv: 300.66330, 2013.
Alexander Shapiro. Asymptotic properties of statistical estimators in stochastic programming. The Annals of Statistics, pp. 458- 1989.
Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log log-ccecave sampling. Advances in Neural Information Processing Systems, 32, 2019.
Marina Sheshukova, Denis Belomestny, Alain Durmus, Eric Moulines, Alexey Naumov, and Sergey Samsonov. Nonasymptotic analysis of stochastic gradient descent with the richardson-romber extrapolation. arXiv preprint arXiv:2410.05106. 2024.
Chaobing Song. Zhengyuan Zhou, Yichao Zhou, Yong Jiang, and Yi Ma. Optimistic dual extrapo- lation for coherent non-monotone variational inequalities. Advances in Neural Information Pro- cessing Systems, 4440--133144 2020.
Matthew Staib, Sashank Reddi, Satyen Kale, Sanjiv Kumar, and Suvrit Sra. Escaping saddle points with adaptive gradient methods. In International Conference Machine Learning, pp. 5956- 5965. PMLR, 2019.
Guido Stampacchia. Formes bilineaires coercitives sur les ensembles convexes. Comptes Rendus Hebdomadaires Des Seances De I. Academie Des Sciences, 258 258(18):4413, 1964.
16

---

Published as a conference paper at ICLR 2026
Weijie Su and Yuancheng Zhu. Statistical inference for online learning and stochastic approximation via hierarchical incremental gradient descent. arXiv preprint arXiv:1802.04887, pp. 3, 2018.
Ruo-Yu Yu Sun. Optimization for Deep Learning: An Overview. Journal of the Operations Research
Society of China, 8(2):299-294, Jun 2020. ISSN 2194-6698. doi: 10. 10.11077/000000--22-0000066666
Vasilis Syrgkanis, Alekh Agarwal, Haipeng Luo, and Robert E Schapire. Fast convergence of regu- larized learning in games. Advances in Neural Information Processing Systems, 28, 2015.
Denis Talay and Luciano Tubaro. Expansion of the global error for numerical schemes solving stochastic differential equations. Stochastic analysis and applications, 8(4):483-509, 1990.
Yan Shuo Tan and Roman Vershynin. Online stochastic gradient descent with arbitrary initialization solves non-smooth, non-convex phase retrieval. Journal of Machine Learning Research, 24(58): 1-47, 2023.
Panos Toulis and Edoardo M Airoldi. Asymptotic and finite-sample properties of estimators based on stochastic gradients. 2017.
Nilesh Tripuraneni, Nicolas Flammarion, Francis Bach, and Michael Jordan. Averaging stochastic gradient descent riemannian manifolds. In Conference On Learning Theory, pp. 650-687. PMLR, 2018.
Sharan Vaswani, Francis Bach, and Mark Schmidt. Fast and faster convergence of sgd for over- parameterized models and accelerated perceptron. In The 22nd international conference on artificial intelligence and statistics, pp. 1195-1204. PMLR, 2019.
Emmanouil Vasileios Vlatakis- Sanaaiii-GGarrakouunii, Angeliki Giannou, Yudong Chen, and Qiaomin Xie. Stochastic methods in variational inequalities: Ergodicity, bias and refinements. In International Conference Artificial Intelligence and Statistics, pp. 4123-4131. PMLR, 2024.
Jingkang Wang, Tianyun Zhang, Sijia Liu, Pin- Yu Chen, Jiacen Xu, Makan Fardad, and Bo Li. Ad- versarial attack generation empowered by min-max optimization. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (eds.), NeurlPS, 2021.
Wikipedia. http:///en.wikipedii.rrr.. Sikkippdiiorg/wik///aassstttttt metric, Accessed: 2025-08-28.
Wasserstein metric wikipedia, the free encyclopedia, 2023.
URL
Lu Yu, Krishnakumar Balasubramanian, Stanislav Volgushev, and Murat A Erdogdu. An analysis of constant step size sgd in the non-convex regime: Asymptotic normality and bias. Advances in Neural Information Processing Systems, 34:4234-444 -2444444.,.... 2021.
Kaiqing Zhang, Zhuoran Yang. and Tamer Başar. Multi-agent reinforcement learning: A selective overview of theories and algorithms. Handbook of reinforcement learning and control, pp. 321- 384, 2021.
Yixuan Zhang and Qiaomin Xie. Constant stepsize learning:: Distributional convergence, bias and extrapolation. arXiv preprint arXiv:2401.13884, 2024.
Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Gradient descent optimizes over- parameterized deep relu networks. Machine learning, 109(3):467-492, 2020.
17

---

Published as a conference paper at ICLR 2026
SUPPLEMENTAL MATERIAL
L
CONTENTS
Appendix
18
A Additional Related Literature
19
Proof Roadmap
23
B.1 Warm-Up:Usefull Inequalities
23
C Proof of Convergence Rate (MSE) of Perturbed SGD-RR (Theorem 4.1)
24
C.1 Preparatory Lemmas & Propositions
24
C.2 Assembling the Lemmas: Proof of Theorem 4.1
26
D Proof of Ergodic Properties and Limit Theorems (Theorem 4.3) & (Theorem 4.4)
28
D.I Proof of Continuous State Time Homogenious Markov Chain (Lemma 4.2)
28
D.2 Proof of Foster-Lyuapunov inequality
29
D.3 Proof of Minorization Property
D.4 Proof of Irreducibility, Aperiodicity and Harris and Positive Recurrence
32
D.5 Proof of Existence of Unique Invariant Distribution at Epoch-level (Theorem 4.3)
32
D.6 Proof of Limit theorems of Epoch-level iterates (Theorem 4.3)
34
Bounding 4th-order Moments of Error Distance
35
E.1 Higher order version of Proposition A.2 (Emmanouilidis et al., 2024) Bound of Kyrtosis for Lipschtiz Operators
35
E.2 A Combinatorial Tool for Higher- Order Moments In Sampling without Replacement
36
F Proof of Bias for RR, and RR RR, Methods (Lemma 4.5) and (Theorem 4.6)
42
F.1 Jacobian bound for one-pass map
42
F.2
Higher-Order Terms of RR, bias (Lemma 4.5)
44
F33 Proof of Bias Refinements of @RR (Theorem 4.6)
45
On Experiments
47
18

---

Published as a conference paper at ICLR 2026
A ADDITIONAL RELATED LITERATURE
Owing to its central role in optimization and machine learning, stochastic gradient descent (SGD) and its numerous variants have generated an extensive body of literature that spans several decades. A complete survey is well beyond the scope of this paper, so we restrict ourselves to highlighting the most relevant threads and pointers.
The only aspects we emphasize at this point are the phenomena most pertinent to our analysis:
1. classical stochastic approximation and asymptotic normality results;
2. constant step-size schemes viewed through the lens of Markov chains and diffusion ap- proximations;
3. the widespread use of random reshuffling and its still-developing theoretical guarantees; and
4. challenges that arise in max and variational inequality settings.
These themes form the backbone of our extended discussion in the appendix, where we provide a more comprehensive account of prior work.
From classical stochastic approximation to modern SGD. The study of stochastic approxima- tion predates machine learning by decades, beginning with the foundational work Robbins & Monro (1951) and Kiefer & Wolfowitz (1952). Early analyses focused vanishing step-sizes obeying the classic L summability rules, and developed the ODE method to describe limiting dynamics; see, e.g., Ljjung (1978; 2003); Benaïm (2006); Bertsekas & Tsitsiklis (2000b). In paral- lel, a rich line of results examined the almost-sure hehavior of stochastic approximation, including avoidance of saddle points and convergence locally stable equilibria (Pemantle, 1990; Brandière & Duflo, 1996; Benaïm & Hirsch, 1995; Hsieh et al., 2021; 2023; Jordan et al., 1998; Mertikopoulos al., 2020; 2024; Staib et al., 2019; Antonakopoulos et al., 2022).
Asymptotic normality and statistical inference. A complementary thread established central limit theorems for stochastic approximation: classical milestones include Chung (1954); Sacks (1958); Fabian (1968); Ruppert (1988); Shapiro (1989), culminating in the Polyak-Juditsky av- eraging principle (Polyak & Juditsky, 1992). Under suitable decaying step-sizes, the averaged SGD iterate is asymptotically normal and attains the Cramér-Rai optimal variance. This statistical per- spective has been leveraged to construct confidence intervals and inference procedures for SGD- based estimators (Tripuraneni al., 2018; Su & Zhu, 2018; 2018).
Toulis & Airoldi, 2017; Fang et al.,
Constant step sizes: bias, speed, and Markovian viewpoints. Constant step-size policies, now standard in large-scale learning, trade a nonvanishing asymptotic error for fast initial progress and robust practical performance. Their benefits in over-parameterized regimes are well documented (Schmidt & Roux, 2013; Needell et al., 2014; Ma et al., 2018; Vaswani et al., 2019). The Markov chain viewpoint provides a unifying language for analyzing such constant-step schemes: schemeess early devel- opments used dynamical-systen properties (Kifer, 1988; Benaim, 1996; Priouret & Veretenikov, 1998; Fort & Pages, 1999; Aguech et al., 2000), with recent refinements quantifying convergence behavior and variance (Dieuleveut et al., 2020; Chee & Toulis, 2018; Tan & Vershynin. 2023). In parallel, diffusion-basee analyses and Langevin-type discretizations connect SGD to MCMC methodology, yielding non-asymptotic
and Markov-process: techniques to establish stability and ergodic
guarantees and sharp mixing rates in log-concave and beyond-log-concave settings (Dalalyan, 2017;
Durmus& & Moulines, 2017; Cheng et al., 2018b; Dalalyan & Karagulyan, 2019; Brosse al., 2017; Cheng et al., 2018a; Bubeck et al., 2018; Dwivedi et al., 2019; Dalalyan & Riou-Durand 2020; Li al., 2019; Shen & Lee, 2019; Erdogdu & Hosseinzadeh, 2021).
Random reshuffting vs. with-replacement sampling. Among finite-sum methods, random reshuffling (RR) occupies a special place: each epoch processes every component exactly once in random order, in contrast to classical with-replacement SGD. RR is ubiquitous in practice- it improves cache locality (Bengio, 2012), often converges faster than with-replacement sampling (Bottou, 2009; Recht& & Ré, 2013), and is the default in deep learning pipelines (Sun, 2020). The
19

---

Published as a conference paper at ICLR 2026
success of RR contrasts with the mature theory for with-replacement SGD, which enjoys tight up- per/lower bounds in many regimes (Rakhlin et al., 2012; Drori & Shamir, 2019; Nguyen et al,, 2019). A key analytical hurdle is bias: within epoch, conditional expectations are no longer unbiased gradients, so classical SGD proofs do not transfer verbatim. Early attempts leveraging the noncommutative crrthnetic-gommetine mined when the conjecture was disproved (Lai & Lim, 2020). More recent works establish rates for twice-smooth and smooth objectives and highlight gaps between theory and prevalent heuristics (Gürbüzzbalaban et al., 2019; Haochen & Sra, 2019; Nagaraj et al., 2019; Safran & Shamir, 2020; Rajput et al., 2020).
mean conjecture (Recht & Ré, 2012) were later under-
Incremental/ordered passes and sensitivity to permutations. Before RR became the default, incremental gradient (IG) methods with fixed orderings were widely used in neural network training (Luo, 1991; Grippo, 1994), with asymptotic convergence known since early work (Mangasarian & Solodov, 1994; Bertsekas & Tsitsiklis, 2000a). However, their performance can depend strongly on the chosen ordering (Nedić & Bertsekas, 2001; Bertsekas, 2011). By randomizing the order every epoch, RR mitigates this sensitivity and-under smoothness-caa outperform both with- replacement SGD and deterministic IG (Gürbüzbalaban al., 2019; Haochen & Sra, 2019), with refined lower/upper bounds developed in follow-up studies (Nagaraj et al., 2019; Safran & Shamir, 2020: Rajput et al., 2020).
Min-max problems and variational inequalities. In large-scale saddle-point and VIP settings, most theoretical analyses assume with-replacement sampling for convenience, whereas implemen- tations overwhelmingly adopt without-replacement sampling (Bottou, 2012a). A growing literature is closing this gap: for minimization problems, several works show (sometimes provably faster) RR rates in finite-sum regimes (Mishchenko al., 2020a; Ahn et al., 2020; Gürbüzbalaban al., 2021; Cai et al., 2023). For min-max and VIPs, guarantees remain comparatively sparse: Chen & Rockafellar (1997) and Korpelevich (1976) initiated the study of stochastic and extragradient-type methods, with modern analyses for SEG and optimistic variants (Gorbunov et al., 2022a;b; Hsich al., 2019 Choudhury et al., 2023). For RR specifically. Das et al. (2022) derive guarantees for SGDA and PPM under strong structural conditions, and Cho & Yun (2023) extend to certain non- monotone settings. Nevertheless, classical SGDA can diverge even in simple monotone bilinear games, while proximal methods are implicit and less practical; filling this theoretical-praat -praciccl gap remains an active direction.
Overparameterization and global convergence phenomena. Finally, SGD training dynamics in overparameterized neural networks reveal regimes where global convergence can emerge from struc- tural properties such as width, depth, and initialization (Du et al., 2019; Zou et al., 2020; Nguyen & Mondelli, 2020; Liu et al., 2023). These results are powerful but specialized: they rely on problem- specific structure (c.g., width scaling or tailored initializations). Our focus is orthogonal seek guarantees for general non-convex or non-mone monoonet landscapes under stochastic approximation, in- dependent of architectural assumptions. For completeness, we refer the reviewer for the related work of the aforementioned work for further surveys about these SGD & overparameterization results in more detail.
Comparison to Prior Work and Overview of Our Contributions. Before introducing the intu- ition behind our algorithmic design, we briefly contrast our results with those of Emmanouilidis et al. (2024), who study RR improvements for the Stochastic Extragradient (SEG) method. Our analysis uncovers fundamentally different phenomenon: the joint use of RR, + RR₂ produces bias cancellation mechanism that eliminates the leading O()) term while preserving the condi- tion number and asymptotic behavior of SGDA. The key distinctions are summarized Table 1. Achieving the best of both worlds- bias order together with tight condition number, as SEG without reshuffling enjoys-remaing an interesting direction for future work.
Summary. To summarize, there is mature theory for with-replacement SGD (both asymptotic and non-asymptotic), well-developed statistical limits via averaging, and powerful diffusion/Markov perspectives for constant-step schemes. RR, despite being the practical default, poses distinctive analytical challenges due to its within-epoch bias, especially in min-max max and VIP settings. Re- cent advances begin to bridge this gap, but a comprehensive understanding of how classical heuris-
20

---

Published as a conference paper at ICLR 2026
ties (constant steps, reshuffling, extrapolation) interact remains incomplete precisely the juncture where our work contributes.
OUR MODEL'S ASSUMPTIONS.
While variational incqualities provide unifying language for optimization, learning, and game dy- namics, no single structural assumption can capture the full complexity of all modern nonconvex- nonconcave problems. From computational standpoint, full generality- -bing tightly connected to Nash equilibria (Papadimitriou et al., 2022; Goldberg & Katzman, 2022), linear complementarity (IEOR, 2011), and constrained saddle-point prob- lems (Daskalakis et al., 2021). Consequently, much of the theoretical literature adopts structured assumptions (strong convexity, quasi-strong monotonicity, quasar or weak convexity, PL/KL condi- tions, Minty conditions, error bounds, etc.), each expressive in specific regimes but not universal.
smooth VIs intractable in
Our work is based on quast-strong monotonicity which falls squarely within this class: it captures stabilizing behaviors of many smooth systems, while remaining far more permissive than strong convexity or global monotonicity. not meant universal model for all adversarial or fully nonconvex-nonconee modern ML applications- including GANs, adversarial robustness, and multi-agent RL-can ex- hibit fundamentally unstable or rotational dynamics (Jin et al., 2020; Han et al., 2023; Kim & Seo, 2022; Bukharin et al., 2023), where even local monotonicity surrogates fail. As such, our theoretical guarantees should be viewed pertaining regimes where minimum amount of local structure is present, rather than to the most adversarial or unstructured cases."
the same time, it is helpful clarify that this assumption
settings. Certain
LIMITATIONS
Limitations of the structural assumption. A central structural assumption in analysis is that the operator F satisfies X-weak j--uasi-stronn monotonicity. This condition broad enough to include several meaningful non-monotone problem classes-such as dissipative dynamical systems, weakly convex optimization, and locally contractive variational inequalities-and is standard in modern analyses of stochastic fixed-point and operator-slliting methods (e.g., Hsieh et al., 2019; Mertikopoulos and Zhou, 2019; Chavdarova et al., 2021). However, it is important to emphasize the following limitations.
1.
Not applicable to general nonconvex-noncoocaax min-max problems. The assump- tion does not hold for arbitrary adversarial problems such as deep GANs, multi-agent RL environments, smooth non-monotone games with persistent rotational dynamics. These settings may lack even local stability (e.g.. Daskalakis et al. 2018; 8; Fiez et al. 2020). Accord- ingly, our theoretical guarantees should not be interpreted applying to fully adversarial or worst-case min- max formulations.
2. Local nature of the assumption. Quasi-strong monotonicity is inherently a local regular- ity condition: smooth operators that are monotone in neighborhood of a solution satisfy it on that region (Lemma A.4, Hsieh et al. 2019). This requires smoothness and regularity that may not hold in problems involving discontinuities, clipping, piecewise-linee linear losses, or hard constraints. In such cases, the assumption may fail even locally.
3. Not capturing highly oscillatory or anti-monotone operators. Allowing À per- mits controlled violation of monotonicity, but the assumption still does not model strongly anti-monotone or highly oscillatory operators. Extending our analysis to Minty variational inequalities, hypomonotone operators, or other generalized monotonicity classes remains an interesting direction for future work.
'A complementary and key fact for our setting, established in (Hsich al., 2019, Lemma A.4), that any smooth VI operator is locally quasi-strone stuiily monotone in a neighborhood of a regular solution. Combined with our Harkov-chain recurrence result whhihhwhchhhhh ensures that the iterates remain in such neighborhoods with probability 1hsh provides natural and widely adopted stability regime in which the RR; and RR2 dehiasing mcchanisms both theoretically justified and practically meaningful.
21

---

Published as a conference paper at ICLR 2026
Despite these limitations, believe the assumption remains meaningful for broad set of struc- tured, smooth VIPs where local stability is present. We hope that this explicit discussion helps avoid any misunderstanding about the scope of upplicability of our results.
22

---

Published as a conference paper at ICLR 2026
B PROOF ROADMAP
Our main theorem relies on several technical components, developed across different parts of the appendix. In Appendix C, we establish the mean-squared convergence rate of Perturbed SGD- RRj While this result is of independent interest-as it extends prior analyses to noisy setting-it primarily serves construct the Lyapunov function that underpins our Markov chain treatment of the algorithm. Armed with this Lyapunov function, the properties of the perturbation, and the epoch- level viewpoint, Appendix D shows that SGD-RS forms a geometrically ergodic Markov chain with all standard consequences: existence of an invariant measure, law of large numbers, and a central limit theorem.
Appendix E develops higher- momeett control. In particular, part E.2 introduces a new combinatorial lemma on fourth moments of finite-sum subsets of vectors-a technically challenging step, moti- vated by the fact that most extrapolation analyses (e.g.. Dieuleveut et al. (2020)) require bounded fourth moments of the reshuffing estimator. With this tool in hand, Appendix F shows that no change in the step-size order is required to accommodate the extrapolation trick: we able to con- trol the higher moments of the Jacobian of the reshuffled biased gradient estimator. To the best of knowledge, this is the first such result. The last parts of Appendix F then contain the full proofs of our main theorem.
Finally, Appendix G presents additional experiments demonstrating the practical gains of our method, which originally motivated this study.
Appendix C Perturbed SGD-RR MSE rote & Lyapunov
Appendix
Appendix D measure, LLN,, CET
Epoch-level Markov
Minorizacion
Appendix Experiments
Foster-Lyapunon
Appendix E.2 4th ith-moment control combinstorial
Appendix Spectral analysis
Appendix F Blus analysis SGD-RR
Appendix Synergy RRresh + Extrap. Bias Olyy)
full--poss operntor
Figure 4: Dependency graph of Appendix results. Solid lines: main logical flow. Dashed lines:
auxiliary inputs.
B.1 WARM-UP:USEFUL INEQUALITIES
We start our technical appendix by providing inequalities that will be useful in proofs
Σ
02
12
2
(2)
iml
n
34
<
3
||0;||
(3)
1
|la
611
>
||a||
||611
(4)
2
(a, b)
=
||a|| 2 + [|6|1
||a -b1
(5)
2
-=
x, Vo 2 0
(6)
1
|/a b112
S
|1a12 +
|1b12 Vt (0.
23

---

Published as a conference paper at ICLR 2026
C PROOF OF CONVERGENCE RATE (MSE) OF PERTURBED SGD-RO R (THEOREM 4.1)
Our first result concerns the Perturbed SGD-RR, variant for X-weak u-quasi strongly monotone VIPs.
Theorem C.1 (Restatement of Theorem 4.1). Let Assumptions 2.1-2.3 hold. Then the iterates of Perturbed SGD-RRj satisfy for S Tmazz
[1 2z []
7mp
1++1
|20
x²
x²112
8n,
L
man
+
82
2
4²
LL662 Liman
where o. = E [x(x") and Ymar
min
In
12n
We first provide some notation that will be necessary for establishing the proof of Theorem 4.1. Consider the epoch-wise update rule of Perturbed SGD-RRJ
n-1
wg (x1)) -Uk
Ib
= 22 Gw, - UUk
(8)
where Gwa (x%) =
Σ Pwj (x1) denotes the epoch-wise operator used update the epoch-level
iterates (rk)
C.I PREPARATORY LEMMAS & PROPOSITIONS
With this notation at hand, we proceed in proving two Lemmas that are necessary for deriving the rate of convergence of the Theorem 4.1. In the first lemma, following the high-level intuition that one epoch random reshuffling with step size progresses the underlying dynamics approximately equal one step of the deterministic GD with step size y' ny, in the first lemma we bound the "progress" that the deterministic algorithm makes in one step.
Lemma C.2. Let Assumptions 2.1- hold. For any x² E X², the iterates of Perturbed SGD-RR, satisfy that
|k+1 - - -mFF((()) Fi
[[1 2mp)² + 22 mar || - x=11² 2 - 2ynA
Proof. For any fixed x² E X it holds that
|| - x² - -mF((x))1
| x. 111 22n -
- x², F(xxx) + +
VA IA
|| x²1 22mpp1114- x²1| 2ym 72n
(1 27mp) |00 x=11 + 2 | F(xx))12
Assumption
(1 2ym + + 1.2aa) x: x11 2ynA
(9)
Taking expectation condition on the filtration Fk. gives
E [| x= ->mF(xx)]1 [F
< (1 2ymp 2,,
- x1112 1
2ynA
I (1 - 72p)) + 22²² L || - x²1 2ynd
Having an expression for the progress made by the deterministic counterpart of Perturbed SGD-RR), we next aim to bound how large the deviation of the two algorithms becomes inside an epoch. To do so, we bound the sum of the distances of the iterates obtain by Perturbed SGD-RR, from the start of the epoch, which corresponds to the fictitious iterate of our comparator deterministic counterpart. The following lemma provides an upper bound dependent on the distance of the current epoch-level iterate from the solution and the variance at the optimum.
24

---

Published as a conference paper at ICLR 2026
Lemma C.3. Let Assumptions 2.1, 2.3 hold. If Perturbed SGD-RR, is ran with step size S then it holds that
V3n (n- (5n(n=))mmas lmma
2 EE
F
I
, maz | - x²1|²
-x"1|" + 2n²
Proof. From the epoch-level update (8), it holds
i-1
| -
x211
72,
Σ
F
(xx)
-=
IA@
i-1
2
ト (= ((x)))())))))))))
3721
Σ
Fog (xx)
Fug Fug
+
5=0
+3, =
Assumption 2.3
37 maxd
E
+ 3
Σ
F Fag F(22)
1=0
+3
where at the last step we have used the Lipschitz property of the operators F,, Vi expectation condition on the filtration Fk. we get
[n]. Taking
E || - =ill FR
I
3 L2
E
Σ
x
2
Fr
man
2
1 Fug(()))) F(x2) (((b))
2 E
F (337)FF(0))1111
From Lemma A.3 in (Emmanouilidis et al., 2024), it holds for A 1 E [((("))] I and vi E [n] that
2 Σ L. of
i(n 1)
E
Σ
Fug (x F(22)
Fr
n (Alla) x²11 x=112+20 +20?)
(11)
From inequality (11) and (10), thus, obtain
| 2
211
<
L2 2maz
Σ 1
=1 2
Fx
+ 3. 33,22(1 i(n - (( i)
n n- All - x²11
1=0
-, 2 i(n
a2 + 3y 3 F(
(12)
n
By summing over 0 Sis 1 we have that
Σ a - =021 Fb
S
12 maa
n2n1 1) n-1
2 Σ
F
||
2 [
22 n(n +
10 - x." x²1
2 + +1)03 +
2nn(n
1)(2n 1)
||
(13)
2
where we used the facts
n-1
Σ n - 1)
n-1
n(n - 1)(2n - 1)
a-
i(n - 1)
n(n + 1)
Σ
12
Σ
2
6
n 1
6
1=0
1=0
25

---

Published as a conference paper at ICLR 2026
For 7 $
rearranging the terms in (13) obtain
/in(n-I)/mma 1)Lmas
n-
Σ
||ا - =211 =
7 7nn(n + 1)A|| - =²1 2n(n + 1)
+n(n - 12 1) 1)7211 F(
Assumption 2.3 S
2y (A+ + n - x² =²1² +2n2n00
AS2L
6n 2 mmaz || x²1|² 2n2
In the preceding subsection, we established a series of preparatory lemmas. We now combine these ingredients into a unified elementwise argument prove Theorem 4.1.
C.2 ASSEMBLING THE LEMMAS: PROOF OF THEOREM 4.1
In this section, we provide the proof of Theorem 4.1, establishing linear convergence of Perturbed SGD-RR; to neighbourhood of the solution. The proof technique leverages the interpretation that one epoch of Perturbed SGD-RR, with sufficiently small step size > 0 is equivalent to one step of the gradient descent with step size 7 ny, as the iterates of Perturbed SGD-RR, inside the epoch do not change drastically. To account for the deviation of the iterates from the initial state z inside the epoch. we have upper bounded the sum of the corresponding distances in Lemma C.3. Thus, using the combining the bound on the progress made by gradient descent from Lemma C.2 with the potential "deviation" between the two algorithms we establish the rate convergence of the method.
Proof. Using the update rule of Perturbed SGD-RR, we have that:
2 k+1
ox n-1
yFutt
(xx-1)
UU
2 Fug (xx) UUx
n-1
2x (mF - Σ
Fug (xt) (xx()))) - W3
(14)
where the last step we used the fact that ymF = ? n-1
Σ Fut (x2) and the finite-sum structure
of the operator F. It holds, thus, that
2
| x²112
ok x² YmFFFFF
Σ
(Fag
pk (x2) U0
(15)
==
From Young's 's inequality, the right-hand side (RHS) (15) can bounded as follows
xx x² (2)(
m-1
|| - x01
+
?
Σ (Fut (0t)
Fry (xf)) + U
1 1 1mm mm
1241
16
0
x² ->FF(a2)11
27
Σ Fug (xf) Fug (22)
+ 27
|| 112
-mpp
mz
ni
IAB
x x² -mFF(22)11
27 n-
Σ
Fag (af) - Fut (0/))
+ 2 | 11 nj
+
1 1 1-mm44 7m144
/1
0=0
26

---

Published as a conference paper at ICLR 2026
Applying the Lipschitz property of the operators, obtain
| x01
| x - -nnF(x))1
2 L max
Σ
1
=211²
27
||u&|1 (16)
mmp
np
Taking expectation condition on the filtration F& (history of z²)) and using the fact that the noise Up Nf (0, we get
11 x² - -mF(22)1
|F
|| =111
2 maz
E
Σ i=1
1
28
1
|Fa
1 ymp
2n
(17)
/
To complete the proof, it suffices to bound each term on the right-hand side of (17). From Lem- mas C.2, C.3, holds for & /3n(n-I)Lmar 1)Lmax
that
F 1 x² - (xF(2)))1 |F S [1--7140) + 2= [.aa]] |ate x²12 + 2yn (18)
|ot - 2 Fx
S
4y2n² L - xz
(19)
Substituting (18) and (19) into (17), we obtain
(1-714- + 7² Limaz
8n 2 Lmaa
[110221-22111 - -=*1 s
1 x²1
1 7pp
4n L Limae
2n7A
a² + 1 7mp
(20)
11662 H
Selecting the stepsize 7 S min
we have that
2np
12nT.2nux
I 2
ymp
and (1 (1-7mp +
?2.2122 maa
32. max
11
7ml
1 mmp
2
and thus substituting in (20) we get
E |102 - x=1 12
I (1
7m 2
||
==112
4n² L2
o² + 4nyA
(21)
Taking expectation both sides and using the tower property of expectations, have that:
E1102 211
(1
2mp
|
x= 11
+
3 Lmex
L2 (1 qm
o² 8nyA
mmp
A+1
|122
x11
4nn Lmax 02 + 8ny2
mp
1++1
1x0 x x11
8nn L
max
+
8A
2
112
27

---

Published as a conference paper at ICLR 2026
D PROOF OF ERGODIC PROPERTIES AND LIMIT THEOREMS (THEOREM 4.3) & (THEOREM 4.4)
We start by proving a series of properties that the induced Markov Chain satisfies, that will be necessary for proving the Theorem 4.3.
Proposition D.1 (Proposition 5.5.3 (Meyn & Tweedie, 2012)). If a set C E B(R") is Vm -smal,,, then it is v.... petite for some 6m 0.
Intuition. The notions small and petite sets are technical tools in Markov chain theory that help verify stability properties. A set C is called small if, starting from C, the chain has a uniform positive chance of reaching any region of the state space within a fixed number of steps. A petite set weaker concept: instead of requiring such uniformity in single time step, it allows the chance of hitting any region to be distributed over random number of steps (via a probability distribution over times). Thus, every small set is automatically petite, but the reverse not true. Intuitively, small sets guarantee "uniform mixing after fixed horizon," while petite sets guarantee the same effect "on average over time."
D.I PROOF OF CONTINUOUS STATE TIME HOMOGENIOUS MARKOV CIIIIN (LEMMA 4.2)
Lemma D.2 (Epoch-level homogeneity and kernel). Fix 7 and N. Then Perturhed-SGD can be described at each epoch k as: Draw we uniformly from G,, and set
Fx+1 = H(r + Uk: Ux ~ ~O, Σ)
where H Hx,, w) denotes the endpoint of one reshuffled pass started at æ with permutation w (i.e., the map induced by n inner updates with step size 7). Then (rx)220 time-homogeneous Markov chain on Rd with transition kernel
P(s, EEE
1
J
Σ
Ays (I,WW, Σ) dy,
A A
n! WEE,
where 6(;m,)) is the d-variate Gaussian density with mean m and covariance Σ.
Proof. Fix 7 > 0 and n E N. For any Rd and w Gn, define the inner-epoch recursion
[0 (I. () = x, and the (measurable) epoch map
x (x, = 2it (x, - Fuli (zobi) (z. w)), j = 0,... n 1,
H(x,w) := Σ Puls (xbl(x
By construction, at epoch & the algorithm updates as
Fk+1 H (z wk) Uk-
where (wx )22 are i.i.d. uniform on G, and (Uk)a20 are i.i.d. with law (00 Ela), independent of (wk)k20 and of øb given the present state.
Markov property. Let A E B(Rd). Using the tower property and the independence of wk, Up from the past given Ik.
Pr(x441 E A 2o To. IK) = Pr(H (r(H wx) + Up E A Ik, wx)
2x = E[PPH (2x, w) + e A)],
where the outer expectation is over Unif((,,) and U ~ (O0, El,), independent. Thus
Pr(se+1 E A ;o. Fx) Pr(x++1 E A mx) P(*,, A),
so (zz))220 is Markov chain.
28

---

Published as a conference paper at ICLR 2026
Tumo-homo and kernel. Since the joint law of (Wk, Uk) does not depend on k. the transition kernel P is time-invariant By conditioning on w and integrating over the Gaussian U,
1 P(z, A)
Σ ((]H(Z,))+US TO! WEE..
1 Σ WEG.
/ olys (I,,w), Eld) dy,
where o(: m, Ela) denotes the d-variate Gaussian density with mean m and covariance Elg- This yields the stated expression for P and establishes time-homogeneity Rd
Augmented formulation (for reference). On the product space Rd x G,, define the next permutation w' Unif (G.) independently of (z,w) and U. Then the augmented chain ((r,, satisfies
(x',w))
w'),
and the associated kernel
|
(((x,w)AAx B) =
oly Hx, Ela) dy
which is manifestly homogeoooge.euus.
We, next, show that there exists energy function that describes the iterates of the Markov chain.
D.2 PROOF OF FOSTER-LYUAPUNOT INEQUALITY
central tool for proving stability and ergodicity of Markov chains is the Doster-Lyapune inequal- ity. The idea is construct "energy" or "Lyapumov" function Exx,x") that tracks the distance of the chain's state from equilibrium. If this function decreases on average outside a bounded region, it ensures that the process cannot drift to infinity and will instead return frequently to compact set. This property, when combined with the minorization condition, implies positive Harris recurrence and geometric ergodicity (Meyn & Tweedie, 2012).
In our case, a natural candidate for such an energy is the squared distance to a solution x", up an additive constant. The following corollary verifies that this choice indeed satisfies Fos- ter Lyapunoo inequality for Perturbed SGD-RR,, showing that the expected energy after one epoch contracts linearly up to fixed additive term.
Corollary D.3. Let Assumptions 2.1-2.3 hold. The function E(x x = |a x=1 + satisfies for any x² E X the inequality
E
k+ E[E(1b*1
,x*)
C x²) + c2.
where
20 and 62
2n
8mn 8 2 more o²
82
Proof. From inequality (21) of Theorem 4.1, we have that
El1lop 101 |
(1- mm4 k+1
11=0 x²112
8n
max
0, +
8A
n
p
Adding in both sides one and using the definition of E(x x*), we obtain
E [1124 - x + 1 1|
S (1- 7mp
1--2mn) (11a - 11 +1 yoy
8ny 2 1 8A
max
0 +
p
E(E(xbt) ,,,,)
C3E(2t x²) 02
(22)
where at the last step we have let cj =
7mp 2
and C2 = 2n
8ng 2 2
o² + 8A
2
29

---

Published as a conference paper at ICLR 2026
Lemma D.4. Let Assumptions 2.1-2.3 hold. If y Imaz, then for any fixed xt E X the functions E (zz x²) E(z, x² (x, x²) x() satisfy the geometric drift property for the iterates of Perturbed SGD-RR), i.e, Vi E [1,2] there exist measurable set C, constants as 0, a oo such that Vx Rd
xE,(x,²") = (x, + +
(23)
where
(x) = SEERS P(z,dr")((x")) - Ex(x) and the constant Yiaa
/146p2 Li
min
3nbman
12n
Proof. In order to prove that the geometric drift property is satisfied, we need to show that there exist
function & : Rd - [1, +00], measurable set C, and constants as > 0, d; < oo such that (23) holds. From Corollary D.3, we have that the function E R + [1, +oo with & x*) = | x x=11 +
satisfies along the iterates of Perturbed SGD-RR, that
E EI FR {-x ===x)) x}} I c1,1(1,x")) +
(24)
704 and 02 = 2 + 2 o² + BA Additionally, for the epoch-level iterates 2
where =1
Jy of Perturbed SGD-RRF the definition of A is
A
J (zz,dz')&1(z/11 - x2x,
[(11(T4+1) - x²) F {Ix =1}
(25)
From (24) and (25), we have that
E E1 x
f = 2}} s (11((()) + Fs Fs:
x*) - x²) F (xx at s
(1
AEi(xx²') S - -1 111--188 8 (") + 02
(26)
Let C1 = fo E Rd (((x,1)) (1-9)
2cg
}. We have that
1 1-
A6(1xx²) - -11 - c E (x) + Ic(2)e2 + (x)
El (1,xx²)
2
S 1 1 c1 3, (z, x + +c)
(27)
where at the last step we used the fact that 1cy (x) < 1 and e₁ E (0, 1). From (27) we conclude that
(1, 2²) satisfies the geometric drift property for the set C1 x E Rd Ex) I 202 and
(1-c1)
with constants 1-c1 C2.
For the Eg(x, 2(( by Jensen's 's inequality it holds that
E[VE(TR+1, Fi (xx = m3}
F [(((0411,0t) Fa fark 2)
< (e +
v V V02
Thus, there exist constants di Voi, d2 = Voz such that it holds
((E2(T4+1, Fx for = m3} S d1E2(1, + d2,
(28)
Since it holds that
Ax, [ ERd
P(x, dr' -((((,,"))__((xx=))) = E [[ x0) --2((0) F {or = x3
(29)
30

---

Published as a conference paper at ICLR 2026
have that
E( E2 x²) 1 (ax = =x} -
- -(1
-111---),((,,, +
7 (62 S -(1 di (1)8((1) x*) + d2
(30)
Let C2 fs E Rd x2(xx 1 (1-d1)
2dz } We have that
1 1-dp d1
AE2(1,))) v -(1 - di ( x + 1c((2) d2 + 1lgg (2)
2 (2(
E2(2, x²) + 1c(()) d2
2 where at the last step we used the fact that (x) < 1 and d1 E (0, 1). Hence, we conclude that (2 satisfies the geometric drift property for the C2 x E Rd E2(x) s [1-))] and with constants 022
1-d1 d2.
D.3 PROOF OF MINORIZATION PROPERTY
The next step in establishing ergodicity is to verify a minorization condition. Intuitively, this prop- erty guarantees that whenever the chain is in certain "small set" C, its one-step transition kernel dominates a fixed nontrivial distribution v, uniformly with probability > 0. In other words, starting from any æ E C. the algorithm has a baseline chance moving into any region of the state space according to v. This is the key ingredient that, together with Lyapunov-Foste drift condition, yields geometric ergodicity of the Markov chain. The following lemma formalizes this property for the iterates of Perturbed SGD-RR).
Lemma D.5 (Minorization property). Let Assumptions 2.1-2.3 hold. If Tmaxx then the iterates of Perturbed SCD--RR, satisfy the minorization condition: there exist constant > 0, probability measure v on (Rd, ((R")) and a set C Rd such that v(C) = 1, ((C)) O, and
P(x, A) 2 6 vc(x)v(A),
Vo E Rd, A E B(RR),
(31)
116 mar
where P(x, A) = Pr(xx+1 E A Fy = x)) and "/max = min
Inkmux
12n
Proof.
Consider the Lyapunov candidate E(x) || xx11 2 + 1 for some x" E X Its sublevel sets
C(r) for Rd E(x) r B(x Vr
r>
bounded, hence suitable for applying small/petite set arguments. At each epoch, the update of Perturbed SGD-RR, can be described by 1k+1 H(E,, wx) UR, where wh is uniform on Gn and UR No0,, Ela), independent of wi and æx. Thus, for any A E B(R*).
P(z. A)
Σ n! WER,
/ Ay: H(,,)), Σ) dy,
Since Φ(v; m. Σla) > 0 for all Rd, the kereel has strictly positive support everywhere.
Now fix ro and restrict C(ro). Define the reference measures for any B(Rd)
(AA = Leb(AnC(ro)) Leb(C(ro)) and Leb(A) = / inf
Σ dy: (x),,, 2) dy.
#EC(ra) n! WEG,
ie., the uniform probability distribution over C(ro). Clearly v(C(ro)) 1 and v(C(ro))) 9) 0. Finally, by continuity of & and compactness of C(ro). there exists 8 2 Leb(C(ra)) > such that P(x, A) > dv(A), Vr E C(ro), A C(ro).
I & C(ro) or A g C(ro), the hand side of (31) is zero and the inequality is trivially satisfied. Hence the minorization condition (31) holds.
31

---

Published as a conference paper at ICLR 2026
D.4 PROOF OF IRREDUCIBILITY, APERIODICITY AND HARRIS AND POSITIVE RECURRENCE Lemma D.6. The Markov chain (zx)k20 of Perturbed SGD RR1 is
1. / iiirreducible for some non-zero o-finite measure & on R& over the Borel o-algebra of Rd.
2. strongly aperiodic.
3. Harris and positive recurrent with an invariant measure.
Proof. We prove each of the three properties in turn.
Irreducibility.
From Lemma D.5, the Markov kernel of Perturbed SGD--RR; is
A)x r2! WEB,
Σ
/ oly: Hxx,w), Σla) dy. A E B(R)).
where <6:;m,, ΣTa) is Gaussian density with strictly positive support. Hence, for any measurable set A of positive Lebesgue measure, P(x, A) > 0. Taking v to be the Lebesgue measure, we conclude that the chain is vi-rrreducible.
Strong Aperiodicity. such that
By Lemma D.5, there exist 8 > 0, probability measure and set C' C Rd
A) 8 Lc()))(A),
Væ Rd A E B(RR).
Since C has positive Lebesgue measure and v(C) 1, v(C²) and given that the sets C(r) in the proof of the Lemma D.5 are small and of positive measure, we get that the Markov chain is strongly aperiodic.
Harris and Positive Recurrence. By Proposition D.I, the small set C' Lemma D.5 is also petite. Combined with the Foster-Lyapunov drift condition of Lemma D.4, the Geometric Ergodic Theorem (Theorem 15.0.1 in (Meyn & Tweedie, 2012)) guarantees that the chain is positive recurrent and admits an invariant probability measure.
Finally, from Theorem 9.1.8 of (Meyn & Tweedie, 2012), the existence of Lyapunov function unbounded off petite sets, satisfying AE < 0 together with @i-irreducibiltty, implics Harris recur- rence.
D.S PROOF OF EXISTENCE OF UNIQUE INVARIANT DISTRIBUTION AT EPOCH-LEVEL (THEOREM 4.3)
By verifying irreducibility, aperiodicity, and positive Harris recurrence (Meyn & Tweedie, 2012), establish unique invariant distribution Ty, geometric convergence in total variation to it, and concentration of scalar observables (admissible test functions) around x=".
Theorem D.7 (Restatement of Theorem 4.3). Under Assumptions 2.1-2.3, run Perturbed SGD- RR, with > S Tmax. Then (Ix)a20 admits a unique stationary distribution N. P2(Rd)) and additionally:
(i) [[[(I=)] - E [!(x)]] & c(1 - pp*
Ve [<()))) L2 + ||x1)),
(ii) [ [[(x)] - f(xx)) < L.VC
Ve: L, - Lipschitz functions,
C for some < oc, p E (0, C) @(MSE(SGD RR()) and 7max defined in Theorem 4.1
Proof. From Lemma D.6, we have that the underlying Markov Chain has an invariant probability measure. Since from Lemma D.4 the induced Markov Chain satisfies the geometric drift property, according to the Strong Ergodic Theorem (Meyn & Tweedie, 2012) we conclude that the measure finite and unique. From the invariant property of my, we have that for so W, the iterates satisfy also that (xx)sso Ky. From Corollary D.3, we have that for an arbitrary fixed x" the iterates of Perturbed SGD-RR, with step size 7 S satisfy for e₁ E (0, 1), 02 > 0 that
E [| - x² =11 + 11
< Ci (I
- x 112 + 1) + C2.
32

---

Published as a conference paper at ICLR 2026
Taking expectation with respect to the invariant measure ny and using the tower law of expectation, we get
Ex ||| x" 211 1 o (max()) < +00.
(32)
-
Combining the above inequality with the fact that ||ا+|| R by Assumption 2.1, we conclude that the invariant measure n, E P2((R)), where P2(RR) is the set of distributions supported in Rd with finite second moment.
We, next, proceed with proving the second statement of the Theorem. By assumption, we have that the test function satisfies Vr R that
|((x))
x Le(1 + ||=]]
< Le(1 + ||>* I| + 1|-- x²1l)
S Le(1 + R | x²1)
S (1 R)L || =²ll)
(33)
where we have used the triangle inequality and the fact that || S R. Applying Cauchy Schwwaz inequality, we can further upper bound [[{(x)]
|((x)] & v2(1 + R)LEV | x=ll s max (1, 22(1 + RR)Le) V((x,z*)
(34)
Letting = max (1, + L1) and E( x*) = eVE((, we have that
[((x)) & x0)
From Lemma D.4 we have that & (xr, E2 (xxx*)) satisfy the geometric drift property and since c 2 1 we have that E(x,x*)= = cE2(2, x²) satisfics also the geometric drift property. According to Theorem 16.0.1 in (Meyn & Tweedie, 2012) Perturbed SGD-RR, is E-uniformly ergodic and there exists (0, 1) and R E (0, +20) such that
ph Exwar, [e(x))] < R(1- p)* (( [
(35)
Letting R (xx, =*){ we have proven the inequality in the statement of the theorem. In order to show that the epoch-level iterates converge under the total variation distance it suffices consider only functions é Rd -> IR that are bounded by 1. In this case, there are constants ò (0, 1) and R E (0, +oo) independent of such that holds
sup pk p=(((x)== - Em [[(x)] | 51
<< R(1 - ((
implying according to the dual representation of Radon metric for bounded initial conditions (Wikipedia, Accessed: 2025-08-28) the geometric convergence under the total variation distance.
In order to prove the third statement of the theorem, apply linearity of expectation and the Lipschitz property of the test function and obtain
E Ezaz, [((x)) E(x*)
EzNz, [[(xx) {(x:))]
F Ex
[Te || x"
Applying Cauchy-Schward inequality and using inequality (32), we obtain that
Eowa, [[(x)] {(x*)) S Le Lev Ezwn, [||x - x |l| < L
where D x
max (7:)) according to (32).
33

---

Published as a conference paper at ICLR 2026
We conclude with the establishment of Law of Large Numbers (LLN) and the corresponding Centra limit Theorem (CLT) that describe the epoch-level iterates of Perturbed SGD-RR)-
D.6 PROOF OF LIMIT THEOREMS OF EPOCH-LEVEL ITERATES (THEOREM 4.3)
Theorem D.8 (Restatement of Theorem 4.4). Suppose Assumptions 2.1-2.3 hold and run Perturbed SRD-RR; y S Timax. (cf. Theorem 4.1).
Let e Idd + R be any test function such that [e(x)] S Le(1 + || and E [[(x)] < co.
Then for the epoch-level iterates, it holds that:
T-1
((x))) - ER [e(x))
(xx+) - F [((x))) 4 o (0),
(LLN)
(CLT)
where o2, (f) limp oc 7 1571 and S7 Σ (E(z1) Ex [(xx)))²
Proof. We show that the Markov Chain induced by the epoch-level iterates of Perturbed SGD-RR, is Harris positive recurrent, it has an invariant measure and satisfies E-uniform ergodicity, and hence by Theorem 17.0.1 in (Meyn & Tweedie, 2012) the stated Law of Large Numbers and Central Limit Theorem hold.
From Lemma D.6, we have that the Markov Chain Harris positive recurrent with an invariant measure. It suffices, thus, to show that the chain is E-unfform ergodic by proving that there exists potential function E-) such that the chain satisfies the geometric drift property of Meyn Tweedie
A+1 xx*) |F for any fixed x" X Accord- (2012) and [[( < E(z). Let (xr, x²) F {(
ing to Lemma D.4, Exx, x²) satisfies the geometrie drift property. Additionally, since l has linear growth it holds that
10(x))2 < L2 + ||
< L3 + ||>||| + | x x1)
s [2 + R || x 11)²
< L + R1)
x= 11)
(36)
From Cauchy-Schwand inequality, it holds that
+ || x='|
< v2v 1 + || - 12
= (1 | -,==0110² 11)
< 2(1 + ||-- x1 113)
(1 || x 11)
s 28(x,x")
(37)
Thus, combining (37) and (36), we obtain
1x)) <
2T2 + Rr)(x,=²))
(38)
Thus, Exx, x²) satisfies the geometric drift property and it holds that [0(x))² < E(x, x*) and hence the chain is E-uniform ergodic, completing the proof.
O
34

---

Published as a conference paper at ICLR 2026
E BOUNDING 4TH-ORDER MOMENTS OF ERROR DISTANCE
E.1 HIGHER ORDER VERSION OF PROPOSITION A.2 (EMMANOUILIDIS ET AL., 2024) BOUND KYRTOSIS FOR LIPSCHTIZ OPERATORS Proposition E.1. Let Assumption 2.3 hold. For any R and any reference point x² E R it
holds that
1 Σ iml
[F(()) -F())) & 128
LA
1 x=11" + 128
where o. = [B(*")!!
Proof. Define
A = F1(2) - Fi(x", A Σ A F(x) - F(x*), : f((x)- F(x"). j=1
Then Pi(x) F(x) = (A, A) + S, Using (a + 6)t - 8(a* + 6*) we have that
[[](()) F())1* 2 A| A 8||
Applying the same inequality once more to A A. we obtain
| A 811 + ||A||").
Averaging over and using Jensen's 's inequality for the convex map ++ ||u|]
[1((()) - F(z))) & 128 ( Σ (1 + Σ IK&.1T):
Moreover, ||&:|| [1 F(x*)| - [|[(x*) + |F (x=")) )ll hence by (a + 6) S 8(a² + b2)
and
Jensen,
n Σ ||&l|* & (6 Σ
|=(x))|* = 160
1=1
By Lipschitz continuity, we have
||ll|| [[F(z)) - F((x')) < Lell= --1
E +1 ( 10-
Combine the last two inequalities to obtain
[F(()) - ((())-F(>>1*+ < 128
4+) =| z"11" + 128
35

---

Published as a conference paper at ICLR 2026
E.2 A COMBINATORIAL TooL FOR HIGHER ORDEER MOMENTS IN SAMPLING WITHOUT REPLACEMENT
As part of our analysis, require bounds on the fourth moment of empirical averages when the underlying data are sampled withou! replacement. While this result is of independent combinatorial saaresstappparinng naturally in the study of randomization effects and variance reduction-it also plays a technical role in controlling higher-order error terms in our proofs. The following lemma provides a clean upper bound in terms of simple population statistics such as Σ and S4-
Lemma E.2 (Fourth moment of sample mean; upper bound). Let X1, Xa E Rd be fixed vectors, let
n
X1.
r r: X; X,
$ re
Σ
n1
n 1=1
and define the population sums
S4
|mill
U2: Σ 1n :2 -(rTr)))
Draw a size-k simple random sample without replacement, with indices (wi,..., wn) uniformly chosen among all k-subsets, and set X
Σi= Xw,- Then
1 k S4
9k(A: 1)
n(n - 1) (n²(tr 1222 sa)
E F| x11
(39)
K n
Proof. Write S = Σ Twn so that X
8 1s and 1 -111
Using the
Frobenius inner product (A, F) = tr(A B),
k
k
||||
=
Two
Twe Twu
T Two
was
Twu F²
1,u=1
t=1
1=
and hence
k
k
||S||
Σ
Σ t,u,s,0=1
(r rww rww (r Tww.)
Tw
Twu
f,u=1
Taking expectation and using the inclusion probabilities for simple random sampling without re- placement,
P(we = i)
P(wa = we =
n (),
we may group terms by the equality pattern among the positions (t. u, 8, v) (Hoeffiing/U-sttatiiiccc enumeration). Only three patterns survive:
(P1) Diagonal-diagonal: (t = u) and (s = v). This contributes
k:k 1) n 1) U2.
S4
+
(P2) Diagonal-off-diagnna (or vice versa): exactly one of the pairs (t, u) or (s) is diagonal and the other is off-diagonal. Counting gives coefficient 4(()) (R-1), leading to the contribution
4 (0) (7
To 6k(k 1) T2.
n(n - 1)
(P3) Off-diagonal-off-aaagona both pairs off- off-diagonal but correspond to the same unordered pair of distinct sampled units. This yields
2. A(R
2k(k: 1) n(n - 1)
U2
U2
U2.
nin
36

---

Published as a conference paper at ICLR 2026
Summing the three contributions obtain the exact identity
k S1
3k(k
6k(k 1)
E||X
x1I
U2
T2
(40)
b
n n(n
1)
n(n - 1)
Next, by
(r (r,Tr;)² ||: ||r| for all i j.
hence T2 S U2. Plugging this into (40) gives
1 E| [ - -II - k
33 + 6) -1)
1
l S1 +
9k(k 1)
S4
+
U2 =
U2
n(n- -
n(n 1)
Finally, observe the exact identity
U2 1r 17 11,
|m:11
|mll (t) S1.
because 1riel1 n tr Substituting this into the previous display yields the claimed bound (39).
Lemma E.3. Let Assumptions 2.1, 2.3 hold. If Perturbed SGD-RR, is run with step size 7 S
then it holds that
I 547"C1124 xx'1" 3456 1)o + 9727 n(nn + 1) (0?22 22ll [Fa]
E
1
where C 64m + 3n?(n + 1)L2 + 20
Proof. From the epoch-level update (8), it holds
||
74 4
(
3=0
277 4:3 Σ
Fag (x+) - Fo((o
+ +77 4,4
F F(x
1=0
+27 4:4 1
Assumption 2.3
277 4,3 LLaaz
Σ
-110k--0 - =i6| + 27774p0 4,1
Σ Fup (x2) - F(zb)
3=0
+27 4x4 =
where the last step have used the Lipschitz property of the operators F,, Vi expectation condition on the filtration Fs. we get
[n]. Taking
[|la - <<
277483
LA max
E
Σ
2H8|
FF
+ +27 4,4 E
Σ
F (x2) - F(22)
Fr
+ 27-444 1 [[(())) (41)
1=0
37

---

Published as a conference paper at ICLR 2026
Substituting in Lemma E.2 Xw, == Fug (xg) and le = n, get that
i-1 Σ j=0
( F(22) Fxx
S
S +
9i(i - 1)
(n²(tr si)]
1
n
nnn - 1)
St +
9in(i
(tr {
(42)
n
n 1
where S2 = -
(o))
F( 2 0 and f = E [F((x{) - F( We, next,
upper bound the terms S4 and ir St. From Lemma F.2, we have that
S4 S
128
10c x"1 128
(43)
n
i=1
Using Proposition A.2 (Emmanouilidis et al., 2024), we have that
1 n-1
F.(x2) - 2(2 s Al| - x²11 20
(44)
72
since each F; is Lipschitz. Substituting (43), (44) into (41), we obtain that
it E
E (o - -
Frx
1281
Σ
L - x² ==11** + 128i
n2
==1
9i(i 1)
+
(Allarx - -xx11 =2111 + 20
n 1
128i
128i
Aallark - xx ==111== +
0
n²
n
+ 18i(i n
1) A |atk x. ||
36i(i 1) (o5)
n
where we have let A = Ex= L for brevity. From inequality (45) and (41), thus, we obtain
||at - 2Q 221 Fi S
277413 L max
E
| - =b11
+27-A A|
128i
A4
I - x x=11"= +
128i
02
+
1) 18i(i
A A - x =1
36i(i 1)
n2 1 (o.)²)
n2
n 1
+277
2, max E <
Σ |1 - -=001') ==21|
+277 [ (1281
A4 + 18i(i - 1) 2) 1::ee - 1=110= 2 +
128i
367 - 1) (0?)²]
0 +
n²
n-
1
+27 [(
(46)
By summing over 0 <is n 1, we have that
n2(n - 1)2 n-1
Σ - =21
< 277 1
Σ E |||at =021
4
+2771
Ag 18n(n + 1) ||00 - =²11 + 1728 1)o
n
972n(n + 1) +-1
nnn 1)(2n 1) 3n - -1111111 (1(())))) (.
(o?)² + 277
(47)
n1
-1
30
38

---

Published as a conference paper at ICLR 2026
where we used the facts
(n - 1)
i(n i)
n(n + 1)
Σ
i(i - 1) n(n + 1)(n - -1),
Σ
2
n-
6
i=0
i=0
i=0
11(11 1)
n-1
nnn + 12n + 1)(m2 + 3n 1)
Σ
Σ
1d
30
k=0
We, thus, have that by rearranging the terms in (47) holds that
[ 277 AI4 n²(n 1)2
|| - x2 =211" Fia
4
1=0
I 27yd 64(n - 1) A4 + 18n(n + 1)A || x²11" ==11"+1287'00 17 1)0
ta* 972m(n + 1)
--1 (02) + 277 n( + 1)(2n + 1)(3? + +3 -1 F(x2))1
30
For S 3nt.mux
1
we have that (1 277 Iimax (m
2 1 2. and thus we obtain
EE[1 - [F s 547" ( (64(n-1) -1)
A4 + 18n(n + 1)A) 1 2"1""
n
+34567 1)0, 77 9722mn + +1)
(0.)²
n n
+54y n(n+ + 1)(2n + 1)(3n² + +3 1)
|F(2))|
30
2.3
n n²(n (n + 1)2 + +1)(")
S 547² (64A4 + 18n(n + 1)A
| =x"1'" x011
10
+3456-4(n - -1)0² + 7d 972n(n + 1)
(02)²
n²(n + 1)2 + +))L*)
I 547* (6nLL + 36n (n 1)
|| 2"11"
10
+34567 1)01 + 47- 9972n(n +
(02)²
n-1
where at the last step we used the fact that A1 = Ei= L < nL and A = 2Σ L
S
2nI2 max-
To establish Theorem 4.6, we will first prove upper bounds on the higher moments of the distance of the epoch-level iterates from the optimum, as well as the connection of the Lipschitz property of
the per-step operators F,, [n), and the Lipschitz constant of the epoch-level operator Gun. We start by providing the bound on the higher moments in the following Lemma.
Lemma E.4. Let Assumption 2.1-2.3 hold and À = 0. Then, the iterates of Perturbed SGD RR,
with stepsize 7 S Tmaz satisfy
E |1:2+1 - x..13 = o 2
F |1102-11 x1.11
O (rr)
v 1+6 1.
3/5
where Ymaz min
12n
3/5
Proof. We, first, provide the bound the third moment. From Holder's inequality, have that
[E [| xx. 22.111 -2_1 $
1
20.1p]]
<
39

---

Published as a conference paper at ICLR 2026
From Theorem 4.1, it holds that
E ||b - x² x²1[]]
= O (177),
and thus substituting we obtain
- 18-41-1-11
11]
< O (3
where we have used the fact that the fourth moment of the stochastic oracles bounded from As- sumption 2.4. Taking the limit of the Markov chain we obtain the bound on the third moment
Ellab+1 - x=118] S O (n2g²²
We, next, bound the fourth moment E ||02+1 - x2.11
We have that
|10bt1
2²11"
(1=0-1 k- x²" 2
(16)
|0 x² - (x
27 L
n-1
Σ
|| =011+=
27
|| |u&l?
2
- 7mp
ny
Letting T
| x² ->F((x))1l
T2 Em
| ==2ll and T3 ||u111 for brevity,
have
| 1²1 (1-mpp
Ti
27 TT
27 Ts) 1
2
+
(48)
mmp
H
From (9) in Lemma C.2, we can bound the term Ti by
Ti [(1 7m44)² + 2 1 1.maz]
| xx x²1 + 2ym
Substituting the bound of Ti into (48), we get
(1- (1p)) + + L2 max
Ih - x1 + 2 T2 + 2%
21
| k+1 - =x111 <
1- 7TU
11p1
(1-3mM))+ + Ln Limax for brevity and performing Young's inequality (7) with t C1s 1-gons
Letting C1
we obtain that
1100 A+I 1²1
<
(allax - 10012=--'11++
27
Ts
2
2 T2
n4
c1p2
412
4c1
2²1 2x I 112
alla 2²11"
T3
Tallax nu
T
12220
Taking expectation with respect to the permutation We and condition on the filtration F (history of x as well as using the fact that the noise UR is independent of the stochasticity of wk. we get
422 12 T3
2,211 L L
[102-11~0111 Fi
s
x 11
401 ny
x²11
E [23 Fi
(49)
Using Lemma E.3, we can bound the term
nnn + 1) (o.)2
E T Fi
< 547*C112x x= 34567 (n 1)0 972y
n n
where 64 2
3n + 1) +
n 1((==1)((2++))++ and thus substituting into (49) we
10
obtain that
E [1lok41 - F=
108-.L L < (al +
Liaa C) 11:00 - x01 +
472 T3 +
4c17
TulTallus x² 12
112
12
np
6912,, L
1)0 + 1944. L
n(n + 1) (o2)²
(n
(50)
112
112
40

---

Published as a conference paper at ICLR 2026
The next step involves taking expectation both sides with respect to the randomness of the Uk and will require a bound on the terms E
[3
E T3 Since U2 ~ N (0, z'u
d - we obtain
?
?2n
[[ = tr
03la = =??n00²
d
and
E || = ((d+2) (1 0² {nn)n((0))
2
as ||Ux||// (22 xa and F [(x2)) d² + 2d.
Thus, taking expectation both sides of (50) and substituting the bounds E have that
[Ts] E [E[3]]
E|11241 - x²11"|
$ c1
c) 1108
x=11+++ +
472
E
E
+
4er7 E T3
|| - x²1
122
22n2c1
nu
6912, 1
1944, 6 I m.max n(n + 1)
max - (n-1)
(02)²
1²
123
n2
+ (c1
108~ L L C) )1100 - =x111 +
1, 12201
+ 40 40 0 11:00 - x x=1
(0
112
6912, 6 max (n 1)o² 1944, 6 Lmaz 112
n(n (0
n 1 Taking expectation on both sides, using the tower law of expectation and the fact that from Theo- rem 4.1 IE |1a8 x²1
122
[|]]
O (n7²) of obtain
E E110241 =x"1=
(1 +
L
C
E El100 x x²11a]
4n
μ22c (o?)2
112
max 19447 L
4c o [
(02)²
ElIEEA x 112
112
11
6912, max
1944, Linax
ne + (02)
(n-1)o
112
n n-1
108 maz
4,012
01 +
E E|1000 - xx x²111]
μ²c (0?)²
122
1944,
(02)2+ + 4cr72n n O (m)²) (0?2²
mmax
112
H
6912,,
L
(1) + 1944, Lina
12(12 + 1)
max
(n
n 1 (02)2
112
122
108- 608-6T.4 c) E E|1108 - x² x²115] -==1|=++++ (1
nnn
cnnys
P + + (o2)
+
112
6912, I
(n 11
(51)
112
Selecting the stepsize 7 S min { |3nmmax'
1/5
{. we have that
3n
Sn 3/3
7 2,,1 1 max
7mpp 2
Ci
ymp +
mmp
and
108- 1
TA
Tmax C
ymp
108 6 I.
6 +
1
+
C
42
2
412
C<Gin' Loor
1
7m8 2
108- L
64n 5 T.A T.max
2
7mp
(52)
4
41

---

Published as a conference paper at ICLR 2026
Thus, substituting (52) in (51)
E [[102 -110211-x"11 - -011
(1
7M4
E |10
x²11*]
+ ( (o
6912, max
(n 1)0
112
Taking expectation with respect to the invariant measure πy and using Assumption 2.4 and the fact that o < +oc, we obtain
nn Exnny - - x²1l* < O (>*1)
and thus we have
E ExYTTY || - x²1 1 (
F PROOF OF BIAS FOR RR1 AND R METHODS (LEMMA 4.5) AND (THEOREM 4.6)
F.I JACOBIAN BOUND FOR ONE-PASS MAP
Lemma F.1 (Jacobian bound for one-pass map). Assume each component Rd + Rd is C' near xt with Fv(x' (==xx)))) lop I L; and let Lmaa max; L, For permutation w G,, define the
inner one-step map
(2)((x) = ------a (F)))( (x),
i .....,
its composition over one reshuffled pass
4 (x) = 4 Pwnn- .... p,,0(=),
and the epoch map Hxxww)
Σ Fun (z),),, (x,(), x, x1it)l(s,w) (2)
Then, at xz it holds that
IV=" S 1+ Σ (nlmaxx)
Consequently, the spectral radius of G,Z""" w) is at most 1 + (Taaax
Our first lemma aims to bound the maximum eigenvalue of the matrix VzH(w, x²) with respect to the known Lipschitz constants of the operators F,, Vi E [n].
Lemma E2. The maximum cigenvalue of the operator VzG(x*, w) is L = 1+ Σ (YLmaz)
Proof. Let twes (II, x -P, Owo (I, Own- 1 (K (),2)...) with 2.e0 (x,x) and obtain that
a (z). Define, also, the k-step operator (x, 2)
abt1 H( w) V. V.G(sh,w) = -- vol (26,
since G()g,w)= == DI-L Rag (r)).
The gradient of G(-,w) w) is computed by deriving first the partial derivatives of 00 (n) (zx 2) with respect and z. We prove by induction that
42

---

Published as a conference paper at ICLR 2026
Veill (x)
Vfu (->)"VPu, (ol (1,,2)) Vfw, (00 (1,,))
Vooll (x,2) =
n-
Σ
E(->)VE (x, (oL" (x,,
1=0
(x,)))
Viun-sti @w
For = 1, we have that 25 (2.2) = ow, (x, = =x-7Fu, (2) and it holds that
V "o (x, ->VFun (0)0 (I,)) I
thus the inductive hypothesis holds for 1. Assuming that holds for n 1, we, next. prove that it holds for n. We have that
V. T (x, 2) = Vedw (z, oi 1) (I, 2) (x) 2)
(-7)VFu ( (ola-1) (2, 2)) (->)"~'VF 'Vuu, (0l ... VFui (dl
(->)"VFun 00 (on-1)(x, z)) V.inni (x5, .... Vlw, (0 (Ir 2))
We, next, compute the gradient with respect to x and get
V (z 2) = Vapu(z, ol o(n (x, z)) + Vedw (z, (z, (1 2)) Ve0l5 (I (I, (53)
(I, 2)) V dw (z, d olm (Ir (z) VF (ol- (x) 2))
Using the fact that Ved ow and the inductive hypothesis for V Va
(I,
(m-1) 1) (I, 2), we obtain
V. v 2) = I- Fn, 12r 2)) Σ
(->) (z,=))) VFan (phn-1-0(x,, (z,, x))
1=0
n-2
= II+
V F (dfni V (8!n 2)) ..... ......unn (6hn (x,22)
1+]
(o."
(z, ... .... (x,
(z.2) (25x-21(I,22) ..... V. (1,
1=0
Thus, in order to compute V=G(w,x*) = I vol") (x", ,²",, we first compute V (m) x²).
Since x² is stationary point, it is a fixed point of the operator o6 (x*,x") = x,, Vj > 0. From chain rule, we have that
V (olm) (a*,z*)
= (x", x²) + V V20(") (1) x²)
(54)
Fu,, (x*) + n-1
= (->)"
III ==
(1>V-- Funn (x"))
(55)
j=1
n Σ1
( V (x*))
(56)
=
43

---

Published as a conference paper at ICLR 2026
In order find the maximum eigenvalue of the operator Vol multiplicative property of the operator norm to get
(n) (x* ,x)), (x*xx,),
we apply the sub-
8
12
Σ i=
y VF (x)),(())))
F FFWVFu (xn))())))
EI
7 F - (20))
1=1
(57)
( [mmax)
(58)
1=1
where L, is the maximum eigenvalue of V (x((")) and Lmaa is the maximum over all eigenvalues of (F, ), [n]- Since VG(w, x²) 1 VoLa (wh (xf, using the submultiplicative property of the operator norm we have that |V x²) ) llop I + 1VVl
(m) (2f. ll) and thus the maximum
eigenvalue of VG(w,>") is T.G Imma + Σ (7 mmaa)'
We, next, provide the theorem establishing that the combination of the two heuristics lead a refined bias of the order O (7*).
F.2 HIGHER-ORDER TERMS OF RR, BIAS (LEMMA 4.5)
Lemma E3 (Extended version of Lemma 4.5). Let and Assumptions 2.1-2.4 hold. If Perturbed SGD-RR, is run with < Ymaz, it holds that
biasPeettr SGD-RR; = limsup || [[ - ==1== +(+) k
$
Exy [r] = x. + A + (77))
where A V
H ) SRS C(x) (dx), C E[ wip
L 1
Dmazi M V HW
xw)
& I I
V
H(W
H (w, x²) 2 V (w, and the
1+6 1+6p2
3/5
maximum step size is Ymax
min
3n 3n
12n
8r L mea
Proof. From third order Taylor expansion of G around v., we have that
V,H(W, x²) v V² (x -0)==x))0 02 RI(a), Rd
(59)
where the reminder R1 (I:) satisfies supzER { ----0' | t00. Taking expectation with respect to the invariant distribution T., and using the fact that Ex,, (w,, x)] = 0, we get
(x)]
(ws, x)) (x - x 2002 + R((x)) = E v (W, x* (x x²) + 2
(60)
From Lemma E.4 and using Holder inequality and the fact that supnerd obtain
{
Mi (x))
)
+00, we
1 2 )(H(x"))) / Rd (I - 12)=2")02 = (2+3)
V. (w,x") Ex, fr x²1 2
(61)
44

---

Published as a conference paper at ICLR 2026
Taking the second order Taylor of G around x², we have that
0 x= = zō x= "V
x²) + (x( + +
(62)
with R2 the second order reminder satisfying SuP {
Ra(a)) a--2" x
+00. From the second order
moment of equation (62), the unbiasedness of the noise UR Vi, k E N, and Theorem E.4, we have that
Jec (x
= [1- (x-
+22 (()))((()) + (77))
Rearranging the terms, we get
M f (x xe)= [y(dx) 7 2/ Rd
+ -
(63)
where M VzH(w,,=*) 8
V VV
@ V (,H(w,,²),
We, next, show that the operator M is invertible for the selected step size by proving that it is
symmetric and positive definite. Let A, Vi E [d], be the eigenvalues V (w,x") with Hejd the corresponding eigenvectors. Note that I 2 V HH,, has eigenvalues (1 }x) 0 and thus for 7 < Amor (V. it symmetric positive definite on the same basis [eedg Hence we can factor the operator M as
M V, x @ DI+IC V.H 1 -V
& V,H(w,x")
VaH(w, s (T V
((w, V, (w,,a))
Thus, the vectors u, & jj Vi,j E [d] diagonalize M with eigenvalues Ma. A1 7dg) + 2 - (1,), ji,j E [d]. From Lemma F.2, we have that the maximum eigenvalue V (w: is I 1 Σ (IIm) and hence for y
have that TImax I and hence L
iG =: 1 n. Selecting the stepsize such that 7 < N it holds that pas >> O, Vi,j E [d] and thus M is positive definite and invertible. Thus, multiplying (63) with M -1 from the left, we get
| (z x=-x") (0-=z")²2T,,00 (dxx) = 7M-1 / ga
((z)T,(dx) + O (7²)
(64)
gd
Substituting (64) into (61) and rearranging the terms, we obtain
V HW, x: E.,, [r x²1 - ((, ( (dz)
-2 (,H((,,')) 22 ([W, x²) M J ((a)m
+ Ex, [r - x =
(65)
Letting 1 = - - v H (w, 2 H x² M JR C(x) (dr). we obtain
E [x] = x. + > 0 (7²)
(66)
F.3 PROOF OF BIAS REFINEMENTS OF RR 0RR (THEOREM 4.6) Theorem E.4 (Restatement of Theorem 4.6). Under the assumptions of Lemma 4.5, Algorithm output satisfies
Last-iterate version (line 9):
[[wx]] x01l
S c(1 p)* + 0(?2)
Averaged- -terrtt version (line 10):
E
1
Fm - x
< q/p + +073)
m
45

---

Published as a conference paper at ICLR 2026
where P E (0,1), oo (cf. Theorem 4.3).
Proof. From Lemma F.3, we have that the iterates Iy,k of Perturbed SGD-RR,will step size satisty
Ex,nory [x] I. + 7A O (>²*)
(67)
Similarly the iterates (#27,)) of SGD-RR with step size 2n satisfy
E [127]] x. 2 O (7²)
(68)
Thus, from (67), (68) compute the Richardson Romberg iterates
(Ez,,~,, [2x] Ex2,~12, - O (7²)
(69)
Consider the test function ((x) z. The function satisfies the assumptions both Theo- rem D.7, 4.4. Combining (69) with the rate that the iterates of the method tend the limiting
invariant distribution and the corresponding Central Limit Theorem from Theorems D.7, 4.4, we obtain
x*1
I c(1 p)* + 0(T)))
(70)
k
x²¹
c/e t + 0
(71)
m=
where p c0 (0, +00).
46

---

Published as a conference paper at ICLR 2026
G ON EXPERIMENTS
In this section, we provide additional details on the experimental setting used for the conducted experiments. We consider the setup of strongly monotone quadratic min-max problems
1 72 b
A x B 2/2 C,a² + of
min max fz 22) ZiERS agERd
x2.
The matrices A, are sampled by first sampling orthogonal matrix P and then sampling a diagonal matrix D; with elements in the diagonal uniformly sampled from the interval [/, L]. The selected parameters Jr, L correspond the strong monotonicity parameter in Assumption 2.2 and the Lips- chitz parameter of the underlying problem respectively. We acquire the matrices A, as the product A, PD, p We sample the matrices B2, C, similarly sampling the matrices A, with the only difference that the elements of the diagonal matrices D; lie in the interval [0. 0.1] and [M, L] respec- tively. The vectors a,, C₁ are follow the normal distribution I). In all experiments, we use = 100, d = 100, while we specify the values of μr, T. in each experiment independently as they differ.
We provide additional experiments on the effect of each heuristic in the convergence of the al- gorithm. More specifically, we compare the classical with-replacement SGDA algorithm, the RRI variant, the RRzvariant and the algorithm utilizing both heuristics. We run the experiment for multiple stepsizes 7 = {10-3 10 10-4, 10 $} and multiple condition numbers x {1, 5, 10).
Strongly Monotone
Strongly Monetone Problem
Stronply Monotone Frostem Ix
503a
5004
S0D4
12-1
50.3--4 481
5C2A-132
10-1
$224-8122
5006 RA1
10-
50D4 HR1
Rebatuve Erter
25
50 (112)
Raaiative Ertur
soib (812-4k))
Resative Eraor
Gobd (882-841
22-
12
13-1
Epochs (110))
Epochs Ix10')
Epachs (110))
Strongly (x
Strongly
20.3
5004
5638-311
$224-811
5224 RR1
503a
5000 8K2
S004 292
Deelaaiion Eltor
Relatie Error
22-1
10-
16-
10-
12"1
Epochs 0110")
Epechs tx10")
Epochs (x10")
Strongly Monotone Problem (x
Strongly Monctone Problem (xm
Strongly Monatone Problem 1xm100)
563a
$25A
SCEA
3030-141 AXL
23
10
5004 m1
Relkrroee Error
Relatve Erer
22-
12-1
12-1
Epochs (X10%)
Epachs kx10?)
Epochs (x10P)
Figure 5: Relative Error of the different variants of SGDA. Each row corresponds to a strongly monotone problem with condition number {1, 5, 10) and each row corresponds to different step size y = {10-3, 100 1,- est relative error comparison the other methods.
10 53. The combination of both heuristics RRy@RR achieves the small-
In Figure 5, we observe that all variants converge linearly to neighbourhood of the solution. Demonstrably, the variant leveraging both heuristics outperforms the other variants, reducing faster the relative error and validating the theoretical results established far.
We, next, provide an ablation study on the effect of the proposed heuristic in a variety of common algorithms used in VI and machine learning settings.
47

---

Published as a conference paper at ICLR 2026
Wasserstein GANs. We train Wasserstein GAN (WGAN) (Arjovsky et al., 2017) for learning the mean multivariate Gaussian and consider the effects of each heuristic in the training of the GAN. In a Wasserstein GAN, the optimization objective is two-player zero-sum game between the generator G()) and the discriminator D).), given by
inf sup ((..)) ((,, N(,, [w., + 0)]
(72)
The discriminator consists a linear classifier (xx; w) = (w.x) x)) of the input x E Rd, while the generator returns noisy estimation G(z; 0) = 2 ++0 of the learned parameter 0 E Rd after sampling noise vector (O I). In our experiments, the aim of the generator is to learn the mean p of the true Gaussian distribution with mean = [3, 4|" and covariance > To
Di23ncc to Optmal Paraaa Puran
10t
SEG
10
1.0
Distance to optimum
Distance to optimum
500 R1
SEG MM/L
20-4
S20 R33
500 R12-2
25n
10-2
30
ZOND
10-4
10-2
5OMD 1912
(X.09)
Rerations 2-10')
Rerations
Figure 6: Wasserstein GAN trained with different heuristics on top of a base algorithm. For all base algorithms, the generator trained with the combination of both heuristies RRy@RR, converges closer to the optimal parameters than the generator trained with any other algorithmic variant.
We examine the effect of the heuristics in a variety of different training algorithms and report the distance from the generator's optimal parameters for each experiment. Similar Emmanouilidis et al. (2024), we, first, consider the Stochastic Extragradient (SEG) method as the main algorithmic template for training and use each one of the 4 variants (SEG, SEG-RR1, SEG-RRg, SEG-RRg + RR1) to train GAN. We use the same constant step size for the generator and discriminator as in Emmanouilidis et al. (2024) and double the step size of the variants that implement Richardson- Romberg extrapolation. Figure 6 shows that the generator trained with SEG-RR2 @ RR, is able to converge closer to the optimal parameters than the generator trained with any other variant and thus the synergy of both heuristics (RR2 e RR₁) is beneficial in training.
Following Daskalakis et al. (2017), we train a WGAN with the use of Optimistic Mirror Descent (OMD). Aiming to the effect of each heuristic even for this algorithm, use the classical OMD method, the RR₁ variant, the random reshuffling (RR2) and the combination of both (RR₂ @ RR)). We let the step size of the generator and the discriminator be na 0.02, 7D 0.01 respectively. As shown in Figure 6, the RR₂ @ RR, outperforms all other variants, indicating that the advantages of this heuristic remain apparent even for the OMD algorithm.
Lastly, we test lightweight second order methods common in the literature of VIs (Mescheder et al., 2017; Loizou et al., 2020). More specifically, Stochastic Consensus Optimization (SCO) is al- gorithm that can be seen as combination of the SGDA algorithm and the Stochastic Hamiltonian (SHMD) method (Loizou al., 2020), where a regularizer articulates the contribution of SHMD that is being introduced in the update rule. Given that the SCO method is related to the SGDA but requires Jacobian vector products, thus being a lightweight second order method, we have run experiments to examine whether the RR₃ ① RR; provides benefits in higher-order methods. Accord- ing to Figure 6, the generator trained with the heuristic RR2 ① RR, converges closer to the optimal parameters than the generator trained with plain SCO or any other variant.
On Single-Run Experiments & Variance of Observed Behaviour. For completeness, we report the variability of our experimental results over single runs, establishing more refined description of the effect each heuristic empirically. More specifically, in Figure 7 we plot the mean and standard deviations for each of the 4 variants over 5 runs. As shown in Figure 7, the RR₂ O RR, variant outperforms all other heuristics even in single trials.
Wall-clock time Comparison of the different heuristics. We, next, compare the wall-clock time of the different heuristics. All 4 variants have the same per iteration cost in terms of gradient eval- uations, since the only difference between the classical SGD algorithm and the RR, variant is the
48

---

Published as a conference paper at ICLR 2026
Strongly Monotone Problem (K 1)
10"
SGDA-RRI
SGDA-RR2
SGDA-IAR2+RAL)
Re Relative Emror
10-3
10
10
Epochs (x10²)
Figure 7: Mean and standard deviation of each heuristic over 5 trials. The combination of both heuristics accheees the smallest relative enror in comparison to the other methods.
way that the mini-batch gradients are sampled, while in RR2 and RR2@RR, two chains of the Richardson extrapolation can be run in parallel.
Interestingly, SGDA with random reshuffling typically runs faster in wall-clock time than SGDA with with-replacement sampling. The following factors explain why random reshufiling can be run faster, as observed in our experiments:
RR1 performs only random operation per epoch. With with-replacemmnt sampling. each iteration requires a random draw is ~ Uniiorm(1..........
3R, calls the PRNG once per epoch (through randperm(n) or equivalent), after which all iterations sequential. This eliminates thousands of PRNG calls and reduces interpreter overhead.
We have reproduced the same experiment as in Figure and have reported the wall-clock time needed for each method. According to table 2. the RR2 + RRj heuristic runs in half the time required for the plain SGDA variant and comparable time with respect to random reshuffling. Hence, thanks to parallelization one obtain the benefits from the synergy of the two heuristics without the need of a higher wall-clock time.
Table 2: Wall-clock time comparison of SGDA variants.
Method
Time (sec)
SGDA
107.59
SGDA-RR
50.92
SGDA-RR2
107.59
SGDA-RR, + RRj
51.94
49


Reviews:
[{"review_id": "bD8x8SoNh8", "reviewer": "Reviewer_WTpW", "summary": "This work studies variant of stochastic gradient methods with fixed step-size for weak quasi-strongly monotone variational inequality problems, which generalizes optimization and min-max problems, showing that combining Random Reshuffling with Richardson–Romberg extrapolation removes the usual bias of $\\gamma^{3/2}$ (for SGD-RR2), $\\gamma^{2}$ (for SGD-RR1) and achieves a cubic-order asymptotic refinement $\\gamma^3$. The paper provides the first theoretical analysis of this combination in quasi-strong monotone VIs, using Markov chain and spectral tensor tools. Experiments confirm the theoretical predictions, demonstrating significant empirical speedups.", "questions": "Please see the weaknesses section. Additionally:\n1. Regarding Thm 3.6, the results are a little bit counterintuitive. Under quasi-strong monotonicity, one can expect a linear rate to a neighborhood of the solution. It is also known that in VIs, average iterates usually converge faster (at least in the monotone deterministic case) than the last iterate, while in the second result of Thm 3.6 the rate for the average iterate is sublinear. Can you please elaborate on this?\n 2. The only assumption you make is on the existence of the solution. Does Theorem 3.6 hold for all solutions $x* \\in X*$ , or only for the projection of $ x_k $ onto $ X^* $?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Overall, paper is well written, have strong theoretical contribution which can be applied in practice. \n1. Authors proposed to use both RR1 and RR2 simultaneously to reduce the bias of SGD for solving finite-sum VIs.\n2. The main result in Theorem 3.6 provides tightens the neighborhood of convergence from previously known $\\gamma^2$ to $\\gamma^3$.\n3. The proof technique seems to be new and not trivial. \n4. Their experimental results support the theory.", "weaknesses": "I did not find any major weakness in the paper and list minor corrections below:\n\n1. Please define PreProcess. It is possible to guess from the next paragraph that it is a Gaussian perturbation, but it would be nice to define it similarly to StochOracle in line 137.\n\n2. The assumptions on the operator of the VI problem are relatively strong. In particular, all $F_i$ need to be Lipschitz. Would the analysis be more involved under the assumption of $F$ being Lipschitz? Would it be possible to relax the quasi-strong monotonicity assumption?\n\n3. I think Theorem 3.3 (lines 346–352) should be explained in more detail, especially the functions $l(x)$ and $L_l(x)$ should be defined for easier readability.\n\n4. I believe it would also be interesting to see plots for a single run of SGD, SGD-RR1, SGD-RR2, and SGD-RR1+RR2, since this is more common in practice. Do you observe that SGD-RR1+RR2 outperforms vanilla, RR1, and RR2 for a single run?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qoqzaBj6OD", "reviewer": "Reviewer_3rkQ", "summary": "The article explores random permutations and Richardson–Romberg extrapolation for constant step-size SGD/SGDA in variational inequality problems (VIPs). It states that the combination of these two techniques reduces the bias order of the stationary distribution to $O(\\gamma^3)$ for quasi-strongly monotone operators. The authors develop a detailed theoretical framework using the analysis of Markov chains at the epoch level, Foster–Lyapunov drift arguments, and tensor–spectral expansions, and confirm their conclusions empirically on synthetic quadratic problems.", "questions": "See my weeknesses. I would be grateful if the authors would resolve my doubts about the additional factor $\\frac{L}{\\mu}$ in the estimate, and about the small values of stepsize. Can the authors provide an experimental confirmation of the revealed effect on realistic ML tasks?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. **Clear and well-motivated question.**\nIt is important to understand how the combination of random reshuffling with Richardson–Romberg extrapolation affects bias. This is a theoretically complex issue.\n\n2. **Technically sophisticated analysis.**\nThe paper presents a nontrivial combination of tools $-$ epoch-level Markov modeling, spectral tensor expansions, and moment control. The analytical pipeline is rigorous and, to my knowledge, novel for this combination of RR and extrapolation.\n\n3. **Theoretical results.**\nLemma 3.5 and Theorem 3.6 provide explicit asymptotic bias rates up to third-order terms in $\\gamma$, which is a novel rate.", "weaknesses": "1. **Comparison to work [1] and convergence rate concerns.**\nWhile the paper claims an $O(\\gamma^3)$ bias improvement, it is not immediately clear that this is conceptually distinct from previous works, e.g. Stochastic Extra-Gradient with Random Reshuffling (SEG-RR) by Emmanouilidis et al. [1].\nTheir work already provides higher-order dependence on $\\gamma$ in convergence bounds for affine and strongly monotone VIPs. The difference lies in the metric (stationary bias vs. convergence rate) and in the algorithmic structure (plain SGD + extrapolation vs. SEG). The convergence estimate itself also raises a concern. In the theoretical analysis of convergence, the factor $\\gamma^2$ appears in the bias term instead of the classical $\\gamma$, but an additional factor $\\frac{L}{\\mu}$ also emerges. This is a problem, since such an extra factor indicates a deterioration of the bound, especially in real-world problems where $\\frac{L}{\\mu}$ is large. Moreover, it seems that the improvement in the estimate of bias is achieved precisely due to such a deterioration in the condition number.\n2. **Unclear practical relevance of $O(\\gamma^3)$.**\nThe bias order is $O(\\gamma^3)$; however, the allowed stepsize $\\gamma_{\\text{max}}$ vary depending on task parameters $n, L_{\\text{max}}, \\mu$. It is possible that the “bias improvement” manifests only at unrealistically small $\\gamma$. The paper should provide numerical examples of typical acceptable values of $\\gamma$ to clarify whether this is only a theoretical improvement or a practically relevant regime.\n\n3. **Dependence on Gaussian pre-processing.**\nThe algorithm’s PREPROCESS step introduces Gaussian smoothing, but the paper does not convincingly argue for its necessity in practice.\n\n4. **Limited empirical validation.**\nExperiments are purely synthetic (quadratic min–max tasks). There are no tests on realistic ML or game-theoretic problems (e.g., GANs or robust regression).\nThis makes it unclear whether the $O(\\gamma^3)$ bias scaling provides any tangible benefit beyond theoretical neatness.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "iJAG7blK82", "reviewer": "Reviewer_Vvn1", "summary": "The paper studies constant–stepsize stochastic methods for finite-sum VIPs and shows that composing random reshuffling (RR1) with Richardson–Romberg extrapolation (RR2) cancels all lower-order bias terms, yielding $O(\\gamma^3)$ asymptotic bias for quasi-strongly monotone smooth VIPs. The analysis blends an epoch-level Markov chain view of reshuffling with a refined bias expansion and a spectral/tensor argument for the extrapolation step. Experiments on saddle-point games support the theory. I am satisfied with this manuscript and believe it will make a valuable contribution to the community.", "questions": "1. **About (c_2) (Line 1494).** You state $(c_2 \\in (0,1))$. I see (c_2) defined around Line 1475, but I could not follow how the admissible range is derived from that definition. Could you clarify the step that enforces (c_2\\in (0,1))?\n2. **Practical problem classes/benchmarks.** Beyond the synthetic/quadratic games, are there data-level benchmarks that exactly satisfy your assumptions? Even if the assumptions are only approximately met, are there realistic tasks (or public benchmarks) where SGDA–RR2⊕RR1 shows clear empirical gains over baselines?\n\n### What would be change my opinion?\nAt the very least, if the revision includes a well-developed Limitations section or any equavalent clarification, I would be inclined to vote for acceptance of this paper. If additional comparison baselines are included or the method’s effectiveness is demonstrated on other ML problems/benchmarks, I will argue even more strongly for the strength of this paper’s contributions.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Modeling one reshuffled pass as a time-homogeneous Markov kernel (with Harris recurrence and LLN/CLT) gives clean asymptotic control of per-epoch iterates, a technically elegant way to handle the permutation-induced bias. This framing clarifies why RR1 can sharpen bias/MSE and provides a principled scaffold for the outer-loop extrapolation analysis. \n- The drift/minorization arguments (via Lyapunov–Foster) yield uniqueness of an invariant distribution and geometric convergence, which in turn legitimizes using stationary expansions for bias cancellation in the RR2 step. The explicit LLN/CLT for scalar observables further explains the empirical stability of epoch-averaged statistics. \n- The paper isolates that RR1 alone yields a linear-plus-cubic bias expansion and then shows RR2 cancels the linear term, leaving $O(\\gamma^3)$. This surpasses known rates for either heuristic in isolation and is obtained under quasi-strong monotonicity. \n- The debiasing proof uses a spectral view of the full-pass operator (akin to multi-step extragradient) to bound eigenstructure under reshuffling noise, enabling a refined Taylor-type expansion where all sub-cubic terms vanish for the composed algorithm. \n- Algorithm 1 aligns with practical pipelines (RR1 inside, RR2 outside), and the plots consistently show the RR2⊕RR1 curve settling to a smaller bias neighborhood than baselines across condition numbers. This match between theory and experiment strengthens the paper’s practical message.", "weaknesses": "* **[Major] Scope of Assumption 2.2.** The λ-weak μ-quasi strong monotonicity condition, while weaker than strong monotonicity, still excludes many ML-relevant VI problems, e.g., GANs, adversarial training, and multi-agent RL (typically modeled as nonconvex–nonconcave games). However, these problems are mentioned in Introduction section as VIPs. To avoid misleading readers, I recommend softening the tone in the Introduction and adding a dedicated **Limitations** section that clearly delineates where the guarantees apply/do not apply.\n  - Jin, Chi, Praneeth Netrapalli, and Michael Jordan. \"What is local optimality in nonconvex-nonconcave minimax optimization?.\" International conference on machine learning. PMLR, 2020. \n  - Han, Andi, et al. \"Nonconvex-nonconcave min-max optimization on Riemannian manifolds.\" Transactions on Machine Learning Research (2023). \n  - Kim, Beomsu, and Junghoon Seo. \"Semi-Implicit Hybrid Gradient Methods with Application to Adversarial Robustness.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2022. \n  - Bukharin, Alexander, et al. \"Robust multi-agent reinforcement learning via adversarial regularization: Theoretical foundation and stable algorithms.\" Advances in neural information processing systems 36 (2023): 68121-68133.\n\n* **[Major] Wall-clock comparisons missing.** Figures 1 and 3 plot error vs. epochs, but SGDA–RR2⊕RR1 likely incurs higher per-iteration computation than SGDA/SGDA-RR1/SGDA-RR2. To substantiate practical gains, please add plots against **wall-clock time** (or FLOPs), ideally with breakdowns for inner/outer steps, so readers can judge time-to-accuracy fairly. \n* **[Minor] Baseline coverage.** Given Loizou et al. (2021) and related work, two additional methods seem relevant and it would be good for them to be included as baselines:\n  - [SCO] stochastic consensus optimization (SCO): Mescheder, Lars, Sebastian Nowozin, and Andreas Geiger. \"The numerics of gans.\" Advances in neural information processing systems 30 (2017). \n  - [SHGD] stochastic Hamiltonian gradient descent: Loizou, Nicolas, et al. \"Stochastic hamiltonian gradient methods for smooth games.\" International Conference on Machine Learning. PMLR, 2020.\n  - Loizou, Nicolas, et al. \"Stochastic gradient descent-ascent and consensus optimization for smooth games: Convergence analysis under expected co-coercivity.\" Advances in Neural Information Processing Systems 34 (2021): 19095-19108.\n  A short ablation with these would contextualize the gains of RR2⊕RR1.\n* Line 1587: (C^{\\circ}) appears undefined in context. It likely should be (C^{c}) (the complement)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)