# Optimal selection and tracking of generalized Nash equilibria in monotone
games

Emilio Benenati, Wicak Ananduta, and Sergio Grammatico The authors are with
the Delft Center of Systems and Control (DCSC), TU Delft, the Netherlands.
E-mail addresses: {e.benenati, w.ananduta<EMAIL_ADDRESS>This work
was partially supported by the ERC under research project COSMOS (802348).

###### Abstract

A fundamental open problem in monotone game theory is the computation of a
specific generalized Nash equilibrium (GNE) among all the available ones, e.g.
the optimal equilibrium with respect to a system-level objective. The existing
GNE seeking algorithms have in fact convergence guarantees toward an
arbitrary, possibly inefficient, equilibrium. In this paper, we solve this
open problem by leveraging results from fixed-point selection theory and in
turn derive distributed algorithms for the computation of an optimal GNE in
monotone games. We then extend the technical results to the time-varying
setting and propose an algorithm that tracks the sequence of optimal
equilibria up to an asymptotic error, whose bound depends on the local
computational capabilities of the agents.

###### Index Terms:

Multi-agent systems, Nash equilibrium seeking, Optimization

##  I Introduction

#### Motivation

Numerous engineering systems of recent interest, such as smart electrical
grids [1, 2], traffic control systems [3], and wireless communication systems
[4, 5, 6] can be modelled as a _generalized game_ , that is, a system of
multiple agents aiming at optimizing their individual, inter-dependent
objectives, while satisfying some common constraints. A typical operating
point for these systems is the Generalized Nash Equilibrium (GNE), where no
agent can unilaterally improve their objective function [7].

The recent literature has witnessed the development of theory and algorithms
for computing a variational GNE (v-GNE) [8, 7, 9], which exhibits desirable
properties of fairness and stability. Semi-decentralized GNE seeking
algorithms, where a reliable central coordinator gathers and broadcasts
aggregate information, have been proposed for strongly monotone [10, 11] and
merely monotone games [12, 13, 14]. A breakthrough idea in [15], later
generalized for non-strongly monotone games [16, 17, 18], enables a
distributed computation of GNEs by exploiting a suitable consensus protocol
[19], thus requiring a peer-to-peer information exchange.

Existing results present, however, two fundamental shortcomings that might
limit their practical application. First, unless strong assumptions are
considered (namely, strong monotonicity of the pseudogradient), a game may
have infinitely many v-GNEs _and the vast majority of the existing algorithms
provide no characterization of the equilibrium computed._ For instance, a Nash
equilibrium can be arbitrarily inefficient with respect to system-level
efficiency metrics (e.g., overall social cost) [20]. Such uncertainty on the
obtained equilibrium is often unacceptable. A notable exception is the
Tikhonov regularization algorithm [12], which guarantees convergence to the
minimum-norm solution. In addition, the method in [21, 22] seeks the (not
necessarily variational) GNE closest to a desired strategy via a double-layer
algorithm. Second, decision-making agents often operate in a time-dependent
environment and, due to the limited computation capabilities and to the time
required to exchange information, it can be impossible to ensure a time-scale
separation between the environment and the algorithm dynamics. This results in
non-constant objectives and constraints between the discrete-time algorithmic
iterations, as discussed in [23], and the references therein, for the
particular case of optimization problems. Only few works, e.g., [24, 25],
consider this setting in the case of game equilibrium problems and only with a
strong monotonicity assumption on the game pseudogradient mapping.

#### Optimal equilibrium selection and tracking

We can formulate the first issue, identified in the seminal work [7, Sect. 6],
as an _optimal GNE selection_ problem, that is, the problem of computing a GNE
of a game (among the potentially infinitely many) that satisfies a selection
criterion. This criterion characterizes the desired equilibrium and can be
formalized as a system-level _selection function_ to be optimized over the set
of GNEs. For example, the system-level objective of an electricity market can
be to minimize the deviation from an efficient operating set-point [23]; for
multiple autonomous vehicles, it can be to minimize the overall travel time of
the network. Meanwhile, the second issue can be cast as an _optimal GNE
tracking_ problem, i.e., the problem of tracking the sequence of optimal GNEs
of a time-varying game, with finite computation time and limited information
on the future instances of the game available. As the GNE set is in general
not a singleton, the tracking objective should be again chosen by means of a
(time-varying) selection function. These problems, although of high practical
interest, have never been addressed in the literature.

Under mild assumptions on the selection function, the optimal GNE selection
problem in a monotone game is a special case of a Variational Inequality
($\operatorname{VI}$) [26] defined over the set of v-GNEs. On the other hand,
as shown in [17, 13, 14], operator splitting techniques [27] can be leveraged
to characterize v-GNEs as the zeros of a monotone operator and, in turn, as
the fixed-point set of a suitable operator. Therefore, here we can cast the
problem as that of fixed-point selection [28]. In the literature, e.g., [28,
29, 30], the latter can be solved by the Hybrid Steepest Descent Method
(HSDM), whose iterations depend on the fixed-point operator (whose definition
depends on the primitives of the game) and the monotone operator that defines
the VI, namely the gradient of the selection function in our setting.

#### Contributions

In the first part of the paper (Sections III and IV), we propose the first
distributed algorithms for solving the optimal GNE selection problem. Our
method employs the Forward-Backward-Forward (FBF) operator [17] combined with
the HSDM. We show that the proposed algorithm guarantees convergence to the
_optimal_ v-GNE set in monotone games. Moreover, for a special class of
monotone games, namely cocoercive games with affine coupling constraints, we
also show that the preconditioned Forward-Backward (pFB) [13] can be paired
with the HSDM to derive optimal GNE selection algorithms. Technically, our
contribution is to show that these operators fulfill special properties that
guarantee the convergence of the HSDM toward the solution set of the
corresponding fixed-point selection $\operatorname{VI}$. Compared to the
methods in [12, 21, 22], our proposed algorithms significantly generalize the
class of selection functions and, being single-layer, they provide a
considerable advantage in computational and communication burden compared to
[21, 22].

In the second part of the paper (Section V), we formalize the optimal GNE
tracking problem as a time-varying fixed-point selection problem. Thus, as a
solution framework, we propose the _restarted HSDM_ , which adapts its
operators when the problem changes. In line with the results in the time-
varying optimization literature [31, 32], we show convergence up to a tracking
error which depends on the problem data and can be controlled by a suitable
tuning of the algorithm parameters. Similarly to the equilibrium selection
problem, the restarted HSDM works with the aforementioned fixed-point
operators to solve the optimal GNE tracking problem for the corresponding
classes of monotone games.

#### Paper organization

In Section II, we survey the required mathematical background and present a
generalization of the class of operators that comply with the conditions for
applying the HSDM. In Section III, we formalize the optimal GNE selection
problem and we explain our FBF-based algorithm for general monotone games,
while Section IV discusses the pFB-based algorithm for cocoercive games. In
Section V, we formalize the optimal GNE tracking problem and we present the
performance properties of the restarted HSDM algorithm. Finally, Section VI
illustrates the advantages of our methods on a peer-to-peer electricity market
case study.

##  II Mathematical preliminaries

##### Notation

The set of real numbers is denoted by $\mathbb{R}$. The vector of all $1$ (or
$0$) with dimension $n$ are denoted by $\mathbf{1}_{n}$ ($\bm{0}_{n}$). We
omit the subscript when the dimension is clear from the context. The operator
$\operatorname{col}(\cdot)$ stacks the arguments column-wise. For a group of
vectors $x_{i}$, $i\in\mathcal{I}=\\{1,2,\dots,N\\}$, we use the bold symbol
to denote their column concatenation, i.e.,
$\bm{x}:=\operatorname{col}(\\{x_{i}\\}_{i\in\mathcal{I}})$. The cardinality
of a set is denoted by $|\cdot|$. The operator $\langle x,y\rangle$ denotes
the inner product. We denote by $\|\cdot\|$ the Euclidean norm and by
$\|\cdot\|_{p}$ the $p$-norm. Let $P\succ 0$ be symmetric. For
$x,y\in\mathbb{R}^{n}$, $\langle x,y\rangle_{P}=\langle x,Py\rangle$ denotes
the $P$-weighted Euclidean inner product. The graph of an operator
$A\colon\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is denoted by
$\operatorname{gph}(A)$. $\operatorname{zer}(A)$ defines the set of zeros of
operator $A$, i.e., $\operatorname{zer}(A):=\\{x\in\operatorname{dom}(A)\mid
0\in A(x)\\}$ whereas $\operatorname{fix}(A)$ defines the set of fixed points
of $A\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$, i.e.,
$\operatorname{fix}(A):=\\{x\in\operatorname{dom}(A)\mid A(x)=x\\}$.

##### Convex functions

A continuously differentiable function $f\colon\mathbb{R}^{n}\to\mathbb{R}$ is
$\sigma$-strongly convex with respect to a $p$-norm, with $\sigma>0$, if, for
all $x,x^{\prime}\in\operatorname{dom}f$, $f(x^{\prime})\geq
f(x)+\langle\nabla
f(x),x^{\prime}-x\rangle+\frac{\sigma}{2}\|x^{\prime}-x\|_{p}^{2}.$
Additionally, $f$ is convex if the previous inequality hold for $\sigma=0$.
The projection onto a closed convex set $C$ is denoted by
$\operatorname{proj}_{C}^{\Psi}(x)=\operatorname{argmin}_{z\in
C}\|x-z\|_{\Psi}$, where $\Psi\succ 0$. For a convex function $f$ with
subdifferential $\partial f$ and $\Psi\succ 0$, the operator
$\mathrm{prox}_{\partial
f}^{\Psi}(x):=\operatorname{argmin}_{z}f(z)+\tfrac{1}{2}\|z-x\|_{\Psi}^{2}$
[27, Def. 12.23]. For example, for the indicator function of a closed convex
set $C$, $\iota_{C}$, where $\partial\iota_{C}=\mathrm{N}_{C}$ being the
normal cone operator [27, Ex. 1.25, 16.13],
$\mathrm{prox}_{\iota_{C}}^{\Psi}(x)=\operatorname{proj}_{C}^{\Psi}(x)$ [27,
Ex. 12.25].

##### Operator theoretic definitions

An operator $A\colon\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is
_monotone_ if, for any $(x,y)\in\operatorname{gph}(A)$ and
$(x^{\prime},y^{\prime})\in\operatorname{gph}(A)$, $\langle
y-y^{\prime},x-x^{\prime}\rangle\geq 0$ [27, Def. 20.1], and $\beta$-_strongly
monotone_ if $A-\beta\mathrm{Id}$, where $\mathrm{Id}$ is the identity
operator, is monotone. Let $C$ be a nonempty subset of $\mathbb{R}^{n}$. A
single-valued operator $\mathcal{T}\colon C\to\mathbb{R}^{n}$ is _Lipschitz
continuous_ if there exists a constant $L>0$, such that, for all
$x,x^{\prime}\in\mathbb{R}^{n}$,
$\|\mathcal{T}(x)-\mathcal{T}(x^{\prime})\|\leq L\|x-x^{\prime}\|$ [27, Def.
1.47]. In particular, the operator $\mathcal{T}$ is (i) _nonexpansive_ if
$L=1$, (ii) _attracting nonexpansive_ if $\mathcal{T}$ is nonexpansive with
$\operatorname{fix}(\mathcal{T})\neq\varnothing$ and
$\|\mathcal{T}(x)-z\|<\|x-z\|$, for all $z\in\operatorname{fix}(\mathcal{T})$
and all $x\notin\operatorname{fix}(\mathcal{T})$; and (iii) _quasi-
nonexpansive_ if $\operatorname{fix}(\mathcal{T})\neq\varnothing$ and
$\|\mathcal{T}(x)-z\|\leq\|x-z\|$, for all
$z\in\operatorname{fix}(\mathcal{T})$ and $x\in\mathbb{R}^{n}$. Moreover,
$\mathcal{T}$ is $\alpha$-averaged nonexpansive, for $\alpha\in(0,1)$, if
there exists a nonexpansive operator $\mathcal{R}\colon C\to\mathbb{R}^{n}$
such that $\mathcal{T}=(1-\alpha)\mathrm{Id}+\alpha\mathcal{R}$. If
$\mathcal{T}$ is averaged nonexpansive with
$\operatorname{fix}(\mathcal{T})\neq\varnothing$, then $\mathcal{T}$ is
attracting [28, Sec. 2.A]. Additionally $\mathcal{T}$ is $\beta$-_cocoercive_
if
$\langle\mathcal{T}(x)-\mathcal{T}(y),x-y\rangle\geq\beta\|\mathcal{T}(x)-\mathcal{T}(y)\|$.

Now, let $C$ be a non-empty, closed, and convex subset of $\mathbb{R}^{n}$,
$\mathcal{T}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ be quasi-nonexpansive under
the $\Psi$-induced norm $\|\cdot\|_{\Psi}$ for some positive definite matrix
$\Psi$, i.e., $\|\mathcal{T}(x)-z\|_{\Psi}\leq\|x-z\|_{\Psi}$, for all
$z\in\operatorname{fix}(\mathcal{T})\neq\varnothing$ and $x\in\mathbb{R}^{n}$.
We define the distance of a point $x\in\mathbb{R}^{n}$ to $C$ by
$\mathrm{dist}_{\Psi}(x,C):=\inf_{z\in C}\|x-z\|_{\Psi},$. For $r\geq 0$, we
define the set

$\displaystyle C_{\geq
r}^{\Psi}:=\\{x\in\mathbb{R}^{n}\mid\mathrm{dist}_{\Psi}(x,C)\geq r\\}.$ (1)

Furthermore, let us define the function

$D_{\Psi}(r):=\begin{cases}{\inf}\ \
\left(\mathrm{dist}_{\Psi}(x,\operatorname{fix}(\mathcal{T}))-\mathrm{dist}_{\Psi}(\mathcal{T}(x),\operatorname{fix}(\mathcal{T}))\right),\\\
\operatorname{s.t.}\ {x\in(\operatorname{fix}(\mathcal{T}))_{\geq
r}^{\Psi}\cap C},\ \text{if }(\operatorname{fix}(\mathcal{T}))_{\geq
r}^{\Psi}\cap C\neq\varnothing,\\\ +\infty,\quad\text{otherwise.}\end{cases}$
(2)

For $\Psi=I$, we omit the subscript of $D$. We sometimes refer to $D_{\Psi}$
as the _shrinkage function_ under the norm $\|\cdot\|_{\Psi}$.

The function $D_{\Psi}$ has the properties stated next in Proposition 1 (see
[33, Prop. 2.6] for the case $\Psi=I$).

######  Proposition 1.

Let $\Psi$ be positive definite. For the function $D_{\Psi}$ defined in (2),
it holds that:

  1. (i)

$D_{\Psi}$ is positive semidefinite and non-decreasing;

  2. (ii)

$D_{\Psi}(\mathrm{dist}(x,\operatorname{fix}(\mathcal{T})))\leq\|x-\mathcal{T}(x)\|_{\Psi}$
for all $x\in C$.

######  Definition 1 (Quasi-shrinking [28]).

A quasi-nonexpansive operator
$\mathcal{T}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is _quasi-shrinking_ on a
non-empty, closed, and convex set $C\subseteq\mathbb{R}^{n}$ if
$\operatorname{fix}(\mathcal{T})\cap C\neq\varnothing$ and
$D(r)=0\Leftrightarrow r=0$, where $D(r)$ is defined as in (2).

######  Remark 1.

Suppose that a quasi-nonexpansive operator $\mathcal{T}$ is quasi-shrinking on
C, i.e., $D(r)=0\Leftrightarrow r=0$. Then, it also holds that
$D_{\Psi}(r)=0\Leftrightarrow r=0$, for any $\Psi\succ 0$.

######  Example 1.

The Euclidean projection onto $C$, $\operatorname{proj}_{C}$ is quasi-
shrinking and its shrinkage function (defined in (2)) is

$D(r)=\inf_{\\{u|\mathrm{dist}(u,C)\geq
r\\}}\mathrm{dist}(u,C)-\underbrace{\mathrm{dist}(\operatorname{proj}_{C}(u),C)}_{=0}=r.$

Finally, we identify a class of quasi-shrinking operators, as formally stated
in Lemma 1, which generalizes the result in [33, Prop. 2.11] and is useful for
our analysis.

######  Definition 2 (Demiclosed operator [27, Def. 4.26]).

Let $C\subseteq\mathbb{R}^{{n}}$ be a closed set. An operator
$\mathcal{T}\colon C\to\mathbb{R}^{{n}}$ is demiclosed at
$u\in\mathbb{R}^{{n}}$ if $\mathcal{T}(\bm{\omega}^{\infty})=u$, for any
sequence $(\bm{\omega}_{k})_{k\in\mathbb{N}}\in C$ such that
$\lim_{k\xrightarrow[]{}\infty}\bm{\omega}_{k}=\bm{\omega}^{\infty}$ and
$\lim_{k\xrightarrow[]{}\infty}\mathcal{T}(\bm{\omega}_{k})=u$.

######  Lemma 1.

Let $\mathcal{T}$ be quasi-nonexpansive, with
$\operatorname{fix}(\mathcal{T})\neq\varnothing$. Let ${\mathcal{T}_{2}}$ be
an operator such that ${\mathrm{Id}}-\mathcal{T}_{2}$ is demiclosed at $0$ and
such that
$\operatorname{fix}({\mathcal{T}_{2}})\subseteq\operatorname{fix}(\mathcal{T})$.
Assume that for any $\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T})$,

$\|\mathcal{T}(\bm{\omega})-\bm{\omega}^{\star}\|^{2}_{\Psi}\leq\|\bm{\omega}-\bm{\omega}^{\star}\|^{2}_{\Psi}-\gamma\|\bm{\omega}-\mathcal{T}_{2}(\bm{\omega})\|^{2}_{\Psi},$
(3)

for some $\gamma>0$ and $\Psi\succ 0$. Then, $\mathcal{T}$ is quasi-shrinking
on any compact convex set $C$ such that
$C\bigcap\operatorname{fix}(\mathcal{T})\neq\varnothing$.

###### Proof.

See Appendix A. ∎

##  III Optimal selection of generalized Nash equilibria

###  III-A Generalized Nash equilibrium problem

Let us consider $N$ agents, denoted by the set
$\mathcal{I}:=\\{1,2,\dots,N\\}$, with inter-dependent optimization problems:

$\displaystyle\underset{x_{i}\in\mathcal{X}_{i}}{\operatorname{min}}\quad$
$\displaystyle J_{i}(\bm{x}):={\ell_{i}(x_{i})+}f_{i}(\bm{x})$ (4a)
$\displaystyle\operatorname{s.t.}\quad$
$\displaystyle\sum_{j\in\mathcal{I}}g_{j}(x_{j})\leq 0,$ (4b)

