# Bandit Online Learning of Nash Equilibria in Monotone Games

Tatiana Tatarenko and Maryam Kamgarpour, IEEE Member T. Tatarenko
<EMAIL_ADDRESS>is with the Control Methods and
Robotics Lab Technical University Darmstadt, Darmstadt, Germany 64283, M.
Kamgarpour<EMAIL_ADDRESS>is with Electrical and Computer Engineering,
UBC, Vancouver, Canada. M. Kamgarpour gratefully acknowledges ERC Starting
Grant CONENE.

###### Abstract

We address online bandit learning of Nash equilibria in multi-agent convex
games. We propose an algorithm whereby each agent uses only obtained values of
her cost function at each joint played action, lacking any information of the
functional form of her cost or other agents’ costs or strategies. In contrast
to past work where convergent algorithms required strong monotonicity, we
prove that the algorithm converges to a Nash equilibrium under mere
monotonicity assumption. The proposed algorithm extends the applicability of
bandit learning in several games including zero-sum convex games with possibly
unbounded action spaces, mixed extension of finite-action zero-sum games, as
well as convex games with linear coupling constraints.

###### Index Terms:

learning in games, monotone bandit learning

##  I Introduction

Game theory is a powerful framework to address optimization and learning of
multiple interacting agents referred to as players. Such multi-agent problems
arise in several application domains including traffic networks, internet,
auctions and adversarial learning. In a multi-agent setting, the notion of
Nash equilibrium captures a desirable solution as it exhibits stability.
Namely, no player has an incentive to unilaterally deviate from this solution.
Nash equilibrium is also consistent with a notion referred to as rationality -
each player aims to optimize her own cost function. Given these theoretical
justifications for Nash equilibria, from the viewpoint of learning the
question is whether players can learn their Nash equilibrium strategies with
limited information about the game. In particular, in several application
domains, each player might not know the functional form of her objective. For
example, the travel times of edges in a traffic network or the market
constraints in an auction are unknown a priori and depend in non-trivial ways
on actions and objectives of other players. However, by playing the game, a
player can have access to the so-called online bandit information feedback,
namely, she can receive payoff information of her objective (zeroth-order
oracle) for any feasible joint actions taken by all the players. Thus, we
focus on how do players learn Nash equilibria given online bandit information.

Bandit learning in games has been mainly explored in finite action settings.
It is known that if each player uses a no-regret algorithm, the time average
of the sequence of played actions converges to a coarse-correlated equilibrium
- a relaxed notion of equilibrium which encompasses Nash equilibria as well as
possibly non-rationalizable strategies. The convergence of time-averaged
sequence of plays in no-regret algorithms to a mixed strategy Nash equilibrium
can be established for special games such as two-player zero-sum game.
However, the convergence of time-averaged actions does not imply the actual
sequence of plays also converges to a Nash equilibrium. This issue was
recently re-examined in light of progresses in bandit convex optimization. In
particular, [13] showed that if each player uses a fairly general model of no-
regret algorithm (FTRL) in continuous-time, the actual sequence of plays may
not converge to a mixed strategy Nash equilibrium in a two-player zero-sum
game and it might exhibit a non-vanishing cycle111The FTRL algorithms explored
had access to more than bandit feedback as each player could observe the
gradient of its mixed extension cost function at each time step.. The work of
[2] analyzed this non-convergence through the lens of Hamiltonian and
potential of a game and proposed convergent algorithms assuming access to
exact first and second order oracles (gradients and Hessians) of the players’
objectives. More recently, motivated by the Hamiltonian analogy of zero-sum
games, [1] extended the analysis of [13] to discrete-time setting and showed
that sequential gradient descent can overcome the divergence of discretized
simultaneous gradient descent in these games.

The mixed extension of a finite action game falls into category of convex
games. In such games, each player’s objective is convex in her own decision
variable for any fixed action of other players (in mixed extension of finite
action games, this objective is linear in each player’s action). Furthermore,
the strategy sets are convex (simplexes in mixed extension games). It is known
that the Nash equilibria in convex differentiable games correspond to the
solution set of a variational inequality problem. Hence, one can use
algorithms for solving variational inequalities to compute Nash equilibria.
This connection has been used in two lines of works addressing bandit learning
of Nash equilibria in convex games. The works in [18, 4] showed that no-regret
learning can converge to Nash equilibria in certain class of convex games.
Both works leveraged the idea of one-point estimation of the gradient of the
player’s cost function using the bandit payoff information. Whereas [4]
performed this gradient estimation by perturbing the played actions with a
point sampled from the uniform distribution on the unit sphere, motivated by
Stoke’s theorem [6, 9], the approach in [18] formed these gradient estimates
using Normal distribution and smoothing motivated by [20, 14]. In both cases,
the convergence was proven by appropriately choosing the stepsizes to tradeoff
the bias and variance of the resulting estimation terms in the stochastic
approximation procedure.

The above works on learning Nash equilibria rely on strict monotonicity of the
game mapping. The game mapping is the stacked vector of gradient of each
player’s objective with respect to her actions. In case the game is potential,
the game mapping is symmetric and thus, the game mapping corresponds to
gradient of a single function, the so-called potential function of the game.
Hence, finding equilibria can be cast as a single-objective optimization and
bandit algorithms can be readily applied to learn Nash equilibria. However, in
general noncooperative games, this game mapping need not be symmetric. In such
cases, a necessary and sufficient condition for convergence of no-regret
algorithms such as those explored in [18, 4, 13] is strict monotonicity of the
game mapping.

The strict monotonicity of the game mapping is a stringent condition. The non-
convergence issue and in particular cyclic behavior discussed in [13, 2, 1] is
due to the fact that the game mapping in several games do not satisfy strict
monotonicity. Zero-sum games for example are only monotone. This class of
games have been widely adopted in robust optimization and control, in models
of perfect competition and more recently in adversarial training and deep
learning. In addition to zero-sum games, generalized Nash equilibria problems,
that is, Nash equilibria with coupling constraints can at best have monotone
extended game mappings due to the coupling constraints. The generalized Nash
equilibria problems arise in several domains where a resource must be shared
between agents and this can only be formulated as a hard rather than soft
constraints [5]. The relevance of games with merely monotone mapping, motivate
our paper on learning Nash equilibria under bandit feedback in this game
class.

Given the connection of Nash equilibria of convex games with solution set of
variational inequality (VI) problem, a natural starting point for learning
Nash equilibria of merely monotone games is to search for algorithms for
solving the corresponding VIs. This topic is well-explored in [16, Chapter 12]
and several approaches including extra-gradient and Tikhonov regularization
for finding solution of merely monotone or pseudo-monotone VIs are proposed.
All the proposed approaches though require at least first-order (gradient)
feedback. The challenge with generalizing these approaches to zeroth-order
(bandit) feedback is that they require coordination among the players and are
not suitable in the bandit feedback setting because a player cannot sample her
objective functions at different points while ensuring actions of other
players remains fixed. For example, the extra-gradient algorithm in [12]
remedies the cyclic behavior of trajectories in monotone game. However, it is
only applicable to exact first-order feedback oracle. On the other hand,
standard Tikhonov regularization requires a double iterative procedure, where
the players would solve a regularized VI corresponding to a regularized game
mapping in each inner iteration. Here, solving the inner VI itself would
require an iterative algorithm and it is not clear how the players should
coordinate stopping time for the inner algorithm in addition to setting the
regularization parameter.

Motivated by generalizing the algorithms for learning Nash equilibria, in our
recent work [19] we proposed an approach to learn Nash equilibria in merely
monotone games. Our approach was inspired by the single timescale Tikhonov
regularization algorithm for solving stochastic variational inequalities [11].
In contrast to the above work that assumed access to noisy gradients, we
considered the single-point estimation of gradients using the bandit feedback.
The main contribution was to appropriately control the bias and variance
introduced in the single-point estimation along with the Tikhonov
regularization parameter to ensure convergence. Our proposed algorithm was not
applicable to online learning because the played actions were not bound to lie
in the feasible set. Our current paper fills this gap by providing a
convergent algorithm while ensuring the query points do lie in the feasible
set.

Our contributions in this paper are as follows. We develop, to our best
knowledge, the first bandit approach for online learning Nash equilibria in
convex games with merely monotone game mappings. In doing so, we propose a
novel single time-scale double-regularization - this double regularization
corresponds to both the Tikhonov regularization and at the same time the
projection of the actions onto a shrunk feasible set. By properly managing the
interplay between the choice of the regularization sequence, the radius of
shrinkage for the feasible set and the stepsize we ensure that bias and
variance of the resulting stochastic approximation vanish with appropriate
rates. The choice of parameters are stated in Assumption 6, whereas the main
convergence result is stated in our Theorem 2. In terms of the proof
techniques, there are few main novelties leading to this theorem: first,
showing that the doubly regularized Tikhonov sequence and a single-time scale
approach to solve it, remain bounded and converges to the least-norm solution
of the variational inequality (see Theorem 3 and Proposition 1, respectively);
second, showing that iterates of our algorithm remain bounded (see Lemma 5).
These results enable us to use well-established results on boundedness and
convergence of stochastic processes to our setup. In summary, we prove
convergence in probability of the sequence of actions to a Nash equilibria in
monotone games.

Notations. The set $\\{1,\ldots,N\\}$ is denoted by $[N]$. Boldface is used to
distinguish between vectors in a multi-dimensional space and scalars. Given
$N$ vectors $\boldsymbol{x}^{i}\in\mathbb{R}^{d}$, $i\in[N]$,
$(\boldsymbol{x}^{i})_{i=1}^{N}:=({\boldsymbol{x}^{1}}^{\top},\ldots,{\boldsymbol{x}^{N}}^{\top})^{\top}\in\mathbb{R}^{Nd}$;
$\boldsymbol{x}^{-i}:=({\boldsymbol{x}^{1}},\ldots,{\boldsymbol{x}^{i-1}},{\boldsymbol{x}^{i+1}},\ldots,{\boldsymbol{x}^{N}})\in\mathbb{R}^{(N-1)d}$.
$\mathbb{R}^{d}_{+}$ and $\mathbb{Z}_{+}$ denote, respectively, vectors from
$\mathbb{R}^{d}$ with non-negative coordinates and non-negative whole numbers.
The standard inner product on $\mathbb{R}^{d}$ is denoted by $(\cdot,\cdot)$:
$\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$, with associated norm
$\|\boldsymbol{x}\|:=\sqrt{(\boldsymbol{x},\boldsymbol{x})}$. Given some
matrix $A\in\mathbb{R}^{d\times d}$, $A\succeq(\succ)0$, if and only if
$\boldsymbol{x}^{\top}A\boldsymbol{x}\geq(>)0$ for all $\boldsymbol{x}\neq 0$.
We use the big-$O$ notation, that is, the function
$f(x):\mathbb{R}\to\mathbb{R}$ is $O(\mathbf{g}(x))$ as $x\to a$, $f(x)$ =
$O(g(x))$ as $x\to a$, if $\lim_{x\to a}\frac{|f(x)|}{|g(x)|}\leq K$ for some
positive constant $K$. And with a slight abuse of notation, we write $f(x)\leq
O(g(x))$ as we estimate certain bounds. We say that a function
$f(\boldsymbol{x})$ grows not faster than a function $g(\boldsymbol{x})$ as
$\boldsymbol{x}\to\infty$, if there exists a positive constant $Q$ such that
$f(\boldsymbol{x})\leq g(\boldsymbol{x})$ $\forall\boldsymbol{x}$ with
$\|\boldsymbol{x}\|\geq Q$. For $x\in\mathbb{R}^{n}$ and convex closed set
$C\subset\mathbb{R}^{n}$, $\mbox{Proj}_{C}x$ denotes the projection of $x$
onto $C$.

###### Definition 1

A mapping $\boldsymbol{M}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is _monotone_ over
$X\subseteq\mathbb{R}^{d}$, if
$(\boldsymbol{M}(\boldsymbol{x})-\boldsymbol{M}(\boldsymbol{y}),\boldsymbol{x}-\boldsymbol{y})\geq
0$ for every $\boldsymbol{x},\boldsymbol{y}\in X$.

##  II Problem Formulation

Consider a game $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$ with $N$ players, the sets
of players’ actions $A_{i}\subseteq\mathbb{R}^{d}$, $i\in[N]$, and the cost
(objective) functions $J_{i}:\boldsymbol{A}\to\mathbb{R}$, where
$\boldsymbol{A}=A_{1}\times\ldots\times A_{N}$ denotes the set of joint
actions. We restrict the class of games as follows.

###### Assumption 1

The game under consideration is _convex_. Namely, for all $i\in[N]$ the set
$A_{i}$ is convex and closed, the cost function
$J_{i}(\boldsymbol{a}^{i},\boldsymbol{a}^{-i})$ is defined on
$\mathbb{R}^{Nd}$, continuously differentiable in $\boldsymbol{a}$ and convex
in $\boldsymbol{a}^{i}$ for fixed $\boldsymbol{a}^{-i}$.

###### Assumption 2

The mapping $\boldsymbol{M}:\mathbb{R}^{Nd}\to\mathbb{R}^{Nd}$, referred to as
the _game mapping_ , defined by

