# A priori data-driven robustness guarantees on strategic deviations from
generalised Nash equilibria

George Pantazis<EMAIL_ADDRESS>Filiberto Fele<EMAIL_ADDRESS>Kostas
Margellos<EMAIL_ADDRESS>Delft Center for Systems and Control,
TU Delft, 2628 CD Delft, Netherlands  Department of Systems Engineering and
Automation, University of Sevilla, Spain  Department of Engineering Science,
University of Oxford, OX1 3PJ, UK

###### Abstract

In this paper we focus on noncooperative games with uncertain constraints
coupling the agents’ decisions. We consider a setting where bounded deviations
of agents’ decisions from the equilibrium are possible, and uncertain
constraints are inferred from data. Building upon recent advances in the so
called scenario approach, we propose a randomised algorithm that returns a
nominal equilibrium such that a _pre-specified_ bound on the probability of
violation for yet unseen constraints is satisfied for an entire region of
admissible deviations surrounding it—thus supporting neighbourhoods of
equilibria with probabilistic feasibility certificates. For the case in which
the game admits a potential function, whose minimum coincides with the social
welfare optimum of the population, the proposed algorithmic scheme opens the
road to achieve a trade-off between the guaranteed feasibility levels of the
region surrounding the nominal equilibrium, and its system-level efficiency.
Detailed numerical simulations corroborate our theoretical results.

###### keywords:

Generalised equilibrium problem; Uncertain games; Randomized methods;
Stochastic control and game theory; Multi-agent systems.

††thanks: This work was partially supported by the UK Engineering and Physical
Sciences Research Council under grant agreement EP/P03277X/1.

, ,

##  1 Introduction

The study of noncooperative games plays a significant role in a panoply of
applications ranging from smart-grids [41] to communication [43] and social
networks [2]. In these setups, agents can be modelled as self-interested
entities that interact with each other and make decisions based on possibly
misaligned individual criteria, while being subject to constraints (local or
global) that restrict the domain of their choices. Even though a variety of
systems can be analysed by means of deterministic game-theoretic tools [43,
23, 35], in many applications the decision making procedure is affected by
uncertainty. A number of results in the literature have explicitly addressed
uncertainty in a noncooperative setting, on the basis of specific assumptions
on the probability distribution [14, 44] and/or the geometry of the sample
space of the uncertainty [3, 25]; however, such assumptions can restrict the
applicability of these methods.

To circumvent these limitations, recent developments adopt a data-driven
perspective, focusing on the connection of game theory with the so called
scenario approach [8]. This is based on the idea that an optimisation problem
with constraints parametrised by an uncertain parameter—with fixed but
possibly unknown support set and probability distribution—can be approximated
by drawing samples of that parameter and solving the problem subject to the
constraints generated by those samples only; this approximation is known as
the _scenario program_. Standard results in the scenario approach [10, 12, 7]
provide certificates on the probability that a new yet unseen constraint will
violate the randomised solution obtained by solving the scenario program.

While the aforementioned results apply to uncertain _convex_ optimisation
problems, the works [11] and [13] paved the way towards the provision of data-
driven robustness guarantees to solutions of more general nonconvex problems.
In [21, 22, 34], these theoretical advancements were leveraged for the first
time in a game-theoretic context, for the formulation of distribution-free
probabilistic feasibility guarantees for randomised Nash equilibria. These
works provide guarantees for one specific equilibrium point (often assumed to
be unique); this was extended in [36, 37], by providing _a posteriori_
feasibility guarantees for the entire domain. Besides the game-theoretic
context, alternative methodologies for set-oriented probabilistic feasibility
guarantees have been proposed in the seminal works [15, 5], which a priori
characterise probabilistic feasibility regions constructed out of sampled
constraints using statistical learning theoretic results. More recently, the
so called probabilistic scaling [4, 31] has been proposed to obtain a
posteriori guarantees on the probability that a polytope generated out of
samples is a subset of some chance-constrained feasibility region. Following
an approach similar to [36], the works [16, 17] deliver tighter probabilistic
feasibility guarantees by focusing on variational-inequality (VI) solution
sets.

The results above follow a standard approach in the game-theoretic literature,
where a strict behavioural assumption—the so called _rationality_ —is imposed
on the players’ decision making. Namely, the players are viewed as rational
agents wishing to maximize their profit (expressed by some given cost
function). However, studies have shown that this is unrealistic in practice
[38, 9, 27, 45] and that agents usually exhibit a _boundedly rational_
behaviour [39], i.e., their decisions can deviate from rationality due to
individual biases, behavioural inertia, restricted computational power/time,
etc. The consequences of this become relevant in engineering applications, as
the human role in technical systems evolves beyond mere users and consumers to
active agents, operators, decision-makers and enablers of efficient, resilient
and sustainable infrastructures [30].

To bridge this gap between real-world applications and the cognate literature,
here we study games with uncertain constraints, where deviations from a
_nominal_ equilibrium are explicitly considered. We follow a randomised
approach to approximate the coupling constraints by means of data. In this
more general setting, where deviations are considered, providing guarantees
for a single solution is devoid of any meaning: indeed, repetition of the game
might lead to a different solution in a neighbourhood around the nominal
equilibrium, irrespective of the employed dataset. Technically speaking, this
renders the identification of the data samples that support the solution (cf.
sample compression [32]) a challenging task. Focussing on the class of
generalised Nash equilibrium seeking problems [18], we contribute to the
provision of data-driven robustness guarantees for the collection of possible
deviations from the equilibrium as follows:

  1. 1.

Adopting a scenario-theoretic paradigm, we establish a methodology for the
provision of _a posteriori_ probabilistic feasibility guarantees for a domain
around the randomised equilibrium of the game under study.

  2. 2.

We design a data-driven equilibrium-seeking algorithm that converges to a
solution that meets an _a priori_ defined level of probabilistic feasibility
for a fixed admissible region surrounding this equilibrium. This can model
possible deviations from a nominal equilibrium that the designer wishes to
take into account when incentivising a certain operation profile. The strength
of the provided feasibility guarantees depends on a prespecified quantity,
which in turn can affect the location of the nominal equilibrium and the
volume of the region for which these probabilistic guarantees hold.
Furthermore, when the game under study admits a potential function—whose
minimum coincides with some social welfare optimum—our methodology provides a
new perspective for trading off the probabilistic feasibility of the region
surrounding the nominal equilibrium and its system-level efficiency.

The rest of the paper is organized as follows. In Section 2 we provide
fundamentals of game theory and the scenario approach which will be used as
main ingredients for the subsequent developments. In Section 3.1 we show how
the feasibility guarantees for a region around the game solution can be a
posteriori quantified. In Section 3.2 we propose a data driven algorithm and
prove its convergence to an equilibrium such that the considered neighbourhood
of strategic deviations can satisfy prespecified probabilistic feasibility
requirements. An illustrative example in Section 4 corroborates our
theoretical analysis. Section 5 concludes the paper and presents future
research directions. To streamline the presentation of our results, some
proofs are deferred to the Appendix.

##  2 Preliminaries

_Notation_ : All vectors are column unless otherwise indicated. We denote by
$\mathbb{R}^{n}_{+}$ the nonnegative orthant in $\mathbb{R}^{n}$. When a
matrix $A$ is positive definite we write $A\succ 0$; similarly, positive semi-
definitess is denoted as $A\succeq 0$. Note that our definition of
(semi-)definiteness does not require the matrix to be symmetric. We denote by
$\mathbf{0}_{q\times r}$ a $q\times r$ matrix full of zeros, by $I_{r}$ the
$r\times r$ identity matrix, and by $\mathbf{1}_{r}$ the vector of $r$ ones;
dimensions can be omitted when clear from the context. $e_{q}$ is the unit
vector whose $q$-th element is 1 and all other elements are 0, $\|\cdot\|_{p}$
the $p$-norm operator, and $(\cdot)_{r}$ denotes the $r$-th component of its
vector argument.
$\mathbb{B}_{p}(x,\rho)=\\{y\in\mathbb{R}^{d}:\|y-x\|_{p}<\rho\\}$ is the open
$p$-normed ball centred at $x$ with radius $\rho$; when $p$ is omitted, any
choice of norm is valid. For a set $S$, $|S|$ denotes its cardinality, while
$2^{S}$ denotes its power set, i.e., the collection of all subsets of $S$.
Finally, given $D\succ 0$,
$\text{proj}_{K,D}[x]:=\operatorname*{arg\,min}_{y\in
K}(y-x)\raisebox{1.0pt}{${}^{\intercal}$}D(y-x)$ is the skewed projection of
$x$ onto the set $K$.

###  2.1 Games with uncertain constraints

We consider a population of agents with index set
$\mathcal{N}=\\{1,\dots,N\\}$. The decision vector $x_{i}$ of each agent
$i\in\mathcal{N}$ takes value in the set $X_{i}\subseteq\mathbb{R}^{n}$, while
$x=(x_{i})_{i=1}^{N}\in X=\prod_{i=1}^{N}X_{i}\subseteq\mathbb{R}^{nN}$ is the
global decision vector that is formed by concatenating the decisions of the
entire population. The vector $x_{-i}\in\mathbb{R}^{n(N-1)}$ comprises all
agents’ decisions except for those of agent $i$. In our setup, the cost
incurred by agent $i\in\mathcal{N}$ is expressed by a real-valued function
$J_{i}(x_{i},x_{-i})$ that depends on local decisions as well as on the
decisions from other agents $j\in\mathcal{N}\setminus\\{i\\}$. In the
following, with a slight abuse of notation, we can exchange $x$ for
$(x_{i},x_{-i})$ to single out agent $i$’s decision from the global decision
vector. Furthermore, we consider _uncertain_ constraints coupling the agents’
decisions. These can be expressed in the form111This formulation can describe
deterministic and/or local constraints as special cases.

$C_{\delta}=\\{x\in X:g(x,\delta)\leq 0\\},\;\delta\in\Delta,$ (1)

where $g:\mathbb{R}^{nN}\times\Delta\rightarrow\mathbb{R}$ depends on some
uncertain parameter $\delta$ taking values in a support set $\Delta$ according
to a probability measure $\mathbb{P}$.

Feasible collective decisions under this setup can be found by letting every
agent $i\in\mathcal{N}$ solve the following optimization program, where the
decisions $x_{-i}$ of all other agents are given,

$\left.\begin{aligned} G:\;&\min_{x_{i}\in X_{i}}J_{i}(x_{i},x_{-i})\\\
&\text{ subject to
}x_{i}\in\bigcap_{\delta\in\Delta}C^{i}_{\delta}(x_{-i})\end{aligned}\right\\}\;\forall
i\in\mathcal{N},$ (2)

where $C^{i}_{\delta}(x_{-i})=\\{x_{i}\in X_{i}:g(x_{i},x_{-i},\delta)\leq
0\\}$ denotes the projection of the coupling constraint on $X_{i}$ for fixed
$x_{-i}$ and uncertain realization $\delta\in\Delta$. The collection of
coupled optimization programs in (2) for all $i\in\mathcal{N}$ constitutes an
_uncertain noncooperative game_ ; we denote it as $G$.

Note that $G$ follows a worst-case paradigm, taking into account all possible
coupling constraints that can be realised by variations of the uncertain
parameter $\delta\in\Delta$. This can render the solutions of $G$ rather
conservative. Furthermore, it is in general not possible to compute a solution
for $G$ without an accurate knowledge of and/or additional assumptions on the
support set $\Delta$ and the probability distribution $\mathbb{P}$. To
circumvent these limitations, we follow a data-driven paradigm and approximate
$G$ by means of a finite number of samples drawn from $\Delta$, namely the
$K$-multisample
$\boldsymbol{\delta}_{K}=(\delta_{1},\dots,\delta_{K})\in\Delta^{K}$. In the
remainder of this document, we hold on to the standing assumption that these
samples are independent and identically distributed (i.i.d.). Apart from this,
no other knowledge on the support set $\Delta$ and the probability
distribution $\mathbb{P}$ of the uncertain parameter is required. Then, for
given multi-sample $\boldsymbol{\delta}_{K}$ and other agents’ decisions
$x_{-i}$, agent $i\in\mathcal{N}$ solves the optimization program

$\left.\begin{aligned} G_{K}:\;&\min_{x_{i}\in X_{i}}J_{i}(x_{i},x_{-i})\\\
&\text{ subject to }\
x_{i}\in\bigcap_{k=1}^{K}C^{i}_{\delta_{k}}(x_{-i})\end{aligned}\right\\}\;\forall
i\in\mathcal{N}.$ (3)

Instead of considering all possible uncertainty realizations $\delta\in\Delta$
as in (2), we let the data encoded in $\boldsymbol{\delta}_{K}$ lead agents to
their decision by solving (3). We refer to the collection of coupled
optimization programs in (3) as the _scenario game_ $G_{K}$, where the
subscript $K$ implies dependence on the drawn multi-sample
$\boldsymbol{\delta}_{K}$. Under standard assumptions, a solution to the
scenario game $G_{K}$ exists and the problem is, in contrast to $G$, tractable
using state-of-the-art equilibrium seeking algorithms.

###  2.2 Variational inequalities and game equilibria

Notably—under certain assumptions detailed next—solutions to the game $G_{K}$
can be retrieved as solutions to a variational inequality (VI), for specific
choices of the mapping $F:\,X\to\mathbb{R}^{nN}$ [18, Thm 3.9]:

