# A relaxed-inertial forward-backward-forward algorithm for Stochastic
Generalized Nash equilibrium seeking

Shisheng Cui1 Barbara Franci2 Sergio Grammatico2 Uday V. Shanbhag1 and Mathias
Staudigl3 *This work was partially supported by NWO under research projects
OMEGA (613.001.702) and P2P-TALES (647.003.003), and by the ERC under research
project COSMOS (802348). Mathias Staudigl’s research benefited from the
support of the FMJH Program PGMO and from the support of EDF.1 Shisheng Cui
and Uday V. Shanbhag are with the Department of Industrial and Manufacturing
Engineering, Pennsylvania State University, University Park, PA 16802, USA
<EMAIL_ADDRESS>Barbara Franci and Sergio Grammatico are with the
Delft Center for System and Control, TU Delft, The Netherlands,
<EMAIL_ADDRESS>Mathias Staudigl is with the Department of
Data Science and Knowledge Engineering, Maastricht University, The
Netherlands<EMAIL_ADDRESS>
###### Abstract

In this paper we propose a new operator splitting algorithm for distributed
Nash equilibrium seeking under stochastic uncertainty, featuring relaxation
and inertial effects. Our work is inspired by recent deterministic operator
splitting methods, designed for solving structured monotone inclusion
problems. The algorithm is derived from a forward-backward-forward scheme for
solving structured monotone inclusion problems featuring a Lipschitz
continuous and monotone game operator. To the best of our knowledge, this is
the first distributed (generalized) Nash equilibrium seeking algorithm
featuring acceleration techniques in stochastic Nash games without assuming
cocoercivity. Numerical examples illustrate the effect of inertia and
relaxation on the performance of our proposed algorithm.

##  I Introduction

A stochastic generalized Nash equilibrium problem (SGNEP) describes a subclass
of competitive multi-agent optimization problems in which local unilateral
minimization of an agent-specific expectation-valued cost function subject to
system-wide shared coupling constraints. Due to the presence of the
uncertainty and the shared constraints, computing a SGNE is generally rather
challenging. However, these problems have recently received the attention of
the system and control community, especially because of their applicability
[1, 2, 3, 4, 5] to relevant problems in the engineering sciences. An important
class of models in this context is that of networked Cournot games with market
capacity constraints and uncertainty in demand and capacity [6]. Instances of
these models arise in transportation systems, where the drivers’ perception of
travel-time is a possible source of uncertainty [7], electricity markets where
companies dispatch electricity without an a priori knowledge of actual demand
[8], and natural gas markets where the companies participate in a bounded
capacity market [9].

Typically, two key concerns arise in any attempt to deal with uncertainty.
First, often the distribution of the random noise is not known to the agent,
thus making the computation of the cost function impossible. Second, even if
the distribution of the stochastic uncertainty is known or predictable from,
say, historical data, a key complication arises when trying to compute the
expected value (and its gradients). Costly simulation-based integration
techniques required employment each time an agent is asked to compute its
decision variable, imposing significant computational burden on each agent. A
versatile alternative to such approaches is provided by stochastic
approximation theory (SA), under which agents’ draw fresh samples at each
iteration.

With the aim of boosting the performance of distributed Nash equilibrium
seeking algorithms, Yi and Pavel [10] introduced a preconditioned forward-
backward splitting with inertial effects in a completely deterministic
environment where agents receive perfect feedback information. The possibility
of including inertia in the basic forward-backward scheme has received some
attention in the field already before (see e.g. [11, 12, 13, 14]). The common
motivation of all these contributions is to exploit momentum effects to
accelerate the numerical schemes, inspired by Nesterov’s accelerated gradient
method [15] for convex optimization. However, in the context of distributed
computation of Nash equilibria, the role of inertial and acceleration effects
is not well understood. This applies in particular to situations where the
game data are subject to stochastic uncertainty so that agents have only noisy
information available in their decision-making process. Even in the most
general problem where one’s aim is to solve a stochastic monotone inclusion
[16, 17, 18, 19], standard acceleration techniques have not received much
attention. Our aim is to shed some light on this highly understudied question
and prove some interesting properties about accelerated game dynamics. This
paper departs from recent progress made in the field of splitting algorithms
for stochastic variational problems, summarized in [20, 21] and [19], which
contain new asymptotic and non-asymptotic results on stochastic sampling-based
algorithms under weaker hypothesis than usually assumed in the computational
game theory literature. In these seminal contributions, stochastic versions of
Tseng’s modified extragradient (hitherto forward-backward-forward) algorithm
[22, 23] have been introduced. The importance of this alternative splitting
technique in the context of distributed Nash equilibrium seeking has been
emphasized in [24]. This work extends all these seminal contributions via an
explicit study of the effects of acceleration parameters. The numerical scheme
presented in this paper is provably convergent (in an almost sure sense)
without assuming co-coercivity of the game operator, and can be implemented
via a disciplined mini-batch stochastic approximation technology, distributed
over a network of competing agents. The main result of this work gives a
precise set of parameter sequences ensuring convergence of the game play to
the set of _variational equilibria_ , an important subset of generalized Nash
equilibria with a clear economic interpretation [25]. Our work extends the
recent results reported in [26], reliant on forward-backward splitting ideas,
and thus require co-coercivity, as well as the stochastic extragradient method
introduced in [27], where no joint coupling constraints are considered.

###  I-A Basic Notation

${\mathbb{R}}$ denotes the set of real numbers and
$\bar{\mathbb{R}}={\mathbb{R}}\cup\\{+\infty\\}$.
$\langle\cdot,\cdot\rangle:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{R}}$
denotes the standard inner product and $\lVert\cdot\rVert$ represents the
associated Euclidean norm. We indicate a (symmetric and) positive definite
matrix $A$, i.e., $x^{\top}Ax>0$, with $A\succ 0$. Given a matrix $\Phi\succ
0$, we define the $\Phi$-induced inner product as $\langle
x,y\rangle_{\Phi}=\langle\Phi x,y\rangle$ and the norm as $\lVert
x\rVert_{\Phi}=\sqrt{\langle\Phi x,x\rangle}$. $A\otimes B$ indicates the
Kronecker product between matrices $A$ and $B$. ${\bf{0}}_{m}$
(${\bf{1}}_{m}$) indicates the vector with $m$ entries all equal to $0$ ($1$).
Given $x_{1},\ldots,x_{N}\in{\mathbb{R}}^{n}$,
$\bm{x}:=\operatorname{col}\left(x_{1},\dots,x_{N}\right)=\left[x_{1}^{\top},\dots,x_{N}^{\top}\right]^{\top}.$

Let $T:{\mathbb{R}}^{n}\rightrightarrows{\mathbb{R}}^{n}$ be a set-valued
operator. The domain of $T$ are defined by
$\operatorname{dom}T=\\{x\in{\mathbb{R}}^{n}\mid T(x)\neq\varnothing\\}.$ The
set of zeros of $T$ is $\mathsf{Zer}(T)=\\{x\in{\mathbb{R}}^{n}\mid 0\in
T(x)\\}$. The resolvent of the operator $T$ is
$\mathrm{J}_{T}=(\operatorname{Id}+T)^{-1}$, where $\operatorname{Id}$
indicates the identity operator. An operator $T$ is monotone if $\langle
T(x)-T(y),x-y\rangle\geq 0$ and it is Lipschitz continuous if, for some
$\beta>0$, $\|T(x)-T(y)\|\leq\beta\|x-y\|$ for all
$x,y\in\operatorname{dom}T$. A monotone operator is maximally monotone if its
graph is not properly contained in the graph of another monotone operator.

Given a proper, lower semi-continuous, and convex function $g$, the
subdifferential is the operator $\partial g(x):=\\{u\in\Omega\mid(\forall
y\in\Omega):\langle y-x,u\rangle+g(x)\leq g(y)\\}$. The proximal operator is
defined as
$\operatorname{prox}_{g}(v):=\operatorname{argmin}_{u\in\Omega}\\{g(u)+\tfrac{1}{2}\lVert
u-v\rVert^{2}\\}=\mathrm{J}_{\partial g}(v)$. $\iota_{C}$ is the indicator
function of the set $C$, i.e., $\iota_{C}(x)=1$ if $x\in C$ and
$\iota_{C}(x)=0$ otherwise. The set-valued mapping
$\mathrm{N}_{C}:{\mathbb{R}}^{n}\rightrightarrows{\mathbb{R}}^{n}$ denotes the
normal cone operator for the the set $C$ , i.e.,
$\mathrm{N}_{C}(x)=\varnothing$ if $x\notin
C,\left\\{v\in{\mathbb{R}}^{n}|\sup_{z\in C}v^{\top}(z-x)\leq 0\right\\}$
otherwise.

All randomness is modeled on a complete probability space
$(\Omega,\mathcal{F},\mathbb{P})$, endowed with a filtration
${\mathbb{F}}=(\mathcal{F}_{k})_{k\geq 0}$.

##  II Mathematical Setup

###  II-A Generalized Nash equilibrium problems

We consider a game where each agent $i\in\mathcal{I}=\\{1,\ldots,N\\}$ chooses
an action $u_{i}\in\mathbb{R}^{d_{i}}$. Let ${\bf
u}=\operatorname{col}(u_{1},\ldots,u_{N})$ and $d\equiv\sum_{i=1}^{N}d_{i}$.
Each agent $i$ has a local cost function
${\mathbb{J}}_{i}:\mathbb{R}^{d}\to\bar{\mathbb{R}}$ of the form