$\displaystyle\boldsymbol{M}(\boldsymbol{a})$
$\displaystyle=(\nabla_{\boldsymbol{a}^{i}}J_{i}(\boldsymbol{a}^{i},\boldsymbol{a}^{-i}))_{i=1}^{N}=(\boldsymbol{M}_{1}(\boldsymbol{a}),\ldots,\boldsymbol{M}_{N}(\boldsymbol{a}))^{\top},$
$\displaystyle\mbox{where
}\boldsymbol{M}_{i}(\boldsymbol{a})=(M_{i,1}(\boldsymbol{a}),\ldots,M_{i,d}(\boldsymbol{a}))^{\top},\mbox{
and}$ $\displaystyle M_{i,k}(\boldsymbol{a})$ $\displaystyle=\frac{\partial
J_{i}(\boldsymbol{a})}{\partial
a^{i}_{k}},\quad\boldsymbol{a}\in\boldsymbol{A},\quad i\in[N],\quad k\in[d],$
(1)

is _monotone on $\boldsymbol{A}$_ (see Definition 1).

We consider a _Nash equilibrium_ in game $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$
as a stable solution outcome because it represents a joint action from which
no player has any incentive to unilaterally deviate.

###### Definition 2

A point $\boldsymbol{a}^{*}\in\boldsymbol{A}$ is called a _Nash equilibrium_
if for any $i\in[N]$ and $\boldsymbol{a}^{i}\in A_{i}$

$J_{i}(\boldsymbol{a}^{i*},\boldsymbol{a}^{-i*})\leq
J_{i}(\boldsymbol{a}^{i},\boldsymbol{a}^{-i*}).$

Our goal is to learn such a stable action in a game through designing a
payoff-based algorithm. To do so, we first connect existence of Nash
equilibria for $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$ with solution set of a
corresponding variational inequality problem.

###### Definition 3

Consider a mapping $\boldsymbol{T}(\cdot)$: $\mathbb{R}^{d}\to\mathbb{R}^{d}$
and a set $Y\subseteq\mathbb{R}^{d}$. The _solution set
$SOL(Y,\boldsymbol{T})$ to the variational inequality problem_
$VI(Y,\boldsymbol{T})$ is the set of vectors $\boldsymbol{y}^{*}\in Y$ such
that
$(\boldsymbol{T}(\boldsymbol{y}^{*}),\boldsymbol{y}-\boldsymbol{y}^{*})\geq
0$, $\forall\boldsymbol{y}\in Y$.

###### Theorem 1

([16, Proposition 1.4.2]) Given a game $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$
with game mapping $\boldsymbol{M}$, suppose that the action sets $\\{A_{i}\\}$
are closed and convex, the cost functions $\\{J_{i}\\}$ are continuously
differentiable in $\boldsymbol{a}$ and convex in $\boldsymbol{a}^{i}$ for
every fixed $\boldsymbol{a}^{-i}$ on the interior of $\boldsymbol{A}$. Then,
some vector $\boldsymbol{a}^{*}\in\boldsymbol{A}$ is a Nash equilibrium in
$\Gamma$, if and only if $\boldsymbol{a}^{*}\in
SOL(\boldsymbol{A},\boldsymbol{M})$.

It follows that under Assumptions 1 and 2 for a game with mapping
$\boldsymbol{M}$, any solution of $VI(\boldsymbol{A},\boldsymbol{M})$ is also
a Nash equilibrium in such games and vice versa. While
$\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$ under Assumptions 1 and 2 might admit a
Nash equilibrium, these two assumptions alone do not guarantee uniqueness of a
Nash equilibrium. More restrictive assumptions for uniqueness are needed, for
example, strong monotonicity of the game mapping or compactness of the action
sets [16]. Here, we do not restrict our attention to such cases. However, to
have a meaningful discussion, we do assume existence of at least one Nash
equilibrium in the game.

###### Assumption 3

The set $SOL(\boldsymbol{A},\boldsymbol{M})$ is not empty.

###### Corollary 1

Let $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$ be a game with game mapping
$\boldsymbol{M}$ for which Assumptions 1, 2, and 3 hold. Then, there exists at
least one Nash equilibrium in $\Gamma$. Moreover, any Nash equilibrium in
$\Gamma$ belongs to the set $SOL(\boldsymbol{A},\boldsymbol{M})$.

The following additional assumptions are needed for convergence of the
proposed payoff-based algorithm to a Nash equilibrium (see proofs of Lemma 5
and Theorem 2).

###### Assumption 4

Each element $\boldsymbol{M}_{i}$ of the game mapping
$\boldsymbol{M}:\mathbb{R}^{Nd}\to\mathbb{R}^{Nd}$, defined in Assumption (2)
is Lipschitz continuous on $\mathbb{R}^{d}$ with a Lipschitz constant $L_{i}$.

###### Assumption 5

Each cost function $J_{i}(\boldsymbol{a})$, $i\in[N]$ grows at most
polynomially in $\boldsymbol{a}$ as $\|\boldsymbol{a}\|\to\infty$. Moreover,
in the case of unbounded joint action set $\boldsymbol{A}$, each continuous
cost function $J_{i}(\boldsymbol{a})$, $i\in[N]$ grows at most linearly in
$\boldsymbol{a}$ as $\|\boldsymbol{a}\|\to\infty$.

###### Remark 1

Note that if the set $\boldsymbol{A}$ is unbounded, Assumption 5 is equivalent
to each cost function $J_{i}(\boldsymbol{a})$, $i\in[N]$, being Lipschitz
continuous on $\mathbb{R}^{Nd}$ with some constant $l_{i}$ Thus, in both
bounded and unbounded $\boldsymbol{A}$, we denote $l=\max_{i\in[N]}l_{i}$ as
the uniform upper bound of the mapping $\boldsymbol{M}$ over $\boldsymbol{A}$.

For the development and analysis of our algorithms, we use the following well-
established and easy to verify result.

###### Lemma 1

Consider a mapping $\boldsymbol{T}(\cdot):\mathbb{R}^{d}\to\mathbb{R}^{d}$ and
a convex closed set $Y\subseteq\mathbb{R}^{d}$. Given $\theta>0$,

$\displaystyle\boldsymbol{y}^{*}\in\text{SOL}(Y,\boldsymbol{T})\iff\boldsymbol{y}^{*}=\mbox{Proj}_{Y}(\boldsymbol{y}^{*}-\theta\boldsymbol{T}(\boldsymbol{y}^{*})).$
(2)

##  III Payoff-Based Algorithm

Given online payoff-based information, also referred to as online bandit or
zeroth-order oracle information, each agent has access to its current action,
referred to as its state and denoted by
$\mathbf{x}^{i}(t)=(x^{i}_{1}(t),\ldots,x^{i}_{d}(t))^{\top}\in\mathbb{R}^{d}$,
and plays the action
$\boldsymbol{a}^{i}(t)=\mbox{Proj}_{A_{i}}(\mathbf{x}^{i}(t))$ at iteration
$t$. After that the cost value $\hat{J}_{i}(t)$ at the joint action
$\boldsymbol{a}(t)=(\boldsymbol{a}^{1}(t),\ldots,\boldsymbol{a}^{N}(t))\in\boldsymbol{A}$,
$\hat{J}_{i}(t)=J_{i}(\boldsymbol{a}(t))$ is revealed to each agent $i$. Given
these pieces of information, in the proposed algorithm each agent $i$ “mixes”
its next state $\mathbf{x}^{i}(t+1)$. Namely, it chooses $\mathbf{x}^{i}(t+1)$
randomly according to the multidimensional normal distribution
$\mathcal{N}(\boldsymbol{\mu}^{i}(t+1)=(\mu^{i}_{1}(t+1),\ldots,\mu^{i}_{d}(t+1))^{\top},\sigma(t+1))$
with the following density function:

$\displaystyle p_{i}$
$\displaystyle(\boldsymbol{x}^{i};\boldsymbol{\mu}^{i}(t+1),\sigma_{t+1})=p_{i}(x^{i}_{1},\ldots,x^{i}_{d};\boldsymbol{\mu}^{i}(t+1),\sigma_{t+1})$
$\displaystyle=\frac{1}{(\sqrt{2\pi}\sigma_{t+1})^{d}}\exp\left\\{-\sum_{k=1}^{d}\frac{(x^{i}_{k}-\mu^{i}_{k}(t+1))^{2}}{2\sigma^{2}_{t+1}}\right\\}.$

The initial value of the means $\boldsymbol{\mu}^{i}(0)$, $i\in[N]$, can be
set to any finite value. The successive means are updated as follows:

$\displaystyle\boldsymbol{\mu}^{i}(t+1)=\mbox{Proj}_{(1-r_{t})A_{i}}\big{[}\boldsymbol{\mu}^{i}(t)$
(3)
$\displaystyle-\gamma_{t}\sigma^{2}_{t}\left({\hat{J}_{i}(t)}\frac{{\mathbf{x}^{i}(t)}-\boldsymbol{\mu}^{i}(t)}{\sigma^{2}_{t}}+\epsilon_{t}\boldsymbol{\mu}^{i}(t)\right)\big{]}.$
(4)

In the above, $(1-r_{t})A_{i}=\\{\boldsymbol{x}\in
A_{i}:\mbox{dist}(\boldsymbol{x},\partial A_{i})\geq r_{t}\\}$ and
$0<r_{t}<1$, is a time-dependent shrinkage parameter, $\gamma_{t}$ is the
stepsize parameter and $\epsilon_{t}>0$ is a regularization parameter. The
convergence of the algorithm is dependent on the interplay of these parameters
and the variance term $\sigma_{t}>0$.

The difference between the proposed approach and that of [18] is due to the
additional term $\epsilon_{t}$ in (3). In the absence of this term the
algorithm would converge only if the game mapping is _strictly_ monotone (see
[18, Theorem 2] and counterexamples in [7, 13]). Moreover, in distinction from
[19], in the bandit online feedback considered here, players can only evaluate
their costs over their feasible action set $A_{i}$ and not over the whole
$\mathbb{R}^{Nd}$, necessitating the additional projection term
$\boldsymbol{a}^{i}(t)=\mbox{Proj}_{A_{i}}(\mathbf{x}^{i}(t))$ and the
shrinkage radius $r_{t}$. As such, the previous convergence analysis does not
apply.

Before stating the convergence result, let us provide insight into the
procedure defined by Equation (3) by deriving an analogy to a regularized
stochastic gradient algorithm.

Let
$p(\boldsymbol{x};\boldsymbol{\mu},\sigma)=\prod_{i=1}^{N}p_{i}(x^{i}_{1},\ldots,x^{i}_{d};\boldsymbol{\mu}^{i},\sigma)$
denote the density function of the joint distribution of agents’ states. Given
$\sigma>0$, for any $i\in[N]$ define
$\tilde{J}_{i}:\mathbb{R}^{Nd}\rightarrow\mathbb{R}$ as

$\displaystyle\tilde{J}_{i}$
$\displaystyle(\boldsymbol{\mu}^{1},\ldots,\boldsymbol{\mu}^{N},\sigma)=\int_{\mathbb{R}^{Nd}}J_{i}(\boldsymbol{x})p(\boldsymbol{x};\boldsymbol{\mu},\sigma)d\boldsymbol{x}.$
(5)

Thus, $\tilde{J}_{i}$, $i\in[N]$ is the $i$th player’s cost function in mixed
strategies. Let
$\boldsymbol{\mu}(t)=(\boldsymbol{\mu}^{1}(t),\ldots,\boldsymbol{\mu}^{N}(t))$
and for $i\in[N]$, define
$\tilde{\boldsymbol{M}}_{i}(\cdot)=(\tilde{M}_{i,1}(\cdot),\ldots,\tilde{M}_{i,d}(\cdot))^{\top}$
as the $d$-dimensional mapping with the following elements:

$\displaystyle\tilde{M}_{i,k}(\boldsymbol{\mu},\sigma)=\frac{\partial{\tilde{J}_{i}(\boldsymbol{\mu},\sigma)}}{\partial\mu^{i}_{k}},\mbox{
for $k\in[d]$}.$ (6)

Our first lemma below shows that the second term inside the projection in (3)
is a sample of the gradient of agent $i$’s cost function in mixed strategies.

###### Lemma 2

Under Assumptions 1 and 5

$\displaystyle\tilde{M}_{i,k}(\boldsymbol{\mu}(t),\sigma_{t})=\mathrm{E}\\{J_{i}(\mathbf{x}^{1}(t),\ldots,\mathbf{x}^{N}(t))\frac{x^{i}_{k}(t)-\mu^{i}_{k}(t)}{\sigma^{2}_{t}}|$
$\displaystyle\qquad\qquad
x^{i}_{k}(t)\sim\mathcal{N}(\mu_{k}^{i}(t),\sigma_{t}),i\in[N],k\in[d]\\}.$
(7)

Moreover, $\tilde{M}_{i,k}(\boldsymbol{\mu},\sigma)$ is bounded for any
$\boldsymbol{\mu}\in\boldsymbol{A}$.

The proof of this Lemma is very similar to that of [19] and is provided in
Appendix -B.

The lemma above implies that had we used the term
$J_{i}(\mathbf{x}^{1}(t),\ldots,\mathbf{x}^{N}(t))\frac{x^{i}_{k}(t)-\mu^{i}_{k}(t)}{\sigma^{2}_{t}}$
in (3) we could perform a one-point estimation of the gradient of the cost
functions in mixed strategies. In the bandit setting considered here, however,
we use the term
$J_{i}(\boldsymbol{a}^{1}(t),\ldots,\boldsymbol{a}^{N}(t))\frac{x^{i}_{k}(t)-\mu^{i}_{k}(t)}{\sigma^{2}_{t}}$.
Despite this difference, in the analysis (see (68), (73)) we show that the
difference between these two terms converges to zero due to the shrinkage
radius selection. Thus, algorithm (3) can be interpreted as a doubly
regularized (due to $r_{r}$ and $\epsilon_{t}$) stochastic projection
algorithms.