$\begin{split}\text{VI}_{K}:&\ \text{Find}\ x^{*}\in\Pi_{K}\ \text{such
that}\\\ &\ (x-x^{*})\raisebox{1.0pt}{${}^{\intercal}$}F(x^{*})\geq 0\
\text{for any}\ x\in\Pi_{K},\end{split}$ (4)

where $\Pi_{K}:=X\cap\bigcap_{k=1}^{K}C_{\delta_{k}}$ denotes the problem
domain. A classic game solution concept, which encounters wide application in
the literature, is the Nash equilibrium (NE) [33]. At a NE, no agent can
decrease their cost by unilaterally changing their decision. Formally, this
can be stated as follows.

######  Definition 1.

A point $x^{*}=(x_{i}^{*},x_{-i}^{*})\in\Pi_{K}$ is called a generalised Nash
equilibrium (GNE) of $G_{K}$ if, for all $i\in\mathcal{N}$,

$J_{i}\left(x_{i}^{*},x_{-i}^{*}\right)\leq J_{i}(y_{i},x_{-i}^{*}),\;\forall
y_{i}\in X_{i}\cap\bigcap_{k=1}^{K}C^{i}_{\delta_{k}}(x_{-i}^{*}).$

For our analysis, we rely on the following conditions:

######  Assumption 1.

For all $i\in\mathcal{N}$, $J_{i}(x_{i},x_{-i})$ is convex and continuously
differentiable in $x_{i}$ for any fixed $x_{-i}$.

######  Assumption 2.

  1. 1.

For any multi-sample $\boldsymbol{\delta}_{K}\in\Delta^{K}$, the domain
$\Pi_{K}$ is non-empty.

  2. 2.

The set $X=\prod_{i=1}^{N}X_{i}$ is compact, polytopic and convex.

  3. 3.

For any $\delta\in\Delta$, $g$ is an affine function of the form
$g(x,\delta)=a(\delta)\raisebox{1.0pt}{${}^{\intercal}$}x-b(\delta)$, where
$a:\Delta\rightarrow\mathbb{R}^{nN}$ and $b:\Delta\rightarrow\mathbb{R}$.

Under these assumptions, we can determine a GNE as in Definition 1 by solving
(4) with

$F(x)=F_{\mathrm{NE}}(x):=\begin{bmatrix}\nabla_{x_{1}}J_{1}(x_{1},x_{-1})\\\
\vdots\\\ \nabla_{x_{N}}J_{N}(x_{N},x_{-N})\end{bmatrix}.$ (5)

A class of problems of common interest can be modelled by the so called
_aggregative_ games [28, 26, 1], where the cost incurred by agents depends on
some aggregate measure—typically the average—of the decision of the entire
population. Formally, such a cost can be expressed in (3) by the real-valued
function $J_{i}(x_{i},\sigma(x))$, where the aggregate
$\sigma:\mathbb{R}^{nN}\rightarrow\mathbb{R}^{n}$ is defined as the mapping
$x\mapsto\frac{1}{N}\sum_{i=1}^{N}x_{i}$. A solution frequently linked to this
class of games is the Wardrop equilibrium (WE), a concept akin to the NE but
specifically defined in the context of transportation networks [6]. The
variational WEs of $G_{K}$ can be expressed by using
$F(x)=F_{\mathrm{WE}}(x):=[\nabla_{x_{i}}J_{i}(x_{i},z)_{|z=\sigma(x)}]_{i\in\mathcal{N}}$;
notice that in this case the second argument of $J_{i}$ is fixed and set to
$\sigma(x)$, consistently with the notion of WE where agents neglect the
impact of their decision on others.

We restrict the considered class of variational mappings as follows:

######  Assumption 3.

The mapping $F$ is

  1. 1.

$\alpha$-strongly monotone, i.e.,
$(x-y)\raisebox{1.0pt}{${}^{\intercal}$}(F(x)-F(y))\geq\alpha\|x-y\|^{2}$ for
any $x,y\in X$,

  2. 2.

$L_{F}$-Lipschitz continuous., i.e., $\|F(x)-F(y)\|\leq L_{F}\|x-y\|$ for any
$x,y\in X$.

Assumptions 1 and 3 are standard in the game-theoretic literature [19, 43].
Assumption 2 is relatively mild; the affine form of the constraints is
exploited in the proposed algorithm (see Section 3) for the convergence to an
equilibrium bearing the desired robustness properties.

We point out that in general only a subset of solutions to $G_{K}$ can be
retrieved through (4): these are referred to as _variational equilibria_ , and
enjoy favourable properties over nonvariational ones, as with the former the
coupling constraints’ burden is equally split among agents [24, 29]. The
following lemma, adapted from [19, Thm. 2.3.3], formalises the connection
between the solutions to VIK and the GNEs (or GWEs) of $G_{K}$.

######  Lemma 1.

Under Assumptions 1, 2 and 3, VIK has a unique solution that is also an
equilibrium of $G_{K}$.

For the considered class of VIs, several algorithms from the literature can be
employed to obtain a variational equilibrium of $G_{K}$; see, e.g., [18, 35].
We remark that, even if not explicitly shown for ease of notation, any
solution $x^{*}$ to $G_{K}$ is itself a function of the drawn multisample
$\boldsymbol{\delta}_{K}\in\Delta^{K}$. Probabilistic feasibility guarantees
for the unique solution of VIK can then be provided both in an _a priori_ and
_a posteriori_ fashion by resorting to the results in [21, 22, 34]. However,
these results are tailored to the provision of probabilistic feasibility
guarantees for a single point (namely the solution of a VI): any strategic
deviation from the equilibrium is not covered by such guarantees. We cover
this issue in Section 3. First, we provide some background on the scenario
approach.

###  2.3 Basic concepts in the scenario approach

A fundamental notion in the scenario approach is the _probability of
violation_ of an uncertain constraint.

######  Definition 2.

  1. 1.

The probability of violation $\mathcal{V}:\mathbb{R}^{nN}\rightarrow[0,1]$ of
a point $x\in\Pi_{K}$ is defined as the probability that a new yet unseen
sample $\delta\in\Delta$ will give rise to a constraint $C_{\delta}$ (as
defined in (1)) such that $x\notin C_{\delta}$, i.e.,

$\displaystyle\mathcal{V}(x):=\mathbb{P}\\{\delta\in\Delta:x\notin
C_{\delta}\\}.$

  2. 2.

The probability of violation $\mathbb{V}:2^{\mathbb{R}^{nN}}\rightarrow[0,1]$
of a set $S\subseteq\Pi_{K}$ is defined as the worst-case $\mathcal{V}$ among
all the points in $S$, i.e.,

$\displaystyle\mathbb{V}(S)=\sup\limits_{x\in
S}\mathbb{P}\\{\delta\in\Delta:x\notin C_{\delta}\\}.$

Formally, a data-driven decision-making process can be characterized by a
mapping—the _algorithm_ —that takes as input the data encoded by the samples
and returns a set of decisions.

######  Definition 3.

An algorithm is a function $\mathscr{A}:\Delta^{l}\rightarrow
2^{\mathbb{R}^{nN}}\times\mathbb{R}^{nN}$ that takes as input an
$l$-multisample and returns the pair $(S_{l}^{*},x^{*})$, namely, a solution
set $S_{l}^{*}$ and a point $x^{*}\in S_{l}$.

In the following, we interpret the above definition as context-dependent, in
that the size $l$ of the input multisample is admitted to vary—all else
remaining fixed for a given algorithm $\mathscr{A}$. A key notion, strongly
linked to that of algorithm, is the _minimal compression set_ [32]. This
concept springs from the observation that typically only a subset of the
sampled data is relevant to a decision or set of decisions, and all other
samples are redundant.

######  Definition 4 (Compression set).

Consider an algorithm $\mathscr{A}$ as in Definition 3. A subset of samples
$I\subseteq\boldsymbol{\delta}_{K}$ is called a _compression_ for
$\mathscr{A}(\boldsymbol{\delta}_{K})$ if
$\mathscr{A}(I)=\mathscr{A}(\boldsymbol{\delta}_{K})$.222With some abuse of
notation, in the remainder the symbol $\boldsymbol{\delta}_{K}$ is interpreted
as either the i.i.d. sample vector $\boldsymbol{\delta}_{K}\in\Delta^{K}$, or
the set comprising its components, i.e.,
$\boldsymbol{\delta}_{K}=\\{\delta_{1},\ldots,\delta_{K}\\}\subseteq\Delta$,
depending on the context. As multiple subset of samples can exist that fulfil
this property, the ones with the minimal cardinality are called _minimal_
compression sets.

If we feed the algorithm with the set of samples corresponding to a
compression, then the same decision (in our case a set and a point in the set)
will be returned as when we feed the algorithm with the entire multi-sample.
The compression set is related to the notion of support samples that is
presented below.

######  Definition 5 (Support sample).

A sample $\delta_{i}\in\boldsymbol{\delta}_{K}$ is a _support sample_ if its
removal changes the decision returned by $\mathscr{A}$, i.e., if
$\mathscr{A}(\boldsymbol{\delta}_{K})\neq\mathscr{A}(\boldsymbol{\delta}_{K}\setminus\delta_{i})$.

Notice that since the output of $\mathscr{A}$ is a pair, it follows directly
from the previous definition that if a support sample is removed, at least one
of the two elements of the pair will have to change for the algorithm’s output
to change as well. A sample whose removal changes the $x^{*}$ component of
$\mathscr{A}$ is called a support sample for $x^{*}$. We link the two previous
concepts by imposing the following _non-degeneracy_ assumption.

######  Assumption 4.

With $\mathbb{P}^{K}$-probability 1, if $I^{*}$ is the set of support samples
of $(S_{K}^{*},x^{*})=\mathscr{A}(\boldsymbol{\delta}_{K})$, then
$\mathscr{A}(\boldsymbol{\delta}_{K})=\mathscr{A}(I^{*})$.

In other words, non-degeneracy implies that the solution returned by algorithm
$\mathscr{A}$ when fed with all the drawn samples is the same as the one
obtained by using only the support samples of the problem. Under such
assumption, the set of support samples and the minimal compression set
coincide—and the latter is unique. This avoids degenerate cases where there is
an accumulation of constraints active at the solution.

##  3 Probabilistic feasibility of sets around equilibria

###  3.1 A first a posteriori result

Returning to the scenario game $G_{K}$ in (3), we now consider a more general
setup where agents are allowed to deviate from $x^{*}$ following, e.g.,
unmodelled changes in their cost functions; while we suppose that these
deviations are feasible with respect to the local constraints, we want to
study the feasibility as regards the coupling constraints obtained through
sampling. Specifically, the region in which agents’ strategies can deviate
from the nominal equilibrium is assumed to lie within a predefined open ball
$\mathbb{B}(x^{*},\rho)$, where $\rho>0$ is a fixed radius that denotes the
maximum possible distance of agents’ deviations from $x^{*}$; the latter is
assumed to be unique as per Lemma 1. As such, the region of interest is
$S_{K}^{*}=\Pi_{K}\cap\mathbb{B}(x^{*},\rho)$.

This is pictorially illustrated in Figure 1 using the ${\infty}$-norm (note
that any other norm could have been used instead): an algorithm $\mathscr{A}$
(see Sec. 2.3) takes as input a multi-sample $\boldsymbol{\delta}_{K}$ and
returns the region $S_{K}^{*}$ around the solution $x^{*}\in\mathbb{R}^{2}$ of
a game with two players whose decisions are defined as scalar quantities. We
assume in this case that $\Pi_{K}$ is shaped exclusively by sampled coupling
constraints: only two constraints (in solid blue) are of support for the set,
as the removal of any of the constraints depicted by dashed red lines will not
change $S_{K}^{*}$. However, depending on the given multi-sample
$\boldsymbol{\delta}_{K}$ and the location of the resulting equilibrium, more
constraints could intersect the surrounding region $\mathbb{B}(x^{*},\rho)$,
thus increasing the cardinality of the support sample set.

\begin{overpic}[trim=0.0pt 40.0pt 0.0pt
85.35826pt,clip,width=184.28981pt]{final_affine_M_2}
\put(40.0,30.0){$\Pi_{K}$} \put(34.0,21.0){$S_{K}^{*}$}
\put(27.0,19.0){$\rho$} \put(28.5,10.7){$x^{*}$} \end{overpic} Figure 1:
Illustration of region $S_{K}$ (in dark green) as the intersection of the set
of deviations $\mathbb{B}_{\infty}(x^{*},\rho)$ around the unique equilibrium
$x^{*}$ (red dot) with the domain $\Pi_{K}$.

We can quantify the number of samples that support $S_{K}^{*}$ in an _a
posteriori_ fashion as established in Theorem 1. To this end, for a fixed
confidence $\beta\in(0,1)$, let the violation level be defined as a function
$\epsilon:\\{0,\dots,K\\}\rightarrow[0,1]$ satisfying

$\epsilon(K)=1\;\text{and}\;\sum_{i=0}^{K-1}{\binom{K}{i}}(1-\epsilon)^{K-i}=\beta,$
(6)

where $K$ is the size of the multisample; see, e.g., [13, Eq. (7)].

######  Theorem 1.

Under Assumptions 1–4, consider some algorithm $\mathscr{A}$ returning a pair
$(S^{*}_{K},x^{*})$ where $S^{*}_{K}$ is parametrised by $x^{*}$. Fix a
confidence parameter $\beta\in(0,1)$ and a violation level
$\epsilon:\\{0,\dots,K\\}\rightarrow[0,1]$ that satisfies (6). We have that

$\mathbb{P}^{K}\left\\{\boldsymbol{\delta}_{K}\in\Delta^{K}:~{}\mathbb{V}(S_{K}^{*})>\epsilon(s^{*}+M)\right\\}\leq\beta,$

