$\displaystyle\|X_{k+1}-p\|^{2}+(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda_{k})}-1\right)\|X_{k+1}-X_{k}\|^{2}-\alpha_{k}\|X_{k}-p\|^{2}$
$\displaystyle\leq\|X_{k}-p\|^{2}+(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda_{k})}-1\right)\|X_{k}-X_{k-1}\|^{2}-\alpha_{k}\|X_{k-1}-p\|^{2}+\Delta
M_{k}+\Delta N_{k}(p,p^{*})$
$\displaystyle+\underbrace{\left(2\alpha_{k}^{2}+(1-\alpha_{k})\left(1-\tfrac{3(1-\alpha_{k})}{2\rho_{k}(1+L\lambda_{k})}\right)\right)}_{\tiny\leq
0}\|X_{k}-X_{k-1}\|^{2}$
$\displaystyle\leq\|X_{k}-p\|^{2}+(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda_{k})}-1\right)\|X_{k}-X_{k-1}\|^{2}-\alpha_{k}\|X_{k-1}-p\|^{2}+\Delta
M_{k}+\Delta N_{k}(p,p^{*}).$

We define $\beta_{k+1}$ as

$\beta_{k+1}\triangleq(1-\alpha_{k})\left(\frac{3}{2\rho_{k}(1+L\lambda_{k})}-1\right)-(1-\alpha_{k+1})\left(\frac{3}{2\rho_{k+1}(1+L\lambda_{k+1})}-1\right),$

and similarly with (4.19), we can show $\\{\beta_{k}\\}$ is non-increasing by
choosing
$\rho_{k}=\tfrac{3(1-\bar{\alpha})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+L\lambda_{k})}$
and $\lambda_{k}\equiv\lambda$. Thus,
$(1-\alpha_{k+1})\left(\tfrac{3}{2\rho_{k+1}(1+L\lambda_{k+1})}-1\right)\leq(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda_{k})}-1\right)$.
Together with $\alpha_{k+1}\geq\alpha_{k}$, the last inequality gives

$\displaystyle\|X_{k+1}-p\|^{2}+(1-\alpha_{k+1})\left(\tfrac{3}{2\rho_{k+1}(1+L\lambda)}-1\right)\|X_{k+1}-X_{k}\|^{2}-\alpha_{k+1}\|X_{k}-p\|^{2}$
$\displaystyle\leq\|X_{k}-p\|^{2}+(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda)}-1\right)\|X_{k}-X_{k-1}\|^{2}-\alpha_{k}\|X_{k-1}-p\|^{2}+\Delta
M_{k}+\Delta N_{k}(p,p^{*}).$

Recall that $\Delta N_{k}(p,p^{\ast})=\Delta
N_{k}(p,0)+2\rho_{k}\lambda\langle p^{\ast},p-Y_{k}\rangle$. Hence, after
setting $\Delta N_{k}(p,0)=\Delta N_{k}(p)$, rearranging the expression given
in the previous display shows that

$\displaystyle 2\rho_{k}\lambda\langle p^{\ast},Y_{k}-p\rangle$
$\displaystyle\leq\left(\|X_{k}-p\|^{2}+(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda)}-1\right)\|X_{k}-X_{k-1}\|^{2}-\alpha_{k}\|X_{k-1}-p\|^{2}\right)$
$\displaystyle-\left(\|X_{k+1}-p\|^{2}+(1-\alpha_{k+1})\left(\tfrac{3}{2\rho_{k+1}(1+L\lambda)}-1\right)\|X_{k+1}-X_{k}\|^{2}-\alpha_{k+1}\|X_{k}-p\|^{2}\right)$
$\displaystyle+\Delta M_{k}+\Delta N_{k}(p).$

Summing over $k=1,\ldots,K$, we obtain

$\displaystyle\sum_{k=1}^{K}2\rho_{k}\lambda\langle p^{\ast},Y_{k}-p\rangle$
$\displaystyle\leq\sum_{k=1}^{K}\left[\left(\|X_{k}-p\|^{2}+(1-\alpha_{k})\left(\tfrac{3}{2\rho_{k}(1+L\lambda)}-1\right)\|X_{k}-X_{k-1}\|^{2}-\alpha_{k}\|X_{k-1}-p\|^{2}\right)\right.$
$\displaystyle\left.-\left(\|X_{k+1}-p\|^{2}+(1-\alpha_{k+1})\left(\tfrac{3}{2\rho_{k+1}(1+L\lambda)}-1\right)\|X_{k+1}-X_{k}\|^{2}-\alpha_{k+1}\|X_{k}-p\|^{2}\right)\right]$
$\displaystyle+\sum_{k=1}^{K}\Delta M_{k}+\sum_{k=1}^{K}\Delta N_{k}(p)$
$\displaystyle\leq\|X_{1}-p\|^{2}+(1-\alpha_{1})\left(\tfrac{3}{2\rho_{1}(1+L\lambda)}-1\right)\|X_{1}-X_{0}\|^{2}-\alpha_{1}\|X_{0}-p\|^{2}$
$\displaystyle+\sum_{k=1}^{K}\Delta M_{k}+\sum_{k=1}^{K}\Delta N_{k}(p)$
$\displaystyle=(1-\alpha_{1})\|X_{1}-p\|^{2}+\sum_{k=1}^{K}\Delta
M_{k}+\sum_{k=1}^{K}\Delta N_{k}(p),$

where we notice $X_{1}=X_{0}$ in the last inequality. $\blacksquare$
$\blacksquare$

Next, we derive a rate statement in terms of the gap function, using the
averaged sequence

$\bar{X}_{K}\triangleq\tfrac{\sum_{k=1}^{K}\rho_{k}Y_{k}}{\sum_{k=1}^{K}\rho_{k}}.$
(4.56)

######  Theorem 4.14 (Rate and oracle complexity under monotonicity of $V$).

Consider the sequence $(X_{k})_{k\in\mathbb{N}}$ generated RISFBF. Suppose
Assumptions 1-5 hold. Suppose $m_{k}\triangleq\lfloor k^{a}\rfloor$ and
$\lambda_{k}=\lambda\in(0,1/(2L))$ for every $k\in\mathbb{N}$ where $a>1$.
Suppose $(\alpha_{k})_{k\in\mathbb{N}}$ is a non-decreasing sequence such that
$0<\alpha_{k}\leq\bar{\alpha}<1$,
$\rho_{k}=\tfrac{3(1-\bar{\alpha})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+L\lambda)}$
for every $k\in\mathbb{N}$. Then the following hold for any $K\in\mathbb{N}$:

  * (i)

$\mathbb{E}[\mathsf{Gap}(\bar{X}_{K}|\mathsf{C})]\leq\mathcal{O}\left(\tfrac{1}{K}\right).$

  * (ii)

Given $\varepsilon>0$, define
$K_{\varepsilon}\triangleq\\{k\in\mathbb{N}|\mathbb{E}[\mathsf{Gap}(\bar{X}_{k}|\mathsf{C})]\leq\varepsilon\\}$,
then
$\sum_{k=1}^{K_{\varepsilon}}m_{k}\leq\mathcal{O}\left(\tfrac{1}{\varepsilon^{1+a}}\right).$

The proof of this Theorem builds on an idea which is frequently used in the
analysis of stochastic approximation algorithms, and can at least be traced
back to the robust stochastic approximation approach of [61]. In order to
bound the expectation of the gap function, we construct an auxiliary process
which allows us to majorize the gap via a quantity which is independent of the
reference points. Once this is achieved, a simple variance bound completes the
result.

###### Proof of Theorem 4.14.

We define an auxiliary process $(\Psi_{k})_{k\in\mathbb{N}}$ such that

$\Psi_{k+1}\triangleq\Psi_{k}+\rho_{k}\lambda_{k}W_{k},\quad\Psi_{1}\in\operatorname{dom}(T).$
(4.57)

Then,

$\displaystyle\lVert\Psi_{k+1}-p\rVert^{2}$
$\displaystyle=\lVert(\Psi_{k}-p)+\rho_{k}\lambda_{k}W_{k}\rVert^{2}=\lVert\Psi_{k}-p\rVert^{2}+\rho^{2}_{k}\lambda_{k}^{2}\lVert
W_{k}\rVert^{2}+2\rho_{k}\lambda_{k}\langle\Psi_{k}-p,W_{k}\rangle,$

so that

$2\rho_{k}\lambda_{k}\langle
W_{k},p-\Psi_{k}\rangle=\lVert\Psi_{k}-p\rVert^{2}-\lVert\Psi_{k+1}-p\rVert^{2}+\rho^{2}_{k}\lambda_{k}^{2}\lVert
W_{k}\rVert^{2}.$

Introducing the iterate $Y_{k}$, the above implies

$\displaystyle 2\rho_{k}\lambda_{k}\langle W_{k},p-Y_{k}\rangle$
$\displaystyle=2\rho_{k}\lambda_{k}\langle
W_{k},p-\Psi_{k}\rangle+2\rho_{k}\lambda_{k}\langle
W_{k},\Psi_{k}-Y_{k}\rangle$
$\displaystyle=\lVert\Psi_{k}-p\rVert^{2}-\lVert\Psi_{k+1}-p\rVert^{2}+\rho^{2}_{k}\lambda_{k}^{2}\lVert
W_{k}\rVert^{2}+2\rho_{k}\lambda_{k}\langle W_{k},\Psi_{k}-Y_{k}\rangle.$

As $\Delta N_{k}(p)=2\rho_{k}\lambda_{k}\langle W_{k},p-Y_{k}\rangle$, this
implies via a telescopian sum argument