${\mathbb{J}}_{i}(u_{i},{\bf u}_{-i})=f_{i}(u_{i},{\bf
u}_{-i})+g_{i}(u_{i}).\vspace{-.15cm}$ (1)

where ${\bf u}_{-i}=\text{col}(\\{u_{j}\\}_{j\neq i})$ is the vector of all
decision variables except for $u_{i}$, and
$g_{i}:{\mathbb{R}}^{d_{i}}\to\bar{\mathbb{R}}$ is a local idiosyncratic cost
function. The function ${\mathbb{J}}_{i}$ in (1) has the typical splitting
into smooth and non-smooth parts.

###### Standing Assumption 1

For each $i\in\mathcal{I}$, the function $g_{i}$ in (1) is proper, convex and
lower semi-continuous and
$\operatorname{dom}(g_{i})=\mathsf{U}_{i}\subseteq{\mathbb{R}}^{d_{i}}$ is
(nonempty) compact and convex. $\square$

Examples for the nonsmooth part are indicator functions to enforce local
constraints, or penalty functions that promote sparsity, or other desirable
structure.

We assume that the function $f_{i}(u_{i},{\bf u}_{-i})$ depends on the own
action $u_{i}$ and a subset of the others actions
$\\{u_{j}\\}_{j\in\mathcal{N}_{i}^{A}}$, where the set
$\mathcal{N}_{i}^{A}\subset\mathcal{I}$ is the _interaction neighborhood_ of
agent $i$. Furthermore, we assume convexity and differentiability, as usual in
the generalized Nash equilibrium problem (GNEP) literature [28, 29, 1].

###### Standing Assumption 2

For each $i\in\mathcal{I}$ and for all ${\bf u}_{-i}$, the function
$f_{i}(\cdot,{\bf u}_{-i})$ in (1) is convex and continuously differentiable.
$\square$

We assume that the game displays joint convexity with affine coupling
constraints defining the collective feasible set

$\bm{\mathcal{C}}=\\{{\bf u}\in\mathsf{U}\mid\;D{\bf
u}-b\leq\mathbf{0}_{m}\\},$ (2)

where $\mathsf{U}=\mathsf{U}_{1}\times\dots\times\mathsf{U}_{N}$,
$D=[D_{1}\mid\dots\mid D_{N}]\in\mathbb{R}^{m\times d}$ and
$b=\sum_{i=1}^{N}b_{i}\in\mathbb{R}^{m}$. Each matrix
$D_{i}\in\mathbb{R}^{m\times d_{i}}$ defines how agent $i$ is involved in the
coupling constraints. Given the strategies of all other agents ${\bf u}_{-i}$,
the set of feasible actions of player $i$ is defined as the following set-
valued map.

$\mathcal{C}_{i}({\bf u}_{-i})=\\{u_{i}\in\mathsf{U}_{i}\mid
D_{i}u_{i}-b_{i}\leq\sum_{j\neq i}^{N}(b_{j}-D_{j}u_{j})\\}.\vspace{-.15cm}$
(3)

###### Standing Assumption 3

The global feasible set $\bm{\mathcal{C}}$ in (2) satisfies Slater’s
constraint qualification. $\square$

For $i\in\mathcal{I}$, the $i$th agent solves the following parametrized
optimization problem.