where $s^{*}$ is the number of samples that support the equilibrium $x^{*}$
and $M$ the number of facets of $\Pi_{K}$ that intersect $S_{K}^{*}$.

_Proof_ : Let $(S_{K}^{*},x^{*})$ be the solution returned by $\mathscr{A}$
for some given $\boldsymbol{\delta}_{K}$, according to Definition 3. We aim at
determining a compression set for $\mathscr{A}(\boldsymbol{\delta}_{K})$, as
this would in turn allow us to obtain the theorem’s conclusion by means of
Theorem 2 in [37]. Under Assumption 4, this is equivalent to determining an
upper bound on the number of support samples for
$\mathscr{A}(\boldsymbol{\delta}_{K})$ (see Definition 5). Since by
construction $S_{K}^{*}$ is parametrised by $x^{*}$, the support samples of
$\mathscr{A}(\boldsymbol{\delta}_{K})$ are given by the union of (i) the
samples that support the equilibrium $x^{*}$, and (ii) the samples that do not
support $x^{*}$, but whose removal can still lead to a change of the region
$S_{K}^{*}$. In the first case, removing a support sample for $x^{*}$ would
move the entire region (e.g., by its centre, as in Fig. 1); we denote the
number of such samples as $s^{*}$. On the other hand, the removal of a sample
corresponding to the second case yields a larger region $S^{*}_{K}$, leaving
$x^{*}$ unaffected; the number of such samples can be upper bounded by the $M$
facets of $\Pi_{K}$ that intersect $S_{K}^{*}$. Hence, the number of samples
that form a compression set for $\mathscr{A}(\boldsymbol{\delta}_{K})$ is
bounded by $s^{*}+M$. Existence of a compression set $I$ with a bound on its
cardinality is sufficient for the application of Theorem 2 in [37]. The fact
that for the minimal compression set $|I^{*}|\leq|I|\leq s^{*}+M$ always holds
leads then to the statement of this theorem. $\blacksquare$

It is important to stress that the application of Theorem 1 is agnostic on the
choice of the equilibrium seeking algorithm. To use the result of Theorem 1,
one needs to quantify (an upper bound of) the number of samples $s^{*}$ that
support the randomised equilibrium $x^{*}$ and (an upper bound of) the number
$M$ of coupling constraints that correspond to facets of $S_{K}^{*}$. While
$s^{*}\leq nN$ under Assumptions 1–3,333By arguments similar to those in [37]
it can be shown that a tighter bound $s^{*}\leq n$ holds for the game $G_{K}$
in case coupling constraints only concern the aggregate variable; see also
[42]. an upper bound for $M$ can in general only be achieved _a posteriori_ ,
i.e., once $\boldsymbol{\delta}_{K}$ is sampled. Then an important question
naturally arises: _When we allow for bounded deviations about some
equilibrium, is it possible to devise an algorithm that converges to a
solution such that an a priori defined robustness for all points in the region
of admissible deviations is achieved?_ We aim at answering this question next.

###  3.2 A priori probabilistic certificates

Consider the scenario game $G_{K}$ and suppose that bounded deviations from
the solution are allowed. We model such deviations as a ball of radius $\rho$
around the equilibrium, as in Section 3.1. In contrast to the a posteriori
nature of the result therein, our goal here is to achieve an a priori bound.
Namely, we aim at establishing the main statement of Theorem 1 with a
_prespecified_ violation level, which does not depend on the given multi-
sample $\boldsymbol{\delta}_{K}$. In other words, we seek a statement—holding
with known confidence—of the form
$\mathbb{V}(\Pi_{K}\cap\mathbb{B}(x^{*},\rho))<\bar{\epsilon}$, with
$\bar{\epsilon}\in(0,1)$ a priori fixed.

To achieve this, we build upon the previous conclusions, which expose a link
between the probability of constraint violation and the number $M$ of
uncertainty samples constituting facets of $\Pi_{K}$ that
$\mathbb{B}(x^{*},\rho)$ intersects. In particular, a monotonic relationship
follows from (6): the smaller the number $M$ of facets of the uncertain domain
intersected by $\mathbb{B}(x^{*},\rho)$, the better, i.e., less conservative,
the theoretical feasibility guarantees on constraint violation for the
strategies belonging to the feasible region $S_{K}^{*}$ surrounding the
equilibrium. Moreover, in certain cases (as illustrated next) a smaller value
of $M$ is associated with a larger region for which the guarantees of Theorem
1 hold—due to a smaller portion of $\mathbb{B}(x^{*},\rho)$ being cut off by
the intersection with $\Pi_{K}$. This motivates us to study the role of $M$ as
a modulating parameter for the robustness of the feasibility certificates
offered for the region $S_{K}^{*}$, as well as the extent of deviation from
the nominal equilibrium covered by such certificates. These concepts are
leveraged in the algorithm proposed next.

####  3.2.1 GNE-seeking algorithm with a priori robustness guarantees

We consider an iterative scheme to determine a solution of VIK in (4). In
particular, since we address a NE problem characterised by coupling
constraints, we build our Algorithm 1 upon a primal-dual scheme, where
constraint satisfaction is achieved by the use of Lagrange multipliers. As
deterministic constraints do not play a role in the evaluation of the
robustness guarantees, suppose for ease of exposition that $\Pi_{K}$ only
comprises uncertain coupling constraints. Let $A\in\mathbb{R}^{m\times nN}$
and $b\in\mathbb{R}^{m}$ such that

$\displaystyle\Pi_{K}=\\{x\in X:Ax\leq b\\},$ (7a)
$\displaystyle\|a_{\ell}\|_{2}=1,\;\ell=1,\ldots,m,$ (7b)

where $a_{\ell}\raisebox{1.0pt}{${}^{\intercal}$}$ denotes the $\ell$-th row
of $A$. Eq. (7) is the irredundant $H$-representation of the polytopic
feasibility region $\Pi_{K}$ defined in (4), where the rows of matrix $A$ are
unit vectors. Property (7b) is key to the second statement in Lemma 2. It
entails no loss of generality, since for any $A,b$ forming an equivalent
$H$-representation of $\Pi_{K}$, (7) can be obtained by normalising each row
of $A$ and the corresponding component of $b$ by the row-vector norm. Thus,
the pair $(A,b)$ encodes the set of randomised coupling constraints that
constitute facets of $\Pi_{K}$. Formally,
$A:\Delta^{K}\rightarrow\mathbb{R}^{m\times nN}$ and
$b:\Delta^{K}\rightarrow\mathbb{R}^{m}$ are mappings from the $K$-multisample
to the space of real $m\times nN$ matrices and $m$-dimensional vectors,
respectively. Although we do not reflect this in the notation for simplicity,
we will consider this to be implicit in the remainder of the paper.

Now, to understand the mechanism underlying Algorithm 1, first note that
tightening an affine constraint by a distance corresponding to the radius
$\rho$ is equivalent to preventing the original constraint from intersecting
$\mathbb{B}(x^{*},\rho)$. Based on this and the idea illustrated in previous
paragraphs, we leverage the information carried by Lagrange multipliers to
control the intersection between $\Pi_{K}$ and $\mathbb{B}(x^{(\kappa)},\rho)$
while the algorithm is running. In particular, all but the coupling
constraints that correspond to the $M$ largest multipliers are tightened in
Algorithm 1 by a distance $\rho$; as a result, the number of facets of
$\Pi_{K}$ intersecting $\mathbb{B}(x^{*},\rho)$ is at most $M$. This enables
us to obtain an _a priori_ estimate on the number of support samples, that in
turn allows us to provide a priori bounds on
$\mathbb{V}(\Pi_{K}\cap\mathbb{B}(x^{*},\rho))$ (see Theorem 3).

It is worth pointing out that the constraints corresponding to the $M$ largest
multipliers contribute the most in reducing the nominal cost if $x^{*}$ is
allowed to lie at their boundaries. It can therefore be understood that while
smaller values for $M$ can result in a more robust—and possibly larger—region
$S_{K}^{*}$, they can also move the location of the nominal equilibrium
$x^{*}$ to a less efficient point towards the interior of $\Pi_{K}$. As we
will demonstrate numerically in the sequel, this is indeed the case with
_potential_ games [20].

Algorithm 1 A priori robust GNE seeking algorithm

1:$y^{(0)},\rho\in\mathbb{R}_{+},\boldsymbol{\delta}_{K}\in\Delta^{K},M\leq
m,\xi\geq 0$

2:$\kappa\leftarrow 0$

3:repeat

4:
$y^{(\kappa+1)}\leftarrow\text{proj}_{X\times\mathcal{M},D}\left[y^{(\kappa)}-D^{-1}T(y^{(\kappa)},\rho,M)\right]$

5: $\kappa\leftarrow\kappa+1$

6:until $\|y^{(\kappa+1)}-y^{(\kappa)}\|\leq\xi$

7:$y^{*}\leftarrow y^{(\kappa+1)}$

8:return $y^{*}$ and $\Pi_{K}\cap\mathbb{B}(x^{*},\rho)$

We build upon the class of asymmetric projection algorithms (APA) [19, Ch. 12]
to seek a solution of the game $G_{K}$. This involves an iterative scheme,
whose trajectory is characterised by the dynamics in line 4 of Algorithm 1,
where
$y:=(x\raisebox{1.0pt}{${}^{\intercal}$},\mu\raisebox{1.0pt}{${}^{\intercal}$})\raisebox{1.0pt}{${}^{\intercal}$}\in\mathbb{R}^{nN+m}$
is the concatenation of the global decision vector $x$ and the Lagrange
multipliers
$\mu=(\mu_{\ell})_{\ell=1}^{m}\in\mathcal{M}\subseteq\mathbb{R}_{+}^{m}$,
$D\succ 0$ is the so called _asymmetric projection_ matrix, and the mapping
$T:\mathbb{R}^{nN}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{N}\rightarrow\mathbb{R}^{nN+m}$
is given by

$T(y,\rho,M)=\begin{bmatrix}F(x)+A\raisebox{1.0pt}{${}^{\intercal}$}\mu\\\
-\left(Ax-b+Q(\mu,M)\boldsymbol{\rho}\right)\end{bmatrix},$ (8)

where $\boldsymbol{\rho}=c\rho\mathbf{1}_{m}$, and $c,\rho,M$ are fixed during
the execution of the algorithm. Note that $T$ follows from the primal-dual
conditions of the game solution; see [18, Sec. 4.2], [19, Sec. 1.4.1]. In (8),
$F$ is defined as in Section 2.2, and $A,b$ are as in (7). Moreover, Algorithm
1 relies on the mapping
$Q:\mathbb{R}^{m}_{+}\times\mathbb{N}\rightarrow\\{0,1\\}^{m\times m}$ which
allows convergence to possibly different nominal solutions $x^{*}$, according
to the specified parameters $M$ and $\rho$. $Q$—formally introduced in the
next subsection—allows to perform constraint tightening along the trajectory
of the algorithm, following the concepts illustrated above.

As illustrated in Figure 2, $\mathcal{M}$ is the union of a finite number of
disjoint convex sets, hence the projection in line 3 of Algorithm 1 is
computationally viable (e.g., projecting on each convex subset of the union
and then setting $y^{(\kappa+1)}$ to be the solution among these projections
that results in the minimum distance from
$y^{(\kappa)}-D^{-1}T(y^{(\kappa)},\rho,M)$). Finally, it should be noted that
Algorithm 1 could be performed in a distributed manner, although this is
outside the scope of the present work.

####  3.2.2 Constraint tightening via mapping $Q$

We define the mapping $Q$ as

$Q(\mu,M)=P\raisebox{1.0pt}{${}^{\intercal}$}(\mu)R(M),$ (9)

where

  * •

$P:\mathbb{R}^{m}\rightarrow\\{0,1\\}^{m\times m}$ returns a permutation
matrix such that $P(\mu)\mu$ is the vector composed by the elements of $\mu$
arranged in decreasing order.

  * •

$R:\mathbb{N}\rightarrow\\{0,1\\}^{m\times m}$ takes as input the number of
coupling constraints $M\leq m$ we allow $\mathbb{B}(x^{*},\rho)$ to intersect
with and returns as output the matrix

$R(M)=\begin{bmatrix}\begin{matrix}\mbox{\large 0}_{m\times
M}\end{matrix}&\vline&\begin{matrix}0_{M\times m-M}\\\
I_{m-M}\end{matrix}\end{bmatrix}.$ (10)

Compatibly with the definition of $P(\cdot)$,
$R(M)P(\mu)\boldsymbol{\rho}=(\mathbf{0}_{M}\raisebox{1.0pt}{${}^{\intercal}$}\;c\rho\mathbf{1}_{m-M}\raisebox{1.0pt}{${}^{\intercal}$})\raisebox{1.0pt}{${}^{\intercal}$}=R(M)\boldsymbol{\rho}$,
where the last equality holds since all components of $\boldsymbol{\rho}$ are
equal.