where $x_{i}\in\mathbb{R}^{n_{i}}$ is the decision variable of agent $i$
whereas
$\bm{x}:=\operatorname{col}((x_{i})_{i\in\mathcal{I}})\in\mathbb{R}^{n}$ is a
concatenated vector of the decision variables of all agents. Let us use
$\bm{x}_{-i}=\operatorname{col}(\\{x_{j}\\}_{j\in\mathcal{I}\backslash\\{i\\}})$
to denote the concatenated decision variables of all agents except agent $i$.
Let $\mathcal{X}_{i}\subseteq\mathbb{R}^{n_{i}}$ denote the local feasible set
of $x_{i}$ and $J_{i}\colon\mathbb{R}^{n}\to\mathbb{R}$ denote the cost
function of agent $i$ that depends on the decision variables of other agents.
Moreover, (4b) represents a separable coupling constraint with the function
$g_{j}\colon\mathbb{R}^{n_{j}}\to\mathbb{R}^{m}$ associated with agent $j$.

We denote the collective feasible set of the game in (4) by

${\bm{\Omega}:=\prod_{i\in\mathcal{I}}\mathcal{X}_{i}\cap\Big{\\{}\bm{x}\mid\sum_{j\in\mathcal{I}}g_{j}(x_{j})\leq
0\Big{\\}}.}$ (5)

Here, we look for equilibrium solutions to (4) where no agent has the
incentive to unilaterally deviate, namely, GNE:

######  Definition 3.

A set of strategies
$\bm{x}^{\ast}:=\operatorname{col}((x_{i}^{\ast})_{i\in\mathcal{I}})$ is a
generalized Nash equilibrium (GNE) of the game in (4) if
$\bm{x}^{\ast}\in{\bm{\Omega}}$ and, for each $i\in\mathcal{I}$,

$J_{i}(\bm{x}^{\ast})\leq J_{i}(x_{i},\bm{x}_{-i}^{\ast}),$ (6)

for any $x_{i}\in\mathcal{X}_{i}\cap\\{y\mid
g_{i}(y)\leq-\sum_{j\in\mathcal{I}\backslash\\{i\\}}g_{j}(x_{j}^{\ast})\\}$.
$\hfill\Box$

Furthermore, we focus on the class of jointly convex GNEP and hence, consider
the following assumptions on Problem (4) [13, Assms 1–2]. We note that [15,
16, 17, 18, 14] consider the case of affine constraint functions.

######  Assumption 1.

In (4), for each $i\in\mathcal{I}$, the functions $f_{i}(\cdot,\bm{x}_{-i})$,
for any $\bm{x}_{-i}$, and $g_{i}(\cdot)$ are component-wise convex and
continuously differentiable; $\ell_{i}$ is convex and lower semicontinuous.
For each $i\in\mathcal{I}$, the set $\mathcal{X}_{i}$ is nonempty, compact,
and convex. The global feasible set $\bm{\Omega}$ defined in (5) is non-empty
and satisfies Slater’s constraint qualification [27, Eq. (27.50)].

######  Assumption 2.

The mapping

$F(\bm{x}):=\operatorname{col}(({\nabla_{x_{i}}f_{i}}(\bm{x}))_{i\in\mathcal{I}}),$
(7)

with $f_{i}$ as in (4a), is monotone.

As in [13, 15, 16, 17, 18, 14], we can formulate the problem of finding a GNE
of the game in (4) as that of a monotone inclusion. To this end, we introduce
the dual variable $\lambda_{i}\in\mathbb{R}^{m}$, for each $i\in\mathcal{I}$,
to be associated with the coupling constraint (4b). Furthermore, we focus on a
subset of GNEs, namely variational GNE (v-GNE), indicated by equal optimal
dual variables, $\lambda_{i}^{\ast}=\lambda^{\ast}$, for all
$i\in\mathcal{I}$. As discussed in [7, 9], a v-GNE enjoys several desirable
properties, such as fairness and larger social stability than non-variational
ones. Under Assumptions 1–2, the set of v-GNEs of the game in (4) is non-empty
[34, Prop. 12.11]. The Karush-Kuhn-Tucker (KKT) optimality conditions of a
v-GNE of the game in (4), denoted by $\bm{x}^{\ast}$, are:

$\displaystyle\mathbf{0}\in\mathrm{N}_{\mathcal{X}_{i}}(x_{i}^{\ast})+\partial_{x_{i}}J_{i}(\bm{x}^{\ast})+\langle\nabla
g_{i}(x_{i}^{\ast}),{\lambda}^{\ast}\rangle,$ (8a)
$\displaystyle\mathbf{0}\in\mathrm{N}_{\mathbb{R}^{m}_{\geq
0}}({\lambda}^{\ast})-\sum_{j\in\mathcal{I}}g_{j}(x_{j}^{\ast}).$ (8b)

To obtain a v-GNE via a fully distributed algorithm, we incorporate a
consensus scheme on the dual variables. In the full information case, one
typically assumes that there exists a communication network over which the
agents exchange information to update their dual variables. Let us represent
this communication network as an undirected graph
$\mathcal{G}^{\lambda}=(\mathcal{I},\mathcal{E}^{\lambda})$ and assume that
$\mathcal{G}^{\lambda}$ is connected. Furthermore, we denote the Laplacian of
$\mathcal{G}^{\lambda}$ by $L$ and the neighbors of agent $i$ in
$\mathcal{G}^{\lambda}$ by $\mathcal{N}_{i}^{\lambda}$, i.e.,
$\mathcal{N}_{i}^{\lambda}:=\\{j\in\mathcal{I}\mid(i,j)\in\mathcal{E}^{\lambda}\\}$.
Additionally, let $\mathcal{N}_{i}^{J}$ denote the set of agents whose
decision variable $x_{j}$ influences the cost function $J_{i}$. For
simplicity, we assume that
$\mathcal{N}_{i}^{J}\subseteq\mathcal{N}_{i}^{\lambda}$.

Now, let us denote $\nu_{i}\in\mathbb{R}^{m}$ as the consensus variable of
agent $i$, and
$\bm{\omega}=(\bm{x},\bm{\lambda},\bm{\nu})\in\mathbb{R}^{n_{\omega}}$, where
$\bm{\lambda}=\operatorname{col}(\\{\lambda_{i}\\}_{i\in\mathcal{I}})$,
$\bm{\nu}=\operatorname{col}(\\{\nu_{i}\\}_{i\in\mathcal{I}})$, and
$n_{\omega}=n+2Nm$. Then, we can define the operators
$\mathcal{A}\colon{\bm{\mathcal{X}}}\times\mathbb{R}^{Nm}_{\geq
0}\times\mathbb{R}^{Nm}\to\mathbb{R}^{n_{\omega}}$,
$\mathcal{B}\colon\mathbb{R}^{n_{\omega}}\to\mathbb{R}^{n_{\omega}}$, and
$\mathcal{C}\colon\mathbb{R}^{n}\times\mathbb{R}_{\geq
0}^{Nm}\times\mathbb{R}^{Nm}\to\mathbb{R}^{n_{\omega}}$, as follows:

$\displaystyle\mathcal{A}(\bm{\omega})$
$\displaystyle:=\prod_{i\in\mathcal{I}}(\mathrm{N}_{\mathcal{X}_{i}}{+\partial\ell_{i}})(x_{i})\times\mathrm{N}_{\mathbb{R}^{Nm}_{\geq
0}}(\bm{\lambda})\times\\{\mathbf{0}_{Nm}\\},$ (9)
$\displaystyle\mathcal{B}(\bm{\omega})$
$\displaystyle:=\operatorname{col}({F(\bm{x})},(L\otimes
I_{m})\bm{\lambda},\mathbf{0}_{Nm}),$ (10)
$\displaystyle\mathcal{C}(\bm{\omega})$
$\displaystyle:=\operatorname{col}((\langle\nabla
g_{i}(x_{i}),\lambda_{i}\rangle)_{i\in\mathcal{I}},-(g_{i}(x_{i}))_{i\in\mathcal{I}}-(L\otimes
I_{m})\bm{\nu},$ $\displaystyle\qquad\quad(L\otimes I_{m})\bm{\lambda}).$ (11)

In turn, we can translate the GNEP in (4) as a monotone inclusion problem,
i.e.,

$\text{find }\bm{\omega}\text{ such that
}{\bm{\omega}\in\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C}).}$
(12)

Similarly to [15, Thm. 2], we can show that for any $\bm{\omega}$ such that
(12) holds, we obtain the pair $(\bm{x},\lambda)$ that satisfies the KKT
conditions in (8) if Assumptions 1-2 hold (see Appendix B for details).
Furthermore, due to the maximal monotonicity of
$(\mathcal{A}+\mathcal{B}+\mathcal{C})$ (Lemma 6 in Appendix B),
$\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$ is convex [27, Prop.
23.39]. Additionally, since the set of v-GNE of the game is bounded as it is a
subset of $\bm{\mathcal{X}}$, the set of solutions of the inclusion in (8) and
the set $\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$ are bounded
[35, Prop. 3.3].

###  III-B Optimal equilibrium selection problem

The inclusion problem in (12) may have multiple solutions. In this section, we
want to find an equilibrium solution that minimizes a selection function,
denoted by $\phi\colon\mathbb{R}^{n_{\omega}}\to\mathbb{R}$, i.e.,

$\left\\{\begin{array}[]{rl}\underset{\bm{\omega}}{\operatorname{argmin}}&\phi(\bm{\omega})\\\
\operatorname{s.t.}&\bm{\omega}\in\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C}).\end{array}\right.$
(13)

For example, we can consider the selection function

${\phi_{\mathrm{ex}}}(\bm{\omega})=\|Q\bm{\omega}-\bm{\omega}^{\mathrm{ref}}\|,$
(14)

for some $Q\succcurlyeq 0$. When $Q=I$ and $\omega^{\mathrm{ref}}=\mathbf{0}$,
the objective is to find a minimum norm v-GNE. The vector
$\omega^{\mathrm{ref}}$ can be any desired strategy of the agents, and thus
the objective is to find the v-GNE closest to this strategy, as discussed in
[21, 22]. In some engineering applications, such as electrical networks, (14)
can represent system level objectives (see Section VI). In the remainder of
the paper, we consider the following technical assumption on the selection
function, which, together with the convexity of
$\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$, guarantees that the
optimization problem in (13) is convex.

######  Assumption 3.

The function $\phi$ in (13) is continuously differentiable, convex, and has
$L_{\phi}$-Lipschitz continuous gradient.

As a first step towards computing an optimal variational GNE, we leverage
existing results to derive operators $\mathcal{T}$ with the property that

$\bm{\omega}\in\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})\Leftrightarrow\bm{\omega}\in\operatorname{fix}(\mathcal{T}),$
(15)

and such that the Banach-Picard iteration of $\mathcal{T}$ [27, Sect. 5.2]
guarantees convergence to a solution of the inclusion in (12). For instance,
for cocoercive generalized games, a preconditioned forward-backward (pFB)
operator presents the desired characteristics[15], whereas the forward-
reflected-backward (FRB) operator [36] or the forward-backward-forward (FBF)
operator [37] meets these requirements even for general monotone games.
Furthermore, we require that the operator $\mathcal{T}$ in (15) can be
evaluated in a distributed manner. Therefore, by (15) and Assumption 3, the
optimal equilibrium selection problem in (13) can be cast as a fixed-point
selection $\operatorname{VI}$:

$\text{find }\bm{\omega}^{\star}\text{ s.t.
}\inf_{\bm{\omega}\in\operatorname{fix}(\mathcal{T})}\left\langle{\bm{\omega}-\bm{\omega}^{\star}},{{\nabla\phi}(\bm{\omega}^{\star})}\right\rangle\geq
0.$ (16)

###  III-C Distributed optimal equilibrium selection algorithm

With the aim of solving the $\operatorname{VI}$ in (16), we consider a fixed-
point selection algorithm called the hybrid steepest descent method (HSDM)
[28], which is defined by the following discrete-time dynamical system or
iteration:

$\bm{\omega}^{(k+1)}=\mathcal{T}(\bm{\omega}^{(k)})-\beta^{(k)}{\nabla\phi}(\mathcal{T}(\bm{\omega}^{(k)})).$
(17)

The HSDM can solve Problem (16) when $\mathcal{T}$ is quasi-nonexpansive and
quasi-shrinking with bounded $\operatorname{fix}(\mathcal{T})$, as formally
stated next.

######  Assumption 4.

The step size of the HSDM $\beta^{(k)}$ satisfies:

  1. (i)

$\lim_{k\to\infty}\beta^{(k)}=0$, $\sum_{k\geq 1}\beta^{(k)}=\infty$;

  2. (ii)

$\sum_{k\geq 1}(\beta^{(k)})^{2}<\infty$. $\hfill\Box$

######  Remark 2.

The sequence $\beta^{(k)}={\beta_{0}}/{k^{p}}$, for any $\beta_{0}>0$ and
$p\in({1/2},1]$, satisfies Assumption 4. $\hfill\Box$

######  Assumption 5.

$\mathcal{T}$ is quasi-nonexpansive. $\hfill\Box$

######  Assumption 6.

There exists a nonempty bounded closed convex set $C$ on which $\mathcal{T}$
is quasi-shrinking. $\hfill\Box$

######  Lemma 2 (From [28, Thm. 5]).

Let Assumption 3 hold and $\Omega^{\star}$ be the set of solutions of the
$\operatorname{VI}$ in (16), with non-empty and bounded
$\operatorname{fix}(\mathcal{T})$. Suppose that $\mathcal{T}$ satisfies
Assumptions 5 and 6 and that $(\bm{\omega}^{(k)})_{k\geq 0}\subset C$. If the
step size $\beta^{(k)}$ satisfies Assumption 4.i, then the HSDM in (17)
generates a sequence $(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ such that

$\lim_{k\to\infty}\mathrm{dist}(\bm{\omega}^{(k)},\Omega^{\star})=0.$

Therefore, our main technical task is to find a suitable operator
$\mathcal{T}$ that can be evaluated in a distributed manner and that satisfies
both (15) and Assumptions 5–6, required for the convergence of the HSDM
sequence.

Under mere monotonicity of the pseudogradient mapping (Assumption 2), perhaps
the most obvious choice is the FRB splitting, which, however, is not quasi-
nonexpansive111The FRB iteration does not generate a Fejér monotone sequence
[36, Prop. 2.3], implying that it is not quasi-nonxepansive and violates
Assumption 5. (and, thus, it is not quasi-shrinking). Another viable option is
the FBF splitting method [37], which works for v-GNE seeking in monotone games
satisfying Assumptions 1–2, as shown in [13, 17]. As our first technical
result, we show that the FBF algorithm satisfies both the desired property in
(15) and Assumptions 5–6. To that end, firstly, we compactly state the FBF
operator for (12), as follows:

$\displaystyle\mathcal{T}_{\mathrm{FBF}}(\bm{\omega})$
$\displaystyle:=((\mathrm{Id}-\Psi^{-1}(\mathcal{B}+\mathcal{C}))(\mathrm{Id}+{\Psi^{-1}\mathcal{A}})^{-1}$
$\displaystyle\quad\cdot(\mathrm{Id}-\Psi^{-1}(\mathcal{B}+\mathcal{C}))+\Psi^{-1}(\mathcal{B}+\mathcal{C}))(\bm{\omega}),$
(18)

where $\Psi\succ 0$ is a diagonal positive definite matrix. The FBF requires
the forward operator, which is $(\mathcal{B}+\mathcal{C})$, to be Lipschitz
continuous. A sufficient condition for this requirement is given in Assumption
7 (see Lemma 7 in Appendix B). Under maximal monotonicity and Lipschitz
continuity, it holds that
$\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})=\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$
(see Lemma 9 in Appendix C). In addition, by denoting $L_{B}$ as the Lipschitz
constant of $\mathcal{B}+\mathcal{C}$, we define the step-size matrix $\Psi$
in Assumption 8, which guarantees the convergence of the sequence generated by
the fixed-point iteration with $\mathcal{T}_{\mathrm{FBF}}$ toward a point in
$\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$.

######  Assumption 7.

The mapping $F(\bm{x})$ in (7) is $L_{F}$-Lipschitz continuous. Furthermore,
for each $i\in\mathcal{I}$, the function $g_{i}$ in (4b) has a bounded and
$L_{g}$-Lipschitz continuous gradient.

######  Assumption 8.

It holds that $|\Psi^{-1}|\leq 1/L_{B}$, where $L_{B}>0$ is the Lipschitz
constant of $\mathcal{B}+\mathcal{C}$ and
$\Psi=\operatorname{diag}(\rho^{-1},\sigma^{-1},\tau^{-1})$, where
$\rho=\operatorname{diag}(\\{\rho_{i}I_{n_{i}}\\}_{i\in\mathcal{I}})$,
$\tau=\operatorname{diag}(\\{\tau_{i}I_{m}\\}_{i\in\mathcal{I}})$, and
$\sigma=\operatorname{diag}(\\{\sigma_{i}I_{m}\\}_{i\in\mathcal{I}})$.

We are now ready to present the distributed FBF for seeking an optimal
variational GNE based on the selection function $\phi(\bm{\omega})$ via the
HSDM as shown in Algorithm 1.

To have a convergence guarantee as stated in Lemma 2, the FBF operator must
satisfy Assumptions 5 and 6. Let us show that this is the case in the
following lemma.

######  Lemma 3.

Let Assumptions 1, 2, and 7 hold. The operator $\mathcal{T}_{\mathrm{FBF}}$ in
(18), where $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ are defined in
(9)–(11) and $\Psi$ is defined in Assumption 8, is quasi-nonexpansive and
quasi-shrinking on any compact convex set $C$ such that
$C\cap\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})\neq\varnothing$ (see
Definition 1).

###### Proof.

See Appendix C-A. ∎

Furthermore, we observe that the HSDM sequence generated by using
$\mathcal{T}_{\mathrm{FBF}}$ is bounded, as formally stated next.

######  Lemma 4.

Let Assumptions 1–4 and 7–8 hold. Then, the sequence
$(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ generated by the HSDM method in (17)
with $\mathcal{T}=\mathcal{T}_{\mathrm{FBF}}$ in (18), where $\mathcal{A}$,
$\mathcal{B}$, and $\mathcal{C}$ are defined in (9)–(11) and $\Psi$ is defined
in Assumption 8, is bounded, i.e., for any arbitrary
$\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$, it
holds that $\|\bm{\omega}^{(k)}-\bm{\omega}^{\star}\|\leq
R(\bm{\omega}^{\star}),$ for some positive finite $R(\bm{\omega}^{\star})$.

###### Proof.

See Appendix C-B. ∎

Algorithm 1 Optimal v-GNE selection via FBF and HSDM

Initialization. Set $x_{i}^{(0)}\in\mathcal{X}_{i}$,
$\lambda_{i}^{(0)}\in\mathbb{R}_{\geq 0}^{m}$, and
$\nu_{i}^{(0)}\in\mathbb{R}^{m}$, for all $i\in\mathcal{I}$.

Iteration of each agent $i\in\mathcal{I}$.

  1. 1.

Receives $x_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{J}$ and
$\lambda_{j}^{(k)},\nu_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{\lambda}$.

  2. 2.

Updates:

$\displaystyle\tilde{x}_{i}^{(k)}$
$\displaystyle={\mathrm{prox}_{\ell_{i}+\iota_{\mathcal{X}_{i}}}^{\rho_{i}}}\left(x_{i}^{(k)}-\rho_{i}\big{(}\nabla_{x_{i}}{f}_{i}(\bm{x}^{(k)})+\nabla
g_{i}(x_{i}^{(k)})^{\top}\lambda_{i}^{(k)}\big{)}\right),$
$\displaystyle\tilde{\lambda}_{i}^{(k)}$
$\displaystyle=\operatorname{proj}_{\geq
0}\Big{(}\Big{.}\lambda_{i}^{(k)}+\tau_{i}\Big{(}\Big{.}g_{i}(x_{i}^{(k)})$
$\displaystyle\qquad\quad\qquad\Big{.}\Big{.}+\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\nu_{i}^{(k)}-\nu_{j}^{(k)}-\lambda_{i}^{(k)}+\lambda_{j}^{(k)}\right)\Big{)}\Big{)},$
$\displaystyle\tilde{\nu}_{i}^{(k)}$
$\displaystyle=\nu_{i}^{(k)}-\sigma_{i}\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\lambda_{i}^{(k)}-\lambda_{j}^{(k)}\right).$

  3. 3.