Following the above interpretation, our main result is Theorem 2 below where
we show that by appropriately choosing the algorithm parameters, we can bound
the bias and variance terms of the stochastic projection and consequently
establish convergence of the iterates $\boldsymbol{\mu}(t)$ in (3) to a Nash
equilibrium.

###### Assumption 6

Let $\beta_{t}=\gamma_{t}\sigma^{2}_{t}$ and choose
$\gamma_{t}=\frac{1}{t^{a}}$, $\sigma_{t}=\frac{1}{t^{b}}$,
$\epsilon_{t}=\frac{1}{t^{c}}$ and $r_{t}=\frac{1}{t^{d}}$ $a,b,c,d>0$
respectively, such that

i) $\sum_{t=0}^{\infty}\beta_{t}=\infty$,
$\sum_{t=0}^{\infty}\beta_{t}\epsilon_{t}=\infty,$

ii)
$\sum_{t=0}^{\infty}\frac{(\epsilon_{t}-\epsilon_{t-1})^{2}}{\beta_{t}\epsilon^{3}_{t}}+\frac{(r_{t}-r_{t-1})^{2}}{\beta_{t}\epsilon^{6}_{t}}<\infty,$

iii) $\sum_{t=0}^{\infty}\gamma^{2}_{t}<\infty,$
$\sum_{t=0}^{\infty}\beta_{t}\sigma_{t}<\infty$,

iv) $\lim_{t\to\infty}\frac{r_{t}}{\sigma_{t}}=\infty$,
$\lim_{t\to\infty}\frac{r_{t}}{\epsilon_{t}}=0$.

As an example for existence of parameters to satisfy Assumption 6, let
$a_{1}=\frac{5}{9}$, $a_{2}=\frac{5}{27}$, $a_{3}=\frac{1}{54}$,
$a_{4}=\frac{1}{6}$.

###### Theorem 2

Let the players in game $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$ choose the states
$\\{\mathbf{x}^{i}(t)\\}$ at time $t$ according to the normal distribution
$\mathcal{N}(\boldsymbol{\mu}^{i}(t),\sigma_{t})$, where the mean
$\boldsymbol{\mu}^{i}(0)$ is arbitrary and $\boldsymbol{\mu}^{i}(t)$ is
updated as in (3). Under Assumptions 1-6, as $t\to\infty$, the mean vector
$\boldsymbol{\mu}(t)$ converges almost surely to a Nash equilibrium
$\boldsymbol{\mu}^{*}=\boldsymbol{a}^{*}$ of the game $\Gamma$ and the joint
action $\boldsymbol{a}(t)$ converges in probability to $\boldsymbol{a}^{*}$.

##  IV Analysis of the Algorithm

To prove Theorem 2 we need to first establish boundedness of the iterates
$\boldsymbol{\mu}(t)$ for the cases in which the action space is unbounded.
Having established the boundedness, we can show that the limit of the iterates
$\boldsymbol{\mu}(t)$ exists and is the minimum norm Nash equilibrium of the
problem. This convergence is proven using existing results on convergence of a
sequence of random variables (Lemma 10 on page 49 in [17]). For ease of
reference, we provide the supporting statements used for boundedness [15,
Theorem 2.5.2] and for convergence [17, Lemma 10 on page 49] of the iterates
in Appendix -A.

###  IV-A Characterizing the terms in the algorithm

We first show that algorithm (3) can be interpreted within the framework of
well-studied Robbins-Monro stochastic approximations procedures [3], where the
iterates are updated according to stochastic gradient descent. In our case,
bias of the game mapping arises due to each player’s one-point estimation of
its gradient. However, in contrast to a stochastic approximation procedure,
the game mapping is in general not gradient of a single function (as its
derivative is not symmetric) unless the game is potential. Furthermore, there
are additional terms in the algorithm iterates due to the projection of the
query points onto the shrunk feasible set and the regularization required to
handle the non-strictly monotone game mapping. Let us specify all these terms
below.

Using the notation
$\boldsymbol{M}_{i}(\cdot)=(M_{i,1}(\cdot),\ldots,M_{i,d}(\cdot))$, we can
rewrite the algorithm step in (3) in the following form:

$\displaystyle\boldsymbol{\mu}^{i}(t+1)=\mbox{Proj}_{(1-r_{t})A_{i}}[\boldsymbol{\mu}^{i}(t)-\gamma_{t}\sigma^{2}_{t}$
(8)
$\displaystyle\times\big{(}\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))+\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t),\sigma_{t})+\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})$
(9)
$\displaystyle\qquad\qquad+\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})+\epsilon_{t})\boldsymbol{\mu}^{i}(t)\big{)}],$
(10)

for all $i\in[N]$, where $\boldsymbol{Q}_{i}$, $\boldsymbol{R}_{i}$,
$\boldsymbol{P}_{i}$ are

$\displaystyle\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t),\sigma_{t})=\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu}(t),\sigma_{t})-\boldsymbol{M}_{i}(\boldsymbol{\mu}(t)),$
$\displaystyle\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})=\boldsymbol{F}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})-\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu}(t),\sigma_{t}),$
$\displaystyle\boldsymbol{F}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})=J_{i}(\mathbf{x}(t))\frac{\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}(t)}{\sigma^{2}_{t}},$
$\displaystyle\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})=\frac{\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}(t)}{\sigma^{2}_{t}}(J_{i}(\boldsymbol{a}(t))-J_{i}(\mathbf{x}(t))).$

Above,
$\boldsymbol{M}(\boldsymbol{\mu}(t))=(\boldsymbol{M}_{1}(\boldsymbol{\mu}(t)),\ldots,\boldsymbol{M}_{N}(\boldsymbol{\mu}(t)))$
corresponds to the gradient term in stochastic approximation procedures. The
mapping $\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu}(t))$ evaluated at
$\boldsymbol{\mu}(t)$ is equivalent to the game mapping in mixed strategies
[19]. That is,

$\displaystyle\tilde{\boldsymbol{M}}_{i}($
$\displaystyle\boldsymbol{\mu}(t))=\int_{\mathbb{R}^{Nd}}{\boldsymbol{M}_{i}}(\boldsymbol{x})p(\boldsymbol{x};\boldsymbol{\mu}(t),\sigma_{t})d\boldsymbol{x}.$
(11)

Thus, the
$\boldsymbol{Q}(\boldsymbol{\mu}(t),\sigma_{t})=(\boldsymbol{Q}_{1}(\boldsymbol{\mu}(t),\sigma_{t}),\ldots,\boldsymbol{Q}_{N}(\boldsymbol{\mu}(t),\sigma_{t}))$
can be interpreted as disturbance of the gradient. Furthermore, according to
(2), we have222The notation $\mathrm{E}_{\mathbf{x}(t)}\\{\cdot\\}$ further is
used to emphasize that the expectation is taken in respect to $\mathbf{x}(t)$
which has the normal distribution with the mean $\boldsymbol{\mu}(t)$ and the
covariance matrix $\sigma(t)I$, where $I$ is the identity matrix.

$\displaystyle\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})=\boldsymbol{F}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})$
(12)
$\displaystyle-\mathrm{E}_{\mathbf{x}(t)}\\{\boldsymbol{F}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})\\},\;\;i\in[N].$

Thus, $\boldsymbol{R}$ below is a martingale difference

$\displaystyle\boldsymbol{R}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})=($
$\displaystyle\boldsymbol{R}_{1}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t}),\ldots,$
$\displaystyle\boldsymbol{R}_{N}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})).$

Finally, due to the projection of query points, the term

$\displaystyle\boldsymbol{P}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})=($
$\displaystyle\boldsymbol{P}_{1}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t}),\ldots,$
$\displaystyle\boldsymbol{P}_{N}(\mathbf{x}(t),\boldsymbol{\mu}(t),\sigma_{t})),$

is the vector of the difference between the gradient estimation based on the
state $\mathbf{x}(t)\in\mathbb{R}^{Nd}$ and the played action
$\boldsymbol{a}(t)\in\boldsymbol{A}$.

Our goal is to bound each of the terms above to ensure convergence and
boundedness of the iterates. However, first, we need to account for the
regularization terms $\epsilon_{t}$, $r_{t}$.

###  IV-B Analyzing the modified Tikhonov sequence

In contrast to stochastic approximation algorithms and the proof in [18], we
have an addition term $\epsilon_{t}\boldsymbol{\mu}(t)$ to be able to address
merely monotone game mappings. As such, to bound $\boldsymbol{\mu}(t)$ we also
relate the variations of the sequence $\boldsymbol{\mu}(t)$ to those of the
_modified Tikhonov sequence_ defined below. Let
$\boldsymbol{y}(t)=(\boldsymbol{y}^{1}(t),\ldots,\boldsymbol{y}^{N}(t))$
denote the solution of the variational inequality
$VI((1-r_{t})\boldsymbol{A},\boldsymbol{M}(\boldsymbol{y})+\epsilon_{t}\boldsymbol{y})$,
namely

$\displaystyle\boldsymbol{y}(t)\in
SOL((1-r_{t})\boldsymbol{A},\boldsymbol{M}(\boldsymbol{y})+\epsilon_{t}\boldsymbol{y}).$
(13)

The Tikhonov sequence corresponds to the solution of the variational
inequality above with the $r_{t}=0$. Thus, the sequence
$\\{\boldsymbol{y}(t)\\}$ can be considered the _modified_ Tikhonov sequence.
Similar to the Tikhonov sequence, $\boldsymbol{y}(t)$ enjoys the following
important property.

###### Theorem 3

Under Assumptions 2, 3, and 4, $\boldsymbol{y}(t)$ defined in (13) exists and
is unique for each $t$. Moreover, for $\epsilon_{t}\downarrow 0$ and $r_{t}\to
0$ and given $\lim_{t\to\infty}\frac{r_{t}}{\epsilon_{t}}=0$,
$\boldsymbol{y}(t)$ is uniformly bounded and converges to the least norm
solution of $VI(\boldsymbol{A},\boldsymbol{M})$.

The significance of the above theorem is that if we can establish
$\boldsymbol{\mu}(t)$ converges to $\boldsymbol{y}(t)$, then from 1 we can
establish convergence to a Nash equilibrium for the game. To prove Theorem 3,
first we establish some useful properties of projecting onto sets
$(1-r_{t})\boldsymbol{A}$, $t=1,2,\ldots$.

###### Lemma 3

For any $\mathbf{x}\in\mathbb{R}^{Nd}$ the following holds:

$\|\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}\mathbf{x}-\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}\mathbf{x}\|=O(|r_{t-1}-r_{t}|).$

Please see the Appendix -C for the proof of above Lemma.

###### Proof:

(of Theorem 3) Let $\boldsymbol{a}$ be the least norm solution of
$VI(\boldsymbol{A},\boldsymbol{M})$. Moreover, let $\boldsymbol{a}^{p}$ be the
projection of $\boldsymbol{a}$ onto the set $(1-r_{t})\boldsymbol{A}$. Next,
let $\boldsymbol{y}(t)$ be the unique solution of the doubly regularized
inequality, namely $\boldsymbol{y}(t)\in
SOL((1-r_{t})\boldsymbol{A},\boldsymbol{M}+\epsilon_{t}I)$. Thus, we conclude
that

$(\boldsymbol{M}(\boldsymbol{a}),\boldsymbol{y}(t)-\boldsymbol{a})\geq 0,$
$(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t),\boldsymbol{a}^{p}-\boldsymbol{y}(t))\geq
0.$

Thus, taking into account monotonicity of $\boldsymbol{M}$, we obtain

$\displaystyle 0\leq$
$\displaystyle(\boldsymbol{M}(\boldsymbol{a}),\boldsymbol{y}(t)-\boldsymbol{a})$
$\displaystyle+(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t),\boldsymbol{a}-\boldsymbol{y}(t))$
$\displaystyle+(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t),\boldsymbol{a}^{p}-\boldsymbol{a})$
$\displaystyle=$
$\displaystyle-(\boldsymbol{M}(\boldsymbol{a})-\boldsymbol{M}(\boldsymbol{y}(t)),\boldsymbol{a}-\boldsymbol{y}(t))+\epsilon_{t}(\boldsymbol{y}(t),\boldsymbol{a}-\boldsymbol{y}(t))$
$\displaystyle+(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t),\boldsymbol{a}^{p}-\boldsymbol{a})$
$\displaystyle\leq$
$\displaystyle\epsilon_{t}(\boldsymbol{y}(t),\boldsymbol{a})-\epsilon_{t}\|\boldsymbol{y}(t)\|^{2}$
$\displaystyle+(\boldsymbol{M}(\boldsymbol{y}(t)),\boldsymbol{a}^{p}-\boldsymbol{a})+\epsilon_{t}(\boldsymbol{y}(t),\boldsymbol{a}^{p}-\boldsymbol{a}).$

Hence,