$\forall
i\in\mathcal{I}:\quad\left\\{\begin{array}[]{cl}\min\limits_{u_{i}\in{\mathbb{R}}^{d_{i}}}&{\mathbb{J}}_{i}(u_{i},{\bf
u}_{-i})\\\ \text{ s.t. }&u_{i}\in\mathcal{C}_{i}({\bf
u}_{-i}).\end{array}\right.\vspace{-.15cm}$ (4)

The usual solution concept for the game with coupling constraints in (4) is
that of _generalized Nash equilibrium_ (GNE) [29, 30], i.e., an $N$-tuple
${\bf
u}^{\ast}=\textnormal{col}(u_{1}^{\ast},\ldots,u_{N}^{\ast})\in\mathsf{U}$
such that for all $i\in\mathcal{I}$,

${\mathbb{J}}_{i}(u_{i}^{\ast},{\bf
u}^{\ast}_{-i})\leq\inf\\{{\mathbb{J}}_{i}(u_{i},{\bf
u}^{\ast}_{-i})\,\mid\,u_{i}\in\mathcal{C}_{i}({\bf u}_{-i})\\}.$

Our computational approach for solving the GNEP in (4) makes use of the
Karush-Kuhn-Tucker (KKT) conditions characterizing the unilateral optimization
of the agents. To achieve a numerically tractable framework, we impose some
conditions on the model concerning the monotonicity and Lipschitz continuity
of the mapping that collects the local pseudogradients of the agents.

###### Standing Assumption 4

The pseudogradient mapping

$F({\bf u})=\mathrm{col}\left(\nabla_{u_{1}}f_{1}({\bf
u}),\ldots,\nabla_{u_{N}}f_{N}({\bf u})\right)$ (5)

is monotone and $\ell$-Lipschitz continuous. $\square$

The KKT conditions corresponding to the game in (4) are necessary and
sufficient for characterizing a tuple of strategies to be a GNE. Among all
possible GNEs of the game, we focus on the computation of _variational
equilibria_ (v-GNE), i.e. a GNE in which all agents share consensus on the
dual variables [29, Theorem 3.1], [31, Theorem 3.1] which is, in turn, a
solution of the variational system

$\forall
i\in\mathcal{I}:\begin{cases}\mathbf{0}_{d_{i}}\in\nabla_{u_{i}}f_{i}(u^{\ast}_{i},{\bf
u}^{\ast}_{-i})+\partial g_{i}(u^{*}_{i})+D_{i}^{\top}\lambda^{\ast}\\\
\mathbf{0}_{m}\in\operatorname{N}_{\mathbb{R}^{m}_{\geq
0}}(\lambda^{\ast})-(D{\bf u}^{\ast}-b).\\\ \end{cases}\vspace{-.15cm}$ (6)

for some $\lambda^{*}\in{\mathbb{R}}^{M}_{\geq 0}$.

###  II-B Distributed GNE via operator splitting

A key challenge one faces in any computational approach in Nash equilibrium
problems is to resolve the question how players access the decision variables
of the other agents. An attractive approach for resolving this issue is the
distributed operator splitting approach pioneered in [10].

We allow each agent to have information on his own local problem data only,
i.e., ${\mathbb{J}}_{i}$, $\mathsf{U}_{i}$, $D_{i}$ and $b_{i}$. Moreover,
each agent $i$ controls its local decision $u_{i}$ and a local copy
$\lambda_{i}\in\mathbb{R}^{m}_{\geq 0}$ of dual variables, as well as a local
auxiliary variable $\mu_{i}\in\mathbb{R}^{m}$ used to enforce consensus of the
dual variables. To reach such consensus, we let the agents exchange
information via an undirected weighted communication graph represented by its
weighted adjacency matrix $\mathbf{W}=[w_{i,j}]\in\mathbb{R}^{N\times N}$. We
assume $w_{ij}>0$ iff $(i,j)$ is an edge in the communication graph. The set
of neighbors of agent $i$ in the communication graph is
$\mathcal{N}^{\lambda}_{i}=\\{j\mid w_{i,j}>0\\}$.

###### Standing Assumption 5

The adjacency matrix $\mathbf{W}$ of the communication graph is symmetric and
irreducible. $\square$

Let us define the weighted Laplacian as
$\mathbf{L}=\operatorname{diag}\left\\{(\mathbf{W}\mathbf{1}_{N})_{1},\dots,(\mathbf{W}\mathbf{1}_{N})_{N}\right\\}-\mathbf{W}$.
It holds that $\mathbf{L}^{\top}=\mathbf{L}$,
$\operatorname{null}(\mathbf{L})=\\{a\bm{1}_{N},a\in{\mathbb{R}}\\}$ and that,
given Standing Assumption 5, $\mathbf{L}$ is positive semi-definite with real
and distinct eigenvalues $0=s_{1}<s_{2}\leq\ldots\leq s_{N}$. Moreover, given
the maximum (weighted) degree of the graph,
$\Delta:=\max_{i\in\mathcal{I}}(\mathbf{W}\mathbf{1}_{N})_{i}$, it holds that
$\Delta\leq s_{N}\leq 2\Delta$. Denoting by $\kappa=\lvert\mathbf{L}\rvert$,
it holds that $\kappa\leq 2\Delta$ [32]. We define the tensorized Laplacian as
the matrix $\bar{\mathbf{L}}=\mathbf{L}\otimes\mathbf{I}_{m}$. We set
$\bar{b}=(b_{1},\ldots,b_{N})^{\top}$, ${\bf
u}=\text{col}(u_{1},\ldots,u_{N})$ and similarly ${\bm{\mu}}$ and
${\bm{\lambda}}$. As the state variable, we consider the triple ${\bf x}=({\bf
u},{\bm{\mu}},{\bm{\lambda}})\in\mathsf{X}:=\mathbb{R}^{n}\times\mathbb{R}^{mN}\times\mathbb{R}^{mN}$
and endow $\mathsf{X}$ with the product topology. Let
$\mathbf{D}=\operatorname{diag}\\{D_{1},\ldots,D_{N}\\}$. Then, we define the
maximally monotone operators

$\displaystyle V({\bf x})$ $\displaystyle=\left[\begin{array}[]{c}F({\bf
u})+\mathbf{D}^{\top}{\bm{\lambda}}\\\ \bar{\mathbf{L}}{\bm{\lambda}}\\\
\bar{b}+\bar{\mathbf{L}}({\bm{\lambda}}-{\bm{\mu}})-\bar{\mathbf{D}}{\bf
u}\end{array}\right],$ (10) $\displaystyle T({\bf x})$ $\displaystyle=G({\bf
u})\times\\{\mathbf{0}_{Nm}\\}\times\operatorname{N}_{\mathbb{R}^{mN}_{\geq
0}}({\bm{\lambda}}),$ (11)

where $G({\bf u})=\partial g_{1}(u_{1})\times\cdots\times\partial
g_{N}(u_{N})$. Let us summarize the properties of the operators above.

###### Lemma 1

The following statements hold:

  * (i)

$V:\mathsf{X}\to\mathsf{X}$ is maximally monotone and
$\ell_{V}=(\ell+2\kappa+\lvert\mathbf{D}\rvert)$-Lipschitz continuous.

  * (ii)

$T:\mathsf{X}\rightrightarrows\mathsf{X}$ is maximally monotone.

###### Proof:

(i) We can split the operator $V$ into the parts $V_{1}({\bf
x})=\text{col}(F({\bf
u}),{\bf{0}}_{mN},\bar{b}+\bar{\mathbf{L}}{\bm{\lambda}})$ and $V_{2}({\bf
x})=\text{col}(\mathbf{D}^{\top}{\bm{\lambda}},\bar{\mathbf{L}}{\bm{\lambda}},-\mathbf{D}{\bf
u}-\bar{\mathbf{L}}{\bm{\mu}})$, which are maximally monotone by [23, Prop.
20.23], [23, Cor. 20.28]. Furthermore, similarly to [24, Lemma 1] $V_{1}$ is
$\ell_{1}=(\ell+\kappa)$-Lipschitz continuous and $V_{2}$ is
$\ell_{2}=(|\mathbf{D}|+\kappa)$-Lipschitz continuous. Hence, $V$ is
$\ell_{1}+\ell_{2}=\ell_{V}$-Lipschitz continuous.

(ii) It follows from [10, Lemma 5], [24, Lemma 1]. ∎

The splitting $V+T$ encodes a distributed version of the KKT conditions for
v-GNE (4). In particular, it can be shown that the zeros of the maximally
monotone inclusion $V+T$ are in correspondence with variational equilibria of
the Nash game.

###### Proposition 1

The set $\operatorname{zer}(V+T)$ coincides with the set of v-GNE of the game
satisfying the KKT conditions in (6).

###### Proof:

This follows from [10, Thm. 2] or [24, Lemma 3]. ∎

###  II-C Stochastic GNEPs

Stochastic uncertainty affecting the decision problem of agent $i$ is modeled
by a random variable $\xi_{i}:\Omega\to\Xi_{i}$, where
$\Xi_{i}\subset\mathbb{R}^{q_{i}}$ is a given measurable set. We assume that
the uncertainty enters the model in the smooth part of the agents’
optimization problem, i.e., for each $i\in\mathcal{I}$

${\mathbb{J}}_{i}(u_{i},{\bf u}_{-i})=\mathbb{E}[\hat{f}_{i}({\bf
u},\xi_{i})]+g_{i}(u_{i}).$ (12)

It follows that the local optimization problems in (4) describes a stochastic
programming problem, parameterized by the decisions of the opponents ${\bf
u}_{-i}$.

Let $k\in\mathbb{N}$ denote the iteration count of our computational
procedure. We assume that at each round $k$ agent $i$ is able to generate a
random sample $\xi_{i,k}=(\xi_{i,k}^{(t)})_{t=1}^{S_{k}}$, consisting of i.i.d
copies of the random element $\xi_{i}$. This sample is used to construct an
agent-specific gradient estimator of the form

$\hat{F}_{i,k}({\bf
u},\xi_{i,k})=\frac{1}{S_{k}}\sum_{t=1}^{S_{k}}\nabla_{u_{i}}\hat{f}_{i}({\bf
u},\xi^{(t)}_{i,k}),$ (13)

where $S_{k}\geq 1$ is the size of the data sample. (13) is an example of a
mini-batch estimator, which interpolates between cheap sampling and precision.
The degree of precision is regulated via the batch size sequence
$\\{S_{k}\\}_{k}$. Dynamically adjusting the size of the batch simulates an
online variance reduction mechanism, which plays a key role in our convergence
analysis of the distributed operator splitting algorithm to come. Mini-batch
samples are prominent in simulation-based optimization, where taking repeated
samples of stochastic gradients is computationally cheap [33, 27, 21, 19, 34].

###### Standing Assumption 6

The batch size $(S_{k})_{k\geq 1}$ is increasing and such that
$\sum_{k\in{\mathbb{N}}}\frac{1}{S_{k}}<\infty.$ $\square$

Under the prevailing i.i.d. assumption, it holds true that
$\mathbb{E}[\hat{F}_{i,k}({\bf u},\xi_{i,k})|{\bf u}]=F_{i}({\bf u})$ for all
$i\in\mathcal{I}$ and all ${\bf u}\in\mathbb{R}^{d}$. Hence, the random
variable (13) is an _unbiased_ estimator of the individual payoff gradient at
each action profile ${\bf u}$. Upon defining the random operator

$\hat{V}_{k}({\bf x},\xi_{k})=\left[\begin{array}[]{c}\hat{F}_{k}({\bf
u},\xi_{k})+\mathbf{D}^{\top}{\bm{\lambda}}\\\
\bar{\mathbf{L}}{\bm{\lambda}}\\\
\bar{b}+\bar{\mathbf{L}}({\bm{\lambda}}-{\bm{\mu}})-\mathbf{D}{\bf
u}\end{array}\right],\\\ $ (14)

with $\xi_{k}=\operatorname{col}(\xi_{i,k})_{i\in\mathcal{I}}$, we see that
$\mathbb{E}[\hat{V}_{k}({\bf x},\xi_{k})|{\bf x}]=V({\bf x})$ for all ${\bf
x}=({\bf u},{\bm{\mu}},{\bm{\lambda}})\in\mathsf{X}$.

Fundamental to the analysis of stochastic approximation algorithms is the
control of the stochastic error, defined for all $k\in{\mathbb{N}}$ as

$\varepsilon_{k}({\bf x},\xi_{k})=\hat{V}_{k}({\bf x},\xi_{k})-V({\bf
x})\quad\forall{\bf x}\in\mathsf{X}.$ (15)

###### Standing Assumption 7

There exists $\sigma>0$ such that for all $k\in{\mathbb{N}}$, the stochastic
error is such that the following hold $\mathbb{P}$-a.s. .

$\displaystyle\mathbb{E}_{\mathbb{P}}[\varepsilon_{k}({\bf x},\xi_{k})\mid{\bf
x}]=0,\text{ and }$ (16)
$\displaystyle\mathbb{E}_{\mathbb{P}}[\lVert\varepsilon_{k}({\bf
x},\xi_{k})\rVert^{2}\mid{\bf x}]\leq\frac{\sigma^{2}}{S_{k}}.$ (17)

$\square$

###### Remark 1

Assumption 7 is rather mild and standard in stochastic optimization [35, 26]
while Condition (16) means that the random operator $\hat{V}_{k}({\bf
x},\xi_{k})$ is a conditionally unbiased estimator of $V({\bf x})$. Note that
(17) can be satisfied if the sequence of martingale difference errors
$\hat{F}_{k}({\bf u},\xi)-F({\bf u})$ satisfies a uniform variance bound [35,
26].

##  III A distributed algorithm

Algorithm 1 Distributed Relaxed Inertial Stochastic Forward Backward Forward
(RISFBF)

Initialization: $u_{i,0}\in\mathbb{R}^{d_{i}},\lambda_{i,0}\in\mathbb{R}_{\geq
0}^{m},$ and $\mu_{i,0}\in\mathbb{R}^{m}.$
Iteration $k$: Agent $i$
($1$) Perform inertia step:

$\displaystyle u^{in}_{i,k}$ $\displaystyle=u_{i,k}+\alpha(u_{i,k}-u_{i,k-1})$
$\displaystyle\mu^{in}_{i,k}$
$\displaystyle=\mu_{i,k}+\alpha(\mu_{i,k}-\mu_{i,k-1})$
$\displaystyle\lambda^{in}_{i,k}$
$\displaystyle=\lambda_{i,k}+\alpha(\lambda_{i,k}-\lambda_{i,k-1}).$

($2$) Receives $u^{in}_{j,k}$ for $j\in\mathcal{N}_{i}^{A}$,
$\lambda^{in}_{j,k}$ and $\mu^{in}_{j,k}$ for $j\in\mathcal{N}_{i}^{\lambda}$
and update

$\displaystyle
u^{md}_{i,k}=\operatorname{prox}_{\gamma_{i}g_{i}}[u^{in}_{i,k}-\gamma_{i}(\hat{F}_{i,k}({\bf
u}^{in}_{k},\xi_{k})+D_{i}^{\top}\lambda_{i,k})]$
$\displaystyle\mu^{md}_{i,k}=\mu^{in}_{i,k}+\sigma_{i}\sum\nolimits_{j}w_{i,j}(\lambda^{in}_{j,k}-\lambda^{md}_{i,k})$
$\displaystyle\lambda^{md}_{i,k}=\Pi_{\mathbb{R}^{m}_{\geq
0}}\\{\lambda^{in}_{i,k}+\tau_{i}(D_{i}u^{in}_{i,k}-b_{i})$
$\displaystyle\quad\quad+\tau\sum\nolimits_{j}w_{i,j}[(\mu^{in}_{i,k}-\mu^{in}_{j,k})-(\lambda^{in}_{i,k}-\lambda^{in}_{j,k})]\\}$

($3$) Receives $u^{md}_{j,k}$ for $j\in\mathcal{N}_{i}^{A}$,
$\lambda^{md}_{j,k}$ and $\mu^{md}_{j,k}$ for $j\in\mathcal{N}_{i}^{\lambda}$
and performs a relaxation step:

$\displaystyle
u_{i,k+1}=(1-\rho_{k})u^{in}_{i,k}+\rho_{k}[u^{md}_{i,k}+\gamma_{i}(\hat{F}_{i,k}({\bf
u}^{in}_{k},\xi_{i,k})+$ $\displaystyle\quad\quad\quad-\hat{F}_{i,k}({\bf
u}^{md}_{k},\eta_{i,k}))+\gamma_{i}D_{i}^{\top}(\lambda^{in}_{i,k}-\lambda^{md}_{i,k})$
$\displaystyle\mu_{i,k+1}=(1-\rho_{k})\mu_{i,k}^{in}+\rho_{k}\\{\mu^{md}_{i,k}+$
$\displaystyle\quad\sigma_{i}\sum\nolimits_{j}w_{i,j}[(\lambda^{in}_{i,k}-\lambda^{in}_{j,k})-(\lambda^{md}_{i,k}-\lambda^{md}_{j,k})]\\}$
$\displaystyle\lambda_{i,k+1}=(1-\rho_{k})\lambda^{in}_{i,k}+\rho_{k}\\{[\lambda^{md}_{i,k}+\tau_{i}D_{i}(u^{in}_{i,k}-u^{md}_{i,k})$
$\displaystyle\quad\quad\quad-\tau_{i}\sum\nolimits_{j\in\mathcal{N}_{i}^{\lambda}}w_{i,j}[(\mu^{md}_{i,k}-\mu^{md}_{j,k})-(\mu^{in}_{i,k}-\mu^{in}_{j,k})]$
$\displaystyle\quad\quad\quad+\tau_{i}\sum\nolimits_{j\in\mathcal{N}_{i}^{\lambda}}w_{i,j}[(\lambda^{in}_{i,k}-\lambda^{in}_{j,k})-(\lambda^{md}_{i,k}-\lambda^{md}_{j,k})]\\}$

With the intent of boosting the convergence of distributed Nash seeking
algorithms, we propose a relaxed inertial forward-backward-forward algorithm
(RISFBF), presented in Algorithm 1. Using operator-theoretic notation, the
numerical scheme can be restated compactly as

$\left\\{\begin{array}[]{l}Z_{k}=X_{k}+\alpha_{k}(X_{k}-X_{k-1}),\\\
Y_{k}=\mathrm{J}_{\Psi^{-1}T}(Z_{k}-\Psi^{-1}\hat{V}_{k}(Z_{k},\xi_{k})),\\\
X_{k+1}=(1-\rho_{k})Z_{k}+\\\
\qquad\rho_{k}[Y_{k}-\Psi^{-1}(\hat{V}_{k}(Y_{k},\eta_{k})-\hat{V}_{k}(Z_{k},\xi_{k}))],\end{array}\right.$
(18)

where $Z_{k}=({\bf u}^{in}_{k},{\bm{\mu}}^{in}_{k},{\bm{\lambda}}_{k}^{in})$,
$X_{k}=({\bf u}_{k},{\bm{\mu}}_{k},{\bm{\lambda}}_{k})$ and $Y_{k}=({\bf
u}^{md}_{k},{\bm{\mu}}^{md}_{k},{\bm{\lambda}}_{k}^{md})$. The random sequence
$\eta_{k}=(\eta_{i,k})_{i\in\mathcal{I}}$ is another i.i.d. random sample,
generated independently by each agent after the first updating step in
Algorithm 1 is completed. The preconditioning matrix

$\Psi=\operatorname{diag}({\bm{\gamma}}^{-1},{\bm{\sigma}}^{-1},{\bm{\tau}}^{-1})$
(19)

collects all agent-specific step sizes, so that
${\bm{\gamma}}=\operatorname{diag}\\{\gamma_{1}\mathbf{I}_{d_{1}},\ldots,\gamma_{N}\mathbf{I}_{d_{N}}\\}$
is a block-diagonal matrix with $\gamma_{i}>0$ (analogously, ${\bm{\sigma}}$
and ${\bm{\tau}}$).

The iterations involve first an inertial step in the primal-dual space. Then,
there is a proximal step corresponding to a gradient-based update given the
stochastic estimate of the pseudogradient and the local estimate of the dual
variable, followed by a consensus-enforcing estimate merging the values of the
dual variables of neighboring agents, and a dual update in the spirit of
Lagrangian methods. The last step is a weighted average between the inertial
iterate $Z_{k}$ and the forward update $Y_{k}$. The algorithm is distributed
and involves communication only in terms of dual variables. This fact makes
the scheme very attractive for decentralized implementations in large
networked game-theoretic settings.

Standing Assumption 4 and Lemma 1 imply that $V$ is monotone and
$\ell_{V,\Psi}=\ell_{V}/\lambda_{\min}(\Psi)$-Lipschitz continuous in the
$\Psi$-induced norm [24].

###### Theorem 1

Suppose $\nu$ is a positive scalar where $0<\nu<1$. Let
$\lambda_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\min}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Psi})\in\left(0,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tfrac{1-\nu}{2\ell_{V}}}\right)$,
$0<\alpha_{k}\leq\bar{\alpha}<1$ and
$\rho_{k}={\tfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(3-\nu)}(1-\bar{\alpha})^{2}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}(2\alpha_{k}^{2}-\alpha_{k}+1)(1+\ell_{V,\Psi})}}$.
Then, the sequence $({\bf u}_{k})_{k\geq 1}$ generated by Algorithm 1
converges almost surely to a v-GNE of the game in (4).