Receives $\tilde{x}_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{J}$ and
$\tilde{\lambda}_{j}^{(k)},\tilde{\nu}_{j}^{(k)}$ from agent
$j\in\mathcal{N}_{i}^{\lambda}$.

  4. 4.

Updates:

$\displaystyle\accentset{\circ}{x}_{i}^{(k)}$
$\displaystyle=\tilde{x}_{j}^{(k)}-\rho_{i}\Big{(}\Big{.}\nabla_{x_{i}}f_{i}(\tilde{\bm{x}}^{(k)})-\nabla_{x_{i}}f_{i}(\bm{x}^{(k)})$
$\displaystyle\quad+\nabla
g_{i}(\tilde{x}_{i}^{(k)})^{\top}\tilde{\lambda}_{i}^{(k)}-\nabla
g_{i}(x_{i}^{(k)})^{\top}\lambda_{i}^{(k)}\Big{.}\Big{)},$
$\displaystyle\accentset{\circ}{\lambda}_{i}^{(k)}$
$\displaystyle=\tilde{\lambda}_{i}^{(k)}+\tau_{i}\Big{(}\Big{.}g_{i}(\tilde{x}_{i}^{(k)})-g_{i}(x_{i}^{(k)})$
$\displaystyle\quad+\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\Big{(}\Big{.}\tilde{\nu}_{i}^{(k)}-\nu_{i}^{(k)}-\tilde{\nu}_{j}^{(k)}+\nu_{j}^{(k)}\Big{.}\Big{)}$
$\displaystyle\quad-\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\Big{(}\Big{.}\tilde{\lambda}_{i}^{(k)}-\lambda_{i}^{(k)}-\tilde{\lambda}_{j}^{(k)}+\lambda_{j}^{(k)}\Big{.}\Big{)}\Big{.}\Big{)},$
$\displaystyle\accentset{\circ}{\nu}_{i}^{(k)}$
$\displaystyle=\tilde{\nu}_{i}^{(k)}-\sigma_{i}\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\tilde{\lambda}_{i}^{(k)}-\lambda_{i}^{(k)}-\tilde{\lambda}_{j}^{(k)}+\lambda_{j}^{(k)}\right).$

  5. 5.

Sends
$(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})$
to a coordinator and receives back
$\nabla_{\omega_{i}}\phi(\accentset{\circ}{\bm{x}}^{(k)},\accentset{\circ}{\bm{\lambda}}^{(k)},\accentset{\circ}{\bm{\nu}}^{(k)})$,
where $\omega_{i}=(x_{i},\lambda_{i},\nu_{i})$.

  6. 6.

Updates:

$\displaystyle(x_{i}^{(k+1)},\lambda_{i}^{(k+1)},\nu_{i}^{(k+1)})$
$\displaystyle=(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})-\beta^{(k)}\nabla_{\omega_{i}}\phi(\accentset{\circ}{\bm{x}}^{(k)},\accentset{\circ}{\bm{\lambda}}^{(k)},\accentset{\circ}{\bm{\nu}}^{(k)}).$
(19)

Thus, we can now show that Algorithm 1 generates a sequence that converges
toward the solution set of the optimal GNE selection problem in (16).

######  Theorem 1.

Let Assumptions 1–4 and 7–8 hold. Let $\Omega^{\star}$ be the set of solutions
to Problem (16) with $\mathcal{T}=\mathcal{T}_{\mathrm{FBF}}$ defined in (18),
where $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ are defined in (9)–(11).
Furthermore, let $(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$, where
$\bm{\omega}^{(k)}=(\bm{x}^{(k)},\bm{\lambda}^{(k)},\bm{\nu}^{(k)})$, be the
sequence generated by Algorithm 1. Then,
$\lim_{k\to\infty}\mathrm{dist}(\bm{\omega}^{(k)},\Omega^{\star})=0,$ and
$(\bm{x}^{(k)})_{k\in\mathbb{N}}$ converges to an optimal v-GNE of the game in
(4).

###### Proof.

See Appendix C-C. ∎

######  Remark 3.

A central coordinator and step 5 of Algorithm 1 are not needed if $\phi$ is a
separable function, i.e.,
$\phi(\bm{\omega})=\sum_{i\in\mathcal{I}}\phi_{i}(\omega_{i})$. In this case,
step 6 can be immediately executed by using local information
$(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})$
only, as long as each agent $i$ knows the gradient $\nabla\phi_{i}$.

##  IV Optimal equilibrium selection in cocoercive games

In this section, we discuss a special class of monotone games, namely
cocoercive games with affine coupling constraints. These games arise as a
generalization of the widely studied class of strongly monotone games
[15],[11]. Differently from the strong monotonicity assumption, however,
cocoercivity alone does not guarantee the uniqueness of the v-GNE.

######  Assumption 9 ([14, Assm. 5]).

The mapping $F$ in (7) is $\eta$-cocoercive.

######  Assumption 10 ([14, Eq. (3)]).

For each $i\in\mathcal{I}$, the function $g_{i}$ in (4b) is affine, i.e.,
$g_{i}(x_{i}):=A_{i}x_{i}-b_{i}$, for some matrix $A_{i}\in\mathbb{R}^{m\times
n_{i}}$ and vector $b_{i}\in\mathbb{R}^{m}$.

For this particular class of games, the preconditioned forward-backward (pFB)
splitting [15] can efficiently compute a variational GNE. We note that,
although [15] considers games with strongly monotone pseudogradient, the FB
splitting only requires cocoercivity of the forward operator[27, Thm. 26.14].
Compared with the FBF, the pFB has the advantages of only having one
communication round per iteration (as opposed to two) and larger step size
bounds. A numerical performance comparison is provided in [17].

Given the particular structure of the coupling constraint as stated in
Assumption 10, we can rewrite the operators in (12) as follows:

$\displaystyle\mathcal{A}(\bm{\omega})$
$\displaystyle:=\prod_{i\in\mathcal{I}}(\mathrm{N}_{\mathcal{X}_{i}}{+\partial\ell_{i}})(x_{i})\times\mathrm{N}_{\mathbb{R}^{Nm}_{\geq
0}}(\bm{\lambda})\times\\{\mathbf{0}_{Nm}\\},$ (20)
$\displaystyle\mathcal{B}(\bm{\omega})$
$\displaystyle:=\operatorname{col}({F(\bm{x})},(L\otimes
I_{m})\bm{\lambda}+\bm{b},\mathbf{0}_{Nm}),$ (21)
$\displaystyle\mathcal{C}(\bm{\omega})$
$\displaystyle:=\operatorname{col}(\bm{A}^{\top}\bm{\lambda},-\bm{A}\bm{x}-(L\otimes
I_{m})\bm{\nu},(L\otimes I_{m})\bm{\lambda}),$ (22)

where $\bm{A}=\operatorname{diag}(\\{A_{i}\\}_{i\in\mathcal{I}})$ and
$\bm{b}=\operatorname{col}(\\{b_{i}\\}_{i\in\mathcal{I}})$. Thus, the pFB
operator for the monotone inclusion in (12) based on the operators
$\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ in (20)–(22) is given by [15,
Eq. (24)]:

$\mathcal{T}_{\mathrm{pFB}}(\bm{\omega}):=(\mathrm{Id}+\Phi^{-1}(\mathcal{A}+\mathcal{C}))^{-1}(\mathrm{Id}-\Phi^{-1}\mathcal{B})(\bm{\omega}),$
(23)

where $\Phi\succ 0$ is a symmetric positive definite preconditioning matrix,
defined as

$\Phi:=\begin{bmatrix}\rho^{-1}&-\bm{A}^{\top}&0\\\
-\bm{A}&\tau^{-1}&-L\otimes I_{m}\\\ 0&-L\otimes
I_{m}&\sigma^{-1}\end{bmatrix},$

with $\rho,\theta,\tau\in\mathbb{R}_{>0}^{N}$ being step sizes similarly
defined as those of the FBF algorithm. Then, we can have an extension of the
pFB for the v-GNE optimal selection of cocoercive games, as stated in
Algorithm 2 with step size rules given in Assumptions 4 and 11. Finally, we
formally state the convergence property of Algorithm 2 in Theorem 2.

######  Assumption 11 ( [15, Eq. (27) and Thm. 3]).

It holds that
$\rho_{i}\leq(\max_{j=1,\dots,n_{i}}\sum_{k=1}^{m}|[A_{i}]_{jk}|\delta)^{-1}$,
$\tau_{i}\leq(\max_{j=1,\dots,n_{i}}\sum_{k=1}^{m}|[A_{i}]_{jk}|+2|\mathcal{N}_{i}^{\lambda}|+\delta)^{-1}$,
and $\sigma_{i}\leq(2|\mathcal{N}_{i}^{\lambda}|+\delta)^{-1}$, for all
$i\in\mathcal{I}$, where
$\delta>1/(\min(\eta,(2\max_{i\in\mathcal{I}}|\mathcal{N}_{i}^{\lambda}|)^{-1}))$.

######  Theorem 2.

Let Assumptions 1–4, 7, and 9–11 hold. Let $\Omega^{\star}$ be the set of
solutions to Problem (16) with $\mathcal{T}=\mathcal{T}_{\mathrm{pFB}}$
defined in (23), where $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ are
defined in (20)–(22). Furthermore, let $(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$,
where $\bm{\omega}^{(k)}=(\bm{x}^{(k)},\bm{\lambda}^{(k)},\bm{\nu}^{(k)})$, be
the sequence generated by Algorithm 2. Then,
$\lim_{k\to\infty}\mathrm{dist}(\bm{\omega}^{(k)},\Omega^{\star})=0,$ and
$(\bm{x}^{(k)})_{k\in\mathbb{N}}$ converges to an optimal v-GNE of the game in
(4).

Algorithm 2 Optimal v-GNE selection via pFB and HDSM for linearly coupled
cocoercive games

Initialization. Set $x_{i}^{(0)}\in\mathcal{X}_{i}$,
$\lambda_{i}^{(0)}\in\mathbb{R}_{\geq 0}^{m}$, and
$\nu_{i}^{(0)}\in\mathbb{R}^{m}$, for all $i\in\mathcal{I}$.

Iteration of each agent $i\in\mathcal{I}$.

  1. 1.

Receives $x_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{J}$ and
$\lambda_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{\lambda}$.

  2. 2.

Updates:

$\displaystyle\accentset{\circ}{x}_{i}^{(k)}$
$\displaystyle={\mathrm{prox}_{\ell_{i}+\iota_{\mathcal{X}_{i}}}^{\rho_{i}}}\left(x_{i}^{(k)}-\rho_{i}(\nabla_{x_{i}}{f_{i}}(\bm{x}^{(k)})+A_{i}^{\top}\lambda_{i}^{(k)})\right),$
$\displaystyle\accentset{\circ}{\nu}_{i}^{(k)}$
$\displaystyle=\nu_{i}^{(k)}-\sigma_{i}\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\lambda_{i}^{(k)}-\lambda_{j}^{(k)}\right).$

  3. 3.

Receives $\tilde{\nu}_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{\lambda}$.

  4. 4.

Updates:

$\displaystyle\accentset{\circ}{\lambda}_{i}^{(k)}=\operatorname{proj}_{\geq
0}\Big{(}\Big{.}\lambda_{i}^{(k)}+\tau_{i}\Big{(}\Big{.}A_{i}(2x_{i}^{(k+1)}-x_{i}^{(k)})-b_{i}$
$\displaystyle\qquad\quad+\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(2\nu_{i}^{(k+1)}-2\nu_{j}^{(k+1)}-\nu_{i}^{(k)}+\nu_{j}^{(k)}\right)$
$\displaystyle\qquad\quad+\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\lambda_{i}^{(k)}-\lambda_{j}^{(k)}\right)\Big{.}\Big{)}\Big{.}\Big{)}.$

  5. 5.

Sends
$(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})$
to a coordinator and receives back
$\nabla_{\omega_{i}}\phi(\accentset{\circ}{\bm{x}}^{(k)},\accentset{\circ}{\bm{\lambda}}^{(k)},\accentset{\circ}{\bm{\nu}}^{(k)})$,
where $\omega_{i}=(x_{i},\lambda_{i},\nu_{i})$.

  6. 6.

Updates:

$\displaystyle(x_{i}^{(k+1)},\lambda_{i}^{(k+1)},\nu_{i}^{(k+1)})$
$\displaystyle=(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})-\beta^{(k)}\nabla_{\omega_{i}}\phi(\accentset{\circ}{\bm{x}}^{(k)},\accentset{\circ}{\bm{\lambda}}^{(k)},\accentset{\circ}{\bm{\nu}}^{(k)}).$
(24)

###### Proof.

See Appendix D-A. ∎

##  V Online tracking of optimal generalized Nash equilibria

###  V-A Online optimal equilibrium tracking problem

In the second part of this paper, we consider the online GNE selection
problem. Specifically, let us introduce the time-varying game:

$\displaystyle\underset{x_{i}\in\mathcal{X}_{i,t}}{\min}\quad$ $\displaystyle
J_{i,t}(\bm{x})$ (25a) $\displaystyle\operatorname{s.t.}\quad$
$\displaystyle\sum_{j\in\mathcal{I}}g_{j,t}(x_{j})\leq 0,$ (25b)

where $t$ denotes the time index. The problem is time-varying in the sense
that the objective functions of the agents, as well as the constraints, may
vary over time. We assume that each instance of the games in (V-A) satisfies
Assumptions 1 and 2. The time-varying GNE selection problem thus concerns the
tracking of the sequence $(\bm{\omega}^{\star}_{t})_{t\in\mathbb{N}}$:

$\displaystyle\underset{\bm{\omega}}{\operatorname{argmin}}~{}$
$\displaystyle\phi_{t}(\bm{\omega})$ (26a)
$\displaystyle\operatorname{s.t.}~{}$
$\displaystyle\bm{\omega}\in\operatorname{zer}(\mathcal{A}_{t}+\mathcal{B}_{t}+\mathcal{C}_{t}).$
(26b)

The problems in (V-A) and (26) are a sequence in time of instances of (4) and
(13), respectively. The operators $\mathcal{A}_{t}$, $\mathcal{B}_{t}$, and
$\mathcal{C}_{t}$ are defined in (9)–(11), for the game in (V-A) at time step
$t$. The agents need to compute the action $\bm{\omega}_{t+1}$, having only
access to the game formulation up to time $t$. This setup describes the case
in which the agents act in a variable environment with limited computation
capabilities, so that they cannot compute the exact optimal selection before
changes in the problem (either in the selection function or in the game)
occur.

For every $t{\in\mathbb{N}}$, and under a suitable choice of operator
$\mathcal{T}_{t}$, such that

$\bm{\omega}\in\operatorname{zer}(\mathcal{A}_{t}+\mathcal{B}_{t}+\mathcal{C}_{t})\Leftrightarrow\bm{\omega}\in\operatorname{fix}(\mathcal{T}_{t}),$

$\bm{\omega}^{\star}_{t}$ (26) can be equivalently found as the solution of
the time-varying fixed-point selection problem

$\inf_{\bm{\omega}\in\operatorname{fix}(\mathcal{T}_{t})}\langle\bm{\omega}-\bm{\omega}^{\star}_{t},\nabla\phi_{{t}}(\bm{\omega}^{\star}_{t})\rangle\geq
0.$ (27)

The sequence $(\bm{\omega}^{\star}_{t})_{t\in\mathbb{N}}$ is well defined
when, for each $t\in\mathbb{N}$, the solution of (26) is unique. Let us then
introduce the following assumptions, which guarantee uniqueness.

######  Assumption 12.

The selection function $\phi_{t}\colon\mathbb{R}^{n_{\omega}}\to{\mathbb{R}}$
in (27) is continuously differentiable, $\sigma$-strongly convex, and has
$L_{\phi}$-Lipschitz continuous gradient for all $t\in\mathbb{N}$.

######  Assumption 13.

The operator $\mathcal{T}_{t}$ in (27) is quasi nonexpansive with
$\operatorname{fix}(\mathcal{T}_{t})\neq\varnothing$ for all $t\in\mathbb{N}$.

Under Assumptions 12 and 13, by [28, Prop. 1], we find
$\operatorname{fix}(\mathcal{T}_{t})$ to be closed and convex for all $t$. By
[26, Thm. 2.3.3], the problem in (27) has a unique solution for all $t$. These
assumptions also guarantee, by [38, Thm. 2A.7], that the solutions of (27)
coincide with the solutions of (26). In the remainder of this section, we
build upon the results of Section III-C to derive an HSDM-inspired algorithm
for tracking $(\bm{\omega}^{\star}_{t})_{t\in\mathbb{N}}$.

###  V-B Online fixed point tracking via the restarted Hybrid Steepest Descent
Method

The existing results on the HSDM algorithm study the asymptotic behavior with
vanishing step size $(\beta^{(k)})_{k\in\mathbb{N}}$ (see Assumption 4).
However, in online scenarios, decision makers may not have the computational
capability to exactly compute the fixed point of the algorithm, since that
would require an infinite amount of iterations in a limited time span before a
new instance of the problem becomes available. Thus, we propose and study the
(approximate) convergence properties of an algorithm that only performs a
finite number of HSDM iterations per time step. Consequently, the sequence of
step sizes becomes truncated and a sequence of vanishing step sizes, which is
required for the convergence of the HSDM, cannot be defined. We therefore
simplify the analysis by considering a constant sequence of step sizes.

Let us introduce the restarted HSDM algorithm. Given an initial state
$\bm{\omega}_{1}$, for each $t\in\mathbb{N}$, we propose the following:

$\displaystyle\bm{y}^{(k+1)}$
$\displaystyle:=\begin{cases}\bm{\omega}_{t},&\text{for}~{}k=1,\\\
\mathcal{T}_{t}(\bm{y}^{(k)})-\beta\nabla\phi_{t}(\mathcal{T}_{t}(\bm{y}^{(k)})),&\text{for}~{}k=2,...,K,\end{cases}$
$\displaystyle\bm{\omega}_{t+1}$ $\displaystyle:=\bm{y}^{(K+1)}.$ (28)