It can be seen from (8) and the above definition that $Q(\cdot,M)$ allows to
tighten the constraints corresponding to the smallest $m-M$ multipliers. In
particular, following the discussion in Section 3.2.1, we consider tightening
these constraints by an amount equal to the radius of the sphere that
circumscribes $\mathbb{B}({x^{*}},\rho)$, for any choice of norm. This amount
is equal to $\boldsymbol{\rho}_{\ell}=c\rho\|a_{\ell}\|_{2}=c\rho$, where the
last equality is due to (7b); $c=1$ if $\mathbb{B}({\cdot},\rho)$ is expressed
by a $p$-norm with $p\leq 2$, $c=\sqrt{n}$ otherwise. Conversely, at most $M$
constraints can intersect $\mathbb{B}(x^{*},\rho)$ upon convergence of the
algorithm. Let $\mathcal{L}(M)\subseteq\\{1,\ldots,m\\}$ contain the indices
of the $M$ largest multipliers. Then,
$\ell\in\mathcal{L}(M)\Leftrightarrow(Q(\mu,M)\boldsymbol{\rho})_{\ell}=0$,
and the second block row of $T$ in (8) expresses

$\begin{cases}a_{\ell}\raisebox{1.0pt}{${}^{\intercal}$}x\leq
b_{\ell}&\text{if $(Q(\mu,M)\boldsymbol{\rho})_{\ell}=0$},\\\
a_{\ell}\raisebox{1.0pt}{${}^{\intercal}$}x\leq b_{\ell}-c\rho&\text{if
$(Q(\mu,M)\boldsymbol{\rho})_{\ell}=c\rho$}.\end{cases}$ (11)

The following example further clarifies this concept.

##### Illustrative example:

Suppose $\Pi_{K}$ is composed of 3 uncertain coupling constraints and that we
allow the region of strategic deviations $\mathbb{B}({\cdot},\rho)$ to
intersect at most $M=1$ of them. It follows from (10) that
$R(M)=\left[\begin{smallmatrix}0&0&0\\\ 0&1&0\\\
0&0&1\end{smallmatrix}\right]$. Now, suppose that at some iteration $\kappa$
of Algorithm 1 the multiplier vector is
$\mu^{(\kappa)}=(\mu^{(\kappa)}_{\ell})_{\ell=1}^{3}$ with

$\mu^{(\kappa)}_{2}>\mu^{(\kappa)}_{1}>\mu^{(\kappa)}_{3}.$ (12)

Then, $P(\mu^{(\kappa)})=\left[\begin{smallmatrix}0&1&0\\\ 1&0&0\\\
0&0&1\end{smallmatrix}\right]$ is a permutation matrix such that
$P(\mu^{(\kappa)})\mu^{(\kappa)}=(\mu^{(\kappa)}_{2}\;\mu^{(\kappa)}_{1}\;\mu^{(\kappa)}_{3})\raisebox{1.0pt}{${}^{\intercal}$}$.
Finally, we have
$Q(\mu^{(\kappa)},M)\boldsymbol{\rho}=P\raisebox{1.0pt}{${}^{\intercal}$}(\mu^{(\kappa)})R(M)\boldsymbol{\rho}=(c\rho\;0\;c\rho)\raisebox{1.0pt}{${}^{\intercal}$}$,
where $P\raisebox{1.0pt}{${}^{\intercal}$}(\cdot)$ applies the correct
ordering to the vector $R(M)\boldsymbol{\rho}$. Therefore, in case (12)
continues to hold for all $j\geq\kappa$, the region
$\mathbb{B}({x^{(j)}},\rho)$ will intersect the constraint associated to the
largest multiplier from some iteration $\hat{j}\geq\kappa$ until convergence.

###  3.3 Convergence analysis and main result

In this section we show convergence for Algorithm 1. Before proceeding, we
provide details on the set $\mathcal{M}\subseteq\mathbb{R}_{+}^{m}$ where the
Lagrange multipliers updates are projected (line 4, Algorithm 1). Let
$\mathcal{Z}:=[\zeta,+\infty)\cup\\{0\\}$, for some small $\zeta>0$, i.e.,
$\mathcal{Z}\subset\mathbb{R}$ contains all nonnegative scalars which take
value greater than $\zeta$ when nonzero.

######  Assumption 5.

$\mathcal{M}$ is compact and admits the form

$\mathcal{M}:=\left\\{\mu\in\mathbb{R}^{m}:(P(\mu)\mu)_{\ell+1}<(P(\mu)\mu)_{\ell}-\zeta,\right.\\\
\left.\forall\ell=1,\dots,m-1\right\\}\cap\mathcal{Z}^{m}.$ (13)

Recalling that $P(\mu)\mu$ rearranges the multipliers in descending order, the
set $\mathcal{M}$ contains all vectors where the difference between every pair
of strictly positive components—and the distance of the smallest of these from
zero—exceeds $\zeta$. We note that (13) can be expressed as the disjoint union
of $q=m!+m+1$ convex subsets of $\mathbb{R}^{m}_{+}$, each of which we denote
in the following as $\mathcal{M}_{j}$, i.e.,
$\mathcal{M}=\bigcup_{j=1}^{q}\mathcal{M}_{j}$. Figure 2 provides an
illustration of this set for the case $m=3$.

Figure 2: Illustration of the domain $\mathcal{M}$ of the Lagrange multipliers
associated to the coupling constraints, for the case $\zeta=0.2$ and $m=3$.
This results in $q=10$ convex subsets, including the origin and a portion of
the axes.

Restricting the Lagrange multipliers to $\mathcal{M}$ facilitates convergence
of the algorithm in cases where the VI$(X\times\mathcal{M},T)$ has multiple
solutions. In particular, it (i) alleviates discontinuity issues of the
mapping $T$, and (ii) endows the latter with desired nonexpansiveness
properties (see proof of Lemma 3.2). In the numerical implementation of the
algorithm, ensuring $\mu\in\mathcal{M}$ can possibly introduce small
perturbations in the multipliers—compared to standard formulations where
$\mu\in\mathbb{R}^{m}_{+}$—which in turn could produce a slight violation of
the constraints (this can be controlled through the magnitude of $\zeta$). We
further note that assuming compactness of $\mathcal{M}$ is only required for
the proof of Theorem 2 below; the assumption can be numerically satisfied by
intersecting (13) with an arbitrarily large compact set. While necessary for a
formal statement of convergence, we note that in many practical instances,
including our numerical study, the desired solution can be attained without
the need of imposing such an assumption on $\mathcal{M}$.

We are now ready to introduce the following lemmas as an intermediate step
towards proving convergence of Algorithm 1.

######  Lemma 3.1.

Define $T$ as in (8)–(9), where $A,b$ satisfy (7). Then, for any
$\mu,\mu^{\prime}\in\mathcal{M}$, $\mu\neq\mu^{\prime}$, there exists an
integer $0\leq h\leq M$ such that

$(\mu-\mu^{\prime})\raisebox{1.0pt}{${}^{\intercal}$}(Q(\mu,M)-Q(\mu^{\prime},M))\boldsymbol{\rho}\leq-h\zeta
c\rho.$ (14)

The construction of $h$ can be found in the proof of Lemma 3.1; nonetheless,
its exact value does not play a role in the subsequent analysis.

######  Lemma 3.2.

Consider $T$ as in (8)–(9), where $A,b$ satisfy (7) and
$\mathcal{M}=\bigcup_{j=1}^{q}\mathcal{M}_{j}$ as in (13). For each
$j=1,\ldots,q$, let VI$(X\times\mathcal{M}_{j},T)$ denote the VI problem
defined by the map $T$ restricted to the subdomain $X\times\mathcal{M}_{j}$.
Under Assumptions 1–3 the following holds:

  1. 1.

$T$ is continuous on $X\times\mathcal{M}$.

  2. 2.

Let

$D=\begin{bmatrix}\frac{1}{\tau}I_{nN}&0\\\
-2A&\frac{1}{\tau}I_{m}\end{bmatrix},$ (15)

and set $\tau>0$ such that

$\tau<\min\bigg{\\{}\frac{-L_{F}^{2}\\!+\\!\sqrt{L_{F}^{4}+4\alpha^{2}\|A\|^{2}}}{2\alpha\|A\|^{2}},\\\
\frac{-\rho(1\\!+\\!\|A\|^{2})\\!+\\!\sqrt{\rho^{2}(1\\!+\\!\|A\|^{2})^{2}+16\zeta^{2}\|A\|^{2}}}{4\zeta\|A\|^{2}}\bigg{\\}}.$
(16)

Then, for any $j=1,\ldots,q$, Algorithm 1 converges to a solution of
VI$(X\times\mathcal{M}_{j},T)$, when the gradient step in line 4 is projected
on the corresponding subdomain, for any $y^{(0)}\in X\times\mathcal{M}_{j}$.

The second part of Lemma 3.2 provides an admissible range of values for $\tau$
such that Algorithm 1 converges to a solution of
VI$(X\times\mathcal{M}_{j},T)$ if at each iteration the projection in line 4
is performed on the (convex) subdomain $\mathcal{M}_{j}\subset\mathcal{M}$,
$j\in\\{1,\ldots,q\\}$. However, we are interested in establishing convergence
on the entire domain $\mathcal{M}$, so at each iteration the projected
solution might belong to a different subdomain $\mathcal{M}_{j}$,
$j\in\\{1,\ldots,q\\}$. This does not trivially follow from the second part of
Lemma 3.2; therefore, we capitalize on Lemmas 3.1 and 3.2 to establish an
additional condition on $\tau$ such that Algorithm 1 retrieves a solution of
VI$(X\times\mathcal{M},T)$. This is formally established in the following
theorem.

######  Theorem 2.

Consider Assumptions 1, 2, 3, and 5. Fixed $0\leq M\leq m$, assume the domain
$\Pi_{K}$ is nonempty for any of the ${m\choose M}$ combinations of
constraints tightened as in (11). Let $D\succ 0$ be defined as in (15), where
$\tau$ satisfies (16) and

$\tau<\frac{-(\bar{C}+\bar{R})+\sqrt{(\bar{C}+\bar{R})^{2}+2\zeta\bar{R}}}{2\bar{R}},$
(17)

where $\bar{R}=\max\big{\\{}\sup_{x\in
X}\sup_{\mu\in\mathcal{M}}\|2A(F(x)+A\raisebox{1.0pt}{${}^{\intercal}$}\mu)\|,\\\
\sup_{x\in X}\|Ax-b\|\big{\\}}$ and $\bar{C}=c\rho\sqrt{m-M}$.

Then Algorithm 1 converges to a solution of VI$(X\times\mathcal{M})$ for any
given multi-sample $\boldsymbol{\delta}_{K}\in\Delta^{K}$ and initial
condition $y^{(0)}\in X\times\mathcal{M}$.

Note that as $\mu^{(\kappa)}\to\mu^{*}$, we have $Q(\mu^{(\kappa)})\to
Q(\mu^{*})=:Q^{*}$. Then, the solution returned by Algorithm 1 is the
equilibrium of a variant of $G_{K}$ with $m-M$ tightened constraints—as it
follows from (11) with $Q(\mu)$ replaced by $Q^{*}$. The next result
accompanies the region $S_{K}^{*}=\Pi_{K}\cap\mathbb{B}({x^{*}},\rho)$ of
strategic deviations from the equilibrium $x^{*}$ with _a priori_
probabilistic feasibility guarantees that can be tuned by means of $M$. It
should be noted that Theorem 2 establishes that there exists a choice of
$\tau$ to guarantee convergence of Algorithm 1. The exact admissible range for
$\tau$—albeit explicit via (16) and (17)—can be difficult to quantify due to
$\bar{R}$. Numerical evidence suggests that selecting a small enough value is
sufficient to achieve convergence.

######  Theorem 3.

Consider Assumptions 1–5. Let $x^{*}$ and
$S_{K}^{*}=\Pi_{K}\cap\mathbb{B}({x^{*}},\rho)$ be returned by Algorithm 1;
fix $\overline{\epsilon}\in(0,1)$ and $M$. We then have that

$\mathbb{P}^{K}\Big{\\{}\boldsymbol{\delta}_{K}\in\Delta^{K}:\;\mathbb{V}(S_{K}^{*})\leq\overline{\epsilon}\Big{\\}}\\\
\geq
1-\sum_{i=0}^{nN+M-1}\binom{K}{i}\overline{\epsilon}^{i}(1-\overline{\epsilon})^{K-i}.$
(18)

By Definition 2, Theorem 3 guarantees that for any point in $S_{K}^{*}$, the
probability of constraint violation is bounded by $\bar{\epsilon}$, with
confidence at least
$1-\sum_{i=0}^{nN+M-1}\binom{K}{i}\overline{\epsilon}^{i}(1-\overline{\epsilon})^{K-i}$.
The dependence of this term on $M$ gives us an additional degree of freedom in
trading the robustness of the solution for its associated probabilistic
confidence. The choice of $M$ can also have an effect on the size of
$S_{K}^{*}$, as well as on the location of $x^{*}$, thus resulting in a trade-
off between performance and robustness.

For the case in which the coupling constraints concern exclusively the
aggregate variable, it can be shown that the support rank for all points in
the region $S_{K}^{*}$ is upper-bounded by $n+M-1$. This allows to state (18)
with a much higher confidence of
$1-\sum_{i=0}^{n+M-1}\binom{K}{i}\overline{\epsilon}^{i}(1-\overline{\epsilon})^{K-i}$;
for details, we refer the reader to [37] and [42].

##  4 Numerical example

Consider a game with $N$ agents whose decisions are subject to deterministic
local constraints and uncertain coupling constraints on the aggregate
decision:

$\left.\begin{aligned} &\min_{x_{i}\in X_{i}}\
x_{i}\raisebox{1.0pt}{${}^{\intercal}$}(C\sigma(x)+d)\\\ &\text{ subject to }\
\underline{b}_{\delta_{k}}\leq\sigma(x)\leq\overline{b}_{\delta_{k}},\\\
&\qquad\qquad\qquad\qquad\qquad k=1,\dots,K\end{aligned}\right\\}\;\forall
i\in\mathcal{N},$ (19)