###### Proof:

See Section V-B. ∎

###### Remark 2

The classic SFBF algorithm [21] is obtained as a special case by taking
$\alpha=0$ and $\rho_{k}=1$. Under noise-free feedback, this scheme would
coincide with the operator-splitting approach of [24].

##  IV Numerical Results

$M_{1}$$M_{2}$$M_{3}$$M_{4}$$M_{5}$$M_{6}$$M_{7}$$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$$C_{8}$$C_{9}$$C_{10}$
Figure 1: Networked Cournot game: an edge from $C_{i}$ to $M_{j}$ means that
company $i$ sells energy in market $j$.

In this section, we report the results of some numerical simulations to
illustrate the improved performance of the RISFBF algorithm (Algorithm 1)
compared to the classic SFBF [35, 24] and with the preconditioned SFB [26,
10].

Let us consider a networked Cournot problem with market capacity constraints
as, for instance, the electricity market or the gas market, inspired by [36].
We suppose that there are $N=10$ firms selling energy in $m=7$ markets. Not
every company sells quantities on each market. Instead, we let
$\mathcal{M}_{i}\subseteq\\{1,\ldots,m\\}$ denote the subset of markets firm
$i$ is active on. Each company has a cost of production
$c_{i}(u_{i})=c_{i}^{\top}u_{i}$ where $c_{i}\in{\mathbb{R}}^{d_{i}}$ is
chosen according to a truncated normal distribution, i.e.,
$[c_{i}]_{j}=\max(N(2,1),0.6)$. Moreover, each market $j$ has an inverse
demand function $P_{j}({\bf u},\xi)=q_{j}+p_{j}(\xi)[S_{j}({\bf u})]^{\sigma}$
where $q_{j}=400$ and $p_{j}(\xi)$ depends on the unknown random variable,
e.g, the overall demand. The values of $p_{j}(\xi)$ are randomly generated
with a normal distribution with mean $0.02$ and bounded variance. The variable
$S_{j}({\bf u})=\sum_{i\in\mathcal{I}}[u_{i}]_{j}$ couples the actions of the
companies and it represents the total energy sold in market $j$. Hence, the
cost function of each company is
${\mathbb{J}}_{i}(u_{i},\bm{u}_{-i})=c_{i}(u_{i})-\sum_{j\in\mathcal{M}_{i}}{\mathbb{E}}[P_{j}({\bf
u},\xi)[u_{i}]_{j}]$. The corresponding pseudogradient mapping is monotone,
according to [36, Section 4], for $1<\sigma\leq 3$. Therefore, we fix
$\sigma=1.2$. Moreover, we suppose that the companies have a limited
production, i.e., $0\leq[u_{i}]_{j}\leq\theta_{i,j}$ with
$\theta_{i,j}=\max(N(250,50),0)$ for $j\in\mathcal{M}_{i}$. This can be
incorporated by setting
$g_{i}(u_{i})=\sum_{j\in\mathcal{N}_{i}^{A}}\iota_{[0,\theta_{i,j}]}([u_{i}]_{j})$.
Similarly, the markets have a bounded capacity $b_{j}\in[5,10],j=1,\ldots,m,$
and the coupling between the companies can be retrieved from Figure 1.