In words, at each time step $t$ the auxiliary variable $\bm{y}^{(k)}$, with
$k=1,...K$, is updated with $K$ iterations of the HSDM. Then, the decision
variable at time step $t+1$ is obtained as $\bm{\omega}_{t+1}=\bm{y}^{(K+1)}$.
The algorithm is then restarted when the information on the selection function
and game for the next time step becomes available. Next, let us postulate the
following technical assumptions:

######  Assumption 14.

There exists a compact set $\mathcal{Y}$ such that
$\bm{\omega}^{\star}_{t}\in\mathcal{Y}$ for all $t\in\mathbb{N}.$

######  Assumption 15.

There exists $U\geq 0$ such that
$\sup_{\bm{\omega}\in\bigcup_{{\tau}\in\mathbb{N}}\mathrm{Im}(\mathcal{T}_{{\tau}}),t\in\mathbb{N}}\|{\nabla}\phi_{t}(\bm{\omega})\|\leq
U.$

Assumption 14 is practically reasonable, since we can assume that we do not
aim at tracking a divergent sequence. Assumption 15 specifies an upper bound
for the gradient of the selection function and is in line with the online
optimization literature (see [39, Assm. 5], [40, Assm. 5], among others).

As shown in Section III-C, the HSDM method converges to the solution of a
selection problem over the fixed point set of a quasi-shrinking operator. In
the online scenario, assuming the operator $\mathcal{T}_{t}$ to be quasi-
shrinking for all $t$ is not enough, as the quasi-shrinking property might not
hold asymptotically. Thus, we also postulate the following technical
assumption:

######  Assumption 16.

(Uniformly quasi-shrinking operator) For any closed convex set $C$ such that
$C\cap\operatorname{fix}(\mathcal{T}_{t})\neq\varnothing$, there exists
$D:\mathbb{R}_{\geq 0}\to\mathbb{R}$ positive semidefinite such that
$D_{t}(r)\geq D(r)$ for all $t{\in\mathbb{N}}$ and for all $r\geq 0$, where
$D_{t}(\cdot)$ is the shrinkage function of $\mathcal{T}_{t}$ defined as in
(2).

######  Remark 4.

Assumption 16 implies that $\mathcal{T}_{t}$ is quasi-shrinking on any closed,
convex set $C$ such that
$\mathcal{C}\cap\operatorname{fix}(\mathcal{T}_{t})\neq\varnothing$, $\forall
t\in\mathbb{N}$.

The next lemma outlines a contraction property of the restarted HSDM to the
solution sequence of Problem (27) up to an additive error, which can be
controlled by an appropriate choice of the step size $\beta$ and the number of
iterations $K$.

######  Lemma 5.

Let Assumptions 12–16 hold. For any $t\in\mathbb{N}$, let $\bm{\omega}_{t+1}$
be generated by the restarted HSDM algorithm in (28). For any $\gamma>0$,
there exist $K,\beta>0$, such that

$\|\bm{\omega}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\left(1-{\tau(\beta)}\right)^{K}\|\bm{\omega}_{t}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+\gamma,$
(29)

where $\tau(\beta):=1-\sqrt{1-\beta(2\sigma-\beta L_{\phi}^{2})}\in(0,1).$

###### Proof.

See Appendix E-B. ∎

######  Remark 5.

For decreasing values of the tolerable error $\gamma$, the stepsize $\beta$
has to be decreased and the number of iterations $K$ has to be increased (see
the proof of Lemma 5).

We now proceed to show how the property in (29) can be exploited to derive an
error bound on the trajectory tracking of the solution sequence of the problem
in (27) via the restarted HSDM (28). Thus, we introduce the following
assumption:

######  Assumption 17.

There exist scalars $\delta_{1},\delta_{2}\geq 0$ such that

  1. (i)

$\sup_{t\in\mathbb{N}}\|\bm{\bm{\omega}^{\star}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|\leq\delta_{1}$;

  2. (ii)

$\sup_{t\in\mathbb{N}}\mathrm{dist}(\bm{\omega}^{\star}_{t},(\operatorname{fix}(\mathcal{T}_{t+1}))\leq\delta_{2}$.

Assumption 17.i is standard in online optimization (e.g. [23, Assm. 1], [32,
Assm. 3.1], and [40, Assm. 3]). We note that Assumption 17.i implies
Assumption 17.ii. The latter is nevertheless introduced to distinguish the
effects of the time variation of $\mathcal{T}_{t}$ (which influences both
$\delta_{1}$ and $\delta_{2}$) from the one of $\phi_{t}$ (which only
influences $\delta_{1}$).

######  Remark 6.

If $\mathcal{T}_{t}=\mathcal{T}$, for all $t\in\mathbb{N}$, and the time
dependence can be expressed through a parametrization, that is,
$\phi_{t}(\bm{\omega})=\phi(\bm{\omega},t)$, then an estimate for $\delta_{1}$
can be found. In fact, if $\phi(\bm{\omega},t)$ is continuously
differentiable, we find by [38, Thm. 2F.7] that the solution mapping, that is,
the mapping from $t$ to the solution of
$\mathrm{VI}(\operatorname{fix}(\mathcal{T}),\nabla_{x}\phi(\cdot,t))$, is
Lipschitz continuous in a neighbourhood of any ${t}$ with Lipschitz constant
$\sigma^{-1}|\nabla_{t}\phi(\bm{\omega}^{\star}_{{t}},{t})|$. Thus, if the
time variation between two consecutive time steps $t_{1}$ and $t_{2}$ is small
enough, $\delta_{1}$ can be estimated as
$\sigma^{-1}|\nabla_{t}\phi(\bm{\omega}^{\star}_{t_{1}},{t_{1}})|(t_{2}-t_{1})$.
The solution mapping is in general discontinuous when $\mathcal{T}_{t}$ is
time-varying; thus, a similar estimate cannot be found in the general case.

######  Theorem 3.

Let Assumptions 12–17 hold. Let the sequence
$(\bm{\omega}_{t})_{t\in\mathbb{N}}$ be generated by the restarted HSDM in
(28). For any $\gamma>0$, there exist
$\beta\in(0,\frac{2\sigma}{L^{2}_{{\phi}}})$ and $\bar{K}$, such that, for all
$K\geq\bar{K}$, the sequence $(\bm{\omega}_{t})_{t\in\mathbb{N}}$ is bounded
and

$\limsup_{t\rightarrow\infty}\|\bm{\omega}_{t}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\frac{({\gamma}+\delta_{1}^{2})}{1/2-\alpha},$
(30)

where $\alpha=(1-\tau(\beta))^{K}<\frac{1}{2}$.

###### Proof.

See Appendix E-C. ∎

######  Remark 7.

In Theorem 3, $\gamma$ is derived from the additive error in (29). Thus, to
control the approximation error in (30), $\beta$ must be chosen small so to
obtain small values of ${\gamma}$, as pointed out in Remark 5. However, the
value $\tau(\beta)$ tends to $0$ for small values of $\beta$. This leads to
the denominator in (30) to be small for small stepsizes, unless the number of
iterations $K$ is increased. Therefore, a smaller step size leads to a better
approximation error only if it is shouldered by an increase in the number of
iterations of the algorithm per time step.

In summary, we find that the restarted HSDM (28) asymptotically tracks the
solutions trajectory of the online fixed point selection problem in (27), with
an asymptotic error that can be controlled up to the variability of the
problem $\delta_{1}$, via an appropriate choice of $\beta$, $K$, as shown in
Theorem 3. Additionally, we emphasize that the results hold for a more general
problem, i.e., one could replace $\nabla\phi_{t}$ in Problem (27) with a
strongly monotone operator to obtain an extension to the fixed-point selection
problem in [28]. In the next section, we use the restarted HSDM to solve the
online GNE tracking problem in (26).

###  V-C Distributed optimal equilibrium tracking algorithm for monotone games

We recall from Section III-C that the set of variational GNEs for a monotone
game can be characterized as the set of fixed points of the operator
$\mathcal{T}_{\mathrm{FBF}}$ defined in (18). Thus, for the time-varying game
in (V-A) at time $t$, let $\mathcal{T}_{\mathrm{FBF},t}$ be the FBF operator
defined as:

$\displaystyle\mathcal{T}_{\mathrm{FBF},t}(\bm{\omega})$
$\displaystyle:=((\mathrm{Id}-\Psi^{-1}(\mathcal{B}_{t}+\mathcal{C}_{t}))(\mathrm{Id}+{\Psi^{-1}\mathcal{A}_{t}})^{-1}$
$\displaystyle\quad\cdot(\mathrm{Id}-\Psi^{-1}(\mathcal{B}_{t}+\mathcal{C}_{t}))+\Psi^{-1}(\mathcal{B}_{t}+\mathcal{C}_{t}))(\bm{\omega}),$
(31)

where $\mathcal{A}_{t}$, $\mathcal{B}_{t}$, and $\mathcal{C}_{t}$ are those in
Problem (26) and associated with the game in (V-A) at time $t$. The solutions
of the time-varying GNE selection problem in $\eqref{eq:online_GNE_selection}$
are equivalent to the solutions of (27), with
$\mathcal{T}_{t}=\mathcal{T}_{\mathrm{FBF},t}$ for all $t$. By Lemma 3,
$\mathcal{T}_{\mathrm{FBF},t}$, for each $t$, is a quasi-nonexpansive, quasi-
shrinking operator. Therefore, the restarted HSDM algorithm in (28) can be
employed for tracking the solution trajectory, with an asymptotic tracking
error given by Theorem 3. We introduce an assumption for the GNE selection
problem, which is equivalent to Assumption 17.ii:

######  Assumption 18.

There exists a scalar $\delta_{2}\geq 0$ such that
$\sup_{t\in\mathbb{N}}\mathrm{dist}(\bm{\omega}^{\star}_{t},\mathrm{zer}(\mathcal{A}_{t+1}+\mathcal{B}_{t+1}+\mathcal{C}_{t+1}))\leq\delta_{2}$.

######  Corollary 1.

Let us consider the online GNE tracking problem in (26) for the time-varying
game in (V-A) that satisfies Assumptions 1, 2, 7, for each $t\in\mathbb{N}$.
Suppose that Assumptions 12, 14, 15, 17, 18 hold. Let
$\mathcal{T}_{t}=\mathcal{T}_{\mathrm{FBF},t}$ satisfy Assumption 16. Then,
for any $\gamma>0$ there exist $\beta\in(0,\frac{2\sigma}{L_{\phi}^{2}})$ and
$\bar{K}$ such that, for any $K\geq\bar{K}$, the asymptotic tracking error of
Algorithm 3 is given by (30).

###### Proof.

See Appendix E-D. ∎

######  Remark 8.

In Corollary 1, Assumption 16 is satisfied for example when at every time step
$t$, the feasible set of Problem (26) is selected among the GNE sets of
finitely many games. That is, consider a finite set of operators

$\mathcal{A}_{h},\mathcal{B}_{h},\mathcal{C}_{h},~{}~{}\text{with}~{}h\in\\{1,...,H\\},$

and for each $h$, the associated FBF operator
$\mathcal{T}_{\mathrm{FBF}}^{h}$. Defining a mapping from the time step $t$ to
the indexes of the operators $\eta:\mathbb{N}\to\\{1,...,H\\},$ Problem (26)
is defined by

$\bm{\omega}^{\star}_{t}:=\left\\{\begin{array}[]{rl}&\operatorname{argmin}_{\bm{\omega}}\phi_{t}(\bm{\omega})\\\
&\text{s.t.}~{}~{}\bm{\omega}\in\operatorname{zer}(\mathcal{A}_{\eta(t)}+\mathcal{B}_{\eta(t)}+\mathcal{C}_{\eta(t)}).\end{array}\right.$

Let us denote with $D^{h}(\cdot)$ the shrinkage function of
$\mathcal{T}_{\mathrm{FBF}}^{h}$. By Lemma 1, $\mathcal{T}_{\mathrm{FBF}}^{h}$
is quasi-shrinking and, therefore, $D^{h}(\cdot)$ is positive semidefinite.
Assumption 16 is then satisfied with $D(r)=\min_{h\in\\{1,...,H\\}}D^{h}(r)$.
This problem class includes the case when only the selection function ${\phi}$
varies, i.e., $H=1$.

Algorithm 3 Optimal v-GNE tracking via FBF and HSDM

Initialization. Set $x_{i,0}\in\mathcal{X}_{i}$,
$\lambda_{i,0}\in\mathbb{R}_{\geq 0}^{m}$, and $\nu_{i,0}\in\mathbb{R}^{m}$,
for all $i\in\mathcal{I}$.
Iteration at time $t\in\mathbb{N}_{0}$ of each agent $i\in\mathcal{I}$:

  1. 1.

Receives $J_{i,t}(\cdot)$, $g_{i,t}(\cdot)$, and $\mathcal{X}_{i,t}(\cdot)$.

  2. 2.

Assigns $\hat{x}_{i}^{(1)}\leftarrow x_{i,t}$,
$\hat{\lambda}_{i}^{(1)}\leftarrow\lambda_{i,t}$, and
$\hat{\nu}_{i}^{(1)}\leftarrow\nu_{i,t}$.

  3. 3.

For $k=1,...,K$:

  1. (i)

Receives $\hat{x}_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{J}$ and
$\hat{\lambda}_{j}^{(k)},\hat{\nu}_{j}^{(k)}$ from agent
$j\in\mathcal{N}_{i}^{\lambda}$.

  2. (ii)

Updates:

$\displaystyle\tilde{x}_{i}^{(k)}$
$\displaystyle={\mathrm{prox}_{\ell_{i,t}+\iota_{\mathcal{X}_{i,t}}}^{\rho_{i}}}\left(\hat{x}_{i}^{(k)}-\rho_{i}(\nabla_{x_{i}}{f}_{i,t}(\hat{\bm{x}}^{(k)})\right.$
$\displaystyle\left.\quad+\nabla
g_{i,t}(\hat{x}_{i}^{(k)})^{\top}\hat{\lambda}_{i}^{(k)})\right),$
$\displaystyle\tilde{\lambda}_{i}^{(k)}$
$\displaystyle=\operatorname{proj}_{\geq
0}\Big{(}\Big{.}\hat{\lambda}_{i}^{(k)}+\tau_{i}\Big{(}\Big{.}g_{i,t}(\hat{x}_{i}^{(k)})$
$\displaystyle\quad+\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\hat{\nu}_{i}^{(k)}-\hat{\nu}_{j}^{(k)}-\hat{\lambda}_{i}^{(k)}+\hat{\lambda}_{j}^{(k)}\right)\Big{.}\Big{)}\Big{.}\Big{)},$
$\displaystyle\tilde{\nu}_{i}^{(k)}$
$\displaystyle=\hat{\nu}_{i}^{(k)}-\sigma_{i}\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\hat{\lambda}_{i}^{(k)}-\hat{\lambda}_{j}^{(k)}\right).$

  3. (iii)

Receives $\tilde{x}_{j}^{(k)}$ from agent $j\in\mathcal{N}_{i}^{J}$ and
$\tilde{\lambda}_{j}^{(k)},\tilde{\nu}_{j}^{(k)}$ from agent
$j\in\mathcal{N}_{i}^{\lambda}$.

  4. (iv)

Updates:

$\displaystyle\begin{split}\accentset{\circ}{x}_{i}^{(k)}=&\tilde{x}_{j}^{(k)}-\rho_{i}\Big{(}\nabla_{x_{i}}f_{i,t}(\tilde{{\bm{x}}}^{(k)})-\nabla_{x_{i}}f_{i,t}(\hat{\bm{x}}^{(k)})+\\\
&\nabla g_{i,t}(\tilde{x}_{i}^{(k)})^{\top}\tilde{\lambda}_{i}^{(k)}-\nabla
g_{i,t}(\hat{x}_{i}^{(k)})^{\top}\hat{\lambda}_{i}^{(k)}\Big{)},\\\
\accentset{\circ}{\lambda}_{i}^{(k)}=&\tilde{\lambda}_{i}^{(k)}+\tau_{i}\Big{(}\Big{.}g_{i,t}(\tilde{x}_{i}^{(k)})-g_{i,t}(\hat{x}_{i}^{(k)})+\\\
&\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\Big{(}\Big{.}\tilde{\nu}_{i}^{(k)}-\hat{\nu}_{i}^{(k)}-\tilde{\nu}_{j}^{(k)}+\hat{\nu}_{j}^{(k)}\Big{.}\Big{)}-\\\
&\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\Big{(}\Big{.}\tilde{\lambda}_{i}^{(k)}-\hat{\lambda}_{i}^{(k)}-\tilde{\lambda}_{j}^{(k)}+\hat{\lambda}_{j}^{(k)}\Big{.}\Big{)}\Big{.}\Big{)},\\\
\accentset{\circ}{\nu}_{i}^{(k)}=&\tilde{\nu}_{i}^{(k)}-\sigma_{i}\textstyle\sum_{j\in\mathcal{N}_{i}^{\lambda}}\left(\tilde{\lambda}_{i}^{(k)}-\hat{\lambda}_{i}^{(k)}-\tilde{\lambda}_{j}^{(k)}+\hat{\lambda}_{j}^{(k)}\right).\end{split}$

  5. (v)

Sends
$(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})$
to a coordinator and receives
$\nabla\phi^{t}_{\bm{\omega}_{i}}(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})$,
where $\bm{\omega}_{i}=(x_{i},\lambda_{i},\nu_{i})$.

  6. (vi)

Updates:

$\displaystyle\begin{split}&(\hat{x}_{i}^{(k+1)},\hat{\lambda}_{i}^{(k+1)},\hat{\nu}_{i}^{(k+1)})\\\
&=(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)})-\beta\nabla\phi^{t}_{\bm{\omega}_{i}}(\accentset{\circ}{x}_{i}^{(k)},\accentset{\circ}{\lambda}_{i}^{(k)},\accentset{\circ}{\nu}_{i}^{(k)}).\end{split}$

End For

  1. 4)

Assigns ${x}_{i,t}\leftarrow\hat{x}_{i}^{(K+1)}$,
${\lambda}_{i,t}\leftarrow\hat{\lambda}_{i}^{(K+1)}$,
${\nu}_{i,t}\leftarrow\hat{\nu}_{i}^{(K+1)}$.

######  Remark 9.

The result of this section holds similarly if we substitute the FBF operator
with the pFB operator in (23), which is quasi-shrinking (see the proof of
Theorem 2), for cocoercive games with affine coupling constraints.

##  VI Illustrative example