where $C\succ\alpha I_{n}$, for some $\alpha>0$, and $d\in\mathbb{R}^{n}$. We
impose no knowledge of $\Delta$ and $\mathbb{P}$; we rely instead on a
scenario-based approximation of the game, whereby each sample
$\delta_{k}\in\boldsymbol{\delta}_{K}$ gives rise to
$\underline{b}_{\delta_{k}},\overline{b}_{\delta_{k}}$. Eq. (19) is an
_aggregative_ game in the form of (3).

In this instance, we assume each agent’s action has negligible effect on the
aggregate, and accordingly consider a WE-seeking problem. Following the
definition of $F_{\mathrm{WE}}$ (Sec. 2.2), we get
$F(x)=F_{\mathrm{WE}}=[C\sigma(x)+c]_{i\in\mathcal{N}}$. It can be verified
that $F$ is Lipschitz continuous and strongly monotone with respect to
$\sigma$: by [19, Thm. 2.3.3], (19) admits a unique aggregate equilibrium
$\sigma^{*}=\sigma(x^{*})$.444We note that this case slightly transcends the
conditions in Theorem 2, as $F$ does not comply with Assumption 3-(1).
Convergence of Algorithm 1 (following from the nonexpansiveness of $T$ on each
subdomain $\mathcal{M}_{j}$) can still be ensured on this example due to the
affine structure of $F$; cf. [19, Sec. 12.5.1].

We employ Algorithm 1 to seek a WE $x^{*}$ such that, by fixing $M$, a
prespecified theoretical violation level is guaranteed for the set
$\Pi_{K}\cap\mathbb{B}(x^{*},\rho)$. Due to uniqueness of $\sigma^{*}$, all
sets $\mathbb{B}(x^{*},\rho)$—parametrised by any $x^{*}$ solving (19)—are
projected into the unique ball $\mathbb{B}(\sigma^{*},\rho/N)$ in the
aggregate space. Also note that by definition of $\sigma$, at most $n$ non-
redundant samples will contribute to define the domain
$\Pi^{\sigma}_{K}:=\tfrac{1}{N}(\mathbf{1}_{N}\raisebox{1.0pt}{${}^{\intercal}$}\otimes
I_{n})X\cap\left(\bigcap_{k=1}^{K}C_{\delta_{k}}\right)$ in (19). For the
derivation of the robustness guarantees, we can thus restrict our attention to
$S_{K}^{*}=\Pi^{\sigma}_{K}\cap\mathbb{B}(\sigma^{*},\rho/N)\subseteq\mathbb{R}^{n}$.
As remarked at the end of Section 3.3, we can apply (18) from Theorem 3 with
tighter confidence, since in this case the support rank for all points
belonging to $S_{K}^{*}$ is bounded by $n+M-1$. For the case $n=2,N=50$, and
for different choices of $M$, Figure 3 depicts the iterates
$\\{\sigma(x^{(\kappa)})\\}$, $\kappa=1,2,\ldots$, generated by our algorithm,
i.e., the projection of the algorithm iterates on the space
$\Pi^{\sigma}_{K}$. It can be observed how the region $S_{K}^{*}$ changes as
the value of $M$ is modified.

It is worth noting that in this case $F(x)$ is _integrable_ —this can be
inferred by [19, Thm. 1.3.1] since the Jacobian of the game is symmetric,
i.e., $\nabla_{x}F(x)=\nabla_{x}F(x)\raisebox{1.0pt}{${}^{\intercal}$}$.
Therefore, a WE $x^{*}$ can also be obtained by solving

$\displaystyle\min_{x\in X}\
\sigma(x)\raisebox{1.0pt}{${}^{\intercal}$}C\sigma(x)+d\raisebox{1.0pt}{${}^{\intercal}$}\sigma(x)$
(20) $\displaystyle\text{subject to}\
\underline{b}_{\delta_{k}}\leq\sigma(x)\leq\overline{b}_{\delta_{k}},\;\
k=1,\dots,K.$

In other words, this game admits a _potential function_
$E(x):=\sigma(x)\raisebox{1.0pt}{${}^{\intercal}$}C\sigma(x)+d\raisebox{1.0pt}{${}^{\intercal}$}\sigma(x)$,
whose minimizers correspond to a WE. $E$ can be interpreted as the total cost
incurred by the population of agents, and its minimization leads to the
optimum social welfare. The contour lines of $E$ are depicted in Figure 3:
since $x^{*}$ minimises $E(\cdot)$, $\sigma^{*}$ lies on the contour
associated to the minimum value of $E$ within the feasible domain. Lower
values of $M$ result in larger regions for which guarantees are provided.
Figure 4 shows how the sequence $\\{E(x^{(\kappa)})\\}_{\kappa=1,2,\ldots}$
converges to the minimum _potential_ within the possibly tightened feasibility
region. It can be observed how in this case the efficiency of the equilibrium
decreases as smaller values of $M$ are chosen. The three panels in Figure 4
show the trade-off between system level efficiency and the guaranteed
robustness levels. The lower the value of $M$, the lower the empirical
constraints violation—corresponding to a better confidence bound in the right-
hand side of (18).

(a) $M=0$

Figure 3: Iterates generated by Algorithm 1 (blue diamonds) for different
choices of $M$. In this numerical instance, $N=50$, $\rho=10$, and
$X_{i}:=\\{x_{i}\in\mathbb{R}^{n}:x_{i}\in[\underline{x}_{i},\overline{x}_{i}]\\}$,
with $\underline{x}_{i}=[0\ 0]\raisebox{1.0pt}{${}^{\intercal}$}$,
$\overline{x}_{i}=[3.5\ 3.5]\raisebox{1.0pt}{${}^{\intercal}$}$. The randomly
generated coupling constraints form the rectangular feasibility region
$\Pi^{\sigma}_{K}$ (delineated by the solid black line). The red-shaded region
represents the intersection between the latter and the ball
$\mathbb{B}_{1}(\sigma^{*},\rho/N)$ around the aggregate equilibrium
$\sigma^{*}$ (red diamond marker). In this instance, its volume increases as
larger values for $M$ are chosen. The value associated to the contour lines of
the potential function $E$ decreases from top-right to bottom-left, showing
that $\sigma^{*}$ is the unique minimiser in the admissible region (shaded in
green) after constraint tightening is performed by the algorithm (see Sec.
3.2.2). Figure 4: Potential function $E(x^{(\kappa)})$ evaluated along the
iterations of Algorithm 1. Lower values of $M$ yield better confidence on the
theoretical robustness certificates for the considered region (see Thm. 3),
which results in a lower empirical probability of constraint violation. On the
other hand, the system-level efficiency of the equilibrium increases for
higher values of $M$.

##  5 Concluding remarks

This work proposes a data-driven equilibrium-seeking algorithm such that
probabilistic feasibility guarantees are provided for a region surrounding a
game equilibrium. These guarantees are a priori and the region that is
accompanied with such a probabilistic certificate is tunable. For games that
admit a potential function, the proposed scheme is shown to achieve a trade-
off between cost and the level of probabilistic feasibility guarantees. In
fact, our scheme returns the most efficient equilibrium such that the
predefined guarantees are achieved. Proving this conjecture is left for future
work.

##  6 Appendix

###  6.1 Proof of Lemma 3.1

Let $\mu,z$ be arbitrary vectors in $\mathcal{M}$ and, as in the proof of
Lemma 3.2, define $\vec{\mu},\vec{z}$ as the vectors composed by rearranging
the elements of $\mu,z$ in decreasing order. According to this arrangement,
let $\mathcal{I}_{\mu}=\\{i_{1},i_{2},\dots,i_{m}\\}$ be the ordered set of
indices of $\mu$, i.e., $i_{k}:\,\mu_{i_{k}}=\vec{\mu}_{k}$, $k=1,\ldots,m$;
as a result, $i_{1}$ and $i_{m}$ will be the indices of the largest and
smallest components of $\mu$, respectively. Applying a similar definition to
$z$, we denote the corresponding set
$\mathcal{I}_{z}:=\\{j_{1},j_{2},\dots,j_{m}\\}$. Then, the first $M$ indices
in $\mathcal{I}_{\mu}$ and $\mathcal{I}_{z}$, denoted as $\mathcal{L}_{\mu}$
and $\mathcal{L}_{z}$, respectively, are relative to the constraints not
tightened by the application of $Q(\cdot,M)$. In other words, for all
$\ell\in\mathcal{L}_{\mu}$, $(Q(\mu,M)\boldsymbol{\rho})_{\ell}=0$—and
similarly for $z$. Vice versa, the complementary sets
$\mathcal{L}_{\mu}^{c}=\mathcal{I}_{\mu}\setminus\mathcal{L}_{\mu}$ and
$\mathcal{L}_{z}^{c}=\mathcal{I}_{z}\setminus\mathcal{L}_{z}$ are such that
for all $\ell\in\mathcal{L}_{\mu}^{c}$,
$(Q(\mu,M)\boldsymbol{\rho})_{\ell}=c\rho$, and for all
$\ell\in\mathcal{L}_{z}^{c}$, $(Q(z,M)\boldsymbol{\rho})_{\ell}=c\rho$. Let
$q=[Q(\mu,M)-Q(z,M)]\boldsymbol{\rho}$. We distinguish between the following
cases:

  1. 1.

$\ell\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}$: we have
$(Q(\mu,M)\boldsymbol{\rho})_{\ell}=c\rho$ since
$\ell\in\mathcal{L}_{\mu}^{c}$, while $(Q(z,M)\boldsymbol{\rho})_{\ell}=0$ as
$\ell\in\mathcal{L}_{z}$. Then, $q_{\ell}=c\rho$.

  2. 2.

$\ell\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}$: from
$\ell\in\mathcal{L}_{z}^{c}$ we have $(Q(z,M)\boldsymbol{\rho})_{\ell}=c\rho$.
On the other hand, since $\ell\in\mathcal{L}_{\mu}$,
$(Q(\mu,M)\boldsymbol{\rho})_{\ell}=0$. This results in $q_{\ell}=-c\rho$.

  3. 3.

$\ell\in(\mathcal{L}_{\mu}\cap\mathcal{L}_{z})\cup(\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}^{c})$.
If $\ell\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}$ then
$(Q(\cdot,M)\boldsymbol{\rho})_{\ell}=0$ for both $\mu$ and $z$. Therefore,
$q_{\ell}=0$. Conversely, if
$\ell\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}^{c}$, then
$(Q(\cdot,M)\boldsymbol{\rho})_{\ell}=c\rho$ for both $\mu$ and $z$, which
results again in $q_{\ell}=0$.

The sets $\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}$,
$\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}$,
$(\mathcal{L}_{\mu}\cap\mathcal{L}_{z})\cup(\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}^{c})$
are pairwise disjoint and exhaust the set $\\{1,\dots,m\\}$. Hence we can
write

$\begin{split}&U=(\mu-z)\raisebox{1.0pt}{${}^{\intercal}$}(Q(\mu,M)-Q(z,M))\boldsymbol{\rho}=\sum_{\ell=1}^{m}(\mu_{\ell}-z_{\ell})q_{\ell}\\\
&=\sum_{i\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}}\\!\\!\\!(\mu_{i}-z_{i})c\rho+\\!\\!\\!\sum_{j\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}}\\!\\!\\!(\mu_{j}-z_{j})\cdot(-c\rho)\\\
&=\bigl{(}\\!\\!\underbrace{\sum_{i\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}}\\!\\!\\!\mu_{i}\;-\\!\\!\\!\sum_{j\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}}\\!\\!\\!\mu_{j}}_{=:U_{1}}\;+\\!\\!\\!\underbrace{\sum_{j\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}}\\!\\!\\!z_{j}\;-\\!\\!\\!\sum_{i\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}}\\!\\!\\!z_{i}}_{=:U_{2}}\bigr{)}c\rho,\end{split}$
(21)

Now, notice that for any
$i\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}\subseteq\mathcal{L}_{\mu}^{c}$
and $j\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}\subseteq\mathcal{L}_{\mu}$,
we have by definition of $\mathcal{L}_{\mu}$ and $\mathcal{L}_{\mu}^{c}$ that
$\mu_{i}\leq\mu_{j}$ (which by (13) only holds with equality if
$\mu_{i}=\mu_{j}=0$). With analogous reasoning, we have $z_{i}\geq z_{j}$ for
any $i\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}\subseteq\mathcal{L}_{z}$ and
$j\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}\subseteq\mathcal{L}_{z}^{c}$.
Let $h_{1}$ be the cardinality of the set
$\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}$, and $h_{2}$ that of
$\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}$. Then,

$h_{1}=|\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}|\stackrel{{\scriptstyle(a)}}{{=}}|\mathcal{L}_{z}\setminus\mathcal{L}_{\mu}|=|\mathcal{L}_{z}|-|\mathcal{L}_{\mu}\cap\mathcal{L}_{z}|\\\
\stackrel{{\scriptstyle(b)}}{{=}}|\mathcal{L}_{\mu}|-|\mathcal{L}_{z}\cap\mathcal{L}_{\mu}|=|\mathcal{L}_{\mu}\setminus\mathcal{L}_{z}|=|\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}|=h_{2},$