$\sum_{k=1}^{K}\Delta
N_{k}(p)\leq\lVert\Psi_{1}-p\rVert^{2}+\sum_{k=1}^{K}\rho^{2}_{k}\lambda^{2}_{k}\lVert
W_{k}\rVert^{2}+\sum_{k=1}^{K}2\rho_{k}\lambda_{k}\langle
W_{k},\Psi_{k}-Y_{k}\rangle.$ (4.58)

Using Lemma 4.13 and setting $\lambda_{k}\equiv\lambda$, for any
$(p,p^{\ast})\in\operatorname{gr}(F)$ it holds true that

$\displaystyle\sum_{k=1}^{K}2\rho_{k}\lambda\langle p^{\ast},Y_{k}-p\rangle$
$\displaystyle\leq(1-\alpha_{1})\lVert X_{1}-p\rVert^{2}+\sum_{k=1}^{K}\Delta
M_{k}+\sum_{k=1}^{K}\Delta N_{k}(p).$

Define $\mathtt{c}_{1}\triangleq(1-\alpha_{1})\lVert X_{1}-p\rVert^{2}$,
divide both sides by $\sum_{k=1}^{K}\rho_{k}$ and using our definition of an
ergodic average (4.56), this gives

$2\lambda\langle
p^{\ast},\bar{X}_{K}-p\rangle\leq\frac{1}{\sum_{k=1}^{K}\rho_{k}}\left\\{\mathtt{c}_{1}+\sum_{k=1}^{K}\Delta
M_{k}+\sum_{k=1}^{K}\Delta N_{k}(p)\right\\}.$

Using the bound established in eq. (4.58), it follows

$\displaystyle 2\lambda\langle
p^{\ast},\bar{X}_{K}-p\rangle\leq\frac{1}{\sum_{k=1}^{K}\rho_{k}}\left\\{\mathtt{c}_{1}+\sum_{k=1}^{K}\Delta
M_{k}+\lVert\Psi_{1}-p\rVert^{2}+\sum_{k=1}^{K}\rho^{2}_{k}\lambda^{2}\lVert
W_{k}\rVert^{2}+\sum_{k=1}^{K}2\rho_{k}\lambda\langle
W_{k},\Psi_{k}-Y_{k}\rangle\right\\}.$

Choosing $\Psi_{1},p\in\mathsf{C}$ and introducing
$\mathtt{c}_{2}\triangleq\mathtt{c}_{1}+4D^{2}$, we see that the above can be
bounded by a random quantity which is independent of $p$:

$\displaystyle 2\lambda\langle
p^{\ast},\bar{X}_{K}-p\rangle\leq\frac{1}{\sum_{k=1}^{K}\rho_{k}}\left\\{\mathtt{c}_{2}+\sum_{k=1}^{K}\Delta
M_{k}+\sum_{k=1}^{K}\rho^{2}_{k}\lambda^{2}\lVert
W_{k}\rVert^{2}+\sum_{k=1}^{K}2\rho_{k}\lambda_{k}\langle
W_{k},\Psi_{k}-Y_{k}\rangle\right\\}.$

Taking the supremum over pairs $(p,p^{\ast})$ such that $p\in\mathcal{C}$ and
$p^{\ast}\in F(y)$, it follows

$2\lambda\mathsf{Gap}(\bar{X}_{K}|\mathsf{C})\leq\frac{\mathtt{c}_{2}}{\sum_{k=1}^{K}\rho_{k}}+\frac{\sum_{k=1}^{K}\Delta
M_{k}+\sum_{k=1}^{K}\rho^{2}_{k}\lambda^{2}\lVert
W_{k}\rVert^{2}+\sum_{k=1}^{K}2\rho_{k}\lambda_{k}\langle
W_{k},\Psi_{k}-Y_{k}\rangle}{\sum_{k=1}^{K}\rho_{k}}$ (4.59)

In order to proceed, we bound the first moment of the process $\Delta M_{k}$
in the same way as in (4.17), in order to get

$\displaystyle\mathbb{E}[\Delta M_{k}|\mathcal{F}_{k}]$
$\displaystyle\leq\frac{6\rho_{k}\lambda_{k}^{2}}{1+L\lambda_{k}}\mathbb{E}[\lVert
W_{k}\rVert^{2}|\mathcal{F}_{k}]+\lambda^{2}_{k}\left(\frac{6\rho_{k}}{1+L\lambda_{k}}+\frac{\rho_{k}\lambda^{2}_{k}}{2}\right)\mathbb{E}[\lVert
U_{k}\rVert^{2}|\mathcal{F}_{k}]$
$\displaystyle=\frac{\left(\frac{12\rho_{k}\lambda^{2}_{k}}{1+L\lambda_{k}}\sigma^{2}+\frac{\rho_{k}\lambda^{2}_{k}}{2}\sigma^{2}\right)}{m_{k}}\triangleq\frac{\mathtt{a}_{k}\sigma^{2}}{m_{k}}.$

Next, we take expectations on both sides of inequality (4.59), and use the
bound (3.8), and $\mathbb{E}[\langle
W_{k},\Psi_{k}-Y_{k}\rangle]=\mathbb{E}\left[\mathbb{E}\left(\langle
W_{k},\Psi_{k}-Y_{k}\rangle|\hat{\mathcal{F}}_{k}\right)\right]=0.$ This
yields

$\displaystyle
2\lambda\mathbb{E}\left[\mathsf{Gap}(\bar{X}_{K}|\mathsf{C})\right]\leq\frac{\mathtt{c}_{2}}{\sum_{k=1}^{K}\rho_{k}}+\frac{1}{\sum_{k=1}^{K}\rho_{k}}\left(\sum_{k=1}^{K}\frac{\mathtt{a}_{k}\sigma^{2}}{m_{k}}+\sum_{k=1}^{K}\rho^{2}_{k}\lambda^{2}\frac{\sigma^{2}}{m_{k}}\right).$

Since $\alpha_{k}\uparrow\bar{\alpha}\in(0,1)$, we know that
$\rho_{k}\geq\tilde{\rho}\triangleq\frac{3(1-\bar{\alpha}^{2})}{2(1+L\lambda)(2\bar{\alpha}^{2}+1)}$.
Similarly, since $2\alpha^{2}_{k}-\alpha_{k}+1\geq 7/8$ for all $k$, it
follows $\rho_{k}\leq\bar{\rho}\triangleq\frac{12(1-\bar{\alpha})^{2}}{7}$.
Using this upper and lower bound on the relaxation sequence, we also see that
$\mathtt{a}_{k}\leq\lambda^{2}\left(\frac{12\bar{\rho}}{1+L\lambda}+\frac{\bar{\rho}}{2}\right)\equiv\bar{\mathtt{a}}$,
so that

$\displaystyle
2\lambda\mathbb{E}\left[\mathsf{Gap}(\bar{X}_{K}|\mathsf{C})\right]\leq\frac{\mathtt{c}_{2}}{\tilde{\rho}K}+\frac{1}{\tilde{\rho}K}\left(\bar{\mathtt{a}}\sigma^{2}+\bar{\rho}^{2}\lambda^{2}\sigma^{2}\right)\sum_{k=1}^{K}\frac{1}{m_{k}}\leq\frac{\mathtt{c}_{3}}{K}$

where
$\mathtt{c}_{3}\triangleq\frac{\mathtt{c}_{2}}{\tilde{\rho}}+\frac{1}{\tilde{\rho}}\left(\bar{\mathtt{a}}\sigma^{2}+\bar{\rho}\lambda^{2}\sigma^{2}\right)\sum_{k=1}^{\infty}\frac{1}{m_{k}}$.
Hence, defining the deterministic stopping time
$K_{\varepsilon}=\\{k\in\mathbb{N}|\mathbb{E}[\mathsf{Gap}(\bar{X}_{k}|\mathsf{C})]\leq\varepsilon\\}$,
we see
$K_{\varepsilon}\geq\frac{\mathtt{c}_{3}}{2\lambda\varepsilon}=\frac{\mathtt{c}_{4}}{\varepsilon}$.

(ii). Suppose $m_{k}=\lfloor k^{a}\rfloor$, for $a>1$. Then the oracle
complexity to compute an $\bar{X}_{K}$ such that
$\mathbb{E}[\mathsf{Gap}(\bar{X}_{k}|\mathsf{C})]\leq\epsilon$ is bounded as

$\displaystyle\sum_{k=1}^{K}m_{k}$
$\displaystyle\leq\sum_{k=1}^{\lceil(\mathtt{c}_{4}/\varepsilon)\rceil}m_{k}\leq\sum_{k=1}^{\lceil(\mathtt{c}_{4}/\varepsilon)\rceil}k^{a}\leq\int_{k=1}^{(\mathtt{c}_{4}/\varepsilon)+1}x^{a}dx\leq\tfrac{((\mathtt{c}_{4}/\varepsilon)+1)^{a+1}}{a+1}\leq\left(\tfrac{\mathtt{c}_{4}}{\varepsilon^{a+1}}\right).$

$\blacksquare$ $\blacksquare$

######  Remark 4.6.

In the prior result, we employ a sampling rate $m_{k}=\lfloor k^{a}\rfloor$
where $a>1$. This achieves the optimal rate of convergence. In contrast, the
authors in [45] employ a sampling rate, loosely given by $m_{k}=\lfloor
k^{1+a}(\ln(k))^{1+b}\rfloor$ where $a>0,b\geq-1$ or $a=0,b>0$. We observe
that when $a>0$ and $b\geq-1$, the mini-batch size grows faster than our
proposed $m_{k}$ while it is comparable in the other case.

##  5 Applications

In this section, we compare the proposed scheme with its SA counterparts on a
class of monotone two-stage stochastic variational inequality problems (Sec.
5.1) and a supervised learning problem (Sec. 5.2) and discuss the resulting
performance.