We consider a peer-to-peer electricity market clearing problem with
operational constraints of the electrical network, adapted from [2]. We assume
that each bus of a distribution network consists of one agent that has access
to either a storage unit or a dispatchable generation unit. Each agent
$i\in\mathcal{I}$ has decision authority on the power generated
$p^{\text{g}}_{i,h}$, the power bought from the main grid
$p^{\text{mg}}_{i,h}$, the power drawn from the storage unit
$p^{\text{st}}_{i,h}$, the power traded with the trading partners
$p^{\mathrm{tr}}_{(i,j),h},~{}j\in\mathcal{N}_{i}$ and the phase at the bus
$\theta_{i,h}$ over the horizon $h=1,...,H$. Let us denote
$\bm{x}_{i,h}=\mathrm{col}(p^{\text{g}}_{i,h},p^{\text{mg}}_{i,h},p^{\text{st}}_{i,h},\\{p^{\mathrm{tr}}_{(i,j),h}\\}_{j\in\mathcal{N}_{i}},\theta_{i,h})$,
for all $i\in\mathcal{I}$ and $h=1,...,H$, and denote
$\bm{x}_{i}:=\mathrm{col}(\\{\bm{x}_{i,h}\\}_{h=1,...,H})$,
$\bm{x}:=\mathrm{col}(\\{\bm{x}_{i}\\}_{i\in\mathcal{I}})$. Each agent aims at
minimizing its local cost function [2, Eq. (17)]:

$\displaystyle\begin{split}J_{i}(\bm{x})=&\sum_{h=1}^{H}f^{\text{g}}_{i,h}(p^{\text{g}}_{i,h})+f^{\text{tr}}_{i,h}(\\{p^{\text{tr}}_{(i,j),h}\\}_{j\in\mathcal{N}_{i}})\\\
&+f^{\text{mg}}_{i,h}(p^{\text{mg}}_{i,h},p^{\text{mg}}_{-i,h}),\end{split}$
(32)

where $f^{\text{tr}}_{i,h}$ encodes the cost or revenue of the trading with
other agents and $f^{\text{mg}}_{i,h}$ encodes the cost of purchasing energy
from the main grid as in [2, Eq. (11)], while $f^{\text{g}}_{i,h}$ is a linear
function which encodes the cost of power generation. The local feasible sets
$\mathcal{X}_{i},~{}i=1,...,N$ include the satisfaction of the power demand at
the bus, as well as the operating constraints of the generators and storage
units. The shared constraints are of the form $g(\bm{x})\leq\bm{0}_{n_{c}},$
with $g$ affine. They include the operating limits of the grid, the trading
reciprocity
$\\{p^{\mathrm{tr}}_{(i,j),h}=-p^{\mathrm{tr}}_{(j,i),h},~{}\forall~{}i\in\mathcal{N},~{}\forall~{}j\in\mathcal{N}_{i}\\}$
and the linearized power flow equations with DC approximation
$\\{p^{\text{g}}_{i,h}+p^{\text{st}}_{i,h}+\iota^{\text{mg}}_{i}\sum_{j\in\mathcal{N}}p^{\text{mg}}_{j,h}+\sum_{j\in\mathcal{B}_{i}}B_{ij}(\theta_{i,h}-\theta_{j,h})=0\\}$,
where $\iota^{\text{mg}}_{i}\in\\{0,1\\}$ is $1$ if and only if $i$ is
connected to the main grid, $\mathcal{B}_{i}$ is the set of buses that are
connected to bus $i$ on the electric grid and $B$ is the susceptance matrix.
We note that the game satisfies Assumptions 1 and 2. In addition, we consider
the IEEE 13-bus distribution feeder for our numerical simulations, performed
in Matlab.
We first simulate the day-ahead market clearing (with 24 hourly time steps)
via the standard FBF-based algorithm, which can obtain a v-GNE, and Algorithm
1, which solves the optimal selection problem of this game. Specifically, we
consider the GNE selection function:

$\displaystyle\begin{split}&\phi(\bm{x})=\sum_{h=1}^{H}\\{\|\bm{p}^{\text{g}}_{h}-\bar{\bm{p}}^{\text{g}}\|_{Q_{\text{d}}}+\|\bm{p}^{\text{mg}}_{h}\|_{Q_{\text{mg}}}+\|\bm{\theta}_{h}-\bar{\bm{\theta}}\|_{Q_{\theta}}\\\
&+\|G\bm{\theta}_{h}\|_{Q_{\text{pf}}}+\|\bm{p}^{\text{tr}}_{h}\|_{Q_{\text{tr}}}+\|\bm{p}^{\text{st}}_{h}\|_{Q_{\text{st}}}\\}+\|\bm{\lambda}\|_{Q_{\lambda}}+\|\bm{\nu}\|_{Q_{\nu}},\end{split}$
(33)

where we denoted in bold the column stack of the respective variables for each
agent and the matrices $Q_{\star}$ are diagonal positive definite. We choose
$\bar{\bm{p}}^{\text{g}}$ to be the column vector of the maximum generation
production for each agent, in order to maximize the renewable energy
production, and $\bar{\bm{\theta}}$ to be a vector which elements are all
equal to the phase of the node connected to the main grid, in order to reduce
the grid imbalances. The cost factors related to
${\bm{p}}^{\text{mg}},{\bm{p}}^{\text{st}},{\bm{p}}^{\text{tr}}$ aim at
reducing the burden on the transmission grid, increasing the lifespan of the
storage units and reducing the load of the trading platform, respectively. The
terms in $\bm{\lambda}$ and $\bm{\nu}$ act as regularization of the dual
variables. Finally, $G$ is a matrix that maps the phase of the nodes to the
power flowing through the lines. In this test, we aim at maximizing the
lifespan of the grid lines by setting the non-zero elements of $Q_{\text{pf}}$
to be large. The solution obtained by Algorithm 1 and that of the standard FBF
are depicted in Figure 1. As expected, since the v-GNE computed by Algorithm 1
minimizes the selection function (33), it has a lower load on the power lines
than that of the standard FBF.

Figure 1: Power flowing through the lines (day-ahead market).

Secondly, we test Algorithm 3 on a real-time market scenario, formulated as a
time-varying game. The horizon is set to 2 hours, with a sampling time of $15$
minutes. The simulation is run over a 24 hour span for different values, thus
resulting in 12 consecutive instances of GNE selection problems. Let us index
these problems with $t=1,...,12$. The cost function of each agent is given by
(32), with an additional term
$f^{\text{st}}_{i,t}(\\{p^{\text{st}}_{i,h}\\}_{h=1,...,H})=\|{x}_{\text{ch},t}^{i}-\bar{x}_{\text{ch},t}^{i}\|^{2}_{2}$.
This term penalizes the deviations of the storage units charge state at the
end of the $t$-th horizon from the charge state planned in the day-ahead
market clearance $\bar{x}_{\text{ch},t}^{i}$. The charge state at the end of
the $t$-th horizon is given by
${x}_{\text{ch},t}^{i}={x}_{\text{ch},t-1}^{i}-\sum_{h=1}^{H}p^{\text{st}}_{i,h}$,
and the initial state ${x}^{i}_{\text{ch},0}$ is known. Because of the
variability along the day of the power demand, the local power balance
constraint defined in [2, Eq. (6)] depends on $t$. The cost functions and
constraints of the game are therefore time-varying, with $t$ representing the
time index. Furthermore, in this scenario, we aim at computing a v-GNE that
minimizes the power flowing on the line connecting buses 632 and 671 during
peak hours. Thus, we consider (33) as the selection function at each $t$ where
the element of $Q_{\text{pf}}$ related to this line is time-varying, i.e., it
is set high between 6AM and 4PM. We note that this setup falls into the case
considered in Remark 8, whilst $\\{\phi_{t}\\}_{t=1,...,12}$ satisfies
Assumption 12. We run the simulation for different values of the parameters
$K$ and $\beta$ and Figure 2 illustrates the results. An increasing $K$
results in a diminishing residual, that is, a better convergence to the GNE
set, and a diminishing power load on the penalized line during peak hours, as
expected from the imposed penalty term in the selection function. A
diminishing $\beta$ implies a slower reduction of the cost function, which
results in a higher cost for small values of $K$, as shown in Figure 2(b).
Figure 3 depicts some particular trajectories of the power flowing through the
penalized line with $\beta=5\cdot 10^{-4}$.

(a) Average residual $\|T_{\text{FBF}}(\bm{\omega}_{t})-\bm{\omega}_{t}\|$.

(b) Average power flow on the line connecting buses 632 and 671 during peak
hours (kW).

Figure 2: Algorithm performance for several restarted HSDM parameters.  Figure
3: Power flowing on the line connecting buses 632 and 671 (real-time market)
with ${\beta=5\cdot 10^{-4}}.$

##  VII Conclusion

The optimal generalized Nash equilibrium selection problem in monotone games
can be solved distributively by combining the hybrid steepest descent method
with an appropriate fixed-point operator. The key requirement to guarantee
convergence to the set of optimal generalized Nash equilibria is the quasi-
shrinking property, which holds true for certain fixed-point operators. The
hybrid steepest descent method can also be modified to track a time-varying
optimal generalized Nash equilibria. The resulting approach is suitable for
real-time decision making in multi-agent dynamic environments.

##  Appendix A Proof of Lemma 1

Let us proceed by contradiction. We assume that there exists $r>0$ such that
$D_{\Psi}(r)=0$. Then, by the definition of $D_{\Psi}(\cdot)$ in (2), there
exists a sequence
$(\bm{\omega}_{k})_{k\in\mathbb{N}}\in(\operatorname{fix}(\mathcal{T})_{\geq
r}^{\Psi})\bigcap C$ such that

$\lim_{k\xrightarrow[]{}\infty}\mathrm{dist}_{\Psi}(\bm{\omega}_{k},\operatorname{fix}(\mathcal{T}))-\mathrm{dist}_{\Psi}(\mathcal{T}(\bm{\omega}_{k}),\operatorname{fix}(\mathcal{T}))=0.$

By the definition of projection, we have

$\displaystyle\mathrm{dist}_{\Psi}(\mathcal{T}(\bm{\omega}_{k}),\operatorname{fix}(\mathcal{T}))$
$\displaystyle=\|\mathcal{T}(\bm{\omega}_{k})-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\mathcal{T}(\bm{\omega}_{k}))\|_{\Psi}$
$\displaystyle\leq\|\mathcal{T}(\bm{\omega}_{k})-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}.$
(34)

By the quasi-nonexpansiveness of $\mathcal{T}$ and the latter inequality,

$\displaystyle\begin{split}0&\leq\underbrace{\|\bm{\omega}_{k}-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}}_{=\mathrm{dist}_{\Psi}(\bm{\omega}_{k},\operatorname{fix}(\mathcal{T}))}-\|\mathcal{T}(\bm{\omega}_{k})-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}\\\
&\leq\mathrm{dist}_{\Psi}(\bm{\omega}_{k},\operatorname{fix}(\mathcal{T}))-\mathrm{dist}_{\Psi}(\mathcal{T}(\bm{\omega}_{k}),\operatorname{fix}(\mathcal{T}))\xrightarrow[]{k\xrightarrow{}\infty}0.\end{split}$

It follows that

$\lim_{k\xrightarrow{}\infty}\|\bm{\omega}_{k}-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}-\|\mathcal{T}(\bm{\omega}_{k})-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}=0.$

By (3), we then have that

$\displaystyle\begin{split}&\|\bm{\omega}_{k}-\mathcal{T}_{2}(\bm{\omega}_{k})\|^{2}_{\Psi}\leq\\\
&\tfrac{1}{\gamma}(\|\bm{\omega}_{k}-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|^{2}_{\Psi}-\|\mathcal{T}(\bm{\omega}_{k})-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|^{2}_{\Psi})\leq\\\
&\tfrac{2d}{\gamma}(\|\bm{\omega}_{k}-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}-\|\mathcal{T}(\bm{\omega}_{k})-\operatorname{proj}^{\Psi}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}_{k})\|_{\Psi}),\end{split}$

where the latter inequality follows from $a^{2}-b^{2}=(a-b)(a+b)$ for
$a,b\in\mathbb{R}$ and where we substituted $d:=\sup_{\bm{\omega}\in
C}\|\bm{\omega}_{k}-\bm{\omega}\|_{\Psi}$, which is finite since the set $C$
is compact. We conclude that

$\lim_{k\xrightarrow{}\infty}\|\bm{\omega}_{k}-\mathcal{T}_{2}(\bm{\omega}_{k})\|^{2}_{\Psi}=0.$
(35)

By the Bolzano-Weierstrass theorem and the boundedness of $\bm{\omega}_{k}$,
there exists a convergent subsequence $(\bm{\omega}_{k_{l}})_{l\in\mathbb{N}}$
with accumulation point $\bm{\omega}^{\infty}$. By (35),
$\lim_{l\xrightarrow{}\infty}\mathcal{T}_{2}(\bm{\omega}_{k_{l}})=\bm{\omega}^{\infty}.$

By the demiclosedness of ${\mathrm{Id}}-\mathcal{T}_{2}$ and by
$\operatorname{fix}(\mathcal{T}_{2})\subset\operatorname{fix}(\mathcal{T})$,
$\bm{\omega}^{\infty}-\mathcal{T}_{2}(\bm{\omega}^{\infty})=0\Rightarrow\bm{\omega}^{\infty}\in\operatorname{fix}(\mathcal{T}_{2})\Rightarrow\bm{\omega}^{\infty}\in\operatorname{fix}(\mathcal{T}).$
However, since $(\operatorname{fix}(\mathcal{T})_{\geq r}^{\Psi})\bigcap C$ is
a closed set, then
$\bm{\omega}^{\infty}\in\operatorname{fix}(\mathcal{T})_{\geq r}^{\Psi}$,
which is in contradiction with
$\bm{\omega}^{\infty}\in\operatorname{fix}(\mathcal{T})$. $\hfill\blacksquare$

##  Appendix B Properties of operators $\mathcal{A}$, $\mathcal{B}$, and
$\mathcal{C}$ in (9)–(11)

######  Lemma 6.

Let Assumption 1 hold. Then, the operators $\mathcal{A}$, $\mathcal{B}$, and
$\mathcal{C}$ in (9)–(11) are maximally monotone. Thus,
$\mathcal{A}+\mathcal{B}+\mathcal{C}$ is also maximally monotone. $\hfill\Box$

###### Proof.

By Assumption 1, $\mathrm{N}_{{\mathcal{X}_{i}}}$ and $\partial\ell_{i}$ are
maximally monotone [27, Thm. 20.25 & Example 20.26]. The operator
$\mathcal{A}$ is thus maximally monotone by [27, Prop. 20.23 & Cor. 25.5]. The
operator $F$ is maximally monotone by Assumption 2 and by continuity in
Assumption 7. Meanwhile $L$ is a linear positive semidefinite operator and,
therefore, it is maximally monotone; thus, the operator $\mathcal{B}$ is
maximally monotone. We can write
$\mathcal{C}=\mathcal{C}_{1}+\mathcal{C}_{2}$, where
$\mathcal{C}_{1}=\operatorname{col}(\\{\langle\nabla_{x_{i}}g_{i}(x_{i}),\lambda_{i}\rangle\\}_{i\in\mathcal{I}},-\\{g_{i}(x_{i})\\}_{i\in\mathcal{I}},\mathbf{0}_{Nm})$
and $\mathcal{C}_{2}=\operatorname{col}(\mathbf{0}_{n},-(L\otimes
I_{m})\bm{\nu},(L\otimes I_{m})\bm{\lambda})$. The operator $\mathcal{C}_{1}$
is maximally monotone by continuity and by noting that, for any
$\bm{\omega},\bm{\omega}^{\prime}\in\mathbb{R}^{n}\times\mathbb{R}^{Nm}_{\geq
0}\times\mathbb{R}^{Nm}$,

$\displaystyle\langle\mathcal{C}_{1}(\bm{\omega})-\mathcal{C}_{1}(\bm{\omega}^{\prime}),\bm{\omega}-\bm{\omega}^{\prime}\rangle$
$\displaystyle=\textstyle\sum_{i\in\mathcal{I}}\langle
g_{i}(x_{i}^{\prime})-g_{i}(x_{i})-\nabla_{x_{i}}g_{i}(x_{i})^{\top}(x_{i}^{\prime}-x_{i}),\lambda_{i}\rangle$
$\displaystyle\quad+\textstyle\sum_{i\in\mathcal{I}}\langle
g_{i}(x_{i})-g_{i}(x_{i}^{\prime})-\nabla_{x_{i}}g_{i}(x_{i}^{\prime})^{\top}(x_{i}-x_{i}^{\prime}),\lambda_{i}^{\prime}\rangle\geq
0,$

where the inequality follows by the convexity of $g_{i}$. As $\mathcal{C}_{2}$
is a linear skew-symmetric operator, it is maximally monotone [27, Example
20.35]. By invoking [27, Cor. 25.5], the result follows. ∎

######  Lemma 7.

Let Assumptions 1 and 7 hold. Then the operators $\mathcal{B}$,
$\mathcal{C}$,and $\mathcal{B}+\mathcal{C}$, defined in (10)–(11), are
Lipschitz continuous.

###### Proof.

Due to Assumption 7, the operator $\mathcal{B}$ is $L_{F}$-Lipschitz
continuous. Lipschitz continuity of $\mathcal{C}$ can be evaluated as follows.
Similarly to the proof of Lemma 6, let us split
$\mathcal{C}=\mathcal{C}_{1}+\mathcal{C}_{2}$. The operator $\mathcal{C}_{2}$
is Lipschitz continuous by linearity, while Lipschitz continuity of
$\mathcal{C}_{1}$ is shown as follows. Let us denote the bound of
$\nabla_{x_{i}}g_{i}(x_{i})$ by $b_{\nabla g_{i}}$, i.e.,
$\|\nabla_{x_{i}}g_{i}(x_{i})\|\leq b_{\nabla g_{i}}$ (c.f. Assumption 7) and
the bound of $\lambda_{i}$ by $b_{\lambda_{i}}$, for all $i\in\mathcal{I}$,
which exists due to [35, Prop. 3.3]. For any
$\bm{\omega},\bm{\omega}^{\prime}\in\mathbb{R}^{n+2Nm}$,