where $(a)$ holds since
$\mathcal{L}_{\mu},\mathcal{L}_{z}\subseteq\\{1,\ldots,M\\}$, and $(b)$
follows from $|\mathcal{L}_{\mu}|=|\mathcal{L}_{z}|=M$. Therefore
$h_{1}=h_{2}=:h$ and $0\leq h\leq M$, which implies $U_{1}\leq 0$ and
$U_{2}\leq 0$ in (21). We can observe that $U_{1}<0$ and $U_{2}<0$ if
$\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}$ and
$\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}$ are nonempty and the corresponding
components of $\mu$ and $z$ are nonzero. In such a case $h\geq 1$ and we can
write

$U_{1}=\sum_{i\in\mathcal{L}_{\mu}^{c}\cap\mathcal{L}_{z}}\\!\\!\\!\mu_{i}\;-\\!\\!\\!\sum_{j\in\mathcal{L}_{\mu}\cap\mathcal{L}_{z}^{c}}\\!\\!\\!\mu_{j}\leq-h\zeta,$
(22)

where the inequality follows from (13) and the above discussion. A similar
reasoning holds for $U_{2}$. Lastly, note that if $\mu\neq z$ and $h\geq 1$,
then at least one of $U_{1}\leq-h\zeta$ and $U_{2}\leq-h\zeta$ will hold. By
(21), we can thus conclude $U\leq-h\zeta c\rho$ for any $\mu,z\in\mathcal{M}$,
$\mu\neq z$. $\blacksquare$

###  6.2 Proof of Lemma 3.2

Part (1): To prove that the mapping $T$ is continuous on its domain, we first
notice that $T$ is by construction continuous on $X\times\mathcal{M}$ when the
operator $Q(\cdot,M)$ is continuous on $\mathcal{M}$ (as the parameter $M$ is
fixed). Therefore, it is sufficient to show that for any $\mu,z\in\mathcal{M}$
and any $\eta>0$, there exists $\delta>0$ such that

$\|\mu-z\|<\delta~{}\Rightarrow~{}\|Q(\mu,M)-Q(z,M)\|\|\boldsymbol{\rho}\|<\eta,$
(23)

where $\boldsymbol{\rho}=c\rho\mathbf{1}_{m}\neq\mathbf{0}$. To this end,
consider any $\mu,z\in\mathcal{M}$ such that $\|\mu-z\|<\tfrac{\zeta}{2}$,
with $\zeta$ as defined in (13).555The proof of this part also holds for
$\mu\in\mathbb{R}^{m}_{+}\supset\mathcal{M}$. Let $\vec{\mu}$ and $\vec{z}$
denote the vectors $\mu$ and $z$ sorted in decreasing order; thus,
$\vec{\mu}_{\ell}$ is the $\ell$-th largest element of $\mu$ (and similarly
for $z$). For any given $\ell$, let $i:\,\mu_{i}=\vec{\mu}_{\ell}$,
$j:\,z_{j}=\vec{z}_{\ell}$, and
$\bar{\ell}:=\min_{\\{1,\ldots,m\\}}\ell:\,i\neq j$. In words, $\bar{\ell}$ is
the smallest index for which the $\ell$-th largest elements of $\mu$ and $z$
do not appear at the same row of their respective vectors. We then let
$\mathcal{I}$ be the set of indices for which the ordering of the elements of
$\mu$ and $z$ agrees, i.e., for all $k\in\mathcal{I}$, there exists
$\ell<\bar{\ell}$ such that $i=j=k$, with $i:\,\mu_{i}=\vec{\mu}_{\ell}$ and
$j:\,z_{j}=\vec{z}_{\ell}$.

We prove our statement by contradiction. Suppose there exists
$i,j\notin\mathcal{I}$ such that $i:\,\mu_{i}=\vec{\mu}_{\ell}$ and
$j:\,z_{j}=\vec{z}_{\ell}$ for some $\ell>\bar{\ell}$, where $\mu_{i}<\mu_{j}$
and $z_{i}>z_{j}$. First, we note that such an instance exists by hypothesis,
as otherwise the only possible case is where $i=j$, which contradicts
$i,j\notin\mathcal{I}$ and implies $Q(\mu,M)=Q(z,M)$. Since $z\in\mathcal{M}$,
it further holds $z_{j}<z_{i}-\zeta$, which by
$\|\mu_{i}-z_{i}\|\leq\|\mu-z\|<\tfrac{\zeta}{2}$ implies

$z_{j}<\mu_{i}+\frac{\zeta}{2}-\zeta.$ (24)

We bound (24) from below by noting $z_{j}>\mu_{j}-\tfrac{\zeta}{2}$, which
holds since $\|\mu_{j}-z_{j}\|<\tfrac{\zeta}{2}$, obtaining

$\mu_{j}-\frac{\zeta}{2}<\mu_{i}+\frac{\zeta}{2}-\zeta,$

or equivalently $\mu_{j}<\mu_{i}$, which contradicts our hypothesis. Hence the
elements of any pair of vectors $\mu,z\in\mathcal{M}$ such that
$\|\mu-z\|<\tfrac{\zeta}{2}$ must follow the same ordering. By definition of
$P(\cdot)$, this implies $P(\mu)=P(z)$ and, in turn, $\|Q(\mu,M)-Q(z,M)\|=0$.
This validates (23) with $\delta=\tfrac{\zeta}{2}$ and any $\eta>0$,
establishing the continuity of $Q(\cdot,M)$ on $\mathcal{M}$ and concluding
the proof of the first part.

Part (2): We show here that the mapping $T$ fulfils certain nonexpansiveness
properties required for the convergence of Algorithm 1, for compatible choices
of $\tau$. In particular, we provide here a sufficient condition for which the
iteration

$y^{(\kappa+1)}=\text{proj}_{X\times\mathcal{M}_{j},D}\left[y^{(\kappa)}-D^{-1}T(y^{(\kappa)},\rho,M)\right],$
(25)

converges to a solution of VI$(X\times\mathcal{M}_{j},T)$, where
$j\in\\{1,\ldots,q\\}$ is fixed, for any $y^{(0)}\in X\times\mathcal{M}_{j}$.
Notice that in (25) the skew projection is performed on the convex subdomain
$X\times\mathcal{M}_{j}$. Then (25) is the solution of the
VI$(X\times\mathcal{M}_{j},T_{D}^{(\kappa)})$ (see [19, Sec. 12.5.1]), where
$T_{D}^{(\kappa)}(y):=T(y^{(\kappa)},\rho,M)+D(y-y^{(\kappa)})$ is strongly
monotone due to $D\succ 0$ and
$(T(y,\rho,M)-T(y^{\prime},\rho,M))\raisebox{1.0pt}{${}^{\intercal}$}(y-y^{\prime})\geq
0$, for all $y,y^{\prime}\in X\times\mathcal{M}$, which in turn follows from
Assumption 3 and Lemma 3.1. The fixed-point iteration (25) is an instance of
the forward-backward splitting method: we thus resort to standard results in
the literature to prove its convergence. Following the notation in [19, Sec.
12.5.1], let $\tilde{D}:=D_{s}^{{-}1/2}(D-D_{s})D_{s}^{{-}1/2}$, where
$D_{s}:=\frac{D+D\raisebox{1.0pt}{${}^{\intercal}$}}{2}$. Also, let
$\mathcal{U}_{j}:=\\{D_{s}^{{}1/2}y:\,y\in X\times\mathcal{M}_{j}\\}$,
$\mathcal{U}=\bigcup_{j=1}^{q}\mathcal{U}_{j}$, and
$\tilde{T}(w):=D_{s}^{{-}1/2}T(D_{s}^{{-}1/2}w,\rho,M)$, for all
$w\in\mathcal{U}$. To ease notation, we drop the dependence of $\tilde{T}$ and
$\tilde{T}_{D}$ on $\rho,M$, as they remain fixed throughout the proof.
According to [19, Thm. 12.5.2] (see also [46, Sec. 4.3]), to ensure
convergence of (25) to a solution of the VI$(X\times\mathcal{M}_{j},T)$ it is
sufficient to show that $\tilde{T}_{D}=\tilde{T}-\tilde{D}$ is
$\beta$-cocoercive on $\mathcal{U}_{j}$, i.e.,

$(\tilde{T}_{D}(v)-\tilde{T}_{D}(w))\raisebox{1.0pt}{${}^{\intercal}$}(v-w)\geq\beta\|\tilde{T}_{D}(v)-\tilde{T}_{D}(w)\|^{2},$
(26)

for some $\beta>\tfrac{1}{2}$ and all $v,w\in\mathcal{U}_{j}$,
$j\in\\{1,\ldots,q\\}$. In fact, we will go a step forward and demonstrate
here that $\tilde{T}_{D}$ is co-coercive on $\mathcal{U}$ with
$\beta>\tfrac{1}{2}$. Due to the particular saddle problem structure of the
mapping in (8), we adopt the procedure in [19, Prop. 12.5.4] and define $D$ as
in (15) (see also [35]). It then follows from the above definitions that
$\tilde{T}_{D}(w)$, for any $w\in\mathcal{U}$, reduces to

$\tilde{T}_{D}(D_{s}^{{}1/2}y)=D_{s}^{{-}1/2}\begin{bmatrix}F(x)\\\
b-Q(\mu,M)\boldsymbol{\rho}\end{bmatrix},\;\forall y\in X\times\mathcal{M},$
(27)

which can be easily seen by rewriting (8) as

$T(y,\rho,M)=\begin{bmatrix}F(x)\\\
0\end{bmatrix}+\underbrace{\begin{bmatrix}0&A\raisebox{1.0pt}{${}^{\intercal}$}\\\
-A&0\end{bmatrix}}_{D-D_{s}}y+\begin{bmatrix}0\\\
b-Q(\mu,M)\boldsymbol{\rho}\end{bmatrix}.$

Define
$W:=(D_{s}^{{-}1/2})\raisebox{1.0pt}{${}^{\intercal}$}D_{s}^{{-}1/2}=D_{s}^{-1}$,
and let $\vec{Q}(\cdot)$ be a shorthand for $Q(\cdot,M)\boldsymbol{\rho}$ (as
$M$ is a fixed parameter). Then, for any $w_{a},w_{b}\in\mathcal{U}$, we can
expand (26) by using (27), obtaining

$\begin{split}(w_{a}&-w_{b})\raisebox{1.0pt}{${}^{\intercal}$}(\tilde{T}_{D}(w_{a})-\tilde{T}_{D}(w_{b}))-\beta\|\tilde{T}_{D}(w_{a})-\tilde{T}_{D}(w_{b})\|^{2}\\\
={}&(D_{s}^{{-}1/2}w_{a}-D_{s}^{{-}1/2}w_{b})\raisebox{1.0pt}{${}^{\intercal}$}\begin{bmatrix}F(x_{a})-F(x_{b})\\\
\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\end{bmatrix}\\\ &\hskip
60.00009pt-\beta\left\|D_{s}^{{-}1/2}\begin{bmatrix}F(v_{a})-F(v_{b})\\\
\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\end{bmatrix}\right\|^{2}\\\
=&\begin{bmatrix}x_{a}-x_{b}\\\
\mu_{a}-\mu_{b}\end{bmatrix}^{\intercal}\begin{bmatrix}F(x_{a})-F(x_{b})\\\
\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\end{bmatrix}\\\
&-\beta\begin{bmatrix}F(x_{a})-F(x_{b})\\\
\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\end{bmatrix}^{\intercal}W\begin{bmatrix}F(x_{a})-F(x_{b})\\\
\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\end{bmatrix},\end{split}$ (28)

for all $y_{a},y_{b}\in X\times\mathcal{M}$, where the last equality follows
from the definition of $\mathcal{U}_{j}$ and by expanding the norm. Matrix $W$
can be written as $W=\left[\begin{smallmatrix}W_{11}&W_{12}\\\
W_{21}&W_{22}\end{smallmatrix}\right]$, where $W_{11}\in\mathbb{R}^{nN\times
nN}$, $W_{12}\in\mathbb{R}^{nN\times m}$, $W_{33}\in\mathbb{R}^{m\times m}$
are:

$\begin{split}W_{11}&=\tau(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1},\\\
W_{12}&=W_{21}\raisebox{1.0pt}{${}^{\intercal}$}=\tau^{2}(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}A\raisebox{1.0pt}{${}^{\intercal}$},\\\
W_{22}&=\tau
I_{m}+\tau^{3}A(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}A\raisebox{1.0pt}{${}^{\intercal}$}.\end{split}$

Expanding the inner product in (28) with respect to the matrix blocks
$W_{11},W_{12},W_{21},W_{33}$ we obtain

$\begin{split}\beta&(F(x_{a})-F(x_{b}))\raisebox{1.0pt}{${}^{\intercal}$}\big{[}\tfrac{1}{\beta}(x_{a}-x_{b})\\\
&-W_{11}(F(x_{a})-F(x_{b}))-2W_{12}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{1}))\big{]}\\\
&+\beta(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\raisebox{1.0pt}{${}^{\intercal}$}\big{[}\tfrac{1}{\beta}(\mu_{a}-\mu_{b})\\\
&\hskip 60.00009pt-W_{22}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\big{]}\\\
={}&(F(x_{a})-F(x_{b}))\raisebox{1.0pt}{${}^{\intercal}$}(x_{a}-x_{b})\\\
&-\beta(F(x_{a})-F(x_{b}))\raisebox{1.0pt}{${}^{\intercal}$}W_{11}(F(x_{a})-F(x_{b}))\\\
&-2\beta(F(x_{a})-F(x_{b}))\raisebox{1.0pt}{${}^{\intercal}$}W_{12}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\\\
&+(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\raisebox{1.0pt}{${}^{\intercal}$}(\mu_{a}-\mu_{b})\\\
&-\beta(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\raisebox{1.0pt}{${}^{\intercal}$}W_{22}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})).\end{split}$