###  5.1 Two-stage stochastic variational inequality problems

In this section, we describe some preliminary computational results obtained
from the (RISFBF) method when applied to a class of two-stage stochastic
variational inequality problems, recently introduced by Rockafellar and Wets
[70].

Consider an imperfectly competitive market with $N$ firms playing a two-stage
game. In the first stage, the firms decide upon their capacity level
$x_{i}\in[l_{i},u_{i}]$, anticipating the expected revenues to be obtained in
the second stage in which they compete by choosing quantities à la Cournot.
The second-stage market is characterized by uncertainty as the per-unit cost
$h_{i}(\xi_{i})$ is realized on the spot and cannot be anticipated. To compute
an equilibrium in this game, we assume that each player is able to take
stochastic recourse by determining production levels $y_{i}(\xi)$, contingent
on random convex costs and capacity levels $x_{i}$. In order to bring this
into the terminology for our problem, let use define the feasible set for
capacity decisions of firm $i$ as
$\mathcal{X}_{i}\triangleq[l_{i},u_{i}]\subset\mathbb{R}_{+}$. The joint
profile of capacity decisions is denoted by an $N$-tuple
$x=(x_{1},\ldots,x_{N})\in\mathcal{X}\triangleq\prod_{i=1}^{N}\mathcal{X}_{i}=\mathcal{X}$.
The capacity choice of player $i$ is then determined as a solution to the
parametrized problem (Play${}_{i}(x_{-i})$)

$\displaystyle\min_{x_{i}\in\mathcal{X}_{i}}\,c_{i}(x_{i})-\left(p(X)x_{i}-\mathbb{E}_{\xi}[\mathcal{Q}_{i}(x_{i},\xi)]\right),$
(Play${}_{i}(x_{-i})$)

where $c_{i}:\mathcal{X}_{i}\to\mathbb{R}_{+}$ is a $\tilde{L}^{c}_{i}$-smooth
and convex cost function and $p(\cdot)$ denotes the inverse-demand function
defined as $p(X)=d-rX$, $d,r>0$. The function $\mathcal{Q}_{i}(\cdot,\xi)$
denotes the optimal cost function of firm $i$ in scenario $\xi\in\Xi$,
assuming a value $\mathcal{Q}_{i}(x_{i},\xi)$ when the capacity level $x_{i}$
is chosen. The recourse function
$\mathbb{E}_{\xi}[\mathcal{Q}_{i}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cdot},\xi)]$
denotes the expectation of the optimal value of the player $i$’s second stage
problem and is defined as

$\displaystyle\mathcal{Q}_{i}(x_{i},\xi)$
$\displaystyle\triangleq\min\\{h_{i}(\xi)y_{i}(\xi)|y_{i}(\xi)\in[0,x_{i}]\\}$
(Rec${}_{i}(x_{-i})$) $\displaystyle=\max\\{\pi_{i}(\xi)x_{i}|\pi_{i}(\xi)\leq
0,h_{i}(\xi)-\pi_{i}(\xi)\geq 0\\}.$

A Nash equilibrium of this game is given by a tuple
$(x^{\ast}_{1},\cdots,x^{\ast}_{N})$ where $x^{\ast}_{i}\mbox{ solves
}(\mbox{Play}_{i}(x_{-i}^{\ast}))$ for each $i=1,2,\ldots,N$. A simple
computation shows that $Q_{i}(x_{i},\xi)=\min\\{0,h_{i}(\xi)x_{i}\\}$, and
hence it is nonsmooth. In order to obtain a smoothed variant, we introduce
$\mathcal{Q}_{i}^{\epsilon}(\cdot,\xi_{i})$, defined as

$\mathcal{Q}^{\epsilon}_{i}(x_{i},\xi)\triangleq\max\\{x_{i}\pi_{i}(\xi)-\tfrac{\epsilon}{2}(\pi_{i}(\xi))^{2}|\pi_{i}(\xi)\leq
0,\pi_{i}(\xi)\leq h_{i}(\xi)\\},\quad\epsilon>0.$

This is the value function of a quadratric program, requiring the maximization
of an $\epsilon$-strongly concave function. Hence,
$\mathcal{Q}_{i}^{\epsilon}(x_{i},\xi)$ is single-valued and
$\nabla_{x_{i}}\mathcal{Q}^{\epsilon}_{i}(\cdot,\xi)$ is
$\tfrac{1}{\epsilon}$-Lipschitz and $\epsilon$-strongly monotone [71, Prop.
12.60] for all $\xi\in\Xi$. The latter is explicitly given by

$\nabla_{x_{i}}\mathcal{Q}^{\epsilon}_{i}(x_{i},\xi)\triangleq\operatorname*{argmax}\\{x_{i}\pi_{i}(\xi)-\tfrac{\epsilon}{2}(\pi_{i}(\xi))^{2}|\pi_{i}(\xi)\leq
0,\pi_{i}(\xi)\leq h_{i}(\xi)\\}.$

Employing this smoothing strategy in our two-stage noncooperative game yields
the individual decision problem

$\displaystyle(\forall i\in\\{1,\ldots,N\\}):\min_{x_{i}\in\mathcal{X}_{i}}\
c_{i}(x_{i})-p(X)x_{i}+\mathbb{E}_{\xi}[\mathcal{Q}^{\epsilon}_{i}(x_{i},\xi)].$
(Play${}^{\epsilon}_{i}(x_{-i})$)

The necessary and sufficient equilibrium conditions of this
$\epsilon$-smoothed game can be compactly represented by the inclusion problem
(SGEϵ)

$\displaystyle 0\in F^{\epsilon}(x)\triangleq V^{\epsilon}(x)+T(x),\mbox{
where
}V^{\epsilon}(x)=C(x)+R(x)+D^{\epsilon}(x),T(x)=\operatorname{\mathsf{N}}_{\mathcal{X}}(x),$
(SGEϵ)

and $C$, $R$, and $D^{\epsilon}$ are single-valued maps given by

$C(x)\triangleq\left(\begin{array}[]{c}c_{1}^{\prime}(x_{1})\\\ \vdots\\\
c_{N}^{\prime}(x_{N})\end{array}\right),\quad R(x)\triangleq
r(X\text{{1}}+x)-d,\text{ and
}D^{\epsilon}(x)\triangleq\left(\begin{array}[]{c}\mathbb{E}_{\xi}[\nabla_{x_{1}}\mathcal{Q}_{1}^{\epsilon}(x_{1},\xi)]\\\
\vdots\\\
\mathbb{E}_{\xi}[\nabla_{x_{N}}\mathcal{Q}_{N}^{\epsilon}(x_{N},\xi)]\end{array}\right).$

We note that the interchange between the expectation and the gradient operator
can be invoked based on smoothness requirements (cf. [76, Th. 7.47]). The
problem (SGEϵ) aligns perfectly with the structured inclusion (MI), in which
$T$ is a maximal monotone map and $V$ is an expectation-valued maximally
monotone map. In addition, we can quantify the Lipschitz constant of $V$ as
$L_{V}=L_{C}+L_{R}+L_{D}^{\epsilon}$ ,where $L_{C}=\max_{1\leq i\leq
N}\tilde{L}^{c}_{i}$,
$L_{R}=r\lVert\operatorname{Id}+\text{{1}}\text{{1}}^{\top}\rVert_{2}=r(N+1)$
and $L_{D}^{\epsilon}=\tfrac{1}{\epsilon}$. Here, $\operatorname{Id}$ is the
$N\times N$ identity matrix, and 1 is the $N\times 1$ vector consisting only
of ones.

Problem parameters for 2-stage SVI. Our numerics are based on specifying
$N=10$, $r=0.1$, and $d=1$. We consider four problem settings of $L_{V}$
ranging from $10,\cdots,10^{4}$ (See Table 1). For each setting, the problem
parameters are defined as follows.

  1. (i)

Specification of $h_{i}(\xi)$. The cost parameters
$h_{i}(\xi_{i})\triangleq\xi_{i}$ where $\xi_{i}\sim{\mathtt{Uniform}}[-5,0]$
and $i=1,\cdots,N$.

  2. (ii)

Specification of $L_{V},L_{R},$ $L_{D}^{\epsilon}$, $L_{C}$, and
$\hat{b}_{1}$. Since $\lVert\operatorname{Id}+11^{T}\rVert_{2}=11$ when
$N=10$,
$L_{R}=r\lVert\operatorname{Id}+\text{{1}}\text{{1}}^{\top}\rVert=1.1$. Let
$\epsilon$ be defined as $\epsilon=\frac{10}{L_{V}}$ and
$L_{D}^{\epsilon}=\tfrac{1}{\epsilon}=\tfrac{L_{V}}{10}$. It follows that
$L_{C}=L_{V}-L_{R}-L_{D}^{\epsilon}$ and $\hat{b}_{1}=L_{C}$.

  3. (iii)

Specification of $c_{i}(x_{i})$. The cost function $c_{i}$ is defined as
$c_{i}(x_{i})=\tfrac{1}{2}\hat{b}_{i}x_{i}^{2}+a_{i}{x_{i}}$ where
$a_{1},\cdots,a_{N}\sim{\mathtt{Uniform}}[2,3]$ and
$\hat{b}_{2},\cdots,\hat{b}_{N}\sim{\mathtt{Uniform}}[0,\hat{b}_{1}].$
Further, $a\triangleq[a_{1},\dots,a_{N}]^{T}\in\mathbb{R}^{N}$ and
$B\triangleq\operatorname{diag}(\hat{b}_{1},\dots,\hat{b}_{N})$ is a diagonal
matrix with nonnegative elements.

