arXiv CSopen access

Single-Loop Stochastic Projected Damped Extragradient Methods for Stochastic Nonconvex--(Strongly) Concave Minimax Optimization

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

Single-Loop Stochastic Projected Damped Extragradient Methods for Stochastic Nonconvex–(Strongly) Concave Minimax Optimization Huiling Zhang1† Minhao Zhang2† Zi Xu2∗

arXiv:2609.21747v1 [math.OC] 18 Sep 2026

LSEC, ICMSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China [email protected]

Department of Mathematics, College of Sciences, Shanghai University, Shanghai 200444, China [email protected]; [email protected]

Abstract We develop single-loop stochastic projected damped extragradient methods for stochastic nonconvex–(strongly) concave minimax optimization, with complexity guarantees for both game stationarity (GS) and optimization stationarity (OS). Our approach combines a stochastic projected damped extragradient (SPDE) method with a recursive variance-reduced variant, VR-SPDE, both of which retain a single-loop structure. Under an unbiased stochastic gradient oracle with uniformly bounded variance, SPDE finds an ε-game-stationary point with stochastic first-order oracle (SFO) complexities of O(κε−4 ) and O(ε−5 ) in the nonconvex–strongly concave and nonconvex–concave settings, respectively, where κ = L/µ. Under an additional mean-square Lipschitz condition on the stochastic gradients, VR-SPDE improves these GS complexities to O(κ3/2 ε−3 ) and O(ε−9/2 ), respectively. For an ε-optimization-stationary point, SPDE achieves SFO complexities of O(κε−4 ) and O(ε−6 ), while VR-SPDE achieves O(κ3/2 ε−3 ) and O(ε−6 ), in the two settings, respectively. These OS guarantees match the best-known bounds achieved by multi-loop methods while preserving a single-loop implementation. To the best of our knowledge, our results provide the best-known SFO complexity guarantees among single-loop stochastic first-order methods for the respective stationarity criteria and problem classes.

Keywords: stochastic nonconvex–(strongly) concave minimax optimization; stochastic projected damped extragradient method; variance reduction; single-loop algorithms.

Introduction

We consider the stochastic minimax problem min max F (x, y) := Eω∼D [f (x, y; ω)], x∈X y∈Y

(1)

where X ⊆ Rn is closed and convex, Y ⊆ Rp is compact and convex, and F is smooth. The objective may be nonconvex in x and is concave or strongly concave in y. The random variable ω follows an unknown distribution D, and E denotes expectation. Problems of this form arise in distributed nonconvex optimization [4], wireless systems [1], and statistical learning [12]. Their growing scope of applications motivates stochastic first-order methods that combine a simple iteration structure with strong oracle complexity guarantees. † ∗

The first two authors contributed equally to this paper. Corresponding author.

Two stationarity criteria are central to nonconvex minimax optimization: game stationarity (GS) and optimization stationarity (OS). Game stationarity concerns first-order conditions for the two players at a primal–dual pair, whereas optimization stationarity concerns the minimization of the primal value function ϕ(x) := max F (x, y). y∈Y