$\displaystyle\|\mathcal{C}_{1}(\bm{\omega})-\mathcal{C}_{1}(\bm{\omega}^{\prime})\|^{2}$
$\displaystyle\overset{\\{1\\}}{\leq}\textstyle\sum_{i\in\mathcal{I}}\Big{(}\Big{.}2\|\nabla_{x_{i}}g_{i}(x_{i})^{\top}(\lambda_{i}-\lambda_{i}^{\prime})\|^{2}+\|g_{i}(x_{i})-g_{i}(x_{i}^{\prime})\|^{2}$
$\displaystyle\qquad+2\|(\nabla_{x_{i}}g_{i}(x_{i})-\nabla_{x_{i}}g_{i}(x_{i}^{\prime}))^{\top}\lambda_{i}^{\prime}\|^{2}\Big{.}\Big{)}$
$\displaystyle\overset{\\{2\\}}{\leq}\textstyle\sum_{i\in\mathcal{I}}\Big{(}\Big{.}2\|\nabla_{x_{i}}g_{i}(x_{i})^{\top}\|^{2}\|\lambda_{i}-\lambda_{i}^{\prime}\|^{2}+b_{\nabla
g_{i}}^{2}\|x_{i}-x_{i}^{\prime}\|^{2}$
$\displaystyle\qquad+2\|\lambda_{i}^{\prime}\|^{2}\|\nabla_{x_{i}}g_{i}(x_{i})-\nabla_{x_{i}}g_{i}(x_{i}^{\prime})\|^{2}\Big{.}\Big{)}$
$\displaystyle\overset{\\{3\\}}{\leq}\textstyle\sum_{i\in\mathcal{I}}\Big{(}2b_{\nabla
g_{i}}^{2}\|\lambda_{i}-\lambda_{i}^{\prime}\|^{2}+(2b_{\lambda_{i}}^{2}L_{g}^{2}+b_{\nabla
g_{i}}^{2})\|x_{i}-x_{i}^{\prime}\|^{2}$
$\displaystyle\leq\textstyle\sum_{i\in\mathcal{I}}{\max(2b_{\nabla
g_{i}}^{2},2b_{\lambda_{i}}^{2}L_{g}^{2}+b_{\nabla
g_{i}}^{2})}\|\bm{\omega}_{i}-\bm{\omega}_{i}^{\prime}\|^{2},$

where $\\{1\\}$ follows by adding and subtracting the term
$\nabla_{x_{i}}g_{i}(x_{i})^{\top}\lambda_{i}^{\prime}$ and by the bound
$\|a+b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$; $\\{2\\}$ is obtained by the Cauchy-
Schwartz inequality and by the fact that $g_{i}$ is Lipschitz since it has a
bounded gradient; $\\{3\\}$ is obtained by the Lipschitz continuity of
$\nabla_{x_{i}}g_{i}$. Hence, $\mathcal{C}_{1}$ is
$L_{\mathcal{C}_{1}}$-Lipschitz continuous, where
$L_{\mathcal{C}_{1}}=\max_{i\in\mathcal{I}}(\max(2b_{\nabla
g_{i}},\sqrt{2b_{\lambda_{i}}^{2}L_{g}^{2}+b_{\nabla g_{i}}^{2})})$. Since the
sum of Lipschitz continuous operators is Lipschitz continuous, the result
follows. ∎

######  Lemma 8.

Let $(\bm{x}^{\star},\bm{\lambda}^{\star})$ be a solution to the monotone
inclusion in (12). Then, $(\bm{x}^{\star},\bm{\lambda}^{\star})$ is also a
solution to the monotone inclusion in (8). $\hfill\Box$

###### Proof.

The proof follows that of [15, Thm. 2(i)]. ∎

##  Appendix C Results and Proofs of Section III

The following lemma shows the equivalence between
$\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$ and
$\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$.

######  Lemma 9.

Let Assumptions 1, 2, 7, and 8 hold. Furthermore, let
$\mathcal{T}_{\mathrm{FBF}}$ be defined by (18) while $\mathcal{A}$,
$\mathcal{B}$, and $\mathcal{C}$ be defined in (9)–(11). Then,
$\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})=\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$.
$\hfill\Box$

###### Proof.

The proof is analogous to that of [17, Prop. 1]. ∎

The following lemma is used to prove the quasi-shrinking property of the FBF
operator (18).

######  Lemma 10.

Let $\mathcal{A}$ and $\mathcal{B}$ maximally monotone and $\mathcal{B}$
continuous. Let

$\mathcal{T}={(\mathrm{Id}+\Psi^{-1}\mathcal{A})^{-1}}(\mathrm{Id}-\Psi^{-1}\mathcal{B}).$

Then $\mathrm{Id}-\mathcal{T}$ is demiclosed at $0$. $\hfill\Box$

###### Proof.

Let us consider a sequence $(v_{k})_{k\in\mathbb{N}}$ such that

$\lim_{k\xrightarrow{}\infty}v_{k}=v,\qquad\lim_{k\xrightarrow{}\infty}\mathcal{(}\mathrm{Id}-\mathcal{T})(v_{k})=0.$

We want to prove that $v-\mathcal{T}(v)=0$ or, equivalently,
$v\in\operatorname{fix}(\mathcal{T})$. Let us define
$u_{k}:=\mathcal{(}\mathrm{Id}-\mathcal{T})(v_{k}).$ Then,

$\displaystyle\begin{split}&v_{k}-u_{k}=(\mathrm{Id}+\Psi^{-1}\mathcal{A})^{-1}(\mathrm{Id}-{\Psi^{-1}}\mathcal{B})(v_{k})\\\
&\Leftrightarrow(\mathrm{Id}-{\Psi^{-1}}\mathcal{B})(v_{k})\in(\mathrm{Id}+{\Psi^{-1}}\mathcal{A})(v_{k}-u_{k})\\\
&\Leftrightarrow
v_{k}-{\Psi^{-1}}\mathcal{B}(v_{k})+u_{k}-v_{k}\in{\Psi^{-1}}\mathcal{A}(v_{k}-u_{k})\\\
&\Leftrightarrow\mathcal{-}\mathcal{B}(v_{k})+\Psi
u_{k}\in\mathcal{A}(v_{k}-u_{k}).\end{split}$

By the continuity of $\mathcal{B}$ and [27, Fact 1.19], we conclude that
$\lim_{k\xrightarrow{}\infty}-{\mathcal{B}(v_{k})}+\Psi
u_{k}=-\mathcal{B}(v).$ By [27, Prop. 20.37],
$\operatorname{gph}(\mathcal{A})$ is closed. Therefore, since
$\lim_{k\xrightarrow{}\infty}v_{k}-u_{k}=v$, we conclude that
$-\mathcal{B}(v)\in\mathcal{A}(v).$ By [27, Prop. 26.1(iv)], we obtain
$v\in\operatorname{fix}(\mathcal{T})$. ∎

###  C-A Proof of Lemma 3

By Lemmas 6 and 7, the operator $\mathcal{A}$ is maximally monotone whereas
the operator $\mathcal{B}+\mathcal{C}$ is maximally monotone and Lipschitz
continuous with Lipschitz constant denoted by $L_{B}$. Then, [17, Cor. 1]
shows that $\mathcal{T}_{\mathrm{FBF}}$ is quasi-nonexpansive when the step
size matrix $\Psi$, satisfy Assumption 8. Specifically, it holds that [17,
Prop. 2]:

$\|\mathcal{T}_{\mathrm{FBF}}(\bm{\omega})-\bm{\omega}^{\star}\|_{\Psi}^{2}\leq\|\bm{\omega}-\bm{\omega}^{\star}\|_{\Psi}^{2}-(L_{B}/\mu_{\min}(\Psi))^{2}\|\tilde{\bm{\omega}}-\bm{\omega}\|_{\Psi}^{2},$
(36)

where $\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$,
$\mu_{\min}(\Psi)$ is the smallest eigenvalue of $\Psi$ and
$\tilde{\bm{\omega}}=(\mathrm{Id}+\Psi^{-1}\mathcal{A})^{-1}(\mathrm{Id}-\Psi^{-1}(\mathcal{B}+\mathcal{C}))(\bm{\omega})$.
Finally, we prove that $\mathcal{T}_{\mathrm{FBF}}$ is quasi-shrinking by
invoking Lemma 1. Specifically, we choose
$\mathcal{T}_{2}=(\mathrm{Id}+\Psi^{-1}\mathcal{A})^{-1}(\mathrm{Id}-\Psi^{-1}(\mathcal{B}+\mathcal{C}))$.
As shown in the proof of [17, Prop. 1],
$\operatorname{fix}(\mathcal{T}_{2})=\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})=\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$.
Moreover, Lemma 10 shows that $\mathrm{Id}-\mathcal{T}_{2}$ is demiclosed at 0
and (36) is indeed the inequality in (3) for $\mathcal{T}_{\mathrm{FBF}}$.
$\hfill\blacksquare$

######  Remark 10.

Although [17, Cor. 1] shows quasi-nonexpansiveness of
$\mathcal{T}_{\mathrm{FBF}}$ and [17, Prop. 2] shows the inequality in (36)
for Problem (4) with a linear coupling constraint, these results also holds
for nonlinear functions $g_{i}(x_{i})$, for all $i\in\mathcal{I}$, as long as
Assumption 7 holds, since the operator $\mathcal{C}$ in (11) remains Lipschitz
continuous.

###  C-B Proof of Lemma 4

Firstly, we show that, for an arbitrary
$\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$,

$\displaystyle\|\mathcal{T}_{\mathrm{FBF}}(\bm{\omega})-\bm{\omega}^{\star}\|_{\Psi}^{2}<\|\bm{\omega}-\bm{\omega}^{\star}\|_{\Psi}^{2},$
(37)

for all $\bm{\omega}\notin\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$. To
this end, let us recall the inequality (36) in the proof of Lemma 3:

$\|\mathcal{T}_{\mathrm{FBF}}(\bm{\omega})-\bm{\omega}^{\star}\|_{\Psi}^{2}\leq\|\bm{\omega}-\bm{\omega}^{\star}\|_{\Psi}^{2}-(L_{B}/\mu_{\min}(\Psi))^{2}\|\tilde{\bm{\omega}}-\bm{\omega}\|_{\Psi}^{2},$

which holds for any
$\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$ and
$\bm{\omega}\in\operatorname{dom}(\mathcal{T}_{\mathrm{FBF}})$. Furthermore,
we consider any
$\bm{\omega}\notin\operatorname{fix}({\mathcal{T}_{\mathrm{FBF}}})$. Since
$\tilde{\bm{\omega}}=(\mathrm{Id}+{\Psi^{-1}\mathcal{A}})^{-1}(\mathrm{Id}-\Psi^{-1}(\mathcal{B}+\mathcal{C})))(\bm{\omega})$,
when $\tilde{\bm{\omega}}=\bm{\omega}$, it holds that

$\displaystyle\Psi(\tilde{\bm{\omega}}-\bm{\omega})-(\mathcal{B}+\mathcal{C})(\bm{\omega})+(\mathcal{B}+\mathcal{C})(\tilde{\bm{\omega}})$
$\displaystyle\in(\mathcal{A}+\mathcal{B}+\mathcal{C})(\tilde{\bm{\omega}})$
$\displaystyle\Leftrightarrow 0$
$\displaystyle\in(\mathcal{A}+\mathcal{B}+\mathcal{C})(\tilde{\bm{\omega}}),$

implying that
$\bm{\omega}=\tilde{\bm{\omega}}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$.
Hence, $\tilde{\bm{\omega}}\neq\bm{\omega}$ if
$\bm{\omega}\notin\operatorname{fix}({\mathcal{T}_{\mathrm{FBF}}})$. We
observe from the preceding inequality that when
$\tilde{\bm{\omega}}\neq\bm{\omega}$, the inequality (37) holds.

By using the inequality (37) and the fact that
$\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$ is bounded, we can then show
that for any arbitrary fixed point
$\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$, there
exists $R>0$ satisfying $\inf_{\|\bm{\omega}-\bm{\omega}^{\star}\|\geq
R}(\|\bm{\omega}-\bm{\omega}^{\star}\|-\|\mathcal{T}_{\mathrm{FBF}}(\bm{\omega})-\mathcal{T}_{\mathrm{FBF}}(\bm{\omega}^{\star})\|)>0$.
The proof of the previous statement is analogous to that of [41, Lem. 1],
which claims a similar inequality for attracting non-expansive operators.
Finally, we follow the proof of [41, Thm. 2], which claims the boundedness of
the HSDM sequence with an attracting non-expansive operator $\mathcal{T}$,
since, for some $R>0$, the inequality
$\inf_{\|\bm{\omega}-\bm{\omega}^{\star}\|\geq
R}(\|\bm{\omega}-\bm{\omega}^{\star}\|-\|\mathcal{T}_{\mathrm{FBF}}(\bm{\omega})-\mathcal{T}_{\mathrm{FBF}}(\bm{\omega}^{\star})\|)>0$
holds not only for attracting non-expansive operators but also
$\mathcal{T}_{\mathrm{FBF}}$, $\nabla\phi$ is monotone and Lipschitz
continuous (Assumption 3), and the step size $\beta^{(k)}$ is non-summable but
square summable (Assumption 4). $\hfill\blacksquare$

###  C-C Proof of Theorem 1

Let
$\tilde{\bm{\omega}}^{(k)}=(\tilde{\bm{x}}^{(k)},\tilde{\bm{\lambda}}^{(k)},\tilde{\bm{\nu}}^{(k)})$
and
$\accentset{\circ}{\bm{\omega}}^{(k)}=(\accentset{\circ}{\bm{x}}^{(k)},\accentset{\circ}{\bm{\lambda}}^{(k)},\accentset{\circ}{\bm{\nu}}^{(k)})$,
where
$\tilde{\bm{x}}^{(k)}=\operatorname{col}(\\{\tilde{x}_{i}\\}_{i\in\mathcal{I}})$
and the other variables are defined similarly. The updates of
$\tilde{\bm{\omega}}^{(k)}$ in Step 2 of Algorithm 1 can be compactly written
as

$\tilde{\bm{\omega}}^{(k)}=(\mathrm{Id}+{\Psi^{-1}\mathcal{A}})^{-1}(\mathrm{Id}-\Psi^{-1}(\mathcal{B}+\mathcal{C}))(\bm{\omega}^{(k)}),$

whereas the updates of $\accentset{\circ}{\bm{\omega}}^{(k)}$ in Step 4 of
Algorithm 1 can be compactly written as
$\accentset{\circ}{\bm{\omega}}^{(k)}=\tilde{\bm{\omega}}^{(k)}-\Psi^{-1}(\mathcal{B}+\mathcal{C})(\tilde{\bm{\omega}}^{(k)}-\bm{\omega}^{(k)}),$
implying that
$\accentset{\circ}{\bm{\omega}}^{(k)}=\mathcal{T}_{\mathrm{FBF}}(\bm{\omega}^{(k)})$
and the updates in (19) is compactly written as

$\bm{\omega}^{(k+1)}=\mathcal{T}_{\mathrm{FBF}}(\bm{\omega}^{(k)})-\beta^{(k)}\nabla\phi(\mathcal{T}_{\mathrm{FBF}}(\bm{\omega}^{(k)})),$
(38)

which is the HSDM applied to $\mathcal{T}_{\mathrm{FBF}}$.

We can then invoke Lemma 2 to claim the hypothesis. By Lemma 9,
$\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})=\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})$;
therefore $\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$ is non-empty and
bounded. Moreover, by Assumption 4, the step size $\beta^{(k)}$ meets the
conditions in Lemma 2. Lemma 3 shows that $\mathcal{T}_{\mathrm{FBF}}$ is
quasi-nonexpansive and quasi-shrinking on any bounded closed convex set, $C$
such that
$C\cap\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})\neq\varnothing$. On the
other hand, Lemma 4 shows that the FBF-HSDM sequence
$(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ obtained by the iterations in (38) is
bounded, i.e., for any
$\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$, there
exists a positive finite $R(\bm{\omega}^{\star})$ such that
$\|\bm{\omega}^{(k)}-\bm{\omega}^{\star}\|\leq R(\bm{\omega}^{\star})$.
Therefore, for an arbitrarily chosen
$\bm{\omega}^{\star}\in\operatorname{fix}(\mathcal{T}_{\mathrm{FBF}})$, we can
construct the following bounded closed set
$\mathfrak{B}(\bm{\omega}^{\star}):=\\{x\in\operatorname{dom}(\mathcal{T}_{\mathrm{FBF}})\mid\|x-\bm{\omega}^{\star}\|\leq
R(\bm{\omega}^{\star})\\},$ on which the sequence
$(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ lies. Moreover, we can observe that
indeed
$\mathfrak{B}\cap\operatorname{fix}({\mathcal{T}_{\mathrm{FBF}}})\neq\varnothing$,
since $\bm{\omega}^{\star}\in\mathfrak{B}$ is a fixed point of
$\mathcal{T}_{\mathrm{FBF}}$. Hence, $\mathcal{T}_{\mathrm{FBF}}$ is quasi-
shrinking on $\mathfrak{B}$, which completes the proof. $\hfill\blacksquare$

##  Appendix D Proofs of Section IV

###  D-A Proof of Theorem 2

First, we observe that in Algorithm 2,
$\accentset{\circ}{\bm{\omega}}^{(k)}=(\accentset{\circ}{\bm{x}}^{(k)},\accentset{\circ}{\bm{\lambda}}^{(k)},\accentset{\circ}{\bm{\nu}}^{(k)})$
is updated by using $\mathcal{T}_{\mathrm{pFB}}$ in (23), i.e.,
$\accentset{\circ}{\bm{\omega}}^{(k)}=\mathcal{T}_{\mathrm{pFB}}(\bm{\omega}^{(k)})$
[15, Section 4, Algorithm 1]. Hence, we can see that $\bm{\omega}^{(k)}$ is
updated via the HSDM method, i.e.,

$\bm{\omega}^{(k+1)}=\mathcal{T}_{\mathrm{pFB}}(\bm{\omega}^{(k)})-\beta^{(k)}\nabla\phi(\mathcal{T}_{\mathrm{pFB}}(\bm{\omega}^{(k)})),$
(39)

Similarly to the proof of Theorem 1, due to the boundedness of
$\operatorname{fix}(\mathcal{T}_{\mathrm{pFB}})=\operatorname{zer}(\mathcal{A}+\mathcal{B}+\mathcal{C})\neq\varnothing$
and the step size rule of $\beta^{(k)}$ in Assumption 4, we can invoke Lemma
2. Specifically, the operator $\mathcal{T}_{\mathrm{pFB}}$ is averaged
nonexpansive when Assumptions 1, 2, 7, and 9–11 hold [15, Thm. 3]. Therefore,
$\mathcal{T}_{\mathrm{pFB}}$ is also quasi-nonexpansive [27, Section 4.1]. By
[27, Prop. 4.35 (iii)], the condition in (3) holds with
$\mathcal{T}_{2}=\mathcal{T}$. By [27, Thm. 4.27],
$\mathrm{Id}-\mathcal{T}_{\mathrm{pFB}}$ is demiclosed at 0. Therefore, by
Lemma 1, $\mathcal{T}_{\mathrm{pFB}}$ is quasi-shrinking on any closed bounded
convex set whose intersection with
$\operatorname{fix}(\mathcal{T}_{\mathrm{pFB}})$ is nonempty. Furthermore,
since $\mathcal{T}_{\mathrm{pFB}}$ is averaged nonexpansive,
$\mathcal{T}_{\mathrm{pFB}}$ is attracting. Therefore, by [41, Thm. 2] and due
to the choice of the step size $\beta^{(k)}$ in Assumption 4, the sequence
generated by (39) is bounded. $\hfill\blacksquare$

##  Appendix E Proofs of Section V

###  E-A Preliminary results

First, we show a series of preliminary results in Lemmas 11–13 that lead to
the proofs of Lemma 5 and Theorem 3. The proofs of this section are provided
in the standard Euclidean norm for ease of notation. However, the case for any
$\Psi$-induced norm, with $\Psi\succ 0$, follows verbatim. First, Lemma 11
shows the convergence of a particular sequence and can be regarded as a
finite-iteration version of [28, Lem. 1].

######  Lemma 11.

Let $\psi:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ be non-decreasing and
non-negative. Let a sequence $(b^{(k)})_{k\in\mathbb{N}}$ be non-increasing,
non-negative. Let $(a^{(k)})_{k\in\mathbb{N}}\subset[0,\infty)$ satisfy

$a^{(k+1)}\leq a^{(k)}-\psi(a^{(k)})+b^{(k+1)}.$ (40)

Let $K\in\mathbb{N}$. If there exists $\xi>0$ such that
$\psi(\xi)\geq\max\\{2b^{(1)},\frac{2}{K-1}a^{(1)}\\}$, then

$a^{(k)}\leq\xi+b^{(k)},~{}~{}~{}~{}\forall~{}k\geq K.$ (41)

###### Proof.

Let us first show that there exists an $M\in\mathbb{N}$, $M\leq K$ such that
$a^{(M)}\leq\xi$. We proceed by contradiction, assuming that $a^{(k)}>\xi$
$\forall k=1,...,K$. Then, by noting that $\psi(\cdot)$ is non-decreasing and
that $\psi(\xi)\geq 2b^{(k)}$ for all $k\in\mathbb{N}$, we have

$\displaystyle\begin{split}a^{(k+1)}&\leq a^{(k)}-\psi(a^{(k)})+b^{(k+1)}\\\
&\leq
a^{(k)}-\psi(\xi)+\tfrac{1}{2}\psi(\xi)=a^{(k)}-\tfrac{1}{2}\psi(\xi).\end{split}$

By iterating the latter relation and recalling that
$\psi(\xi)\geq\frac{2}{K-1}a^{(1)}$, we find that

$\displaystyle\begin{split}&a^{(k+1)}\leq a^{(1)}-\tfrac{k}{2}\psi(\xi)\leq
a^{(1)}-\tfrac{k}{K-1}a^{(1)}.\end{split}$

For $k=K$, we then obtain the contradiction $a^{(K+1)}<0$. Thus, there exists
$M\leq K$ such that $a^{(M)}\leq\xi$. We then proceed by induction to prove
(41). Let us prove that, if $a^{(k)}\leq\xi+b^{(k)}$ then
$a^{(k+1)}\leq\xi+b^{(k+1)}$ for all $k\geq M$. We distinguish two cases:
1) Case $a^{(k)}<\xi$. Then, by (40) and by the non-negativity of
$\psi(\cdot)$, $a^{(k+1)}\leq a^{(k)}+b^{(k+1)}<\xi+b^{(k+1)}.$
2) Case $\xi\leq a^{(k)}\leq\xi+b^{(k)}$. Then, by the non-decreasing property
of $\psi$, $a^{(k)}\geq\xi\Rightarrow\psi(a^{(k)})\geq\psi(\xi)$. By the
assumptions, $\psi(\xi)\geq 2b^{(1)}$ and by the non-incresing property of
$(b^{k})_{k\in\mathbb{N}}$, $2b^{(1)}\geq b^{(k)}+b^{(k+1)}$. We thus obtain
$\psi(a^{(k)})\geq b^{(k)}+b^{(k+1)}$. Substituting into (40) leads to