Algorithm specifications. We compare (RISFBF) with a standard stochastic
approximation (SA) scheme and a stochastic forward-backward-forward (SFBF)
scheme. Solution quality is compared by estimating the residual function
$\mathsf{res}(x)=\|x-\Pi_{\mathcal{X}}(x-\lambda V^{\epsilon}(x))\|$. All of
the schemes were implemented in MATLAB on a PC with 16GB RAM and 6-Core Intel
Core i7 processor (2.6GHz).

(i) (SA): Stochastic approximation scheme. The (SA) scheme utilizes update
(SA).

$\displaystyle
X_{k+1}:=\Pi_{\mathcal{X}}\left[X_{k}-\lambda_{k}\widehat{V}^{\epsilon}(X_{k},\xi_{k})\right],$
(SA)

where
$V^{\epsilon}(X_{k})=\mathbb{E}_{\xi}[\widehat{V}^{\epsilon}(X_{k},\xi)]$ and
$\lambda_{k}=\tfrac{1}{\sqrt{k}}$. The operator $\Pi_{\mathcal{X}}[\cdot]$
means the orthogonal projection onto the set $\mathcal{X}$. Note that $x_{0}$
is randomly generated in $[0,1]^{N}$.

(ii) (SFBF): Variance-reduced stochastic modified forward-backward scheme.

$\left\\{\begin{array}[]{l}Y_{k}=\Pi_{\mathcal{X}}[X_{k}-\lambda_{k}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}A_{k}(X_{k})}],\\\
X_{k+1}=Y_{k}-\lambda_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B_{k}(Y_{k})-A_{k}(X_{k})}).\end{array}\right.$
(SFBF)

where
$A_{k}(X_{k})=\frac{1}{m_{k}}\sum_{t=1}^{m_{k}}\widehat{V}^{\epsilon}(X_{k},\xi_{k})$,
$B_{k}(Y_{k})=\frac{1}{m_{k}}\sum_{t=1}^{m_{k}}\widehat{V}^{\epsilon}(Y_{k},\eta_{k})$.
We choose a constant $\lambda_{k}\equiv\lambda=\tfrac{1}{4L_{V}}$. We assume
$m_{k}=\lfloor k^{1.01}\rfloor$ for merely monotone problems and
$m_{k}=\lfloor 1.01^{k}\rfloor$ for strongly monotone problems.

(iii) (RISFBF): Relaxed inertial stochastic forward-backward-forward scheme.
We choose a constant steplength $\lambda_{k}\equiv\lambda=\tfrac{1}{4L_{V}}$.
In merely monotone settings, we utilize an increasing sequence
$\alpha_{k}=\alpha_{0}(1-\tfrac{1}{k+1})$, where $\alpha_{0}=0.1$, the
relaxation parameter sequence $\rho_{k}$ defined as
$\rho_{k}=\tfrac{3(1-\alpha_{0})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+L_{V}\lambda)}$,
and $m_{k}=\lfloor k^{1.01}\rfloor$. In strongly monotone regimes, we choose a
constant inertial parameter $\alpha_{k}\equiv\alpha=0.1$, a constant
relaxation parameter $\rho_{k}\equiv\rho=1$, and $m_{k}=\lfloor
1.01^{k}\rfloor$.

In Fig. 1, we compare the three schemes under maximal monotonicity and strong
monotonicity, respectively and examine their senstivities to inertial and
relaxation parameters. Both sets of plots are based on selecting
$L_{V}=10^{2}$.

Figure 1: Trajectories for (SA), (SFBF), and (RISFBF) (left: monotone, right:
s-monotone) Table 1: Comparison of (RISFBF) with (SA) and (SFBF) under various
Lipschitz constant merely monotone, 20000 evaluations
---
$L_{V}$ | RISFBF | SFBF | SA
error | time | CI | error | time | CI | error | time | CI
1e1 | 2.2e-4 | 2.7 | [2.0e-4,2.5e-4] | 1.6e-3 | 2.6 | [1.3e-3,1.8e-3] | 5.3e-2 | 2.7 | [5.0e-2,5.7e-2]
1e2 | 2.7e-4 | 2.7 | [2.5e-4,3.0e-4] | 1.9e-3 | 2.6 | [1.6e-3,2.1e-3] | 6.1e-2 | 2.7 | [5.8e-2,6.4e-2]
1e3 | 6.9e-4 | 2.7 | [6.7e-3,7.1e-4] | 2.2e-3 | 2.6 | [2.0e-3,2.5e-3] | 7.6e-2 | 2.5 | [7.3e-2,7.9e-2]
1e4 | 2.7e-3 | 2.7 | [2.5e-3,3.0e-3] | 5.9e-3 | 2.6 | [5.4e-3,6.2e-3] | 9.4e-2 | 2.6 | [9.0e-1,9.7e-1]
strongly monotone, 20000 evaluations
$L_{V}$ | RISFBF | SFBF | SA
error | time | CI | error | time | CI | error | time | CI
1e1 | 1.5e-6 | 2.6 | [1.3e-6,1.7e-6] | 1.5e-5 | 2.6 | [1.2e-5,1.7e-5] | 2.9e-2 | 2.5 | [2.7e-2,3.1e-2]
1e2 | 3.7e-6 | 2.6 | [3.5e-6,3.9e-6] | 3.6e-5 | 2.5 | [3.3e-5,3.9e-5] | 4.1e-2 | 2.5 | [3.8e-2,4.4e-2]
1e3 | 4.5e-6 | 2.6 | [4.3e-6,4.7e-6] | 5.6e-5 | 2.5 | [4.2e-6,4.7e-6] | 5.5e-2 | 2.4 | [5.2e-2,5.7e-2]
1e4 | 1.4e-5 | 2.6 | [1.1e-5,1.7e-5] | 7.4e-5 | 2.5 | [7.1e-5,7.7e-5] | 6.0e-2 | 2.5 | [5.7e-2,6.3e-2]

Key insights. Several insights may be drawn from Table 1 and Figure 1.

(a) First, from Table 1, one may conclude that on this class of problems,
(RISFBF) and (SFBF) significantly outperform (SA) schemes, which is less
surprising given that both schemes employ an increasing mini-batch sizes,
leading to performance akin to that seen in deterministic schemes. We should
note that when $\mathcal{X}$ is somewhat more complicated, the difference in
run-times between SA schemes and mini-batch variants becomes for more
pronounced; in this instance, the set $\mathcal{X}$ is relatively simple to
project onto and there is little difference in run-time across the three
schemes.

(b) Second, we observe that while both (SFBF) and (RISFBF) schemes can contend
with poorly conditioned problems, as seen by noting that as $L_{V}$ grows,
their performance does not degenerate significantly in terms of empirical
error; However, in both monotone and strongly monotone regimes, (RISFBF)
provide consistently better solutions in terms of empirical error over (SFBF).
Figure 1 displays the range of trajectories obtained for differing relaxation
and inertial parameters and in the instances considered, (RISFBF) shows
consistent benefits over (SFBF).

(c) Third, since such schemes display geometric rates of convergence for
strongly monotone inclusion problems, this improvement is reflected in terms
of the empirical errors for strongly monotone vs monotone regimes.

###  5.2 Supervised learning with group variable selection

Our second numerical example considers the following population risk
formulation of a composite absolute penalty (CAP) problem arising in
supervised statistical learning [88]

$\min_{w\in\mathcal{W}}\tfrac{1}{2}\mathbb{E}_{(a,b)}[(a^{\top}w-b)^{2}]+\eta\sum_{\textsl{g}\in\mathcal{S}}\|w_{\textsl{g}}\|_{2},$
(CAP)

where the feasible set $\mathcal{W}\subseteq\mathbb{R}^{d}$ is a Euclidean
ball with $\mathcal{W}\triangleq\\{w\in\mathbb{R}^{d}\mid\|w\|_{2}\leq D\\}$,
$\xi=(a,b)\in\mathbb{R}^{d}\times\mathbb{R}$ denotes the random variable
consisting of a set of predictors $a$ and output $b$. The parameter vector $w$
is the sparse linear hypothesis to be learned. The sparsity structure of $w$
is represented by group $\mathcal{S}\in 2^{\\{1,\dots,l\\}}$. When the groups
in $\mathcal{S}$ do not overlap,
$\sum_{\textsl{g}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{S}}}\|w_{\textsl{g}}\|_{2}$
is referred to as the group lasso penalty [33, 46]. When the groups in
$\mathcal{S}$ form a partition of the set of predictors, then
$\sum_{\textsl{g}\in\mathcal{S}}\|w_{\textsl{g}}\|_{2}$ is a norm afflicted by
singularities when some components $w_{\textsl{g}}$ are equal to zero. For any
$\textsl{g}\in\\{1,\cdots,l\\}$, $w_{\textsl{g}}$ is a sparse vector
constructed by components of $x$ whose indices are in g, i.e.,
$w_{\textsl{g}}:=(w_{i})_{i\in\textsl{g}}$ with few non-zero components in
$w_{\textsl{g}}$. Here, we assume that each group $\textsl{g}\in\mathcal{S}$
consists of $k$ elements. Introduce the linear operator
$L:\mathbb{R}^{d}\to\mathbb{R}^{k}\underbrace{\times\cdots\times}_{l-\text{times}}\mathbb{R}^{k}$,
given by $Lw=[\eta w_{g_{1}},\ldots,\eta w_{g_{l}}]$. Let us also define

$\displaystyle
Q=\mathbb{E}_{\xi}[aa^{\top}],q=\mathbb{E}_{\xi}[ab],c=\frac{1}{2}\mathbb{E}_{\xi}[b^{2}],$
$\displaystyle h(w)\triangleq\frac{1}{2}w^{\top}Qw-w^{\top}q+c,\text{ and
}f(w)\triangleq\delta_{\mathcal{W}}(w),$