Figure 2: Residual distance of the primal variable form the solution.

The plot in Figure 2 shows the performance, in terms of the residual, of our
proposed algorithm in comparison with the SFBF and SFB algorithms. The
residual mapping is defined as $\operatorname{res}(x^{k})=\lVert
x^{k}-\operatorname{proj}_{\mathcal{C}}(x^{k}-F(x^{k}))\rVert$ and it measures
the distance of the primal variable from being a Nash equilibrium. The thick
line indicates the average performance and the transparent area is the
variability over 10 simulations. The acceleration parameter is updated
according to $\alpha_{k}=\bar{\alpha}(1-\frac{1}{k+1})$ with
$\bar{\alpha}=0.1$ and the relaxation parameter is
$\rho_{k}=\frac{(3-\nu)(1-\bar{\alpha})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+\ell_{V,\Psi})}$,
with $\nu=0.01$. Figure 3 shows how the performance changes varying the
inertial parameter $\bar{\alpha}$ while $\rho_{k}=1$ is fixed. For the sake of
comparison, we also include the performance with the same parameters as in
Figure 2, the updating rule for $\rho_{k}$ as in Theorem 1 and the SFBF
($\alpha_{k}=0$, $\rho_{k}=0$).

##  V Analysis

###  V-A Preparatory facts

To simplify the analysis, let us define the random processes
$\mathcal{A}_{k}:=\hat{V}_{k}(Z_{k},\xi_{k})$ and
$\mathcal{B}_{k}:=\hat{V}_{k}(Y_{k},\eta_{k})$. Define the sub-sigma algebra
$\mathcal{F}_{k}:=\sigma({\bf
x}_{0},\xi_{0},\ldots,\xi_{k-1},\eta_{0},\ldots,\eta_{k-1})$, and
$\mathcal{G}_{k}:=\sigma\left(\mathcal{F}_{k}\cup\sigma(\xi_{k})\right)$. We
introduce the centered error processes
$U_{k}:=\mathcal{A}_{k}-\mathbb{E}[\mathcal{A}_{k}\mid\mathcal{F}_{k}]$ and
$W_{k}:=\mathcal{B}_{k}-\mathbb{E}[\mathcal{B}_{k}\mid\mathcal{G}_{k}]$. Note
that Standing Assumption 7 implies that
$\mathbb{E}[U_{k}\mid\mathcal{F}_{k}]=\mathbb{E}[W_{k}\mid\mathcal{F}_{k}]=0$
and that $\left(\mathbb{E}[\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}\mid\mathcal{F}_{k}]\right)_{k\geq 1}$ and
$\left(\mathbb{E}[\lVert
W_{k}\rVert^{2}_{\Psi^{-1}}\mid\mathcal{F}_{k}]\right)_{k\geq 1}$ are summable
sequences.

Define the residual function for the monotone inclusion as $r_{\Psi}(x)=\lVert
x-\mathrm{J}_{\Psi^{-1}T}(x-\Psi^{-1}V(x))\rVert.$ For every $\Psi\succ 0$,
$x\in\operatorname{zer}(V+T)\Leftrightarrow r_{\Psi}(x)=0$.

###### Lemma 2

For $x,y\in\mathsf{X}$ and $\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}.$

Figure 3: Residual distance of the primal variable form the solution for the
RISFBF algorithm varying the inertial parameter.

######  Lemma 3 (Robbins-Siegmund)

[37, Lemma 11, page 50]. Let
$(\Omega,\mathcal{F},{\mathbb{F}}=(\mathcal{F}_{k})_{k\geq 0},\mathbb{P})$ be
a discrete stochastic basis. Let $(\alpha_{k})_{k\in{\mathbb{N}}}$,
$(\theta_{k})_{k\in{\mathbb{N}}}$, $(\eta_{k})_{k\in{\mathbb{N}}}$ and
$(\chi_{k})_{k\in{\mathbb{N}}}$ be non-negative processes such that
$\sum_{k}\eta_{k}<\infty$, $\sum_{k}\chi_{k}<\infty$ and let

$\forall
k\in{\mathbb{N}},\quad{\mathbb{E}}[\alpha_{k+1}|\mathcal{F}_{k}]+\theta_{k}\leq(1+\chi_{k})\alpha_{k}+\eta_{k}\quad
a.s.$

Then $\sum_{k}\theta_{k}<\infty$ and $(\alpha_{k})_{k\in{\mathbb{N}}}$
converges a.s. to a non negative random variable.

###### Lemma 4

For all $k\geq 1$ we have

$-\lVert Z_{k}-Y_{k}\rVert^{2}_{\Psi}\leq\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}-\tfrac{1}{2}r^{2}_{\Psi}(Z_{k}).$ (20)

###### Proof:

By definition

$\displaystyle\tfrac{1}{2}r_{\Psi}^{2}(Z_{k})=\tfrac{1}{2}\lVert
Z_{k}-\mathrm{J}_{\Psi^{-1}T}(Z_{k}-\Psi^{-1}V(Z_{k})\rVert^{2}$
$\displaystyle\leq\lVert Z_{k}-Y_{k}\rVert^{2}_{\Psi}$
$\displaystyle+\lVert\mathrm{J}_{\Psi^{-1}T}(Z_{k}-\Psi^{-1}\mathcal{A}_{k})-\mathrm{J}_{\Psi^{-1}T}(Z_{k}-\Psi^{-1}V(Z_{k}))\rVert^{2}_{\Psi}$
$\displaystyle\leq\lVert Z_{k}-Y_{k}\rVert^{2}_{\Psi}+\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}$

where the last inequality uses the non-expansiveness of the resolvent operator
$J_{\Psi^{-1}T}$ under the norm $\lVert\cdot\rVert_{\Psi}$. ∎

###  V-B Convergence analysis of RISFBF algorithm

Define the stochastic processes

$\displaystyle\Delta
M_{k}:=\textstyle{\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(3-\nu)}\rho_{k}}{1+\ell_{V,\Psi}}\lVert\mathtt{e}_{k}\rVert^{2}+\textstyle{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu}\rho_{k}}\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}},$ (21) $\displaystyle\Delta
N_{k}(p):=2\rho_{k}\langle W_{k},p-Y_{k}\rangle,$ (22)

with ${\mathtt{e}_{k}}:=W_{k}-U_{k}$. We start proving the following
fundamental inequality.

######  Lemma 5 (Fundamental Recursion)

Fix $p\in\operatorname{zer}(V+T)$ arbitrary, and set
$R_{k}=Y_{k}+\Psi^{-1}(\mathcal{A}_{k}-\mathcal{B}_{k})$. For all $k\geq 0$,
it holds true that

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}\leq(1+\alpha_{k})\lVert
X_{k}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert X_{k-1}-p\rVert^{2}_{\Psi}$
$\displaystyle+\Delta M_{k}+\Delta
N_{k}(p)-\textstyle{\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu}\rho_{k}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}}r^{2}_{\Psi}(Z_{k})$
$\displaystyle+\alpha_{k}\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}\textstyle{\left(2\alpha_{k}+\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}3-\nu}(1-\alpha_{k})}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}\rho_{k}(1+\ell_{V,\Psi})}\right)}$
$\displaystyle-(1-\alpha_{k})\textstyle{\left(\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}3-\nu}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}\rho_{k}(1+\ell_{V,\Psi})}-1\right)}\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}.$

###### Proof:

Start by observing

$\displaystyle\lVert Z_{k}-p\rVert^{2}_{\Psi}$ $\displaystyle=\lVert
Z_{k}-Y_{k}+Y_{k}-R_{k}+R_{k}-p\rVert^{2}_{\Psi}$ $\displaystyle=\lVert
Z_{k}-Y_{k}\rVert^{2}_{\Psi}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\
-\ }\lVert Y_{k}-R_{k}\rVert^{2}_{\Psi}+\lVert R_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+2\langle Z_{k}-R_{k},Y_{k}-p\rangle_{\Psi}.$

Since

$\displaystyle\lVert Y_{k}-R_{k}\rVert^{2}_{\Psi}$
$\displaystyle=\lVert\Psi^{-1}(\mathcal{A}_{k}-\mathcal{B}_{k})\rVert^{2}_{\Psi}$
$\displaystyle=\lVert V(Y_{k})-V(Z_{k})+W_{k}-U_{k}\rVert^{2}_{\Psi^{-1}}$
$\displaystyle\leq 2\ell^{2}_{V,\Psi}\lVert
Y_{k}-Z_{k}\rVert^{2}_{\Psi}+2\lVert W_{k}-U_{k}\rVert^{2}_{\Psi^{-1}}$
$\displaystyle=2\ell^{2}_{V,\Psi}\lVert
Y_{k}-Z_{k}\rVert^{2}_{\Psi}+2\lVert{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(23)

Hence,

$\displaystyle\lVert Z_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overset{\eqref{bd-
YkRk}}{\geq}}(1-2\ell^{2}_{V,\Psi})\lVert
Z_{k}-Y_{k}\rVert^{2}_{\Psi}-2\lVert{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathtt{e}_{k}}\rVert^{2}_{\Psi^{-1}}$
$\displaystyle+\lVert R_{k}-p\rVert^{2}_{\Psi}+2\langle
Z_{k}-R_{k},Y_{k}-p\rangle_{\Psi}.$