Setting
$p_{\tau}:=(I-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1/2}(F(x_{a})-F(x_{b}))$
and
$q_{\tau}:=\tau(I-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1/2}A\raisebox{1.0pt}{${}^{\intercal}$}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))$
above we obtain

$\begin{split}(&F(x_{a})-F(x_{b}))\raisebox{1.0pt}{${}^{\intercal}$}(x_{a}-x_{b})\\\
&+(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\raisebox{1.0pt}{${}^{\intercal}$}(\mu_{a}-\mu_{b})\\\
&-\beta\tau(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\raisebox{1.0pt}{${}^{\intercal}$}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\\\
&-\beta\tau(p_{\tau}+q_{\tau})\raisebox{1.0pt}{${}^{\intercal}$}(p_{\tau}+q_{\tau})\\\
&\geq\alpha\|x_{a}-x_{b}\|^{2}+2h\zeta c\rho\\\ &\hskip 30.00005pt-2\beta\tau
h(c\rho)^{2}-2\beta\tau(p_{\tau}\raisebox{1.0pt}{${}^{\intercal}$}p_{\tau}+q_{\tau}\raisebox{1.0pt}{${}^{\intercal}$}q_{\tau}),\end{split}$
(29)

where for the last inequality we used, in order, $(i)$ strong monotonicity of
$F$ (cf. Assumption 3), $(ii)$ Lemma 3.1, $(iii)$
$\|\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\|^{2}\leq 2h(c\rho)^{2}$—which follows
from the same arguments used in the proof of Lemma 3.1—and $(iv)$
$(p_{\tau}+q_{\tau})\raisebox{1.0pt}{${}^{\intercal}$}(p_{\tau}+q_{\tau})\leq
2(p_{\tau}\raisebox{1.0pt}{${}^{\intercal}$}p_{\tau}+q_{\tau}\raisebox{1.0pt}{${}^{\intercal}$}q_{\tau})$.
Expanding the term containing $p_{\tau},q_{\tau}$ in (29) we get

$\begin{split}&\alpha\|x_{a}-x_{b}\|^{2}+2h\zeta c\rho-2\beta\tau
h(c\rho)^{2}\\\
&-2\beta\tau(F(x_{a})\\!-\\!F(x_{b}))\raisebox{1.0pt}{${}^{\intercal}$}(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}(F(x_{a})\\!-\\!F(x_{b}))\\\
&-2\beta\tau^{3}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\raisebox{1.0pt}{${}^{\intercal}$}\\\
&\hskip 70.0001pt\cdot
A(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}A\raisebox{1.0pt}{${}^{\intercal}$}(\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a}))\\\
&\overset{(a)}{\geq}\alpha\|x_{a}-x_{b}\|^{2}+2h\zeta c\rho-2\beta\tau
h(c\rho)^{2}\\\
&\quad-2\beta\tau\|F(x_{a})\\!-\\!F(x_{b})\|^{2}\cdot\|(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}\|\\\
&\quad-2\beta\tau^{3}\|\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\|^{2}\cdot\|(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}\|\cdot\|A\|^{2}\\\
&\overset{(b)}{\geq}(\alpha-2\beta\tau
L_{F}^{2}\|(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}\|)\|x_{a}-x_{b}\|^{2}\\\
&\quad+2h\zeta c\rho-2\beta\tau
h(c\rho)^{2}\bigg{(}1+\frac{2\tau^{2}}{1-\tau^{2}\|A\|^{2}}\|A\|^{2}\bigg{)},\end{split}$
(30)

where $(a)$ is obtained by applying the Cauchy-Schwarz inequality, and in
$(b)$ we use the Lipschitz continuity of $F$ (cf. Assum. 3),
$\|\vec{Q}(\mu_{b})-\vec{Q}(\mu_{a})\|^{2}\leq 2h(c\rho)^{2}$, and the
triangle inequality. Now, notice that for the last term in (30),

$2\beta\tau
h(c\rho)^{2}\bigg{(}1+\frac{2\tau^{2}}{1-\tau^{2}\|A\|^{2}}\|A\|^{2}\bigg{)}\\\
=2\beta\tau h(c\rho)^{2}\frac{1+\tau^{2}\|A\|^{2}}{1-\tau^{2}\|A\|^{2}}\leq
2\beta\tau h(c\rho)^{2}\frac{1+\|A\|^{2}}{1-\tau^{2}\|A\|^{2}},$ (31)

holds for any choice of
$\tau\in\big{(}0,\max\big{\\{}\tfrac{1}{\|A\|},1\big{\\}}\big{)}$. Recall that
by invoking [19, Thm. 12.5.2], our objective is to show that (26) holds for
some $\tau>0$ and $\beta>\tfrac{1}{2}$. Then, by inspecting (30) and using
(31), to achieve this it is sufficient to guarantee

$\displaystyle\alpha-\tau
L_{F}^{2}\|(I_{n}-\tau^{2}A\raisebox{1.0pt}{${}^{\intercal}$}A)^{-1}\|>0,$
$\displaystyle 2h\zeta c\rho-\tau
h(c\rho)^{2}\frac{1+\|A\|^{2}}{1-\tau^{2}\|A\|^{2}}>0,\;\text{if $1\leq h\leq
M$}.$

Solving the quadratic expressions above with respect to $\tau$ results in the
admissible range of values in (16) (these are also satisfying
$\tau\in\big{(}0,\max\big{\\{}\tfrac{1}{\|A\|},1\big{\\}}\big{)}$, required
for (31) to hold). Therefore, for any $\tau$ satisfying this condition,
$\tilde{T}_{D}$ is co-coercive with $\beta>\tfrac{1}{2}$ on the entire domain
$\mathcal{U}$, which in turn implies that co-coercivity of $\tilde{T}_{D}$
holds on each subdomain $\mathcal{U}_{j}$, $j=1,\ldots,q$, with the same
modulus. By [19, Thm. 12.5.2], this is sufficient to guarantee the convergence
of (25) to a solution of the VI$(X\times\mathcal{M}_{j},T)$, thus concluding
the proof. $\blacksquare$

###  6.3 Proof of Theorem 2

Fix any $\tau$ satisfying the conditions of Lemma 3.2 and (17). The sequence
of iterates $\\{y^{(\kappa)}\\}_{\kappa=1,2,\ldots}$ (where
$y^{(\kappa)}=(x^{(\kappa)},\mu^{(\kappa)})$) generated by Algorithm 1 lives
in a bounded set since $X$ and $\mathcal{M}$ are assumed to be compact (see
Assumption 5). As such there exist convergent subsequences, or in other words,
the set

$\Omega:=\Big{\\{}\bar{y}=(\bar{x},\bar{\mu})\colon\;\exists\text{ subsequence
}\\{\kappa_{i}\\}_{i\in\mathbb{N}}\\\ \text{such that
}\lim_{i\to\infty}\kappa_{i}=\infty,\,\lim_{i\to\infty}y^{(\kappa_{i})}=\bar{y}\Big{\\}},$
(32)

containing the limit points of $\\{y^{(\kappa)}\\}$ is non-empty; see, e.g.,
[40, p. 48]. In particular, we will show that $\Omega$ is a singleton for any
$\tau$ satisfying (16)–(17), which implies that the iterates generated by
Algorithm 1 have a unique limit point, hence they converge. To achieve this,
we assume for the sake of contradiction that there exist two cluster points
$\bar{y}_{1},\bar{y}_{2}\in\Omega$, where
$\bar{y}_{1}=(\bar{x}_{1},\bar{\mu}_{1})$ and
$\bar{y}_{2}=(\bar{x}_{2},\bar{\mu}_{2})$. Moreover, we assume that
$\bar{\mu}_{1}\in\mathcal{M}_{i}$, and $\bar{\mu}_{2}\in\mathcal{M}_{j}$, with
$i\neq j$. Note that if this were not the case, then we would be in a trivial
case where $\bar{y}_{1}=\bar{y}_{2}$, due to co-coercivity of $T$ (see Lemma
3.2)—by which Algorithm 1 converges to a unique solution when restricted to
any convex subdomain $X\times\mathcal{M}_{j}$, $j=1,\ldots,q$. To ease the
notation in the remainder of the proof, we assume without loss of generality
that $\bar{\mu}_{1}\in\mathcal{M}_{1}$, $\bar{\mu}_{2}\in\mathcal{M}_{2}$ (see
Fig. 5).

Figure 5: Domain $\mathcal{M}$ of the Lagrange multipliers associated to the
coupling constraints, for the case $m=2$. Notice the minimum distance $\zeta$
between any two subdomains of $\mathcal{M}$ involves the origin as one of the
subdomains.

By (32) there exist an infinite-length subsequence
$\\{\kappa_{i}\\}_{i\in\mathbb{N}}$ of the iterates generated by Algorithm 1
whose elements get arbitrarily close to $\bar{\mu}_{1}$ while staying in
$\mathcal{M}_{1}$ where this cluster point belongs (similarly for
$\bar{\mu}_{2}$). We then have that for any $\delta>0$, there exists
$\tilde{\kappa}$ such that for all $\kappa_{i}\geq\tilde{\kappa}$,
$\|y^{(\kappa_{i})}-\bar{y}_{1}\|\leq\delta$; this implies
$\|x^{(\kappa_{i})}-\bar{x}_{1}\|\leq\delta$ and
$\|\mu^{(\kappa_{i})}-\bar{\mu}_{1}\|\leq\delta$.

Due to our contradiction hypothesis (recall that
$\\{\kappa_{i}\\}_{i\in\mathbb{N}}$ is a subsequence), the sequence of
iterates generated by Algorithm 1 would be leaving $\mathcal{M}_{1}$ towards
$\mathcal{M}_{2}$ infinitely often. Denote then by
$\overline{\kappa}>\tilde{\kappa}$ the smallest index of the subsequence such
that $\mu^{(\overline{\kappa})}\in\mathcal{M}_{1}$, but
$\mu^{(\overline{\kappa}+1)}\in\mathcal{M}_{2}$, i.e., after the
$\overline{\kappa}$-th iterate the original sequence would jump to
$\mathcal{M}_{2}$ for the first time after $\tilde{\kappa}$.

For this jump to occur, the unprojected solution for the Lagrange multipliers
must be “closer” to $\mathcal{M}_{2}$ than to any other sub-domain of
$\mathcal{M}$. To see this more formally, let $D_{\mu}^{-1}$ denote the lower
block-row of $D^{-1}=\left[\begin{smallmatrix}\tau I_{nN}&0\\\ 2A\tau^{2}&\tau
I_{m}\end{smallmatrix}\right]$, corresponding to the Lagrange multiplier
update in line 4 of Algorithm 1. By definition of $\mathcal{M}$, such a jump
requires the distance between the unprojected gradient step at
$\overline{\kappa}+1$ and $\mu^{(\overline{\kappa})}$ to satisfy

$\|\mu^{(\overline{\kappa})}-D_{\mu}^{-1}T(y^{(\overline{\kappa})},\rho,M)-\mu^{(\overline{\kappa})}\|>\zeta/2,$
(33)

where $\zeta$ is the worst-case (minimum) distance between any two subdomains
of $\mathcal{M}$—corresponding to the case where one of these is the origin
(in this case $\mathcal{M}_{2}$). Figure 5 provides a pictorial illustration
of this construction. However, we have that

$\displaystyle\|$
$\displaystyle\mu^{(\overline{\kappa})}-D_{\mu}^{-1}T(y^{(\overline{\kappa})},\rho,M)-\mu^{(\overline{\kappa})}\|$
$\displaystyle=\tau\|-2\tau
A(F(x^{(\overline{\kappa})})+A\raisebox{1.0pt}{${}^{\intercal}$}\mu^{(\overline{\kappa})})$
$\displaystyle\hskip
11.38092pt+Ax^{(\overline{\kappa})}-b+Q(\mu^{(\overline{\kappa})},M)\boldsymbol{\rho}\|$
$\displaystyle=\tau\|-2\tau
A(F(x^{(\overline{\kappa})})-F(\bar{x}_{1})+A\raisebox{1.0pt}{${}^{\intercal}$}(\mu^{(\overline{\kappa})}-\bar{\mu}_{1}))$
$\displaystyle\hskip
11.38092pt+A(x^{(\overline{\kappa})}-\bar{x}_{1})+Q(\mu^{(\overline{\kappa})},M)\boldsymbol{\rho}$
$\displaystyle\hskip 11.38092pt-2\tau
A(F(\bar{x}_{1})+A\raisebox{1.0pt}{${}^{\intercal}$}\bar{\mu}_{1})+(A\bar{x}_{1}-b)\|$
$\displaystyle\leq\tau^{2}\|2A(F(\bar{x}_{1})+A\raisebox{1.0pt}{${}^{\intercal}$}\bar{\mu}_{1})\|+\tau\|A\bar{x}_{1}-b\|$
$\displaystyle\hskip
11.38092pt+\tau\|Q(\mu^{(\overline{\kappa})},M)\boldsymbol{\rho}\|+\tau\|A\|\|x^{(\overline{\kappa})}-\bar{x}_{1}\|$
$\displaystyle\hskip
11.38092pt+2\tau^{2}\big{(}\|A(F(x^{(\overline{\kappa})})-F(\bar{x}_{1}))\|+\|AA\raisebox{1.0pt}{${}^{\intercal}$}(\mu^{(\overline{\kappa})}-\bar{\mu}_{1})\|\big{)}$
$\displaystyle\leq(\tau^{2}+\tau)\bar{R}+\tau c\rho\sqrt{m-M}$
$\displaystyle\hskip
11.38092pt+\tau\delta\big{(}2\tau(L_{F}\|A\|+\|AA\raisebox{1.0pt}{${}^{\intercal}$}\|)+\|A\|\big{)}$
(34)