where $\delta_{\mathcal{W}}(\cdot)$ denotes the indicator function with
respect to the set $\mathcal{W}$. Then (CAP) becomes

$\displaystyle\min_{w\in\mathbb{R}^{d}}\
\\{h(w)+g(Lw)+f(w)\\},\quad\text{where
}g(y_{1},\ldots,y_{l})\triangleq\sum_{i=1}^{l}\lVert y_{i}\rVert.$

This is clearly seen to be a special instance of the convex programming
problem (1.2). Specifically, we let $\mathsf{H}_{1}=\mathbb{R}^{d}$ with the
standard Euclidean norm, and
$\mathsf{H}_{2}=\mathbb{R}^{k}\underbrace{\times\cdots\times}_{l-\text{times}}\mathbb{R}^{k}$
with the product norm

$\lVert(y_{1},\ldots,y_{l})\rVert_{\mathsf{H}_{2}}\triangleq\sum_{i=1}^{l}\lVert
y_{i}\rVert_{2}.$

Since

$g^{\ast}(v_{1},\ldots,v_{l})=\sum_{i=1}^{l}\delta_{\mathbb{B}(0,1)}(v_{i})\qquad\forall
v=(v_{1},\ldots,v_{l})\in\mathbb{R}^{k}\underbrace{\times\cdots\times}_{l-\text{times}}\mathbb{R}^{k},$

the Fenchel-dual takes the form (1.3). Accordingly, a primal-dual pair for
(CAP) is a root of the monotone inclusion (MI) with

$V(w,v)=(\nabla h(w)+L^{\ast}v,-Lw)\text{ and }T(w,v)\triangleq\partial
f(w)\times\partial g^{\ast}(v)$

involving $d+kl$ variables.

Problem parameters for (CAP): We simulated data with $d=82$, covered by 10
groups of 10 variables with 2 variables of overlap between two successive
groups: $\\{1,\dots,10\\},\\{9,\dots,18\\},\dots,\\{73,\dots,82\\}$. We assume
the nonzeros of $w_{\rm true}$ lie in the union of groups 4 and 5 and sampled
from i.i.d. Gaussian variables. The operator $V(w,v)$ is estimated by the
mini-batch estimator using $m_{k}$ iid copies of the random input-output pair
$\xi=(a,b)\in\mathbb{R}^{d}\times\mathbb{R}$. Specifically, we draw each
coordinate of the random vector $a$ from the standard Gaussian distribution
${\mathtt{N}}(0,1)$ and generate $b=a^{\top}w_{\rm true}+\varepsilon$, for
$\varepsilon\sim{\mathtt{N}}(0,\sigma^{2}_{\varepsilon})$. In the concrete
experiment reported here, the error variance is taken as
$\sigma_{\varepsilon}=0.1$. In all instances, the regularization parameter is
chosen as $\eta=10^{-4}$. The accuracy of feature extraction of algorithm
output $w$ is evaluated by the relative error to the ground truth, defined as

$\frac{\|w-w_{\rm true}\|_{2}}{\|w_{\rm true}\|_{2}}.$

Algorithm specifications. We compare (RISFBF) with stochastic extragradient
(SEG) and stochastic forward-backward-forward (SFBF) schemes and specify their
algorithm parameters. Again, all the schemes are run on MATLAB 2018b on a PC
with 16GB RAM and 6-Core Intel Core i7 processor (2.6$\times$8GHz).

(i) (SEG): Variance-reduced stochastic extragradient scheme. Set
$\mathcal{X}\triangleq\mathcal{W}\times\operatorname{dom}(g^{\ast})$. The
(SEG) scheme [45] utilizes update (SEG).

$\displaystyle\begin{aligned}
Y_{k}&:=\Pi_{\mathcal{X}}\left[X_{k}-\lambda_{k}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}A_{k}(X_{k})}\right],\\\
X_{k+1}&:=\Pi_{\mathcal{X}}\left[X_{k}-\lambda_{k}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B_{k}(Y_{k})}\right],\end{aligned}$
(SEG)

where $A_{k}(X_{k})=\frac{1}{m_{k}}\sum_{t=1}^{m_{k}}V(X_{k},\xi_{k})$,
$B_{k}(Y_{k})=\frac{1}{m_{k}}\sum_{t=1}^{m_{k}}V(Y_{k},\eta_{k})$. In this
scheme, $\lambda_{k}\equiv\lambda$ is chosen to be
$\tfrac{1}{4{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}L_{V}}}$
(${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}L_{V}}$
is the Lipschitz constant of $V$). We assume
$m_{k}=\left\lfloor\tfrac{k^{1.1}}{n}\right\rfloor$.

(ii) (SFBF): Variance-reduced stochastic modified forward-backward scheme. We
employ the algorithm paramters employed in (i). Specifically, we choose a
constant
$\lambda_{k}\equiv\lambda=\tfrac{1}{4{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}L_{V}}}$
and $m_{k}=\left\lfloor\tfrac{k^{1.1}}{n}\right\rfloor$.

(iii) (RISFBF): Relaxed inertial stochastic forward-backward-forward scheme.
Here, we employ a constant steplength
$\lambda_{k}\equiv\lambda=\tfrac{1}{4{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}L_{V}}}$,
an increasing sequence $\alpha_{k}=\alpha_{0}(1-\tfrac{1}{k+1})$, where
$\alpha_{0}=0.85$, a relaxation parameter sequence
$\rho_{k}=\tfrac{3(1-\alpha_{0})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}L_{V}}\lambda)}$,
and assume $m_{k}=\left\lfloor\tfrac{k^{1.1}}{n}\right\rfloor$.

Table 2: The comparison of the RISFBF, SFBF and SEG algorithms in solving
(CAP)

Iteration |  | RISFBF |  |  | SFBF |  |  | SEG |
---|---|---|---|---|---|---|---|---|---
$N$ |  | Rel. error | CPU |  | Rel. error | CPU |  | Rel. error | CPU
400 |  | 5.4e-1 | 0.1 |  | 34.6 | 0.1 |  | 34.7 | 0.1
800 |  | 8.1e-3 | 0.5 |  | 1.1e-1 | 0.5 |  | 1.5e-1 | 0.5
1200 |  | 6.0e-3 | 1.1 |  | 2.4e-2 | 1.1 |  | 2.4e-2 | 1.1
1600 |  | 5.2e-3 | 2.0 |  | 2.0e-2 | 2.0 |  | 1.9e-2 | 2.0
2000 |  | 4.6e-3 | 3.1 |  | 1.6e-2 | 3.1 |  | 1.5e-2 | 3.1

  * •

The relative error and CPU time in the table is the average results of 20 runs

Insights. We compare the performance of the schemes in Table 2 and observe
that (RISFBF) outperforms its competitors others in extracting the underlying
feature of the datasets. In Fig. 2, trajectories for (RISFBF), (SFBF) and
(SEG) are presented where a consistent benefit of employing (RISFBF) can be
seen for a range of choices of $\alpha_{0}$.

Figure 2: Trajectories for (SEG), (SFBF), and (RISFBF) for problem (CAP).

##  6 Conclusion

In a general structured monotone inclusion setting in Hilbert spaces, we
introduce a relaxed inertial stochastic algorithm based on Tseng’s forward-
backward-forward splitting method. Motivated by the gaps in convergence claims
and rate statements in both deterministic and stochastic regimes, we develop a
variance-reduced framework and make the following contributions: (i)
Asymptotic convergence guarantees are provided under both increasing and
constant mini-batch sizes, the latter requiring somewhat stronger assumptions
on $V$; (ii) When $V$ is monotone, rate statements provided in terms of a
restricted gap function, inspired by the Fitzpatrick function for inclusions,
show that the expected gap of an averaged sequence diminishes at the rate of
$\mathcal{O}(1/k)$ and oracle complexity of computing an $\epsilon$-solution
is $\mathcal{O}(1/\epsilon^{1+a})$ where $a>1$; (iii) When $V$ is strongly
monotone, a non-asymptotic linear rate statement can be proven with an oracle
complexity of $\mathcal{O}(\log(1/\epsilon))$ of computing an
$\epsilon$-solution. In addition, a perturbed linear rate is also developed.
It is worth emphasizing that the rate statements in the strongly monotone
regime accommodate the possibility of a biased SO. Unfortunately, the growth
rates in batch-size may be onerous in some situations, motivating the analysis
of a polynomial growth rate in sample-size which is easily modulated. This
leads to an associated polynomial rate of convergence.

Various open questions arise from our analysis. First, we exclusively focused
on a variance reduction technique based on increasing mini-batches. From the
point of view of computations and oracle complexity, this approach can become
quite costly. Exploiting different variance reduction techniques, taking
perhaps special structure of the single-valued operator $V$ into account (as
in [66]), has the potential of improving the computational complexity of our
proposed method. At the same time, this will complicate the analysis of the
variance of the stochastic estimators considerably and consequently, we leave
this as an important question for future research.

Second, our analysis needs knowledge about the Lipschitz constant $L$. While
in deterministic regimes, line search techniques have obviated such a need,
such avenues are far more challenging to adopt in stochastic regimes. Efforts
to address this in variational regimes have centered around leveraging
empirical process theory [44]. This remains a goal of future research. Another
avenue emerges in applications where we can gain a reasonably good estimate
about this quantity via some pre-processing of the data (see e.g. Section 6 in
[39]). Developing such an adaptive framework robust to noise is an important
topic for future research.

## Acknowledgments

The authors thank Radu I. Bot̨ and Robert E. Csetnek for valuable discussions
on this topic. Special thanks to Patrick Johnstone for helping us to clarify
the discussion in Section 4.2.

##  Appendix A Technical Facts

######  Lemma A.1.