$\displaystyle\epsilon_{t}\|\boldsymbol{y}(t)\|^{2}$
$\displaystyle\leq\epsilon_{t}(\boldsymbol{y}(t),\boldsymbol{a})+(\boldsymbol{M}(\boldsymbol{y}(t)),\boldsymbol{a}^{p}-\boldsymbol{a})$
$\displaystyle\qquad+\epsilon_{t}(\boldsymbol{y}(t),\boldsymbol{a}^{p}-\boldsymbol{a})$
$\displaystyle\leq\epsilon_{t}\|\boldsymbol{y}(t)\|\|\boldsymbol{a}\|+l\|\boldsymbol{a}^{p}-\boldsymbol{a}\|+\epsilon_{t}\|\boldsymbol{y}(t)\|\|\boldsymbol{a}^{p}-\boldsymbol{a}\|$
$\displaystyle=\epsilon_{t}\|\boldsymbol{y}(t)\|\|\boldsymbol{a}\|+lO(r_{t})+\epsilon_{t}\|\boldsymbol{y}(t)\|O(r_{t}),$

where in the first inequality we used Remark 1 and in the second one we
applied Lemma 3. Hence,

$\|\boldsymbol{y}(t)\|^{2}\leq\|\boldsymbol{y}(t)\|\|\boldsymbol{a}\|+lO\left(\frac{r_{t}}{\epsilon_{t}}\right)+\|\boldsymbol{y}(t)\|O(r_{t}).$

By taking the upper limit $t\to\infty$ in the inequality above and due to the
settings for $\epsilon_{t}$ and $r_{t}$, we obtain

$\displaystyle\overline{\lim}_{t\to\infty}[\|\boldsymbol{y}(t)\|^{2}]\leq$
$\displaystyle\|\boldsymbol{a}\|\overline{\lim}_{t\to\infty}\|\boldsymbol{y}(t)\|$
$\displaystyle+\overline{\lim}_{t\to\infty}\|\boldsymbol{y}(t)\|\overline{\lim}_{t\to\infty}O(r_{t}).$

It implies that
$\overline{\lim}_{t\to\infty}\|\boldsymbol{y}(t)\|\leq\|\boldsymbol{a}\|$,
and, thus, the sequence $\|\boldsymbol{y}(t)\|$ is upper bounded. Moreover,
any accumulation point of this sequence is bounded above by the Euclidean norm
of the least-norm solution of $VI(\boldsymbol{A},\boldsymbol{M})$. Hence,
according to the fact that the function
$\mbox{Proj}_{(1-r)\boldsymbol{A}}(\mathbf{x})$ is continuous in both $r$ and
$\mathbf{x}$, $\boldsymbol{y}(t)$ converges to the least norm solution of
$VI(\boldsymbol{A},\boldsymbol{M})$. ∎

Since our goal now is to relate $\boldsymbol{\mu}(t)$ to $\boldsymbol{y}(t)$,
aligned with procedure (3), we will now design a one-time scale approach to
solving (13) as per (14) below:

$\displaystyle\mathbf{z}(t+1)=\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\mathbf{z}(t)-\beta_{t}(\boldsymbol{M}(\mathbf{z}(t))+\epsilon_{t}\mathbf{z}(t))],$
(14)

where $\beta_{t}$ is defined in Assumption 6. We show that the procedure above
is a one time-scale algorithm and similarly to $\boldsymbol{y}(t)$, it
converges to the least norm solution of $VI(\boldsymbol{A},\boldsymbol{M})$.

###### Proposition 1

Under conditions of Theorem 3 and given Assumption 6, the sequence
$\mathbf{z}(t)$ defined by (14) converges to the least norm solution of
$VI(\boldsymbol{A},\boldsymbol{M})$.

To prove the result above, we bound $\|\mathbf{z}(t+1)-\boldsymbol{y}(t)\|$ in
terms of the previous terms in the sequence, namely,
$\|\mathbf{z}(t)-\boldsymbol{y}(t-1)\|$ and show that [17, Lemma 10, page 49]
on convergence of a random sequence applies. To do so though, first we need to
bound the variations of $\boldsymbol{y}(t)$ as below.

###### Lemma 4

Under Assumptions 2, 4, and 6, the Tikhonov sequence $\boldsymbol{y}(t)$
defined in (13) satisfies

$\displaystyle\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}=$ $\displaystyle
O\left(\frac{(\epsilon_{t}-\epsilon_{t-1})^{2}}{\epsilon^{2}_{t}}+\frac{(r_{t}-r_{t-1})^{2}}{\epsilon^{5}_{t}}\right).$

Please see Appendix -D for the proof.

In summary, the results of this section enabled us to prove Proposition 1.
This proposition serves as the main new result in comparison to non-
regularized stochastic gradient procedures in order to show almost-sure
boundedness of $\|\boldsymbol{\mu}(t)\|$ and the convergence of the algorithm
to a Nash equilibrium. Please see Appendix -E for the proof.

###  IV-C Boundedness of the iterates

###### Lemma 5

Let Assumptions 2-6 hold in $\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$ and
$\boldsymbol{\mu}(t)$ be the vector updated in the run of the payoff-based
algorithm (3). Then, $\Pr\\{\sup_{t\geq
0}\|\boldsymbol{\mu}(t)\|<\infty\\}=1$.

###### Proof:

If the set $\boldsymbol{A}$ is compact, then, according to the updated in (8),
the norm of the vector $\boldsymbol{\mu}_{t}$ is bounded for all $t$. So, let
us consider the case of the unbounded $\boldsymbol{A}$.

Define $V(t,\boldsymbol{\mu})=\|\boldsymbol{\mu}-\mathbf{z}(t)\|^{2}$, where
$\mathbf{z}(t)$ is given in (13). Consider the generating operator of the
Markov process $\boldsymbol{\mu}(t)$

$\displaystyle
LV(t,\boldsymbol{\mu})=E[V(t+1,\boldsymbol{\mu}(t+1))\mid\boldsymbol{\mu}(t)=\boldsymbol{\mu}]-V(t,\boldsymbol{\mu}).$

We aim to show that $LV(t,\boldsymbol{\mu})$ satisfies the following decay

$\displaystyle
LV(t,\boldsymbol{\mu})\leq-\alpha(t+1)\psi(\boldsymbol{\mu})+\phi(t)(1+V(t,\boldsymbol{\mu})),$
(15)

where $\psi$, $\phi$ and $\alpha$ are terms arising from $\boldsymbol{Q}$,
$\boldsymbol{R}$ and $\boldsymbol{P}$ in (8). Our goal is to show that
$\psi\geq 0$ on $\mathbb{R}^{Nd}$, $\phi(t)>0$, $\forall t$,
$\sum_{t=0}^{\infty}\phi(t)<\infty$, $\alpha(t)>0$,
$\sum_{t=0}^{\infty}\alpha(t)=\infty$. This combined with the boundedness of
the iterates $\mathbf{z}(t)$ stated in Proposition 1 enable us to apply
Theorem 2.5.2 in [15] to conclude almost sure boundedness of
$\boldsymbol{\mu}(t)$.

From now on, for simplicity in notation, we omit the argument $\sigma_{t}$ in
the terms $\tilde{\boldsymbol{M}}$, $\boldsymbol{Q}$, and $\boldsymbol{R}$. In
certain derivations, for the same reason we omit the time parameter $t$ as
well.

Let us analyze each term $i=1,\dots,N$ in

$V(t+1,\boldsymbol{\mu}(t+1))=\sum_{i=1}^{N}\|\boldsymbol{\mu}^{i}(t+1)-\mathbf{z}^{i}(t+1)\|^{2}.$

From the procedures for the update of $\boldsymbol{\mu}(t)$ and
$\mathbf{z}(t)$ and the non-expansion property of the projection operator, we
obtain

$\displaystyle\|$
$\displaystyle\boldsymbol{\mu}^{i}(t+1)-\mathbf{z}^{i}(t+1)\|^{2}$ (16)
$\displaystyle\leq\|\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t)-\beta_{t}\big{[}\epsilon_{t}(\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$
(17)
$\displaystyle+(\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\mathbf{z}(t))+\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t))+\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))$
(18)
$\displaystyle\qquad\qquad+\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\big{]}\|^{2}$
(19) $\displaystyle=\|\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t)\|^{2}$ (20)
$\displaystyle\qquad-2\beta_{t}(\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\mathbf{z}(t)),\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$
(21)
$\displaystyle\qquad-2\beta_{t}\epsilon_{t}(\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t),\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$
(22)
$\displaystyle\qquad-2\beta_{t}(\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t))+\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)),\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$
(23)
$\displaystyle\qquad-2\beta_{t}(\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)),\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$
(24)
$\displaystyle\qquad+\beta^{2}_{t}\|\boldsymbol{G}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2},$
(25)

where, for ease of notation, we have defined

$\displaystyle\boldsymbol{G}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))=$
$\displaystyle\epsilon_{t}(\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$ (26)
$\displaystyle+\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\mathbf{z}(t))$
(27)
$\displaystyle+\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t))+\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))$
(28) $\displaystyle+\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)).$
(29)

Note that the terms in
$\|\boldsymbol{G}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}$ are given as

$\displaystyle\|\boldsymbol{G}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}={\epsilon^{2}(t)\|\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t)\|^{2}}$
(30)
$\displaystyle+\|\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\mathbf{z}(t))\|^{2}$
(31)
$\displaystyle+\|\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t))\|^{2}+\|\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}+\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}$
(32)
$\displaystyle+2(\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t)),\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))$
(33)
$\displaystyle+2(\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)),\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))$
(34)
$\displaystyle+2(\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t)),\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))$
(35)
$\displaystyle+2\epsilon_{t}(\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\mathbf{z}(t)),\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))$
(36)
$\displaystyle+2(\epsilon_{t}(\boldsymbol{\mu}^{i}(t)-\mathbf{z}^{i}(t))+\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\mathbf{z}(t)),$
(37)
$\displaystyle\qquad\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t))+\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))+\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))),$
(38)

Thus, accounting for the above, for (12), which implies
$\mathrm{E}\\{\|\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\||\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}=0$
for any $\boldsymbol{\mu}$, and for the Cauchy-Schwarz inequality, we get from
(16)

$\displaystyle\mathrm{E}\\{\|\boldsymbol{\mu}^{i}(t+1)-\mathbf{z}^{i}(t+1)\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(39)
$\displaystyle\leq(1-2\beta_{t}\epsilon_{t})\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|^{2}$
(40)
$\displaystyle\quad-2\beta_{t}(\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t)),\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t))$
(41)
$\displaystyle\quad-2\beta_{t}(\boldsymbol{Q}_{i}(\boldsymbol{\mu}),\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t))$
(42)
$\displaystyle\quad+2\beta_{t}\mathrm{E}\\{\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\||\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|$
(43)
$\displaystyle\quad+\beta^{2}_{t}\mathrm{E}\\{\|\boldsymbol{G}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(44)
$\displaystyle\leq(1-2\beta_{t}\epsilon_{t})\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|^{2}$
(45)
$\displaystyle\quad-2\beta_{t}(\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t)),\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t))$
(46)
$\displaystyle\quad+2\beta_{t}\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t))\|$
(47)
$\displaystyle\quad+2\beta_{t}\mathrm{E}\\{\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\||\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|$
(48)
$\displaystyle{+\beta^{2}_{t}\epsilon^{2}(t)\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|^{2}}$
(49)
$\displaystyle+\beta^{2}_{t}[\|\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t))\|^{2}$
(50) $\displaystyle\quad+\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|^{2}$ (51)
$\displaystyle\quad+\mathrm{E}\\{\|\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}+\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(52)
$\displaystyle\quad+2\mathrm{E}\\{(\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)),\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(53)
$\displaystyle\quad+2\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|\mathrm{E}\\{\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))\||\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(54)
$\displaystyle\quad+2(\epsilon_{t}\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|+\|\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t))\|)$
(55)
$\displaystyle\quad\times(\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|+\mathrm{E}\\{\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\||\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\})].$
(56)

We proceed estimating the terms in the inequality above. Due to Assumption 4,
we conclude that

$\displaystyle\|\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t))\|^{2}\leq
L^{2}_{i}\|\boldsymbol{\mu}-\mathbf{z}(t)\|^{2}=O(V(t,\boldsymbol{\mu}))$ (57)
$\displaystyle(\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t)),\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t))$
(58)
$\displaystyle\leq\|\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t))\|\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|$
(59) $\displaystyle\leq
L_{i}\|\boldsymbol{\mu}-\mathbf{z}(t)\|\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|=O(V(t,\boldsymbol{\mu})).$
(60)

Let us analyze the terms containing the disturbance of gradient, namely
$\boldsymbol{Q}_{i}$, in Equation (30). Since
$\boldsymbol{Q}_{i}(\boldsymbol{\mu}(t))=\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu}(t))-\boldsymbol{M}_{i}(\boldsymbol{\mu}(t))$,
due to Assumption 2 and Equation (11), we obtain

