Variance-Reduced Forward-Reflected-Backward Splitting Methods for Nonmonotone Generalized Equations
Authors: Quoc Tran-Dinh
Venue: ICML 2025 poster (2025)
Keywords: variance-reduction; forward-reflected-backward splitting method; nonmonotone generalized equations; SVRG; SAGA
Primary area: optimization->stochastic
TL;DR: We develop two novel stochastic variance-reduced optimistic gradient-type methods to solve a class of nonmontone generalized equations.

Abstract: We develop two novel stochastic variance-reduction methods to approximate solutions of a class of nonmonotone [generalized] equations. Our algorithms leverage a new combination of ideas from the forward-reflected-backward splitting method and a class of unbiased variance-reduced estimators. We construct two new stochastic estimators within this class, inspired by the well-known SVRG and SAGA estimators. These estimators significantly differ from existing approaches used in minimax and variational inequality problems. By appropriately choosing parameters, both algorithms achieve state-of-the-art oracle complexity of $\mathcal{O}(n + n^{2/3} \epsilon^{-2})$ for obtaining an $\epsilon$-solution in terms of the operator residual norm for a class of nonmonotone problems, where $n$ is the number of  summands and $\epsilon$ signifies the desired accuracy. This complexity aligns with the best-known results in SVRG and SAGA methods for stochastic nonconvex optimization.  We test our algorithms on some numerical examples and compare them with existing methods. The results demonstrate promising improvements offered by the new methods compared to their competitors.

Variance-Reduced ForwarddReefectee--Bcckww Splitting Methods for Nonmonotone Generalized Equations
Quoc Tran-Dinh
Abstract
ing and robust optimization has brought renewed interest to generalized equations and their special case: minimax problem. They serve as powerful tools for handling Nash's equilibria and minimax models in generative adversarial nets, adversarial training, and robust learning, see (Arjovsky et al., 2017: Goodfellow al., 2014; Madry et al., 201 8; Namkoong & Duchi. 2016). Notably, most problems arising from these applications are nonmonotone, nonsmooth, and large-scale. This paper develops new and simple stochastic algorithms with variance reduction for solving this class of problems, equipped with rigorous theoretical guarantees.
We develop two novel stochastic variance- reduction methods to approximate solutions of a class of nonmonotone [generalized] equations. Our algorithms leverage a new combination of ideas from the sewwardreflected--accwaa split- ting method and class of unbiased variance- reduced estimators. We construct two new stochastic estimators within this class, inspired by the well-known SVRG and SAGA estimators. These estimators significantly differ from existing approaches used in minimax and variational in- equality problems. By appropriately choosing pa- rameters, both algorithms achieve state-of-the-ar oracle complexity of O(n + n n=2 ing an e-solution in terms of the operator residual norm for a class of nonmonotone problems, where n is the number of summands and signifies the desired accuracy. This complexity aligns with the best-known results in SVRG and SAGA methods for stochastic nonconvex optimization. We test algorithms some numerical examples and compare them with existing methods. The results demonstrate promising improvements offered by the new methods compared to their competitors.
1.1. Nonmonotone finite-sum generalized equations The central problem we study in this paper the following [possibly nonmonotone] generalized equation (also known composite inclusion) (Rockafellar & Wets, 1997):
?) for obtain-
Find x² E RP such that: 0 Gx - Tx" (NI)
where G TRP - RP is a single-valued operator, possibly nonlinear, and T R = 2 is multivalued mapping from E' to 2d (the set of all subsets of R?).. In addition, we assume that G is given in the following large finite-sum:
Gx ==
G.
(I)
where G, : R + Rp are given operators for all i [n]] {1, .n] and n >>>>>> 1. This structure often arises from statistical learning, [generative] machine learning, networks, distributed systems, and data science. For simplicity of notation, we denote V = G+T T and dom (V) := dom(G)) dom(T), (T). where dom R) is the domain of R.
1. Introduction
[Non]linear equations and inclusions comerstones of computational mathematics, finding applications in diverse fields like engineering, mechanics, economics, statistics, optimization, and machine learning, see, e.g., (Bauschke Combettes, 2017; Burachik & lusem, 2008; Facchinei & Pang, 2003; Phelps, 2009; Ryu Yin, 2022; Ryu Boyd, 2016). These problems, known as generalized equations (Rockafellar & Wets, 1997), are equivalent to fixed-point problems. The recent revolution in modern machine learn-
We highlight that the methods developed in this paper be straightforwardly extended to tackle Gr Eş G(z, E) as the expectation of a stochastic operator G involving random vector & defined on a probability space (52, IP, 2).
1.2. Equivalent forms and special cases
DDepartment of Statistics and Operations Research, The Univer- sity North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA. Correspondence to: Quoc Tran- <quoetd@ email.unc.co und.
The model (NI) covers many fundamental problems in opti- mization and related fields, including the following ones.
(a) [Non]linear equation. If 0. then (NI) reduces to the following [non]linear equation:
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
Find x² E dom(G) such that: Gz 0.
(NE)

---