For $x,y\in\mathsf{H}$ and scalars $\alpha,\beta\geq 0$ with $\alpha+\beta=1$,
it holds that

$\lVert\alpha x+\beta y\rVert^{2}=\alpha\lVert x\rVert^{2}+\beta\lVert
y\rVert^{2}-\alpha\beta\lVert x-y\rVert^{2}.$ (A.1)

We recall the Minkowski inequality: For $X,Y\in
L^{p}(\Omega,\mathcal{F},\mathbb{P};\mathsf{H}),\mathcal{G}\subseteq\mathcal{F}$
and $p\in[1,\infty]$,

$\mathbb{E}[\lVert X+Y\rVert^{p}|\mathcal{G}]^{1/p}\leq\mathbb{E}[\lVert
X\rVert^{p}|\mathcal{G}]^{1/p}+\mathbb{E}[\lVert
Y\rVert^{p}|\mathcal{G}]^{1/p}.$ (A.2)

In the convergence analysis, we use the Robbins-Siegmund Lemma [68, Lemma 11,
pg. 50].

######  Lemma A.2 (Robbins-Siegmund).

Let $(\Omega,\mathcal{F},{\mathbb{F}}=(\mathcal{F}_{n})_{n\geq 0},\mathbb{P})$
be a discrete stochastic basis. Let $(v_{n})_{n\geq 1},(u_{n})_{n\geq
1}\in\ell^{0}_{+}({\mathbb{F}})$ and $(\theta_{n})_{n\geq
1},(\beta_{n})_{n\geq 1}\in\ell^{1}_{+}({\mathbb{F}})$ be such that for all
$n\geq 0$,

$\mathbb{E}[v_{n+1}|\mathcal{F}_{n}]\leq(1+\theta_{n})v_{n}-u_{n}+\beta_{n}\qquad\mathbb{P}-\text{a.s.
}.$

Then $(v_{n})_{n\geq 0}$ converges a.s. to a random variable $v$, and
$(u_{n})_{n\geq 1}\in\ell^{1}_{+}({\mathbb{F}})$.

The next technical lemma will be needed in deriving a linear convergence
result.

######  Lemma A.3.

Let $z\geq 0$ and $0<q<p<1$. Then, if $D\geq\frac{1}{\exp(1)\ln(p/q)}$, it
holds true that $zq^{z}\leq Dp^{z}$ for all $z\geq 0$.

###### Proof.

We want to find a positive constant $D_{\min}>0$ such that
$D_{\min}\exp(z\ln(p))=z\exp(z\ln(q))$ for all $z>0$. Choosing $D$ larger than
this, gives a valid value. Rearranging, this is equivalent to
$D=z\left(\frac{q}{p}\right)^{z}\geq 0$ for all $z\geq 0$, or, which is still
equivalent to $\ln(D)-\ln(z)-z\ln(q/p)=0.$ Define the extended-valued function
$f:[0,\infty)\to[-\infty,\infty]$ by $f(z)=\ln(D)-\ln(z)-\ln(q/p)$ if $z>0$,
and $f(z)=\infty$ if $z=0$. Then, for all $z>0$, simple calculus show
$f^{\prime}(z)=-1/z-\ln(q/p)$ and $f^{\prime\prime}(z)=1/z^{2}$. Hence,
$z\mapsto f(z)$ is a convex function with a unique minimum
$z_{\min}=\frac{1}{\ln(p/q)}>0$ and a corresponding function value
$f(z_{\min})=\ln(D)+\ln(\ln(p/q))+1$. Hence, for $D\geq
D_{\min}=\frac{1}{\exp(1)\ln(p/q)}$, we see that $f(z_{\min})>0$, and thus
$zq^{z}\leq Dp^{z}$ for all $z\geq 0$. $\blacksquare$ $\blacksquare$