$\displaystyle\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|$
$\displaystyle=\|\int_{\mathbb{R}^{Nd}}[\boldsymbol{M}_{i}(\boldsymbol{x})-\boldsymbol{M}_{i}(\boldsymbol{\mu})]p(\boldsymbol{x};\boldsymbol{\mu},\sigma_{t})d\boldsymbol{x}\|$
(61)
$\displaystyle\leq\int_{\mathbb{R}^{Nd}}\|\boldsymbol{M}_{i}(\boldsymbol{x})-\boldsymbol{M}_{i}(\boldsymbol{\mu})\|p(\boldsymbol{x};\boldsymbol{\mu},\sigma_{t})d\boldsymbol{x}$
(62)
$\displaystyle\leq\int_{\mathbb{R}^{Nd}}L_{i}\|\boldsymbol{x}-\boldsymbol{\mu}\|p(\boldsymbol{x};\boldsymbol{\mu},\sigma_{t})d\boldsymbol{x}$
(63)
$\displaystyle\leq\int_{\mathbb{R}^{Nd}}L_{i}\left(\sum_{i=1}^{N}\sum_{k=1}^{d}|x^{i}_{k}-\mu^{i}_{k}|\right)p(\boldsymbol{x};\boldsymbol{\mu},\sigma_{t})d\boldsymbol{x}$
$\displaystyle=O(\sigma_{t}),$ (64)

where the last equality is due to the fact that the first central absolute
moment of a random variable with a normal distribution
$\mathcal{N}(\mu,\sigma)$ is $O(\sigma)$. The estimation above imply, in
particular, that for any $\boldsymbol{\mu}\in\boldsymbol{A}$

$\displaystyle\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t)\|=O(\sigma_{t})(1+V(t,\boldsymbol{\mu}))$
(65)
$\displaystyle\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|\|M_{i}(\boldsymbol{\mu})-M_{i}(\mathbf{z}(t))\|\leq
L_{i}\|\boldsymbol{Q}_{i}(\boldsymbol{\mu})\|\|\boldsymbol{\mu}-\mathbf{z}(t)\|$
$\displaystyle=O(\sigma_{t})(1+V(t,\boldsymbol{\mu})).$ (66)

We bound the martingale term
$\|\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}$.

$\displaystyle\mathrm{E}\\{\|\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
$\displaystyle\leq\sum_{k=1}^{d}\int_{\mathbb{R}^{Nd}}{J_{i}}^{2}(\boldsymbol{x})\frac{(x^{i}_{k}-\mu^{i}_{k})^{2}}{\sigma^{4}(t)}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\leq\frac{f_{i}(\boldsymbol{\mu},\sigma_{t})}{\sigma^{4}(t)}=\frac{O(1+V(t,\boldsymbol{\mu}))}{\sigma^{4}(t)},$
(67)

where the first inequality is due to the fact that
$\mathrm{E}(\xi-\mathrm{E}\xi)^{2}\leq\mathrm{E}\xi^{2}$ and taking into
account (12), the second inequality is due to Assumption 5, with
$f_{i}(\boldsymbol{\mu},\sigma_{t})$ being a quadratic function of
$\boldsymbol{\mu}$ and $\sigma_{t}$, $i\in[N]$ (see Appendix -F for more
details).

We proceed estimating the terms containing
$P_{t}(\mathbf{x}_{t},\boldsymbol{\mu})$. For any
$\boldsymbol{\mu}\in\boldsymbol{A}$ we have

$\displaystyle\mathrm{E}\\{\|\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(68)
$\displaystyle=\mathrm{E}\frac{\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}|J_{i}(\mathbf{x}(t))-J_{i}(\boldsymbol{a}(t))|^{2}}{\sigma^{4}_{t}}$
(69)
$\displaystyle=\Pr\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\boldsymbol{A}\\}\mathrm{E}\frac{\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}|J_{i}(\mathbf{x}(t))-J_{i}(\boldsymbol{a}(t))|^{2}}{\sigma^{4}_{t}}$
(70)
$\displaystyle\leq\Pr\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\boldsymbol{A}\\}\mathrm{E}l^{2}\frac{\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}\|\mathbf{x}(t)-\boldsymbol{\mu}\|^{2}}{\sigma^{4}_{t}}$
(71)
$\displaystyle=k_{0}\Pr\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\boldsymbol{A}\\},\,\mbox{for
some $k_{0}>0$},$ (72)

where the inequality is due to Assumption (5) implying
$|J_{i}(\mathbf{x}(t))-J_{i}(\boldsymbol{a}(t))|^{2}\leq
l^{2}\|\mathbf{x}(t)-\boldsymbol{a}(t)\|^{2}$, and furthermore because
$\|\mathbf{x}(t)-\boldsymbol{a}(t)\|^{2}\leq\|\mathbf{x}(t)-\boldsymbol{\mu}\|^{2}$
for $\boldsymbol{a}(t)=\mbox{Proj}_{\boldsymbol{A}}\mathbf{x}(t)$.

Next, let us estimate
$\Pr\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\boldsymbol{A}\\}$. The idea is
that due to the fact that $\mathbf{x}(t)$ is sampled from a Gaussian
distribution with mean $\boldsymbol{\mu}(t)$, $\mathbf{x}(t)$ concentrates
around its mean $\boldsymbol{\mu}(t)$ with high probability. Since the mean is
projected onto a shrunk version of the set $\boldsymbol{A}$, namely,
$(1-r_{t})\boldsymbol{A}$, by appropriately tuning $r_{t}$ and $\sigma_{t}$ we
can ensure that $\mathbf{x}(t)$ stays within the original feasible set with
high probability. Let
$\mathcal{O}_{r_{t}}(\boldsymbol{\mu})=\\{\boldsymbol{y}\in\mathbb{R}^{Nd}|\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}<r^{2}_{t}\\}$
denote the $r_{t}$-neighborhood of the point $\boldsymbol{\mu}$. Hence,
$\sup_{\boldsymbol{y}\notin\mathcal{O}_{r_{t}}(\boldsymbol{\mu})}-\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}=-r^{2}_{t}$
. Then, taking into account the fact that
$\mathcal{O}_{r_{t}}(\boldsymbol{\mu})$ is contained in $\boldsymbol{A}$ and
$r_{t}<1$, we obtain that for any $t$ and any bounded $\sigma>\sigma_{t}$:

$\displaystyle\Pr$
$\displaystyle\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\boldsymbol{A}\\}\leq\Pr\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\mathcal{O}_{r_{t}}(\boldsymbol{\mu})\\}$
(73)
$\displaystyle=\int_{\boldsymbol{y}\notin\mathcal{O}_{r_{t}}(\boldsymbol{\mu})}\frac{1}{(2\pi)^{Nd/2}\sigma_{t}^{Nd}}\exp\left\\{-\frac{\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}}{2\sigma_{t}^{2}}\right\\}d\boldsymbol{y}$
(74)
$\displaystyle=\int_{\boldsymbol{y}\notin\mathcal{O}_{r_{t}}(\boldsymbol{\mu})}\exp\left\\{-\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}\left(\frac{1}{2\sigma_{t}^{2}}-\frac{1}{2\sigma^{2}}\right)\right\\}$
(75)
$\displaystyle\qquad\qquad\times\frac{\sigma^{Nd}}{\sigma_{t}^{Nd}}\frac{1}{(2\pi)^{{Nd}/2}\sigma^{Nd}}\exp\left\\{-\frac{\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}}{2\sigma^{2}}\right\\}d\boldsymbol{y}$
(76)
$\displaystyle\leq\exp\left\\{-r_{t}\left(\frac{1}{2\sigma_{t}^{2}}-\frac{1}{2\sigma^{2}}\right)\right\\}\frac{\sigma^{Nd}}{\sigma_{t}^{Nd}}$
(77)
$\displaystyle\quad\times\int_{\boldsymbol{y}\notin\mathcal{O}_{r_{t}}(\boldsymbol{\mu})}\frac{1}{(2\pi)^{Nd/2}\sigma^{Nd}}\exp\left\\{-\frac{\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}}{2\sigma^{2}}\right\\}d\boldsymbol{y}$
(78) $\displaystyle\leq
k_{2}\frac{e^{-\frac{r^{2}_{t}}{2\sigma_{t}^{2}}}}{\sigma_{t}^{Nd}}$ (79)

for some finite $k_{2}>0$. The last inequality holds because

$\int_{\boldsymbol{y}\notin\mathcal{O}_{r_{t}}(\boldsymbol{\mu})}\frac{1}{(2\pi)^{{Nd}/2}\sigma^{Nd}}\exp\left\\{-\frac{\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}}{2\sigma^{2}}\right\\}d\boldsymbol{y}\leq
1$

and, thus, there exists $0<k_{3}<\infty$:

$\int_{\boldsymbol{y}\notin\mathcal{O}_{r_{t}}(\boldsymbol{\mu})}\frac{e^{\frac{r^{2}_{t}}{2\sigma^{2}}}\sigma^{n}}{(2\pi)^{n/2}\sigma^{Nd}}\exp\left\\{-\frac{\|\boldsymbol{y}-\boldsymbol{\mu}\|^{2}}{2\sigma^{2}}\right\\}d\boldsymbol{y}\leq
k_{3}.$

From (39) it now remains to bound the term

$\mathrm{E}\\{(\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)),\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}.$

According to definitions of $\boldsymbol{P}_{i}$ and $\boldsymbol{R}_{i}$,
Remark 1, and the Cauchy-Schwarz inequality,

$\displaystyle\mathrm{E}\\{(\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)),\boldsymbol{P}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t)))|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(80)
$\displaystyle=\mathrm{E}(J_{i}(\mathbf{x}(t))\frac{\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}}{\sigma^{2}_{t}}-\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu},\sigma_{t}),$
$\displaystyle\quad\qquad\frac{\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}}{\sigma^{2}_{t}}(J_{i}(\boldsymbol{a}(t))-J_{i}(\mathbf{x}(t))))$
$\displaystyle\leq\|\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu},\sigma_{t})\|\mathrm{E}\|P_{i}(\mathbf{x}_{t},\boldsymbol{\mu})\|$
$\displaystyle\,-\mathrm{E}\left\\{\frac{J_{i}(\mathbf{x}(t))(J_{i}(\boldsymbol{a}(t))-J_{i}(\mathbf{x}(t)))\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}}{\sigma^{4}_{t}}\right\\}$
$\displaystyle\leq\|\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu},\sigma_{t})\|\mathrm{E}\|P_{i}(\mathbf{x}_{t},\boldsymbol{\mu})\|$
$\displaystyle\quad\qquad+\Pr\\{\mathbf{x}(t)\in\mathbb{R}^{Nd}\setminus\boldsymbol{A}\\}\times$
$\displaystyle\times\mathrm{E}\left\\{\frac{|J_{i}(\mathbf{x}(t))|\|\mathbf{x}(t)-\boldsymbol{\mu}\|\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}}{\sigma^{4}_{t}}\right\\}.$

Note that $\|\tilde{\boldsymbol{M}}_{i}(\boldsymbol{\mu},\sigma_{t})\|$ is
bounded from Lemma 2), and that

$\displaystyle\mathrm{E}\left\\{\frac{|J_{i}(\mathbf{x}(t))|\|\mathbf{x}(t)-\boldsymbol{\mu}\|\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}}{\sigma^{4}_{t}}\right\\}=\frac{h_{i}(\boldsymbol{\mu},\sigma_{t})}{\sigma^{4}_{t}},$
(81)

where $h_{i}(\boldsymbol{\mu},\sigma_{t})$ is a quadratic function of
$\boldsymbol{\mu}$ and $\sigma_{t}$, $i\in[N]$ (see Appendix -F for more
details). Hence, due to the choice of the parameters $r_{t}$ and $\sigma_{t}$
(in particular, Assumption 6 d)) and the estimations in (68)- (80), we
conclude that the terms containing $\boldsymbol{P}_{i}$ are dominated by other
terms in the inequality in (39). Thus, by inserting (57)-(IV-C) into (39), we
obtain

$\displaystyle\mathrm{E}\\{\|$
$\displaystyle\boldsymbol{\mu}^{i}(t+1)-\mathbf{z}^{i}(t+1)\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}$
(82)
$\displaystyle\leq(1-2\beta_{t}\epsilon_{t})\|\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t+1)\|^{2}$
(83)
$\displaystyle\quad-2\beta_{t}(\boldsymbol{M}_{i}(\boldsymbol{\mu})-\boldsymbol{M}_{i}(\mathbf{z}(t)),\boldsymbol{\mu}^{i}-\mathbf{z}^{i}(t))$
(84) $\displaystyle\quad+2\beta_{t}O(\sigma_{t})(1+V(t,\boldsymbol{\mu}))$
(85) $\displaystyle\quad{+\beta^{2}_{t}\epsilon^{2}(t)V(t,\boldsymbol{\mu})}$
(86) $\displaystyle\quad+O(\gamma^{2}_{t})(1+V(t,\boldsymbol{\mu})),$ (87)

where in the last inequality we used the fact that $\epsilon_{t}\to 0$
(Assumption 6 a)), $\gamma_{t}\to 0$, and $\sigma_{t}\to 0$ for all $i\in[N]$
as $t\to\infty$ (Assumption 6 c), d)). Thus, taking into account Assumption 6
c), d) and (82), we obtain

$\displaystyle\mathrm{E}[$
$\displaystyle\|\boldsymbol{\mu}(t+1)-\mathbf{z}(t+1)\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}]$
(88) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\mathrm{E}[\|\boldsymbol{\mu}^{i}(t+1)-\mathbf{z}^{i}(t+1)\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}]$
(89) $\displaystyle\leq$
$\displaystyle(1-2\epsilon_{t}\beta_{t})\|\boldsymbol{\mu}-\mathbf{z}(t)\|^{2}-2\beta_{t}(\boldsymbol{M}(\boldsymbol{\mu})-\boldsymbol{M}(\mathbf{z}(t)),\boldsymbol{\mu}-\mathbf{z}(t))$
(90)
$\displaystyle+O(\beta_{t}\sigma_{t}+\gamma^{2}_{t})(1+V(t,\boldsymbol{\mu})).$
(91)