These criteria capture different aspects of approximate solutions, and their complexity guarantees must be distinguished. This distinction is particularly relevant in the nonconvex–concave setting, where ϕ may be nonsmooth. We state the precise GS and OS criteria in Section 2 and analyze both criteria for the methods developed in this paper. For stochastic nonconvex–strongly concave problems, representative multi-loop methods obtain strong stochastic first-order oracle (SFO) complexity guarantees by approximately solving an inner maximization problem or a regularized saddle-point subproblem. Luo et al. [8] proposed SREDA, which combines nested ascent with SPIDER recursion and requires O(κ3 ε−3 ) stochastic gradient evaluations, where κ = L/µ, L is the gradient Lipschitz constant, and µ is the strong concavity parameter. Zhang et al. [15] developed SAPD+, which combines an inexact proximalpoint outer scheme with an accelerated primal–dual inner solver. Their method requires O(Lκε−4 ) oracle calls, while its variance-reduced variant requires O(Lκ2 ε−3 ) calls. Among single-loop methods, Lin et al. [7] analyzed two-timescale stochastic gradient descent–ascent (SGDA) and established an O(κ3 ε−4 ) stochastic gradient complexity for value-function stationarity. Huang et al. [5] introduced Acc-MDA with a STORM-type estimator and obtained complexities of e 9/2 ε−3 ) without a large batch and O(κ e 5/2 ε−3 ) with a batch size of order κ4 . O(κ For stochastic nonconvex–concave problems, the possible nonsmoothness of the value function presents an additional challenge. Multi-loop methods often address this difficulty through proximal regularization or regularized saddle-point subproblems. Zhang et al. [15] applied SAPD+ to this setting and established an O(ε−6 ) oracle complexity. Among single-loop methods, Lin et al. [7] established an O(ε−8 ) stochastic gradient complexity for SGDA under value-function stationarity. For game stationarity, Zhang and Xu [14] proposed FORMDA, which combines vanishing dual regularization with recursive momentum and achieves an iteration complexity of e −13/2 ). O(ε These results highlight the challenge of obtaining strong complexity guarantees within a single-loop implementation. Multi-loop methods achieve sharper dependence on the target accuracy in several settings, but their outer iterations may require inner maximization, proximal, or regularized saddle-point subproblems to be solved to prescribed accuracies. The resulting inner-loop lengths and parameter schedules can depend on the target accuracy. This motivates the following question: Can single-loop stochastic first-order methods achieve improved game-stationarity guarantees and match the best-known multi-loop optimization-stationarity guarantees for both nonconvex–strongly concave and nonconvex–concave problems? We address this question through a stochastic projected damped extragradient method and its recursive variance-reduced variant. Both methods retain a single-loop structure, and their analysis provides complexity guarantees for both GS and OS. Contributions. Our contributions concern the algorithmic design and the complexity guarantees under the two stationarity criteria. (i) Single-loop stochastic damped extragradient methods. We propose a stochastic projected damped extragradient (SPDE) method under an unbiased stochastic gradient oracle with uniformly bounded variance. We further incorporate a recursive variancereduction mechanism to obtain VR-SPDE under an additional mean-square Lipschitz condition on the stochastic gradients. Both algorithms use a single-loop structure without nested iterative solvers for maximization, proximal, or saddle-point subproblems. 2

Table 1: Representative complexity guarantees for stochastic nonconvex minimax optimization. Bounds are total SFO complexities unless marked by †, which denotes an iteration complexity. The dependence on L is suppressed in this comparison. Method

Setting

Loop VR Criterion

Complexity

SGDA [7] Acc-MDA [5] SREDA [8] SAPD+ [15] VR-SAPD+ [15] SPDE (this paper) VR-SPDE (this paper)

NC–SC NC–SC NC–SC NC–SC NC–SC NC–SC NC–SC

Single Single Multi Multi Multi Single Single

No Yes Yes No Yes No Yes

OS OS OS OS OS GS/OS GS/OS

O(κ3 ε−4 ) e 9/2 ε−3 ) O(κ O(κ3 ε−3 ) O(κε−4 ) O(κ2 ε−3 ) O(κε−4 ) O(κ3/2 ε−3 )

SGDA [7] SAPD+ [15] SPDE (this paper) VR-SPDE (this paper)

NC–C NC–C NC–C NC–C

Single Multi Single Single

No No No Yes

OS OS OS OS

O(ε−8 ) O(ε−6 ) O(ε−6 ) O(ε−6 )

FORMDA [14] SPDE (this paper) VR-SPDE (this paper)

NC–C NC–C NC–C

Single Yes Single No Single Yes

GS GS GS

e −13/2 )† O(ε O(ε−5 ) O(ε−9/2 )

(ii) Complexity guarantees for game stationarity. SPDE finds an ε-game-stationary point with total SFO complexities of O(κε−4 ) and O(ε−5 ) in the nonconvex–strongly concave and nonconvex–concave settings, respectively. The corresponding iteration complexities √ are O( κε−2 ) and O(ε−5/2 ). Under the additional mean-square Lipschitz condition, VRSPDE improves the GS oracle complexities to O(κ3/2 ε−3 ) and O(ε−9/2 ), respectively, while preserving the single-loop structure. (iii) Complexity guarantees for optimization stationarity. For an ε-optimizationstationary point, SPDE requires O(κε−4 ) and O(ε−6 ) SFO calls in the nonconvex–strongly concave and nonconvex–concave settings, respectively. VR-SPDE achieves the corresponding complexities of O(κ3/2 ε−3 ) and O(ε−6 ). These OS guarantees match the best-known multi-loop dependence on the target accuracy in both settings while retaining a single-loop implementation. To the best of our knowledge, these results provide the best-known total SFO complexity guarantees among single-loop stochastic first-order methods for the respective stationarity criteria and problem classes. In particular, the same single-loop algorithms support both GS and OS guarantees. In the strongly concave setting, the two criteria have the same stated SFO bounds for each method. In the merely concave setting, the GS bounds are O(ε−5 ) for SPDE and O(ε−9/2 ) for VR-SPDE, whereas both methods attain an OS bound of O(ε−6 ). Table 1 summarizes representative results and identifies the stationarity criterion associated with each bound. In Table 1, NC–SC and NC–C denote the nonconvex–strongly concave and nonconvex–concave settings, respectively. “VR” indicates whether a method uses recursive variance reduction, and “GS” and “OS” denote game stationarity and optimization stationarity. The Acc-MDA entry reports the bound without a large batch. Each result is stated under the assumptions of the corresponding reference; comparisons must account for both the stationarity criterion and the stochastic oracle assumptions. Organization and notation. Section 2 presents the preliminaries and defines the GS and OS criteria. Section 3 introduces SPDE and its stochastic oracle model, develops a unified convergence analysis, and establishes the GS complexity bounds for the nonconvex–strongly concave and nonconvex–concave settings. Section 4 introduces the paired oracle, recursive estimator, and VR-SPDE algorithm, and derives the corresponding variance-reduced GS complexity bounds

with all oracle calls counted. Section 5 establishes the OS guarantees for the original primal value function. Section 6 concludes the paper. Throughout, ∥·∥ denotes the Euclidean norm, ⟨·, ·⟩ denotes the corresponding inner product, and ΠC and NC denote the Euclidean projection onto and the normal cone of a closed convex set C, respectively. The notation O(·) suppresses constants independent of the target accuracy ε e and of any parameters whose dependence is displayed explicitly. The notation O(·) additionally suppresses logarithmic factors.

Preliminaries and stationarity criteria

We first state the standing assumptions and then introduce the stationarity criteria used in our complexity analysis. Throughout, product spaces are equipped with the Euclidean norm. Assumption 1. (i) The set X ⊆ Rn is nonempty, closed, and convex. The set Y ⊆ Rp is nonempty, compact, and convex, with DY := max∥y∥ > 0. y∈Y

(ii) For almost every ω, the sample loss f (·, ·; ω) is differentiable on a neighborhood of X × Y . (iii) The function F is continuously differentiable on a neighborhood of X × Y , and there exists L > 0 such that, for all (x, y), (x′ , y ′ ) ∈ X × Y , ∥∇F (x′ , y ′ ) − ∇F (x, y)∥ ≤ L∥(x′ − x, y ′ − y)∥.

(2)

(iv) For every x ∈ X, the function F (x, ·) is concave on Y . The primal value function satisfies ϕ(x) := max F (x, y), y∈Y

x ∈ X,

ϕinf := inf ϕ(x) > −∞. x∈X

(3)

Compactness of Y and continuity of F ensure that the maximum in (3) is attained for every x ∈ X. The quantity DY bounds the norm of points in Y ; it is not the diameter of Y . In the nonconvex–strongly concave setting, we additionally assume that F (x, ·) is µ-strongly concave on Y , uniformly in x ∈ X, for some µ > 0. Equivalently, for every x ∈ X and y, y ′ ∈ Y , F (x, y ′ ) ≤ F (x, y) + ⟨∇y F (x, y), y ′ − y⟩ −

µ ′ ∥y − y∥2 . 2

We write κ := L/µ in this setting. The stochastic oracle assumptions are stated separately with the corresponding algorithms. Projection and normal cone.

For a nonempty closed convex set C ⊆ Rd , define (

{η ∈ Rd : ⟨η, a − u⟩ ≤ 0 for all a ∈ C}, u ∈ C, ∅, u∈ / C. (4) The projection is uniquely defined and satisfies 1 ΠC (w) := argmin ∥u − w∥2 , u∈C 2

u = ΠC (w)

⇐⇒

NC (u) :=

w − u ∈ NC (u),

∥ΠC (w) − ΠC (w′ )∥ ≤ ∥w − w′ ∥.

For a nonempty set A, we use dist(v, A) := inf a∈A ∥v − a∥.

(5)

Game stationarity. For a feasible pair (x, y) ∈ X × Y , define the game-stationarity residual of the original minimax problem by R(x, y)2 := dist2 0, ∇x F (x, y) + NX (x) + dist2 0, −∇y F (x, y) + NY (y) .

(6)

Thus R(x, y) = 0 if and only if 0 ∈ ∇x F (x, y) + NX (x),

0 ∈ −∇y F (x, y) + NY (y).

These inclusions express the first-order conditions for the minimization and maximization variables, respectively. By concavity of F (x, ·), the second inclusion is equivalent to y ∈ arg maxv∈Y F (x, v). We call a random feasible pair (xout , yout ) ε-game-stationary in expectation if E R(xout , yout )2 ≤ ε2 .

(7)

By the Cauchy–Schwarz inequality, this condition also implies

E R(xout , yout ) ≤ ε. Optimization stationarity. Optimization stationarity concerns the constrained minimization of the primal value function ϕ. To incorporate the constraint explicitly, define (

ϕ̄(x) :=

ϕ(x), x ∈ X, +∞, x ∈ / X.

For each y ∈ Y , Assumption 1 implies that F (·, y) is L-weakly convex on X; that is, x 7−→ F (x, y) +

L ∥x∥2 2

is convex on X. Taking the pointwise maximum over y ∈ Y shows that ϕ is also L-weakly convex on X. Moreover, continuity of F and compactness of Y imply that ϕ is continuous relative to X. Consequently, ϕ̄ is a proper, lower semicontinuous, L-weakly convex function on Rn . For z ∈ Rn , define the Moreau envelope of ϕ̄ with parameter 1/(2L) and its associated proximal point by n

x⋆0 (z) := argmin ϕ(x) + L∥x − z∥2 .

p0 (z) := min ϕ(x) + L∥x − z∥2 , x∈X

(8)

x∈X

Since 1/(2L) < 1/L, the standard Moreau-envelope theorem for weakly convex functions guarantees that x⋆0 (z) exists and is unique for every z ∈ Rn . Furthermore, p0 is continuously differentiable on Rn , with ∇p0 (z) = 2L z − x⋆0 (z) . (9) We therefore use the optimization-stationarity measure SOS (z) := ∥∇p0 (z)∥.

(10)

A random output zout ∈ Rn is called ε-optimization-stationary in expectation if E SOS (zout )2 ≤ ε2 .

This condition also implies E[SOS (zout )] ≤ ε. The GS criterion measures first-order residuals at a feasible primal–dual pair, whereas the OS criterion measures the gradient of a Moreau envelope of the constrained primal value function. The envelope is defined for every z ∈ Rn , and its proximal point x⋆0 (z) always belongs to X. Although the two criteria are related, their approximate guarantees are not interchangeable without further analysis. We establish the corresponding complexity bounds separately. 5

A stochastic projected damped extragradient method

We develop a stochastic projected damped extragradient (SPDE) method for nonconvex–(strongly) concave minimax optimization. The method combines a projected predictor–corrector step, a damped dual momentum recursion, and a relaxed update of a primal proximal center. These updates are performed within a single loop: each iteration uses a fixed sequence of projections and stochastic gradient queries, without an inner iterative solver for a regularized saddle-point subproblem. In this section, we establish the game-stationarity guarantees; the corresponding optimization-stationarity guarantees are developed in Section 5. Regularized objective.

For a center z ∈ X and a regularization parameter 0 < τ ≤ L, define

τ g(x, y; z, ω) := f (x, y; ω) + L∥x − z∥2 − ∥y∥2 , 2 τ G(x, y; z) := Eω∼D [g(x, y; z, ω)] = F (x, y) + L∥x − z∥2 − ∥y∥2 . 2

(11)

Under Assumption 1, F (·, y) is L-weakly convex on X. Consequently, for each fixed z, G(·, ·; z) is L-strongly convex in x and τ -strongly concave in y. If F (x, ·) is µ-strongly concave, the dual strong concavity modulus of G is µ + τ . These properties concern the population objective G; they need not hold for individual sample objectives g. The regularized objective guides the algorithmic updates. The stationarity criteria remain those of the original problem: the residual R in (6) and the Moreau-envelope measure SOS in (10). The effect of dual regularization on these criteria is accounted for in the subsequent complexity analysis. Stochastic gradient estimates. For a feasible query point u = (x, y) and a batch Ω = (ω1 , . . . , ωm ) of independent samples from D, define b (u; Ω) := 1 ∇F

m X

m i=1

∇f (x, y; ωi ).

(12)

The regularization gradients are evaluated exactly, so stochastic estimation is needed only for ∇F . At iteration t, SPDE queries three fresh batches sequentially: at the current iterate, at a projected predictor, and at the corrected iterate. We use the notation (0)

b ((xt , yt ); Ω ), := ∇F t

(1)

b ((x et , yet ); Ωt ), := ∇F

(2)

b ((xt+1 , yt+1 ); Ω ), := ∇F t

gbt gbt gbt (r)

(r)

(0) (1)

(13) (2)

(r)

and write gbt = (gbx,t , gby,t ). Each batch is drawn independently of all information available before that batch is queried. Single-loop update. SPDE maintains the primal–dual iterate (xt , yt ), the proximal center zt , normal-cone vectors ξt and nt , and a dual momentum vector vt . The first stochastic gradient et , yet ). The second estimate produces the corrected estimate generates the projected predictor (x iterate and the associated normal-cone vectors. It also forms the intermediate momentum v̄t used in the dual correction. After the corrected point has been computed, a third fresh batch updates the dual momentum at that point. This sampling order makes the endpoint gradient error conditionally mean zero given the corrected iterate and its normal-cone vectors. Finally, the relaxed update zt+1 = (1 − β)zt + βxt+1 moves the proximal center toward the corrected primal iterate. Thus the center evolves together with the predictor–corrector and momentum updates, without a separate outer loop. 6

The coupling of the normal-cone corrections with the damped dual momentum is a central feature of SPDE. In particular, the scaling of nt+1 by 1 + k is part of this coupling, and the same vector enters the next momentum update. Algorithm 1 gives the complete procedure. The convergence results below specify admissible parameter choices within the ranges stated in the algorithm. Algorithm 1 Stochastic projected damped extragradient (SPDE) Input: Deterministic x0 ∈ X and y0 ∈ Y ; batch size m ≥ 1 and iteration budget T ≥ 1; 0 < τ ≤ L, h > 0, 0 < α ≤ 1, k ≥ 0, and 0 < β ≤ 1. 1: Set z0 = x0 , ξ0 = 0 ∈ Rn , and n0 = v0 = 0 ∈ Rp . 2: for t = 0, . . . , T − 1 do (0) (0) 3: Draw a fresh batch Ωt of m samples and compute gbt by (13). (0) et = ΠX xt − h gbx,t + 2L(xt − zt ) + ξt . 4: x

(0)

yet = ΠY yt + h gby,t − τ yt − nt + vt

Draw a fresh batch Ωt of m samples and compute gbt (1) v̄t = αvt + k gby,t − τ yet .

(1)

by (13).

e t − zt ) . xt+1 = ΠX xt − h gbx,t + 2L(x xt − xt+1 (1) et − zt ). − gbx,t − 2L(x 9: ξt+1 = h (1) 10: yt+1 = ΠY yt + h gby,t − τ yet + v̄t . 1 yt − yt+1 (1) 11: nt+1 = + gby,t − τ yet + v̄t . 1+k h (2) (2) 12: Draw a fresh batch Ωt of m samples and compute gbt by (13). (2) 13: vt+1 = αvt + k gby,t − τ yt+1 − nt+1 . 14: zt+1 = zt + β(xt+1 − zt ). 15: end for 16: Draw J ∼ Unif{0, . . . , T − 1} independently of all oracle samples. Output: (xout , yout ) = (xJ+1 , yJ+1 ). 8:

The output rule is used for the GS guarantees in this section. It requires no evaluation of the stationarity residual. The output used for the OS guarantees is specified in Section 5.

3.1

Stochastic oracle and batch errors

We next state the stochastic oracle assumptions and describe the information available at each query. These assumptions connect sample gradients to the population gradient; differentiability of the sample losses alone does not assert unbiasedness. Assumption 2. Let A be the sigma-algebra representing the information available before an oracle query, and let u = (x, y) ∈ X ×Y be an A-measurable query point. The oracle draws a fresh sample ω ∼ D, independently of A, and returns a measurable sample gradient ∇f (u; ω) ∈ Rn+p satisfying h i E[∇f (u; ω) | A] = ∇F (u), E ∥∇f (u; ω) − ∇F (u)∥2 A ≤ σ 2 . (14) Here σ ≥ 0 is independent of the query point, iteration index, and regularization parameter. Within each batch, the samples are independent draws from D and are independent of the pre-query information. Each evaluation of ∇f (u; ω) counts as one stochastic first-order oracle (SFO) call.

Conditional independence and unbiasedness imply that the mini-batch estimator satisfies b (u; Ω) | A] = ∇F (u), E[∇F i

b (u; Ω) − ∇F (u)∥2 A ≤ E ∥∇F

(15)

σ2 . m

Indeed, the cross terms between distinct sample-gradient errors have zero conditional expectation. No smoothness or concavity assumption is imposed on individual sample losses beyond the differentiability in Assumption 1. All smoothness and concavity properties used in this section concern the population objective F . In particular, SPDE does not require the mean-square Lipschitz condition used later for VR-SPDE. (1)

Oracle accounting. Every occurrence of gbt within iteration t uses the same computed vector. The third batch is used for the endpoint momentum update, and the first batch of iteration t + 1 is drawn afresh. With three batches of size m per iteration, Algorithm 1, as written, uses 3mT SFO calls. Only the y component of the third estimate enters the updates; if an implementation evaluates only that component, we still charge one full-gradient SFO call per sample. Thus 3mT is also a valid upper bound under this convention. Filtration and conditional errors. Let Ft be the sigma-algebra generated by all oracle samples drawn before iteration t, together with the deterministic initialization. Define (1)

(0)

(2)

:= σ(Ft , Ωt ),

(1)

(2)

Ft+1 := σ(Ft , Ωt ).

:= σ(Ft , Ωt ),

(1)

et , yet ) is Ft -measurable, and The state (xt , yt , zt , ξt , nt , vt ) is Ft -measurable. The predictor (x (2) the corrected quantities xt+1 , yt+1 , ξt+1 , nt+1 and v̄t are Ft -measurable. The center zt+1 is also (2) Ft -measurable, because its update does not depend on the third batch. Define the batch errors by (0)

at := gbt

− ∇F (xt , yt ),

(1) et , yet ), bt := gbt − ∇F (x (2) ct := gbt − ∇F (xt+1 , yt+1 ).

(16)

Equivalently, the stochastic estimates admit the decompositions (0)

gbt

= ∇F (xt , yt ) + at ,

(1) et , yet ) + bt , gbt = ∇F (x (2) gbt = ∇F (xt+1 , yt+1 ) + ct .

(17)

Applying (15) at the three successive query points gives E[at | Ft ] = 0, E[∥at ∥2 | Ft ] ≤

σ2 , m

(1)

(2)

E[bt | Ft ] = 0, (1)

E[∥bt ∥2 | Ft ] ≤

E[ct | Ft ] = 0,

σ2 , m

(2)

E[∥ct ∥2 | Ft ] ≤

σ2 . m

(18)

Writing ct = (cx,t , cy,t ), the endpoint momentum update therefore takes the form

vt+1 = αvt + k ∇y F (xt+1 , yt+1 ) − τ yt+1 − nt+1 + kcy,t ,

(2)

E[cy,t | Ft ] = 0.

This identity states precisely the conditional centering supplied by the third batch. The errors at , bt , and ct need not be mutually independent, since their query points depend on preceding batches; the analysis uses their conditional moment properties in (18). 8

The deterministic saddle-point, envelope, sensitivity, and normal-cone arguments below concern the population objective. For the stochastic updates, we first substitute (17) into the algorithm and derive inequalities for the realized iterates and errors. We then take conditional expectations with respect to the appropriate sigma-algebras. In particular, conditional mean-zero identities are used only when the other factors are measurable with respect to the corresponding pre-query information. Lemma 1. Suppose that Assumptions 1 and 2 hold, and let the iterates be generated by Algorithm 1. Then, almost surely, for every t = 0, . . . , T , xt , zt ∈ X,

yt ∈ Y,

ξt ∈ NX (xt ),

nt ∈ NY (yt ).

(19)

et , yet ) ∈ X × Y for t = 0, . . . , T − 1. At every finite iteration, all state The predictors satisfy (x components, predictors, intermediate momentum vectors, and batch gradient estimates have finite second moments.

Proof. The initialization is feasible, and the zero vector belongs to the normal cone at every point of a nonempty closed convex set. The projection steps ensure that the predictors and corrected iterates are feasible. Since 0 < β ≤ 1, zt+1 = (1 − β)zt + βxt+1 ∈ X whenever zt , xt+1 ∈ X. By the projection optimality condition in (5), the primal correction satisfies (1)

et − zt ) − xt+1 = hξt+1 ∈ NX (xt+1 ). xt − h gbx,t + 2L(x

Similarly, the dual correction gives (1)

yt + h gby,t − τ yet + v̄t − yt+1 = h(1 + k)nt+1 ∈ NY (yt+1 ). Normal cones are cones, and h > 0 and 1 + k > 0. Dividing by these positive scalars proves (19) by induction. To establish the second-moment claim, fix uref ∈ X × Y . The Lipschitz bound (2) implies ∥∇F (u)∥2 ≤ 2∥∇F (uref )∥2 + 2L2 ∥u − uref ∥2 ,

u ∈ X × Y.

For any square-integrable feasible query point u, conditional unbiasedness and (15) yield h

b (u; Ω)∥2 A ≤ ∥∇F (u)∥2 + E ∥∇F

σ2 . m

Hence the corresponding batch estimate has a finite second moment. Moreover, for any nonempty closed convex set C and any fixed uC ∈ C, nonexpansiveness gives ∥ΠC (w) − uC ∥ = ∥ΠC (w) − ΠC (uC )∥ ≤ ∥w − uC ∥. Thus projection preserves finite second moments. Starting from the deterministic initialization, apply these observations successively to the first batch and predictor, the second batch and correction, and the third batch and momentum update. All remaining updates are affine combinations with fixed finite coefficients. Induction establishes the claimed finite second moments at every finite iteration.

3.2

Convergence analysis

We now introduce the regularized saddle point attached to a fixed center and the potential that will telescope across iterations. The first three results establish sensitivity, nonnegativity, and an accuracy-independent initial budget; none of them uses stochastic independence. For a fixed center z ∈ Rn , define p(z) := min max G(x, y; z), x∈X y∈Y

⋆

(20)

⋆

(x (z), y (z)) := the unique saddle point of G(·, ·; z). Existence follows from strong convexity and coercivity in x, compactness of Y , and the convex– concave minimax theorem. Strong convexity and strong concavity give uniqueness. More explicitly, for fixed y ∈ Y , the function G(x, y; z) has a quadratic lower bound. Thus maxy G(x, y; z) has bounded sublevel sets and attains its minimum. The convex–concave minimax theorem can be applied on a compact convex subset containing all relevant minimizers; coercivity then recovers the problem over X. Lemma 2. Suppose that Assumption 1 holds and 0 < τ ≤ L. Then, for all z, z ′ ∈ Rn , s ⋆

′

⋆

∥x (z ) − x (z)∥ ≤ 2∥z − z∥,

⋆

′

⋆

∥y (z ) − y (z)∥ ≤

2L ′ ∥z − z∥, τ

(21)

and ∇p(z) = 2L(z − x⋆ (z)),

∥∇p(z ′ ) − ∇p(z)∥ ≤ 6L∥z ′ − z∥.

(22)

Proof. Define T (x, y) = (∇x F (x, y) + 2Lx, −∇y F (x, y) + τ y). Adding the two pairs of first-order strong convexity and strong concavity inequalities for G gives ⟨T (x′ , y ′ ) − T (x, y), (x′ − x, y ′ − y)⟩ ≥ L∥x′ − x∥2 + τ ∥y ′ − y∥2 . Add the variational inequalities at the two saddle points and write δx = x⋆ (z ′ ) − x⋆ (z) and δy = y ⋆ (z ′ ) − y ⋆ (z). Then L∥δx ∥2 + τ ∥δy ∥2 ≤ 2L⟨z ′ − z, δx ⟩ ≤

L ∥δx ∥2 + 2L∥z ′ − z∥2 . 2

Rearranging proves (21). Let ϕτ (x) = maxy∈Y {F (x, y)−τ ∥y∥2 /2}. Then p(z) = minx∈X {ϕτ (x)+ L∥x − z∥2 }. Uniqueness of the minimizer and envelope differentiation give ∇p(z) = 2L(z − x⋆ (z)). Together with (21), this yields the Lipschitz bound 2L(1 + 2) = 6L. The potential used in the remaining concave analysis employs the shorthand √ s := 2Lτ .

(23)

Write wt = (xt , yt , ξt , nt , vt ). For w = (x, y, ξ, n, v) with ξ ∈ NX (x), n ∈ NY (y), define s 1 Ez (w) := − 1 − [G(x, y; z) − p(z)] + ∥∇x G(x, y; z) + ξ∥2 16L L 1 s τ + ∥−∇y G(x, y; z) + n − v∥2 + ⟨v, y − y ⋆ (z)⟩ + ∥y − y ⋆ (z)∥2 , L 16L 256 s Hz (w) := Ez (w) + ∥v∥2 , 256L2 V(z, w) := p(z) − ′inf n p(z ′ ) + Hz (w), Vt := V(zt , wt ).

z ∈R

(24)

Lemma 3. Suppose that Assumption 1 holds and 0 < τ ≤ L. For any state as above, let P = ∇x G(x, y; z) + ξ and Q = −∇y G(x, y; z) + n − v. Then L τ ∥x − x⋆ (z)∥2 + ∥y − y ⋆ (z)∥2 , 2 2

(25)

2 s 1 s τ 1 Q − (y − y ⋆ (z)) ≥ 0. ∥P ∥2 + ∥x − x⋆ (z)∥2 + ∥y − y ⋆ (z)∥2 + 2L 32 2 L 32

(26)

⟨−∇y G(x, y; z) + n, y − y ⋆ (z)⟩ + G(x, y; z) − p(z) ≥ and Ez (w) ≥

Consequently, V(z, w) ≥ 0. Proof. Saddle-point optimality, strong convexity, and strong concavity give L ∥x − x⋆ (z)∥2 , 2 τ G(x, y ⋆ (z); z) ≤ G(x, y; z) − ⟨∇y G(x, y; z), y − y ⋆ (z)⟩ − ∥y − y ⋆ (z)∥2 . 2

G(x, y ⋆ (z); z) − p(z) ≥

(27) (28)

Adding these inequalities and using ⟨n, y − y ⋆ (z)⟩ ≥ 0 proves (25). Write Y0 = y − y ⋆ (z). Completing the square in (24) yields s 1 s Ez (w) = − 1 − (G − p) + ∥P ∥2 + ⟨−∇y G + n, Y0 ⟩ 16L L 16L 1 τ + ∥Q − sY0 /32∥2 + ∥Y0 ∥2 L 512 1 1 s [G(x, y ⋆ (z); z) − p(z)] + ∥P ∥2 + ∥Q − sY0 /32∥2 . ≥ − (G − p) + 16L L L

The last line follows from (28) after discarding nonnegative terms. For every a ∈ X, strong convexity and the normal-cone condition also give G(a, y; z) ≥ G(x, y; z) + ⟨P, a − x⟩ +

∥P ∥2 L ∥a − x∥2 ≥ G(x, y; z) − . 2 2L

Since mina G(a, y; z) ≤ G(x⋆ (z), y; z) ≤ p(z) − τ ∥Y0 ∥2 /2, we obtain ∥P ∥2 τ − (G(x, y; z) − p(z)) ≥ ∥Y0 ∥2 . 2L 2 Substituting this inequality and (27) into the completed-square lower bound proves (26). Moreover, ∥y∥ ≤ DY implies inf z p(z) = inf x ϕτ (x) ≥ ϕinf − τ DY2 /2 > −∞. Thus Ez , Hz , and V in (24) are all nonnegative. Define the initial quantity, independent of τ, m, T, ε, by ∥∇x F (x0 , y0 )∥2 + 2∥∇y F (x0 , y0 )∥2 1 B := ϕ(x0 ) − ϕinf + 2DY ∥∇y F (x0 , y0 )∥ + + 3+ LDY2 . L 64 (29)

Lemma 4. Suppose that Assumption 1 holds and 0 < τ ≤ L. Then the initialization in Algorithm 1 satisfies 0 ≤ V0 ≤ B. Proof. Using 0 ≤ ϕ(x) − ϕτ (x) ≤ τ DY2 /2 and z0 = x0 , we have τ p(x0 ) − inf p(z) ≤ ϕ(x0 ) − ϕinf + DY2 . z 2

Substituting ξ0 = n0 = v0 = 0 into (24) gives s τ Ex0 (w0 ) = 1 − p(x0 ) − F (x0 , y0 ) + ∥y0 ∥2 16L 2 2 ∥∇x F (x0 , y0 )∥ + ∥−∇y F (x0 , y0 ) + τ y0 ∥2 τ + + ∥y0 − y ⋆ (x0 )∥2 . L 256

By concavity and the radius bound on Y , the bracketed expression is at most 2DY ∥∇y F (x0 , y0 )∥ + τ DY2 /2. Replace the bracket by this nonnegative upper bound and use 0 < 1 − s/(16L) < 1 together with ∥−∇y F (x0 , y0 ) + τ y0 ∥2 ≤ 2∥∇y F (x0 , y0 )∥2 + 2τ 2 DY2 . This gives V0 ≤ ϕ(x0 ) − ϕinf + 2DY ∥∇y F (x0 , y0 )∥ +

∥∇x F (x0 , y0 )∥2 + 2∥∇y F (x0 , y0 )∥2 L

2τ 2 τ + τ+ + DY2 ≤ B. L 64 The last step uses τ ≤ L; nonnegativity follows from Lemma 3. The next lemma is the deterministic algebraic core of the analysis. It isolates one predictor– corrector step at a fixed center and keeps the two gradient errors explicit. In the SPDE proof (0) (1) these errors are the realized sample-gradient errors generated by Ωt and Ωt ; in the VR-SPDE proof they are the errors of the recursive sample estimator. Throughout this subsection, the center z is fixed. All gradients refer to the same function G(·, ·; z), so neither conditional expectations nor center updates arise here. The next lemma makes the concrete step-size and damping choices at their first use. Consider two feasible states w = (x, y, ξ, n, v) and w+ = (x+ , y + , ξ + , n+ , v + ), and assume normal cone membership at both endpoints. Specifically, ξ ∈ NX (x), ξ + ∈ NX (x+ ), n ∈ NY (y), and n+ ∈ NY (y + ). For compact notation, define P = ∇x G(x, y; z) + ξ,

P + = ∇x G(x+ , y + ; z) + ξ + ,

R = −∇y G(x, y; z) + n,

R+ = −∇y G(x+ , y + ; z) + n+ ,

Q = R − v,

Q+ = R + − v + ,

X + = x+ − x⋆ (z),

Y + = y + − y ⋆ (z).

(30)

For any quantity defined at both endpoints, ∆ denotes its endpoint value minus its initial value; for example, ∆P = P + − P and ∆x = x+ − x. Also set I = ∥∆P ∥2 + ∥∆Q∥2 ,

e = (ex , ey ),

∥e∥2 = ∥ex ∥2 + ∥ey ∥2 .

(31)

We use the following full dissipation associated with the potential function (24): Dz (w+ , w) =

h + 2 3hs 7hs ∥P ∥ + ∥Q+ ∥2 + ∥v + ∥2 4 512L 4096L hτ s + 2 hs(16L − s) 1 + ∥Y ∥ + ∥X + ∥2 + I. 128 1024 4L

(32)

Here v + may be the virtual endpoint momentum constructed below. The argument requires the exact recursion stated next, but does not require ex , ey to vanish.

Lemma 5. Suppose that Assumption 1 holds. Let 0 < τ ≤ L, let s be given by (23), and set LG := 3L + τ,

h :=

hs α := 1 + 16

1 , 64LG

−1

d := hLG =

1 , 64

(33)

h(8L − s) k := . 16 + hs

These choices satisfy s ≤ 2L,

hL ≤

1 , 192

α−1 = 1 +

hs , 16

k h(8L − s) = , α 16

0≤k≤

hL . 2

(34)

Under the preceding notation, suppose that ∆x = −h(P + + ex ), +

∆y = −h(Q+ + ey ), +

v = αv + k ∇y G(x , y ; z) − n and that, for some u ≥ 0,

(35)

(36)

√ ∥e∥ ≤ 2d I + u.

(37)

hs 17 Hz (w+ ) ≤ −Dz (w+ , w) + u2 . 32 L

(38)

Then Hz (w+ ) − Hz (w) +

Proof. We suppress the fixed-center subscript z below. Set ρ = (16L−s)/32 and H = G(x, y; z)− p(z) and H + = G(x+ , y + ; z) − p(z). Using (34) and dividing the momentum recursion by α gives ∆v s 8L − s + L s = − v+ − R = − R + + Q+ . h 16 16 2 16

(39)

Step 1: Endpoint inequalities and a three-component expansion. For any state, write Y = y − y ⋆ (z) and define three local components 1 A = −ρH + ∥P ∥2 , 2

1 B = ∥Q∥2 , 2

C0 =

s Lτ ⟨v, Y ⟩ + ∥Y ∥2 . 32 512

(40)

The definition of the potential function gives A + B + C0 = (L/2)Ez (w). The kinetic energy is excluded from C0 for now and will be added in the final step. Since ξ ∈ NX (x) and n+ ∈ NY (y + ), ⟨ξ, ∆x⟩ ≤ 0 and ⟨n+ , ∆y⟩ ≥ 0. At the mixed point + (x , y), applying L-strong convexity and τ -strong concavity gives L ∥∆x∥2 , 2 τ G(x+ , y; z) ≤ G(x+ , y + ; z) + ⟨R+ , ∆y⟩ − ∥∆y∥2 . 2 G(x+ , y; z) ≥ G(x, y; z) + ⟨P, ∆x⟩ +

Consequently, τ L ∥∆x∥2 + ∥∆y∥2 . (41) 2 2 To establish the strong monotonicity used below, write the following inequalities at the four endpoints: H + − H ≥ ⟨P, ∆x⟩ − ⟨R+ , ∆y⟩ +

L ∥∆x∥2 , 2 L G(x, y + ; z) − G(x+ , y + ; z) ≥ −⟨∇x G(x+ , y + ; z), ∆x⟩ + ∥∆x∥2 , 2 τ + G(x, y ; z) − G(x, y; z) ≤ ⟨∇y G(x, y; z), ∆y⟩ − ∥∆y∥2 , 2 τ + + + + + G(x , y; z) − G(x , y ; z) ≤ −⟨∇y G(x , y ; z), ∆y⟩ − ∥∆y∥2 . 2 G(x+ , y; z) − G(x, y; z) ≥ ⟨∇x G(x, y; z), ∆x⟩ +

The sum of the first two left-hand sides equals the sum of the last two. Comparing the corresponding right-hand sides and rearranging yields ⟨∆(∇x G), ∆x⟩ + ⟨∆(−∇y G), ∆y⟩ ≥ L∥∆x∥2 + τ ∥∆y∥2 . The normal cone definition also gives ⟨∆ξ, ∆x⟩ ≥ 0 and ⟨∆n, ∆y⟩ ≥ 0. Adding these inequalities gives ⟨∆P, ∆x⟩ + ⟨∆R, ∆y⟩ ≥ L∥∆x∥2 + τ ∥∆y∥2 . (42) For the primal part of A, the difference-of-squares identity gives exactly 1 Lρ (∥P + ∥2 − ∥P ∥2 ) − ρ⟨P, ∆x⟩ − ∥∆x∥2 2 2 1 L(128L − τ ) = ⟨P + , ∆P ⟩ − ρ⟨P + , ∆x⟩ − ∥∆P − ρ∆x∥2 − ∥∆x∥2 . 2 1024

(43)

The last term is nonpositive. Multiplying (41) by −ρ and applying (43), we obtain ∆A ≤ − ρ⟨P + , ∆x⟩ + ρ⟨R+ , ∆y⟩ + ⟨P + , ∆P ⟩ 1 ρτ − ∥∆P − ρ∆x∥2 − ∥∆y∥2 . 2 2

(44)

Moreover, 1 ∆B = ⟨Q+ , ∆Q⟩ − ∥∆Q∥2 . 2 Substituting (35) into (42) and using ∆R = ∆Q + ∆v yields 1 ⟨P + , ∆P ⟩ + ⟨Q+ , ∆Q⟩ ≤ − L∥P + + ex ∥2 − τ ∥Q+ + ey ∥2 h 1 1 − ⟨ex , ∆P ⟩ + ⟨ey , ∆Q⟩ − ⟨Q+ + ey , ∆v⟩. h h

(45)

(46)

By (39), the last inner product equals 1 L s − ⟨Q+ + ey , ∆v⟩ = ⟨Q+ + ey , R+ ⟩ − ⟨Q+ + ey , Q+ ⟩. h 2 16 Combining (44)–(46) and using L/2 − ρ = s/32, we obtain the full two-component increment bound ∆A + ∆B ≤ ρ∥P + ∥2 − L∥P + + ex ∥2 − τ ∥Q+ + ey ∥2 + ρ⟨P + , ex ⟩ h s s + ⟨R+ , Q+ + ey ⟩ − ⟨Q+ , Q+ + ey ⟩ 32 16 1 1 − ⟨ex , ∆P ⟩ + ⟨ey , ∆Q⟩ − ∥∆P − ρ∆x∥2 h 2h 1 ρτ − ∥∆Q∥2 − ∥∆y∥2 . 2h 2h

(47)

For the third component, start with the two exact identities ⟨v + , Y + ⟩ − ⟨v, Y ⟩ = ⟨∆v, Y + ⟩ + ⟨v + , ∆y⟩ − ⟨∆v, ∆y⟩, ∥Y + ∥2 − ∥Y ∥2 = 2⟨Y + , ∆y⟩ − ∥∆y∥2 . Substituting (35) and (39), and then using v + = R+ − Q+ and s2 /512 = Lτ /256, gives the termwise expansion ∆C0 s s s = − ⟨R+ , Q+ + ey ⟩ + ∥Q+ ∥2 + ⟨Q+ , ey ⟩ h 32 32 32 Ls + + Lτ s Lτ − ⟨R , Y ⟩ − ⟨Y + , ey ⟩ − ⟨∆v, ∆y⟩ − ∥∆y∥2 . 64 256 32h 512h 14

(48)

The two ⟨Q+ , Y + ⟩ terms cancel because s2 /512 = Lτ /256. In (47) and (48), the ⟨R+ , Q+ + ey ⟩ terms also cancel exactly. Step 2: Add the endpoint weight and bound the current-point terms. The endpoint weight expands as Ls s(16L − s) + s Ez (w+ ) = − H + (∥P + ∥2 + ∥Q+ ∥2 ) 64 1024 64 Lτ + + Lτ s + ⟨v , Y ⟩ + ∥Y + ∥2 . 512 16384

(49)

After writing v + = R+ − Q+ , the coefficient of ⟨R+ , Y + ⟩ is −

Ls Lτ s(16L − s) + =− . 64 512 1024

The saddle-point growth inequality (25) therefore gives −

s(16L − s) + H + ⟨R+ , Y + ⟩ 1024 Ls(16L − s) τ s(16L − s) + 2 ≤− ∥X + ∥2 − ∥Y ∥ . 2048 2048

(50)

Adding (47), (48), and (49), and applying (50), yields L hs Ez (w+ ) − Ez (w) + Ez (w+ ) ≤ C + J , 2h 32

(51)

where the current-point and increment terms are, respectively, s 32L − s + 2 ∥P ∥ − L∥P + + ex ∥2 − τ ∥Q+ + ey ∥2 − ∥Q+ ∥2 64 64 s Lτ Lτ + ρ⟨P + , ex ⟩ − ⟨Q+ , ey ⟩ − ⟨Y + , ey ⟩ − ⟨Q+ , Y + ⟩ 32 256 512 τ s(16L − s) Ls(16L − s) Lτ s − ∥Y + ∥2 − ∥X + ∥2 , + 16384 2048 2048 1 1 ρτ Lτ J = − ∥∆P − ρ∆x∥2 − ∥∆Q∥2 − + ∥∆y∥2 2h 2h 2h 512h 1 s − ⟨∆v, ∆y⟩ − ⟨ex , ∆P ⟩ + ⟨ey , ∆Q⟩ . 32h h C=

(52)

(53)

These expressions include every term in the sum; no term involving u has been discarded. For the primal terms, expand exactly and then use s ≤ 2L: 32L − s + 2 ∥P ∥ − L∥P + + ex ∥2 + ρ⟨P + , ex ⟩ 64 32L + s + 2 48L + s + =− ∥P ∥ − ⟨P , ex ⟩ − L∥ex ∥2 64 32 L ≤ − ∥P + ∥2 + 4L∥ex ∥2 . 4

(54)

The last step uses 2L∥P + ∥∥ex ∥ ≤ (L/4)∥P + ∥2 + 4L∥ex ∥2 and drops the remaining nonpositive terms. Young’s inequality bounds the other three cross terms as follows: s s s |⟨Q+ , ey ⟩| ≤ ∥Q+ ∥2 + ∥ey ∥2 , 32 256 16 Lτ s Lτ s + + + 2 |⟨Q , Y ⟩| ≤ ∥Q ∥ + ∥Y + ∥2 , 512 256 8192 Lτ Lτ s s |⟨Y + , ey ⟩| ≤ ∥Y + ∥2 + ∥ey ∥2 . 256 16384 32 15

(55) (56) (57)

The constants follow directly from s2 = 2Lτ . The resulting coefficient of ∥Q+ ∥2 is −s/128, while the coefficient of the distance term satisfies Lτ s τ s(16L − s) 27Lτ s Lτ s − ≤− ≤− . 4096 2048 4096 256 Dropping −τ ∥Q+ + ey ∥2 and using 3s/32 ≤ 4L, we obtain L + 2 s Lτ s + 2 ∥P ∥ − ∥Q+ ∥2 − ∥Y ∥ 4 128 256 Ls(16L − s) − ∥X + ∥2 + 4L∥e∥2 . 2048

C≤ −

(58)

Step 3: Bound the increment terms and absorb the explicit errors. The inequality ∥a − b∥2 ≥ ∥a∥2 /2 − ∥b∥2 and (35) give −

1 1 ∥∆P − ρ∆x∥2 − ∥∆Q∥2 2h 2h I ≤− + hρ2 (∥P + ∥2 + ∥ex ∥2 ). 4h

(59)

Since ∆v = −∆(∇y G) + ∆n − ∆Q and ⟨∆n, ∆y⟩ ≥ 0, applying LG -Lipschitz continuity only to the true gradient gives s LG s q s ⟨∆v, ∆y⟩ ≤ ∥∆x∥2 + ∥∆y∥2 ∥∆y∥ + ∥∆Q∥∥∆y∥ 32h 32h 32h 3LG s Lτ 1 LG s ∥∆x∥2 + + ∥∆y∥2 + ∥∆Q∥2 ≤ 64h 64h 64h 32h hLG s 3hLG s hLτ I ≤ (∥P + ∥2 + ∥ex ∥2 ) + + (∥Q+ ∥2 + ∥ey ∥2 ) + . 32 32 32 32h √ The intermediate step uses a2 + b2 b ≤ a2 /2 + 3b2 /2 and −

s∥∆Q∥∥∆y∥ ≤ ∥∆Q∥2 +

(60)

s2 ∥∆y∥2 . 4

The bounds s ≤ 2L ≤ 2LG and τ ≤ L imply hρ2 +

hLG s 5Ld ≤ , 32 16

3hLG s hLτ ds + ≤ . 32 32 8

(61)

The first inequality uses ρ ≤ L/2, and the second uses Lτ ≤ LG s. Both right-hand sides are at most L/2. We retain the explicit cost of u in the error inner products. By (37), √ 1 ∥e∥ I 2d + 1/32 8 − ⟨ex , ∆P ⟩ + ⟨ey , ∆Q⟩ ≤ ≤ I + u2 , (62) h h h h ∥e∥2 ≤ 8d2 I + 2u2 . (63) √ The first inequality uses u I ≤ I/32 + 8u2 . Substituting (59)–(62) into (53) and dropping the explicit nonpositive ∥∆y∥2 term gives J ≤

5Ld + 2 ds + 2 L ∥P ∥ + ∥Q ∥ + ∥e∥2 16 8 2 7/32 − 2d − 1/32 8 2 − I+ u . h h

(64)

Adding (58) bounds the right-hand side of (51) by (4 − 5d)L + 2 (1 − 16d)s + 2 Lτ s + 2 ∥P ∥ − ∥Q ∥ − ∥Y ∥ 16 128 256 Ls(16L − s) 9L 7/32 − 2d − 1/32 8 − ∥X + ∥2 + ∥e∥2 − I + u2 . 2048 2 h h −

(65)

Substituting (63) and using Lh ≤ d, the magnitude of the negative increment coefficient is at least 1 1 7 3 − 2d − − 36d . h 32 32 At d = 1/64, the three required numerical margins are 4 − 5d 251 1 = ≥ , 16 1024 8

1 − 16d 3 1 = ≥ , 128 512 256

7 1 10231 1 − 2d − − 36d3 = ≥ . (66) 32 32 65536 8

Multiplying both sides of (51) by 2h/L therefore gives hs Ez (w+ ) 32 h hs hτ s + 2 hs(16L − s) ≤ − ∥P + ∥2 − ∥Q+ ∥2 − ∥Y ∥ − ∥X + ∥2 4 128L 128 1024 1 16 − I + 18h + u2 . 4L L

Ez (w+ ) − Ez (w) +

(67)

Since hL ≤ 1/64, the last positive coefficient is at most 17/L. Step 4: Add the kinetic-energy dissipation and endpoint weight. Equation (39) and R+ = Q+ + v + give the exact filter identity hL + h(8L − s) + 1+ v =v− Q . 2 16

(68)

Set a = h(8L − s)/16. Then 0 ≤ a ≤ hL/2 and v − v + = (hL/2)v + + aQ+ . The difference-ofsquares identity and Young’s inequality yield ∥v∥2 − ∥v + ∥2 = 2⟨v + , v − v + ⟩ + ∥v − v + ∥2 ≥ hL∥v + ∥2 + 2a⟨v + , Q+ ⟩ hL + 2 hL + 2 ≥ ∥v ∥ − ∥Q ∥ . 2 2

(69)

Multiplying by s/(256L2 ), rearranging, and adding the endpoint weight of the kinetic energy gives s hs s (∥v + ∥2 − ∥v∥2 ) + ∥v + ∥2 2 256L 32256L2 hs hs s ≤ ∥Q+ ∥2 − 1− ∥v + ∥2 512L 512L 16L hs 7hs ≤ ∥Q+ ∥2 − ∥v + ∥2 . 512L 4096L

(70)

The last step uses s ≤ 2L. Adding (70) to (67) and using Hz (w) = Ez (w) + s∥v∥2 /(256L2 ), the coefficient of the dual direction becomes −hs/(128L) + hs/(512L) = −3hs/(512L). All remaining coefficients agree with (32), proving (38). The entire argument is pathwise and therefore requires neither independence nor unbiasedness of ex , ey .

The step-size and damping choices in (33) remain in force throughout the rest of the concave analysis. The actual endpoint momentum contains the third stochastic estimate. To take conditional expectations without losing its linear error terms, we compare it with a virtual endpoint state in which only that final estimate is replaced by its conditional mean. For the analysis, define the virtual endpoint momentum before drawing the third batch by 0 vt+1 := αvt + k ∇y F (xt+1 , yt+1 ) − τ yt+1 − nt+1 ,

0 0 wt+1 := (xt+1 , yt+1 , ξt+1 , nt+1 , vt+1 ). (71)

(2)

0 . This quantity is F The algorithm need not compute vt+1 t -measurable, whereas the actual momentum satisfies 0 vt+1 = vt+1 + kcy,t . (72)

Lemma 6. Suppose that Assumption 1 holds, and let the iterates be generated by Algorithm 1 with the parameter choices in (33). Apply the notation of the preceding subsection with fixed 0 . Let I 0 and D 0 denote the corresponding sum of center zt , initial state wt , and endpoint wt+1 t t (0) (1) squared increments and dissipation, respectively. For every realization of Ωt and Ωt , define (0)

(1)

Ut (Ωt , Ωt ) := 2d∥at ∥ + 2∥bt ∥,

(73)

where at and bt are the explicit sample averages in (16). Then 0 Hzt (wt+1 ) − Hzt (wt ) ≤ −Dt0 +

17 2 U , L t

(0)

(1)

Ut := Ut (Ωt , Ωt ).

(74)

Proof. Fix the two realized batches. By (16), their errors are m 1 X (0) at = [∇f (xt , yt ; ωt,i ) − ∇F (xt , yt )], m i=1

bt =

m 1 X (1) et , yet ; ωt,i ) − ∇F (x et , yet )]. [∇f (x m i=1

No expectation is taken in this proof. The correction step and (71) yield the dynamics xt+1 − xt = −h(P + + ex,t ), yt+1 − yt = −h(Q+ + ey,t ),

(75)

et , yet ; zt ) − ∇x G(xt+1 , yt+1 ; zt ) + bx,t , ex,t = ∇x G(x

et , yet ; zt ) − by,t . ey,t = (1 + k) ∇y G(xt+1 , yt+1 ; zt ) − ∇y G(x

(76)

0 . The second Here P + = ∇x G(xt+1 , yt+1 ; zt ) + ξt+1 and Q+ = −∇y G(xt+1 , yt+1 ; zt ) + nt+1 − vt+1 dynamical identity uses 0 et , yet ; zt ) + by,t − ∇y G(xt+1 , yt+1 ; zt ) + nt+1 v̄t − vt+1 = k ∇y G(x

and the definition of nt+1 , with the normal-cone coefficient equal to k − (1 + k) = −1. Write et = (ex,t , ey,t ). Subtracting (75) from the unprojected predictor and applying the nonexpansiveness of the projection gives et − xt+1 , yet − yt+1 )∥ ≤ h ∥(x

It0 + ∥et ∥ + ∥at ∥ .

Thus, setting r = d(1 + k), we obtain from (76) that et − xt+1 , yet − yt+1 )∥ + (1 + k)∥bt ∥ ∥et ∥ ≤ (1 + k)LG ∥(x

≤r

It0 + ∥et ∥ + ∥at ∥

+ (1 + k)∥bt ∥.

Since k ≤ 1/384 and d = 1/64, we have r < 1, r/(1 − r) ≤ 2d, and (1 + k)/(1 − r) ≤ 2. Rearranging yields q ∥et ∥ ≤ 2d It0 + Ut .

(77)

The virtual momentum satisfies the exact recursion required in the preceding subsection. Substituting (77) into the fixed-center descent inequality and dropping the nonnegative weighted 0 ) gives (74). endpoint term (hs/32)Hzt (wt+1 We next restore the center update. The center sensitivity bound converts its motion into a controlled perturbation of the fixed-center potential, after which the three fresh batches yield a one-step conditional descent inequality. Lemma 7. Suppose that Assumption 1 holds and 0 < τ ≤ L. For any state w = (x, y, ξ, n, v) satisfying (19) and any z ′ = z + δ, we have s 1 ∥y − y ⋆ (z)∥ ∥δ∥ H (w) − Hz (w) ≤ 2L∥x − x (z)∥ + 4∥∇x G(x, y; z) + ξ∥ + ∥v∥ + 8 128 (78) 1025L 2 + ∥δ∥ . 128

⋆

z′

Proof. First, (11) and (22) give G(x, y; z ′ ) − G(x, y; z) = −2L⟨x − z, δ⟩ + L∥δ∥2 , |p(z ′ ) − p(z) − ⟨∇p(z), δ⟩| ≤ 3L∥δ∥2 . Using 2L(x−z) = 2L(x−x⋆ (z))−∇p(z), the increment of the function-value term in the potential function is at most 2L∥x − x⋆ (z)∥∥δ∥ + 4L∥δ∥2 . Next, ∇x G(x, y; z ′ ) = ∇x G(x, y; z) − 2Lδ, so 1 ∥∇x G(x, y; z ′ ) + ξ∥2 − ∥∇x G(x, y; z) + ξ∥2 ≤ 4∥∇x G(x, y; z) + ξ∥∥δ∥ + 4L∥δ∥2 . L p

Finally, (21) and s 2L/τ = 2L imply s 1 |⟨v, y ⋆ (z) − y ⋆ (z ′ )⟩| ≤ ∥v∥∥δ∥, 16L 8 τ s L ∥y − y ⋆ (z ′ )∥2 − ∥y − y ⋆ (z)∥2 ≤ ∥y − y ⋆ (z)∥∥δ∥ + ∥δ∥2 . 256 128 128 The remaining terms are independent of z. Adding the four increments proves (78). Lemma 8. Suppose that Assumption 1 holds, and let the iterates be generated by Algorithm 1 with the parameter choices in (33). Choose the center relaxation β := (0)

hs , 4096

0<β<

1 . 12

(79)

(1)

For every realization of Ωt and Ωt , the following inequality holds before any expectation is taken: 1 β 17 0 V(zt+1 , wt+1 ) − Vt ≤ − Dt0 − ∥∇p(zt )∥2 + Ut2 . (80) 2 8L L Proof. Let Xt+ = xt+1 − x⋆ (zt ), Yt+ = yt+1 − y ⋆ (zt ), and Pt+ = ∇x G(xt+1 , yt+1 ; zt ) + ξt+1 . Define the dissipation budget Bt :=

hs + 2 7hs hτ s + 2 7hLs 0 ∥2 + ∥Pt ∥ + ∥vt+1 ∥Y ∥ + ∥Xt+ ∥2 . 8L 4096L 128 t 512

(81)

Since s ≤ 2L and 16L − s ≥ 14L, comparison with each term of the fixed-center dissipation gives 0 ≤ Bt ≤ Dt0 . The weighted Cauchy–Schwarz inequality yields 2 1 0 s 432L 2L∥Xt+ ∥ + 4∥Pt+ ∥ + ∥vt+1 ∥+ ∥Yt+ ∥ ≤ Bt . 8 128 hs

(82)

The constant follows by summing the four coefficients: 16L hs

128 4 1 +8+ + 7 7 1024

432L . hs

The center increment δt = zt+1 − zt satisfies

δt = β

∇p(zt ) Xt+ − 2L

∥δt ∥2 ≤ 2β 2 ∥Xt+ ∥2 +

β2 ∥∇p(zt )∥2 . 2L2

(83)

Substituting (82) into (78) and using ab ≤ a2 /4 + b2 and hs ≤ 32 gives 1 704L 0 0 ) − Hzt (wt+1 ) ≤ Bt + Hzt+1 (wt+1 ∥δt ∥2 4 hs 1408Lβ 2 352β 2 1 ∥Xt+ ∥2 + ∥∇p(zt )∥2 . ≤ Bt + 4 hs Lhs

(84)

Moreover, p is 6L-smooth, so (83) and ⟨∇p(zt ), Xt+ ⟩ ≤ ∥∇p(zt )∥2 /(8L) + 2L∥Xt+ ∥2 imply β 3 3 − β ∥∇p(zt )∥2 + βL(2 + 6β)∥Xt+ ∥2 L 8 2 5 β ≤ − ∥∇p(zt )∥2 + βL∥Xt+ ∥2 . 4L 2

p(zt+1 ) − p(zt ) ≤ −

(85)

Add (74), (84), and (85). Since β = hs/4096, we have 5 1408Lβ 2 7hLs βL + ≤ , 2 hs 2048 β 352β 2 β 1 352 β − = − ≥ . 4L Lhs L 4 4096 8L The right-hand side of the first line is the coefficient in Bt /4 multiplying ∥Xt+ ∥2 . Hence, the sum of the right-hand sides is at most 1 β 17 −Dt0 + Bt − ∥∇p(zt )∥2 + Ut2 . 2 8L L Using Bt ≤ Dt0 now gives (80). The center choice in (79) remains in force throughout the rest of the concave analysis. Proposition 1. Suppose that Assumption 1 holds, and let the iterates be generated by Algo(0) (1) (2) rithm 1 with the parameter choices in (33) and (79). For the three realized batches Ωt , Ωt , Ωt , define (2)

2k 0 ⟨−∇y G(xt+1 , yt+1 ; zt+1 ) + nt+1 − vt+1 , cy,t ⟩ L sk sk + ⟨cy,t , yt+1 − y ⋆ (zt+1 )⟩ + ⟨v 0 , c ⟩ 2 t+1 y,t 16L 128L 1 s + + k 2 ∥cy,t ∥2 , L 256L2

Λt (Ωt ) := −

(86)

where cy,t =

m 1 X (2) [∇y f (xt+1 , yt+1 ; ωt,i ) − ∇y F (xt+1 , yt+1 )]. m i=1

Then, before taking any expectation, β 17 1 (1) (2) (0) ∥∇p(zt )∥2 + Ut (Ωt , Ωt )2 + Λt (Ωt ). Vt+1 − Vt ≤ − Dt0 − 2 8L L

(87)

Proof. For fixed z, x, y, ξ, n, the function Hz (w) is quadratic in v, with quadratic coefficient 0 1/L + s/(256L2 ). Substitute the realized third-batch identity vt+1 = vt+1 + kcy,t from (72) and expand the square. This gives the exact identity (2)

0 Vt+1 − V(zt+1 , wt+1 ) = Λt (Ωt ).

(88)

Adding this identity to the sample-path virtual descent (80) proves (87). Theorem 1. Suppose that Assumptions 1 and 2 hold. Then Algorithm 1, with the parameter choices in (33) and (79), satisfies 1 β 300σ 2 E[Vt+1 | Ft ] ≤ Vt − E[Dt0 | Ft ] − ∥∇p(zt )∥2 + . 2 8L Lm

(89)

Proof. We now take conditional expectations in the sample-path inequality (87). This is the first point in the descent argument at which an expectation is taken. In the first three inner (2) products, all factors other than cy,t are Ft -measurable. Thus, (18) gives the exact conditional expectation identity (2) (2) E[Λt (Ωt ) | Ft ] (90) s 2σ 2 1 (2) 2 2 + k E[∥c ∥ | F ] ≤ . = y,t t L 256L2 Lm The last step uses s ≤ 2L and k ≤ 1; this loose constant suffices. For the first and second batches, we do not invoke unbiased cancellation involving subsequent iterates. Since Ut = 2d∥at ∥ + 2∥bt ∥, for every sample realization we have Ut2 ≤ 8d2 ∥at ∥2 + 8∥bt ∥2 . Taking conditional expectations and using the tower property gives 17 136(1 + d2 )σ 2 E[Ut2 | Ft ] ≤ . L Lm

(91)

Substitute (90) and (91) into the conditional expectation of (87). Since 136(1+1/4096)+2 < 300, this proves (89). Lemma 1, the compactness of Y , and the at most quadratic growth of p, G ensure the integrability of these random quantities, so the conditional expectations and the tower property are well defined.

3.3

Stationarity and complexity results

The descent inequality controls regularized residuals along the trajectory. The following budget lemma averages those residuals, removes the artificial dual regularization, and converts the result into the original criterion defined in (6). Lemma 9. Suppose that Assumptions 1 and 2 hold, and let the iterates be generated by Algorithm 1 with the parameter choices in (33) and (79). Let MT := B + 21

300T σ 2 . Lm

(92)

Then

TX −1 t=0

β 1 EDt0 + E∥∇p(zt )∥2 ≤ MT . 2 8L

(93)

The random output satisfies ER(xJ+1 , yJ+1 )2 ≤

64LB 19200σ 2 + + 2τ 2 DY2 . βT βm

(94)

Proof. Taking total expectations in (89) and summing from t = 0 to T − 1, we obtain (93) from VT ≥ 0 and V0 ≤ B. Write Pt+ = ∇x G(xt+1 , yt+1 ; zt ) + ξt+1 , 0 Q+ t = −∇y G(xt+1 , yt+1 ; zt ) + nt+1 − vt+1 ,

Xt+ = xt+1 − x⋆ (zt ). Retaining the squared terms in Dt0 individually and using hs = 4096β and 16L − s ≥ 14L, we obtain X

8MT , h

E∥Pt+ ∥2 ≤

t<T

2 E∥Q+ t ∥ ≤

t<T

MT E∥Xt+ ∥2 ≤ , 28βL t<T X

LMT , 12β

0 E∥vt+1 ∥2 ≤

t<T

2LMT , 7β

8LMT E∥∇p(zt )∥ ≤ . β t<T X

(95)

Define the following true-gradient certificate solely for the analysis: St := ∥∇x F (xt+1 , yt+1 ) + ξt+1 ∥2 + ∥−∇y F (xt+1 , yt+1 ) + τ yt+1 + nt+1 ∥2 .

(96)

The certificate St contains neither the actual nor the virtual momentum. By (22), ∇x F (xt+1 , yt+1 ) + ξt+1 = Pt+ − 2LXt+ + ∇p(zt ), 0 −∇y F (xt+1 , yt+1 ) + τ yt+1 + nt+1 = Q+ t + vt+1 .

Hence 2 0 2 St ≤ 3∥Pt+ ∥2 + 12L2 ∥Xt+ ∥2 + 3∥∇p(zt )∥2 + 2∥Q+ t ∥ + 2∥vt+1 ∥ .

Summing and substituting (95) gives X t<T

ESt ≤

24β 3 1 4 + + 24 + + hL 7 6 7

LMT 32LMT ≤ . β β

(97)

The last inequality uses β/(hL) = s/(4096L) ≤ 1/2048. √ Normal-cone membership and the triangle inequality give, pathwise, R(xt+1 , yt+1 ) ≤ St + τ ∥yt+1 ∥, and thus R(xt+1 , yt+1 )2 ≤ 2St + 2τ 2 DY2 . The independent uniform index J converts the finite average into the output expectation. Substituting (97) and (92) yields (94). 3.3.1

Nonconvex–strongly concave setting

We first specialize the mini-batch analysis to the case in which the original inner objective is strongly concave. This setting does not require artificial dual regularization and therefore has no regularization-bias term.

Assumption 3. There exists µ > 0 such that, for every x ∈ X and y, y ′ ∈ Y , F (x, y ′ ) ≤ F (x, y) + ⟨∇y F (x, y), y ′ − y⟩ −

µ ′ ∥y − y∥2 . 2

(98)

Define the condition number κ := L/µ. The joint L-smoothness in Assumption 1 implies 0 < µ ≤ L and hence κ ≥ 1. Define the unregularized proximal saddle function and its value by Gsc (x, y; z) := F (x, y) + L∥x − z∥2 ,

psc (z) := min max Gsc (x, y; z). x∈X y∈Y

(99)

The function Gsc is L-strongly convex in x, µ-strongly concave in y, and has a 3L-Lipschitz full gradient. The concrete strongly-concave specialization is stated in the first result that uses it. Lemma 10. Suppose that Assumptions 1, 2, and 3 hold. Set LG,sc := 3L,

ssc :=

hsc :=

2Lµ,

1 1 = , 64LG,sc 192L

hsc ssc −1 hsc (8L − ssc ) , ksc := , 16 16 + hsc ssc √ 2 hsc ssc βsc := = κ−1/2 . 4096 786432

(100)

αsc := 1 +

Run Algorithm 1 with every explicit dual-regularization term −τ yt , −τ yet , or −τ yt+1 deleted and with the parameters in (100). This prescription specifies all updates and introduces neither an inner solve nor a restart. Then SPDE satisfies ER(xJ+1 , yJ+1 )2 ≤

64LB 19200σ 2 + . βsc T βsc m

(101)

Proof. We give the exact substitution map from the preceding mini-batch analysis. Let V sc denote the potential in (24) after replacing G and p by Gsc and psc , deleting every explicit dual-regularization gradient, and replacing the strong-concavity parameter in the analytical p coefficients by µ. In particular, the sensitivity factor becomes 2L/µ, the potential uses ssc ⟨v, y − y ⋆ (z)⟩/(16L) and µ∥y − y ⋆ (z)∥2 /256, and the dissipation uses the same formulas with (s, τ, h, β) replaced by (ssc , µ, hsc , βsc ). Assumption 3 and (100) give 0 < µ ≤ L,

ssc ≤ 2L,

hsc LG,sc =

1 , 64

βsc =

hsc ssc . 4096

These are precisely the scalar relations used in the fixed-center, center-movement, and mini-batch error estimates. Strong concavity (98) supplies every dual quadratic-growth term, while the algorithmic dual gradients are those of Gsc and contain no regularization term. Thus the proof of Theorem 1 gives the corresponding strongly-concave drift inequality with the same constants. The initial potential is bounded by the same B in (29). Indeed, z0 = x0 and the unregularized saddle function give V0sc ≤ ϕ(x0 ) − ϕinf + 2DY ∥∇y F (x0 , y0 )∥ +

∥∇x F (x0 , y0 )∥2 + ∥∇y F (x0 , y0 )∥2 µDY2 + L 64

≤ B, where the first inequality follows from concavity and the radius bound on Y , and the second uses µ ≤ L and (29). Finally, the analytical certificate contains −∇y F (xt+1 , yt+1 ) + nt+1 itself. Its conversion to R is exact, so the term 2τ 2 DY2 in (94) is absent. The remaining two terms give (101). 23

Theorem 2. Suppose that Assumptions 1, 2, and 3 hold. Fix a known bound B ≥ B and set (

128LB T = max 1, βsc ε2

(

38400σ 2 m = max 1, βsc ε2

(102)

Then SPDE returns a point satisfying (7), and !

38400σ 2 1+ βsc ε2 ! √ (LB + σ 2 ) κ LBσ 2 κ =O 1+ + . ε2 ε4

128LB Nsc ≤ 3 1 + βsc ε2

For fixed L, B and σ > 0, this gives √ √ T = O( κε−2 ), m = O( κε−2 ), √ If σ = 0, then m = 1 and Nsc = O( κε−2 ).

(103)

Nsc = O(κε−4 ).

(104)

Proof. Substituting (102) into (101) gives ε2 64LB ≤ , βsc T 2

19200σ 2 ε2 ≤ . βsc m 2

Hence (7) holds. three batches per iteration give the first line of (103). The identity √ The √ −1 βsc = (786432/ 2) κ from (100) yields its second line and (104). Theorem 2 shows that strong concavity removes the regularization bias and fixes the center timescale at βsc = Θ(κ−1/2 ). The resulting product of the iteration and batch orders gives the O(κε−4 ) stochastic-oracle bound. 3.3.2

Nonconvex–concave setting

Without strong concavity, we set the auxiliary dual regularization at the accuracy-dependent scale τ = O(ε). This makes its contribution to the original game-stationarity residual at most order ε, while slowing the center by the factor β = Θ(ε1/2 ). Theorem 3. Suppose that Assumptions 1 and 2 hold. Choose a known upper bound independent of ε, τ satisfying B ≥ B, where B is defined in (29). Given ε > 0, set ε τ = min L, 2DY

(105)

choose h, α, k by (33) and β by (79), and set (

256LB T = max 1, βε2

(

76800σ 2 m = max 1, βε2

(106)

Then the output of Algorithm 1 satisfies (7). The number N of stochastic full-gradient oracle calls satisfies the explicit bound 256LB N ≤ 3mT ≤ 3 1 + βε2 Let

 

Kε := max 1, 

76800σ 2 1+ βε2

 2LDY  . ε 

(107)

(108)

Then

(LB + σ 2 )Kε LBσ 2 Kε2 N =O 1+ + ε2 ε4

(109)

For a fixed problem, initialization, B, and σ > 0, as ε ↓ 0, the iteration count is T = O(ε−5/2 ) and the batch size is m = O(ε−5/2 ), so N = O(ε−5 ). Proof. By (105), τ DY ≤ ε/2 and 0 < τ ≤ L. Substituting (106) into (94), the three error terms satisfy 64LB 64LB ε2 19200σ 2 ε2 ε2 ≤ ≤ , ≤ , 2τ 2 DY2 ≤ . (110) βT βT 4 βm 4 2 When σ = 0, the second error term is zero and m = 1; no division by σ is needed. Adding the three bounds gives ER2 ≤ ε2 . For any a ≥ 0, max{1, ⌈a⌉} ≤ 1 + a. Thus the three batch queries per iteration give (107). By (33) and (79), √

2Lτ 786432 β= , √ 262144(3L + τ ) 2

L 1 1048576 ≤ ≤ √ τ β 2

L , τ

L = Kε . τ

(111)

Expanding (107) and substituting (111) yields (109). In particular, for 0 < ε ≤ 2LDY ,

N = O 1 + (LB + σ 2 ) LDY ε−5/2 + L2 BDY σ 2 ε−5 . p

(112)

This proves the stated sample complexity. The quantity B is used only in the analysis; the budgets use the supplied B and require no additional true-gradient oracle calls to estimate B. Remark 1. The theorem provides an upper bound for a general unbiased oracle with bounded variance, using a fixed batch size while retaining the single-loop structure. The O(ε−5 ) bound counts individual sample-gradient calls; the iteration count is O(ε−5/2 ). The guarantee concerns the expected squared game-stationarity residual of the original problem.

A variance-reduced SPDE method

We develop a variance-reduced variant of SPDE, termed VR-SPDE, that preserves its single-loop structure. The method maintains one recursive gradient estimator along the entire sequence of current, predictor, and corrected query points. This estimator replaces the three independent mini-batch estimates in Algorithm 1, while the projected steps, normal-cone corrections, dual momentum recursion, and proximal-center update retain the same formulas. Periodic refreshes update only the gradient estimator and require no inner optimization procedure. The recursive estimator introduces a distinct analytical difficulty. In the mini-batch analysis, the fresh endpoint estimate has conditionally mean-zero error, which eliminates the corresponding linear noise terms in (88). A recursive estimate generally retains the error from the preceding query, so this cancellation is no longer available. We instead control the endpoint error pathwise and bound the accumulated estimator variance through the motion of the query points. The resulting variance terms are then absorbed into the descent inequality. This approach yields the variance-reduced GS complexity guarantees in this section; the corresponding OS guarantees are established in Section 5.

4.1

Paired oracle and recursive estimator

Variance reduction requires access to stochastic gradient differences evaluated with a common sample. We impose the following additional oracle assumption.

Assumption 4. In addition to Assumption 2, the oracle permits evaluating ∇f (u; ω) and ∇f (u′ ; ω) using the same sample ω. There exists ℓ ≥ L such that, for every pre-query sigmaalgebra A and every pair of A-measurable feasible points u, u′ ∈ X × Y , a fresh sample ω ∼ D, independent of A, satisfies h

E ∥∇f (u; ω) − ∇f (u′ ; ω)∥2 A ≤ ℓ2 ∥u − u′ ∥2 .

(113)

Samples within each newly drawn batch are independent, and each batch is independent of the information available before it is drawn. Each sample-gradient evaluation at one point counts as one SFO call; evaluating a paired difference therefore costs two SFO calls. Write χ := ℓ/L ≥ 1. Assumption 4 controls stochastic gradient differences in mean square. It is stronger than Lipschitz continuity of the population gradient and is used only for VR-SPDE. It does not require every sample loss to have a uniformly Lipschitz gradient or to be concave in the dual variable. The estimator follows the SPIDER principle [3]. Its use here requires accounting for an adaptive query sequence: each new point may depend on the estimates produced at earlier requests. A recursive estimator along the query sequence. j = 0, . . . , 3T − 1 and define u3t = (xt , yt ),

et , yet ), u3t+1 = (x

Index the estimator requests by

u3t+2 = (xt+1 , yt+1 ),

0 ≤ t < T.

(114)

These points are generated sequentially by the algorithm. For integers q, Br , b ≥ 1, set

gj =

 Br   1 X   ∇f (uj ; ωj,i ),   Br

j≡0

b  1X    g + ∇f (uj ; ωj,i ) − ∇f (uj−1 ; ωj,i ) , j−1  

j ̸≡ 0

(mod q),

i=1

b i=1

(115) (mod q).

At each request, the samples appearing in (115) are drawn afresh after the query points have been determined. Within each paired difference, the two gradients use the same sample. Here Br is the refresh batch size, b is the difference batch size, and q is the refresh period measured in estimator requests, rather than algorithmic iterations. In particular, the initial estimate g0 is always a refresh estimate. The same recursion serves all three query locations; no separate estimator is maintained for the predictor or endpoint. For t ≥ 1, the consecutive requests at the boundary between iterations satisfy u3t = u3t−1 = (xt , yt ). If request 3t is not a refresh, every paired difference in (115) is then identically zero. The update can therefore be implemented exactly by setting g3t = g3t−1 without additional oracle calls. If request 3t is a refresh, the prescribed refresh batch is still drawn. The VR-SPDE algorithm. Algorithm 2 inserts the recursive estimator into the SPDE updates. The parameters q, Br , b govern estimator refresh and accuracy; the parameters τ, h, α, k, β retain their roles in Algorithm 1. Their choices for the two problem settings are specified in the complexity results below. Every occurrence of the second estimate within an iteration uses the same vector g3t+1 . A refresh replaces only the gradient estimate: it does not restart the primal, dual, normal-cone, momentum, or center states. Thus periodic refreshing preserves the single-loop structure. As in the SPDE analysis, the displayed output rule is used for GS; the OS output is specified in Section 5. 26

Algorithm 2 Variance-reduced stochastic projected damped extragradient (VR-SPDE) Input: Deterministic x0 ∈ X and y0 ∈ Y ; integers T, q, Br , b ≥ 1; 0 < τ ≤ L, h > 0, 0 < α ≤ 1, k ≥ 0, and 0 < β ≤ 1. 1: Set z0 = x0 , ξ0 = 0 ∈ Rn , and n0 = v0 = 0 ∈ Rp . 2: for t = 0, . . . , T − 1 do 3: Set u3t = (xt , yt ) and compute g3t by (115). (0) et , yet ) using the predictor updates in Algorithm 1. 4: Set gbt = g3t and compute (x et , yet ) and compute g3t+1 by (115). 5: Set u3t+1 = (x (1) b 6: Set gt = g3t+1 and compute v̄t , xt+1 , ξt+1 , yt+1 , nt+1 using the correction updates in Algorithm 1. 7: Set u3t+2 = (xt+1 , yt+1 ) and compute g3t+2 by (115). (2) 8: Set gbt = g3t+2 and compute vt+1 and zt+1 using the momentum and center updates in Algorithm 1. 9: end for 10: Draw J ∼ Unif{0, . . . , T − 1} independently of all oracle samples. Output: (xout , yout ) = (xJ+1 , yJ+1 ). There are

3T 3T − 1 = Nr := 1 + q q refresh requests. Charging two SFO calls for every sample in every nonrefresh difference gives the total bound NSFO ≤ Br Nr + 2b(3T − Nr ).

This bound counts all refresh and difference evaluations. Copying the estimate at repeated query points can only reduce the oracle cost. Conditional estimator errors. Let Hj denote the sigma-algebra containing all information available immediately before request j is processed. The query point uj is Hj -measurable. For j ≥ 1, the preceding point uj−1 and estimate gj−1 are also Hj -measurable. Define ηj := gj − ∇F (uj ),

at := η3t ,

bt := η3t+1 ,

ct := η3t+2 .

(116)

These errors satisfy the algebraic decompositions in (17). Their conditional moments, however, differ from those of independent mini-batch estimates. At a refresh request, Assumption 2 gives E[ηj | Hj ] = 0,

E[∥ηj ∥2 | Hj ] ≤

σ2 . Br

At a nonrefresh request, write ηj = ηj−1 + δj , where δj :=

b 1X ∇f (uj ; ωj,i ) − ∇f (uj−1 ; ωj,i ) − ∇F (uj ) − ∇F (uj−1 ) . b i=1

Conditional unbiasedness of the paired differences and Assumption 4 imply E[δj | Hj ] = 0,

E[∥δj ∥2 | Hj ] ≤

ℓ2 ∥uj − uj−1 ∥2 . b

Since ηj−1 is Hj -measurable, it follows that E[ηj | Hj ] = ηj−1 ,

E[∥ηj ∥2 | Hj ] ≤ ∥ηj−1 ∥2 + 27

ℓ2 ∥uj − uj−1 ∥2 . b

Thus the recursive increment is conditionally centered, whereas the full estimator error generally is not. Accordingly, the conditional zero-mean identities in (18) are not assumed for at , bt , ct in this section. When a nonrefresh request repeats the preceding point, δj = 0 and ηj = ηj−1 exactly. Feasibility and moment bounds. The projection identities in Lemma 1 apply without change, since VR-SPDE modifies only the gradient estimates. In particular, almost surely, xt , zt ∈ X,

yt ∈ Y,

ξt ∈ NX (xt ),

nt ∈ NY (yt ),

and all predictor points belong to X × Y . All states and estimates also have finite second moments at every finite iteration. Indeed, at a refresh request, population smoothness and the bounded-variance oracle assumption give a finite second moment whenever the query point is square-integrable. At a nonrefresh request, the pairedoracle bound gives a finite second moment for the recursive increment whenever the current and preceding query points are square-integrable. The affine updates and nonexpansive projections preserve this property. Induction along the query sequence, starting from the deterministic initialization, therefore establishes the claim. The analysis below combines a pathwise descent estimate with bounds on query motion and accumulated estimator error. Appropriate choices of the refresh period and batch sizes allow the variance terms to be absorbed into the descent bound, yielding the complexity guarantees for the strongly concave and merely concave settings.

4.2

Convergence analysis

Retain the virtual momentum and state in (71), the quantities It0 , Dt0 from Lemma 6, and + + ⋆ Pt+ , Q+ t , Xt defined in the proof of Lemma 9. Also write Yt = yt+1 − y (zt ). The deterministic error bound (77) and virtual joint descent (80) hold for these correlated errors: their proofs use only the update identities and inequalities that hold pathwise. After inserting (115), the errors at , bt , ct are explicit functions of the realized samples {ωj,i }. The next lemma is therefore proved for a fixed realization of the complete recursive estimator path; expectations enter only in Lemma 13 and the subsequent accumulated bounds. Lemma 11. Suppose that Assumption 1 holds, and let the iterates be generated by Algorithm 2 with the parameter choices in (33) and (79). Set C⋆ := 140. For every realization of the samples {ωj,i }0≤j<3T used in (115), 1 β C⋆ Vt+1 − Vt ≤ − Dt0 − ∥∇p(zt )∥2 + (∥at ∥2 + ∥bt ∥2 + ∥ct ∥2 ). 4 16L Lβ

(117)

Proof. Start from (80) and add the exact endpoint expansion (88). Since the y gradient of G is independent of the center, its first inner product contains Q+ t . By (21), with δt = zt+1 − zt , s

∥yt+1 − y ⋆ (zt+1 )∥ ≤ ∥Yt+ ∥ + sk 16L

2L ∥δt ∥, τ

2L kβ kβ ∥δt ∥∥ct ∥ ≤ ∥Xt+ ∥∥ct ∥ + ∥∇p(zt )∥∥ct ∥. τ 8 16L

(118)

Here the second line uses (83) and s2 = 2Lτ . The dissipation (32) has the lower bound Dt0 ≥

24β + 2 7β 0 2 ∥Qt ∥ + ∥v ∥ + 32τ β∥Yt+ ∥2 + 56Lβ∥Xt+ ∥2 . L L t+1

(119)

For clarity, Young’s inequality bounds each linear endpoint term separately: 2k + 24β + 2 k2 ∥Qt ∥∥ct ∥ ≤ ∥Qt ∥ + ∥ct ∥2 , L 16L Lβ sk 7β 0 2 k2 0 ∥v ∥∥c ∥ ≤ ∥v ∥ + ∥ct ∥2 , t t+1 t+1 128L2 16L Lβ sk k2 ∥Yt+ ∥∥ct ∥ ≤ 2τ β∥Yt+ ∥2 + ∥ct ∥2 , 16L Lβ 7 k2 kβ ∥Xt+ ∥∥ct ∥ ≤ Lβ∥Xt+ ∥2 + ∥ct ∥2 , 8 2 Lβ kβ β k2 ∥∇p(zt )∥∥ct ∥ ≤ ∥∇p(zt )∥2 + ∥ct ∥2 . 16L 16L Lβ

(120)

These bounds follow from ab ≤ ra2 +b2 /(4r), s2 = 2Lτ , s ≤ 2L, and β ≤ 1. The four state-square terms in the first four lines sum to at most Dt0 /16 by (119). The quadratic endpoint term satisfies

2k 2 s 1 2 2 k ∥c ∥ ≤ + ∥ct ∥2 . y,t L 256L2 Lβ

Consequently, (88)–(120) imply 0 Vt+1 − V(zt+1 , wt+1 )≤

1 0 β 7k 2 Dt + ∥∇p(zt )∥2 + ∥ct ∥2 . 16 16L Lβ

(121)

Add (80) and (121). Using Ut2 ≤ 8d2 ∥at ∥2 + 8∥bt ∥2 , d, k, β ≤ 1, and 17 · 8 = 136 < C⋆ , we have 1 1 1 − Dt0 + Dt0 ≤ − Dt0 , 2 16 4 17 2 7k 2 C ⋆ U + ∥ct ∥2 ≤ (∥at ∥2 + ∥bt ∥2 + ∥ct ∥2 ). L t Lβ Lβ The gradient coefficients leave −β/(16L), proving (117). Lemma 12. Suppose that Assumption 1 holds, and let the iterates be generated by Algorithm 2 with the parameter choices in (33) and (79). For each iteration, the query points in (114) satisfy ∥u3t+1 − u3t ∥2 + ∥u3t+2 − u3t+1 ∥2 ≤

Dt0 ∥at ∥2 + ∥bt ∥2 + . Lβ L2

(122)

Proof. Use the fixed-center directions P, Q, P + , Q+ in (30) for the actual initial state and the virtual endpoint q state. Projection nonexpansiveness, the predictor update, and ∥(P, Q) −

(P + , Q+ )∥ =

It0 give

∥u3t+1 − u3t ∥ ≤ h ∥u3t+2 − u3t+1 ∥ ≤ h

q q

2 ∥Pt+ ∥2 + ∥Q+ t ∥ +

It0 + ∥at ∥

It0 + ∥et ∥ + ∥at ∥ .

The second inequality is the projection estimate used in the proof of Lemma 6. By (77), ∥et ∥2 ≤ 8d2 It0 + 2Ut2 ,

Ut2 ≤ 8d2 ∥at ∥2 + 8∥bt ∥2 .

Squaring the two motion bounds and using (a + b + c)2 ≤ 3(a2 + b2 + c2 ) yields ∥u3t+1 − u3t ∥2 + ∥u3t+2 − u3t+1 ∥2 2 2 2 0 2 2 2 2 2 ≤ 3h2 (∥Pt+ ∥2 + ∥Q+ t ∥ ) + (6 + 24d )h It + (6 + 48d )h ∥at ∥ + 48h ∥bt ∥ 2 0 2 2 2 ≤ 8h2 (∥Pt+ ∥2 + ∥Q+ t ∥ + It ) + 64h (∥at ∥ + ∥bt ∥ ).

(123)

The last line uses d = 1/64. From (32), ∥Pt+ ∥2 ≤

4 0 D , h t

2 ∥Q+ t ∥ ≤

L 0 D , 24β t

It0 ≤ 4LDt0 .

Their coefficient in (123) is 32h + h2 L/(3β) + 32h2 L. Since hL ≤ 1/64 and β ≤ 1, !

h2 L Lβ 32h + + 32h2 L 3β

32 32 1 + < 1, + 64 3 · 642 642 1 64h2 ≤ 2 . L ≤

Substitution proves (122). Define the cumulative expected error and dissipation by ET :=

3T −1 X j=0

E∥ηj ∥2 =

TX −1

E(∥at ∥2 + ∥bt ∥2 + ∥ct ∥2 ),

AT :=

t=0

TX −1

EDt0 .

(124)

t=0

Lemma 13. Suppose that Assumptions 1, 2, and 4 hold, and let the iterates and estimators be generated by Algorithm 2 with the parameter choices in (33) and (79). Then ET ≤

3T σ 2 qχ2 L qχ2 + AT + ET . Br bβ b

(125)

Proof. At a refresh index j, (14) gives E[∥ηj ∥2 | Hj ] ≤ σ 2 /Br . At a nonrefresh index, define the centered fresh increment ζj :=

b 1X [∇f (uj ; ωj,i ) − ∇f (uj−1 ; ωj,i )] − [∇F (uj ) − ∇F (uj−1 )]. b i=1

Both points and ηj−1 are Hj -measurable. Conditional independence and (113) give ηj = ηj−1 + ζj ,

E[ζj | Hj ] = 0,

E[∥ηj ∥2 | Hj ] = ∥ηj−1 ∥2 + E[∥ζj ∥2 | Hj ] ≤ ∥ηj−1 ∥2 +

ℓ2 ∥uj − uj−1 ∥2 . b

(126)

Let r(j) := q⌊j/q⌋ be the last refresh index. Iterating (126) after that refresh yields j

E∥ηj ∥2 ≤

σ2 ℓ2 X E∥ui − ui−1 ∥2 . + Br b i=r(j)+1

(127)

When (127) is summed over j < 3T , each motion term is counted at most q times. Moreover, u3t = u3t−1 for t ≥ 1. Thus ET ≤

−1 X 3T σ 2 qℓ2 3T + E∥uj − uj−1 ∥2 Br b j=1

3T σ 2 qℓ2 AT ET ≤ + + 2 , Br b Lβ L

where the last line uses (122) and nonnegativity of the endpoint-error squares. Since χ = ℓ/L, this is (125).

Lemma 14. Suppose that Assumptions 1, 2, and 4 hold, and let the iterates be generated by Algorithm 2 with the parameter choices in (33) and (79). Suppose further that b≥

16C⋆ qχ2 , β2

MTvr := B +

6C⋆ T σ 2 . LβBr

(128)

Then Algorithm 2 satisfies −1 1 β TX AT + E∥∇p(zt )∥2 ≤ MTvr , 8 16L t=0

(129)

and its output obeys ER(xJ+1 , yJ+1 )2 ≤

128LB 768C⋆ σ 2 + 2τ 2 DY2 . + βT β 2 Br

(130)

Proof. Since qχ2 /b ≤ β 2 /(16C⋆ ) ≤ 1/2, rearranging (125) gives ET ≤

6T σ 2 2qχ2 L + AT . Br bβ

(131)

Sum (117) in expectation and use VT ≥ 0 and V0 ≤ B to obtain 1 β X C⋆ AT + E∥∇p(zt )∥2 ≤ B + ET 4 16L t<T Lβ 6C⋆ T σ 2 2C⋆ qχ2 + AT LβBr bβ 2 1 ≤ MTvr + AT . 8

≤B+

Moving the last term to the left proves (129). Retaining individual components of (32) gives X t<T

E∥Pt+ ∥2 ≤

32MTvr , h

E∥Xt+ ∥2 ≤

t<T

2 E∥Q+ t ∥ ≤

t<T MTvr

7Lβ

LMTvr , 3β

0 E∥vt+1 ∥2 ≤

t<T

E∥∇p(zt )∥2 ≤

t<T

16LMTvr . β

8LMTvr , 7β (132)

The algebraic certificate St in (96) satisfies the same pointwise bound used to prove (97). Hence X

ESt ≤

t<T

2 16 96β 12 + + 48 + + hL 7 3 7

LMTvr 64LMTvr ≤ . β β

(133)

The last inequality uses β/(hL) ≤ 1/2048. Finally, R2 ≤ 2St + 2τ 2 DY2 and independent uniform selection of J give 128LMTvr ER(xJ+1 , yJ+1 )2 ≤ + 2τ 2 DY2 . βT Substitute (128) to obtain (130).

4.3

Complexity results

The preceding lemmas reduce the remaining work to selecting the refresh batch, difference batch, and refresh period. We first exploit strong concavity, where no bias term is present, and then impose the same accuracy-dependent dual regularization used in the ordinary concave result.

4.3.1

Nonconvex–strongly concave setting

VR-SPDE applies Algorithm 2 to the query stream generated by SPDE: every explicit dualregularization term is deleted and the parameters in (100) are used. Lemma 15. Suppose that Assumptions 1, 2, 4, and 3 hold and that b≥

16C⋆ qχ2 . 2 βsc

(134)

Then VR-SPDE satisfies ER(xJ+1 , yJ+1 )2 ≤

128LB 768C⋆ σ 2 . + 2B βsc T βsc r

(135)

Proof. Apply the strongly-concave substitution established in the proof of Lemma 10 to Lemmas 11–14. The paired-oracle recursion and query-motion estimate use only mean-square smoothness, projection nonexpansiveness, and the scalar relations displayed after (100); hence their constants are unchanged. The initial potential remains bounded by B, and the exact original dual certificate again removes the regularization-bias term. Substituting βsc into (130) therefore gives (135). Theorem 4. Suppose that Assumptions 1, 2, 4, and 3 hold. Fix a known bound B ≥ B and set (

256LB T = max 1, βsc ε2

(

1536C⋆ σ 2 Br = max 1, 2 ε2 βsc

, (

16C⋆ χ2 , Avr,sc := 2 βsc

q = max 1,

Br Avr,sc

(136)

b = ⌈Avr,sc q⌉.

Then VR-SPDE returns a point satisfying (7), and n

Nvr,sc ≤ Br + 18T min Br ,

Avr,sc Br .

For fixed positive σ and χ, the parameter and oracle orders are √ T = O( κε−2 ), Br = O(κε−2 ), q = Θ((χε)−1 ),

(137)

b = O(χκε−1 ),

(138)

Nvr,sc = O κε−2 + χκ3/2 ε−3 . √ If σ = 0, then Br = q = 1 and Nvr,sc ≤ 3T = O( κε−2 ). 2 , so Lemma 15 applies. Moreover, Proof. The batch choice gives b ≥ 16C⋆ qχ2 /βsc

128LB ε2 ≤ , βsc T 2

768C⋆ σ 2 ε2 ≤ , 2B βsc 2 r

which proves (7). There are 3T estimator requests and exactly ⌈3T /q⌉ refreshes. If q ≥ 2, then Br ≥ 4Avr,sc and s q 1 Br ≤ q, b ≤ 2 Avr,sc Br . 2 Avr,sc Charging Br calls per refresh and 2b calls per paired difference gives Nvr,sc ≤ Br +

q 3T Br + 6T b ≤ Br + 18T Avr,sc Br . q

If q = 1, p every request is a refresh and Nvr,sc = 3T Br . In this case Br < 4Avr,sc , so Br ≤ 2 min{Br , Avr,sc Br } and the same bound follows. This proves (137) in both cases. 32

By (100), βsc = Θ(κ−1/2 ). For fixed positive σ and χ, (136) therefore gives Br = Θ(κε−2 ),

Br /Avr,sc = Θ((χε)−1 ).

Avr,sc = Θ(χ2 κ),

Consequently, q = Θ((χε)−1 ), b = O(χκε−1 ), and Nvr,sc ≤ Br + 18T

Avr,sc Br

= O κε−2 + χκ3/2 ε−3 . If σ = 0, then Br = 1 < Avr,sc , hence q = 1; all 3T requests are refreshes of size one and Nvr,sc ≤ 3T . Compared with SPDE, the iteration count is unchanged. Variance reduction replaces a fresh batch at every request by periodic refreshes and paired differences, yielding O(κε−2 + χκ3/2 ε−3 ) oracle calls under the stronger paired mean-square smoothness assumption. 4.3.2

Nonconvex–concave setting

Theorem 5. Suppose that Assumptions 1, 2, and 4 hold. Fix a known bound B ≥ B independent of ε, τ . For ε > 0, choose τ as in (105), choose h, α, k by (33), and choose β by (79). Set (

512LB T = max 1, βε2

(

3072C⋆ σ 2 Br = max 1, β 2 ε2

, (

16C⋆ χ2 , Avr := β2

q = max 1,

Br Avr

, (139)

b = ⌈Avr q⌉.

Then Algorithm 2 satisfies (7). Its number of stochastic full-gradient oracle calls obeys

Nvr ≤ Br

3T q

+ 2b 3T −

≤ Br + 18T min{Br ,

3T q

Avr Br }.

(140)

For fixed problem data, B, and σ > 0, as ε ↓ 0, T = O(ε−5/2 ),

Br = O(ε−3 ),

q = Θ(ε−1 ),

b = O(ε−2 ),

Nvr = O(ε−9/2 ). (141)

If σ = 0, the choices in (139) give Br = q = 1 and Nvr ≤ 3T = O(ε−5/2 ). Proof. By (139), b ≥ 16C⋆ qχ2 /β 2 , so Lemma 14 applies. Its three residual terms satisfy 128LB ε2 ≤ , βT 4

768C⋆ σ 2 ε2 ≤ , β 2 Br 4

2τ 2 DY2 ≤

ε2 . 2

(142)

When σ = 0, the second term is zero. Adding (142) proves ER2 ≤ ε2 . There are 3T estimator requests, including exactly ⌈3T /q⌉ refresh requests. Each refresh costs Br calls and each other request costs at most 2b calls, proving the first line of (140). Since Avr ≥ 1, if q ≥ 2 then Br ≥ 4Avr and 1 2

Br ≤q≤ Avr

Br , Avr

Avr Br = Br + 18T min{Br ,

b≤

Avr Br + 1 ≤ 2 Avr Br .

Consequently, Nvr ≤ Br +

3T Br + 6T b q

≤ Br + 18T

Avr Br }.

If q = 1, √ every request is a refresh and Nvr = 3T Br . Here Br < 4Avr , which implies Br ≤ 2 min{Br , Avr Br }; the same upper bound follows. For fixed σ > 0 and sufficiently small ε, τ = ε/(2DY ) and (111) gives β = Θ(ε1/2 ). Thus (139) implies q Br = Θ(ε−3 ), Avr = Θ(ε−1 ), Br /Avr = Θ(ε−1 ). It follows that q = Θ(ε−1 ), b = O(ε−2 ), and

Nvr = O Br + T

Avr Br = O ε−3 + ε−5/2 ε−2 = O(ε−9/2 ).

When σ = 0, Br = 1 and Avr > 1 force q = 1, so the difference batches are never used and Nvr ≤ 3T . Remark 2. The iteration order remains O(ε−5/2 ) in both stochastic methods. In SPDE, every iteration pays for three batches of order ε−5/2 . In VR-SPDE, a refresh of order ε−3 is spread over Θ(ε−1 ) requests, and the difference batches have order ε−2 . The average oracle work per iteration is therefore O(ε−2 ) for fixed positive variance. The stronger oracle in Assumption 4 enables this improvement. The bounds in Theorems 3 and 5 are upper bounds for the same expected squared game-stationarity criterion.

Optimization-stationarity guarantees

The preceding sections use game stationarity, which evaluates a feasible primal–dual pair. We now analyze the optimization-stationarity criterion introduced in (10). The same centerbased algorithmic constructions, with the OS-specific parameter choices given below, provide this guarantee through the center sequences of SPDE and VR-SPDE without any additional stochastic-gradient evaluation. For 0 ≤ τ ≤ L, define the regularized value function and its 1/(2L)-Moreau envelope over X by n o τ pτ (z) := min ϕτ (x) + L∥x − z∥2 . (143) ϕτ (x) := max F (x, y) − ∥y∥2 , x∈X y∈Y 2 Thus ϕ0 = ϕ, p0 agrees with the original envelope in (8), pτ = p in the regularized analysis of Section 3, and p0 = psc in the strongly-concave specialization. Assumption 1 implies that every ϕτ is L-weakly convex. Consequently, the minimizer n

x⋆τ (z) := argmin ϕτ (x) + L∥x − z∥2

(144)

x∈X

is unique and ∇pτ (z) = 2L(z − x⋆τ (z)).

(145)

For the OS results below, the algorithms use the same independent index J ∼ Unif{0, . . . , T − 1} as before and additionally report zout := zJ . Storing and reporting this center changes neither the updates nor the SFO count. Lemma 16. Suppose that Assumption 1 holds. For every 0 ≤ τ ≤ L and z ∈ Rn , 0 ≤ ϕ(x) − ϕτ (x) ≤

τ 2 D , 2 Y

τ DY2 , L √ ∥∇p0 (z) − ∇pτ (z)∥ ≤ 2DY Lτ . ∥x⋆τ (z) − x⋆0 (z)∥2 ≤

x ∈ X,

(146) (147) (148)

In particular, SOS (z)2 ≤ 2∥∇pτ (z)∥2 + 8Lτ DY2 . 34

(149)

Proof. For every x ∈ X, subtracting the nonnegative quadratic from the inner maximization gives ϕτ (x) ≤ ϕ(x). If y0 (x) ∈ argmaxy∈Y F (x, y), compactness of Y and the definition of DY give τ τ ϕτ (x) ≥ F (x, y0 (x)) − ∥y0 (x)∥2 ≥ ϕ(x) − DY2 . 2 2 This proves (146). The lower smoothness inequality for F (·, y) shows that x 7→ F (x, y) + (L/2)∥x∥2 is convex for every y ∈ Y . Taking the pointwise maximum over y preserves convexity, and hence every ϕτ is L-weakly convex. Moreover, ϕτ ≥ ϕinf − τ DY2 /2, so the proximal objectives below are coercive on the closed set X and attain their minima. For fixed z, set τ Pτ,z (x) := ϕτ (x) + L∥x − z∥2 , δτ := DY2 . 2 Since ϕτ is L-weakly convex, Pτ,z is L-strongly convex on X. Applying strong convexity of P0,z at its minimizer and then (146) gives L ⋆ ∥x (z) − x⋆0 (z)∥2 ≤ P0,z (x⋆τ (z)) − P0,z (x⋆0 (z)) 2 τ ≤ Pτ,z (x⋆τ (z)) + δτ − P0,z (x⋆0 (z)) ≤ Pτ,z (x⋆0 (z)) + δτ − P0,z (x⋆0 (z)) ≤ δτ .

(150)

Rearranging (150) proves (147). Using (145), ∥∇p0 (z) − ∇pτ (z)∥ = 2L∥x⋆τ (z) − x⋆0 (z)∥ ≤ 2DY

√ Lτ ,

which is (148). Finally, ∥a + b∥2 ≤ 2∥a∥2 + 2∥b∥2 gives (149). The next lemma extracts the envelope-gradient estimates already contained in the two descent analyses. Lemma 17. Suppose that Assumptions 1 and 2 hold, and let J be independent and uniform on {0, . . . , T − 1}. For SPDE with 0 < τ ≤ L and the parameters in (33) and (79), E∥∇pτ (zJ )∥2 ≤

8LB 2400σ 2 + . βT βm

(151)

If Assumption 4 also holds and b ≥ 16C⋆ qχ2 /β 2 , then VR-SPDE satisfies E∥∇pτ (zJ )∥2 ≤

16LB 96C⋆ σ 2 + 2 . βT β Br

(152)

Under Assumption 3, the corresponding unregularized SPDE and VR-SPDE estimates are obtained from (151) and (152), respectively, by replacing (pτ , β) with (p0 , βsc ). Proof. Independence and uniformity of J give −1 1 TX E∥∇pτ (zJ )∥ = E∥∇pτ (zt )∥2 . T t=0 2

For SPDE, the last estimate in (95) and the definition (92) yield 8L E∥∇pτ (zJ )∥ ≤ βT 2

300T σ 2 B+ Lm 35

8LB 2400σ 2 + . βT βm

For VR-SPDE, the last estimate in (132) and (128) give 16L 6C⋆ T σ 2 E∥∇pτ (zJ )∥ ≤ B+ βT LβBr

16LB 96C⋆ σ 2 . + 2 βT β Br

The strongly-concave drift used in Lemma 10 retains the same envelope term with (p0 , βsc ). Likewise, the strongly-concave substitution in Lemma 15 retains the corresponding VR envelope term. This proves the last assertion.

5.1

Nonconvex–strongly concave setting

Strong concavity removes the artificial dual regularization, so the envelope controlled by the descent inequality is already the envelope of the original value function. Theorem 6. Suppose that Assumptions 1, 2, and 3 hold, with the strong-concavity modulus chosen so that 0 < µ ≤ L; any larger valid modulus may be replaced by min{µ, L}. Let ε > 0 and fix a known bound B ≥ B. For SPDE, use the unregularized specialization and the parameters in (100), and set (

16LB T = max 1, βsc ε2

(

4800σ 2 m = max 1, βsc ε2

(153)

Then ESOS (zout )2 ≤ ε2 , and the number of SFO calls satisfies 16LB Nos,sc ≤ 3 1 + βsc ε2

4800σ 2 1+ βsc ε2

= O(κε−4 )

(154)

for fixed positive L, B, and σ. If Assumption 4 also holds, run the unregularized VR-SPDE specialization with (

32LB T = max 1, βsc ε2 16C⋆ χ2 Aos,sc := , 2 βsc

(

192C⋆ σ 2 Br = max 1, 2 ε2 βsc

, (

q = max 1,

Br Aos,sc

(155)

b = ⌈Aos,sc q⌉ .

Then ESOS (zout )2 ≤ ε2 and Nos,vr,sc ≤ Br + 18T min Br ,

= O κε−2 + χκ3/2 ε−3

Aos,sc Br

(156)

√ for fixed positive σ and χ. If σ = 0, both methods use O( κε−2 ) SFO calls. Proof. For SPDE, the strongly-concave form of (151) and (153) give ESOS (zout )2 ≤

8LB 2400σ 2 ε2 ε2 + ≤ + = ε2 . βsc T βsc m 2 2

The three fresh batches per iteration give the first inequality in (154). Since βsc = Θ(κ−1/2 ), √ both T and m are O( κε−2 ), proving its stated order. 2 . The strongly-concave form of (152) therefore For VR-SPDE, (155) implies b ≥ 16C⋆ qχ2 /βsc gives ESOS (zout )2 ≤

16LB 96C⋆ σ 2 ε2 ε2 + 2 ≤ + = ε2 . βsc T βsc Br 2 2 36

The refresh and paired-difference accounting used in (137) applies without change and proves the first line of (156). Moreover, √ T = O( κε−2 ), Br = O(κε−2 ), Aos,sc = O(χ2 κ), q = Θ((χε)−1 ),

b = O(χκε−1 ),

Aos,sc Br = O(χκε−1 ).

Substitution into the query bound proves the second line of (156). When σ = 0, the mini-batch √ method uses m = 1, while the VR method has Br = q = 1; both costs are O(T ) = O( κε−2 ).

5.2

Nonconvex–concave setting

In the merely concave setting, the regularized envelope pτ differs from the envelope p0 of the original value function. The square-root transfer error in (148) requires a smaller regularization parameter than the one used for game stationarity. Theorem 7. Suppose that Assumptions 1 and 2 hold. Let ε > 0, fix a known bound B ≥ B independent of ε and τ , and set (

ε2 τos := min L, 16LDY2

)

4LDY ε

19200σ 2 m = max 1, βε2

Kεos := max 1,

(157)

Choose h, α, k by (33) and β by (79), with τ = τos . For SPDE, set (

64LB T = max 1, βε2

(

(158)

Then ESOS (zout )2 ≤ ε2 and 64LB Nos,ncc ≤ 3 1 + βε2

19200σ 2 1+ βε2

(LB + σ 2 )Kεos LBσ 2 (Kεos )2 =O 1+ + ε2 ε4

(159)

For fixed problem data, B, and σ > 0, this gives T = m = O(ε−3 ) and Nos,ncc = O(ε−6 ). If Assumption 4 also holds, run VR-SPDE with (

128LB T = max 1, βε2 16C⋆ χ2 Aos,ncc := , β2

(

768C⋆ σ 2 Br = max 1, β 2 ε2

, (

q = max 1,

, (160)

Br Aos,ncc

b = ⌈Aos,ncc q⌉ .

Then ESOS (zout )2 ≤ ε2 and n

Nos,vr,ncc ≤ Br + 18T min Br ,

Aos,ncc Br

= O(ε−4 + χε−6 ) = O(ε−6 )

(161)

for fixed positive problem data, σ, and χ. If σ = 0, both methods use O(ε−3 ) SFO calls. Proof. The choice (157) gives 8Lτos DY2 ≤

ε2 . 2

(162)

For SPDE, (151) and (158) imply E∥∇pτos (zJ )∥2 ≤

8LB 2400σ 2 ε2 ε 2 ε2 + ≤ + = . βT βm 8 8 4

(163)

Taking expectations in (149) and using (162)–(163) gives ESOS (zout )2 ≤ 2E∥∇pτos (zJ )∥2 + 8Lτos DY2 ≤ ε2 . The three batches per iteration give the first line of (159). By (111) and (157), s

1 =Θ β

L τos

L = Kεos . τos

(164)

Expanding the first line and using (164) proves the second line of (159). For sufficiently small ε, Kεos = Θ(ε−1 ), proving the stated orders. For VR-SPDE, (160) ensures b ≥ 16C⋆ qχ2 /β 2 . Hence (152) gives E∥∇pτos (zJ )∥2 ≤

16LB 96C⋆ σ 2 ε2 ε2 ε2 + 2 + = . ≤ βT β Br 8 8 4

(165)

Equations (149), (162), and (165) prove ESOS (zout )2 ≤ ε2 . The query accounting in (140) is independent of the choice of τ and proves the first line of (161). For fixed positive σ and χ, (164) and (160) give T = O(ε−3 ),

Br = O(ε−4 ),

q = Θ((χε)−1 ),

b = O(χε−3 ),

Aos,ncc = O(χ2 ε−2 ), q

Aos,ncc Br = O(χε−3 ).

Substitution gives Nos,vr,ncc = O(ε−4 ) + O(ε−3 )O(χε−3 ) = O(ε−4 + χε−6 ). When σ = 0, SPDE uses m = 1, while VR-SPDE has Br = q = 1. In both cases the SFO count is O(T ) = O(ε−3 ). The NC–C OS and GS parameter choices differ for a structural reason. Game stationarity incurs the direct dual bias O(τ ) and therefore permits τ = Θ(ε). Optimization stationarity √ transfers a Moreau gradient, whose norm bias is O( τ ), and consequently requires τ = Θ(ε2 ). This reduces the center rate from Θ(ε1/2 ) to Θ(ε) and yields the O(ε−6 ) SFO order for both stochastic estimators under the present analysis.

Conclusions

We developed single-loop stochastic projected damped extragradient methods for smooth nonconvex–(strongly) concave minimax optimization, with complexity guarantees for both game stationarity and optimization stationarity. The proposed SPDE method combines predictor– corrector projections, damped dual momentum, normal-cone corrections, and a relaxed proximalcenter update. Its variance-reduced variant, VR-SPDE, maintains a recursive gradient estimator along the same sequence of query points. Both methods update the primal, dual, momentum, normal-cone, and center states within a single loop, without inner iterative solvers for regularized subproblems. For game stationarity, SPDE returns a feasible pair satisfying E[R(xout , yout )2 ] ≤ ε2 with total stochastic first-order oracle (SFO) complexities of O(κε−4 ) and O(ε−5 ) in the nonconvex– strongly concave and nonconvex–concave settings, respectively, where κ = L/µ. These guarantees 38

require an unbiased stochastic gradient oracle with uniformly bounded variance. With paired sample evaluations and an additional mean-square Lipschitz condition on the stochastic gradients, VR-SPDE improves the corresponding GS complexities to O(κ3/2 ε−3 ) and O(ε−9/2 ). For optimization stationarity, the same algorithms, with criterion-specific parameter choices, return a proximal center zout = zJ satisfying E[SOS (zout )2 ] ≤ ε2 , where SOS is defined through the Moreau envelope of the constrained primal value function of the original problem. SPDE achieves SFO complexities of O(κε−4 ) and O(ε−6 ) in the nonconvex–strongly concave and nonconvex– concave settings, respectively, while VR-SPDE achieves O(κ3/2 ε−3 ) and O(ε−6 ). Thus variance reduction improves the stated OS bound in the strongly concave setting, while both methods attain the same O(ε−6 ) bound in the merely concave setting. To the best of our knowledge, these results provide the best-known SFO complexity guarantees among single-loop stochastic first-order methods for the respective stationarity criteria and problem classes. The OS guarantees also match the best-known multi-loop dependence on the target accuracy in both settings. Together, these results establish that a single-loop damped extragradient framework can support improved GS guarantees and competitive OS guarantees under their respective oracle assumptions.

References [1] J. Chen and V. K. N. Lau, Convergence analysis of saddle point problems in time varying wireless systems—control theoretical approach, IEEE Transactions on Signal Processing, 60(1) (2012), pp. 443–452. https://doi.org/10.1109/TSP.2011.2169407. [2] Y. Fan, S. Lyu, Y. Ying, and B.-G. Hu, Learning with average top-k loss, Advances in Neural Information Processing Systems, 30 (2017), pp. 497–505. https://proceedings. neurips.cc/paper/2017/hash/6c524f9d5d7027454a783c841250ba71-Abstract.html. [3] C. Fang, C. J. Li, Z. Lin, and T. Zhang, SPIDER: Near-optimal non-convex optimization via stochastic path-integrated differential estimator, Advances in Neural Information Processing Systems, 31 (2018), pp. 689–699. https://proceedings.neurips.cc/paper/2018/hash/ 1543843a4723ed2ab08e18053ae6dc5b-Abstract.html. [4] G. B. Giannakis, Q. Ling, G. Mateos, I. D. Schizas, and H. Zhu, Decentralized learning for wireless communications and networking, in Splitting Methods in Communication, Imaging, Science, and Engineering, R. Glowinski, S. J. Osher, and W. Yin (eds.), Springer, Cham, 2016, pp. 461–497. https://doi.org/10.1007/978-3-319-41589-5_14. [5] F. Huang, S. Gao, J. Pei, and H. Huang, Accelerated zeroth-order and first-order momentum methods from mini to minimax optimization, Journal of Machine Learning Research, 23(36) (2022), pp. 1–70. https://jmlr.org/papers/v23/20-924.html. [6] X. Jiang, L. Zhu, T. Zheng, and A. M.-C. So, Efficient single-loop stochastic algorithms for nonconvex-concave minimax optimization, arXiv preprint arXiv:2501.05677v2 (2025). https://arxiv.org/abs/2501.05677v2. [7] T. Lin, C. Jin, and M. I. Jordan, On gradient descent ascent for nonconvex-concave minimax problems, in Proceedings of the 37th International Conference on Machine Learning, Proceedings of Machine Learning Research, 119 (2020), pp. 6083–6093. https: //proceedings.mlr.press/v119/lin20a.html. [8] L. Luo, H. Ye, Z. Huang, and T. Zhang, Stochastic recursive gradient descent ascent for stochastic nonconvex-strongly-concave minimax problems, Advances in Neural Information Processing Systems, 33 (2020), pp. 20566–20577. https://proceedings.neurips.cc/ paper/2020/hash/ecb47fbb07a752413640f82a945530f8-Abstract.html. 39

[9] L. M. Nguyen, J. Liu, K. Scheinberg, and M. Takáč, SARAH: A novel method for machine learning problems using stochastic recursive gradient, in Proceedings of the 34th International Conference on Machine Learning, Proceedings of Machine Learning Research, 70 (2017), pp. 2613–2621. https://proceedings.mlr.press/v70/nguyen17b.html. [10] H. Rafique, M. Liu, Q. Lin, and T. Yang, Weakly-convex–concave min–max optimization: Provable algorithms and applications in machine learning, Optimization Methods and Software, 37(3) (2022), pp. 1087–1121. https://doi.org/10.1080/10556788.2021.1895152. [11] S. Sagawa, P. W. Koh, T. B. Hashimoto, and P. Liang, Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization, in International Conference on Learning Representations, 2020. https://openreview.net/ forum?id=ryxGuJrFvS. [12] S. Shafieezadeh-Abadeh, P. Mohajerin Esfahani, and D. Kuhn, Distributionally robust logistic regression, Advances in Neural Information Processing Systems, 28 (2015), pp. 1576–1584. https://proceedings.neurips.cc/paper/2015/hash/ cc1aa436277138f61cda703991069eaf-Abstract.html. [13] J. Yang, A. Orvieto, A. Lucchi, and N. He, Faster single-loop algorithms for minimax optimization without strong concavity, in Proceedings of the 25th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, 151 (2022), pp. 5485–5517. https://proceedings.mlr.press/v151/yang22b.html. [14] H. Zhang and Z. Xu, An accelerated first-order regularized momentum descent ascent algorithm for stochastic nonconvex-concave minimax problems, Computational Optimization and Applications, 90 (2025), pp. 557–582. https://doi.org/10.1007/s10589-024-00638-9. [15] X. Zhang, N. S. Aybat, and M. Gürbüzbalaban, SAPD+: An accelerated stochastic method for nonconvex-concave minimax problems, Advances in Neural Information Processing Systems, 35 (2022), pp. 21668–21681. https://proceedings.neurips.cc/paper_files/ paper/2022/hash/880d8999c07a8efc9bbbeb0c38f50765-Abstract-Conference.html.

Record · ID 1006870 · SHA-256 271bab7e629de61e
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.