## References

  * Attouch and Cabot [2019] Hedy Attouch and Alexandre Cabot.  Convergence of a relaxed inertial forward–backward algorithm for structured monotone inclusions.  _Applied Mathematics & Optimization_, 80(3):547–598, 2019.  doi: 10.1007/s00245-019-09584-z.  URL https://doi.org/10.1007/s00245-019-09584-z.
  * Attouch and Cabot [2020] Hedy Attouch and Alexandre Cabot.  Convergence of a relaxed inertial proximal algorithm for maximally monotone operators.  _Mathematical Programming_ , 184(1):243–287, 2020\.  doi: 10.1007/s10107-019-01412-0.  URL https://doi.org/10.1007/s10107-019-01412-0.
  * Attouch and Maingé [2011] Hedy Attouch and Paul-Emile Maingé.  Asymptotic behavior of second-order dissipative evolution equations combining potential with non-potential effects.  _ESAIM: Control, Optimisation and Calculus of Variations_ , 17(3):836–857, 2011.
  * Attouch and Peypouquet [2019] Hedy Attouch and Juan Peypouquet.  Convergence of inertial dynamics and proximal algorithms governed by maximally monotone operators.  _Mathematical Programming_ , 174(1):391–432, 2019\.  ISSN 1436-4646.  doi: 10.1007/s10107-018-1252-x.  URL https://doi.org/10.1007/s10107-018-1252-x.
  * Attouch et al. [2010] Hédy Attouch, Luis M. Briceno-Arias, and Patrick L. Combettes.  A parallel splitting method for coupled monotone inclusions.  _SIAM Journal on Control and Optimization_ , 48(5):3246–3270, 2021/03/28 2010.  doi: 10.1137/090754297.  URL https://doi.org/10.1137/090754297.
  * Auslender et al. [1974] A. Auslender, M. Gourgand, and A. Guillet.  Resolution numerique d’inegalites variationnelles.  In _Lecture Notes in Economics and Mathematical Systems (Mathematical Economics)_ , 1974.
  * Barty et al. [2007] Kengy Barty, Jean-Sébastien Roy, and Cyrille Strugarek.  Hilbert-valued perturbed subgradient algorithms.  _Mathematics of Operations Research_ , 32(3):551–562, 2020/11/23 2007.  doi: 10.1287/moor.1070.0253.  URL https://doi.org/10.1287/moor.1070.0253.
  * Barty et al. [2009] Kengy Barty, Jean-Sébastien Roy, and Cyrille Strugarek.  A stochastic gradient type algorithm for closed-loop problems.  _Mathematical Programming_ , 119(1):51–78, 2009\.  doi: 10.1007/s10107-007-0201-x.  URL https://doi.org/10.1007/s10107-007-0201-x.
  * Bauschke and Combettes [2016] Heinz H. Bauschke and Patrick L. Combettes.  _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_.  Springer - CMS Books in Mathematics, 2016.
  * Börgens and Kanzow [2021] Eike Börgens and Christian Kanzow.  Admm-type methods for generalized nash equilibrium problems in hilbert spaces.  _SIAM J. Optim._ , pages 377–403, January 2021.  ISSN 1052-6234.  doi: 10.1137/19M1284336.  URL https://doi.org/10.1137/19M1284336.
  * Borwein and Dutta [2016] Jonathan M Borwein and Joydeep Dutta.  Maximal monotone inclusions and fitzpatrick functions.  _Journal of Optimization Theory and Applications_ , 171(3):757–784, 2016.
  * Bot and Csetnek [2016] Radu Ioan Bot and Erno Robert Csetnek.  An inertial forward-backward-forward primal-dual splitting algorithm for solving monotone inclusion problems.  _Numerical Algorithms_ , 71(3):519–540, 2016\.  doi: 10.1007/s11075-015-0007-5.  URL https://doi.org/10.1007/s11075-015-0007-5.
  * Boţ and Csetnek [2016] Radu Ioan Boţ and ErnöRobert Csetnek.  Second order forward-backward dynamical systems for monotone inclusion problems.  _SIAM Journal on Control and Optimization_ , 54(3):1423–1443, 2021/06/24 2016.  doi: 10.1137/15M1012657.  URL https://doi.org/10.1137/15M1012657.
  * Bot et al. [2020] Radu Ioan Bot, Michael Sedlmayer, and Phan Tu Vuong.  A relaxed inertial forward-backward-forward algorithm for solving monotone inclusions with application to gans.  _arXiv preprint arXiv:2003.07886_ , 2020.
  * Boţ et al. [2021] Radu Ioan Boţ, Panayotis Mertikopoulos, Mathias Staudigl, and Phan Tu Vuong.  Minibatch forward-backward-forward methods for solving stochastic variational inequalities.  _Stochastic Systems_ , 2021.  doi: 10.1287/stsy.2019.0064.  URL https://doi.org/10.1287/stsy.2019.0064.
  * Briceño-Arias and Combettes [2011] Luis M. Briceño-Arias and Patrick L. Combettes.  A monotone+skew splitting model for composite monotone inclusions in duality.  _SIAM Journal on Optimization_ , 21(4):1230–1250, 2011.  doi: 10.1137/10081602X.  URL https://doi.org/10.1137/10081602X.
  * Briceno-Arias and Combettes [2013] Luis M Briceno-Arias and Patrick L Combettes.  _Monotone operator methods for Nash equilibria in non-potential games_ , pages 143–159.  Springer, 2013.
  * Burachik and Millán [2020] Regina S. Burachik and R. Díaz Millán.  A projection algorithm for non-monotone variational inequalities.  _Set-Valued and Variational Analysis_ , 28(1):149–166, 2020.  doi: 10.1007/s11228-019-00517-0.  URL https://doi.org/10.1007/s11228-019-00517-0.
  * Chen et al. [2014] Yunmei Chen, Guanghui Lan, and Yuyuan Ouyang.  Optimal primal-dual methods for a class of saddle point problems.  _SIAM Journal on Optimization_ , 24(4):1779–1814, 2021/04/27 2014.  doi: 10.1137/130919362.  URL https://doi.org/10.1137/130919362.
  * Chen et al. [2017] Yunmei Chen, Guanghui Lan, and Yuyuan Ouyang.  Accelerated schemes for a class of variational inequalities.  _Mathematical Programming_ , 2017.  doi: 10.1007/s10107-017-1161-4.  URL https://doi.org/10.1007/s10107-017-1161-4.
  * Combettes and Pesquet [2012] Patrick L Combettes and Jean-Christophe Pesquet.  Primal-dual splitting algorithm for solving inclusions with mixtures of composite, lipschitzian, and parallel-sum type monotone operators.  _Set-Valued and variational analysis_ , 20(2):307–330, 2012.
  * Combettes and Pesquet [2015] Patrick L. Combettes and Jean-Christophe Pesquet.  Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping.  _SIAM Journal on Optimization_ , 25(2):1221–1248, 2020/07/11 2015.  doi: 10.1137/140971233.  URL https://doi.org/10.1137/140971233.
  * Combettes and Pesquet [2019] Patrick L. Combettes and Jean-Christophe Pesquet.  Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping ii: mean-square and linear convergence.  _Mathematical Programming_ , 174(1):433–451, 2019\.  doi: 10.1007/s10107-018-1296-y.  URL https://doi.org/10.1007/s10107-018-1296-y.
  * Cui and Shanbhag [2021a] Shisheng Cui and Uday V Shanbhag.  On the computation of equilibria in monotone and potential stochastic hierarchical game.  _arXiv preprint arXiv:2104.07860_ , 2021a.
  * Cui and Shanbhag [2021b] Shisheng Cui and Uday V Shanbhag.  On the analysis of variance-reduced and randomized projection variants of single projection schemes for monotone stochastic variational inequality problems.  _Set-Valued and Variational Analysis (to appear)_ , 2021b.
  * Davis and Drusvyatskiy [2019] Damek Davis and Dmitriy Drusvyatskiy.  Stochastic model-based minimization of weakly convex functions.  _SIAM Journal on Optimization_ , 29(1):207–239, 2019.
  * Diakonikolas et al. [2021] Jelena Diakonikolas, Constantinos Daskalakis, and Michael Jordan.  Efficient methods for structured nonconvex-nonconcave min-max optimization.  _International Conference on Artificial Intelligence and Statistics_ , pages 2746–2754, 2021.
  * Duvocelle et al. [2018] Benoit Duvocelle, Panayotis Mertikopoulos, Mathias Staudigl, and Dries Vermeulen.  Learning in time-varying games.  _arXiv preprint arXiv:1809.03066_ , 2018.
  * Facchinei and Pang [2003] Francisco Facchinei and Jong-shi Pang.  _Finite-Dimensional Variational Inequalities and Complementarity Problems - Volume I and Volume II_.  Springer Series in Operations Research, 2003.
  * Fitzpatrick [1988] Simon Fitzpatrick.  Representing monotone operators by convex functions.  In _Workshop/Miniconference on Functional Analysis and Optimization_ , pages 59–65. Centre for Mathematics and its Applications, Mathematical Sciences Institute …, 1988.
  * Franci et al. [2020] B. Franci, M. Staudigl, and S. Grammatico.  Distributed forward-backward (half) forward algorithms for generalized nash equilibrium seeking.  In _2020 European Control Conference (ECC)_ , pages 1274–1279, 2020\.  doi: 10.23919/ECC51009.2020.9143676.
  * Friedlander and Schmidt [2012] M. P. Friedlander and M. Schmidt.  Hybrid deterministic-stochastic methods for data fitting.  _SIAM J. Scientific Computing_ , 34(3):A1380–A1405, 2012.
  * Friedman et al. [2001] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  _The elements of statistical learning_ , volume 1.  Springer series in statistics Springer, Berlin, 2001.
  * Friesz et al. [1993] Terry L. Friesz, David Bernstein, Tony E. Smith, Roger L. Tobin, and B. W. Wie.  Variational inequality formulation of the dynamic network user equilibrium.  _Operations Research_ , 41(1):179–191, 1993.
  * Fukushima [1996] Masao Fukushima.  The primal douglas-rachford splitting algorithm for a class of monotone mappings with application to the traffic equilibrium problem.  _Mathematical Programming_ , 72(1):1–15, 1996\.  ISSN 1436-4646.  doi: 10.1007/BF02592328.  URL https://doi.org/10.1007/BF02592328.
  * Gadat et al. [2018] Sébastien Gadat, Fabien Panloup, and Sofiane Saadane.  Stochastic heavy ball.  _Electronic Journal of Statistics_ , 12(1):461 – 529, 2018.  doi: 10.1214/18-EJS1395.  URL https://doi.org/10.1214/18-EJS1395.
  * Geiersbach and Pflug [2019] Caroline Geiersbach and Georg Ch. Pflug.  Projected stochastic gradients for convex constrained problems in hilbert spaces.  _SIAM Journal on Optimization_ , 29(3):2079–2099, 2020/11/23 2019.  doi: 10.1137/18M1200208.  URL https://doi.org/10.1137/18M1200208.
  * Geiersbach and Wollner [2020] Caroline Geiersbach and Winnifried Wollner.  A stochastic gradient method with mesh refinement for pde-constrained optimization under uncertainty.  _SIAM Journal on Scientific Computing_ , 42(5):A2750–A2772, 2021/07/16 2020.  doi: 10.1137/19M1263297.  URL https://doi.org/10.1137/19M1263297.
  * Ghadimi et al. [2016] Saeed Ghadimi, Guanghui Lan, and Hongchao Zhang.  Mini-batch stochastic approximation methods for nonconvex stochastic composite optimization.  _Mathematical Programming_ , 155(1-2):267–305, 2016.
  * Gidel et al. [2018] Gauthier Gidel, Hugo Berard, Gaëtan Vignoud, Pascal Vincent, and Simon Lacoste-Julien.  A variational inequality perspective on generative adversarial networks.  _arXiv preprint arXiv:1802.10551_ , 2018.
  * Gower et al. [2020] R. M. Gower, M. Schmidt, F. Bach, and P. Richtárik.  Variance-reduced methods for machine learning.  _Proceedings of the IEEE_ , 108(11):1968–1983, 2020.  doi: 10.1109/JPROC.2020.3028013.
  * Gunzburger et al. [2014] Max D Gunzburger, Clayton G Webster, and Guannan Zhang.  Stochastic finite element methods for partial differential equations with random input data.  _Acta Numer._ , 23:521–650, 2014.
  * Han et al. [2019] Ke Han, Gabriel Eve, and Terry L. Friesz.  Computing dynamic user equilibria on large-scale networks with software implementation.  _Networks and Spatial Economics_ , 19(3):869–902, 2019.  doi: 10.1007/s11067-018-9433-y.  URL https://doi.org/10.1007/s11067-018-9433-y.
  * Iusem et al. [2019] Alfredo N. Iusem, Alejandro Jofré, Roberto I. Oliveira, and Philip Thompson.  Variance-based extragradient methods with line search for stochastic variational inequalities.  _SIAM Journal on Optimization_ , 29(1):175–206, 2021/07/16 2019.  doi: 10.1137/17M1144799.  URL https://doi.org/10.1137/17M1144799.
  * Iusem et al. [2017] AN Iusem, Alejandro Jofré, Roberto I Oliveira, and Philip Thompson.  Extragradient method with variance reduction for stochastic variational inequalities.  _SIAM Journal on Optimization_ , 27(2):686–724
  * Jacob et al. [2009] Laurent Jacob, Guillaume Obozinski, and Jean-Philippe Vert.  Group lasso with overlap and graph lasso.  _Proceedings of the 26th annual international conference on machine learning_ , pages 433–440, 2009.
  * Jalilzadeh et al. [2018] Afrooz Jalilzadeh, Uday V Shanbhag, Jose H Blanchet, and Peter W Glynn.  Smoothed variable sample-size accelerated proximal methods for nonsmooth stochastic convex programs.  _arXiv preprint arXiv:1803.00718_ , 2018.
  * Jiang and Xu [2008] H. Jiang and H. Xu.  Stochastic approximation approaches to the stochastic variational inequality problem.  _IEEE Transactions on Automatic Control_ , 53(6):1462–1475, 2008.  ISSN 0018-9286.  doi: 10.1109/TAC.2008.925853.
  * Jofré and Thompson [2019] Alejandro Jofré and Philip Thompson.  On variance reduction for stochastic smooth convex optimization with multiplicative noise.  _Mathematical Programming_ , 174(1-2):253–292, 2019.
  * Juditsky et al. [2011] Anatoli Juditsky, Arkadi Nemirovski, and Claire Tauvel.  Solving variational inequalities with stochastic mirror-prox algorithm.  pages 17–58, 2011.  doi: 10.1214/10-SSY011.  URL http://projecteuclid.org/euclid.ssy/1393252123.
  * Kannan and Shanbhag [2019] Aswin Kannan and Uday V Shanbhag.  Optimal stochastic extragradient schemes for pseudomonotone stochastic variational inequality problems and their variants.  _Computational Optimization and Applications_ , 74(3):779–820, 2019.
  * Lan [2020] Guanghui Lan.  _First-order and Stochastic Optimization Methods for Machine Learning_.  Springer Series in the Data Sciences. Springer Nature, 2020.  ISBN 3030395685.
  * Latafat et al. [2019] Puya Latafat, Nikolaos M. Freris, and Panagiotis Patrinos.  A new randomized block-coordinate primal-dual proximal algorithm for distributed optimization.  _IEEE Transactions on Automatic Control_ , 64(10):4050–4065, 2019.  doi: 10.1109/TAC.2019.2906924.
  * Lei and Shanbhag [2019] Jinlong Lei and Uday V. Shanbhag.  Distributed variable sample-size gradient-response and best-response schemes for stochastic nash equilibrium problems over graphs.  _https://arxiv.org/abs/1811.11246_ , 2019.
  * Lei and Shanbhag [2020] Jinlong Lei and Uday V Shanbhag.  Asynchronous variance-reduced block schemes for composite non-convex stochastic optimization: block-specific steplengths and adapted batch-sizes.  _Optimization Methods and Software_ , pages 1–31, 2020.
  * Lorenz and Pock [2015] Dirk A. Lorenz and Thomas Pock.  An inertial forward-backward algorithm for monotone inclusions.  _Journal of Mathematical Imaging and Vision_ , 51(2):311–325, 2015.  doi: 10.1007/s10851-014-0523-2.  URL https://doi.org/10.1007/s10851-014-0523-2.
  * Malitsky [2020] Yura Malitsky.  Golden ratio algorithms for variational inequalities.  _Mathematical Programming_ , 184(1):383–410, 2020\.  doi: 10.1007/s10107-019-01416-w.  URL https://doi.org/10.1007/s10107-019-01416-w.
  * Mertikopoulos and Staudigl [2017] Panayotis Mertikopoulos and Mathias Staudigl.  Convergence to nash equilibrium in continuous games with noisy first-order feedback.  In _56th IEEE Conference on Decision and Control_ , 2017.
  * Mishchenko et al. [2020] Konstantin Mishchenko, Dmitry Kovalev, Egor Shulgin, Peter Richtárik, and Yura Malitsky.  Revisiting stochastic extragradient.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4573–4582. PMLR, 2020.
  * Monteiro and Svaiter [2010] Renato D. C. Monteiro and B. F. Svaiter.  On the complexity of the hybrid proximal extragradient method for the iterates and the ergodic mean.  _SIAM Journal on Optimization_ , 20(6):2755–2787, 2020/11/24 2010.  doi: 10.1137/090753127.  URL https://doi.org/10.1137/090753127.
  * Nemirovski et al. [2009] Arkadi Nemirovski, Anatoli Juditsky, Guanghui Lan, and Alexander Shapiro.  Robust stochastic approximation approach to stochastic programming.  _SIAM Journal on optimization_ , 19(4):1574–1609, 2009.
  * Nesterov [2004a] Y. Nesterov.  _Introductory Lectures on Convex Optimization: A Basic Course_.  Kluwer, Dordrecht, 2004a.
  * Nesterov [1983] Yu. Nesterov.  A method of solving a convex programming problem with convergence rate $o(1/k^{2})$.  _Soviet Mathematics Doklady_ , 27(2):372–376, 1983.
  * Nesterov [2004b] Yurii Nesterov.  _Introductory Lectures on Convex Optimization: A Basic Course_.  Number 87 in Applied Optimization. Kluwer Academic Publishers, 2004b.
  * Nesterov [2007] Yurii Nesterov.  Dual extrapolation and its applications to solving variational inequalities and related problems.  _Mathematical Programming_ , 109(2):319–344, 2007\.  doi: 10.1007/s10107-006-0034-z.  URL https://doi.org/10.1007/s10107-006-0034-z.
  * Palaniappan and Bach [2016] Balamurugan Palaniappan and Francis Bach.  Stochastic variance reduction methods for saddle-point problems.  In _Advances in Neural Information Processing Systems_ , pages 1416–1424, 2016.
  * Polyak [1964] B. T. Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4(5):1–17, 1964.  doi: https://doi.org/10.1016/0041-5553(64)90137-5.  URL http://www.sciencedirect.com/science/article/pii/0041555364901375.
  * Polyak [1987] Boris T. Polyak.  _Introduction to Optimization_.  Optimization Software, 1987.
  * Rockafellar [1974] R. Tyrrell Rockafellar.  _Conjugate Duality and Optimization_.  Society for Industrial and Applied Mathematics, 1974.  doi: 10.1137/1.9781611970524.  URL https://epubs.siam.org/doi/abs/10.1137/1.9781611970524.
  * Rockafellar and Wets [2017] R. Tyrrell Rockafellar and Roger J.-B. Wets.  Stochastic variational inequalities: single-stage to multistage.  _Math. Program._ , 165(1):331–360, 2017.  doi: 10.1007/s10107-016-0995-5.  URL https://doi.org/10.1007/s10107-016-0995-5.
  * Rockafellar and Wets [1998] Tyrrel R. Rockafellar and Roger J-B. Wets.  _Variational Analysis_.  Springer Verlarg - Grundlehren der Mathematischen Wissenschaften, Volume 317, 1998.
  * Rosasco et al. [2016a] Lorenzo Rosasco, Silvia Villa, and Băng Công Vũ.  A stochastic inertial forward–backward splitting algorithm for multivariate monotone inclusions.  _Optimization_ , 65(6):1293–1314, 2016a.  doi: 10.1080/02331934.2015.1127371.  URL https://doi.org/10.1080/02331934.2015.1127371.
  * Rosasco et al. [2016b] Lorenzo Rosasco, Silvia Villa, and Bang Công Vũ.  Stochastic forward–backward splitting for monotone inclusions.  _Journal of Optimization Theory and Applications_ , 169(2):388–406, 2016b.  doi: 10.1007/s10957-016-0893-2.  URL https://doi.org/10.1007/s10957-016-0893-2.
  * Rudin et al. [1992] Leonid I. Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 60(1):259–268, 1992.  ISSN 0167-2789.  doi: https://doi.org/10.1016/0167-2789(92)90242-F.  URL https://www.sciencedirect.com/science/article/pii/016727899290242F.
  * Shanbhag [2013] Uday V. Shanbhag.  _Stochastic Variational Inequality Problems: Applications, Analysis, and Algorithms_ , chapter Chapter 5, pages 71–107.  2013\.  doi: 10.1287/educ.2013.0120.  URL https://pubsonline.informs.org/doi/abs/10.1287/educ.2013.0120.
  * Shapiro et al. [2009] Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczyński.  _Lectures on stochastic programming: modeling and theory_.  SIAM, 2009.
  * Simons and Zalinescu [2004] S. Simons and C. Zalinescu.  A new proof for rockafellar’s characterization of maximal monotone operators.  132(10):2969–2972, 2004.  URL http://www.jstor.org/stable/4097259.
  * Spall [1992] James C Spall.  Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  _IEEE transactions on automatic control_ , 37(3):332–341, 1992.
  * Spall [1997] James C. Spall.  A one-measurement form of simultaneous perturbation stochastic approximation.  _Automatica_ , 33(1):109–112, 1997.
  * Staudigl and Mertikopoulos [2019] Mathias Staudigl and Panayotis Mertikopoulos.  Convergent noisy forward-backward-forward algorithms in non-monotone variational inequalities.  _IFAC-PapersOnLine_ , 52(3):120–125, 2019.
  * Su et al. [2016] Weijie Su, Stephen Boyd, and Emmanuel J. Candes.  A differential equation for modeling nesterov’s accelerated gradient method: Theory and insights.  _Journal of Machine Learning Research_ , 2016.
  * Thong et al. [2021] Duong Viet Thong, Aviv Gibali, Mathias Staudigl, and Phan Tu Vuong.  Computing dynamic user equilibrium on large-scale networks without knowing global parameters.  _Networks and Spatial Economics (forthcoming) arXiv preprint arXiv:2010.04597_ , 2021.
  * Tibshirani et al. [2005] Robert Tibshirani, Michael Saunders, Saharon Rosset, Ji Zhu, and Keith Knight.  Sparsity and smoothness via the fused lasso.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 67(1):91–108, 2019/04/05 2005.  doi: 10.1111/j.1467-9868.2005.00490.x.  URL https://doi.org/10.1111/j.1467-9868.2005.00490.x.
  * Tibshirani and Taylor [2011] Ryan J. Tibshirani and Jonathan Taylor.  The solution path of the generalized lasso.  _Ann. Statist._ , 39(3):1335–1371, 2011.  doi: 10.1214/11-AOS878.  URL https://projecteuclid.org:443/euclid.aos/1304514656.
  * Tseng [2000] P. Tseng.  A modified forward-backward splitting method for maximal monotone mappings.  _SIAM Journal on Control and Optimization_ , 38(2):431–446, 2018/09/13 2000.  doi: 10.1137/S0363012998338806.  URL https://doi.org/10.1137/S0363012998338806.
  * Yi and Pavel [2019] Peng Yi and Lacra Pavel.  An operator splitting approach for distributed generalized nash equilibria computation.  _Automatica_ , 102:111–121, 2019.  doi: https://doi.org/10.1016/j.automatica.2019.01.008.  URL http://www.sciencedirect.com/science/article/pii/S0005109819300081.
  * Yousefian et al. [2017] Farzad Yousefian, Angelia Nedić, and Uday V. Shanbhag.  On smoothing, regularization, and averaging in stochastic approximation methods for stochastic variational inequality problems.  _Mathematical Programming_ , 165(1):391–431, 2017\.  doi: 10.1007/s10107-017-1175-y.  URL https://doi.org/10.1007/s10107-017-1175-y.
  * Zhao et al. [2009] Peng Zhao, Guilherme Rocha, and Bin Yu.  The composite absolute penalties family for grouped and hierarchical variable selection.  _The Annals of Statistics_ , 37(6A):3468–3497, 12 2009.  doi: 10.1214/07-AOS584.  URL https://doi.org/10.1214/07-AOS584.

  *[SFBF]: stochastic forward-backward-forward
  *[SO]: stochastic oracle
  *[SVI]: stochastic variational inequality
  *[RISFBF]: relaxed inertial stochastic forward-backward-forward
  *[SEG]: stochastic extragradient
  *[SFB]: stochastic forward-backward