Using the definition
$Y_{k}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\triangleq\
}\mathrm{J}_{\Psi^{-1}T}(Z_{k}-\Psi^{-1}\mathcal{A}_{k})$, we get

$\displaystyle Y_{k}+\Psi^{-1}T(Y_{k})$
$\displaystyle\ni(Z_{k}-\Psi^{-1}\mathcal{A}_{k})$ $\displaystyle\mbox{ or
}T(Y_{k})$ $\displaystyle\ni\Psi(Z_{k}-Y_{k}-\Psi^{-1}\mathcal{A}_{k})$

Since $p\in\operatorname{zer}(T+V)$, $(p,{\bf 0})\in\mbox{gr}(T+V)$, implying
that ${\bf 0}-V(p)\in T(p)$. Consequently, by monotonicity of $T$, we have
that

$\displaystyle\langle\Psi(Z_{k}-Y_{k}-\Psi^{-1}\mathcal{A}_{k})+V(p),Y_{k}-p\rangle\geq
0$ $\displaystyle\mbox{ or }\langle
Z_{k}-Y_{k}-\Psi^{-1}(\mathcal{A}_{k}-\mathcal{B}_{k}),Y_{k}-p\rangle_{\Psi}$
$\displaystyle\geq\langle
V(Y_{k})-V(p),Y_{k}-p\rangle+\langle\mathcal{B}_{k}-V(Y_{k}),Y_{k}-p\rangle.$

By definition, $R_{k}=Y_{k}+\Psi^{-1}(\mathcal{A}_{k}-\mathcal{B}_{k})$,
$W_{k}=\mathcal{B}_{k}-V(Y_{k})$, and
$V(Y_{k})=\mathbb{E}[\mathcal{B}_{k}\mid\mathcal{G}_{k}]$, we have that

$\displaystyle\langle Z_{k}-R_{k},Y_{k}-p\rangle_{\Psi}$
$\displaystyle\geq\langle V(Y_{k})-V(p),Y_{k}-p\rangle$ $\displaystyle+\langle
W_{k},Y_{k}-p\rangle$

Since $V$ is a monotone operator, this implies $\langle
Z_{k}-R_{k},Y_{k}-p\rangle_{\Psi}\geq\langle W_{k},Y_{k}-p\rangle.$ Whence,

$\displaystyle\lVert Z_{k}-p\rVert^{2}_{\Psi}\geq$
$\displaystyle(1-2\ell^{2}_{V,\Psi})\lVert Y_{k}-Z_{k}\rVert^{2}_{\Psi}+\lVert
R_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle-2\lVert{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathtt{e}_{k}}\rVert^{2}_{\Psi^{-1}}+2\langle
W_{k},Y_{k}-p\rangle.$

Rearranging, we arrive at

$\begin{array}[]{l}\lVert R_{k}-p\rVert^{2}_{\Psi}\leq\lVert
Z_{k}-p\rVert^{2}_{\Psi}+2\lVert{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathtt{e}_{k}}\rVert^{2}_{\Psi^{-1}}\\\
-(1-2\ell^{2}_{V,\Psi})\lVert Y_{k}-Z_{k}\rVert^{2}_{\Psi}+2\langle
W_{k},p-Y_{k}\rangle.\end{array}$ (24)

Next, we use Lemma 2 to arrive at

$\displaystyle\lVert
X_{k+1}-p\rVert^{2}_{\Psi}=\lVert(1-\rho_{k})Z_{k}+\rho_{k}R_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle=(1-\rho_{k})\lVert Z_{k}-p\rVert^{2}_{\Psi}+\rho_{k}\lVert
R_{k}-p\rVert^{2}_{\Psi}$ $\displaystyle-\rho_{k}(1-\rho_{k})\lVert
R_{k}-Z_{k}\rVert^{2}_{\Psi}$ $\displaystyle=(1-\rho_{k})\lVert
Z_{k}-p\rVert^{2}_{\Psi}+\rho_{k}\lVert
R_{k}-p\rVert^{2}_{\Psi}-\textstyle{\textstyle{\frac{1-\rho_{k}}{\rho_{k}}}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}}$ $\displaystyle\leq\lVert
Z_{k}-p\rVert^{2}_{\Psi}-\textstyle{\frac{1-\rho_{k}}{\rho_{k}}}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}+2\lambda^{2}\rho_{k}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}$
$\displaystyle-\rho_{k}(1-2\ell^{2}_{V,\Psi})\lVert
Z_{k}-Y_{k}\rVert^{2}_{\Psi}-2\rho_{k}\langle W_{k},Y_{k}-p\rangle$
$\displaystyle=\lVert
Z_{k}-p\rVert^{2}_{\Psi}-\textstyle{\textstyle{\frac{1-\rho_{k}}{\rho_{k}}}}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}-2\rho_{k}\langle W_{k},Y_{k}-p\rangle$
$\displaystyle-\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(1-\nu)}-2\ell^{2}_{V,\Psi})\lVert
Z_{k}-Y_{k}\rVert^{2}_{\Psi}-\textstyle{\textstyle{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu\rho_{k}}}\lVert
Y_{k}-Z_{k}\rVert^{2}_{\Psi}}$
$\displaystyle+2\rho_{k}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}.$

Using (20), this implies

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}\leq\lVert
Z_{k}-p\rVert^{2}_{\Psi}-\textstyle{\frac{1-\rho_{k}}{\rho_{k}}}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}$
$\displaystyle-\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(1-\nu)}-2\ell^{2}_{V,\Psi})\lVert
Z_{k}-Y_{k}\rVert^{2}_{\Psi}-\textstyle{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\nu\rho_{k}}{2}}}r^{2}_{\Psi}(Z_{k})$
$\displaystyle-2\rho_{k}\langle
W_{k},Y_{k}-p\rangle+2\rho_{k}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}+\textstyle{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu}\rho_{k}}\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}.$

Furthermore,

$\displaystyle\tfrac{1}{\rho_{k}}\lVert X_{k+1}-Z_{k}\rVert_{\Psi}$
$\displaystyle=\lVert R_{k}-Z_{k}\rVert_{\Psi}$
$\displaystyle\leq\lVert\mathcal{B}_{k}-\mathcal{A}_{k}\rVert_{\Psi^{-1}}+\lVert
Y_{k}-Z_{k}\rVert_{\Psi}$ $\displaystyle\leq(1+\ell_{V,\Psi})\lVert
Y_{k}-Z_{k}\rVert_{\Psi}+\lVert\mathtt{e}_{k}\rVert_{\Psi^{-1}},$

which implies

$\displaystyle\tfrac{1}{2\rho^{2}_{k}}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}\leq(1+\ell_{V,\Psi})^{2}\lVert
Y_{k}-Z_{k}\rVert^{2}_{\Psi}+\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}.$

Multiplying both sides by
$\frac{\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi})}{1+\ell_{V,\Psi}}$,
we obtain

$\displaystyle\textstyle{\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi}}{2\rho_{k}}(1+\ell_{V,\Psi})}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}$
$\displaystyle\leq\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi})(1+\ell_{V,\Psi})\lVert
Y_{k}-Z_{k}\rVert^{2}_{\Psi}$
$\displaystyle+\textstyle{\frac{\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi})}{1+\ell_{V,\Psi}}}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}.$

Rearranging terms, and noting that
$({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi})(1+\ell_{V,\Psi})\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell^{2}_{V,\Psi}$,
the above estimate becomes

$\displaystyle-\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell^{2}_{V,\Psi})\lVert
Y_{k}-Z_{k}\rVert^{2}_{\Psi}$
$\displaystyle\leq-\textstyle{\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi}}{2\rho_{k}}(1+\ell_{V,\Psi})}\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}+\textstyle{\frac{\rho_{k}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi}\Psi)}{1+\ell_{V,\Psi}}}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}.$

Substituting this bound into the first majorization of the anchor process
$\lVert X_{k+1}-p\rVert^{2}_{\Psi}$, we see

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}\leq\lVert
Z_{k}-p\rVert^{2}_{\Psi}+\textstyle{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu}\rho_{k}}\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}-\textstyle{\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nu}\rho_{k}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}}r^{2}_{\Psi}(Z_{k})$
$\displaystyle-\left(\textstyle{\frac{1-\rho_{k}}{\rho_{k}}}+\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi}}{2\rho_{k}(1+\ell_{V,\Psi})}\right)\lVert
X_{k+1}-Z_{k}\rVert^{2}_{\Psi}$
$\displaystyle+\rho_{k}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}\textstyle{\left(2+\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1-\nu}-2\ell_{V,\Psi}}{1+\ell_{V,\Psi}}\right)}-2\rho_{k}\langle
W_{k},Y_{k}-p\rangle.$

Observe that