Thus,

$\displaystyle LV(t,\boldsymbol{\mu})\leq$
$\displaystyle-2\beta_{t}(\boldsymbol{M}(\boldsymbol{\mu})-\boldsymbol{M}(\mathbf{z}(t)),\boldsymbol{\mu}-\mathbf{z}(t))$
(92)
$\displaystyle+O(\beta_{t}\sigma_{t}+\gamma^{2}_{t})(1+V(t,\boldsymbol{\mu})).$
(93)

According to Assumption 6 b)-c),
$\sum_{t=0}^{\infty}\beta_{t}\sigma_{t}+\gamma^{2}_{t}<\infty$. Furthermore,
from Assumption 6 a) $\sum_{t=0}^{\infty}\beta_{t}=\infty$. Taking into
account this, (92), and monotonicity of $\boldsymbol{M}$ implying

$\displaystyle(\boldsymbol{M}(\boldsymbol{\mu})-\boldsymbol{M}(\mathbf{z}(t)),\boldsymbol{\mu}-\mathbf{z}(t))\geq
0,\;\forall t,\,\forall\boldsymbol{\mu}\in\boldsymbol{A},$ (94)

we conclude that $LV(t,\boldsymbol{\mu})$ satisfies the decay needed for the
application of Theorem 2.5.2 in [15] and consequently, $\boldsymbol{\mu}(t)$
is finite almost surely for any $t\in\mathbb{Z}_{+}$ irrespective of
$\boldsymbol{\mu}(0)$. ∎

###  IV-D Convergence to Nash equilibrium

We will use the bound estimations in the previous section to prove convergence
of the algorithm. In particular, we use Inequality (88), which bounds the
decay of the sequence
$\mathrm{E}[\|\boldsymbol{\mu}(t+1)-\mathbf{z}(t)\|^{2}|\boldsymbol{\mu}(t)]$
in terms of $\|\boldsymbol{\mu}-\mathbf{z}(t+1)\|^{2}$. We will show that this
decay satisfies the conditions for applying Lemma 10 in [17]. From this, it
can readily be inferred that random variables
$\|\boldsymbol{\mu}(t)-\mathbf{z}(t)\|$ converge to zero.

First, however, let us verify that even in the compact action case,
Inequalities (88), (92) hold.

###### Remark 2

If the set $\boldsymbol{A}$ is compact, due to Assumption 5, the inequality
(IV-C) can be replaced by

$\displaystyle\mathrm{E}\\{\|\boldsymbol{R}_{i}(\mathbf{x}(t),\boldsymbol{\mu}(t))\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}\\}=O\left(\frac{1}{\sigma^{4}_{t}}\right).$

Moreover, the inequalities (68) and (80) hold for the case of the bounded set
$\boldsymbol{A}$. Indeed, due to polynomial behavior of $J_{i}(\mathbf{x}(t))$
for large $\mathbf{x}(t)$, the terms
$\mathrm{E}\frac{\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}|J_{i}(\mathbf{x}(t))-J_{i}(\boldsymbol{a}(t))|^{2}}{\sigma^{4}_{t}}$
and
$\mathrm{E}\left\\{\frac{|J_{i}(\mathbf{x}(t))(J_{i}(\boldsymbol{a}(t))-J_{i}(\mathbf{x}(t)))|\|\mathbf{x}^{i}(t)-\boldsymbol{\mu}^{i}\|^{2}}{\sigma^{4}_{t}}\right\\}$
are upper bounded by some constants. Thus, for the bounded set
$\boldsymbol{A}$, the inequality (88) can be rewritten as

$\displaystyle\mathrm{E}[$
$\displaystyle\|\boldsymbol{\mu}(t+1)-\mathbf{z}(t+1)\|^{2}|\boldsymbol{\mu}(t)=\boldsymbol{\mu}]$
$\displaystyle\leq$
$\displaystyle(1-2\epsilon_{t}\beta_{t})\|\boldsymbol{\mu}-\mathbf{z}(t)\|^{2}-2\beta_{t}(\boldsymbol{M}(\boldsymbol{\mu})-\boldsymbol{M}(\mathbf{z}(t)),\boldsymbol{\mu}-\mathbf{z}(t))$
$\displaystyle+O(\beta_{t}\sigma_{t}+\gamma^{2}_{t}).$

###### Proof:

(of Theorem 2) Note that we can rewrite (88) as:

$\displaystyle\mathrm{E}[\|\boldsymbol{\mu}(t+1)-\mathbf{z}(t+1)\|^{2}|\mathcal{F}_{t}]$
(95) $\displaystyle\leq$
$\displaystyle(1-2\epsilon_{t}\beta_{t})\|\boldsymbol{\mu}(t)-\mathbf{z}(t)\|^{2}+O(\gamma^{2}_{t}+\beta_{t}\sigma_{t}),$
(96)

where $\mathcal{F}_{t}$ is the $\sigma$-algebra generated by the random
variables $\\{\mathbf{x}(k),\boldsymbol{\mu}(k)\\}_{k=0}^{t}$. In (95) we used
(94) and Lemma 5.

From Assumption 6, and the choices of $\gamma_{t}$, $\sigma_{t}$,
$\epsilon_{t}$, we get
$O(\gamma^{2}_{t}+\beta_{t}\sigma_{t})=O(\frac{1}{t^{n}})$,
$\epsilon_{t}\beta_{t}=\frac{1}{t^{m}}$, with $n>1$, $m\leq 1$. Thus,

$\lim_{t\to\infty}\frac{O(\gamma^{2}_{t}+\beta_{t}\sigma_{t})}{\epsilon_{t}\beta_{t}}=0.$

Assumption 6 d), the fact that
$\sum_{t=0}^{\infty}\gamma^{2}_{t}+\beta_{t}\sigma_{t}<\infty$ and the above
result in the decay (95) imply that we can apply Lemma 10 in [17] to the
sequence $\|\boldsymbol{\mu}(t+1)-\mathbf{z}(t+1)\|^{2}$ to conclude its
almost sure convergence to $0$ as $t\to\infty$. Next, by taking into account
Theorem 3 and Theorem 1, we obtain that

$\Pr\\{\lim_{t\to\infty}\boldsymbol{\mu}(t)=\boldsymbol{a}^{*}\\}=1,$

where $\boldsymbol{a}^{*}$ is the least norm Nash equilibrium in the game
$\Gamma(N,\\{A_{i}\\},\\{J_{i}\\})$. Finally, Assumption 6 implies that
$\lim_{t\to\infty}\sigma_{t}=0$. Taking into account that
$\mathbf{x}(t)\sim\mathcal{N}(\boldsymbol{\mu}(t),\sigma_{t})$ and
$\lim_{t\to\infty}\|\boldsymbol{a}(t)-\mathbf{x}(t)\|=0$, we conclude that
$\boldsymbol{a}(t)$ converges weakly to a Nash equilibrium
$\boldsymbol{a}^{*}=\boldsymbol{\mu}^{*}$. Moreover, according to Portmanteau
Lemma [10], this convergence is also in probability.

∎

##  V Discussion

In the proposed algorithm convergence is established under mild conditions as
strict monotonicity of the game mapping is not implied. This significantly
extends the applicability of bandit online learning. For example, the zero-sum
game considered in [13] with an interior Nash equilibrium satisfies the
assumption of our theorem. Whereas all the follow-the-regularized leader
(FTRL) learning approach fails to converge in a simple zero-sum game (such as
matching penny), our doubly regularized approach can resolve this problem. In
general, examples of games that satisfy assumptions above include mixed
extensions of zero-sum games, Cournot competition, continuous-action
congestion games and convex potential games. On the other hand, mixed
extensions of non-zero sum games do not satisfy the monotonicity assumption in
general.

In accordance with the payoff-based information structure, the parameters
$\gamma_{t},\sigma_{t},\epsilon_{t},r_{t}$ are independent of the problem data
including the Lipschitz constant of the game mapping or the constraint sets.
Below, we further specify feasible choices of the parameters to ensure
convergence.

###### Lemma 6

A sufficient condition on $0<a_{1},a_{2},a_{3},a_{4}<1$ for satisfying
Assumption 6 is as follows:

i) $a_{1}+2a_{2}<1$, $a_{1}+2a_{2}+a_{3}<1$.

ii) $a_{1}+2a_{2}+a_{3}<1$, $a_{1}+2a_{2}+6a_{3}-2a_{4}<1$.

iii) $2a_{1}>1$, $a_{1}+3a_{2}>1$.

iv) $a_{3}<a_{4}<a_{2}$.

###### Proof:

The series $\sum_{t=0}^{\infty}1/t^{m}$ converges for $m>1$ and diverges
otherwise. Thus, the statements i), iii), iv) above follow. To show statement
ii), let us consider the term $(\epsilon_{t}-\epsilon_{t-1})^{2}$ in the first
summand of b), namely,
$\sum_{t=0}^{\infty}\frac{(\epsilon_{t}-\epsilon_{t-1})^{2}}{\beta_{t}\epsilon^{3}_{t}}$:

$\displaystyle(\epsilon_{t}-\epsilon_{t-1})^{2}=(t^{-a_{3}}-(t-1)^{-a_{3}})^{2}\quad\text{(multiply
by $\frac{t^{2a_{3}}}{t^{2a_{3}}}$)}$
$\displaystyle={\big{(}(1-1/t)^{-a_{3}}-1\big{)}^{2}}/{t^{2a_{3}}}\quad\quad\text{(do
Taylor approximation)}$
$\displaystyle={(1+a_{3}/t+O(t^{-2})-1)^{2}}/{t^{2a_{3}}}=O(t^{-2-2a_{3}}).$

Combining the above with the denominator $\beta_{t}\epsilon_{t}^{3}$, we
obtain that
$\sum_{t}\frac{(\epsilon_{t}-\epsilon_{t-1})^{2}}{\beta_{t}\epsilon^{3}_{t}}$
converges if $a_{3}1+2a_{2}+a_{3}<1$. Repeating the same analysis for
$\frac{(r_{t}-r_{t-1})^{2}}{\beta_{t}\epsilon^{6}_{t}}$, we obtain
$\sum_{t}\frac{(r_{t}-r_{t-1})^{2}}{\beta_{t}\epsilon^{6}_{t}}$ converges if
$a_{1}+2a_{2}+6a_{3}-2a_{4}<1$ and ii) is verified. ∎

##  VI Conclusions

We designed an algorithm for learning Nash equilibria in convex games with
monotone game mappings using online bandit feedback information. Our algorithm
relied on a suitable double regularization to handle non-strictly monotone
game maps as well as feasibility of the queried actions (online setting). The
implications of our result is that players can learn Nash equilibria in
several monotone games such as finite action zero-sum games, infinite action
zero-sum convex games and convex games with linear coupling constraints.
Several points remain open and are topic of our current study. These include
showing that our algorithm is no-regret, unifying different sampling
approaches to perform one-point estimation of the game mapping for bandit
learning in games, and analyzing the convergence rate of the algorithm.

## References

  * [1] J. P. Bailey, G. Gidel, and G. Piliouras.  Finite regret and cycles with fixed step-size via alternating gradient descent-ascent.  In Conference on Learning Theory, pages 391–407, 2020.
  * [2] D Balduzzi, S Racaniere, J Martens, J Foerster, K Tuyls, and T Graepel.  The mechanics of n-player differentiable games.  In ICML, volume 80, pages 363–372. JMLR. org, 2018.
  * [3] B. Bharath and V. S. Borkar.  Stochastic approximation algorithms: Overview and recent trends.  Sadhana, 24(4):425–452, 1999.
  * [4] M. Bravo, D. Leslie, and P. Mertikopoulos.  Bandit learning in concave n-person games.  In Advances in Neural Information Processing Systems, pages 5661–5671, 2018.
  * [5] F. Facchinei and C. Kanzow.  Generalized Nash equilibrium problems.  4OR, 5(3):173–210, 2007.
  * [6] A. D. Flaxman, A. T. Kalai, and H. B. McMahan.  Online convex optimization in the bandit setting: gradient descent without a gradient.  In Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithms, pages 385–394. Society for Industrial and Applied Mathematics, 2005.
  * [7] S. Grammatico.  Comments on distributed robust adaptive equilibrium computation for generalized convex games (automatica 63(2016) 82-91).  Automatica, 97:186 – 188, 2018.
  * [8] J.-B. Hiriart-Urruty and C. Lemaréchal.  Fundamentals of Convex Analysis.  Grundlehren Text Editions. Springer, 2001.
  * [9] R. D. Kleinberg.  Nearly tight bounds for the continuum-armed bandit problem.  In Advances in Neural Information Processing Systems, pages 697–704, 2005.
  * [10] A. Klenke.  Probability theory: a comprehensive course.  Springer, London, 2008.
  * [11] J. Koshal, A Nedić, and U. Shanbhag.  Single timescale regularized stochastic approximation schemes for monotone nash games under uncertainty.  In IEEE Conference on Decision and Control, pages 231–236, 2010\.
  * [12] P. Mertikopoulos, B. Lecouat, H. Zenati, Ch.-Sh. Foo, V. Chandrasekhar, and G. Piliouras.  Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile.  arXiv preprint arXiv:1807.02629, 2018.
  * [13] P. Mertikopoulos, C. Papadimitriou, and G. Piliouras.  Cycles in adversarial regularized learning.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 2703–2717, 2018.
  * [14] Y. Nesterov and V. Spokoiny.  Random gradient-free minimization of convex functions.  Found. Comput. Math., 17(2):527–566, April 2017.
  * [15] M. B. Nevelson and R. Z. Khasminskii.  Stochastic approximation and recursive estimation [translated from the Russian by Israel Program for Scientific Translations ; translation edited by B. Silver].  American Mathematical Society, 1973.
  * [16] J.-S. Pang and F. Facchinei.  Finite-dimensional variational inequalities and complementarity problems : vol. 1.  Springer series in operations research. Springer, New York, Berlin, Heidelberg, 2003.
  * [17] B. T. Poljak.  Introduction to optimization.  Optimization Software, 1987.
  * [18] T. Tatarenko and M. Kamgarpour.  Learning generalized nash equilibria in a class of convex games.  IEEE Transactions on Automatic Control, 2018.  to appear.  URL: https://arxiv.org/abs/1703.04113.
  * [19] T. Tatarenko and M. Kamgarpour.  Learning nash equilibria in monotone games.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 3104–3109, 2019.
  * [20] A. L. Thathachar and P. S. Sastry.  Networks of Learning Automata: Techniques for Online Stochastic Optimization.  Springer US, 2003.
  * [21] V.A. Zorich and R. Cooke.  Mathematical Analysis II.  Mathematical Analysis. Springer, 2004.