Variance-Reduced
Splitting Methods
Both (NT) and (NE) are also called root-finding problems. Clearly, (NE) is special case of (ND). However, un-
(ii) Optimality certification. Existing stochastic methods often target special cases of (NI) such (NE) and (VIP). In addition, these methods frequently rely monotonic- ity assumption, which excludes many problems of current interest, e.g., (Alacaoglu et al., 2023; Alacaoglu & Malit- sky, 2022; Beznosikov et al., 2023; Gorbunov et al., 2022a; Loizou et al., 2021). Furthermore, existing methods analyze convergence based on a [duality] gap function (Facchinei & Pang. 2003) restricted gap function (Nesterov, 2007). As discussed in (Cai et al., 2024; Diakonikolas, 2020), these metrics have limitations, particularly in nonmonotone set- tings. It is important to note that standard gap functions are not applicable to our settings due to Assumption 1.4. Re- garding oracle complexity, several works, e.g., (Alacaoglu & Malitsky, 2022: Beznosikov et al., 2023; Gorbunov et al., 2022a; Loizou et al., 2021) claim an oracle complexity of (n + vnc to attain e-solution, but this is measured using a restricted gap function. Again, as highlighted in (Cai et al., 2024; Diakonikolas, 2020), this certification does not translate to the operator residual norm and is inapplicable to nonmonotone settings. Therefore, a direct comparison between our results and these previous works is challenging due these methodological discrepancies (see Table 1).
der appropriate assumptions on G and/or T (e.g., using the resolvent of T), one can also transform (NI) to (NE). Let ver(() (V) : [x* E dom( V)= 0 E Vx*"} } and zer(G) == fo* E dom(G) : Gx* = 0} be the solution sets of (NI) and (NE), respectively, which are assumed to be nonempty.
(b) Variational inequality problem (VIP). Nxx, the normal cone of a nonempty, closed, and convex set X' in R", then (NI) reduces to the following VIP:
Find x² E X such that: (Gz" x - x²) 2 O, VI E X. (VIP)
If T = åg, the subdifferential of a convex function g. then (NI) reduces to mixed VIP, denoted by MVIP. Both VIP and MVIP cover many problems in practice, including min- imax problems and Nash's equilibria, see, e.g. (Burachik & Iusem, 2008; Facchinei & Pang, 2003; Phelps, 2009).
(c) Minimax optimization. Another important special case of (NI) (or MVIP) is the following minimax optimization (or (ondle-point) problem, which has found various applications in machine learning and robust optimization:
{ [C(u. v) == p(u)+ + H(u, v) - v(v)). (SP)
min max KERPI VERRS
(iii) New and simple algorithus.. Various existing stochas- tic methods for solving (VIP) and (NI) rely on estab- lished techniques. These include mirror-prox/averaging and extragradicnt-type schemes combined with the clas- sie Robbin-Monn stochastic approximation (Robbins & Monro, 1951) (e-g.. (Cui & Shanbhag, 2021; Iusem et al., 2017: Juditsky et al., 2011; Kannan & Shanbhag, 2019; Kot- salis et al.. 2022; Yousefian et al., 2018)). Some approaches utilize increasing mini-batch sizes for variance reduction (e.g., (lusem et al., 2017; Pethick et al., 2023)). Recent works have explored alternative variance-reduced methods for (NI) and its special cases (e.g.. (Alacaoglu et al., 2023; Alacaoglu & Malitsky, 2022: Bot et al., 2019; Cai et al., 2022: Davis 2022)). However, these methods primarily adapt existing optimization estimators to approximate the operator G without significant differences. Our approach de- parts from directly approximating G. Instead, we construct an intermediate quantity S == Gr - -G as a linear combination of two consecutive evaluations of G (i.e. G and G1 We then develop stochastic variance-reduced estimators specifically for S Note that S alone is not new, but our idea of using it in stochastic methods is new. This idea allows us to design new and simple algorithms with a single loop for both (NE) and (NI) where the state-of-the-art oracle complexity is achieved (ef. Sections and 4).
where H IR:: x R TR is smooth function, and op and vo are proper, closed, and convex. Let us define = [u, E R" as the concatenation of u and e with p= P1 + P2. +2, Ga ;= VV H(u, u), V H(,, r)}, and Tr := [āt(u)), [.] (.)))) Then, the optimality condition of (SP) is written in the form of (NI). Since (VIP), and in particular, (SP) are special cases of (NI), our algorithms for (NI) in the sequel can be specified to solve these problems.
(d) Fixed-point problem. Problem (NE) is equivalent to the following fixed-point problem:
(FP) where F ;= 1- XG with being the identity operator and X O. Since (FP) is equivalent to (NE), our algorithms for (NE) from this paper can also be applied to solve (FP).
Find x² dom (F)) such that: x' Fz"
1.3. Motivation
Our work is mainly motivated by the following aspects.
(i) Recent applications. Both (NE) and (NI) cover the mini- max problem (SP) as special case. This minimax problem, especially in sonconvex-nonconcane settings, has recently gained its popularity as provides a powerful tool to model applications in generative adversarial networks (Arjovsky et al., 2017: Goodfellow et al., 2014), robust and distribu- tionally robust optimization (Ben-Tal et al., 2009; Bertsimas & Caramanis, 2011; Levy et al., 2020), adversarial train- ing (Madry et al., 2018), online optimization (Bhatia & Sridharan, 2020), and reinforcement learning (Azar et al., 2017; Zhang et al., 2021). Our work is motivated by those applications.
1.4. Basic assumptions
We tackle both (NE) and (NI) covered by the following basic assumptions (see (Bauschke & Combettes, 2017) for terminologies and concepts used in these assumptions).
Assumption 1.1. (Well-definedness) zer(VV of (NI) and
2

---

Variance-Reduced Forward-Reflectedd Backkaar Splitting Methods
zer(G) of (NE) are nonempty.
range 7 E (1/2, excludes these classical methods from recovering special cases of S.. Second, since our SVRG and SAGA estimators are designed specifically for Sh: they differ from existing estimators in the literature, including recent works (Alacaoglu et al., 2023; Alacaoglu & Malitsky, 2022; Bot et al., 2019). Third, both proposed algorithms are single-loop and straightforward to implement. Fourth, our algorithm for nonlinear inclusions (NI) significantly differs from existing methods, including deterministic ones, due to the additional term -2y - 1) - =²*). For compre- hensive survey of deterministic methods, we refer to (Tran- Dinh, 2023). Fifth, our oracle complexity estimates rely on the metrie [||| [ or E[|| + 101 for ot E Tx commonly used in nonmonotone settings. Unlike the mono- tone case, this metric cannot be directly converted to a gap function, see, e.g., (Alacaoglu et al., 2023; Alacaoglu & Malitsky, 2022). Our complexity bounds match the best known in stochastic nonconvex optimization using SAGA or SVRG without additional enhancements, e.g., utilizing a nested technique in (Zhou et al., 2018).
Assumption 1.2. [Maximal monotonicity Tf T in (NI) is maximally monotone on dom(T)).
Assumption 1.3. [Lipschitz continuity of G] G in (1) is
L-averaged Lipschitz continuous, i.e. Vx, E dom(G):
1 Ei ||G,z - -ayll - 12 yll?.
(2)
Assumption 1.4. (Weak-Minty solution] There exist solu- tion x² E zer(yV) () and k 2 0 such that (Ga + v. x*) > -<||GG + =|12 for all x E (V) and v E Tx.
While Assumption 1.1 is basic, Assumption 1.2 guarantees the single-valued and well-definiteness of the resolvent Jr of T. In fact, this assumption can be relaxed to some classes of nonmonotone operators T. but we omit this extension. The L-averaged Lipschitz continuity (2) is standard and has been used in most deterministic, randomized, and stochas- tic methods. It is slightly stronger that the I--Lipschitz continuity of the sum G. The star-co-hypomonotoniic in Assumption 1.4 is significantly different from the star- strong monotonicity used in, c.g., (Kotsalis et al., 2022). covers a class of nonmonotone operators (see Supp. Doc. A.2 for a concrete example). It is also weaker than the co-hypomotonicitty used, e.g., in (Cai et al., 2024).
Related work. Since both theory and solution methods for solving (NE) and (NI) are ubiquitous, see, e.g., (Bauschke & Combettes, 2017; Burachik & Iusem, 2008; Facchinei & Pang, 2003; Phelps, 2009; Ryu & Yin, 2022; Ryu & Boyd, 2016), especially under the monotonicity, we only highlight the most recent related works (see more in Supp. Doc. A).
1.5. Contribution and related work
Our primary goal is to develop class of veriance-reduction methods to solve both (NE) and (NI), their special cases such as ( (VIP) and (SP), and equivalent problems, like (FP).
(i) Weak-Minry solution. Assumption 1.4 is known as a weak-Minty solution of (NI) (in particular, of (NE)), which has been widely used in recent works, e.g.. (Böhm, 2022; Diakonikolas et al., 2021; Lee & Kim, 2021; Pethick et al., 2022: Tran-Dinh, 2023a) for deterministic methods and, e.g., (Lee& & Kim, 2021; Pethick et al., 2023; Tran-Dinh & Luo, 2025) for stochastic methods. This weak-Minty solution condition is weaker than the co-hypomonotoniccty (Bauschke et al., 2020), which was used carlier in proximal- point methods (Combettes & Pemnanen, 2004). Diakoniko- las et al. exploited this condition to develop an extragradient variant (called EG+) to solve (NE). Following up works in- clude (Böhm, 2022: Cai & Zheng, 2023; Luo & Tran-Dinh, 2022; Pethick et al., 2022; Tran-Dinh, 2023a). A recent survey in (Tran-Dinh, 2023) provides several deterministic methods that rely on this condition. This assumption covers a class of nonmonotone operators G or G + T.
Our contribution. Our main contribution consists of:
(a) We exploit the variable S in (FRQ) and introduce class of unbiased variance-reduced estimators 8 for St. not for G. covered by Definition 2.1.
(b) We construct two instances of 54 by leveraging the SVRG (Johnson& & Zhang, 2013) and SAGA (Defazio et al., 2014) 4) estimators, respectively that fulfill our Def- inition 2.1. These estimators are also of independent interest, and can be applied to develop other methods. (c) We develop a variance-reduced forward-reflected-yye method (VFR) to solve (NE) required O(n 2/33 evaluations of G; to obtain an e-solution.
(d) We design a novel stochastic pariance-reducee forward- reflected-backwee splitting method (VFRBS) to solve (NI), also required 2 evaluations of G
(ii) Stochastic approximation methods. Stochastic meth- ods for both (NE) and (NI) and their special cases have been extensively developed, see, e.g., (Juditsky et al., 2011; Kotsalis et al., 2022; Pethick et al., 2023). Several meth- ods exploited mirror prox and averaging techniques such as (Juditsky et al., 201 1; Kotsalis et al., 2022), while others relied on projection or extragradient schemes, e.g., (Cui & Shanbhag, 2021; lusem et al., 2017; Kannan & Shanbhag, 2019; Pethick et al., 2023; Yousefian et al., 2018). Many of these algorithms use standard Robbin-Monro stochastic
Table 1 below compares our work and some existing single- loop variance-reduction methods, but in either the co- coercive or monotone settings. Now, let us highlight the following points of our contribution. First, our intermediate quantity S can be viewed as a generalization of the forward- splitting (FRBS) operator (Malitsky & Tam, 2020) or an optimistic gradient operator (Daskalakis et al., 2018) used in the literature. However, the chosen
3

---

Yariance-Reduced
Splitting Methods
Table I: Comparison of recent existing single-loop variance-reductioon methods and our algorithms
Papers
Problem
Assumptions
Estimatoos
Residual Rates
Oracle Complexity
(Davis, 2022)
(NE) and (NI)
Co-coercive/SQM
SVRG & SAGA
linear
O((L/R) log(e
"Tran-Dinh, 2024)
(NE) and (NT)
co-coercive
a class
((1kk"
O(n + n n770
(Cai et al., 2024)
(NE) and (NI)
co-coercive
SARAH
O(1/k")
O(n + 17 log(n)e
(Alacuoglu & Malitsky, 2022)
(VIP)
monotone
SVRG
X O(n + +n"" (
(Alacaoglu et al., 2023)
(VIP)
monotone
SVRG
O
Ours
(NE) and (NT)
weak Minity
class
(1/kk O(n + +
Notes: SQM means "strong quasi monooonioot.,,,"". Residual Rate is the convergence rate on [[]]G=*II]] or [[I|Gr* + [²] for o E Tx"; and a class is a class of variance-reduccd estimators. The complexity of (Alacaoglu & Malitsky, 2022) and (Alacaoglu et al.,
2023) marked by magenta is on a gap function, a different metric than in other works in Table 1. Thus, it is unclcar how to compare them.
approximation with fixed or increasing batch sizes. Some other works generalized the analysis to a general class of algorithms such as (Beznosikov et al., 2023; Gorbunov et al., 2022a; Loizou et al., 2021) covering both standard stochas- tic approximation and variance reduction algorithms.
class of stochastic estimators for it. It also constructs two instances: SVRG and SAGA, and proves their key proper- ties. Section develops an algorithm for solving (NE) and establishes its oracle complexity. Section designs a new algorithm for solving
(NI) and proves its oracle complexity. Section 5 presents two concrete numerical examples. Proofs and additional results are deferred to Sup. Docs. A to E.
(iii) Variance-reduction methods. Variance-reduucioon tech- niques have been broadly explored in optimization, where many estimators were proposed, including SAGA (Defazio et al., 2014) 4), SVRG (Johnson & Zhang, 2013), SARAH (Nguyen et al., 2017), and Hybrid-SGD (Tran-Dinh et al., 2019; 2022), and STORM (Cutkosky& & Orabona, 2019). Re- searchers have adopted these estimators to develop methods for (NE) and (NI). For example, (Davis, 2022) proposed SAGA-type methods for (NE) under [quasi]-strong mono- tonicity. The authors in (Alacaoglu et al., 2023; Alacaoglu & Malitsky, 2022) employed SVRG estimators and developed methods for (VIP). Other works can be found in (Bot et al., 2019; Carmon et al., 2019; Chavdarova et al., 2019; Huang et al., 2022; Palaniappan & Bach, 2016; Yu al., 2022). All of these results different from ours. Some recent works exploited Halpern's 's fixed-poin iterations and develop cor- responding variance-reduced methods, see, e.g., (Cai et al., 2024; 2022). However, varying parameters or incorporat- ing bouble-loop/inexctt methods must be used to achieve improved theoretical oracle complexity. We believe that such approaches may be challenging to select parameters and to implement in practice. Finally, unlike optimization, it has been realized that using biased estimators such SARAH or Hybrid-SGDD SGD/STORM for (NT) (including (SP) and (VIP)) is challenging due to the lack of an objective function, a key metric to prove convergence, and product terms like (et, x+1 x²) in convergence analyses, where * is a bias rendered from 5tt (see Supp. Doc. A).
Quantity and Its 2. Stochastie Variance-Reduced Estimators
We first define forward-reflected quantity (FRQ) for G in (NE) and (NI) using here. Next, we propose a class of unbiased variance-reduced estimators for FRQ. Finally, we construct two instances relying on the two well-known esti- mators: SVRG from (Johnson & Zhang, 2013) and SAGA from (Defazio et al., 2014).
2.1. The forward-reflected quantity Our methods for solving (NE) and (NI) rely on the fol- lowing intermediate quantity constructed from G' via two consecutive iterates x -- and xt controlled by y E [0, 1]:
s == G.* ~GGx
(FRQ)
Here, 7 E (.. 1) plays crucial role in our methods in the sequel. Clearly, i = 2. then we can write S 1 ] + 1/ (Gr* - 11 = 1/4 [2G=* - Gr* used in both the forwardreefectee-backwwa splitting (FRBS) method (Mal- itsky & Tam, 2020) and the optimistic gradient method (Daskalakis et al., 2018). In deterministic unconstrained set- tings (i.e. solving (NE)), see (Tran-Dinh, 2023), FRBS is also equivalent to Popov's past-extragradient method (Popov, 1980), reffecteddforward--bckk rithm (Cevher & Vü, 2021; Malitsky, 2015), and optimistic gradient scheme (Daskalakis et al., 2018). the determin- istic constrained case, i.e. solving (NI), these methods are different. Since 7 (2) 1), our methods below exclude these classical schemes. However, due to a similarity pat- tern of (FRQ) and FRBS, we still term our quantity Sh by the "forward-reffected quantity", abbreviated by FRQ.
splitting algo-
Notation. We use Fx : oxx", a
xx*) to denote the o-algebra generated by z0. ... 2tth up to the iteration k. [. = | denotes the conditional expectation w.r.t.
FR. and E [ is the total expectation. We also use O (-) to characterize convergence rates and oracle complexity. For an operator G. dom(G) := {x: : Gx + 0} denotes its domain, and Je denotes its resolvent.
2.2. Unbiased variance-reduced estimators for FRQ Now, let us propose the following class of stochastic variance-reduced estimators 5 of S
Paper organization. Section 2 introduces St and defines

---

Variance-Reduced Seeird--eefeeeed-aaeeeeee Splitting Methods
Definition 2.1.
A stochastic estimator Stt is said to be stochastic unbiased variance-reduced estimator of S in (FRQ) if there exist constants p E (0, C 2 0 and C 2 0, and nonnegative sequence } such that:
(b) SAGA estimator for St. Given S as in (FRQ) and mini-batch estimator GB, as in (L-SVRG), we construct the following SAGA estimator for Sh::
S = [GB, x* - -GGB (1 - [.]
Ex - | = 0,
(SAGA)
+ 1-7 E G.¹
E [1S - I Ak:
where By is a mini-batch of size b, and G. is updated as
4 S (1- - p) + C. Un + + +--UK-I UK-1,
where U2 1 - =| G:ah-113)) 11-..
Gjak
if E Bx.
G
(5)
Here, 0, x-1 x., and [x[.] and [ are the conditional and total expectations defined earlier, respectively. The condition p > 0 is important to achieve a variance reduction as long as x is close to ===1 and == close to x²²22 Otherwise, 5 may not be a variance- reduced estimator of S Since S is evaluated at both = and _.t our bound for the estimator 54 depends on three consecutive points x²² _ and a², which is different from previous works, including (Alacaoglu et al., 2021; Beznosikov et al., 2023; Davis, 2022; Driggs et al., 2020).
G
otherwise.
To form 5k. Sy. we need to store n components G computed so
far for i [m] in a table To :- 14 C/. Gb. initialized at G. G,z²D for all i E [x]. Clearly, the SAGA estimator requires significant memory to store Ty if and p are both large. We have the following result.
Lemma 2.3. Let St. be defined by (FRQ) and S be gener- ated by the SAGA estimator (SAGA), and
A = E B[IG -G - (1 7)
Now, will construct two variance-reduced estimators satisfying Definition 2.1 by exploiting SVRG (Johnson & Zhang, 2013) and SAGA (Defazio et al., 2014).
Then, Sty satisfies Definition 2.1 with this {Aa), p==
(0, 1]. C= [2(a- b)((()))))) and := 2(n (2)
nb
nub
(a) Loopless-SVRe estimator for St.. Consider a mini- batch B& S [n] := ",,....,,,] with a fixed batch size b- : |Byl. Denote Gg.z == 1/2 E.EB, G₁z for = E dGm We define the following estimator for Sh:
We only provide two instances: (L-SVRG) and (SAGA) covered by Definition 2.1. However, we believe that similar estimators for S relied on, c.g., JacSketch (Gower et al., 2021) or SEGA (Hanzely et al., 2018), among others can fulfill our Definition 2.1.
3 = (1 (Gw GBBW*))
(L-SVRG)
+ GB,** -GBB
3. A Variance- Forward-Reflected Method for [Non]linear Equations
where the snapshot point wh is selected randomly as follows:
We first utilize the class of stochastic estimators in Defi- nition 2.1 develop a variance-reduced forward-reflected
x² with probability p
wh
k+1
(4)
w with probability 1 1-p. p.
(VFR) method for solving (NE) under Ass. 1.3 and 1.4.
The probability p E 1) will appropriately be chosen later by nonuniformly flipping a coin. This estimator is known as a loopless variant (Kovalev et al., 2020) of the SVRG estimator (Johnson & Zhang, 2013). However, it is different from existing ones used in root-finding algorithms, including (Davis, 2022) because we define it for S... not for Gr. In addition, the first term is also damped by a factor 1 - 7 to guarantee the unbiasedness of Sit to s4.
3.1. The VFR method and its convergence guarantee
(a) VER Method. Our method is described as follows. Starting from x E dom(G), at each iteration k 2 0, we construct an estimator Š satisfying Definition 2.1 with parameters p E 10, C 2 0, and Ö 2 0, and then update
x++1 = x n
(VFR)
where v > O and v > 0 are determined below, x x: -2 = x0, and 5. = (1 Gx.
The following lemma shows that 5y satisfies Definition 2.1. Lemma 2.2. Let St be given by (FRQ) and 5tt be generated by the SVRG estimator (L-SVRG) and
At least two estimators St: the Loopless SVRG estimator in (L-SVRG) and the SAGA estimator in (SAGA), can be used in our method (VFR). In terms of per-iteration complexity, each iteration i of VFR, the loopless SVRG variant requires three mini-batch evaluations Gg, w², Gg, 2², and GB, - k-1 of G, and occasionally computes one full evaluation Gw
A = nb ET E[[|Ga* - - -11-1-(11
Then, 54 satisfies Definition 2.1 with this { p :=
P.
ppp and C 22" (23+tt)) C bp
hp
5

---

Variance-Reduced Seaard--Relectee-Baaeeee Splitting Methods
of G with the probability p. It needs one more mini-batch evaluation GB compared to SVRG-type methods in optimization. Similarly, the SAGA estimator also requires two mini-batco evaluations x and GB, f=s,z=1, which is one more mini-batch GB, compared to SAGA-type methods in optimization, see, e.g., (Reddi et al., 2016b). The SAGA estimator can avoid the occasional full-batch evalua- tion Gw from L-SVRG, but as a compensation, we need to store table To == [G9, Ch ................... G*), which requires significant memory in the large-scale regime.
For O, if we choose p -n-1/3 and b = [n2 then
(VFR) requires To: == n+] 4I'L2 "PL"" evaluations of G
1 to attain K+1 E B ||| S where r := 731.
Corollary 3.2 states that the oracle complexity of (VFR) is O (n x2) matching (up to constant) the one of SVRG in nonconvex optimization, see, e.g., (Allen-Zhu & Hazan, 2016; Reddi et al., 201 It improves by a fac- tor O (m (/5) compared to deterministic counterparts. This complexity is known to be the best for SVRG so far without any additional enhancement (e.g., nested techniques (Zhou et al., 2018)) even for a special case of (NE): Gr = Vf(()) in nonconvex optimization.
(b) Convergence guarantee. Fixed E 1)) with p. C, and Ö as in Definition 2.1 we define
x1+55) M == 3(2-11) 3127-13"
1-
C+C and 8 = 27-1 (6)
8VM
Note that 77 can be computed explicitly when b and p given. For example, if n 10000 and we choose n 1/3 = 0.0464 and b = |n 464, then = 0.1456 If p = 0.1, then = 0.3038 choose appropriate p (n3) and b == O (n2/3).
Then, the following theorem states the convergence of (VFR), whose proof is given in Supp. Doc. C.
I
L Note that, in general, we can
Theorem 3.1. Let us fix E (2,1), and define M and 6 as in (6). Suppose that Assumptions 1.1, 1.3, and 1.4 hold
for (NE) with some 2 0 such that Ls S. 8. Let {x²) be generated by (VFR) using a learning rate 1) > 0 such that Then, the following bounds hold:
Alternatively, we can apply Theorem 3.1 to (SAGA).
Corollary 3.3. Suppose that Assumptions 1.1, 1.3, and 1.4 hold for (NE) with k. 2 0 as in Theorem 3.1. Let {a*} be generated by (VFR) using (SAGA), 7 == 3.4 and r ==
BR
Sn S
1 LVM
2--1
1
0111x00 11
B[IGG**1?
$
2
0.149488/2 provided that S b S n2/3, Then nL
K +1
K +1
LVM
K
1
K
E [1a
k-1 112 S
12 P2lla:0 x0111
K ++1
K+1
1
Σ k=
K E B|IG**1|1
489I. 101
bp² + 1)
(9)
1m1
2(1+ where 01 7(-7))02 and O2 == 3(27 1127-100 1)(1 - -))))))
8(1+13n7²)
Moreover, for a given > 0, if we choose b = |n2/3], then
(VFR) requires To, := n+| 3rL R K to achieve x K+I Ek B[I S Se, where I == 2816.
Theorem 3.1 only proves (L/K) convergence rate K+1 E[Gaa]]) and K+1 ΣK E |
evaluations of G
of both 1 k
x² -A11 but does not characterize the oracle complexity of (VFR). If we choose y == 3. then from (6), we have 11(C+() and 0 = I which can simplify the bounds in Theorem 3.1. In addition. it allows k > 0 such that Lx S 5 = (VP).. which means that x can be positive, but depends on vp. This condition allows us to cover a class of nonmonotone operators G, where a weak-Minty solution exists as stated in Assumption 1.4.
Similar to Corollary 3.2, the learning rate ry in Corollary 3.3 can explicitly be computed if we know n and b. For instance, if n 10000, and we choose b = |n2/3| J. then 1) = 0.1603
1'
M = 57 + 3p
24
L
If ς = 0, i.e. G reduces to a star-monotone operator, then we can choose 7 E (5. 1) and r as:
For SVRG, we have n (0, 1 If we choose
p p= -=(nn-1/) and b = (n2/3), then 1) = (++)
=
3.2. Complexity Bounds of VFR with SVRG and SAGA Let first apply Theorem 3.1 to the mini-batch SVRG estimator (L-SVRG) in Section 2. For simplicity, we choose
. For SAGA, we have 1) (0, W If we choose b b= b= then 77 O(+)
y ;= 3 and 17 ;= 1
L/M but any 7 E (g. 1) still works.
Hitherto, the constant factor I' in both corollaries is still rel- atively large, but it can be further improved by refining our technical proofs (e.g., carefully using Young's 's inequality).
Corollary 3.2. Suppose that Assumptions 1.1, 1.3, and 1.4 hold for (NE) with K. 2 0 as in Theorem 3. Let { be
generated by (VFR) using (L-SVRG), 7 == 2. and 1) :=
4. A New Variance-Reduced FRBS Method for Nonmonotone Generalized Equations
1
0.1440/bb, provided that bp² < 1. Then
LVM
In this section, we develop a new stochastic variance- reduced forvarddreffectte--aatt method to solve (NI) under Assumptions 1.2, 1.3, and 1.4.
5 5261. 2 1 - 12
E ||
(8)
splitting (FRBS)
K +1 1
bp² + 1)
0=
6

---

Variance-Reduced
Splitting Methods
4.1. The algorithm and its convergence
(d) Convergence analysis. For simplicity of our presenta-
(a) The variance-reduced FRBS method (VFRBS). Our scheme for solving (NI) is as follows. Starting from x" E (V), at each iteration & 2 O, we generate an estimator S4, satisfying Definition 2.1 and update
tion, for given y (2. 1), with p. C, and Ć in Defini- tion 2.1, we define the following two parameters:
1y M = 472 + 1-7 P
C+C
and 6=
-(2--1)
(11)
(3
xtt x² - -n n(wuttt - -227- - 1)ut), (VERBS)
Then, Theorem 4.1 below states the convergence of (VFRBS), whose proof can be found in Supp. Doc. D.
where 17 0 and 7 > 0 are determined later, whth E Tw = x=² == xº, and 2. == (1 G/"
Theorem 4.1. Let us fix y E (2:1), 1), and define M and â as in (1I). Suppose that Assumptions 1.1, 1.2, 1.3, and 1.4 hold for (NI) for some K. 2 O such that Lx < 5. Let {a"} be generated by (VERBS) using a fixed learning rate lj such
(b) Implementable version. Since u++1 Tx*++1 appears on the RHS of (VERBS), using the resolvent JyTT (I + mmf" of T, we can rewrite (VFRBS) equivalently to
that 13x 1)x <1 S
Then, we have
722
LVM
7 gh+ r*hh - -n + (24-1) (y* xx),
θ R
BII|G* 00113]
(10)
xht1 == Jynr( (on((ttt)..
K
122(K
+1)
(12)
K
22
Here, 3 E dom((V is given, and x0 = x == Jnpr("
B ||
- x 11
S
This an implementable variant of (VFRBS) using the resolvent Jyny. Clearly, if y = 2. then (10) reduces to
K
K +1
k=0
where RB. U1 and Q₂ are respectively given by
+ == ///22) - 1/21)))
R ||00 x 112 + /1 + 111
which can be viewed as a stochastic forward-reflected- backward splitting scheme. However, our y E (2. 1), mak- ing (10) different from existing methods, even in the deter- ministic case.
θ = (1->)[((>-1>>> (3y-1
(37-1)m
Ô == (1=9)(1-MDpp)
4(3--1)
The bounds in Theorem 4.1 are similar Theorem 3.1, but their proof relies on new Lyapunov function. Note that the condition on Lx still depends on p as Lx S 5 = (VD).
Compared to (Alacaoglu & Malitsky, 2022). (10) requires only one Jynp as in (Alacaoglu et al., 2023), while (Ala- caoglu & Malitsky, 2022) needs more than ones. Moreover, our estimator 5tt is also different from (Alacaoglu & Mal- itsky, 2022). Compared to (Beznosikov al., 2023) and (Alacaoglu et al., 2023), the term - (2y - 1)(y* x*) makes it different from SGDA in (Beznosikov et al., 2023)
4.2. Complexity of VFRBS with SVRG and SAGA Similar to Section 3, we can apply Theorem 4.1 for the mini-batch SVRG estimator in Section 2.
Corollary 4.2. Suppose that Assumptions 1.I, 1.2 1.3, and 1.4 hold for (NI) with K. 2 as in Theorem 4.1. Let {a*} be generated by (VFRBS) using the SVRG estima-
and the golden-ratio method in (Alacaoglu et al., 2023), and also other existing deterministic methods.
(c) Approximate solution certification. To certify an ap- proximate solution of (NI), we note that its exact solution xx E zer( ()) satisfies ||Ga* +C***+o" 1|2 = 0 for some o* E Tx* Therefore, if (2**, o*) satisfies E B[I|G2* + =*1 < 2 for some u* E Tack, then we can say that =² is an e-solution of (ND). Alternatively, we can define the following forward-
tor (L-SVRG), ? 12. 1). and n = LVM 2 ovTp with V1-7 2 2/8474775 provided that bp < 1. Then, we have K E E[I + 1
Σ k=0
(13)
K
0 o²bp² (K +1)
backward splitting (FBS) residual of (NI):
where RG = 1| x+ + 7 Gx + 101
Gnt:= = :,-==111<0 - - Jn(x - ,/Gx)),
For given > 0, if we choose p == n 1/3 and b == 1/3/3 then (VFRBS) requires To, :=
for any given 7) > O. It is well-known that xt E zer(T) iff Gjæt* = 0. Hence, if E [1[] < 2. then zt is also called an e-solution of (NI). One can easily prove that || < |!r.* + w²|| for any w E Ta Clearly, the former metric implies the latter one. Therefore, it is sufficient to only certify EE [[]]G* + [8] S e², which
4r1.3 R2/ 2/3
eval-
uations of G2 and Tr = T ft evaluations of Jayp to 6
achieve X Lk K EE [[ + < e2, where I :== at.
Alternatively, we can apply Theorem 4.1 to the mini- batch SAGA estimator (SAGA) in Section 2.
implies E E[11,,,,** [1p] S ?
7

---

Variance-Reduced -efeaardd-Reeeeee--aeeeeee Splitting Methods
Corollary 4.3. Suppose that Assumptions 1.I, 1.2, 1.3, and 1.4 hold for (NI) with K 2 0 as in Theorem 4.1. Let {et} } be generated by (VFRBS) using the SAGA estimator ab²/2
FRBS scheme (VFRBS) in (Alacaoglu al., 2023), and the variance-reduced extragradient algorithm (VEG) in (Ala- caoglu & Malitsky, 2022). We select the parameters sug- gested by our theory, while choosing appropriate parameters for OG, VFRBS, and VEG. The details of this experiment, including generating data and specific choice of parameters, are given in Supp. Doc. E.
(SAGA), 7 E 12) and 1) == LVM 2 aL with a ==
-- 2/ 71100777779
provided that 1 S 5. 1213, Then
1 K K+1
2 1,2 R3
Σ
D[[|G+*+ + + - 02bb((K K + +1)
(14)
The relative residual norm [[/*ll//|//a" against the num- ber of epochs averaged on 10 problem instances is revealed in Figure 1 for two datasets (p, n) = (100, 5000) and (p.n) = (200, 10000).
where R3 = 1|00 - 12 1 2 + 001
For a given e > 0. if we choose == (VERBS) requires TG, == n of Gj and Tr
| then evaluations
a RBna/a g²
Experiment 1:
5000 and
100
operter norm [Gx (100))
T RD 2
evaluations of J to achieve
10
K+I E E[I + 1 =,,, where I = 6
30
7
30
Similar to Subsection 3.2, when 7, n, b, and are given, we can compute concrete values of the theoretical learning rate
10
QG
J0
VFR-arx
LVFR-wng
y in both corollaries. They are larger than the corresponding lower bounds given in these corollaries.
30
VERBS
VEG
0
10
21
40
5. Numerical Experiments
Numbor spects
We provide two examples to illustrate (VFR) and (VFRBS) and compare them with other methods.
Experiment 2: 10000 and 200
Rolative operaator nears [ [1
10-1
5.1. Nonconvex-noncomme
minimax optimization
-0-2
We consider the following sonconvex-nonconece minimax optimization problem as a special case of (SP):
70
VFR-wwng
min max {C(u, v) == v(u)= + H(u, v) - b())), (15) HERPI WERP2
LVERavnn
70-7
VFR-agg
VFRBS
where H(,, v) == 1
n -ー= H4(,v)) 4 E A,u + n
u LiU 0² Biu + b, G o] such that A; E R2ixmL and B, E RP2xp2 are symmetric matrices, L, E I Xpa b; E IRP, and e E RP2; &o and vb are two proper, closed, and convex functions. The optimality of (15) becomes (NI) (see Supp. Doc. E). In our experiments, we choose A, and B, to be not positive semidefinite such that Assumption 1.4 holds. Thus, (15) is nonconvex- mnnconve--noncoeevve.
Number of spocts
Figure 1: Comparison of 6 algorithms to solve the uncon- strained (15) on 2 experiments (The average of 10 runs).
Clearly, with these experiments, three SVRG variants of our method (VFRBS) work well and significantly outperform other competitors. The LVFR- svrg variant of (VFRBS) seems work best, while VFRBS and VEG still cannot beat the deterministic algorithm OG in this example.
We generate A; QF for given orthonormal ma- trix Q; and diagonal matrix D,, where its entries D/d are generated from standard normal distribution and clipped by max{ -0.13. The matrix B, is also generated by the same way, while L,, by and e; are generated from standard normal distribution. Hence, G in (NI) is not symmetric and also not positive semidefinite.
The constrained case. We now adding two simplex con- straints u E Ap, and v E Ap, to (15), where Ap == {u e R : E u4 = 1} is the standard simplex in IRP These constraints are common in bilinear games. To handle these constraints, we set pp(u) == مممة (u) and vb()) := مسمة (v) as the indicators of Ap, and Ap2, respectively.
The unconstrained case. We implement three variants of (VFR): VFR-svrg (double-loop SVRG), LVFR-SVig (loopless SVRG), VFR-saga (using SAGA estimator) in Python to solve (15) when both yp and &i are vanished, i.c. its optimality is a special case of (NE). We also compare our methods with the deterministic optimistic gradient method (OG) in (Daskalakis et al., 2018), the variance-reduced
Again, we run 6 algorithms for solving this constrained case of (15) using the same parameters as the unconstrained case. We report the relative norm of the FBS residual 11,,11/1 against the number of epochs. The re- sults are revealed in Figure 2 for two datasets (p. n) = (100, 5000) and (p. n) = (200, 10000).
8

---

Variance-Reduced
Splitting Methods
Experiment 1: 5000 and 100
we take the nominal feature vector X; and add random noise generated from normal distribution of zero mean and variance of xه = 0.5. In our test, we choose 7 := 10 and m == 10. The relative FBS residual norm 19, (0,1f11/111 l//
Reletive uperator norm [Gx1//6X4
10.2
OG
VFR-20X
LVFR
10
saga
Viien
10-7
VEG
against the epochs is plotted in Figure 3 for both datasets.
10
The Dataset: (n. (32561. 134. 1D)
10
Relative uperator norm [G,' |
OG
VFR-svro
10"11
LVE
10-3
VFR-saga
20
60
VFABS
Number of apochs
VEG
10-2
Experiment 2: n 10000 and p 200
ative [Gx1/|6>
10-1
0G
10-3
ViR-wng
LVFR song
10-1
VERBS
10
10"1
VEG
20
40
60
80
100
Number of epochs
10-11
wea Dataset: (n. p. (45546, 311, 10)
10-11
Rellattv operatar norm 16Gx"1
svrg
10-11
10-1
LVFR-swrg
VFR.caga
80
100
VFRES
Number epochs
Figure 2: The performance of 6 algorithms to solve the constrained (15) on 2 experiments (The average of 10 runs).
10"
Clearly, with these experiments, both SVRG variants of our method (VFRBS) work well and significantly outper- form other competitors. The SVRG variant (VFR-svip) of (VFRBS) seems work best, while our VFR-saça has similar performance as VEG. Again, we also see that VERBS tends to have a similar performance OG.
60
Number of opachs
Figure 3: Comparison of algorithms to solve (16) on two real datasets: a8a and wBa.
As we can observe from Figure 3 that three variants VFR-svrg, LVFR-svrg, and VFR-saga have similar performance and are better than their competitors. Among three competitors, VFRBS still works well, and is much bet- ter than OG and VEG. The deterministic method, OG, is the worst one in terms of oracle complexity. In this test, VEG has a larger learning rate than ours and VERBS.
5.2. Logistic regression with ambiguous features We consider the following minimax optimization problem arising from a regularized logistic regression with ambigu- ous features (see Supp. Doc. E for more details):
min max [(w, === N WER"ZER"
N 1m ET (Xl(( w).u)) vc)
(16)
+ +R((w) - س2nm {
6. Conclusions
where ((T log(1 + cxp(T)) ST is the standard logistic loss, R(w)= = || is an f-norm regularizer, is regularization parameter, and سA... is the indicator of
We develop two new vuriance-reducee algorithms based on the wwwanddrffeected-backwwaa splitting method tackle both root-finding problems (NE) and (NI). These methods encompass both SVRG and SAGA estimators as special cases. By carefully selecting the parameters, our algorithms achieve the state-of-the-art oracle complexity for attain- ing an e-solution, matching the state-of-the-arr complexity bounds observed nonconvex optimization methods us- ing SVRG and SAGA. While the first scheme resembles stochastic variant of the optimistic gradient method, the second one is entirely novel and distinct from existing ap- proaches, even their deterministic counterparts. We have validated our methods through numerical examples, and the results demonstrate promising performance compared to existing techniques under careful parameter selections.
Am to handle the constraint E Am. Then, the optimality condition of (16) can be cast into (NI), where := [w, 2).
We implement three variants of (VFRBS) to solve (16):
VFR- LVFR- -vrr, and VER- VFR-saga. We also com- pare our methods with OG, VERBS, and VEG as in Subsec- tion 5.1. We cary out mannual tuning procedure select appropriate learning rates for all methods. We test these algorithms on two real datasets: a (134 features and 3561 samples) and w8a (311 features and 45546 samples) downloaded from LIBSVM (Chang & Lin, 2011). We first normalize the feature vector X; and add a column of all ones to address the bias term. To generate ambiguous features,

---

Variance-Reduced
Splitting Methods
Impact Statement
Bertsimas, D. Brown, D. and Caramanis, C. Theory and Applications of Robust Optimization. SIAM Review, 53 (3):44-5501, 2011.
This paper proposes new algorithms with rigorous conver- gence guarantees and complexity estimates for solving broad class of large-scale problems. These problems cover many fundamental challenges and applications in optimiza- tion, machine learning, and related fields special cases. We believe that our new algorithmms have the potential to make a significant impact in machine learning and related areas. Additionally, there are various potential societal con- sequences of our work, though none that we feel require specific emphasis at this time.
Beznosikov, A., Gorbunov, E., Berard, H., and Loizou, N. Stochastic gradient descent-ascent: Unified theory and new efficient methods. In International Conference on Artificial Intelligence and Statistics, pp. 172-235. PMLR, 2023.
Bhatia, K. and Sridharan, K. Online learning with dynamics: A minimax perspective. Advances in Neural Information Processing Systems, 33:15020-15030, 2020.
Acknowledgements
Böhm, A. Solving sanconvex-nonconee min-max prob- lems exhibiting weak Minty solutions. Transactions on Machine Learning Research, 2022.
This work was partly supported by the National Sci- ence Foundation (NSF): NSF-RTG grant No. NSF DMS- 2134107 and the Office of Naval Research (ONR), grant No. N00014-23-1-2588.
Bot, R. L., Mertikopoulos, P., Staudigl, M., and Vuong, methods with variance
P. T. reduction for stochastic variational inequalities. arXiv preprint / 2019.
References
Alacaoglu, A. and Malitsky, Stochastic variance reduc- tion for variational inequality methods. In Conference Learning Theory, pp. 778 816. PMLR, 2022.
Burachik, R. S. and Iusem, A. Set-Valued Mappings and En- largements of Monotone Operators. New York: Springer, 2008
Alacaoglu, A., Malitsky, Y., and Cevher, V. Forward- sefeeted-backwaard method with variance reduction. Com- put. Optim. Appl., 80(2):321-344 2021.
Cai, X., Song, C., Guzmán, C., and Diakonikolas, J. A stochastic Halpern iteration with variance reduction for stochastic monotone inclusion problems. In Proceed- ings of the 12th International Conference on Learn- ing Representations (ICLR 2022), 2022. URL https: // //openreveew.ro net wttp:/ffummmddmddd
Alacaoglu, A., Böhm, A., and Malitsky, Y. Beyond the golden ratio for variational inequality algorithms. J. Mach. Learn. Res, 24(172):1-33, 2023.
8TpCf.
Cai, X., Alacaoglu, A., and Diakonikolas, J. Variance duced halpern iteration for finite-sum monotone inclu- sions. In The 12th International Conference on Learning Representations (ICLR), pp. 1-33,2000 2024.
Allen-Zhu, Z. and Hazan, E. Variance reduction for faster non-convex optimization. In International conference on machine learning, pp. 699-707. PMLR, 2016.
Arjovsky, M., Chintala, S., and Bottou, L.. Wasserstein gen- erative adversarial networks. In International Conference on Machine Learning, pp. 214-223, 2017.
Cai, Y. and Zheng, W. Accelerated Single-Call Methods for Constrained Min-Max Optimization. In The 11th International Conference on Learning Representations, ICLR 2023. The Eleventh International Conference Learning Representations, ICLR 2023, 2023.
Azar, M. G., Osband, I., and Munos, R. Minimax regret bounds for reinforcement learning. In International Con- ference on Machine Learning. pp. 263-272. PMLR, 2017.
Carmon, Y., Jin, Y., Sidford, A., and Tian, K. Variance reduction for matrix games. Advances in Neural Informa- tion Processing Systems, 32, 2019.
Bauschke, H. H. and Combettes, P. Convex analysis and monotone operators theory in Hilbert spaces. Springer- Verlag, 2nd edition, 2017.
Cevher, V. and Vū, B. A reflected ting method for monotone inclusions involving Lips- chitzian operators. Ver-Valued and Variational Analysis, 29()):133-774, 2021.
split-
Bauschke, H. H., Moursi, W. M., and Wang, X. Generalized monotone operators and their averaged resolvents. Math. Program., pp. 1 1-20, 2020.
Chang, C.-C. and Lin, C.J. LIBSVM: A library for Sup- port Vector Machines. ACM Transactions on Intelligent Systems and Technology, 2:27:1-27:27, 2011.
Ben-Tal, A., Ghaoui, L., and Nemirovski, A. Robust optimization, Princeton University Press, 2009.
10

---

Variance-Reduced Seardd--eefeeee---aeeeee Splitting Methods
Chavdarova, T., Gidel, G., Fleuret, F., and Lacoste-Julien, S. Reducing noise in gan training with variance reduced tragradient. Advances in Neural Information Processing Systems, 32:393 403, 2019.
Gower, M., Richtárik, P., and Bach, Stochastic quasi- gradient methods: Variance reduction via Jacobian sketch- ing. Math. Program., (88(1):135-19))) 2021.
Hanzely, F., Mishchenko, K., and Richtárik, P. SEGA:
Combettes, P. L. and Pennanen, T. Proximal methods for cohypomonotone operators. SIAM J. Control Optim., 43 (2):731-742,
Variance reduction via gradient sketching. Advances in Neural Information Processing Systems, pp. 2082-2093, 2018.
Cui, S. and Shanbhag. U. On the analysis of variance- reduced and randomized projection variants of single projection schemes for monotone stochastic variational inequality problems. Ser- Valued and Variational Analysis, 29(2):453-44 49155333-9,220 2021.
Huang, K., Wang, N., and Zhang, S. An accelerated vari- ance reduced extra-point approach to finite-sum vi and optimization. arXiv preprint arXiv:22/1.03206 2022.
lusem, A. N., Jofré, A., Oliveira, R. I., and Thompson, Ex- tragradient method with variance reduction for stochastic variational inequalities. SIAM Optim., 27(2):686-724, 2017.
Cutkosky, A. and Orabona, F. Momentum-based variance reduction in non-convex SGD. Advances in Neural Information Processing Systems, pp. 15210-15219, 2019.
Daskalakis, C., Ilyas, A., Syrgkanis, V., and Zeng. H. Train- ing GANs with Optimism. In International Conference on Learning Representations (ICLR 2018), 2018.
Johnson, R. and Zhang, T. Accelerating stochastic gradient descent using predictive variance reduction. NIPS, pp. 315-323, 2013.
Davis, D. Variance reduction for root-finding problems. Math. Prugram., pp. 1-36, 2022.
Juditsky, A., Nemirovski, A., and Tauvel, C. Solving varia- tional inequalities with stochastic mirror-prox algorithm. Stochastic Systems, 1(1):17-58, 2011. 2011.
Defazio, A., Bach, F., and Lacoste-Julien, S. SAGA: A fast incremental gradient method with support for non- strongly convex composite objectives. Advances in Neural Information Processing Systems (NIPS), pp. 1646- 1654, 2014.
Kannan, A. and Shanbhag, U. V. Optimal stochastic extra- gradient schemes for pseudomonotone stochastic varia- tional inequality problems and their variants. Comput. Optim. Appl., 74(3):779 820, 2019.
Diakonikolas, J. Halpern iteration for optimoltimal and parameter-free monotone inclusion and strong solutions to variational inequalities. In Conference Learning Theory, pp. 1428-1451. PMLR, 2020.
Kotsalis, G., Lan, G., and Li, T. Simple and optimal methods for stochastic variational inequalities, i: operator extrapo- lation. SIAM J. Optim.. 32(3):2041-20735 2022.
Diakonikolas, J., Daskalakis, C., and Jordan, M. Efficient methods for structured nonconvex-nonee nonconevve min-max optimization. In International Conference on Artificial Intelligence and Statistics, pp. 2746-2754, PMLR, 2021.
Kovalev, D., Horvath, S., and Richtarik, P. Don't 't jump through hoops and remove those loops: SVRG and Katyusha are better without the outer loop. Algorithmic Learning Theory, pp. 451-467. PMLR, 2020.
Driggs, D.,, Ehrhardt, M. J., and Schönlieb, C.-B. Accelerat- ing variance-reduuce reduced stochastic gradient methods. Math. Program, (online linst):1-45, 2020.
Lee, S. and Kim, D. Fast extra gradient methods for smooth structured imconnee-nnnncome Thirty-fifth Conference on Neural Information Processing Systems (NeurlP:2021), 2021.
minimax problems.
Facchinei, F. and Pang, J.-S. dimeni-aime variational inequalities and complementarity problems, volume 1-2. Springer- Verlag, 2003.
Levy, D., Carmon, Y., Duchi, C., and Sidford, A. Large- scale methods for distributionally robust optimization. Advances in Neural Information Processing Systems, 33: 8847-8860) 2020.
Goodfellow, I., Pouget-Abadic, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. Generative adversarial nets. In Advances in neural information processing systems, pp. 2672-2680, 2014.
Loizou, N., Berard, H., Gidel, G., Mitliagkas, I., and Lacoste-Julien, S. Stochastic gradient descent-ascent and consensus optimization for smooth games: Convergence analysis under expected co-coercivity. Advances in Neu- ral Information Processing Systems, 34:19095-19108, 2021.
Gorbunov, E., Berard, H., Gidel, G., and Loizou, N. Stochas- tic extragradient: General analysis and improved rates. In International Conference on Artificial Intelligence and Statistics, pp. 7865-7901, PMLR, 2022.

---

Variance-Reduced Seeerd--Refeeet--aeeeeeee Splitting Methods
Luo, Y. and Tran-Dinh, Q. Extragradient-type methods for co-monotone root-finding problems. (UNC-STOR Technical Report), 2022.
Reddi, S. J., Sra, S., Póczos, B., and Smola, A. Fast in- cremental method for smooth nonconvex optimization. In 2016 IEEE 55th Conference Decision and Control (CDC), pp. 1971-1977. IEEE, 2016b.
Madry, A., Makelov, A., Schmidt, L., Tsipras, D., and Vladu, A. Towards deep learning models resistant to adversarial attacks. In International Conference on Learn- ing Representations, 2018.
Robbins, H. and Monro, S. A stochastic approximation method. The Annals of Mathematical Statistics, 22(3): 400-407, 1951.
Malitsky, Y. Projected reflected gradient methods for mono- tone variational inequalities. SIAM J. Optim., 25(1):502- 520, 2015.
Rockafellar, and Wets, R. J.-B. Variational Analysis. Springer- Verlag, 1997.
Ryu, E. and Yin, W. Large-scale convex optimization: Algo-
Malitsky, Y. and Tam, M. K. A forward-backward splitting method for monotone inclusions without cocoercivity. SIAM J. Optim., (002)):451-1472 2020.
rithms & analyses via monotone operators. Cambridge University Press, 2022.
Ryu, E. K. and Boyd, S. Primer on monotone operator methods. Appl. Comput. Math, 15(1):3 43, 2016.
Namkoong, H. and Duchi, Stochastic gradient meth- ods for distributionally robust optimization with f- divergences. Advances in neural information processing systems, 29, 2016.
Tran-Dinh, Q. Extragradient-Type Methods with (1/k))- Convergence Rates for Co-Hypomonotone Inclusions. J. Global Optim., pp. 1-25, 2023a.
Nesterov, Y. Dual extrapolation and applications to solv- ing variational inequalities and related problems. Math. 2007.
Tran-Dinh, Q. Sublinear Convergence Rates of Extragradient Methods: A Survey on Classical and Recent Developments. arXiv preprint arXiv:2308.17112, 2023b.
Program., 109(2-3):33
Nguyen, L. M., Liu, J., Scheinberg, K., and Takáč, M. SARAH: A novel method for machine learning problems using stochastic recursive gradient. In Proceedings of the 34th International Conference on Machine Learning, pp. 2613-2621, 2017.
Tran-Dinh, Q. Accelerated Variance- Forward- Reflected Methods for foot-Finding Problems. arXiv preprint arXiv:2406.00 02413,,, 2024.
Palaniappan, B. and Bach, F. Stochastic variance reduc- tion methods for saddle-point problems. In Advances in Neural Information Processing Systems, pp. 1416-1424, 2016.
Tran- Tran-Dinhs Dinh, Q. and Luo, Y. Randomized Clock-Coordinate Optimistic Gradient Algorithms for Root- Root-Finding Prob- lems. Math. Oper. Res., in press, 2025.
Tran-Dinh, Q., Pham, H. N., Phan, T. D., and Nguyen, M. L. Hybrid stochastic gradient descent algorithms for stochas- tie nonconvex optimization. Preprint: MeXiv:/905.009 2019.
Pethick, T., Patrinos, P., Fercoq, O., and Cevher, V. Es- caping limit cycles: Global convergence for constrained munconvex-nonconeee minimax problems. Interna- tional Conference on Learning Representations, 2022.
Fran-DDnh, Q.. Pham, N. H., Phan, D. T., and Nguyen, L. M. A hybrid stochastic optimization framework for stochastic composite nonconvex optimization. Math. Program., 191: 1005-1071, 2022.
Pethick, T., Fercoq, O., Latafat, P., Patrinos, P., and Cevher, V. Solving stochastic weak Minty variational inequalities without increasing batch size. In Proceedings of Interna- tional Conference on Learning Representations (ICLR), pp. 1-34, 2023.
Yousefian, F., Nedić, A., and Shanbhag, U. V. On stochastic mirror-prox algorithms for stochastic cartesian variational inequalities: Randomized block coordinate and optimal averaging schemes. Vel-Vel and Variational Analysis, 26:789 819, 2018.
Phelps, R. R. Convex functions, monotone operutors and differentiability, volume 1364. Springer, 2009.
Popov, I.. D. A modification of the Arrow-Hurwicz method for search of saddle points. Math. notes of the Academy of Sciences of the USSR, 28(5):845-845 1980.
Yu, Y., Lin, T., Mazumdar, E. V., and Jordan, M. Fast dis- tributionally robust learning with variance-reduced min- max optimization. In International Conference on Artifi- cial Intelligence and Statistics, pp. 1219-1250, PMLR, 2022.
Reddi, S. J., Hefny, A., Sra, S., Póczos, B., and Smola, A. Stochastic variance reduction for nonconvex optimization. In ICML., pp. 314-323, 2016a.
12

---

Variance-Reduced
Splitting Methods
Zhang, K., Yang. Z., and Başar, T. Multi-agent reinforce- ment learning: A selective overview of theories and algo- rithms. Handbook of reinforcement learning and control, pp. 321-384, 2021.
Zhou, D., Xu, P., and Gu, Q. Stochastic nested variance reduction for nonconvex optimization. In Proceedings of the 32nd International Conference on Neural Information
Processing Systems, pp. 9925-3936. Curran Associates Inc., 2018.
13

---

Variance-Reduced Reeee-eeeeeeeeccceeeee Splitting Methods
Supplementary Document: Variance-Reduced Barward-Reflected-Bateee Splitting Methods for Nonmonotone Generalized Equations
Due to space limit, some parts of our algorithmic construction and theory are not described in detail and motivated in the main text. This supplementrry document aims at providing more details of the algorithmic construction, motivation, related work, technical proofs, and additional experiments related our methods.
A. Further Discussion of Related Work and Assumptions
Let us further expand our discussion of related work in the main text. Then, we show that our assumptions, Assumptions 1.3 and 1.4, indeed cover nonmonotone problems.
A.1. Further Discussion of Related Work
As we already discussed in the introduction of the main text, both standard stochastic approximation and variance-reduction methods have been broadly studied for (NE) and (NI), including (Juditsky et al., 2011; Kotsalis al., 2022; Pethick et al., 2023). In this section, we further discuss some other related work to (NE) and (NI), their special cases, and equivalent forms.
Beyond monotonicity. Classical methods such as extragradient, prox-mirror, and projective schemes often relax the monotonicity to star-monotonicity, and other forms such as Noor, 2003; Noor & Al-Said, 1999; Tu, 2018). These assumptions are certainly weaker than the monotonicity and can cover some wider classes of problems, including some nonmonotone subclasses. Another extension of monotonicity is the weak-Minty solution condition in Assumption 1.4, which was proposed in early work, perhaps in the most recent one such as (Diakonikolas al., 2021), as an extension of the star-monotonicity and stantwaak- assumptions. Other following-up up works include (Böhm, 2022; Gorbunov et al., 2022b; Luo & Tran-Dinh 2022). A comprehensive survey for extragradient-type methods using the weak-Minty solution condition can be found in (Tran-Dinh, 2023; Tran-Dinh & Nguyen- Trung, 2025a;b). The monotonicity has also been extended to a weak monotonicity, or related, prox-regularity (Rockafellar & Wets, 1997) (in particular, weak-convexity). Other types of hypo-monotonicity or co-monotonicity concepts can be found, c.g., in (Bauschke et al., 2020). These concepts have been exploited to develop algorithms for solving (NE) and (NI) and their special cases. For stochastic methods, extensions beyond monotonicity have been also extensively explored. For instance, some further structures beyond monotonicity such as weak solution were exploited for MVIs in (Song et al., 2020), pseudo-monotonicce
and quasi-monotonicity (Konnov, 2001;
was used in (Bot et al., 2021; Kannan & Shanbhag, 2019) for stochastic VIPs, two-sided Polyak-Lojasiewicz condition was extended to VIP in (Yang et al., 2020) to tackle class on nonconvex-nomoneea minimax problems, an expected co-coercivity was used (Loizou et al., 2021), and strongly star-monotone was further exploited in (Gorbunov et al., 2022a). While these structures are occasionally used in different works, the relation between them is still largely clusive. In addition, their relation to concrete applications is still not well studied.
(b) Further discussion on stochastic methods. Under the monotonicity, several authors have exploited the stochastic approximation approach (Robbins & Monro, 1951) to develop stochastic variants for solving (NE) and (NT) and their special cases. For example, stochastic Mirror-Prox was proposed in (Juditsky et al., 201 1), which has convergence on a gap function, but requires bounded domain assumption. This approach was later extended to the extragradient method under additional assumptions in (Mishchenko et al., 2020). In (Hsieh et al., 2019), the authors discussed several methods for solving MVIs, a special case of (NI), including stochastic methods. They experimented on numerical examples and showed that the norm of the operator can asymptotically converge for unconstrained MVIs with a double learning rate. In the last few years, there were many works focusing on developing stochastic methods for solving (NE) and (NI), and their special cases using different techniques such as single-call stochastic schemes in (Hsich al., 2019), non-accelerated and accelerated variance reduction with Halpern-type iterations in (Cai et al., 2024; 2022), co-coercive structures in (Beznosikov et al., 2023), and bilinear game models in (Li et al., 2022).
(c) The challenge of using biased estimators in algorithms for solving (NI) and related problems. In optimization, especially in nonconvex optimization, stochastic methods using biased estimators such as SARAH, Nested SVRG, Hybrid- SGD, and STORM can achieve better, even "optimal" oracle complexity compared to unbiased ones such as standard SVRG and SAGA, see, e.g., (Cutkosky & Orabona, 2019; Driggs et al., 2022; Pham et al., 2020; Tran-Dinh et al., 2022). However,
14

---

Variance-Reduced Beaerdd-Refeeeee-BBaaete Splitting Methods
it becomes challenging in root-finding algorithms for solving (NT) and related problems such as minimax optimization and VIPs. One main reason for this is that the convergence analysis of these optimization methods relies the objective function as a key metric to prove convergence guarantee and to estimate oracle complexity, However, in (NI) and (NE), we do not have such an object, making it difficult to process the biased terms, including product terms such as (e", _x*11 - ==)) and (e*; _++11_ x²), where et is a bias rendered from the underlying stochastic estimator. This is currently one of the main obstacles to move from using unbiased to biased estimators in root-finding algorithms, including our methods in this paper.
(d) Comparison (Cai al., 2024). Among many existing works, perhaps, (Cai et al., 2024) is one of the most recent works that develops suriance-reduction methods for solving (NI) and achieves the state-of-the- oracle complexity, However, (Cai et al., 2024) explores different approach than ours, which relies on some recent development of the Halpern fixed-point iteration and a biased SARAH estimator. Let us clarify the differences of this work and our paper here. Algorithm 1 in (Cai et al., 2024) is a single-loop and achieves a better oracle complexity. However, it requires much stronger assumption, Assumption 3, which is a co-coercive condition. Note that this assumption excludes the well-known bilinear matrix game, or the synthetic WGAN model (37) below. Section 4 of (Cai et al., 2024) studies both the monotone and the co-hypomonotone cases of (NT). The main idea is to reformulate (NT) into a resolvent equation JNG+Z) = O and then apply deterministic variant of Algorithm I to this equation, where Jy(G+T) is co-coercive. However, exactly evaluating T)) is impractical, one needs to approximate it by an appropriate algorithm. For instance, (Cai al., 2024) suggests to use the variance-reduced FRBS method in (A (Alacaoglu et al., 2023) to approximate this resolvent, leading to a double loop algorithm. Note that this method also relies on a unbiased estimator, namely SVRG. This approach is not a direct variance-reduced method (i.e., the inner loop can be any algorithm) as ours or Algorithm I of (Cai et al., 2024). Morcover, practically implementing as well as rigorously analyzing an inexact double loop algorithm. when the inner loop is also a stochastic method, is often very challenging and technical as it is difficult to conduct a stopping criterion of the inner loop, and to select appropriate parameters. Nevertheless, our algorithms developed in this paper are simple to implement and applicable to both (NE) and (NT) whose weak-Minty solution exists. These problems are broader than the ones in (Cai et al., 2024). We also believe that our oracle complexity in this paper can be further improved by exploiting enhancement techniques such as nested trick or multiple loops as done in (Cai et al., 2024; Zhou et al., 2018).
(e) Randomized coordinate and cyclic coordinate methods for (NE) and (NI). Together with stochastic algorithms for solving (NE) and (NI) and their special cases, randomized coordinate methods have also been proposed to solve these problems, including (Combettes& Eckstein, 2018; Combettes & Pesquet, 2015; Peng et al., 2016). Recent works on randomized coordinate and cyclic coordinate methods can be found, e.g., in (Chakrabarti et al., 2024; Cui & Shanbhag, 2021; Hamedani al., 2018; Song & Diakonikolas, 2023; Tran-Dinh & Luo, 2025; Yousefian al., 2018). These methods are not directly related to our work, but they can be considered as a dual form of stochastic methods in certain settings such as convex-concave minimax problems. Studying relations between randomized coordinate methods and stochastic algorithms for (NE) and (NI) appears to be an interesting research topic.
A.2. A Nonmonotone Example (NT)
As an example of (NI), we can consider the following linear operators
Gz = Gr + g and T= Tr.
where G and T are given square matrices in R and g Rp is a given vector. Clearly, for any G, G is L-Lipschitz continuous with L := || (the operator norm of G). Our goal is to choose G and T such that Vr == Gx + Tx is nonmonotone and satisfies Assumption 1.4.
Clearly, we can choose G and T such that (G + GT) is positive semidefinite and G + GT +T TT) is not positive semidefinite. This shows that is nonmonotone.
Now, Assumption 1.4 holds if for any (x, u), (y, v) E gra(V)
t-co-hypomonotone, i.e. there exists x > such that (u - U. v) =l|u 0 (V). In the linear case, this condition is equivalent
S = (G + G +T+TT) + <GG + T) T T(G (G TT) 2
is positive semidefinite.
However, since 1 (G + G is positive semidefinite, this condition holds if 1 (T T + + ( positive semidefinite. In particular, vi satisfies Assumption 1.4.
For example, given any e > 0, we choose
-
G :=
T =
and k:= > 0.
0
0
15

---

Variance-Reduced
Splitting Methods
Then, is clear that G is symmetric and positive semidefinite, but G + ' is symmetric and not positive semidefiitt. Thus, U/ Gr Tx: is nonmonotone. Moreover, Gr Ga + g is L-LLipschizz continuous with I. == 1.
Next, we check Assumption 1.4. Clearly, G is symmetric and positive semidefinite. Moreover, since G and T are symmetric, we have
= E +
0 0
1
c = E + k( (11+6 +2)
1
-KE
M = 2
0
1
0
1
0
If we choose k ==c e, then one can easily check that M is positive semidefinite. Hence, we can conclude that V is k-Co- hypomonotone with a 6 > 0. In particular, satisfies Assumption 1.4. addition, we can choose sufficiently small such that Lu is sufficiently small, which fulfills the condition LE S 8 in Theorem 4.1.
B. The Proof of Technical Results in Section 2
This supplementary section provides the full proof of Lemma 2.2 and Lemma 2.3.
Further discussion of the FR quantity. Let us recall our quantity S defined by (FRQ) as follows:
s = Gr G1
(FRO)
As we mentioned earlier, " plays a crucial role in our methods as y E (2, 1). If 2. then we can write $1/2 Gx + 2
] (Gr* -G* = [2 - Gr* used in both the forward-reflected beacadd-ad splitting (FRBS) method (Malitsky& & Tam, 2020) and the optimistic gradient method (Daskalakis et al., 2018).
Note that if we write Gx - Gr = Ja(a")(* - 1) by the Value Theorem, where Ja(x*) = fo VG(x*-1+ T<** (ddr, then S (1 - (xx( + (aa("") 1*-1). Clearly, if is small, then 5t can be considered as an approximation of G. augmented by second-order correction term (z - 1) (called Hessian-driven damping term or second-order dissipative term) widely used in dynamical systems for convex optimization, see, e.g., (Adly & Attouch, 2021; Attouch & Cabot, 2020). These two viewpoints motivate the use our new operator S not only in our (VFR) and (VFRBS), but in other methods such as accelerated algorithms. Thus, the results in Section 2 are of independent interest, and can potentially be used develop other methods.
Other possible stochastic estimators for S One natural idea to construct unbiased estimator for S is to use an increasing mini-batch stochastic estimator as S =- 1/2 EiEs, GG, - "G,x*=1,, where By is an increasing mini-batch b 2 bb- see, e.g., (lusem et al., 2017). While this idea may work well for the general in [n]. with by := | 2 I-p2
expectation case Gx Eg
G (=r,))] [) may not be an ideal choice for the finite-sum operator (I) as by S n, which requires
Im1-pp iterations). Other stochastic approximations may also fall into our class in Definition 2.1 such as JacSketch (Gower et al., 2021), SEGA (Hanzely et al., 2018), and quantized and compressed estimators ( see, e.g., (Horváth et al., 2023)).
to stop increasing after finite iterations (i.e. O
In(n))
B.1. Proof of Lemma 2.2: Loopless-SVRG Estimator
Let us further expand Lemma 2.2 in detail as follows and then provide its full proof.
Lemma B.1. Let S = Ga k- be defined by (FRQ) and 5 be generated by (L.SVRR).. We define
A = 1 Στ i= nb
E.
E[
GGaa (1- (1-1)G,w01p].
(17)
Then, we have
E& [
S Gx GGr
E [|| 54
< A |B[lG -Gr (011-7GGi Gw 12] < Aks
(18)
A
< (1 2 A + (4-6p+3p")
Σn E G k -, -13
nbp
+ nóp
(2p
E [[llaarkk- - G 2112..
Consequently, the SVRG estimator št constructed by (L-SVRG) satisfies Definition 2.1 with De in (17), p: = - 10,
+p2 and C --p+p=) bp
C
hp
16

---

Variance-Reduced Hrrward-Reffected-Baekee Splitting Methods
Proof. It is well-known, see, e.g., (Johnson & Zhang, 2013), that 54 is an unbiased estimator of 5 conditioned on Fk. we have Eg [5]] St.
Next, let X3 == G - -,G,ah-1- u for any E [n]. Then, we have B2 [x.] Gr for any i [n]. Since BR is in Fk. using the property of expectation, we can derive
Gw
Ex [|$ 11] (.SYRR))
Ex [11 Σ BE X; |Ga + GGG - (1 -))
Es 3 Eies, [X; - E. [x:]]]]
E EieB, || - EE [x:] [
EB EieB, ||Gr - GG - (((1~)),uu G [13 - t Ex [x.]]
Ei= ||G GGi (1- (1-7) G 112 - 1 [ [x
Here, ① holds due to the i...d property of Bx. and @ holds since E& [|X; - Ex X.] || = Ex [II - (Exx (x.)) This estimate implies the second line of (18) by taking the total expectation B[ both sides and the definition of A from (17).
Now, from (4) and (17), we can show that
(17)
IE| [[[ - -G (1 7) G/
k:
12 1i=1
(4)
P)
Zi=1 18
[[[|G - GGat (1 7)., 11| 1|
mbb
E G| GG - (1 - G),
ki=
①
(1+c)(1 p)
Σ 1= G - (1 - 7) G
mh
(1+c)(1-p
Ei= E G|t
GG |G -G [11|]]
omb
nb Ei=1 E [[[
G..k-111
18 (1+c)(1- p) ΣB E G 1- GG (1 G 1112
nh
(11+)) Σ [||Gra* k -
nèe
no [p+ + 2(1+c)(1 p)
E E[|l G k ||3]
c
= (1 c1 p) A - +
2(++ (1-))> Σ E [|||aaa
21
# pp + 2(++0](1p))
Σ E |||
G - G, 112]
Here, in both inequalities ( and @. have used Young's inequality twice. If we choose c:= 2(1-p) then (1+c) (1p))
(1+c)(1-)) = (1- (1+ p(1-l) = (2-p)(1-p) = 2-3p+p2 and 2(1+c)(1-)) +p 4p Hence, we obtain
1 2'
p
A S (1 -)))
(4- +p+
Dd B[1l G. [12]
Rbp
2- 3p p²
E G.x G.rt-2112)
nbp
This is exactly the last inequality of (18).
B.2. Proof of Lemma 2.3: SAGA estimator
Similarly, we also further expand Lemma 2.3 in detail as follows and then provide its full proof.
Lemma B.2. Let St == Gx - -1 be defined by (FRQ) and 54 be generated by the SAGA estimator (SAGA), and = Skk_ -Sk We consider the following quantity:
As
ub Σ
E
[[I -G, -1-
(19)
17

---

Gariance-Reduced Reserrd-Reeleeee--Baeeae Splitting Methods
Then, we have
Ex [5]]
Sg Gx -G
E [| - | < 4 1 G -GG
(1-1)
Li=1 0:113] I Ak:
n
(20)
A
$ (1 2n
b A +
b
> 1=
B[ Gp [13]
n2
+ 2ín b) +b2²
E.
G,ah
n=b²
Consequently, the SAGA estimator 5 constructed by (SAGA) satisfies Definition 2. with A in (19), == b C= D((n-b)(2++ +b)
2n (0,
+ and C := 2(n (2n=h))2n+)7?
nb2
ab²
Proof. It is well-known, see, e.g., (Defazio et al., 2014), that 5 defined by (SAGA) is an unbiased estimator of Sh. Indeed, we have Eg [G.]] 1 Ei= G.. Es [B Gx and Ex GGs --=1] Gx*-1, Using these relations and the definition of 8% we have
E& [5*] = R& (1-?) -i== [4] (1 - BB [G + E& [GB, >1 [
(1-1) Eim1 1 G (1-2) Ein G + Gx - 7Gxh-1
Gr - Gx
Shh
Hence, 5h is an unbiased estimator of 5
Next, let X, == G -G,f--GGa---- - (1 G for any i E [n]. Then, we have E& [X.] Gx -G 1
(1-7) Di= C
for any i [n]. Therefore, we can derive
Ex [15 - S = Es 11 EiEB, X, |Ga + GG ل Σ 1=
E [1]
EE 11 EieB, X, [
J. E& |E EieBe ||X; Ex [X,] [18]
= 82 Ex EiEB, ||Grot - 7 - () [E x2]
3 Ex ||Gao 7G - G/ [Ex [x.]]²
This implies the second line of (20) by taking the total expectation E both sides. Now, from (5) and (19) and the rule (5), for any c > O, by Young's inequality, we can show that
(19)
nh E ||||G GG - (1- - 7) |-|
(1 (1- b)
E. E. E |/, -G,a- (1 - G 112]
+ a
h.h n E B[
7G,20 - (1 - x ||3]
< u+t (1 - ) Σ G - - (1 - G
+ 0t0 (1- - h) Σ E (/ -Gr (Grz -G 21)
+ -
R [[[ xk G k 112]]
< (1+c) )) 4 k- + 1 (1 2)
--)
2(1+c)
Σ E [| * Giakk
anb
2(1+hn2 (1 b Σ E E[I G.I 2|1 21
cnb
If we choose c:= inequality as
b 2m
E (0, 1), then (1 - 4)
<
Hence, we can further upper bound the last 2n
b.
2m
2n
D < (1
h A +
12( 61(2 b
Dm EE [[[ G
2n
+
2(n (2m 6)
E B[I 21122212]
18

---

Variance-Reduced
Splitting Methods
This is exactly the last inequality of (20).
C. Convergence Analysis of VFR for (NE): Technical Proofs
To analyze our (VFR) scheme, we introduce the following two functions:
Lx = -11 mG
==11 112 + ull=* 1
(21)
E Lx -((+)))00)))
P A + 1((()(1 11 - 26-2112
where y is a given positive parameter, p. C, Õ, and As are given in Definition 2.1, and x = x x.. Clearly, we have Cg 2 and E 2 0 for all k 2 0 a.s.
One key step to analyze the convergence of (VFR) is to prove a descent property of & defined by (21). The following lemma provides such key estimate to prove the convergence of (VFR).
Lemma C.1. Suppose that Assumptions 1.3 and 1.4 hold for (NE). Let {a²} be generated by (VFR) and &g be defined by for any y E [0, 1]. Then, with M == 2(1+a=7) + (++))++C) we have
Ap
E E& EE Ex4+11]
2 1 (1 M 1211²) 18 [1=x - 112
+ 17(1 - -) (m - 1 14) - 2x E |||
(22)
+ 7 7) + [][[1|G 11|
Proof. First, using k+1 = xx -1/4 from (VFR), we can expand
| + vnnG k x=11 (VFR)
| x² + mGG - -115
= || - x011 + 2y x² - x²) + 7 Gr*11
2n rz - x²) - 2x12 2 3tt) + 1
Second, it is obvious to show that
||a + yn -111| 1|a - 112 + 2ym a - x²) + 1
Third, using again x++1 = x² -n from (VFR), we can show that
- - 17 121
Combining three expressions above, and using La from (21), we can establish that
La Lxt1
||* mmGxx=1 - 1=1 1| ||, + mmG - 101 + +lllxt - -/11111 11/1 11 111
-2m ,, - xx*) 2nm ,x xx*) + +>"n"
(23)
- + (3,
Zhk I² - x²) + +2mr (Gx* , 3ik)
+ =1|la* - -- 112 21 + 1||
Next, since E& [5 = St Ga. - - shown in the first line of (3) of Definition 2.1. Moreover, since 54 is conditionally independent of x. - xx and G. w.r.t. the ofiild Fx, we have
E3 (55 y x* - xt) = (Ga* - x²) - (Gat 1 -0")
2IE [(5t, G
= 2|1 - 22
= = [1)>)>l 112 - 7 +1|Ga* - 1112 + 7// 111--112.
19

---

Variance-Reduced Beaerdd-Rfleetee-Bacckee Splitting Methods
Taking the conditional expectation Ex both sides of (23) and using the last two expressions, we can show that
Lx- ER [x++]]
2y1(Ga . x*) 2ym(Grx 2* - x²) + 1 "/=
- + 2nyEEx [ (5 x² - 2²)| + 2 En [G 34)]
- -1 + E) || + po 1,0 - -x1-11²
= =2n(1 - 7) (Gz², "Gtt,,* , - x²) + 27(1 - ->) +)) G ||
+ + 1x² 1 + e) [| + p||l** x
Since 34 is unbiased estimator of S if we denote e² : 5k S then we have E& [et] -0. Hence, can show that
Ex || = E& ||s + ex [|-] = || + 2 [(e*, s + E& |||e* = E& ||| || + ||S 11. Using this relation and
S Gx - -GG we can show that
Ex [||
||s?|12 + E& [||]~113] 1/ - 110 E& [|]]
||Gx²1| - ( + +221| + E& [1 [1]
= (1 - >) 112 7(1 - 7) |G + + - 1r* + Ex [[]cb[]3]-
Substituting this expression into the last estimate, we can show that
Ls Ex [f
2n(1 - ())( at - x²) + 1/2 7) - 11 1))1
+ +1 - 7) H)ll (1 u Gx* 1112
- -1 + 1) ||2 + ul|=** -
Taking the total expectation E[-) both sides of this expression, we get
B[Ax] - B[cx+1] = 2(1 - >)) (Ga+ .² - xx0= + n² (1 + x) [|]
+ + - (2 - 1 11 || + μB || - -111
- + =1
-z 111 21 + [||c 1|
By Young's inequality in ( and (2) of Assumption 1.3, we have
1 1 E G112 1 Ei |G, 11
||Gx* Gz 112
(24)
(2)
2 L 12 x
Utilizing this inequality, (Ga* x - x²) 2 - x|| 11 from Assumption 1.4 with T = 0. and E [||]*1|]] < Aa from (3), we can bound the last expression as
E [[CA] E E[Ck+1] 2 [u-
p- [)) 1|1]
+ (1 + (1) - -2)
(25)
+ (1 (1
- 1-
2x] [[I 21 + n))
By the third line of (3) in Definition 2.1 and again (2), we have
A S (1 ()- CT. [12] +
x
Rearranging this inequality, we get
A S p (AR - AR) +
ŌI. E
[|
k-1 *
112 - E
[1|,,* 2 -1111
(C+0)12 E ||= - 112]
p
20

---

Variance-Reduced
Bhwnard-Reffccte-BBaccck Splitting Methods
Substituting this inequality into (25), show that
B[CN] - E [[ 2 - (1=))
- + 11(1 -1 2
1 x [13]
||
+ (1 1) - -
13770(1+1) EE| 11 ||- x²
||(
||
11
(((+)))-p)) (Ax-1 AR).
Rearranging this inequality and using & from (21), we obtain (22).
Now, we ready prove our first main result, Theorem 3.1 in the main text.
Proof of Theorem 3.1. Let denote by M 2(1+ + (I+u)(++0), Then, to keep the right-hand side of (22) 1 positive, we need to choose the parameters such that 1²2 S AL'R² (24-1-w)² S 12 S M
PP
and 17 2
These two conditions lead to 27-1-H
2x
M
Now, for given y E (g: 1). let choose μ = 3(2y-1) 0. Then, the last condition holds if La S 0 := 2y-1
4
8vM as
stated in Theorem 3.1. In this case, we have M 7(1+57)
3(2--1)
+ 3(2--1)
1+67
C+C as stated in (6). Hence, we can choose
- Sn
1 as claimed in Theorem 3.1. T/M
1+62 Next, utilizing μ + = 2
2 1 and μ = 3(27-1) (22) reduces to 4
E E|
E E( 312--1) (1 M LM-2)))) E[ x
7(1 7), ||
Averaging this inequality from h: = 0 to k == K, we obtain
k
E B[1|
K+I
Ex=
11
E [||I
1 <
K+1
k=
1I(K+)
Finally, since = 20 and S chosen as S == (1 Gzz we have A Do Using this fact, Gx* O.
the Lipschitz continuity of G, p [0, 1), and < 1,
show that
Eol = || + GG(()) -
- 1p + n (1(+))(1-p) Do
s 21 x² [3] + 2
12 + (1+67)n² $ 2 1B - 11] + (1+6))n
Do
4p
10 Do
< 2 (1+ 1.31)))) -
Substituting this upper bound into the above estimates, we get the second bound of (7). For the first bound, we replace - 1 by k, and K by + 1, using ||G=0|2 2| and then multiplying both sides of the result by line of (7).
K+2 to obtain the first K+I
Next, we restate Corollary 3.2 for the case 7 E () instead of y = 3 as in the main text. Then, we derive the proof of Corollary 3.2 from this result by fixing y =
Corollary C.2. Suppose that Assumptions 1.1, 1.3, and 1.4 hold for (NE) with k 2 0 as in Theorem 3.1. Let {a*} be generated by (VFR) using the SVRG estimator (L-SVRG), 7 (+) 1), and
7) LVM
where A == 4(1+7 (0 1p) and M 2(1+5y) 3(2--1) + 3(2-11)
6p
1+6y
A.
(26)
21

---

Variance-Reduced
Splitting Methods
Then, we have 17 2 ovbp L foro: V84497+ 13, 48y5
3(21-1)
and the following bound holds:
K+I ER=
B[IGT++ < 201+2212)3
7[1->)p((-1) where Ro == 1|x00 x*11.
(27)
For given tolerance c > O, if we choose p := n 1/3 and b == [n2/3], then (VFR) requires To, := +
1/3
ITL
Rg
6²
evaluations of G2 to achieve 1 K+1
B [[| [P] S è, where r 2(5- 47144)7(8--8111 3)(8+)444-7--1 1372 +487²)
27 ([]-1))())(1455))
Proof. For the SVRG estimator (L.-SVRG), by Lemma 2.2, we have
B E (0, C ==
1-6p+3p2 and C == up
8(1+72) and M in Theorem 3.1 as Therefore, we can compute == C+C 4(177) -2) 223+2777 pp I bp M := 312 +
hp
bp
+11+57) 1+67
7(1+57) 3(2--1) A v 3(27-1) + 3(27-1)6p² as stated in (26). The estimate (27) is exactly the first line of (7).
8(1+6:)(1+2))
Now, suppose that hp² 1. Then, by (26), have M < 8 19 Therefore, if we choose 1j 3(2y-1)6p) 1)6p"
then "
LVM M
13(2--1/Vbp avbp where v /8440,+ 43, L
satisfies the conditions of Theorem 3.1, provided that Lp < 8. Moreover, we have 17 2
/3(2-11) V8+499+ 13,2
21+L211) From (27), to guarantee K+1 ΣK E [| S e,, we need to impose >(->))"(K+1) S 2²,, where Ro == || x=
However, since + 12 = 1+ 1
7733 and 17 2 avop M S (1+5g)
the last condition holds if we choose K
r
1 R
L +77 19733 (2- 1)(1 (1-7)
bp*
2(5 +77-3 where I == ,- (=72(11))(++5)
Finally, note that, at each iteration k, (VFR) requires 3 mini-batches of size b. and occasionally compute the full batch Gw leading to the cost of np + 3b. The total complexity is
Te = K(mp + 3b) TL (m(+++3) bp'r²
TL n
+ 3. bp p²)
If we choose b == | and := n 1/3 then bp² = 1 and To = one needs to compute Gw,, which requires n evaluations of G,. Hence, the total complexity of the algorithm is To. 2/3 22 as stated.
2 For the SVRG estimator (L-SVRG),
n +
Proof of Corollary 3.2. Since we fix 7 = 3. we can easily compute a = V8+49y 137
3(2--1)
0.144025 2 0.1440
and r
2(3- ()
730.736842 S 731. Therefore, obtain v) > O. 14 and To. == n +
37 1122-10019)4466
Tn 2 where I' : 731. Moreover, (27) reduces to K+I > [[[]G*"]]] < 30.1100Up bp-(K+I)
321 526-1 Ra
Finally, we also restate Corollary 3.3 for the case 7 E (2. 1) and then derive the oracle complexity of Corollary 3.3 from this result by fixing := 3 4.
Corollary C.3. Suppose that Assumptions 1.1, 1.3, and 1.4 hold for (NE) with k 2 0 as in Theorem 3.1. Let {a*} be generated by (VFR) using the SAGA estimator (SAGA), 7 E (b, 1). and
:= LVM
where == 2
4(1+ 2(1 bi((nb))
and M 7(++37)
3(2--1) +
3(2--1J
1467
A.
(28)
0/3/2 nt.
1)
Then, we have 17 2 for 0
and the following bound holds:
v V10+61 100-11-13724480 44.
Lk= B[||Gr+|P]]
S 201+L31118
where
Ro == 111:00 - x111.
(29)
KtI
7[1->90((K11)
Moreover. for a given tolerance e > O, if we choose |/ then (VFR) requires TG, := n
3r1. n2,
evaluations
6
of G2 achieve 1 E S e where :=
= 2(7h (2-1-)))2-1
22

---

Variance-Reduced BerwarddRefleceel-Backky Splitting Methods
2n E (0, 1]. C = [(m-b)e2n and Ö = niz
Proof. Since we use the SAGA estimator (SAGA), we have p ==
2(n bX20+))", In this case, since b 2 I, we can easily show that A == C-C c b()
S + 411+>")n" 6
nh² Hence, M in Theorem 3.1 reduces to
M >(1+67) 3(2--1) +
1+6y 3(2--1)
A v
2+1
8(1+ )
3(2--1)
3(27
Suppose that 1 b S n2/3, Then, one can prove that M <
2+13745-2
((++7
n²
(10+61
448 11447
3(2
3(24-1)
3(2--1) 136
Hence, if we choose 17 obtain (27) from the first line of (7) as before.
0 where a 110+611+11+94488
(3(2-11
then we get 17 2
aba/2
stated. Moreover, we
v
TVM
mL
Now, for 77 LVM
1
2 aba/2 from (27), to guarantee
2(1-12.p)m2 K+1 2k=u E [|]]] S Se, we need to impose
7(I-7)*((+1)
e, where Rg : ||x=0 x² ||. Since 1 1+ 12 = 1+ M - (1+5m)
7y45y²-3 and 1j 2 nt. the last condition holds if we choose (1-- 138
ab3/2
n where := o. (1-+55)))
2 --3)
207y45y2 (1 -1-0)))3)1)))
K == r 6
Finally, at each iteration k, (VFR) requires 3 complexity is
of size b, leading to the cost of 3b per iteration. Hence, the total
3r1.2 8
Te == 3bK =
If we choose = [n2/3], then Te 3T S0TL2R n For the SAGA estimator (SAGA), one needs to compute Gw which
3T n/
requires n evaluations of GG. Hence, the total complexity of the algorithm To. = n +
3
Proof of Corollary 3.3. Since we choose := 4. we have o :=
0.14948 2 0.1494 and :=
V10+6l 13,7
+487") 2815.8 < 2816. Applying the results of Corollary C.3, obtain our conclusions in
Corollary 3.3. Moreover, (29) reduces to K+1
Σx
k=
E
[[
[13] S
32( 32 R2
(p+(K+1)
489 489.1, 1.22
as stated.
D. Convergence Analysis of VFRBS for (NI): Technical Proofs
One key step to analyze the convergence of (VFRBS) to construct an appropriate potential function. For this purpose, we introduce the following function:
Cg == |la + (m + wt)
112 + ul|a* zh-1 + mn >+)
(30)
where / > O is a given parameter and w E Tar is given. This function is then combined with & from (21) to establish the convergence of (VFRBS).
Let us first state and prove Lemma D.1, which provides a key estimate for our convergence analysis of (VFRBS) in Theorem 4.1.
Lemma D.1. Suppose that Assumption 1.3 holds for (NT). Let fx+² be generated by (VFRBS), Ly be defined by (30), and Eg be defined by (21). Then, we have
LE - Es [C+11] 2 2(1 - >)n + ut, - -x*)===>>>> (1 (1-1)1)22 - 1)n??)|Gr' 1/112
+ > - - 1x33-
Gx + 2 (1+ 4) [
(31)
+ [ - 2(1 11)((1 7) 1,2, -
If. additionally, Assumption 1.4 holds for (NI), then we have
B[Ea] E [[Ex+]] 2 1 1- (1 + 2)) 22 2 - -
+ - 1)]r"" 1 + , [2]
(32)
+ (11-)) (1+ 1)n 2x LE [I 1811
23

---

Variance-Reduced
Splitting Methods
Proof. Let us introduce two notations w Gæ + ve and ú* == Gr + ot, where vy Ts We also recall S S/=== = Gr -Ga and ex = 5t st from (FRQ). Then, it is obvious that S S e Gr GG ex.
Now, using 3tx Gx - -GG 1 + +e,, it follows from (VFRBS) that
x
x² n mmu - -27 1)me
= r² - mmGGx + + - (1 ))((G** + + -m((G= + +20)
ne
(33)
= x² - -mputtt - (1 - 7)qo + tomit - nne
Then, using (33) and wh+1 Gr bkt1, can show that
Ti = || 1 + >m(( + x2+1)= x-112 = || - x² + ym
(33) 11 - -mmohh+ +1 (1 ?)p x + Tppot net - x² + 11
= || - -x*112 - 2((1-- - (w0 x²h - x²) + 2ym I - x²) + 17211=*11
+ (1 - 11=)) 1|2 27 (wt, rit)
2n7 x² x²) + 2(1 ?)? (e", w²) 2ym2 (et itt).
Alternatvell,, using th = Gx + v²t we also have
TØ: == ||* + (n(C + u*)
) - 1=1 -==1==113
= || x² + m1
I|** - x=1 + 2ym/tw x² - x²) + +22,7² 11=
Subtracting T from T3) we can show that
Ts = 1| + "n + )²) 2012 ||z + m + xhtl)) - 11?
2(1 >))(/o x x²) (1 1 7) ? (| 122 2y -)) ,
+ 2n1(" x x²) - 2(1 7)n² (et, w²) ) + +ym² ( - 17² 12222
(34)
= -2(1 - >)0 x - x²) + (1 - (2 - 11
+ 7(1 - 7)17 100001 7(1 - 7) - x -
+ 2 k xt x²) - 2(11 7) (eth w+) 2ym (0 wb - 1121221122
Next, using again x++1 = Gx + ++ and (33), we have
Ti4] == || | +1 - a² + Ty + o²11)11== 11+21 x²t + +
(33
r [| - ?) - yúi + 112
(1 -))22121 1|1 1 2yy(1 7)) (n+* wht + 1
+ , + 2( - -7) w/y - (e 2ith)
= -(1 - 727 1),>?|/u*11 + 7(27 - 11)1² 10 l10* 12 + (1 - ?) -
+ + x + 2(1 - (e wh) - 2 20*).
Moreover, by the Cauchy-Schward inequality in @ and Young's inequality in @. we can prove that
11 k-1 vm0 1|
b - 1122 2mm , bb11) + 7,n || |
| - - 1|2- - 2 ( 16 ||| - - 11 + 2 ||000
Ive - 111 - > 2100
-
24

---

Variance-Reduced -eaed--Reeeeeeeeeceeecee Splitting Methods
Combining the last two expressions, we can show that
Tia ||0 r=-1 + mmGG + e )12 |1 - xt + mmGr* + 11
1| x=*1 + 1| - || 1 - x² + 111
2 || 111 (1 - (2) - 1)1²² | 112 1----77 -3 - 110
7(1 1) 20 - s7² 1| - 2(1 - -)) (e", , w²) ) + 2yn (e", 20))
Multiplying T(s) by μ > 0. and adding the result to (34), and using Ls from (30), we have
x²' 1|? Lk - La+1 || + mG + n²t) x 1| - || + mG + xkt1)
+ u||=** - xk-1 + + (Ga + )²² / x² + (G + +
2 2(1 ?)) ,,,,,,<,,,000 x x²) ) + $1 1| + (1 +1)( 7) - -1 1),1(110*11
+ 7((1 - 7) 1137 1)] 1|00
(1 -0
+ (n(e" x² - --))- 2 (1) (c w²
+ 2(1 + (ct, wb*) - (1+ 1)n
Taking the conditional expectation Eg [.] both sides of this expression, and noting that
Ex [ee wt - ==*}] = (Ex [e*]] s - x²) = 0,
Ex (eR; w
(Ex [e*] ,, = 0,
Ex (ett [
(Ex [e*]. 2002 = O,
we obtain
Lx Eg [xt1] 2 2(1 - ")n(w*, 1thh - x²) + 2. 21100 111 + (1 ())
+ 7 - -) - |1(3 - 1)]| 1>)b>"/10*11 (1 + 1) - >)n?(1 10/112
(1 T [>)+
Finally, by the L-Lipschitz continuity of G from (2) of Assumption 1.3, we have I - x |[Gr 1x 1.21 - - 1|2 as shown in (24). Using this inequality into the last estimate, we can show that
L Ex [Ce+1] 2 2(1 - (w x - x²) + (1 #1 - 11
+ 7 9 - p(3 1)]n 110 (1 +1+ B (a
+ 1/2 [a- 2(1 1)7(1 7) - - 1
which proves (31) by recalling w : Gx + + and uh = Gz + -
Taking the full expectation of (31) and using (Gx* + uk, 2k - -x*) 2 -xi -ril + 211 from Assumption 1.4 and E2 [[]] < A from (3), we can bound it as
[[[c]] 16 [C++1 2
[ - 2(1 + +1) 2) LE [12t* - 11-1 - -1 A
+ + 7 - 103--1333111111
- 1))) [ Gr + +
(35)
+ (1 (1->)n((1+)) (1+ 1)m 2x E [| + u11..
By the third line of (3) Definition 2.1 and utilizing again (2), we have
4 s (1- A) CL'E(II*A z + ÓL [1 2
25

---

Variance-Reduced Seeerdd-Reeeeee--Beeeeeee Splitting Methods
Rearranging this inequality, we get
P (AR-) - A&) ŌL
E ||
111
-B11x* x6 1111
(C+0rL E [[1e** - 2
Substituting this inequality into (35), we can show that
B|ck| Cx+1]
2
1 2 - 2(1 2 (h((1 2121P(144(C4C) 11 x² 112]
+ 7|1 - 113 - 1)) [|| 11?
(1 - 7)1 (1 + (2 - 1)1 2x B |||| [|
L' (((x()
P | || -
E8 -
10++1(111)) ml (Ai-1 -+))
Rearranging this inequality and using & from (21), we obtain (32).
Now, we are ready to prove our second main result, Theorem 4. 1 in the main text.
Proof of Theorem 4.1. Since we fix y E (2:1) and μ == 3,11 we have μ 0 and + p
27 3y-1-
Let us denote by
M = 4-2 +
4y C+C as in Theorem 4.1. Then, (32) reduces to
I-
E E[Ea] E E[Ex+1]
())))
E [|1 x² x-
213-11)
(36)
+ 2(1 >)
1(2--1)n
E ElI + v 1
3--1
Let us choose 77 >> 0 such that 7(24-1)m K > 0 and 1 1-M-L" 13 2 0. These two conditions lead to (8y-1)M 3-
< 17 S LVM
7(4-1)
as stated in Theorem 4.1. However, this condition holds if 1,2x² <
122-111 This condition is equivalent to Le S as our
M(3--13
condition in Theorem 4.1, where 5 = 112--11
(3y-1)yM
Averaging (36) from k 0 to K and noting that [[Ex]] 0 for all k 2 0, we get
K
E
[[ + u* x+112 122
(3y- 1(E
K+1
0=
2(1-7)h(27-
1)4 (3 1)
M 2 K+1
Ek k= B [1|a** 2x 111 v 2(3--1)E(Eo)
(K+))
Finally, since x- -=---- - 2- == 20, we have A Ag. However, since 50 (1-7)Gxx SG. we get Ao = -
Using these relations, p E [0, 1] and y < 1, we can show that
((E) [11 + yG( + 1-]
2 (11+))(1-p)) Ao
S 2E [| x² + 2, 2 [| + +001 +
27m""
Ao
(37-1)p
= 2 ||0000 x" [| + 2 [| + =113].
Substituting this upper bound into the above two estimates, we get two lines of (12).
Finally, we prove Corollaries 4.2 and 4.3 in the main text. Unlike Corollaries 3.2 and 3.3 where we fix 7 == 3. here we state these corollaries for any value of y E (1.1).
Proof of Corollary 4.2. For the SVRG estimator (L-SVRG), we have p = B (0, 1]. C = 1-6++3p" ò ==
up
(2-3ppppp) and Ao 0 due to (17) and 2.0
--1 w In this case, we have A = C+C
p
1(1+ (2 2233222""") 7
8(1+7") and thus M in Theorem 3.1 reduces M 4,2 + - < 472 bp bp
32(1+7)
26

---

Variance-Reduced hewaarddRRffcctee-BBcckyy Splitting Methods
Suppose that 6p² & 1. Since A - 8(1+7")
bpp and M- = 4, +
47
S 47 + 32%(1+7")
(1-7)bp² S
1x(++477")) (1-3)bp"
we choose
1-7
1) = LVM then we have 7) V--Vvb
2L(/((+77773))
ovbe L with
2
then it satisfies (3-1)K
7(2-11)
21 LVM
in
Theorem 4.1, provided that Ls I 8. Note that using 1) 2 ovbp in (12) of Theorem 4.1 we obtain the bound (13).
1.
Now, from the first line of (12), to guarantee समा E E || + b S ش we need to impose (((K+I) S e where
efi
gvip the last condition holds if we choose K := r R] = | 2=1 + 12 Gx + v01... Since 2 L
R
where
bp
: -
Finally, at each iteration k, (VFRBS) requires 3 mini batche of size b, and occasionally compute the full Gw Ieading to the cost of np + 3b per iteration. Thus the total complexity is
Te == K(np + 36) r1. (mp++)3)) bppr
r1 4T 2²
7.)
hp
If we choose b : [n2/38 and p n- 1/3, 1133, then bp² = and To
For the SVRG estimator (L-SVRG), one needs to compute Gw which requires evaluations of G4, Hence, the total evaluations of G, is To, + 4rm Moreover, at each iteration, we need one evaluation of Jnny. Therefore, the total evaluations of Jnnp is Tr = K
2²
R r L2AE
r
Proof of Corollary 4.3. Since we use the SAGA estimator (SAGA), we have p == b E (0, 1), ;= == (1
and C := 2(nn +bm In this case, since b 2 1, we get A == C+C = 3/4 + 4(2+>2 ()-
2 + 801 Hence,
nb²
83
M in Theorem 3.1 reduces to
M 47 + A 4n 42(2+ 1
I-7
(1- 11-7,6
Suppose that 1 6 6 S n 2/3 Then, we can show that M S 1-(2+)-
+
32t
n²
40104++777)
where d
n²²
1-y
1-y
(1-7)03
gb2/2
ab3/2
2 +70+ Hence, if we choose 77 := LVM then get 17 2
Note that using nj 2
in (12) of
aL
mL
Theorem 4.1 we obtain the bound (14).
For f = LVM 2
1
0b2/2 from the first line of (12), to guarantee
K+1 E
1
E [l| + 1²1 S e, we need to impose
nL
7(K+I] S 2, where RB = ||~00 x²1 + G + 111 Since 17 2 nL the last condition holds if we choose
ab3/2
K = [
L' L2 where 1 == θ Fee
Finally, at each iteration k, (VFRBS) requires 3 mini-batches of size b, leading to the cost of 3b per iteration. Thus the total complexity is
To == 36K- = 3 RE
6222
If we choose b = [n2/3], then Te 3 2/, For the SAGA estimator (SAGA), one needs to compute Gwo, which requires evaluations of G,. We conclude that (VFRBS) requires To, = n 3T 2/3 evaluations of G;, Moreover, since each iteration, it requires one evaluation of Juyr, we need Tr := K = =
L Rg
evaluations of Ju
Remark D.2. For the SVRG estimator, if we choose y = 3 1. then we have o == 0.0702. Hence, we have 7) > 0.0702.6p
0.1027.6p If we choose b = |// and n -1/3,, then the latter lower However, if we choose := 0.55, then 17 2 L bound becomes 17 > 0.1027
For the SAGA estimator, if we choose = 3 4. then we have a := 0.0753. Hence, we get >
0.075383/2 mL
However, if we set
7 : 0.55, then 7) 2 mL
0.127/83/2 If we choose = | then the latter lower bound becomes v) 2 L
0.1271
Note that these lower bounds of r can be further improved by refining the related parameters in Lemma D.I, and carefully choosing P in the proof of Theorem 4.1.
27

---

Variance-Reduced -eee--RefeeeeeBacceeeee Splitting Methods
E. Details of Experiments and Additional Experiments
Due to space limit, we do not provide the details of experiments in Section 5. In this Supp. Doc., we provide the details of our implementation and experiments. We also add more examples to illustrate our algorithms and compare them with existing methods. All algorithms are implemented in Python, and all the experiments are run on MacBookPro. 2.8GHz Quad-Core Intel Core 17, 16Gb Memory.
E.1. Synthetic WGAN Example
We modify the synthetic example in (Daskalakis al., 2018) built up on WGAN from (Arjovsky al., 2017) as first example. Suppose that the generator is simple additive model Gø(z) - 0 + with the noise input generated from normal distribution. (0,, I), and the discriminator is also a linear function Dp = (K) 3, w) for given matrix K, where E R and / E R and K E TRPx a given matrix. The goal of the generator to find a true distribution = 0", leading to the following loss:
[(), == E=- n) [K3, w))
EE- [KK 0 2)].
Suppose that we have n samples for both w and leading to the following bilinear minimax problem:
inf sup (= 0ERPI BERP2
Σ [K3, -; - -0] (
(37)
Here, we add two convex functions f(0) and (B3) to possibly handle constraints or regularizers associated with 0 and B, respectively.
If we define := 10, 3] 1 Gz [0) (3)) (g))) then the optimality condition of this minimax problem becomes 0 Gr + Tz, which is a special case of (NI) with Gr being linear. The model (37) is different from the one in (Daskalakis et al., 2018) at two points:
[K3, / Σi i- K (w; <T(w;
- 0)} and T :=
It involves a linear operator K, making it more general than (Daskalakis et al., 2018).
It has two additional terms f and g. making it broader also cover constraints or non smooth regularizers.
In our experiments below, consider cases:
Case 1 (Unconstrained setting). We assume that 0 E IRP: and / E IRP2,
Case 2 (Constrained setting). Assume that 0 and 33 stays in an -ball of radius r > 0. leading to f(0) := 6(-r,)p1 (0) and ((3)= := -/-rrj (3), the indicator of the -balls.
E. 1.1. THE UNCONSTRAINED CASE
Algorithms. We implement three variants of (VFR) to solve (37).
The first variant is using a double-loop SVRG strategy (called VFR- svrr), where the full operator Gw" at a snapshot point w is computed at the beginning of each epoch s. Then, we perform |n/b] iterations k to update x² using (VFR), where à is the mini-batch size. Finally, we set the next snapshot point w=+1 := x++1 after finishing the inner loop. The second variant is called a loopless one, LVFR-svrg, where we implement exactly the same scheme (VFR) as in this paper and using the Loopless-SVRG estimator.
The third variant is VaRasaga, where
the SAGA estimator in (VFR).
We also compare our methods with the deterministic optimistic gradient (OG) in (Daskalakis et al., 2018), the variance- reduced FRBS (VFRBS) in (Alacaoglu et al., 2023), and the variance-reduced extragradient (VEG) in (Alacaoglu & Malitsky, 2022).
(b) Input data. For (NE), we generate vector 0* from the standard normal distribution as our true mean in R., Then, we generate i.i.d. samples w; and 2; from normal distribution N/(0", II) and (000, I), respectively for i 1,2, , in RP and TR respectively. We perform two expertiments: Experiment I with n 5000 and p1 = P2 100, and Experiment 2 with 10000 and p1 p2 200. For each experiment, we run 10 times up to 100 epochs, corresponding to 10 problem instances, using the same setting, but different input data (w,, 2;), and then compute the mean of the relative operator norm This mean then plotted.
|G ||/[[/]/]/.
(c) Parameters. For the optimistic gradient algorithm (OG), we choose its learning rate 1) == +. where L is the Lipschitz constant of G. though its theoretical learning rate much smaller. For our methods (VFR), 5000, and we choose b [0.5m2/3 2/3 146, and the probability p n4/
0.1170, then 17 LVM
1905 However, due to the under
28

---

Variance-Reduced
Splitting Methods
estimation of M, we instead use larger learning rate 1j == 2I. for all three variants, and choose mini-batch of size b = J. and a probability p = + for the loopless SVRG variant.
7277
For the rrrward-reflected-aaaeo its learning rate 1) : 1.96(--II-P as suggested by (Alacaoglu et al., 2023). However, we still choose the probability 21.
splitting method with variance reduction (VPRBS) in (Alacaoglu et al., 2023), we choose
p n// and the mini mini-batch size b = [ as methods. These values are much larger the ones suggested in (Alacaoglu et al., 2023), typically p= (1/n).
For the variance reduction extragradient method (VEG) in (Alacaoglu & Malitsky, 2022), choose its learning rate 1 17 and = [0.5n 2/3 in this method, which is the same as ours, though their theoretical results suggest smaller values of p (e.g., p .,). Note that if n 5000, then the batch size b= 150 and the probability p =- 0.062, but if n 10000, then = 239 and = 0.0479.
17 := 0.95y/1a for a == 1 - p from the paper. However, again, we also choose ==
(d) Experiments for K = L. We perform two experiments: Experiment 1 with (n, p) = (5000, 200) and Experiment with (n, p) (10000, 400) as discussed above. We run each experiment with 10 problem instances and compute the mean of the relative residual norm |/|| The results of this test are plotted in Figure
Experiment 1: n 5000 and = 200
Experiment 2: n 10000 and p = 400
Relatvv operator |G/1 [//A//
Relative operator norm |G>A1/((a))
10--
10-2
10
10"
10-9
10
10
OG
DG
10-19
VFRisug LVFR surg
10
VER-sarg ovig VFR-sago
10-13
VFR-saga
10
VFR85
VERBS
10-14
VEG
VEG
Number of epochs
Number of epochs
Figure 4: Performance of 6 algorithms to solve (37) on 2 experiments when K I.
For these particular experiments, our methods highly outperform OG, VFRBS, and VEG. It shows that VFR-svrg is the best overall, while VFRBS and VEG do not perform well in this test and they are much slower than ours and also OG. This is perhaps due to a small learning rate of VERBS although we choose the same mini-batch size b and the same probability p as ours.
and VFR- svrg have a similar performance in both experiments. Both the competitors:
(e) Experiments for K A I. Now, we test these 6 algorithms for the case K + in our extended model (37), where K is generated randomly from the standard normal distribution. Then, we normalize K as K/||K]| to get unit Lipschitz constant L = 1.
Again, we use the same configuration as in Figure 4 and also run our experiments on 10 problems and report the mean results. We perform two experiments: Experiment 1 with n = 5000 and p1 = P2 = p = 100, and Experiment 2 with n 10000 and pl = p2 = = 200. The results are reported in Figure 5.
We still observe that our algorithms work well and outperform their competitors. However, after 100 epochs, these methods can only reach 10-2 accuracy level for an approximate solution.
E. 1.2. THE UNCONSTRAINED CASE - VARYING AND p We can certainly tune the parameters to make our competitors (VFRBS) and (VEG) work better. However, such parameter configurations are far from satisfying the conditions of their theoretical results. For example, if we set p v then both VFRBS and VEG work better. In particular, if 5000, then get p suggested value p-
0.28, which is several times larger than its V
20
2 x 10-4
Let us further experiment other choices of parameters (i.e. the mini-batch size 6 and the probability of flipping a coin) to observe the performance of these algorithms.
(a) Larger Figure reveals the performance of these algorithms when we increase the mini-batch size to a larger value
29

---

Variance-Reduccd
Splitting Methods
Experiment 5000
Experiment 10000 and p 400
OG
10P
Relabve [02"1/100]
VFR-swig LVFR-sung
Relative (G.M)//(G)]
VER-swig
LVFR-owi
saga
VEG
VEG
10-1
10-1
20
40
61
100
20
40
80
100
epochs
spachs
Figure 5: Performance of algorithms to solve (37) on experiments when K
b 10.1n], while keeping the probability p unchanged.
Experiment
Experiment
p
Relaove [[22+ M|G.
10-1
Reletive operator norm [Gr+ [/|G:/|
10-1
10-3
10
10
10
OG
10
06
VFR-swig
10-11
VFR-avrg
LVER.wng
LVER
10-11
VFR-saga
10
VFR-4age
10-11
10-15
VEG
VEG
20
10
80
80
200
201
40
80
80
100
spochs
Number
opachs
Figure 6: Performance of algorithms for a large [0.1n] and unchanged p 173
Note that for n 5000, we have = 500 and p 0.058, and for n 10000, we have b = 1000 and p 0.046. With these large mini-batches, our algorithms still outperform other methods, while VERBS and VEG are significantly slowed down. The double-loop variant of (VFR) with SVRG performs best, while LVER-SVTg and VFR- have a similar performance.
(b) Medium b and larger p. Next, we b medium size of b |0.05n| (corresponding 250 for 5000 and b = 500 for 10000) and increase p (corresponding p 0.119 for 5000 and p 0.1 for 10000).
Then, the results shown in Figure 7.
Experiment 1: 5000 and 200
Experiment 2: 10000 and 400
10-1
Relative operator norm [G**)/))??
10-1
Reistive [Gu*1/]
10
10
10
10-9
10-1
OG
10
10-11
VFR-weg
11
VFR-oveg
LVFRisung
10-13
VFR
10
VFR soga
VFRB5
VFRBS
VEG
VEG
Number of epochs
Number of epochs
Figure 7: Performance of 6 algorithms for medium 6 [0.05m] and larger p 774
30

---

Variance-Reduced -eae---eeeeeeeeeceeeeee Splitting Methods
Then, observe that LVER- svrg and vaga superiorly outperform the others. The performance of the double-loop VFR-svrg is still similar the previous tests since it is not affected by p. In addition, VEG is now comparable with OG, but VERBS remains the slowest one.
(c) Large b and small To see the effect of on our competitors: VFRBS and VEG, as suggested by their theory, we decrease to 17 (corresponding to 0.014 for n 5000 and = 0.01 for n = 10000) and still set b = [0.1n], and the results are plotted Figure 8.
Experiment 5000 and 200
Experiment 10000 and 400
10'
Relative operaor |G/41/|G2"
Relative operator norm [[=*1//G)))
10-1
10-
10-1
10-1
10-4
10-3
DG
10
10-113
VERssng LVFR sung
10-11
VER-sag
10-11
LVER svng
10-11
VFR-snga
10-14
VERBS
VERBS
10-13
VEG
VEG
100
100
Number of spochs
Number of spachs
Figure 8: Performance of algorithms for large b [0.1ml and small p
As we can observed from Figure our methods highly outperform VERBS and VEG, suggesting that these competitors require larger probability select the snap-shot point w for full-batch evaluation. This is certainly not suggested in their theoretical results.
1.3. THE CONSTRAINED CASE
Next, we choose f(0) 61 --rrij (0) and ( == 61-r.)?= 22 (8) as the indicators of the lo -alls of radius r 5, respectively, In this case, implement three variants of (VFRBS): the double-loop (VFR-svrg), the loopless (LVFR-SVrg) and the SAGA (VFR-sad saga) variants solve (NI) and compare against 3 algorithms in the unconstrained case. Using the same data generating procedure as in the unconstrained case, we obtain the results as shown in Figure when K =
Experiment
200
Experiment 10000 p
Relabive |/>11/||aM
10"
Reaaiie operator norm [G**)/||a]]
10
10
10-7
10-7
10
10-9
10-11
VFR-swrg
10-11
VERcovg
LVER ovzz VFR-saga
LVER
10-13
VFR-saga
VFRBS
10-11
10
VEG
VEG
20
40
60
80
100
40
60
80
100
epochs
Number
epachs
Figure 9: Comparison of algorithms to solve constrained instances (37) on experiments when I [[Te average of 10 runs).
We see that the two SVRG variants of our (VFRBS): VFR-svig and LVFR-svig, as well our VFR-saga variant remain working well compared to other methods. They superiorly outperform the three competitors.
Finally, we test our methods and their competitors for the case K & as we done in Figure 5. Our results are plotted in Figure 10, where we observe a similar behavior as in Figure 5.
31

---

Variance-Reduced
Splitting Methods
Experiment 5000
Experiment 10000 and p 400
OG
10"
Relabve operator [a2*1//00]]
VFR-swig LVFR-oung
Relative [G*M)//(c+]
LVFR-ovig
saRa
saga
VEG
VEG
10"3
20
40
100
20
43
80
1001
epochs
spachs
Figure 10: Comparison of algorithms to solve constrained instances of (37) on experiments when K o (hhe average of 10 runs).
E.2. Nonconvex-Nonconcave Quadratic Minimax Problems We recall the noncomvex-noncconeay quadratic minimax optimization problem (15) in this subsection:
min max (C(u, = p(u) 1 wERPI vER*2
u A,u + u² uTuT x² Bi + buTu cq - (())).
(38)
where A; E X and B, E RP2Xp= are symmetric matrices, L, E RPiXPs, b, E IRP, C E IRp, and = sb., and
/ = ريحة are the indicator of standard simplexes in R and RR2,, respectively.
Let us first define x = [u, v] e R as the concatenation of the primal and dual variables and v, where p = Pt + P2. Next, we define
A Lj
u
b;
Auu + L2U bb
04p
G = G + =
+
and T=
-1 B
v
C;
Lau + B,U + t
Oub
A; Li
b
Then, we denote G; =
and gi:
Clearly, Gi(-) is an affine mapping from RP to R", but G2 is
L, B
C3
nonsymmetric. Let Gr = EE G E G.) + -E = Gz + g, where G = Σ G; and := - +
gi. Then, the optimality condition of (38) becomes 0 E + which is exactly in the form (NI). Clearly, A, and/or B, are not positive semidefinite, then (38) possibly covers nonconvex nonconccvvv minimax optimization instances.
1. 1. THE UNCONSTRAINED CASE
We consider the case 0 and ui = 0, leading to an unconstrained setting of (38), i.c. T = 0 as considered in (15) of the main text. Hence, the optimality condition of (38) reduces to Ga: = 0, which is of the form (NE).
(a) How generate data? To run our experiments, we generate synthetic data as follows. First, we fix dimensions p and P2 and the number of components n. We generate A, Q.D.Q for a given orthonormal matrix Qi and diagonal matrix D; diag( D!.... ..... D where its elements negative entries as max (D(, e} for j- 1,... PPI and e := -.... This choice of A; guarantees that A; is symmetric, but possibly not positive semidefinite. The matrix B, also generated by the same way. The pay-off matrix L; an p1 p2 matrix, which is also generated from the standard normal distribution for all E [n]. The vectors b, and c; are generated from the standard normal distribution for i E [n).. With this data generating procedure, G; is not symmetric and possibly not positive semidefinite.
generated from standard normal distribution and clipped its
(b) Algorithms. We again test 6 algorithms: two variants (double-loop SVRG VFR-svrr) and (loopless SVRG (VERg of (VFR), our (VFR) with SAGA estimator (VFR-saqa), VERBS from (Alacaoglu et al., 2023), VEG from (Alacaoglu & Malitsky, 2022), and OG (the standard optimistic gradient method), e.g., from (Daskalakis al., 2018).
32

---

Variance-Reduced Cerward-Refleceed Backterd Splitting Methods
(c) The details of Subsection 5.1 in Section 5. First, we provide the details of Subsection 5.1 in Section 5. The purpose of this example is to verify our theoretical results stated in Corollaries 3.2 and 3.3.
For the SVRG estimator, let us first choose " : 0.75. b = |n3/3|, and := wt7a as suggested by Corollary 3.2. Then, directly compute 77 := LVM where A == 6.25:2-3pp) 3p) 125p² and M 2.375 + A. Clearly, if 5000, then 6p
1
0.146153 Alternatively, if n = 10000, then 0.148934 1. Figure 1.
These learning rates are used in our experiments plotted in 1.
Similarly, for the SAGA estimator, we also choose ? == 0.75 and b == |n2/s] In this case, by Corollary 3.3, we can also 0.145693 0.146153 Alternatively, if n 10000, then 7 L These
directly compute 17 := LVM If n 5000, then 1) L
learning rates are used in VaR-saça.
Note that since the theoretical value of p in VFRBS and VEG is too small, we instead choose p == w.rh and also b = [ as in our methods. Then, we compute the learning rate n of these methods based on the formula given in (Alacaoglu et al., 2023) for VFRBS and (Alacaoglu & Malitsky, 2022) for VEG, respectively.
(d) Results for different set of parameters. Unlike Subsection 5.1 in the main text, we choose the parameters for these algorithms as in Subsection E.1. The 6 algorithms are run on 2 experiments. The first experiment is with n = 5000 and pi = P2 = 50, while the second one is with n = 10000 and p1 P2 = 100. These experiments are run 10 times, corresponding to 10 problem instances, and the average results are reported in Figure 11 in terms of the relative operator norm ||G**I| ||/ against the number of epochs.
Experiment 1: n 5000 and 100
Experiment 2: 10000 and p 200
Relative norn |Gx [
10-1
Relative operator norm [Gx"]/XX"
10-1
10-3
10-3
10
10
10-7
10-7
10
0G
10
06
10-11
VFR svrg LVFR-wing
10-11
VFR sng
IMFR-we
10-13
VER-aaga
10-13
VERaaaaa
VFRBS
VERBS
10-15
VEG
10-15
10
30
40
50
10
20
50
Number of opachs
Number of spochs
Figure 11: The performance of algorithms to solve the unconstrained case of (38) on 2 experiments (The average of 10 runs).
Clearly, under this configuration, both SVRG variants of our methods work well and significantly outperform other competitors. The loopless SVRG variant (VFR-svrg) of (VFR) seems to work best, while our VFR-saga has a similar performance as VEG. We also see that VFRBS has a similar performance as OG.
To improve the performance of these competitors, especially, VERBS and VEG, one can tune their parameters as in Subsection E. where the probability p of updating the snapshot point w is increased. However, with such a choice of p, its value is often greater or equal 0.5, making these methods be closed to deterministic variants. Hence, their theoretical complexity bounds are no longer improved over the deterministic counterparts.
E.2.2. THE CONSTRAINED CASE
We conduct two more experiments for the constrained case of (38) as in the main text when u E Ap, and v E App. where Ap = {u E R. D u, 1} the standard simplex in IR..
We run 6 algorithms for solving the constrained case of (38) using the same parameters Subsection 5.1, but with larger problems. We report the relative norm of the FBS residual || || (,, against the number of epochs. The results are revealed in Figure 12 for two datasets (p. n) (500, 5000) and (p. n) (300, 10000).
With these two additional experiments, both SVRG variants of our method (VFRBS) again work well and significantly outperform other competitors. The loopless SVRG variant (VER-svig) of (VFRBS) tends to work best, while our VFR-saga has relatively similar performance as VEG. We also see that VFRBS has similar performance trend as OG,
33

---

Variance-Reduced Forward-Reflectedd -eee--eeececceecceeeeee Backward Splitting Methods
Experiment 1: 5000 and p 500
Experiment 2: 10000 and p 300
Relative operator norm [Gx* |Gx* 1/1
10-1
Relative operator norm |Gx* [[Ixx/||xx||
10
VFR-swig
VFR-svig
LVFR-svig
LVFR-svyn
10
VER-saga
10
VFR-saga
VERBS
VERBS
10
VEG
10
VEG
10-
10-10
10-13
13
10-16
10
100
100
Number of epochs
Number of epochs
Figure 12: The performance of 6 algorithms to solve the constrained case of (38) on 2 experiments (The average of 10 runs).
but works better.
E.3. The & Re-Reuulrr Logistic Regression with Ambiguous Features
This Supp. Doc. provides the details of Subsection 5.2 Section the main text.
(a) Model. We consider standard regularized logistic regression model associated with given dataset {(X,, 31) H where X, is an i.i.d. sample feature vector and y; E {0, 13 is the associated label of X, Unfortunately, X; ambiguous, ie. belongs one of possible examples { Since do know X, evaluate the loss, we consider the worst-case loss f((w) maxigism (((xij, w), yt) computed from examples, where {(T, s) = log(1 + exp(r))
the standard logistic loss.
Using the fact that maxisssm (-)) ma En= (-)((), where Am is the standard simplex in RR we can model this regularized logistic regression into the following minimax problem:
min max [(w, 2) == Σ Σラ= 1 23 (Xi,, w) 9i) + (w) - سسmm (=) WER"
(39)
where ((T,s) s) log(1 + exp(r)) sT is the standard logistic loss, R(w)= ||l|| is -nomm regularizer, is regularization parameter, and dA.,, is the indicator of Am that handles the constraint Am. This problem is exactly the one stated in (16) of the main text.
First, let us denote x == [w: 2] E IR as the concatenation of w and z with = d + m, and
EØ ((Xij, w), Ni)
{(Xn, w).
and Tx TdR(w)
Gą =
055m (2)
-E((Xim, w),yt) y:)
where (T(T,, 8) exp(r)
[+exp(7) 8. Then, the optimality condition of (39) can be written as (NT): 0 E Gr + Tæ, where Gr n Σm GeD.
(b) Input data. We test our algorithms and their competitors two real datasets: 9a (134 features and 3561 samples) and w8a (311 features and 45546 samples) downloaded from LIBSVM (Chang & Lin, 2011). For a given nominal dataset {(8,, 9:) } we first normalize the feature vector X, such that its column norm one, and then add a column all ones to address the bias term. To generate ambiguous features, we take the nominal feature vector X; and add random noise generated from a normal distribution of zero mean and variance of a = 0.5. In our test, we choose 7 := 10 10-3 and m = 10 for all the experiments.
(c) Algorithms. before, we implement variants of our method (VFRBS): VFR-svrg, LVFR-svrg, and VFR-saga
34

---

Variance-Reduced Beaaadd-Refeetee-Baaaete Splitting Methods
solve (39). We also compare them with OG, VFRBS, and VEG. We choose x 0.5 ones(p) in all experiments. We all the algorithms for 100 epochs and report the relative FBS residual norm 11,1l against the epochs.
(d) Parameters. Since it is very difficult to estimate the Lipschitz constant L of G, we are unable to set correct learning rate 17 in the underlying algorithms. We instead compute an estimation L := || and then set 1J 2. by tuning w for each algorithm. More specifically, after tuning, we obtain the following configuration.
For the three variants of (VERBS): VER- svrg, LVER-SVEE and VFR-saga, we set vj = L for w8a.
25 for a a and rj = L
50 for a9a and 1j 100 for w8a. For OG, set
7.
For VFRBS, we choose 47.5(1 w- for a and 1] = 95(1
for wBa.
21
2L
For VEG, we select 1) = 47.5vI-a for a and 7j 95vT-a for w wBa with a == 1 p.
i
We still choose the mini-batch size b and the probability p of updating the snapshot point w in SVRG variants as b [0.5n2/33 and n n1/33, respectively for all the algorithms.
We conduct two more experiments using the well-known MNIST dataset (n 70000 and = 780) where we want to classify the even and odd numbers into two different classes, respectively. We use the same parameter selection as the experiment with the a9a dataset. In the first experiment, we choose m = 10 and the variance of noise o = 0.3. In the second experiment, we choose m and the variance of noise 0.25. We only run algorithms and leave out the VFRBS method since have not managed find the parameters that make work stably. The result of this experiment shown in Figure 13.
The MNIST Dataset: (n, p. m) (70000, 10)
The MNIST Dataset: (n, m) (70000, 801, 20)
10a
OG
Relative operator norm ||Gx** G/
10a
OG
Relative operator norm ||,,x**/
VFR-svrg
VFR-svrg
LVFR--vrr
LVFR-svig
10-1
VFR-saga VEG
10-1
VFR-saga
VEG
10
10-2
10
10
60
100
80
100
Number epochs
Number of epochs
Figure 13: Comparison of 6 algorithms to solve (16) on the real dataset: MNIST.
We from Figure 13 again that three variants VFR-svrg, LVFR-SvIR and VFR-saga have similar performance and are much better than their two competitors. Here, VEG is better than OG, but both methods are slower than ours.
References
Adly, S. and Attouch, H. First-order inertial algorithms involving dry friction damping. Math. Program., pp. 1-41, 2021. Alacaoglu, A. and Malitsky, Y. Stochastic variance reduction for variational inequality methods. In Conference Learning Theory, pp. 778-816. PMLR, 2022.
Alacaoglu, A., Böhm, A., and Malitsky, Y. Beyond the golden ratio for variational inequality algorithms. J. Mach. Learn Res, 24(172):1-33, 2023.
Arjovsky, M., Chintala, S., and Bottou, I.. Wasserstein generative adversarial networks. In International Conference on Machine Learning, pp. 214-223, 2017.
Attouch, H. and Cabot, A. Convergence of a relaxed inertial proximal algorithm for maximally monotone operators. Math. Program., 184(1):243-287 2020.
35

---

Variance-Reduced BeraarddRefleeeed-Baaaaae Splitting Methods
Bauschke, H. H., Moursi, W. M., and Wang, X. Generalized monotone operators and their averaged resolvents. Math. Program., pp- 1-20, 2020.
Beznosikov, A., Gorbunov, E., Berard, H., and Loizou, N. Stochastic gradient descent-ascent: Unified theory and new efficient methods. In International Conference on Artificial Intelligence and Statistics, pp. 172-235. PMLR, 2023.
Böhm, A. Solving nonconvex-nonconcaay min-max problems exhibiting weak Minty solutions. Transactions on Machine Learning Research, 2022.
Bot, R. I., Mertikopoulos, P., Staudigl, M., and Vuong, P. T. Minibatch stochastic variational inequalities. Stochastic Systems, 11(2):1 12-113, 2021.
forward methods for solving
Cai, X., Song, C.,, Guzmán, C., and Diakonikolas, J. A stochastic Halpern iteration with variance reduction for stochastic monotone inclusion problems. In Proceedings of the 12th International Conference on Learning Representations (ICLR 2022), 2022. URL https:://///m. http::::///ppneeeeeeeee www.://wnwneeooooeeeeeeeeeeeeeee wiw.amodddd----ii-ee-eoooooooose 8TPCf.
Cai, X., Alacaoglu, A., and Diakonikolas, J. Variance reduced halpemn iteration for finite-sum monotone inclusions. In The 12th International Conference on Learning Representations (ICLR), pp. 1-33, 2024.
Chakrabarti, D., Diakonikolas, J., and Kroer, C. Block-coordinate methods and restarting for solving extensive-form games. Advances in Neural Information Processing Systems, 36, 2024.
Chang, C.-C. and Lin, C.-J. LIBSVM: library for Support Vector Machines. ACM Transactions on Intelligent Systems and Technology. 2:27:1-27:20 2011.
Combettes, P. L. and Eckstein, J. Asynchronous block-iterative primal-dual decomposition methods for monotone inclusions. Math. Program., 168(1))6456672, 2018.
Combettes, P. and Pesquet, J.-C. Stochastic quasi-Fejér block-coordinate fixed point iterations with random sweeping. SIAM Optim., 25(2):1221-12488 2015.
Cui, S. and Shanbhag, U. On the analysis of variance-reduced and randomized projection variants of single projection schemes for monotone stochastic variational inequality problems. Sei-Valued and Variational Analysis, 29(2)4553-49, 2021.
Cutkosky, A. and Orabona, Momentum-baat based variance reduction in non-convex SGD. In Advances in Neural Information Processing Systems, pp- 15210-15219, 2019.
Daskalakis, C., Ilyas, A., Syrgkanis, V., and Zeng, H. Training GANs with Optimism. In International Conference Learning Representations (ICLR 2018), 2018.
Defazio, A., Bach, F,, and Lacoste-Julien, S. SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives. In Advances in Neural Information Processing Systems (NIPS), pp. 1646-1654, 2014.
Diakonikolas, J., Daskalakis, C.. and Jordan, M. Efficient methods for structured nonconvex-noncoonan min-max PMLR, 2021.
optimization. In International Conference Artificial Intelligence and Statistics, pp.
Driggs, D., Liang, J., and Schönlieb, C.-B. On biased stochastic gradient estimation. Journal of Machine Learning Research, 23(24):1-43, 2022.
Gorbunov, E., Berard, H., Gidel, G., and Loizou, Stochastic extragradient: General analysis and improved rates. In International Conference on Artificial Intelligence and Statistics, pp. 7865-7901. PMLR. 2022a.
Gorbunov, E., Loizou, N., and Gidel, G. Extragradient method: (l1/k) last-iterate convergence for monotone variational inequalities and connections with cocoercivity. In International Conference on Artificial Intelligence and Statistics, pp. 366 402. PMLR, 2022b.
Gower, R. M., Richtárik, P., and Bach, F. Stochastic quasi-gradient methods: Variance reduction via Jacobian sketching. Math. Program., 188(1):135-199) 2021.
36

---

Variance-Reduced Beserdd-Refeeeeed-Baeeae Splitting Methods
Hamedani, E. Y., Jalilzadeh, A., Aybat, N. S., and Shanbhag, U. V. Iteration complexity of randomized primal-dual methods for convex-concave saddle point problems. arXiv preprint
2018.
Hanzely, F., Mishchenko, K., and Richtárik, SEGA: Variance reduction via gradient sketching. In Advances Neural Information Processing Systems, pp. 2082-2093, 2018.
Horváth, S., Kovalev, D., Mishchenko, K.. Richtárik, P., and Stich, S. Stochastic distributed learning with gradient quantization and voriinc- reduction. Optimization Methods and Software, 38(1):91-106, 2023.
Hsieh, Y.GG., lutzeler, F., Malick, J., and Mertikopoulos, P. On the convergence of single-call stochastic extra- gradient methods. In Advances in Neural Information Processing Systems, pp.
2019.
lusem, A. N., Jofré, A., Oliveira, I., and Thompson, P. Extragradient method with variance reduction for stochastic variational inequalities. SIAM J. Optim., 27(2):686-724, 2017.
Johnson, R. and Zhang, Accelerating stochastic gradient descent using predictive variance reduction. In NIPS, pp. 315-323, 2013.
Juditsky, A., Nemirovski, A., and Tauvel, C. Solving variational inequalities with stochastic mirror-prox algorithm. Stochastic Systems, 1(1):17-58, 2011.
Kannan, A. and Shanbhag, U. V. Optimal stochastic extragradient schemes for pseudomonotone stochastic variational inequality problems and their variants. Comput. Optim. Appl., 14(3):779-820 2019.
Konnov, I. Combined relaxation methods for variational inequalities. Springer-Verlag, 2001.
Kotsalis, G., Lan, G., and Li, T. Simple and optimal methods for stochastic variational inequalities, i: operator extrapolation.
SIAM J. Optim., 32(3):0041-2075 2022.
Li, C. J., Yu, Y., Loizou, N., Gidel, G., Ma, Y., Roux, N. L., and Jordan, M. On the convergence of stochastic extragradient for bilinear games using restarted iteration averaging. In International Conference on Artificial Intelligence and Statistics, pp. 9793-9826. PMLR, 2022.
Loizou, N., Berard, H., Gidel, G., Mitliagkas, I., and Lacoste-Julien, S. Stochastic gradient descent-ascent and consensus optimization for smooth games: Convergence analysis under expected co-coercivity. Advances in Neural Information Processing Systems, 34:19095-19108, 2021.
Luo, Y. and Tran-Dinh, Q. Extragradient-type methods for co-monotone root-finding problems. (UNC-STOR Technical Report), 2022.
Malitsky, Y. and Tam, M. K. A forward-backward splitting method for monotone inclusions without cocoercivity. SIAM J. Optim.. 30(2):1451-1477) 2020.
Mishchenko, K., Kovalev, D., Shulgin, E., Richtárik, P., and Malitsky, Y. Revisiting stochastic extragradient. In International
Conference on Artificial Intelligence and Statistics, pp. 4573-4582. PMLR, 2020.
Noor, M. A. Extragradient methods for pseudomonotone variational inequalities. J. Optim. Theory Appl., 117(3):475-488, 2003.
Noor, M. A. and Al-Said, E. Wiener-Hopf equations technique for quasimonotone variational inequalities. J. Optim. Theory Appl., 103:705-714, 1999.
Peng, Z., Xu, Y., Yan, M., and Yin, W. ARock: algorithmic framework for asynchronous parallel coordinate updates.
SIAM J. Scientific Comput., 3855128851 2016.
Pethick, T., Fercoq, O., Latafat, P., Patrinos, P., and Cevher, V. Solving stochastic weak Minty variational inequalities without increasing batch size. In Proceedings of International Conference on Learning Representations (ICLR). pp. 1-34, 2023.
Pham, H. N., Nguyen, M. L., Phan, T. D., and Tran- Tran-Dinh, Q. ProXSARAH: An efficient algorithmic framework for stochastic composite nonconvex optimization. Mach. Learn. Res., 21:1 21:1-48, 2020.
37

---

Variance-Reduced Beeerd-RReleeeed-Baaaaaee Splitting Methods
Robbins, H. and Monro, S. A stochastic approximation method. The Annals of Mathematical Statistics, 22(3):400-077, 1951.
Rockafellar, R. and Wets, R. J.-B. Variational Analysis. Springer-Vetee 1997.
Song, C. and Diakonikolas. J. Cyclie coordinate dual averaging with extrapolation. SIAM J. Optim., 33(4):23352261,, 2023.
Song, C., Zhou, Z., Zhou, Y., Jiang, Y., and Ma, Y. Optimistic dual extrapolation for coherent non-monotone variational inequalities. Advances in Neural Information Processing Systems, 33:14303-14314, 2020.
Tran-Dinh, Q. Sublinear Convergence Rates of Extragradient-Type Methods: A Survey on Classical and Recent Develop-
ments. arXiv preprint arXiv:3303.7192,, 2023.
Tran-Dinh, Q. and Luo, Y. Randomized Block-Coordinate Optimistic Gradient Algorithms for Root-Finding Problems.
Math. Oper. Res., in press, 2025.
Tran-Dinh, Q. and Nguyen-Trung N. Guriance-Reduced Accelerated Better Convergence Guarantees. Tech. Report UNC-STOR, 2025a.
Methodsfor Generalized Equations:
Tran-Dinh, Q. and Nguyen-Trung, N. Accelerated extragradient-type methods-par 2: Generalization and sublinear convergence rates under co-hypomonotonicctyy arXiv preprint arXiv:2501.0458: 2025b.
Tran-Dinh, Q., Pham, N. H., Phan, D. T., and Nguyen, L. M. A hybrid stochastic optimization framework for stochastic composite nonconvex optimization. Math. Program., 191:1005-100 1001, 2022.
Tu, V. On the weak convergence of the extragradient method for solving pseudo-monotone variational inequalities. J. Optim. Theory Appl., 176(2):399-40) 409399999999220 2018.
Yang. J., Kiyavash, N., and He, N. Global convergence and variance-reduced optimization for a class of nonconvex- nonconcave minimax problems. arXiv preprint arXiv:2002.09621, 2020.
Yousefian, F., Nedić, A., and Shanbhag, U. V. On stochastic mirror polror-ppox algorithms for stochastic cartesian variational inequalities: Randomized block coordinate and optimal averaging schemes. Ser- Valued and Variational Analysis, 26: 789-819, 2018.
Zhou, D., Xu, P., and Gu, Q. Stochastic nested variance reduction for nonconvex optimization. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 3925-399 3936. Curran Associates Inc., 2018
38


Reviews:
[{"review_id": "wLjbdNtlmM", "reviewer": "Reviewer_bar9", "summary": "In this paper, the author proposes two novel algorithms for solving a class of non-monotone equations, building upon the Forward-Reflected Backward Splitting framework and incorporating variance reduction techniques such as SVRG and SAGA. The proposed methods are accompanied by rigorous convergence guarantees and achieve state-of-the-art complexity bounds. To validate the theoretical results, the authors conduct a series of numerical experiments demonstrating the practical effectiveness of the algorithms.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "First of all, we highly acknowledge the reviewer for comments and feedback. Below is our response to each point.\n\nC1: In this paper, the author proposes two novel algorithms for solving a class of non-monotone equations, building upon the Forward-Reflected Backward Splitting framework and incorporating variance reduction techniques such as SVRG and SAGA. The proposed methods are accompanied by rigorous convergence guarantees and achieve state-of-the-art complexity bounds. To validate the theoretical results, the authors conduct a series of numerical experiments demonstrating the practical effectiveness of the algorithms.\n\n>R1: Thank you for a great summary of our work and contribution. \nUnlike existing variance-reduction methods in the literature which primarily address co-coercive or monotone problems, we instead consider a class of problems satisfying a weak-Minty condition, which is possibly non-monotone. Moreover, our methods are new and significantly different from the existing ones due to the use of $S^k$. The complexity $\\mathcal{O}(n + n^{2/3}/\\epsilon)$ is the best-known (but not optimal) for SVRG and SAGA in nonconvex optimization (without using enhancement tricks, like researt or nested loops).\n\nC2: On lines 140–142 of the second column, please ensure that the citation of “Diakonikolas et al.” is formatted correctly in accordance with the appropriate referencing style.\n\n>R1: Thank you. We will correct the issue you identified.", "claims_and_evidence": "The claims in the submission are well-supported by rigorous theoretical analysis and empirical validation, with no evident inconsistencies.", "methods_and_evaluation_criteria": "Yes, the proposed methods and evaluation criteria are appropriate and well-aligned with the problem and application context.", "theoretical_claims": "Yes, I reviewed the proofs and identified issues related to the assumptions and justification of some theoretical claims.", "experimental_designs_or_analyses": "I find the experiments to be well-designed and informative, with no apparent issues in their methodology or analysis.", "relation_to_broader_scientific_literature": "The proposed approach and novel algorithms can be helpful for optimization community to create new algorithms or can be useful for practical tasks.", "essential_references_not_discussed": "No, all essential related works appear to be appropriately cited and discussed.", "other_strengths_and_weaknesses": "The manuscript is well-composed, and the results are presented in a clear and coherent manner.", "other_comments_or_suggestions": "On lines 140–142 of the second column, please ensure that the citation of “Diakonikolas et al.” is formatted correctly in accordance with the appropriate referencing style.", "questions_for_authors": "I do not have any questions to raise at this point.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TWWDwgjpZu", "reviewer": "Reviewer_tBMV", "summary": "This paper studies the forward-reflected operator in two types of variance-reduced estimators: SVRG and SAGA. Using these estimators, the authors propose the Variance-Reduced Forward-Reflected Method and the Variance-Reduced Forward-Reflected-Backward Splitting, which solve nonlinear and generalized equations, respectively. Through Lyapunov analysis, the authors establish an complexity of $O(n+n^{2/3}\\epsilon^{-2})$ to obtain $\\epsilon$-residual error, matching the best-known results. Theoretical findings are further validated through numerical experiments.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "First of all, we highly acknowledge the reviewer for constructive comments and feedback. Below is our response to each point.\n\nQ1: The authors ...  results. I request that authors verify this by providing a table or a comparison with prior works in terms of complexity and assumptions. This is necessary and would enhance the novelty and clarity.\n>R1: Thank you! Yes, we will provide a table to summarize and compare our work with existing results. We will also highlight our novelty using this summarized table as suggested.\n\nQ2: This paper systematically studies variance-reduced estimators, SVRG and SAGA estimators of forward-reflected operator and proposes new algorithms with theoretical complexity matched to the best known results. The paper is well organized, and the experiments demonstrate the theoretical results. Furthermore, it shows that proposed algorithms perform better than ones in prior work. I believe this work provides a valid contribution to the literature on stochastic approximation and nonmonotone inclusion problems, since, although their theoretical results do not improve upon prior work, this paper introduces an interesting new operator with symmetrical variance reduction analyses and demonstrates its practical potential.\n>R2: We appreciate your adequate summary of our paper and contribution. The complexity we achieved in this paper is similar to the one in variance-reduced methods for nonconvex optimization using SVRG and SAGA, and this complexity is unimprovable without using enhancement tricks such as restart or nested loop. \n\nQ3: Do the complexities of the proposed methods in this paper indeed match the best-known results? \n>R3: To our best knowledge, our oracle complexity results are the best-known so far for non-monotone problems (i.e., covered by Assumptions 1.3-1.4) using SVRG and SAGA. We will add a table to compare as suggested.", "claims_and_evidence": "The authors state that their results match the best known results. I request that authors verify this by providing a table or a comparison with prior works in terms of complexity and assumptions. This is necessary and would enhance the novelty and clarity.", "methods_and_evaluation_criteria": "Yes, methods and evaluation criteria makes sense.", "theoretical_claims": "I briefly reviewed the proofs in both the main text and the supplementary material.", "experimental_designs_or_analyses": "Analyses and experiments look sound.", "relation_to_broader_scientific_literature": "As the authors introduced, this work falls within the literature on stochastic approximation and nonmonotone inclusion problems.", "essential_references_not_discussed": "I believe this papers address prior works.", "other_strengths_and_weaknesses": "This paper systematically studies variance-reduced estimators, SVRG and SAGA estimators of forward-reflected operator and proposes new algorithms with theoretical complexity matched to the best known results. The paper is well organized, and the experiments demonstrate the theoretical results. Furthermore, it shows that proposed algorithms perform better than ones in prior work. I believe this work provides a valid contribution to the literature on stochastic approximation and nonmonotone inclusion problems, since, although their theoretical results do not improve upon prior work, this paper introduces an interesting new operator with symmetrical variance reduction analyses and demonstrates its practical potential.", "other_comments_or_suggestions": "Please refer to the Questions for Authors section.", "questions_for_authors": "Do the complexities of the proposed methods in this paper indeed match the best-known results? (As I wrote in Claims and Evidence section, please make table or a comparison with prior works in terms of complexity and assumptions.)", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nLxTc2bjrE", "reviewer": "Reviewer_9taN", "summary": "1. inspired by  SVRG & SAGA, construct two variance-reduced estimators for the forward-reflect operator\n2. show that VFR and VFRBS methods achieve SOTA oracle complexity for non-monotone operator splitting problems", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "First of all, we highly acknowledge the comments and questions from the reviewer. \nBelow is our detailed response.\n\nQ1. \"Does the convergence ... estimator? If I replace SVRG/SAGA ...  converge?\n>R1. The answer is \"no\". In lines 233-240 and 297-304, we have stated that any estimator $S^k$ satisfies Definition 2.1. can be used in our methods, covering a wide class of unbiased estimators. The convergence guarantees are stated in Theorems 3.1 and 4.1 without specifying $S^k$. Therefore, any estimator satisfying Definition 2.1. can be used in our methods to have convergence guarantees as in these theorems. SVRG and SAGA are two concrete instances of $S^k$. However, the oracle complexity depends on the specific estimator since we do need to know the cost of constructing such an estimator (See Subsections 3.2 and 4.2).\n\nQ2. \"Could you provide a table to compare with the existing literature, in terms of problem setup and sample complexity?\"\n>R2.  Thank you for the request. Yes, we will provide a table summarizing existing results and comparing them with ours. \n\nQ3. \"The authors claimed ... very small.\"\n>R3. Our paper primarily focuses on new algorithms, theoretical convergence, and oracle complexity rather than experiments. Our experiments are used to validate the theory and we did not specifically focus on concrete ML applications. Nevertheless, based on your suggestion, we will increase the sizes of our experiments in the revision.\n\n>Regarding nonsmooth and non-monotone problems, Supp. E2 provides nonsmooth and non-monotone problems since the underlying matrix is not necessarily positive semidefinite (non-monotone), and they have constraints (nonsmoothness).\nExample 2 tackles a nonlinear minimax problem, which is nontrivial compared to the quadratic minimax one in Example 1. As mentioned above, we will substantially increase the size of our experiments as you suggested. \n\nQ4: \"Solving non-monotone problem ... are non-monotone.\"\n>R4: Example 1 allows $\\mathbf{G}$ to be non-positive semidefinite. Therefore, it violates the condition $\\langle Gx - Gy, x - y\\rangle \\geq 0$ (since $\\langle Gx - Gy, x - y\\rangle = (x-y)^T\\mathbf{G}(x - y)$). \nThis means that $G$ is non-monotone. We will clearly explain this in the revision.\n\nQ5: \"I understand that it is hard to do a theoretical comparison ... even diverge behavior.\"\n>R5: The work by (Cai et al 2023) relies on Halpern's fixed-point iteration, and indeed requires stronger assumptions than ours. We will conduct numerical examples to compare both methods, at least in the co-coercivity case. We are not sure if it will diverge in the non-monotone case, but there is no theoretical result to guarantee convergence for this method under Assumption 1.4. \n\nQ6: Weakness:\n-- very dense theory and very weak experiments\n-- the motivating ML applications, e.g. GAN and adversarial training, are not the main line for generative modeling and robust deep learning.\n>R6.1: We agree with the reviewer that the paper is dense, and we will try to improve its presentation as suggested. It was significantly challenging for us to conduct a 10-page paper with both strong theories and experiments. Therefore, we chosed to focus on methodology and theoretical results in our submission. However, as we discussed above, additional experiments will be added to the supplementary in the revision.\n\n>R6.2: We will rephrase some inaccurate statements related to robust deep learning and generative modeling as suggested.\n\nQ7: \"I think this paper is too dense ... applications in ML. \nFor example, what kind of non-convex problem satisfies Assumption 1.4?\"\n>R7: We will improve the presentation of our paper as suggested. We will add more experiments to the second example to handle nonconvex regularizers such as CAD, to obtain non-monotone problems. \nRegarding Assumption 1.4., one simple example is to consider a nonconvex quadratic problem $\\min_x \\frac{1}{2}x^TQx + q^Tx$, where $Q$ is symmetric and invertible, but not necessarily positive semidefinite. The optimality condition is $Qx + q = 0$.\nIf we define $Gx = Qx + q$ and $T = 0$, then $\\langle Gx - Gy, x - y\\rangle = (x-y)^TQ(x-y) \\geq \\lambda_{\\min}(Q^{-1})\\Vert Q(x-y)\\Vert^2 =  \\lambda_{\\min}(Q^{-1})\\Vert Gx - Gy\\Vert^2$. Hence, if $\\lambda_{\\min}(Q^{-1}) < 0$, then it satisfies Assumption 1.4. with $\\rho := -\\lambda_{\\min}(Q^{-1}) > 0$.\nNote that one can verify that our examples in Section 5 and SupDoc. E satisfy Assumption 1.4. We will add a paragraph to explain this case.\n\nQ8: Overall Recommendation: 2: Weak reject (i.e., leaning towards reject, but could also be accepted)\n>R8: We hope the reviewer will read our responses above, and re-evaluate our work based on new algorithms (especially, the one for inclusion) and theoretical contributions. We promise to improve our paper and add a table for comparison, larger experiments, and a comparison to [Cai et al 2023] as requested.", "claims_and_evidence": "1. Does the convergence of your splitting algorithm have strong connection with your estimator? If I replace SVRG/SAGA with anothor variance-reduced estimators, will the algorithm converge? \n2. Could you provide a table to compare with the existing literature, in term of problem setup and sample complexity?", "methods_and_evaluation_criteria": "The authors claimed their work is to solve non-monotone, non-smooth, and large-scale problem arising from generative ML, adversarial learning, and robust learning. However, their numerical examples are limited to synthetic or toy example whose scale is very small.", "theoretical_claims": "I didn't check the proofs.", "experimental_designs_or_analyses": "1. Solving non-monotone problem is one of the main differences compared to the existing works. However, those examples didn't explain why they are non-monotone. \n2. I understand that it is hard to do theoretical comparison to those works requiring monotonicity assumption, e.g. [Cai et al., 2023]. Could you provide some empirical comparison? Since the monotonicity assumption is violated for those approaches, one can expect slow convergence or even diverge behavior.", "relation_to_broader_scientific_literature": "Exploring the application of operator splitting in non-monotone and stochastic setup could be useful for broader engineering problems than ML learning.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Strength:\n1. solid theoretic work\n\nWeakness:\n1. very dense theory and very weak experiments\n2. the motivating ML applications, e.g. GAN and adversarial training, are not the main line for generative modeling and robust deep learning.", "other_comments_or_suggestions": "I think this paper is too dense in technique presentation while leaving very little room for its possible applications in ML. For example, what kind of non-convex problem satisfying Assumption 1.4?", "questions_for_authors": "N/A", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HC3ubv0BsB", "reviewer": "Reviewer_Q5J9", "summary": "This paper proposes two stochastic variance-reduction algorithms to solve a class of nonmonotone equations. The key technical tool used in this paper is the intermediate object $S_{\\gamma}$. \nThe authors apply classical variance reduction techniques on $S_{\\gamma}$ instead of the operator $G$, and they show their proposed algorithms obtain convergence rates of $\\mathcal{O}(n+n^{2/3}\\epsilon^{-2})$.\n\n## update after rebuttal\nThe author's reply makes sense, so I keep my score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "First of all, we acknowledge the reviewers for his/her comments and feedback on our work. Below is our response to each point.\n\nP1. Weakness:\n + Q1.1: Paper organization is not good enough. The comparison with (Cai et al. 2023) can be presented in the main text instead. The comparison with existing work can be summarized in a table. The author may consider highlighting the co-coercive condition presented in (Cai et al. 2023) in the main text.\n> R1.1: Thank you for your comment. We will try to improve the organization of the paper as suggested. We will bring the comparison with [Cai et al 2023] back to the main text, and highlight the co-coercivity condition used in [Cai et al 2023] vs. our assumptions. This co-coercivity condition is stronger than monotonicity+Lipschitz, and quite limited since it does not cover bilinear matrix games. We will also add a table summarizing the comparison with existing works. \n+ Q1.2: The convergence rate of the proposed methods is inferior to the methods obtained in (Cai et al. 2023). Could the authors recover the convergence rate if the co-coercive assumption is imposed?\n> R1.2: In terms of convergence rates, our rate is slower than [Cai et al. 2023] ($1/k$ vs. $1/k^2$) when solving co-coercive equations. However, since our assumptions are significantly weaker, it is unclear if this theoretical comparison is fair.\nNote that Halpern's method in [Cai et al 2023] is a type of accelerated method, while our method is not. We believe that $O(1/k)$ is the best rate one can get for non-monotone problems, similar to nonconvex optimization methods, though we do not have a lower bound result to compare at this moment. In nonconvex optimization, this rate is optimal (up to a constant factor). \n+ Q1.3: In Theorem 3.1, the choice of $\\kappa$ seems restrictive.\n> A1.3: So far, the range of $\\kappa$ is indeed restrictive since we did not focus on optimizing $M$ and $\\gamma$ in Theorems 3.1 and 4.1. In deterministic methods, the largest range of $L\\kappa$ is $L\\kappa < 1$. We believe that our range on $L\\kappa$ can also be improved, but it requires substantial changes in the algorithmic design and analysis. \n+ Q1.4: In the experiments, the author may consider comparing with the methods proposed in (Cai et al. 2023).\n> A1.4: Yes, we will add some examples to compare with [Cai et al 2023] for co-coercive problems as suggested.\n\nP2: Other Strengths And Weaknesses:\n+ Q2.1: In both VFR-SVRG and VFR-SAGA algorithms, the minibatch size is quite large. It is of the order $\\mathcal{O}(n^{2/3})$, which is larger than most minibatch sizes of the stochastic methods solving nonconvex optimization.\n> A2.1: This choice of mini-batch (i.e., $b = \\mathcal{O}(n^{2/3})$) is to attain the complexity of $\\mathcal{O}(n + n^{2/3}\\epsilon^2)$. \nOur theory works with any batch size $b \\in [1, n]$. However, the choice of $b$ affects the final oracle complexity. Therefore, in practice, we can choose any batch size, and our method still converges, but the overall oracle complexity will be worse than the best one $\\mathcal{O}(n + n^{2/3}\\epsilon^2)$. We will add a brief discussion on this aspect in the revision.\n\n+ Q2.2: The choice of $\\gamma$ is quite restrictive. It must belong to the range of $(1/2, 1)$, so it cannot recover the classic FRBS and optimistic gradient methods. Could the authors provide some intuition as to why cannot be chosen for the proposed approaches?\n> A2.2: As can be seen from (22) of Lemma C.1., if $\\gamma = 1/2$, then we must impose $\\mu = 0$ and $\\kappa = 0$, leading to star-monotonicity in Assumption 1.4. Our theory still works, but no longer covers Assumption 1.4. with $\\kappa > 0$, but covers monotone+Lipschitz problems. As far as we know, existing FRBS methods require monotonicity, which is consistent with our theory when choosing $\\gamma = 1/2$. We will add a discussion of this special case.\n\nP3: Other Comments Or Suggestions:\n+ Q3.1: In Corollary 3.3, line 324, why there is a  in the denominator? The SAGA estimator no longer needs probability.\n>A3.1.: Thank you. It is indeed a typo. We will correct it.\n+ Q3.2: A table summarization of existing work with different assumptions and convergence criteria could be provided.\n>A3.2: We will add a table to summarize existing works with assumptions and convergence criteria as suggested.\n\nP4: Overall Recommendation: 3: Weak accept (i.e., leaning towards accept, but could also be rejected)\n>A4: We hope we have addressed your concerns. We will implement what we have promised above. \nGiven our novel methods and theoretical contributions, and the effort of addressing your concerns, your re-evaluation is highly appreciated.", "claims_and_evidence": "Strength:\n\n* The authors provide a theoretical analysis of their proposed approaches and show their algorithm enjoys a convergence rate of $\\mathcal{O}(n+n^{2/3}\\epsilon^{-2})$.\n\n* Applying variance reduction on the intermediate object $S_{\\gamma}$ is a novel and interesting idea.\n\nWeakness:\n\n* Paper organization is not good enough. The comparison with (Cai et al. 2023) can be presented in the main text instead. The comparison with existing work can be summarized in a table. The author may consider highlighting the co-coercive condition presented in (Cai et al. 2023) in the main text.\n\n* The convergence rate of the proposed methods is inferior to the methods obtained in (Cai et al. 2023). Could the authors recover the $\\mathcal{O}(n+n^{1/2}\\epsilon^{-2})$ convergence rate if the co-coercive assumption is imposed?\n\n*  In Theorem 3.1, the choice of $\\kappa$ seems restrictive.\n\n* In the experiments, the author may consider comparing with the methods proposed in (Cai et al. 2023).", "methods_and_evaluation_criteria": "See above.", "theoretical_claims": "See above.", "experimental_designs_or_analyses": "See above.", "relation_to_broader_scientific_literature": "It is an interesting work for solving large-scale nonmonotone equations. Overall, the idea of applying variance reduction on the intermediate object is novel and interesting.", "essential_references_not_discussed": "I think the most relevant references have been discussed/", "other_strengths_and_weaknesses": "Some other weaknesses include:\n\n* In both VFR-SVRG and VFR-SAGA algorithms, the minibatch size is quite large. It is of the order $\\mathcal{O}(n^{2/3})$, which is larger than most minibatch sizes of the stochastic methods solving nonconvex optimization.\n\n* The choice of $\\gamma$ is quite restrictive. It must belong to the range of $(1/2, 1)$, so it cannot recover the classic FRBS and optimistic gradient methods. Could the authors provide some intuition as to why $\\gamma$ cannot be chosen $1/2$ for the proposed approaches?", "other_comments_or_suggestions": "* In Corollary 3.3, line 324, why there is a $p$ in the denominator? The SAGA  estimator no longer needs probability $p$.\n\n* A table summarization of existing work with different assumptions and convergence criteria could be provided.", "questions_for_authors": "See above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: For monotone inclusion problems, the paper develops novel stochastic variance-reduction methods, which combine the ideas from the forward-reflected-backward splitting method and unbiased variance-reduced estimators. By proposing a new definition for variance-reduced estimators, a quite generic theoretical analysis is delivered. From the reviews, the reviewers reach a consensus that the paper contains interesting result and novel contribution to the community. In the meantime, the paper contain the following issues
 - Paper presentation. The reviewers raised several points to improve the readability of the paper.
 - Weak numeric experiments, and lacks comparison against some existing methods.

The authors are encouraged to modify the paper accordingly.
Decision: Accept (poster)