$\displaystyle\lVert X_{k+1}-Z_{k}\rVert^{2}_{\Psi}$
$\displaystyle\geq(1-\alpha_{k})\lVert X_{k+1}-X_{k}\rVert^{2}_{\Psi}$
$\displaystyle+(\alpha^{2}_{k}-\alpha_{k})\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi},$ (25) $\displaystyle\lVert
Z_{k}-p\rVert^{2}_{\Psi}$ $\displaystyle=(1+\alpha_{k})\lVert
X_{k}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert X_{k-1}-p\rVert^{2}_{\Psi}$
$\displaystyle+\alpha_{k}(1+\alpha_{k})\lVert X_{k}-X_{k-1}\rVert^{2}_{\Psi}.$
(26)

Choose parameters $\alpha_{k}$ and $\rho_{k}$ such that
$\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}3-\nu}-2\rho_{k}(1+\ell_{V,\Psi})}{2\rho_{k}(1+\ell_{V,\Psi})}>0$.
Then, using both of these relations in the last estimate for $\lVert
X_{k+1}-p\rVert^{2}_{\Psi}$, we arrive at

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}\leq(1+\alpha_{k})\lVert
X_{k}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert X_{k-1}-p\rVert^{2}_{\Psi}$
$\displaystyle+\alpha_{k}(1+\alpha_{k})\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}-2\rho_{k}\langle W_{k+1},Y_{k}-p\rangle$
$\displaystyle-\textstyle{\frac{\nu\rho_{k}}{2}}r^{2}_{\Psi}(Z_{k})+\frac{(3-\nu)\rho_{k}}{1+\ell_{V,\Psi}}\lVert\mathtt{e}_{k}\rVert^{2}_{\Psi^{-1}}+\textstyle{\nu\rho_{k}}\lVert
U_{k}\rVert^{2}_{\Psi^{-1}}$
$\displaystyle-\textstyle{\left(\frac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)}[(1-\alpha_{k})\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$
$\displaystyle+(\alpha^{2}_{k}-\alpha_{k})\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}].$

Using the respective definitions of the stochastic increments $\Delta
M_{k+1},\Delta N_{k}(p)$ in (21) and (22), we arrive at

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}\leq(1+\alpha_{k})\lVert
X_{k}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert X_{k-1}-p\rVert^{2}_{\Psi}$
$\displaystyle+\Delta M_{k}+\Delta
N_{k}(p)-\textstyle{\frac{\nu\rho_{k}}{2}}r^{2}_{\Psi}(Z_{k})$
$\displaystyle+\alpha_{k}\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}\textstyle{\left(2\alpha_{k}+\frac{(3-\nu)(1-\alpha_{k})}{2\rho_{k}(1+\ell_{V,\Psi})}\right)}$
$\displaystyle-(1-\alpha_{k})\textstyle{\left(\frac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)}\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$

∎

Rearranging the fundamental recursion, we see

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert
X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+(1-\alpha_{k})\textstyle{\left(\frac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)}\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$ $\displaystyle\leq\lVert
X_{k}-p\rVert_{\Psi}^{2}-\alpha_{k}\lVert X_{k-1}-p\rVert^{2}_{\Psi}+\Delta
M_{k}+\Delta N_{k}(p)$
$\displaystyle+(1-\alpha_{k})\textstyle{\left(\frac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)}\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}-\tfrac{\nu\rho_{k}}{2}r^{2}_{\Psi}(Z_{k})$
$\displaystyle+\textstyle{\left(2\alpha_{k}^{2}+(1-\alpha_{k})\left(1-\tfrac{(3-\nu)(1-\alpha_{k})}{2\rho_{k}(1+\ell_{V,\Psi})}\right)\right)}\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}.$

Suppose $(\alpha_{k})_{k}$ is a non-decreasing sequence satisfying
$0<\alpha_{k}\leq\bar{\alpha}<1$ and
$\rho_{k}={\tfrac{(3-\nu)(1-\bar{\alpha})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+\ell_{V,\Psi})}}.$
Since
$\rho_{k}\leq\tfrac{(3-\nu)(1-\alpha_{k})^{2}}{2(2\alpha_{k}^{2}-\alpha_{k}+1)(1+\ell_{V,\Psi})}$,
we claim that

$\displaystyle\lVert X_{k+1}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert
X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+(1-\alpha_{k})\left(\tfrac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$ $\displaystyle\geq\lVert
X_{k+1}-p\rVert^{2}_{\Psi}-\alpha_{k}\lVert X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+(1-\alpha_{k})\left(\tfrac{2\alpha_{k}^{2}-\alpha_{k}+1}{(1-\alpha_{k})^{2}}-1\right)\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}\geq 0.$

To see this, observe that for any $\alpha>0$,

$\displaystyle\lVert X_{k+1}-p\rVert_{\Psi}^{2}-\alpha\lVert
X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+(1-\alpha)\textstyle{\left(\frac{2\alpha^{2}-\alpha+1}{(1-\alpha)^{2}}-1\right)}\lVert
X_{k+1}-X_{k}\rVert_{\Psi}^{2}$ $\displaystyle>\lVert
X_{k+1}-p\rVert^{2}_{\Psi}-\alpha\lVert X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+(1-\alpha)\textstyle{\left(\frac{\alpha^{2}-\alpha+1}{(1-\alpha)^{2}}-1\right)}\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$ $\displaystyle=\lVert
X_{k+1}-p\rVert^{2}_{\Psi}-\alpha\lVert X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle-\left(\tfrac{\alpha^{2}-\alpha+1}{1-\alpha}-\tfrac{1-2\alpha+\alpha^{2}}{1-\alpha}\right)\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$ $\displaystyle=(\alpha+(1-\alpha))\lVert
X_{k+1}-p\rVert^{2}_{\Psi}-\alpha\lVert X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+\textstyle{\left(\alpha+\frac{\alpha^{2}}{1-\alpha}\right)}\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$ $\displaystyle\geq\alpha\lVert
X_{k+1}-p\rVert^{2}_{\Psi}+\alpha\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}-\alpha\lVert X_{k}-p\rVert^{2}_{\Psi}$
$\displaystyle+2\alpha\lVert X_{k+1}-p\rVert_{\Psi}\cdot\lVert
X_{k+1}-X_{k}\rVert_{\Psi}$ $\displaystyle=\alpha(\lVert
X_{k+1}-p\rVert_{\Psi}+\lVert X_{k+1}-X_{k}\rVert_{\Psi})^{2}-\alpha\lVert
X_{k}-p\rVert^{2}_{\Psi}$ $\displaystyle\geq\alpha\lVert
X_{k+1}-p+X_{k}-X_{k+1}\rVert^{2}_{\Psi}-\alpha\lVert
X_{k}-p\rVert^{2}_{\Psi}=0$

where the third inequality follows from Young’s inequality. Under this
specific coupling of the inertial and relaxation parameters, it holds that
$2\alpha_{k}^{2}+(1-\alpha_{k})\left(1-\tfrac{(3-\nu)(1-\alpha_{k})}{2\rho_{k}(1+\ell_{V,\Psi})}\right)\leq
0$. Now, let $H_{k}(p)=\lVert X_{k}-p\rVert_{\Psi}^{2}-\alpha_{k}\lVert
X_{k-1}-p\rVert^{2}_{\Psi}+(1-\alpha_{k})\left(\frac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}$, and
$\delta_{k}=\tfrac{\nu\rho_{k}}{2}r^{2}_{\Psi}(Z_{k})-\left(2\alpha_{k}^{2}+(1-\alpha_{k})\left(1-\tfrac{(3-\nu)(1-\alpha_{k})}{2\rho_{k}(1+\ell_{V,\Psi})}\right)\right)\lVert
X_{k}-X_{k-1}\rVert^{2}_{\Psi}.$ Then,
$(1-\alpha_{k+1})\left(\tfrac{3-\nu}{2\rho_{k+1}(1+\ell_{V,\Psi})}-1\right)\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}\leq(1-\alpha_{k})\left(\tfrac{3-\nu}{2\rho_{k}(1+\ell_{V,\Psi})}-1\right)\lVert
X_{k+1}-X_{k}\rVert^{2}_{\Psi}$. Therefore, for all $k\geq 0$, we conclude

$\displaystyle\mathbb{E}[H_{k+1}(p)\mid\mathcal{F}_{k}]\leq
H_{k}(p)-\delta_{k}(p)+\mathbb{E}[\Delta M_{k}\mid\mathcal{F}_{k}].$

Using Standing Assumption 7, we deduce that $\left(\mathbb{E}[\Delta
M_{k}\mid\mathcal{F}_{k}]\right)_{k\in\mathbb{N}}$ is summable, and thus we
can apply Lemma 3 to the above recursion. Hence, we readily deduce the
existence of an a.s. finite limiting random variable $H_{\infty}(p)$ such that
$\mathbb{P}\left(\lim_{k\to\infty}H_{k}(p)=H_{\infty}(p)\right)=1$ and
$\sum_{k\in{\mathbb{N}}}\delta_{k}(p))_{k\in\mathbb{N}}<\infty$. Therefore,
there exists a measurable set $\Omega_{0}\in\mathcal{F}$ with
$\mathbb{P}(\Omega_{0})=1$ such that for each $\omega\in\Omega_{0}$ it holds
true that