$\displaystyle a^{(k+1)}$ $\displaystyle\leq a^{(k)}-\psi(a^{(k)})+b^{(k+1)}$
$\displaystyle\leq
a^{(k)}-b^{(k)}-b^{(k+1)}+b^{(k+1)}=a^{(k)}-b^{(k)}\leq\xi.$

We conclude by induction that $a^{(k)}\leq\xi+b^{(k)}$ for all $k\geq M$ and,
since $M\leq K$, the claim in (41) immediately follows. ∎

######  Lemma 12.

Let $\mathcal{T}$ quasi-nonexpansive and $\mathcal{F}$ strongly monotone, such
that $\|\mathcal{F}(\bm{\omega})\|\leq U$ for all
$\bm{\omega}\in\mathrm{im}(\mathcal{T})$. Let
$(\bm{\omega}^{(k)})_{k\in\mathbb{N}}$ be generated from (17) with constant
stepsize $\beta^{(k)}=\beta>0$ for all $k$. Let $K\in\mathbb{N}$ and let
$\bm{\omega}^{\star}$ be the solution of
$\mathrm{VI}(\mathcal{F},\operatorname{fix}(\mathcal{T}))$. If there exists
$\xi$ such that the shrinkage function $D(\cdot)$ of $\mathcal{T}$, defined in
(2), satisfies $D(\xi)\geq\max\\{2\beta
U,2\frac{\mathrm{dist}(\omega^{(1)},\operatorname{fix}(\mathcal{T}))}{K-1}\\}$,
then the following inequalities hold:

$\displaystyle\sup_{k\geq
K}\mathrm{dist}(\bm{\omega}^{(k)},\operatorname{fix}(\mathcal{T}))$
$\displaystyle\leq\xi+\beta U,$ (42) $\displaystyle\sup_{k\geq
K}\|\mathcal{T}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)}\|$ $\displaystyle\leq
2(\xi+\beta U),$ (43) $\displaystyle\sup_{k\geq
K}\langle\mathcal{T}(\bm{\omega}^{(k)})-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle$
$\displaystyle\leq 3(\xi+\beta U)\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|.$
(44)

_Proof._ (i.) For all $k$, it holds by the definition of distance and by the
algorithm definition in (17) that:

$\displaystyle\mathrm{dist}(\bm{\omega}^{(k+1)},\operatorname{fix}(\mathcal{T}))\leq\|\bm{\omega}^{(k+1)}-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\mathcal{T}(\bm{\omega}^{(k)}))\|=$
$\displaystyle\|\mathcal{T}(\bm{\omega}^{(k)})-\beta\mathcal{F}(\mathcal{T}(\bm{\omega}^{(k)}))-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\mathcal{T}(\bm{\omega}^{(k)}))\|\leq$
$\displaystyle\underbrace{\|\mathcal{T}(\bm{\omega}^{(k)})-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\mathcal{T}(\bm{\omega}^{(k)}))\|}_{=\mathrm{dist}(\mathcal{T}(\bm{\omega}^{(k)}),\operatorname{fix}(\mathcal{T}))}+\beta\|\mathcal{F}(\mathcal{T}(\bm{\omega}^{(k)}))\|\leq$
$\displaystyle\mathrm{dist}(\mathcal{T}(\bm{\omega}^{(k)}),\operatorname{fix}(\mathcal{T}))+\beta
U.$ (45)

Let us define
$a^{(k)}:=\mathrm{dist}(\bm{\omega}^{(k)},\operatorname{fix}(\mathcal{T}))$.
Then, from (45) we find immediately $a^{(k+1)}-\beta
U\leq\mathrm{dist}(\mathcal{T}(\bm{\omega}^{(k)}),\operatorname{fix}(\mathcal{T}))$.
By the definition of shrinkage function in (2) and the latter inequality, we
can write

$\displaystyle\begin{split}&D(a^{(k)})\leq
a^{(k)}-\mathrm{dist}(\mathcal{T}(\bm{\omega}^{(k)}),\operatorname{fix}(\mathcal{T}))\leq\\\
&a^{(k)}-a^{(k+1)}+\beta U\Rightarrow a^{(k+1)}{\leq}a^{(k)}+\beta
U-D(a^{(k)}),\end{split}$

which defines a sequence of the kind in (40) with $\psi(\cdot)=D(\cdot)$ and
$b^{(k)}=\beta U$ for all $k$. By Lemma 11, then
$\mathrm{dist}(\bm{\omega}^{(k)},\operatorname{fix}(\mathcal{T}))\leq\xi+\beta
U$ for all $k\geq K$.
(ii.) By the triangle inequality, we can write
$\|\mathcal{T}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)}\|\leq{\|\mathcal{T}(\bm{\omega}^{(k)})-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)})\|}+\|\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)}\|.$
By quasi-nonexpansiveness of $\mathcal{T}$, we obtain, for all $k\geq K$,

$\displaystyle\begin{split}&{\|\mathcal{T}(\bm{\omega}^{(k)})-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)})\|}\leq\\\
&\|\bm{\omega}^{(k)}-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)})\|=\mathrm{dist}(\bm{\omega}^{(k)},\operatorname{fix}(\mathcal{T}))\\\
&\Rightarrow\|\mathcal{T}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)}\|\leq
2\mathrm{dist}(\bm{\omega}^{(k)},\operatorname{fix}(\mathcal{T})).\end{split}$

Finally, combining the last inequality and (42) yields (43).
(iii) By the Cauchy-Schwarz inequality, we can write

$\displaystyle\langle\mathcal{T}(\bm{\omega}^{(k)})-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle=$
$\displaystyle\langle\mathcal{T}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle+\langle\bm{\omega}^{(k)}-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle\leq$
$\displaystyle\|\mathcal{T}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)}\|\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|+\langle\bm{\omega}^{(k)}-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle.$
(46)

Based on (43), for all $k\geq K$, we can bound the first term on the right-
hand side of (46) by
$\|\mathcal{T}(\bm{\omega}^{(k)})-\bm{\omega}^{(k)}\|\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|\leq
2(\xi+\beta U)\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|$ and rewrite the
second term as

$\displaystyle\begin{split}\langle\bm{\omega}^{(k)}-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle=&{\langle\bm{\omega}^{(k)}-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)}),-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle}\\\
+&{\langle\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)})-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle}.\end{split}$

We observe that the second addend is non-positive by the definition of
$\mathrm{VI}$ solution. By applying the Cauchy-Schwarz inequality, the
definition of projection, and (42), we obtain

$\displaystyle\begin{split}&\langle\mathcal{T}(\bm{\omega}^{(k)})-\bm{\bm{\omega}^{\star}},-\mathcal{F}(\bm{\bm{\omega}^{\star}})\rangle\leq\\\
&2(\xi+\beta
U)\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|+\|\bm{\omega}^{(k)}-\operatorname{proj}_{\operatorname{fix}(\mathcal{T})}(\bm{\omega}^{(k)})\|\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|=\\\
&2(\xi+\lambda
U)\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|+\mathrm{dist}(\bm{\omega}^{(k)},\operatorname{fix}(\mathcal{T}))\|{\mathcal{F}}(\bm{\bm{\omega}^{\star}})\|\leq\\\
&3(\xi+\lambda U)\|\mathcal{F}(\bm{\bm{\omega}^{\star}})\|.\end{split}$

$\hfill\blacksquare$

######  Lemma 13.

Let Assumptions 12–15 hold. For any $t\in\mathbb{N}$, let $\bm{\omega}_{t+1}$
be generated from the step at time $t$ of the restarted HSDM algorithm in
(28). Let $D_{t}(\cdot)$ be the shrinkage function of $\mathcal{T}_{t}$ as
defined in (2). If there exists $\xi>0$ such that

$D_{t}(\xi)\geq\max{\left\\{2\beta
U,2\tfrac{\mathrm{dist}(\bm{\omega}_{t},\operatorname{fix}(\mathcal{T}_{t}))}{K-1}\right\\}},$
(47)

then, the bound in (29), i.e.,

$\|\bm{\omega}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\left(1-{\tau(\beta)}\right)^{K}\|\bm{\omega}_{t}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+\gamma,$

holds with

$\gamma=\tfrac{\beta}{\tau(\beta)}U(6\xi+11\beta U).$ (48)

###### Proof.

Let us define the operator
$\mathcal{T}^{\beta}_{t}(\bm{\omega}):=\mathcal{T}_{t}(\bm{\omega})-\beta\nabla\phi_{t}(\mathcal{T}_{t}(\bm{\omega})).$
By
$\mathcal{T}_{t}(\bm{\bm{\omega}^{\star}}_{t})=\bm{\bm{\omega}^{\star}}_{t}$
and by the definition of the algorithm in (28),
$\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}=\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}$.
We sum and subtract $\beta\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})$ and
substitute $\mathcal{T}^{\beta}_{t}$ to obtain

$\displaystyle\begin{split}&\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\\\
&=\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}_{t}(\bm{\bm{\omega}^{\star}}_{t})+\beta\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})-\beta\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}\\\
&=\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})-\beta\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}.\end{split}$

Expanding the square ${\\{1\\}}$, expanding $\mathcal{T}^{\beta}_{t}$
$\\{2\\}$, and regrouping $\\{3\\}$ leads to

$\displaystyle\begin{split}&\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\\\
&\overset{\\{1\\}}{=}\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+\beta^{2}\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}\\\
&+2\langle\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t}),-\beta\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\rangle\\\
&\overset{\\{2\\}}{=}\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+\beta^{2}\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}-\\\
&2\beta\langle\mathcal{T}_{t}(\bm{y}^{(K)})-\beta\nabla\phi_{t}(\mathcal{T}_{t}(\bm{y}^{(K)}))-\\\
&\mathcal{T}_{t}(\bm{\bm{\omega}^{\star}}_{t})+\beta\nabla\phi_{t}(\mathcal{T}_{t}(\bm{\bm{\omega}^{\star}}_{t})),\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\rangle\\\
&\overset{\\{3\\}}{=}\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+\beta^{2}\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}\\\
&+2\beta\langle\mathcal{T}_{t}(\bm{y}^{(K)})-\bm{\bm{\omega}^{\star}}_{t},-\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\rangle\\\
&+2\beta^{2}\langle\nabla\phi_{t}(\mathcal{T}_{t}(\bm{y}^{(K)}))-\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t}),\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\rangle.\end{split}$
(49)

We note that, by applying the Cauchy-Schwarz, the triangle inequalities and
Assumption 15, we have
$\langle\nabla\phi_{t}(\mathcal{T}_{t}(\bm{y}^{(K)}))-\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t}),\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\rangle\leq\|\nabla\phi_{t}(\mathcal{T}_{t}(\bm{y}^{(K)}))-\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|\leq(U+\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|)\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|.$
By (47), the bounds in Lemma 12 hold. We then substitute in (49), the latter
relation, and the bound in (44) to obtain

$\displaystyle\begin{split}&\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+\\\
&6\beta(\xi+\beta
U)\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|+\beta^{2}(2U+3\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|)\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|.\end{split}$

Applying Assumption 15 and rearranging the terms leads to

$\displaystyle\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+$
$\displaystyle\quad\quad 6\beta(\xi+\beta
U)\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|+\beta^{2}5U\|\nabla\phi_{t}(\bm{\bm{\omega}^{\star}}_{t})\|$
$\displaystyle\leq\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+\beta(6\xi+11\beta
U)U$
$\displaystyle\leq\|\mathcal{T}^{\beta}_{t}(\bm{y}^{(K)})-\mathcal{T}^{\beta}_{t}(\bm{\bm{\omega}^{\star}}_{t})\|^{2}+{\tau(\beta)}\gamma.$
(50)

By quasi-nonexpansiveness of $\mathcal{T}_{t}$ as well as strong monotonicity
and Lipschitz continuity of $\nabla\phi_{t}$, we can apply [28, Lem. 4a] to
obtain
$\|\mathcal{T}^{\beta}_{t}(\bm{\omega})-\mathcal{T}^{\beta}_{t}(\bar{\bm{\omega}})\|\leq(1-\tau(\beta))\|\bm{\omega}-\bar{\bm{\omega}}\|,$
for all
$\bm{\omega}\in\mathrm{dom}(\mathcal{T}^{\beta}_{t}),\bar{\bm{\omega}}\in\operatorname{fix}(\mathcal{T}_{t})$,
which we substitute in (50) to obtain

$\displaystyle\begin{split}&\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\left(1-{\tau(\beta)}\right)^{2}\|\bm{y}^{(K)}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+{\tau(\beta)}\gamma\\\
&\leq\left(1-{\tau(\beta)}\right)\|\bm{y}^{(K)}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+{\tau(\beta)}\gamma.\end{split}$

By iterating, we obtain

$\displaystyle\begin{split}&\|\bm{\bm{\omega}}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\\\
&\left(1-{\tau(\beta)}\right)^{2}\|\bm{y}^{(K-1)}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+\left(1-{\tau(\beta)}\right){\tau(\beta)}\gamma+{\tau(\beta)}\gamma\\\
&\leq\dots\leq(1-{\tau(\beta)})^{K}\|\bm{y}^{(1)}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+\sum_{j=0}^{K-1}\left(1-{\tau(\beta)}\right)^{j}{\tau(\beta)}\gamma\\\
&\leq\left(1-{\tau(\beta)}\right)^{K}\|\bm{y}^{(1)}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}+\sum_{j=0}^{\infty}\left(1-{\tau(\beta)}\right)^{j}{\tau(\beta)}\gamma.\end{split}$

Applying the geometric series convergence and recalling from (28) that
$\bm{y}^{(1)}=\bm{\omega}_{t}$ leads to (29). ∎

###  E-B Proof of Lemma 5

Let us consider $\xi:=\frac{\gamma\sigma}{12U}.$ Since $\mathcal{T}_{t}$ is
quasi-shrinking, the shrinkage function $D_{t}$ of $\mathcal{T}_{t}$ satisfies
$D_{t}(\xi)>0$. Thus, there exist
$\bar{\beta}\in(0,\frac{2\sigma}{L_{\phi}^{2}})$ and $K$ such that, for any
$\beta\in(0,\bar{\beta}]$,

$D_{t}(\xi)\geq\max{\left\\{2\beta
U,2\tfrac{\mathrm{dist}(\bm{\omega}_{t},\operatorname{fix}(\mathcal{T}_{t}))}{K-1}\right\\}}.$
(51)

######  Remark 11.

As $D_{t}(\xi)$ decreases with $\gamma$, for smaller values of $\gamma$ a
smaller stepsize $\beta$ and a larger $K$ are necessary.

It can be verified that $\lim_{\beta\rightarrow
0^{+}}\tfrac{\beta}{\tau(\beta)}=\tfrac{1}{\sigma}$. Then,

$\displaystyle\begin{split}\lim_{\beta\rightarrow
0^{+}}\tfrac{\beta}{\tau(\beta)}(6\xi+11\beta U)U&=\tfrac{6\xi
U}{\sigma}\leq\tfrac{1}{2}\gamma,\end{split}$ (52)

We thus find $\beta\in(0,\bar{\beta}]$ small enough, such that

$\tfrac{\beta}{\tau(\beta)}(6\xi+11\beta U)U\leq\gamma.$ (53)

Hence, the hypothesis holds by invoking Lemma 13. $\hfill\blacksquare$

###  E-C Proof of Theorem 3

We begin the proof by constructing a suitable stepsize $\bar{\beta}$ and
number of iterations $\bar{K}$. We then proceed with proving that the
statement holds for the chosen variables. Let us first define the auxiliary
variable $\xi=\frac{\gamma\sigma}{12U}.$ By (52), we can choose a small enough
$\bar{\beta}\in(0,\min\\{\frac{2\sigma}{L^{2}_{\varphi}},\frac{D(\xi)}{2U}\\})$,
such that

$\tfrac{\bar{\beta}}{\tau(\bar{\beta})}(6\xi+11\bar{\beta}U)U\leq\gamma.$ (54)

We now define $\alpha(K):=(1-\tau(\bar{\beta}))^{K}.$ Since
$\tau(\bar{\beta})\in(0,1)$, $\alpha$ is decreasing with $K$. We can then
choose $K_{1}$, such that $\alpha(K_{1})<\frac{1}{2}$. Then, we define the
mapping ${a}:\mathbb{N}_{\geq K_{1}}\to\mathbb{R}$