where the first equality follows from the definition of $D_{\mu}^{-1}$ and
$T$, and the second one by adding and subtracting $F(\bar{x}_{1})$,
$A\raisebox{1.0pt}{${}^{\intercal}$}\bar{\mu}_{1}$ and $A\bar{x}_{1}$. The
first inequality is due to the triangle inequality, while the last one follows
from the previous one by upper-bounding (i) the first two terms using the
definition of $\bar{R}$; (ii)
$\|Q(\mu^{(\overline{\kappa})},M)\boldsymbol{\rho}\|$ by $c\rho\sqrt{m-M}$
based on its definition; and (iii) the last three terms using
$\|F(x^{(\overline{\kappa})})-F(\bar{x}_{1})\|\leq
L_{F}\|x^{(\overline{\kappa})}-\bar{x}_{1}\|$ by Assumption 3, and
$\|x^{(\overline{\kappa})}-\bar{x}_{1}\|\leq\delta$,
$\|\mu^{(\overline{\kappa})}-\bar{\mu}_{1}\|\leq\delta$.

By (34), and choosing $\tau$ as in (17), we have that

$\displaystyle\|\mu^{(\overline{\kappa})}-D_{\mu}^{-1}T(y^{(\overline{\kappa})},\rho,M)-\mu^{(\overline{\kappa})}\|<\frac{\zeta}{2}+\bar{K}\delta$
(35)

where $\bar{K}$ is a constant, emanating from the coefficient of $\delta$ in
(34) when substituting for $\tau$ the upper-bound in (17). Since $\delta$ is
arbitrary, taking $\lim\sup_{\delta\to 0}$ in (35) establishes a contradiction
with (33). The latter implies then that $\bar{\mu}_{2}$ must belong to the
same subdomain $\mathcal{M}_{1}$ with $\bar{\mu}_{1}$, i.e., all cluster
points should be in the same subdomain of $\mathcal{M}$. However, by Lemma
3.2, due to co-coercivity of $T$ on each subdomain $X\times\mathcal{M}_{j}$,
$j=1,\ldots,q$, we can only have one cluster point, i.e., $\Omega$ is a
singleton, implying that Algorithm 1 converges, thus concluding the proof.
$\blacksquare$

###  6.4 Proof of Theorem 3

The elements of the minimal compression set $I$ of Algorithm 1 can belong to
one or both of the following subsets:

  1. 1.

The subset $I_{1}$ of samples that support the solution $x^{*}$. Note that
since Algorithm 1 converges to the point $(x^{*},\mu^{*})$ for a fixed choice
of $M$, $Q(\mu^{*},M)$ will be a fixed quantity. As such, Algorithm 1 will
converge to a solution of

$\displaystyle\ \text{Find}\ x^{*}\in\widehat{\Pi}_{K}\ \text{such that}$
$\displaystyle F(x^{*})\raisebox{1.0pt}{${}^{\intercal}$}(x-x^{*})\geq 0\
\text{for any}\ x\in\widehat{\Pi}_{K},$ (36)

where $\widehat{\Pi}_{K}$ denotes the polytope obtained from $\Pi_{K}$ by
tightening at most $M$ coupling constraints, as dictated by (11) with
$Q(\mu^{*},M)$. The constraints in (36) are equivalent to
$F(x^{*})\raisebox{1.0pt}{${}^{\intercal}$}x\geq
F(x^{*})\raisebox{1.0pt}{${}^{\intercal}$}x^{*}$ for all
$x\in\widehat{\Pi}_{K}$. Then, $x^{*}$ is the minimiser of

$\displaystyle\min_{x\in\mathbb{R}^{nN}}F(x^{*})\raisebox{1.0pt}{${}^{\intercal}$}x$
$\displaystyle\text{subject to }x\in\widehat{\Pi}_{K},$ (37)

which is unique due to Lemma 1. Since the cost function is linear in $x$ and
$\widehat{\Pi}_{K}$ is convex by Assumption 2, we obtain a scenario program as
in [10]. Applying [10, Prop. 1] to (37), we have that $|I_{1}|\leq nN$, i.e.,
the number of support samples of $x^{*}$ is bounded by the dimension of the
decision vector $nN$.

  2. 2.

The subset $I_{2}$ of samples whose corresponding coupling constraints
intersect $\mathbb{B}(x^{*},\rho)$. By construction of Algorithm 1 we have
that $|I_{2}|\leq M$.

As such, we have that $I=I_{1}\cup I_{2}$ is a compression set with
cardinality $|I|=|I_{1}\cup I_{2}|\leq|I_{1}|+|I_{2}|\leq nN+M$. Then, from
Corollary 2 in [32] it follows that

$\displaystyle\mathbb{P}^{K}\Big{\\{}\delta_{K}\in\Delta^{K}$
$\displaystyle:~{}\mathbb{V}(S_{K}^{*})>\overline{\epsilon}\Big{\\}}$
$\displaystyle\leq\sum_{i=0}^{nN+M-1}{K\choose
i}\overline{\epsilon}^{i}(1-\overline{\epsilon})^{K-i},$ (38)

which concludes the proof. $\blacksquare$

## References

  * [1] D. Acemoglu and M. K. Jensen.  Aggregate comparative statics.  Games and Economic Behavior, 81:27 – 49, 2013.
  * [2] D. Acemoglu and A. Ozdaglar.  Opinion dynamics and learning in social networks.  Dynamic Games and Applications, 1(1):3–49, 2011.
  * [3] M. Aghassi and D. Bertsimas.  Robust game theory.  Math. Program., 107(1-2):231–273, 2006.
  * [4] T. Alamo, V. Mirasierra, F. Dabbene, and M. Lorenzen.  Safe approximations of chance constrained sets by probabilistic scaling.  2019 18th European Control Conference (ECC), pages 1380–1385, 2019\.
  * [5] T. Alamo, R. Tempo, and E. Camacho.  Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems.  IEEE Transactions on Automatic Control, 54(11):2545–2559, 2009\.
  * [6] M. J. Beckmann, C. B. McGuire, and C. B. Winsten.  Studies in the economics of transportation.  Yale University Press, 1956.
  * [7] G. C. Calafiore.  Random convex programs.  SIAM Journal on Optimization, 20(6):3427–3464, 2010.
  * [8] G. C. Calafiore and M. C. Campi.  The scenario approach to robust control design.  IEEE Transactions on Automatic Control, 51(5):742–753, 2006.
  * [9] C. F. Camerer.  Behavioral Game Theory: Experiments in Strategic Interaction.  Princeton University Press, 2003.
  * [10] M. C. Campi and S. Garatti.  The exact feasibility of randomized solutions of uncertain convex programs.  SIAM Journal on Optimization, 19(3):1211–1230, 2008.
  * [11] M. C. Campi and S. Garatti.  Wait-and-judge scenario optimization.  Mathematical Programming, 167(1):155–189, Jan 2018.
  * [12] M. C. Campi, S. Garatti, and M. Prandini.  The scenario approach for systems and control design.  Annual Reviews in Control, 33(2):149–157, 2009.
  * [13] M. C. Campi, S. Garatti, and F. A. Ramponi.  A general scenario theory for nonconvex optimization and decision making.  IEEE Transactions on Automatic Control, 63(12):4067 – 4078, 2018\.
  * [14] P. Couchman, B. Kouvaritakis, M. Cannon, and F. Prashad.  Gaming strategy for electric power with random demand.  IEEE Transactions on Power Systems, 20(3):1283–1292, 2005.
  * [15] D. P. de Farias and B. V. Roy.  The linear programming approach to approximate dynamic programming.  Operations Research, 51:850–865, 2003.
  * [16] F. Fabiani, K. Margellos, and P. J. Goulart.  On the robustness of equilibria in generalized aggregative games.  2020 59th IEEE Conference on Decision and Control (CDC), pages 3725–3730, 2020.
  * [17] F. Fabiani, K. Margellos, and P. J. Goulart.  Probabilistic feasibility guarantees for solution sets to uncertain variational inequalities.  Automatica, 137:110120, 2022.
  * [18] F. Facchinei and C. Kanzow.  Generalized Nash equilibrium problems.  Annals of Operations Research, 175(1):177–211, 2010.
  * [19] F. Facchinei and J.-S. Pang.  Finite-Dimensional Variational Inequalities and Complementarity Problems.  Springer-Verlag New York, 2003.
  * [20] F. Facchinei, V. Piccialli, and M. Sciandrone.  Decomposition algorithms for generalized potential games.  Computational Optimization and Applications, 50(2):237–262, 2011\.
  * [21] F. Fele and K. Margellos.  Probabilistic sensitivity of Nash equilibria in multi-agent games: a wait-and-judge approach.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 5026–5031, Dec 2019.
  * [22] F. Fele and K. Margellos.  Probably approximately correct Nash equilibrium learning.  IEEE Transactions on Automatic Control, 66(9):4238–4245, 2021.
  * [23] S. Grammatico, F. Parise, M. Colombino, and J. Lygeros.  Decentralized convergence to Nash equilibria in constrained deterministic mean field control.  IEEE Transactions on Automatic Control, 61(11):3315–3329, 2016\.
  * [24] P. T. Harker.  Generalized Nash games and quasi-variational inequalities.  European Journal of Operational Research, 54(1):81–94, 1991.
  * [25] S. Hayashi, N. Yamashita, and M. Fukushima.  Robust Nash equilibria and second-order cone complementarity problems.  Journal of Nonlinear and Convex Analysis, 6, 2005.
  * [26] M. K. Jensen.  Aggregative games and best-reply potentials.  Economic Theory, 43(1):45–66, 2010.
  * [27] D. Kahneman and A. Tversky.  Prospect theory: An analysis of decision under risk.  Econometrica, 47(2):263–291, 1979.
  * [28] N. S. Kukushkin.  Best response dynamics in finite games with additive aggregation.  Games and Economic Behavior, 48(1):94–110, 2004.
  * [29] A. A. Kulkarni and U. V. Shanbhag.  On the variational equilibrium as a refinement of the generalized Nash equilibrium.  Automatica, 48(1):45–55, 2012.
  * [30] F. Lamnabhi-Lagarrigue, A. Annaswamy, S. Engell, A. Isaksson, P. Khargonekar, R. M. Murray, H. Nijmeijer, T. Samad, D. Tilbury, and P. Van den Hof.  Systems & control for the future of humanity, research agenda: Current and future roles, impact and grand challenges.  Annual Reviews in Control, 43:1–64, 2017.
  * [31] M. Mammarella, V. Mirasierra, M. Lorenzen, T. Alamo, and F. Dabbene.  Chance-constrained sets approximation: A probabilistic scaling approach.  Automatica, 137:110–108, 2022.
  * [32] K. Margellos, M. Prandini, and J. Lygeros.  On the connection between compression learning and scenario based single-stage and cascading optimization problems.  IEEE Transactions on Automatic Control, 60(10):2716–2721, 2015\.
  * [33] J. F. Nash.  Equilibrium points in $N$-person games.  Proceedings of the National Academy of Sciences of the United States of America, 36(48-49), 1950.
  * [34] D. Paccagnan and M. C. Campi.  The scenario approach meets uncertain game theory and variational inequalities.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 6124–6129, 2019.
  * [35] 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, 64(4):1373–1388, April 2019\.
  * [36] G. Pantazis, F. Fele, and K. Margellos.  A posteriori probabilistic feasibility guarantees for Nash equilibria in uncertain multi-agent games.  IFAC-PapersOnLine, 53(2):3403–3408, 2020.  21st IFAC World Congress.
  * [37] G. Pantazis, F. Fele, and K. Margellos.  On the probabilistic feasibility of solutions in multi-agent optimization problems under uncertainty.  European Journal of Control, 63:186–195, 2022.
  * [38] C. R. Plott and V. L. Smith, editors.  Handbook of Experimental Economics Results, volume 1.  Elsevier, 2008.
  * [39] A. Rubinstein.  Modeling Bounded Rationality.  MIT Press, 1998.
  * [40] W. Rudin.  Principles of mathematical analysis.  McGraw-Hill, Inc., 3rd edition, 1976.
  * [41] W. Saad, Z. Han, H. V. Poor, and T. Basar.  Game-theoretic methods for the smart grid: An overview of microgrid systems, demand-side management, and smart grid communications.  IEEE Signal Processing Magazine, 29(5):86–105, 2012.
  * [42] G. Schildbach, L. Fagiano, and M. Morari.  Randomized solutions to convex programs with multiple chance constraints.  SIAM Journal on Optimization, 23, May 2012.
  * [43] G. Scutari, F. Facchinei, J. Pang, and D. P. Palomar.  Real and complex monotone communication games.  IEEE Transactions on Information Theory, 60(7):4197–4231, July 2014\.
  * [44] V. V. Singh, O. Jouini, and A. Lisser.  Existence of Nash equilibrium for chance-constrained games.  Operations Research Letters, 44(5):640 – 644, 2016.
  * [45] Y. Wang, W. Saad, N. B. Mandayam, and H. V. Poor.  Load shifting in the smart grid: To participate or not?  IEEE Transactions on Smart Grid, 7(6):2604–2614, 2016.
  * [46] D. L. Zhu and P. Marcotte.  Co-coercivity and its role in the convergence of iterative schemes for solving variational inequalities.  SIAM Journal on Optimization, 6(3):714–726, 1996.