$\displaystyle\lim_{k\to\infty}\lVert
X_{k}(\omega)-X_{k-1}(\omega)\rVert_{\Psi}=0\text{ and }$
$\displaystyle\lim_{k\to\infty}\textstyle{\frac{\nu\rho_{k}}{2}}r^{2}_{\Psi}(Z_{k}(\omega))=0.$

Assuming that $\liminf_{k\to\infty}\rho_{k}>0$, we conclude that
$\lim_{k\to\infty}r^{2}_{\Psi}(Z_{k})=0$ $\mathbb{P}$-a.s. Therefore, we
conclude that $(X_{k})_{k\geq 1}$ converges a.s. to a limiting random variable
with values in $\operatorname{zer}(V+T)$.

##  VI Conclusion

In the context of shared constraint variants of stochastic generalized Nash
equilibrium problems, the convergence of the forward-backward-forward
algorithm can be boosted via an integrated acceleration-relaxation procedure.
In the presence of stochastic uncertainty, convergence can be proved assuming
only monotonicity and Lipschitz continuity of the expected-valued operator.
Specifically, our main result is the claim of almost sure global convergence
of the trajectory of actions to the set of variational equilibria. In future
research, we aim to investigate the question how one may relax monotonicity
and Lipschitz continuity assumptions even further, derive rate statements, and
examine how partial information may be introduced into the algorithm.

## References

  * [1] A. A. Kulkarni and U. V. Shanbhag, “On the variational equilibrium as a refinement of the generalized Nash equilibrium,” _Automatica_ , vol. 48, no. 1, pp. 45–55, 2012.
  * [2] U. Ravat and U. V. Shanbhag, “On the characterization of solution sets of smooth and nonsmooth convex stochastic Nash games,” _SIAM Journal on Optimization_ , vol. 21, no. 3, pp. 1168–1199, 2011.
  * [3] M. Staudigl and P. Mertikopoulos, “Convergent noisy forward-backward-forward algorithms in non-monotone variational inequalities,” _IFAC-PapersOnLine_ , vol. 52, no. 3, pp. 120–125, 2019.
  * [4] C.-K. Yu, M. Van Der Schaar, and A. H. Sayed, “Distributed learning for stochastic generalized Nash equilibrium problems,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 15, pp. 3893–3908, 2017.
  * [5] M. Ye and G. Hu, “Game design and analysis for price-based demand response: An aggregate game approach,” _IEEE transactions on cybernetics_ , vol. 47, no. 3, pp. 720–730, 2016.
  * [6] A. Kannan, U. V. Shanbhag, and H. M. Kim, “Addressing supply-side risk in uncertain power markets: stochastic nash models, scalable algorithms and error analysis,” _Optim. Methods Softw._ , vol. 28, no. 5, pp. 1095–1138, 2013. [Online]. Available: https://doi.org/10.1080/10556788.2012.676756
  * [7] D. Watling, “User equilibrium traffic network assignment with stochastic travel times and late arrival penalty,” _European Journal of Operational Research_ , vol. 175, no. 3, pp. 1539–1556, 2006.
  * [8] R. Henrion and W. Römisch, “On m-stationary points for a stochastic equilibrium problem under equilibrium constraints in electricity spot market modeling,” _Applications of Mathematics_ , vol. 52, no. 6, pp. 473–494, 2007\.
  * [9] I. Abada, S. Gabriel, V. Briat, and O. Massol, “A generalized Nash–Cournot model for the Northwestern European natural gas markets with a fuel substitution demand function: The gammes model,” _Networks and Spatial Economics_ , vol. 13, no. 1, pp. 1–42, 2013.
  * [10] P. Yi and L. Pavel, “An operator splitting approach for distributed generalized Nash equilibria computation,” _Automatica_ , vol. 102, pp. 111–121, 2019.
  * [11] F. Alvarez and H. Attouch, “An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping,” _Set-Valued Analysis_ , vol. 9, no. 1, pp. 3–11, 2001.
  * [12] H. Attouch and A. Cabot, “Convergence of a relaxed inertial forward–backward algorithm for structured monotone inclusions,” _Applied Mathematics & Optimization_, vol. 80, no. 3, pp. 547–598, 2019.
  * [13] ——, “Convergence of a relaxed inertial proximal algorithm for maximally monotone operators,” _Mathematical Programming_ , vol. 184, no. 1, pp. 243–287, 2020.
  * [14] D. A. Lorenz and T. Pock, “An inertial forward-backward algorithm for monotone inclusions,” _Journal of Mathematical Imaging and Vision_ , vol. 51, no. 2, pp. 311–325, 2015.
  * [15] Y. Nesterov, “A method of solving a convex programming problem with convergence rate $o(1/k^{2})$.” _Soviet Mathematics Doklady_ , vol. 27, no. 2, pp. 372–376, 1983.
  * [16] P. Bianchi, “Ergodic convergence of a stochastic proximal point algorithm,” _SIAM Journal on Optimization_ , vol. 26, no. 4, pp. 2235–2260, 2016.
  * [17] L. Rosasco, S. Villa, and B. C. Vũ, “A stochastic inertial forward–backward splitting algorithm for multivariate monotone inclusions,” _Optimization_ , vol. 65, no. 6, pp. 1293–1314, 2016.
  * [18] ——, “Stochastic forward–backward splitting for monotone inclusions,” _Journal of Optimization Theory and Applications_ , vol. 169, no. 2, pp. 388–406, 2016.
  * [19] S. Cui and U. V. Shanbhag, “Variance-reduced proximal and splitting schemes for monotone stochastic generalized equations,” _arXiv preprint arXiv:2008.11348_ , 2020.
  * [20] A. Kannan and U. V. Shanbhag, “Optimal stochastic extragradient schemes for pseudomonotone stochastic variational inequality problems and their variants,” _Computational Optimization and Applications_ , vol. 74, no. 3, pp. 779–820, 2019.
  * [21] R. I. Bot, P. Mertikopoulos, M. Staudigl, and P. T. Vuong, “Mini-batch forward-backward-forward methods for solving stochastic variational inequalities,” _Forthcoming: Stochastic Systems_ , 2021.
  * [22] P. Tseng, “A modified forward-backward splitting method for maximal monotone mappings,” _SIAM Journal on Control and Optimization_ , vol. 38, no. 2, pp. 431–446, 2018/09/13 2000. [Online]. Available: https://doi.org/10.1137/S0363012998338806
  * [23] H. H. Bauschke and P. L. Combettes, _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_. Springer - CMS Books in Mathematics, 2016.
  * [24] B. Franci, S. Grammatico, and M. Staudigl, “Distributed forward-backward (half) forward algorithms for generalized Nash equilibrium seeking,” in _European Control Conference (ECC), St. Petersburg, Russia,_ , 2020.
  * [25] A. A. Kulkarni and U. V. Shanbhag, “On the variational equilibrium as a refinement of the generalized nash equilibrium,” _Automatica_ , vol. 48, no. 1, pp. 45–55, 2012. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0005109811004821
  * [26] B. Franci and S. Grammatico, “A distributed forward-backward algorithm for stochastic generalized nash equilibrium seeking,” _IEEE Transactions on Automatic Control_ , 2020.
  * [27] A. Iusem, A. Jofré, R. I. Oliveira, and P. Thompson, “Extragradient method with variance reduction for stochastic variational inequalities,” _SIAM Journal on Optimization_ , vol. 27, no. 2, pp. 686–724
  * [28] F. Facchinei and C. Kanzow, “Generalized nash equilibrium problems,” _4or_ , vol. 5, no. 3, pp. 173–210, 2007.
  * [29] F. Facchinei, A. Fischer, and V. Piccialli, “On generalized Nash games and variational inequalities,” _Operations Research Letters_ , vol. 35, no. 2, pp. 159–164, 2007.
  * [30] F. Facchinei and J.-s. Pang, _Finite-Dimensional Variational Inequalities and Complementarity Problems - Volume I and Volume II_. Springer Series in Operations Research, 2003.
  * [31] A. Auslender and M. Teboulle, “Lagrangian duality and related multiplier methods for variational inequality problems,” _SIAM Journal on Optimization_ , vol. 10, no. 4, pp. 1097–1115, 2000.
  * [32] C. Godsil and G. F. Royle, _Algebraic graph theory_. Springer Science & Business Media, 2013, vol. 207.
  * [33] R. H. Byrd, G. M. Chin, J. Nocedal, and Y. Wu, “Sample size selection in optimization methods for machine learning,” _Mathematical Programming_ , vol. 134, no. 1, pp. 127–155, 2012. [Online]. Available: https://doi.org/10.1007/s10107-012-0572-5
  * [34] J. Lei and U. V. Shanbhag, “Distributed variable sample-size gradient-response and best-response schemes for stochastic nash equilibrium problems over graphs,” _arXiv preprint arXiv:1811.11246_ , 2018.
  * [35] R. Bot, P. Mertikopoulos, M. Staudigl, and P. Vuong, “Mini-batch forward-backward-forward methods for solving stochastic variational inequalities,” _Stochastic Systems_ , 2020.
  * [36] A. Kannan and U. V. Shanbhag, “Distributed computation of equilibria in monotone nash games via iterative regularization techniques,” _SIAM Journal on Optimization_ , vol. 22, no. 4, pp. 1177–1205, 2012.
  * [37] B. T. Polyak, _Introduction to Optimization_. Optimization Software, 1987.