${a}(K)=\max\left\\{\|\bm{\omega}_{1}\|+\sup_{\bm{\omega}\in\mathcal{Y}}\|\bm{\omega}\|,\sqrt{\tfrac{2\alpha(K)\delta_{1}^{2}+\gamma}{1-2\alpha(K)}}\right\\},$
(55)

We can verify that ${a(\cdot)}$ is non-increasing. Consequently, the sequence
$\left(\frac{2({a}(K)+\delta_{2})}{{K}-1}\right)_{K\geq K_{1}}$ is decreasing.
We can then choose any sufficiently large $\bar{K}\geq K_{1}$, such that

$D(\xi)\geq\tfrac{2(\bar{a}+\delta_{2})}{\bar{K}-1},$ (56)

where $\bar{a}:=a(\bar{K})$. We also define $\bar{\alpha}:=\alpha(\bar{K})$.

We now prove by induction that

$\|\bm{\omega}_{t}-\bm{\omega}_{t-1}^{\star}\|\leq\bar{a}~{}~{}~{}~{}\text{for
all}~{}t>1.$ (57)

To that end, we first show that

$\|\bm{\omega}_{t}-\bm{\omega}_{t-1}^{\star}\|\leq\bar{a}\Rightarrow\|\bm{\omega}_{t+1}-\bm{\omega}_{t}^{\star}\|\leq\bar{a}.$
(58)

Let us then write

$\displaystyle\mathrm{dist}(\bm{\omega}_{t},$
$\displaystyle\operatorname{fix}(\mathcal{T}_{t}))\overset{\\{1\\}}{\leq}\|\bm{\omega}_{t}-\operatorname{proj}_{\operatorname{fix}{(\mathcal{T}_{t})}}(\bm{\omega}^{\star}_{t-1})\|$
$\displaystyle\overset{\\{2\\}}{\leq}\|\bm{\omega}_{t}-\bm{\omega}^{\star}_{t-1}\|+\|\bm{\omega}^{\star}_{t-1}-\operatorname{proj}_{\operatorname{fix}{(\mathcal{T}_{t})}}(\bm{\omega}^{\star}_{t-1})\|$
$\displaystyle\overset{\\{3\\}}{\leq}\|\bm{\omega}_{t}-\bm{\omega}^{\star}_{t-1}\|+\delta_{2}\leq\bar{a}+\delta_{2},$
(59)

where $\\{1\\}$ follows from the definition of distance, $\\{2\\}$ from the
triangle inequality and $\\{3\\}$ from Assumption 17.ii. Then, by Assumption
16, by the choice $\bar{\beta}\leq\frac{D(\xi)}{2U}$ and (56),

$\displaystyle\begin{split}D_{t}({\xi})\geq\max\left\\{2\bar{\beta}U,\tfrac{2(\bar{a}+\delta_{2})}{\bar{K}-1}\right\\}\geq\max\left\\{2\bar{\beta}U,\tfrac{2\mathrm{dist}(\bm{\omega}_{t},\operatorname{fix}(\mathcal{T}_{t}))}{\bar{K}-1}\right\\}.\end{split}$
(60)

By Lemma 13 and (54), we then have

$\|\bm{\omega}_{t+1}-\bm{\omega}_{t}^{\star}\|^{2}\leq\bar{\alpha}\|\bm{\omega}_{t}-\bm{\omega}_{t}^{\star}\|^{2}+\gamma.$
(61)

Applying on (61) the triangle inequality, the fact $(a+b)^{2}\leq
2a^{2}+2b^{2}$ and Assumption 17.i leads to

$\displaystyle\|\bm{\omega}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}$
$\displaystyle\leq
2\alpha(\|\bm{\omega}_{t}-\bm{\bm{\omega}^{\star}}_{t-1}\|^{2}+\|\bm{\omega}^{\star}_{t-1}-\bm{\omega}^{\star}_{t}\|^{2})+\gamma$
$\displaystyle\leq
2\bar{\alpha}(\|\bm{\omega}_{t}-\bm{\bm{\omega}^{\star}}_{t-1}\|^{2}+\delta_{1}^{2})+\gamma$
$\displaystyle\leq 2\bar{\alpha}(\bar{a}^{2}+\delta_{1}^{2})+\gamma.$ (62)

Finally, by (55), it holds that

$2\bar{\alpha}(\bar{a}^{2}+\delta_{1}^{2})+\gamma\leq\bar{a}^{2}\Leftrightarrow\bar{a}^{2}\geq\tfrac{2\bar{\alpha}\delta_{1}^{2}+\gamma}{1-2\bar{\alpha}}.$
(63)

Thus, we obtain
$\|\bm{\omega}_{t+1}-\bm{\bm{\omega}^{\star}}_{t}\|^{2}\leq\bar{a}^{2}.$ We
now continue the induction argument by proving

$\|\bm{\omega}_{2}-\bm{\bm{\omega}^{\star}}_{1}\|^{2}\leq\bar{a}^{2}.$ (64)

From the triangle inequality and from (55),
$\|\bm{\omega}_{1}-\bm{\bm{\omega}^{\star}}_{1}\|\leq\|\bm{\omega}_{1}\|+\|\bm{\bm{\omega}^{\star}}_{1}\|\leq\bar{a}.$
From the definition of distance, we obtain

$\mathrm{dist}(\bm{\omega}_{1},\operatorname{fix}(\mathcal{T}_{1}))\leq\|\bm{\omega}_{1}-\bm{\bm{\omega}^{\star}}_{1}\|\leq\bar{a}\leq\bar{a}+\delta_{2}.$
(65)

Then, $D_{t}({\xi})\geq D({\xi})\geq\max\left\\{2\beta
U,\frac{2(\bar{a}+\delta_{2})}{\bar{K}-1}\right\\}\geq\max\left\\{2\beta
U,\frac{2\mathrm{dist}(\bm{\omega}_{1},\operatorname{fix}(\mathcal{T}_{1}))}{\bar{K}-1}\right\\}.$
By Lemma 13 and (54), we find

$\displaystyle\begin{split}&\|\bm{\omega}_{2}-\bm{\omega}_{1}^{\star}\|^{2}\leq\bar{\alpha}\|\bm{\omega}_{1}-\bm{\omega}_{1}^{\star}\|^{2}+\gamma.\end{split}$

By using (65) and (63) to upperbound the right hand side of the last
inequality, we then obtain

$\|\bm{\omega}_{2}-\bm{\omega}_{1}^{\star}\|^{2}\leq\bar{\alpha}\bar{a}^{2}+\gamma\leq\bar{\alpha}(2\bar{a}^{2}+2\delta_{1}^{2})+\gamma\leq\bar{a}^{2}.$

Therefore, combining (58) and (64) leads to
$\sup_{t>1}\|\bm{\omega}_{t}-\bm{\bm{\omega}^{\star}}_{t-1}\|\leq\bar{a}.$
Recalling that, from Assumption 14, $\bm{\omega}^{\star}_{t}\in\mathcal{Y}$
for all $t$, this immediately implies
$\mathrm{dist}(\bm{\omega}_{t},\mathcal{Y})\leq\bar{a}~{}~{}~{}\text{for
all}~{}t>1,$ which proves that the sequence is bounded.
We now proceed with proving (30). We note that the relation in (61) holds for
all $t$. We then observe that, by the triangle inequality, by $(a+b)^{2}\leq
2a+2b$, and by Assumption 17,

$\displaystyle\begin{split}\|\bm{\omega}_{t+1}-\bm{\omega}^{\star}_{t+1}\|^{2}&\leq
2\|\bm{\omega}_{t+1}-\bm{\omega}^{\star}_{t}\|^{2}+2\|\bm{\omega}^{\star}_{t+1}-\bm{\omega}^{\star}_{t}\|^{2}\\\
&\leq
2\|\bm{\omega}_{t+1}-\bm{\omega}^{\star}_{t}\|^{2}+2\delta_{1}^{2}.\end{split}$

By using (61) to upper bound
$\|\bm{\omega}_{t+1}-\bm{\omega}^{\star}_{t}\|^{2}$ and iterating, we find:

$\displaystyle\begin{split}&\|\bm{\omega}_{t+1}-\bm{\omega}^{\star}_{t+1}\|^{2}\leq
2\bar{\alpha}\|\bm{\omega}_{t}-\bm{\omega}^{\star}_{t}\|^{2}+2(\gamma+\delta_{1}^{2})\\\
&\leq(2\bar{\alpha})^{2}\|\bm{\omega}_{t-1}-\bm{\omega}^{\star}_{t-1}\|^{2}+2(\gamma+\delta_{1}^{2})+2\bar{\alpha}(2\gamma+2\delta_{1}^{2})\\\
&\leq\dots\leq(2\bar{\alpha})^{t}\|\bm{\omega}_{1}-\bm{\omega}^{\star}_{1}\|^{2}+\textstyle\sum_{j=0}^{t-1}(2\bar{\alpha})^{j}(2\gamma+2\delta_{1}^{2}).\end{split}$

By taking the limit for $t\rightarrow\infty$ and by applying the convergence
of the geometric sequence, we obtain (30). $\hfill\blacksquare$

###  E-D Proof of Corollary 1

Steps i–vi of Algorithm 3 are analogous to Steps 1–6 of Algorithm 1.
Analogously to the proof of Theorem 1, we see that the variable
$\bm{y}^{(k)}:=(\hat{x}_{i}^{(k)},\hat{\lambda}_{i}^{(k)},\hat{\nu}_{i}^{(k)})$
is updated at each time step by $K$ iterations of the HSDM:

$\bm{y}^{(k+1)}=\mathcal{T}_{\mathrm{FBF},t}(\bm{y}^{(k)})-\beta\nabla\phi_{t}(\mathcal{T}_{\mathrm{FBF},t}(\bm{y}^{(k)})),~{}~{}k=1,...,K.$

Then, the variable $\bm{\omega}_{t}$ is updated as
$\bm{\omega}_{t}=\bm{y}^{(K+1)}.$ Thus, we see that Algorithm 3 is a
particular instance of the restarted HSDM algorithm (28). By Theorem 3,
$\bm{\omega}_{t}$ is bounded, therefore there exists a compact set
$\mathcal{Z}$ such that $(\bm{\omega}_{t})\in\mathcal{Z}$ for all $t$. By
Lemma 3, $\mathcal{T}_{\mathrm{FBF},t}$ is quasi-nonexpansive and quasi-
shrinking on any bounded, closed convex set $C$ such that
$C\cap\operatorname{fix}(\mathcal{T}_{\mathrm{FBF},t})\neq\emptyset$. In
particular, it is quasi-shrinking on any convex set
$C\supset\mathcal{Y}\cup\mathcal{Z}$, where $\mathcal{Y}$ is a compact set
such that $\bm{\omega}_{t}^{\star}\in\mathcal{Y},\forall t\in\mathbb{N}$
(Assumption 14). We then find Assumption 13 to hold and, by Theorem 3, the
tracking error is given by (30). $\hfill\blacksquare$

## References

  * [1] I. Shilov, H. Le Cadre, and A. Busic, “Privacy impact on generalized Nash equilibrium in peer-to-peer electricity market,” Operations Research Letters, vol. 49, no. 5, pp. 759–766, 2021.
  * [2] G. Belgioioso, W. Ananduta, S. Grammatico, and C. Ocampo-Martinez, “Operationally-safe peer-to-peer energy trading in distribution grids: A game-theoretic market-clearing mechanism,” IEEE Transactions on Smart Grid, 2022.  accepted.
  * [3] B. G. Bakhshayesh and H. Kebriaei, “Decentralized equilibrium seeking of joint routing and destination planning of electric vehicles: A constrained aggregative game approach,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–10, 2021.  Early access at https://doi.org/10.1109/TITS.2021.3123207.
  * [4] J.-S. Pang, G. Scutari, F. Facchinei, and C. Wang, “Distributed power allocation with rate constraints in Gaussian parallel interference channels,” IEEE Transactions on Information Theory, vol. 54, no. 8, pp. 3471–3489, 2008.
  * [5] G. Scutari, D. P. Palomar, F. Facchinei, and J.-S. Pang, “Monotone games for cognitive radio systems,” in Distributed Decision Making and Control, pp. 83–112, Springer London, 2012.
  * [6] J. Wang, M. Peng, S. Jin, and C. Zhao, “A generalized Nash equilibrium approach for robust cognitive radio networks via generalized variational inequalities,” IEEE Transactions on Wireless Communications, vol. 13, no. 7, pp. 3701–3714, 2014.
  * [7] F. Facchinei and C. Kanzow, “Generalized Nash equilibrium problems,” Annals of Operations Research, vol. 175, no. 1, pp. 177–211, 2010.
  * [8] 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.
  * [9] 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.
  * [10] D. Paccagnan, B. Gentile, F. Parise, M. Kamgarpour, and J. Lygeros, “Nash and Wardrop equilibria in aggregative games with coupling constraints,” IEEE Transactions on Automatic Control, vol. 64, no. 4, pp. 1373–1388, 2018\.
  * [11] G. Belgioioso and S. Grammatico, “Projected-gradient algorithms for generalized equilibrium seeking in aggregative games are preconditioned forward-backward methods,” in Proceedings of the 2018 European Control Conference (ECC), pp. 2188–2193, 2018.
  * [12] H. Yin, U. V. Shanbhag, and P. G. Mehta, “Nash equilibrium problems with scaled congestion costs and shared constraints,” IEEE Transactions on Automatic Control, vol. 56, no. 7, pp. 1702–1708, 2011.
  * [13] G. Belgioioso and S. Grammatico, “Semi-decentralized Nash equilibrium seeking in aggregative games with separable coupling constraints and non-differentiable cost functions,” IEEE Control Systems Letters, vol. 1, no. 2, pp. 400–405, 2017.
  * [14] G. Belgioioso and S. Grammatico, “Semi-decentralized generalized Nash equilibrium seeking in monotone aggregative games,” IEEE Transactions on Automatic Control, 2021.  Early access at https://doi.org/10.1109/TAC.2021.3135360.
  * [15] P. Yi and L. Pavel, “An operator splitting approach for distributed generalized Nash equilibria computation,” Automatica, vol. 102, pp. 111–121, 2019.
  * [16] P. Yi and L. Pavel, “Distributed generalized Nash equilibria computation of monotone games via double-layer preconditioned proximal-point algorithms,” IEEE Transactions on Control of Network Systems, vol. 6, no. 1, pp. 299–311, 2019.
  * [17] B. Franci, M. Staudigl, and S. Grammatico, “Distributed forward-backward (half) forward algorithms for generalized Nash equilibrium seeking,” in Proceedings of the 2020 European Control Conference (ECC), pp. 1274–1279, IEEE, 2020.
  * [18] G. Belgioioso, A. Nedić, and S. Grammatico, “Distributed generalized Nash equilibrium seeking in aggregative games on time-varying networks,” IEEE Transactions on Automatic Control, vol. 66, no. 5, pp. 2061–2075, 2020\.
  * [19] M. Bianchi, G. Belgioioso, and S. Grammatico, “Fast generalized Nash equilibrium seeking under partial-decision information,” Automatica, vol. 136, no. 110080, 2022.
  * [20] J. R. Marden and T. Roughgarden, “Generalized efficiency bounds in distributed resource allocation,” IEEE Transactions on Automatic Control, vol. 59, no. 3, pp. 571–584, 2014.
  * [21] A. Dreves, “How to select a solution in generalized Nash equilibrium problems,” Journal of Optimization Theory and Applications, vol. 178, no. 3, pp. 973–997, 2018.
  * [22] A. Dreves, “An algorithm for equilibrium selection in generalized Nash equilibrium problems,” Computational Optimization and Applications, vol. 73, no. 3, pp. 821–837, 2019.
  * [23] A. Simonetto, E. Dall’Anese, S. Paternain, G. Leus, and G. B. Giannakis, “Time-varying convex optimization: Time-structured algorithms and applications,” Proceedings of the IEEE, vol. 108, no. 11, pp. 2032–2048, 2020.
  * [24] K. Lu, G. Li, and L. Wang, “Online distributed algorithms for seeking generalized Nash equilibria in dynamic environments,” IEEE Transactions on Automatic Control, vol. 66, no. 5, pp. 2289–2296, 2021.
  * [25] M. Meng, X. Li, Y. Hong, J. Chen, and L. Wang, “Decentralized online learning for noncooperative games in dynamic environments,” arXiv preprint, https://arxiv.org/abs/2105.06200, 2021.
  * [26] F. Facchinei and J.-S. Pang, Finite-dimensional variational inequalities and complementarity problems.  Springer, 2007.
  * [27] H. H. Bauschke and P. L. Combettes, Convex analysis and monotone operator theory in Hilbert spaces.  Springer, 2011.
  * [28] I. Yamada and N. Ogura, “Hybrid steepest descent method for variational inequality problem over the fixed point set of certain quasi-nonexpansive mappings,” Numerical Functional Analysis and Optimization, vol. 25, no. 7-8, pp. 619–655, 2005.
  * [29] H. Xu and T. Kim, “Convergence of hybrid steepest-descent methods for variational inequalities,” Journal of Optimization Theory and Applications, vol. 119, no. 1, pp. 185–201, 2003.
  * [30] A. Cegielski and R. Zalas, “Properties of a class of approximately shrinking operators and their applications,” Fixed Point Theory, vol. 15, no. 2, pp. 399–426, 2014.
  * [31] N. Bastianello, A. Simonetto, and R. Carli, “Primal and dual prediction-correction methods for time-varying convex optimization,” arXiv preprint, arXiv:2004.11709, 2020.
  * [32] A. Simonetto, “Time varying convex optimization via time-varying averaged operators,” arXiv preprint, arXiv:1704.07338, 2017.
  * [33] A. Cegielski, A. Gibali, S. Reich, and R. Zalas, “An algorithm for solving the variational inequality problem over the fixed point set of a quasi-nonexpansive operator in Euclidean space,” Numerical Functional Analysis and Optimization, vol. 34, no. 10, pp. 1067–1096, 2013.
  * [34] D. P. Palomar and Y. C. Eldar, Convex optimization in signal processing and communications.  Cambridge university press, 2010.
  * [35] 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.
  * [36] Y. Malitsky and M. K. Tam, “A forward-backward splitting method for monotone inclusions without cocoercivity,” SIAM journal on optimization, vol. 30, no. 2, pp. 1451–1472, 2020.
  * [37] 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, 2000.
  * [38] A. L. Dontchev and R. T. Rockafellar, Implicit Functions and Solution Mappings.  Springer Science & Business Media, 2014.
  * [39] M. Zinkevich, “Online convex programming and generalized infinitesimal gradient ascent,” Proceedings of the 20th international conference on machine learning, 2003.
  * [40] E. Dall’Anese and A. Simonetto, “Optimal power flow pursuit,” IEEE Transactions on Smart Grid, 2016.
  * [41] N. Ogura and I. Yamada, “Nonstrictly convex minimization over the bounded fixed point set of a nonexpansive mapping,” Numerical Functional Analysis and Optimization, vol. 24, no. 1-2, pp. 129–135, 2003.