The appendix provides supporting theorems and proof of certain lemmas and
statements.

###  -A Supporting Theorems

Let $\\{\mathbf{X}(t)\\}_{t}$, $t\in\mathbb{Z}_{+}$, be a discrete-time Markov
process on some state space $E\subseteq\mathbb{R}^{d}$, namely
$\mathbf{X}(t)=\mathbf{X}(t,\omega):\mathbb{Z}_{+}\times\Omega\to E$, where
$\Omega$ is the sample space of the probability space on which the process
$\mathbf{X}(t)$ is defined. The transition function of this chain, namely
$\Pr\\{\mathbf{X}(t+1)\in\Gamma|\mathbf{X}(t)=\mathbf{X}\\}$, is denoted by
$P(t,\mathbf{X},t+1,\Gamma)$, $\Gamma\subseteq E$.

###### Definition 4

The operator $L$ defined on the set of measurable functions
$V:\mathbb{Z}_{+}\times E\to\mathbb{R}$, $\mathbf{X}\in E$, by

$\displaystyle LV(t,\mathbf{X})$
$\displaystyle=\int{P(t,\mathbf{X},t+1,dy)[V(t+1,y)-V(t,\mathbf{X})]}$
$\displaystyle=E[V(t+1,\mathbf{X}(t+1))\mid\mathbf{X}(t)=\mathbf{X}]-V(t,\mathbf{X}),$

is called a _generating operator_ of a Markov process
$\\{\mathbf{X}(t)\\}_{t}$.

Next, we formulate the following theorem for discrete-time Markov processes,
which is proven in [15], Theorem 2.5.2.

###### Theorem 4

Consider a Markov process $\\{\mathbf{X}(t)\\}_{t}$ and suppose that there
exists a function $V(t,\mathbf{X})\geq 0$ such that $\inf_{t\geq
0}V(t,\mathbf{X})\to\infty$ as $\|\mathbf{X}\|\to\infty$ and

$LV(t,\mathbf{X})\leq-\alpha(t+1)\psi(t,\mathbf{X})+f(t)(1+V(t,\mathbf{X})),$

where $\psi\geq 0$ on $\mathbb{R}\times\mathbb{R}^{d}$, $f(t)>0$,
$\sum_{t=0}^{\infty}f(t)<\infty$. Let $\alpha(t)$ be such that $\alpha(t)>0$,
$\sum_{t=0}^{\infty}\alpha(t)=\infty$. Then, almost surely $\sup_{t\geq
0}\|\mathbf{X}(t,\omega)\|=R(\omega)<\infty$.

The following result related to the convergence of the stochastic process is
proven in Lemma 10 (page 49) in [17].

###### Theorem 5

Let $v_{0},\ldots,v_{k}$ be a sequence of random variables, $v_{k}\geq 0$,
$\mathrm{E}v_{0}<\infty$ and let

$\mathrm{E}\\{v_{k+1}|\mathcal{F}_{k}\\}\leq(1-\alpha_{k})v_{k}+\beta_{k},$

where $\mathcal{F}_{k}$ is the $\sigma$-algebra generated by the random
variables $\\{v_{0},\ldots,v_{k}\\}$, $0<\alpha_{k}<1$,
$\sum_{k=0}^{\infty}\alpha_{k}=\infty$, $\beta_{k}\geq 0$,
$\sum_{k=0}^{\infty}\beta_{k}<\infty$,
$\lim_{k\to\infty}\frac{\beta_{k}}{\alpha_{k}}=0$. Then $v_{k}\to 0$ almost
surely, $\mathrm{E}v_{k}\to 0$ as $k\to\infty$.

###  -B Proof of Lemma 2

###### Proof:

We verify that the differentiation under the integral sign in (5) is
justified. It can then readily be verified that (2) holds, by taking the
differentiation inside the integral. A sufficient condition for
differentiation under the integral is that the integral of the formally
differentiated function with respect to $\mu^{i}_{k}$ converges uniformly,
whereas the differentiated function is continuous (see [21], Chapter 17). By
formally differentiating the function under the integral sign and omitting the
arguments $t$, we obtain

$\displaystyle\frac{1}{\sigma^{2}}\int_{\mathbb{R}^{Nd}}J_{i}(\boldsymbol{x})(x^{i}_{k}-\mu^{i}_{k})p(\boldsymbol{\mu},\boldsymbol{x},\sigma)d\boldsymbol{x}.$
(97)

Given Assumption 1,
$J_{i}(\boldsymbol{x})(x^{i}_{k}-\mu^{i}_{k})p(\boldsymbol{\mu},\boldsymbol{x},\sigma)$
is continuous. Thus, it remains to check that the integral of this function
converges uniformly with respect to any $\boldsymbol{\mu}\in\boldsymbol{A}$.
If $\boldsymbol{A}$ is bounded, then the conclusion follows from the
polynomial behavior of the function $J_{i}$ on the infinity.

Now, we move to the case of the unbounded set $\boldsymbol{A}$. To this end,
we can write the Taylor expansion of the function $J_{i}$ around the point
$\boldsymbol{\mu}(i,k)\in\mathbb{R}^{Nd}$ with the coordinates
$\mu(i,k)^{i}_{k}=\mu^{i}_{k}$ and $\mu(i,k)^{j}_{m}=x^{j}_{m}$ for any $j\neq
i$, $m\neq k$, in the integral (97):

$\displaystyle\int_{\mathbb{R}^{Nd}}J_{i}(\boldsymbol{x})(x^{i}_{k}-\mu^{i}_{k})p(\boldsymbol{\mu},\boldsymbol{x},\sigma)d\boldsymbol{x}$
$\displaystyle=\int_{\mathbb{R}^{Nd}}[J_{i}(\boldsymbol{\mu}(i,k))$
$\displaystyle\quad+\frac{\partial
J_{i}(\boldsymbol{\eta}(\boldsymbol{x},\boldsymbol{\mu}))}{\partial
x^{i}_{k}}(x^{i}_{k}-\mu^{i}_{k})](x^{i}_{k}-\mu^{i}_{k})p(\boldsymbol{\mu},\boldsymbol{x},\sigma)d\boldsymbol{x}$
$\displaystyle=\int_{\mathbb{R}^{Nd}}\frac{\partial
J_{i}(\boldsymbol{\eta}(\boldsymbol{x},\boldsymbol{\mu}))}{\partial
x^{i}_{k}}(x^{i}_{k}-\mu^{i}_{k})^{2}p(\boldsymbol{\mu},\boldsymbol{x},\sigma)d\boldsymbol{x}$
$\displaystyle=\int_{\mathbb{R}^{Nd}}\frac{\partial
J_{i}(\boldsymbol{\eta}_{1}(\boldsymbol{y},\boldsymbol{\mu}))}{\partial
x^{i}_{k}}(y^{i}_{k})^{2}p(\boldsymbol{0},\boldsymbol{y},\sigma)d\boldsymbol{y},$

where
$\boldsymbol{\eta}(\boldsymbol{x},\boldsymbol{\mu})=\boldsymbol{\mu}(i,k)+\theta(\boldsymbol{x}-\boldsymbol{\mu}(i,k))$,
$\theta\in(0,1)$, $\boldsymbol{y}=\boldsymbol{x}-\boldsymbol{\mu}(i,k)$,
${\boldsymbol{\eta}_{1}}(\boldsymbol{y},\boldsymbol{\mu})=\boldsymbol{\mu}(i,k)+\theta\boldsymbol{y}$.
The uniform convergence of the integral above and, in particular, its
boundedness, follows from the fact333see the basic sufficient condition using
majorant [21], Chapter 17.2.3. that under Assumption 5, $\frac{\partial
J_{i}(\boldsymbol{\eta}_{1}(\boldsymbol{y},\boldsymbol{\mu}))}{\partial
x^{i}_{k}}\leq l^{i}_{k}$ for some positive constant $l^{i}_{k}$ and for all
$i\in[N]$ and $k\in[d]$. Thus, $|\frac{\partial
J_{i}(\boldsymbol{\eta}_{1}(\boldsymbol{y},\boldsymbol{\mu}))}{\partial
x^{i}_{k}}(y^{i}_{k})^{2}p(\boldsymbol{0},\boldsymbol{y},\sigma)|\leq
h(\boldsymbol{y})=l(y^{i}_{k})^{2}p(\boldsymbol{0},\boldsymbol{y},\sigma),$where
$\int_{\mathbb{R}^{Nd}}h(\boldsymbol{y})d\boldsymbol{y}<\infty$.∎

###  -C Proof of Lemma 3

###### Proof:

Without loss of generality, assume $\mathbf{x}\notin(1-r_{t-1})\boldsymbol{A}$
(otherwise,
$\|\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}\mathbf{x}-\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}\mathbf{x}\|=0$).
Due to convexity of the set $\boldsymbol{A}\subset\mathbb{R}^{Nd}$ there
exists a convex function $g:\mathbb{R}^{Nd}\to\mathbb{R}$ such that
$\boldsymbol{A}=\\{\mathbf{x}:g(\mathbf{x})\leq 0\\}$, whereas
$(1-r_{t})\boldsymbol{A}=\\{\mathbf{x}:g(\mathbf{x})\leq-r_{t}\\}$ for any
$t$. Moreover, since
$\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}\mathbf{x}=\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}\\{\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}\mathbf{x}\\}$,
we have
$\|\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}\mathbf{x}-\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}\mathbf{x}\|=d,$
where

$\displaystyle d=\min_{\mathbf{y}}$
$\displaystyle\|\mathbf{y}-\mathbf{x}^{\prime}\|,\,\,\mathbf{x}^{\prime}=\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}\mathbf{x},$
s.t. $\displaystyle g(\mathbf{y})=-r_{t-1}.$

The optimization problem has a solution $\mathbf{y}^{*}$ for which the
gradient of the corresponding Lagrangian is zero, namely

$\frac{(\mathbf{y}^{*}-\mathbf{x}^{\prime})}{\|\mathbf{y}^{*}-\mathbf{x}^{\prime}\|}+\lambda\nabla
g(\mathbf{y}^{*})=\boldsymbol{0},$

where $\lambda>0$ is the dual multiplier of the problem under consideration.
Notice that due to Assumption 1 and the choice of $r_{1}$ the Slator’s
condition for the constraints $g(\mathbf{x})\leq-r_{t}$ holds for all $t$.
Hence, for any $\mathbf{x}\in\mathbb{R}^{Nd}$ there exists a constant
$\Lambda>0$ such that $\lambda<\Lambda$ (see [8]). Thus, we conclude that

$\nabla
g(\mathbf{y}^{*})=-\frac{(\mathbf{y}^{*}-\mathbf{x}^{\prime})}{\lambda\|\mathbf{y}^{*}-\mathbf{x}^{\prime}\|}.$

Next, due to convexity of the function $g$,

$\displaystyle g(\mathbf{x}^{\prime})$ $\displaystyle\geq
g(\mathbf{y}^{*})+(\nabla
g(\mathbf{y}^{*}),\mathbf{x}^{\prime}-\mathbf{y}^{*})$
$\displaystyle=-r_{t-1}+\frac{\|\mathbf{y}^{*}-\mathbf{x}^{\prime}\|^{2}}{\lambda|\mathbf{y}^{*}-\mathbf{x}^{\prime}\|}\geq-
r_{t-1}+\frac{\|\mathbf{y}^{*}-\mathbf{x}^{\prime}\|}{\Lambda}.$

Thus, taking into account that $g(\mathbf{x}^{\prime})\leq-r_{t}$, we obtain

$d=\|\mathbf{y}^{*}-\mathbf{x}^{\prime}\|\leq\Lambda(r_{t-1}-r_{t})=O(|r_{t-1}-r_{t}|).$

∎

###  -D Proof of Lemma 4

###### Proof:

Let us use $\theta=\epsilon_{t}^{3/2}$ in Lemma 1 to express
$\boldsymbol{y}(t)$ as
$\boldsymbol{y}(t)=\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\boldsymbol{y}(t)-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t))]$.
Using this equivalence, the triangle inequality and non-expansion of
projection operator we obtain

$\displaystyle\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|$
$\displaystyle\leq\|\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\boldsymbol{y}(t)-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t))]$
$\displaystyle-\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\boldsymbol{y}(t-1)-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t-1))$
$\displaystyle+\epsilon_{t-1}\boldsymbol{y}(t-1))]\|+\|\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\boldsymbol{y}(t-1)$
$\displaystyle-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t-1))+\epsilon_{t-1}\boldsymbol{y}(t-1))]$
$\displaystyle-\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}[\boldsymbol{y}(t-1)-$
$\displaystyle\quad\quad\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t-1))+\epsilon_{t-1}\boldsymbol{y}(t-1))]\|$
$\displaystyle\leq\|(1-\epsilon^{3/2}_{t}\epsilon_{t})(\boldsymbol{y}(t)-\boldsymbol{y}(t-1))$
$\displaystyle-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t))-\boldsymbol{M}(\boldsymbol{y}(t-1)))$
$\displaystyle+\boldsymbol{y}(t-1)\epsilon^{3/2}_{t}(\epsilon_{t-1}-\epsilon_{t})\|$
$\displaystyle+\|\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\tilde{\boldsymbol{y}}(t-1)]-\mbox{Proj}_{(1-r_{t-1})\boldsymbol{A}}[\tilde{\boldsymbol{y}}(t-1)]\|,$

where
$\tilde{\boldsymbol{y}}(t-1)=\boldsymbol{y}(t-1)-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t-1))+\epsilon_{t-1}\boldsymbol{y}(t-1))$.
Next, due to Lemma 3 we have for any $\theta_{t}>0$ and $\kappa_{t}>0$

$\displaystyle\|$
$\displaystyle\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}\leq(1+\theta_{t})\|(1-\epsilon^{3/2}_{t}\epsilon_{t})(\boldsymbol{y}(t)-\boldsymbol{y}(t-1))$
$\displaystyle-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t))-\boldsymbol{M}(\boldsymbol{y}(t-1)))$
$\displaystyle+\boldsymbol{y}(t-1)\epsilon^{3/2}_{t}(\epsilon_{t-1}-\epsilon_{t})\|^{2}+(1+1/{\theta_{t}})O((r_{t}-r_{t-1})^{2})$
$\displaystyle\leq(1+\theta_{t})(1+\kappa_{t})\|(1-\epsilon^{3/2}_{t}\epsilon_{t})(\boldsymbol{y}(t)-\boldsymbol{y}(t-1))$
$\displaystyle-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t))-\boldsymbol{M}(\boldsymbol{y}(t-1)))\|^{2}$
$\displaystyle+(1+1/{\kappa_{t}})\|\boldsymbol{y}(t-1)\|^{2}\epsilon^{3}_{t}(\epsilon_{t-1}-\epsilon_{t})^{2}$
$\displaystyle+(1+1/{\theta_{t}})O((r_{t}-r_{t-1})^{2}).$

Furthermore, there exists $T$ such that for all $t>T$

$\displaystyle\|(1-\epsilon^{3/2}_{t}\epsilon_{t})(\boldsymbol{y}(t)-\boldsymbol{y}(t-1))$
$\displaystyle\qquad\qquad-\epsilon^{3/2}_{t}(\boldsymbol{M}(\boldsymbol{y}(t))-\boldsymbol{M}(\boldsymbol{y}(t-1)))\|^{2}$
$\displaystyle\leq(1-\epsilon^{3/2}_{t}\epsilon_{t})^{2}\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle\qquad\qquad+\epsilon^{3}_{t}\|\boldsymbol{M}(\boldsymbol{y}(t))-\boldsymbol{M}(\boldsymbol{y}(t-1))\|^{2}$
$\displaystyle\leq(1-2\epsilon^{3/2}_{t}\epsilon_{t}+\epsilon^{3}_{t}(\epsilon^{2}_{t}+L^{2}))\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle\leq(1-\epsilon^{3/2}_{t}\epsilon_{t})\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2},$

where the first inequality is due to monotonicity of the mapping
$\boldsymbol{M}$, the second one is as it is Lipschitz continuous, and the
third one is due to the fact that
$\epsilon^{3}_{t}(\epsilon^{2}_{t}+L^{2})\leq\epsilon^{3/2}_{t}\epsilon_{t}$
for sufficiently large $t$ (since $\epsilon_{t}\to 0$, see Assumption 6).
Thus,

$\displaystyle\|$
$\displaystyle\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}\leq(1+\theta_{t})(1+\kappa_{t})(1-\epsilon^{3/2}_{t}\epsilon_{t})\times$
$\displaystyle\quad\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}+(1+1/{\kappa_{t}})M_{\boldsymbol{y}}^{2}\epsilon^{3}_{t}(\epsilon_{t-1}-\epsilon_{t})^{2}$
$\displaystyle+(1+1/{\theta_{t}})O((r_{t}-r_{t-1})^{2}),$

where $M_{\boldsymbol{y}}^{2}$ is the uniform upper bound of the norm of the
sequence $\boldsymbol{y}(t)$. By rearranging the terms above, we obtain

$\displaystyle(1-(1+\theta_{t})$
$\displaystyle(1+\kappa_{t})(1-\epsilon^{3/2}_{t}\epsilon_{t}))\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle\leq(1+1/{\kappa_{t}})M_{\boldsymbol{y}}^{2}\epsilon^{3}_{t}(\epsilon_{t-1}-\epsilon_{t})^{2}$
$\displaystyle+(1+1/{\theta_{t}})O((r_{t}-r_{t-1})^{2}).$

We conclude the proof by noticing that, according to the choice of
$\epsilon^{3/2}_{t}$ and $\epsilon_{t}$ and by taking
$\kappa_{t}=\theta_{t}=\frac{1}{4}\epsilon^{3/2}_{t}\epsilon_{t}$, we obtain,

$\displaystyle(1-(1+\theta_{t})(1+\kappa_{t})(1-\epsilon^{3/2}_{t}\epsilon_{t}))$
$\displaystyle\geq\epsilon^{3/2}_{t}\epsilon_{t}-\theta_{t}-\kappa_{t}-\theta_{t}\kappa_{t}$
$\displaystyle\geq 0.4\epsilon^{3/2}_{t}\epsilon_{t}.$

∎

###  -E Proof of Proposition 1

###### Proof:

We demonstrate that $\|\mathbf{z}(t)-\boldsymbol{y}(t-1)\|\to 0$ as
$t\to\infty$. Then, the result follows from Theorem 3. From Lemma 1, since
$\boldsymbol{y}(t)=\mbox{Proj}_{(1-r_{t})\boldsymbol{A}}[\boldsymbol{y}(t)-\beta_{t}(\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}\boldsymbol{y}(t))]$,
we have

$\displaystyle\|$
$\displaystyle\mathbf{z}(t+1)-\boldsymbol{y}(t)\|^{2}\leq(1-2\epsilon_{t}\beta_{t})\|\mathbf{z}(t)-\boldsymbol{y}(t)\|^{2}$
(98)
$\displaystyle-2\beta_{t}(\mathbf{z}(t)-\boldsymbol{y}(t),\boldsymbol{M}(\mathbf{z}(t))-\boldsymbol{M}(\boldsymbol{y}(t)))$
(99)
$\displaystyle+\beta^{2}_{t}\|\boldsymbol{M}(\mathbf{z}(t))-\boldsymbol{M}(\boldsymbol{y}(t))+\epsilon_{t}(\mathbf{z}(t)-\boldsymbol{y}(t))\|^{2}$
(100) $\displaystyle\leq$
$\displaystyle(1-2\epsilon_{t}\beta_{t}+2L\beta^{2}_{t}+2\epsilon^{2}_{t}\beta^{2}_{t})\|\mathbf{z}(t)-\boldsymbol{y}(t)\|^{2}$
(101) $\displaystyle\leq$
$\displaystyle(1-\epsilon_{t}\beta_{t})\|\mathbf{z}(t)-\boldsymbol{y}(t)\|^{2},$
(102)

where the last inequality is due to the fact that
$\sum_{t=1}^{\infty}\epsilon_{t}\beta_{t}=\infty$,
$\sum_{t=1}^{\infty}\beta^{2}_{t}<\infty$ (see Assumption 6), and, thus, there
exists $t_{0}$ such that
$2L\beta^{2}_{t}+2\epsilon^{2}_{t}\beta^{2}_{t}\leq\epsilon_{t}\beta_{t}$ for
all $t>t_{0}$. Next, for sufficiently large $t$ we obtain

$\displaystyle\|$
$\displaystyle\mathbf{z}(t+1)-\boldsymbol{y}(t)\|^{2}\leq(1+\beta_{t}\epsilon_{t})(1-\epsilon_{t}\beta_{t})\|\mathbf{z}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle+(1+1/{\beta_{t}\epsilon_{t}})(1-\epsilon_{t}\beta_{t})\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle\leq$
$\displaystyle(1-0.5\epsilon_{t}\beta_{t})\|\mathbf{z}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle+(1+1/{\beta_{t}\epsilon_{t}})(1-\epsilon_{t}\beta_{t})\|\boldsymbol{y}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle\leq$
$\displaystyle(1-0.5\epsilon_{t}\beta_{t})\|\mathbf{z}(t)-\boldsymbol{y}(t-1)\|^{2}$
$\displaystyle+O\left(\frac{(\epsilon_{t}-\epsilon_{t-1})^{2}}{\beta_{t}\epsilon^{3}_{t}}+\frac{(r_{t}-r_{t-1})^{2}}{\beta_{t}\epsilon^{6}_{t}}\right),$
(103)

where in the last inequality we used Lemma 4. Next by using Assumption 6 and
the result in Theorem 5 (see Appendix -A) for the sequence
$\|\mathbf{z}(t)-\boldsymbol{y}(t-1)\|^{2}$, we conclude that
$\mathbf{z}(t)\to\boldsymbol{y}(t-1)$ as $t\to\infty$. Thus, according to
Theorem 3, the sequence $\mathbf{z}(t)$ converges to the least norm solution
of $VI(\boldsymbol{A},\boldsymbol{M})$. ∎

###  -F Verification of Equations (IV-C) and (81)

Due to Assumption 5 there exists a compact set
$\mathbb{S}\subset\mathbb{R}^{Nd}$ such that for any
$\boldsymbol{x}\notin\mathbb{S}$

$J_{i}(\boldsymbol{x})\leq(\boldsymbol{c},\boldsymbol{x})+b_{0}$

for some
$\boldsymbol{c}=(c_{1}^{1},\ldots,c_{d}^{1},\ldots,c_{1}^{N},\ldots,c_{d}^{N})\in\mathbb{R}^{Nd}$
and $b_{0}\in\mathbb{R}$. Thus, for some positive $S$, $d_{1}$ and $d_{2}$ we
get

$\displaystyle\int_{\mathbb{R}^{Nd}}{J_{i}}^{2}(\boldsymbol{x})\frac{(x^{i}_{k}-\mu^{i}_{k})^{2}}{\sigma^{4}_{t}}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\leq\int_{\mathbb{S}}{J_{i}}^{2}(\boldsymbol{x})\frac{(x^{i}_{k}-\mu^{i}_{k})^{2}}{\sigma^{4}_{t}}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\qquad+\int_{\mathbb{R}^{Nd}\setminus\mathbb{S}}[d_{1}\|x\|^{2}+d_{2}]\frac{(x^{i}_{k}-\mu^{i}_{k})^{2}}{\sigma^{4}_{t}}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\leq\frac{S}{\sigma^{4}_{t}}+\int_{\mathbb{R}^{Nd}}[d_{1}\|x\|^{2}+d_{2}]\frac{(x^{i}_{k}-\mu^{i}_{k})^{2}}{\sigma^{4}_{t}}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle=\frac{S}{\sigma^{4}_{t}}+\frac{f_{i}(\boldsymbol{\mu},\sigma_{t})}{\sigma^{4}_{t}},$

$f_{i}(\boldsymbol{\mu},\sigma_{t})$ being a quadratic function of
$\boldsymbol{\mu}$ and $\sigma_{t}$. The last equality is due to the fact that

$\displaystyle\int_{\mathbb{R}^{Nd}}(x^{i}_{k}-\mu^{i}_{k})^{2}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}=\sigma^{2}_{t},$
$\displaystyle\int_{\mathbb{R}^{Nd}}(x^{i}_{k})^{2}(x^{i}_{k}-\mu^{i}_{k})^{2}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\quad\leq\int_{\mathbb{R}^{Nd}}[2(x^{i}_{k}-\mu^{i}_{k})^{2}+(\mu^{i}_{k})^{2})](x^{i}_{k}-\mu^{i}_{k})^{2}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\quad=2\sigma^{4}_{t}+(\mu^{i}_{k})^{2})\sigma^{2}_{t},$
$\displaystyle\int_{\mathbb{R}^{Nd}}(x^{j}_{m})^{2}(x^{i}_{k}-\mu^{i}_{k})^{2}p(\boldsymbol{\mu},\boldsymbol{x})d\boldsymbol{x}$
$\displaystyle\quad=(\sigma^{2}_{t}+(\mu^{j}_{m})^{2})\sigma^{2}_{t}|.$

Analogously one obtains the estimation (81).