# An Optimal Distributed Algorithm with Operator Extrapolation for Stochastic
Aggregative Games

Tongyu Wang<EMAIL_ADDRESS>Peng Yi<EMAIL_ADDRESS>Jie Chen
<EMAIL_ADDRESS>[ Department of Control Science and Engineering,
Tongji University, Shanghai, 201804, China;  The Shanghai Research Institute
for Intelligent Autonomous Systems,Shanghai, 201804, China.  The Institute of
Advanced Study, Tongji University, Shanghai, 201804, China  Frontiers Science
Center for Intelligent Autonomous Systems, Tongji University, Shanghai, China
Shanghai Institute of Intelligent Science and Technology, Tongji University,
Shanghai 200092, China

###### Abstract

This work studies Nash equilibrium seeking for a class of stochastic
aggregative games, where each player has an expectation-valued objective
function depending on its local strategy and the aggregate of all players’
strategies. We propose a distributed algorithm with operator extrapolation, in
which each player maintains an estimate of this aggregate by exchanging this
information with its neighbors over a time-varying network, and updates its
decision through the mirror descent method. An operator extrapolation at the
search direction is applied such that the two step historical gradient samples
are utilized to accelerate the convergence. Under the strongly monotone
assumption on the pseudo-gradient mapping, we prove that the proposed
algorithm can achieve the optimal convergence rate of $\mathcal{O}(1/k)$ for
Nash equilibrium seeking of stochastic games. Finally, the algorithm
performance is demonstrated via numerical simulations.

###### keywords:

Aggregative games, Stochastic game, Distributed algorithm, Operator
extrapolation, Optimal convergence rate

t1t1footnotetext: This research is supported by Shanghai Municipal Science and
Technology Major Project No. 2021SHZDZX0100, Shanghai Municipal Commission of
Science and Technology Project No. 19511132101, National Natural Science
Foundation of China under Grant 62003239, Shanghai Sailing Program No.
20YF1453000, and the Fundamental Research Funds for the Central Universities,
China.

url]http://www.elsevier.com

##  1 Introduction

In recent years, the non-cooperative games have been widely applied in
decision-making for networked systems, where selfish networked players aim to
optimize their own objective functions that have coupling with other players’
decisions. Nash equilibrium (NE) is one of the most widely used solution
concept to non-cooperative games, under which no participant can improve its
utility by unilaterally deviating from the equilibrium strategy [1].
Specifically, the aggregative games is an important class of non-cooperative
games, in which each player’s cost function depends on its own strategy and
the aggregation of all player strategies [2]. This aggregating feature emerges
in numerous decision-making problems over networks systems, when the
individual utility is affected by network-wide aggregate. Hence, aggregative
games and its NE seeking algorithm get wide research, with application in
charging coordination for plug-in electric vehicles [3, 4], multidimensional
opinion dynamics[5, 6], communication channel congestion control [7], and
energy management in power grids [8, 9], etc. It also worth pointing out that
stochastic game should be taken as a proper decision-making model when there
are information uncertainties, dating back to [10], while the NE seeking for
stochastic games with various sampling schemes are recently studied [11, 12].
Hence, when the game has both information uncertainties and aggregating-type
objective functions, stochastic aggregative game and its NE seeking are
studied [13, 14, 15].

Generally speaking, NE seeking methods include centralized, distributed, and
semi-decentralized methods. In the centralized methods [16], each player needs
to know the complete information to compute an NE in an introspective manner.
In the semi-decentralized methods [17], there exists a central coordinator
gathering and broadcasting signals to all players. While in distributed
methods [18], each player has its private information and seeks the NE through
local communications with other players. As such, designing distributed NE
seeking algorithm has gained tremendous research interests due to the benefits
owning to decentralization.

Existing works on distributed NE seeking largely resided in best-response(BR)
schemes [19][12][6] and gradient-based approaches [20][21][22]. In BR schemes,
players select a strategy that best responds to its rivals’ current
strategies. For example, [23] proposed an inexact generalization of this
scheme in which an inexact best response strategy is computed via a stochastic
approximation scheme. In gradient-based approaches, the algorithm is easily
implementable with a low computation complexity per iteration. For instance,
[24] designed an accelerated distributed gradient play for a strongly monotone
game and proved its geometric convergence, while it needs to estimate all
players’ strategies. To reduce the communication costs, [25] proposed a
consensus-based distributed gradient algorithm for aggregative games, in which
players only need to estimate the aggregate in each iteration. Furthermore,
[26] proposed a distributed gradient algorithm based on iterative Tikhonov
regularization method to resolve the class of monotone aggregative games. [27]
proposed a fully asynchronous distributed algorithm and rigorously show the
convergence to a Nash equilibrium. Besides, [28, 9] proposed a continuous-time
distributed algorithm for aggregative games.

Literature |  stochastic |  Method |  rate |  Iteration |  communication | sample
---|---|---|---|---|---|---
[6] |  $\times$ |  BR-based |  - |  - |  $v$ rounds |  -
[27] |  $\times$ |  gradient-based |  $1/k$ |  - |  $1$ round |  -
[29] |  ✓ | gradient-based |  $1/k$ |  $\mathcal{O}(\ln(1/\epsilon))$ |  $k+1$ rounds |  $N_{k}$
This work |  ✓ |  gradient-based |  $1/k$ |  Corollary1 |  $1$ round |  1

Table 1: A list of some recent research papers on distributed methods for
resolving aggregative games.

In column 2 “stochastic”, $\checkmark$ means that the studied problem is
stochastic, while $\times$ implies that the studied problem is deterministic;
In columns 4 and 5, dash - implies that the literature has not studied this
property; In column 6 “communication” stands for that the communication round
required at each iteration; In column 7 “sample” stands for the number of
sampled gradients per iteration.

Moreover, fast convergence rate is indispensable bonus for distributed Nash
equilibrium seeking algorithm, since it implies less communication rounds.
Particulary, for stochastic games, distributed NE seeking with a fast
convergence rate is highly desirable since its also implies low sampling cost
of stochastic gradients. Motivated by this, [30] designed an extra-gradient
method to improve the speed of NE seeking, while it needs two steps of
operator evaluations and two step of projections at each iteration. To further
reduce the computation burden, [31] proposed a simpler recursion with one step
of operator evaluation and projection at each iteration, given by
$x_{t+1}=\operatorname{argmin}_{x\in X}\gamma_{t}\left\langle
F\left(x_{t}+\beta_{t}\left(x_{t}-x_{t-1}\right)\right),x\rangle+V\left(x_{t},x\right)\right.$,
called projected reflected gradient method. But the iterate
$x_{t}+\beta_{t}\left(x_{t}-x_{t-1}\right)$ may sit outside the feasible set
$X$, so it needs to impose the strong monotone assumption on the pseudo-
gradient over $\mathbb{R}^{n}$. Recently, [32] proposed an operator
extrapolation (OE) method to solve the stochastic variational inequality
problem, for which the optimal convergence rate can be achieved through one
operator evaluation and a single projection per iteration.

Motivated by the demands for fast NE seeking algorithm of stochastic
aggregative games and inspired by [32], we propose a distributed operator
extrapolation (OE) algorithm via mirror descent and dynamical averaging
tracking. At each stage, each player aligns its intermediate estimate by a
consensus step with its neighbors’ estimates of the aggregate, samples an
unbiased estimate of its payoff gradients, takes a small step via the OE
method, and then mirrors it back to the strategy set. The algorithm achieves
the optimal convergence rate $\mathcal{O}(1/k)$ for the class of stochastic
strongly monotone aggregative games. The numerical experiments on a Nash-
cournot competition problem demonstrate it advantages over some existing NE
seeking methods. In addition, we compare our algorithm with some other
distributed algorithms for aggregative games and summarize those works in
Table 1. Specially, we compare the iteration complexity and communication
rounds to achieve an $\epsilon-\mathrm{NE}$ when the iterate $x$ satisfies
$\mathbb{E}[\|x-$ $\left.x^{*}\|^{2}\right]\leq\epsilon$. The proposed
algorithm can achieve the optimal convergence rate with only one communication
round per iteration in comparison with multiple communication rounds required
by [6] and [29]. Though [27] can achieve convergence rate $1/k$ with fixed
step-sizes for a deterministic game, our work can achieve the optimal
convergence result with diminishing step for a stochastic game, which can
model various information uncertainties.

The rest of paper is organized as follows. In section 2, we give the
formulation of the stochastic aggregate games, and state the assumptions. In
section 3, we propose a distributed operator extrapolation algorithm and
provide convergence results for the class of strongly monotone games. The
proof of main results are given in Section 4. We show the empirical
performance of our algorithm through the Nash-cournot model in section 5, and
give concluding remarks in Section 6.

Notations: A vector $x$ is a column vector while $x^{T}$ denotes its
transpose. $\langle x,y\rangle=x^{T}y$ denotes the inner product of vectors
$x,y.$ $\|x\|$ denotes the Euclidean vector norm, i.e., $\|x\|=\sqrt{x^{T}x}$.
A nonnegative square matrix $A$ is called doubly stochastic if
$A\mathbf{1}=\mathbf{1}$ and $\mathbf{1}^{T}A=\mathbf{1}^{T}$, where
$\mathbf{1}$ denotes the vector with each entry equal 1.
$\mathbf{I}_{N}\in\mathbb{R}^{N\times N}$ denotes the identity matrix. Let
$\mathcal{G}=\\{\mathcal{N},\mathcal{E}\\}$ be a directed graph with
$\mathcal{N}=\\{1,\cdots,N\\}$ denoting the set of players and $\mathcal{E}$
denoting the set of directed edges between players, where
$(j,i)\in\mathcal{E}$ if player $i$ can obtain information from player $j$.
The graph $\mathcal{G}$ is called strongly connected if for any
$i,j\in\mathcal{N}$ there exists a directed path from $i$ to $j$, i.e., there
exists a sequence of edges $(i,i_{1}),(i_{1},i_{2}),\cdots,(i_{p-1},j)$ in the
digraph with distinct nodes $i_{m}\in\mathcal{N},~{}\forall m:1\leq m\leq
p-1$.

##  2 Problem Formulation

In this section, we formulate the stochastic aggregative games over networks
and introduce basic assumptions.

###  2.1 Problem Statement

We consider a set of $N$ non-cooperative players indexed by
$\mathcal{N}\triangleq\\{1,\cdots,N\\}$. Each player $i\in{\mathcal{N}}$
choose its strategy $x_{i}$ from a strategy set $X_{i}\in\mathbb{R}^{m}$.
Denote by
$\mathbf{x}\triangleq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(x_{1}^{T},\cdots,x_{N}^{T})^{T}}\in\mathbb{R}^{mN}$
and $\mathbf{x}_{-i}\triangleq\\{x_{j}\\}_{j\neq i}$ the strategy profile and
the rival strategies, respectively. In an aggregative game, each player $i$
aim to minimize its cost function $f_{i}(x_{i},\sigma(x))$, where
$\sigma(\mathbf{x})\triangleq\sum_{j=1}^{N}x_{j}$ is an aggregate of all
players’ strategies. Furthermore, given
$\sigma(\mathbf{x}_{-i})\triangleq\sum_{j=1,j\neq i}^{N}x_{j}$, the objective
of player $i$ is to minimize its parameterized stochastic local cost:

$\displaystyle\min_{x_{i}\in
X_{i}}f_{i}(x_{i},x_{i}+\sigma(\mathbf{x}_{-i}))\triangleq\mathbb{E}\left[\psi_{i}(x_{i},x_{i}+\sigma(\mathbf{x}_{-i});{\xi_{i}(\omega)})\right],$
(1)

where ${\xi_{i}}:{\Omega}\to{\mathbb{R}^{d_{i}}}$ is a random variable defined
on the probability space $(\Omega,\mathcal{F},\mathbb{P}_{i})$,
$\psi_{i}:\mathbb{R}^{m}\times\mathbb{R}^{m}\times\mathbb{R}^{d_{i}}\to\mathbb{R}$
is a scalar-valued function, and $\mathbb{E}[\cdot]$ denotes the expectation
with respect to the probability measure $\mathbb{P}_{i}$.

A Nash equilibrium (NE) of the aggregative game (1) is a tuple
$\mathbf{x}^{*}=\\{x_{i}^{*}\\}_{i=1}^{N}$ such that for each
$i\in\mathcal{N}$,

$\displaystyle f_{i}(x_{i}^{*},\sigma(\mathbf{x}^{*}))\leq
f_{i}(x_{i},x_{i}+\sigma(\mathbf{x}_{-i}^{*})),\quad\forall x_{i}\in X_{i}.$

We consider the scenario that each player $i\in\mathcal{N}$ knows the
structure of its private function $f_{i}$, but have no access to the aggregate
$\sigma(\mathbf{x}).$ Instead, each player may communicate with its neighbors
over a time varying graph $\mathcal{G}_{k}=\\{\mathcal{N},\mathcal{E}_{k}\\}$.
Define $W_{k}=[\omega_{ij,k}]_{i,j=1}^{N}$ as the adjacency matrix, where
$\omega_{ij,k}>0$ if and only if $(j,i)\in\mathcal{E}_{k}$, and
$\omega_{ij,k}=0$, otherwise. Denote by
$N_{i,k}\triangleq\\{j\in\mathcal{N}:(j,i)\in\mathcal{E}_{k}\\}$ the
neighboring set of player $i$ at time $k$.

###  2.2 Assumptions

We impose the following conditions on the time-varying communication graphs
$\mathcal{G}_{k}=\\{\mathcal{N},\mathcal{E}_{k}\\}$.

###### Assumption 1

(a) $W_{k}$ is doubly stochastic for any $k\geq 0$;
(b) There exists a constant $0<\eta<1$ such that

$\omega_{ij,k}\geq\eta,\quad\forall j\in\mathcal{N}_{i,k},~{}~{}\forall
i\in\mathcal{N},~{}\forall k\geq 0;$

(c) There exists a positive integer $B$ such that the union graph
$\big{\\{}\mathcal{N},\bigcup_{l=1}^{B}\mathcal{E}_{k+l}\big{\\}}$ is strongly
connected for all $k\geq 0$.

We define a transition matrix $\Phi(k,s)=W_{k}W_{k-1}\cdots W_{s}$ for any
$k\geq s\geq 0$ with $\Phi(k,k+1)=\mathbf{I}_{N}$, and state a result that
will be used in the sequel.

###### Lemma 1

[33, Proposition 1] Let Assumption 1 hold. Then there exist
$\theta=(1-\eta/(4N^{2}))^{-2}>0$ and $\beta=(1-\eta/(4N^{2}))^{1/B}$ such
that for any $k\geq s\geq 0,$

$\displaystyle\left|\left[\Phi(k,s)\right]_{ij}-1/N\right|\leq\theta\beta^{k-s},\quad\forall
i,j\in\mathcal{N}.$ (2)

We require the player-specific problem to be convex and continuously
differentiable.

###### Assumption 2

For each player $i\in\mathcal{N},$ (a) the strategy set $X_{i}$ is closed,
convex and compact; (b) the cost function $f_{i}(x_{i},\sigma(x))$ is convex
in $x_{i}\in X_{i}$ for every fixed $x_{-i}\in X_{-i}\triangleq\prod_{j\neq
i}X_{j}$; (c) $f_{i}(x_{i},\sigma)$ is continuously differentiable in
$(x_{i},\sigma)\in X_{i}\times\mathbb{R}^{m}$.

For any $x\in X\triangleq\prod_{i=1}^{N}X_{i}$ and $z\in\mathbb{R}^{m}$,
define

$\displaystyle
F_{i}(x_{i},z)\triangleq\big{(}\nabla_{x_{i}}f_{i}(\cdot,\sigma)+\nabla_{\sigma}f_{i}(x_{i},\cdot)\big{)}\mid_{\sigma=z},$
$\displaystyle\phi_{i}(\mathbf{x})\triangleq\nabla_{x_{i}}f_{i}(x_{i},\sigma(\mathbf{x}))=F_{i}(x_{i},\sigma(\mathbf{x})),\mbox{~{}and~{}}\phi(\mathbf{x})\triangleq\left(\phi_{i}(\mathbf{x})\right)_{i=1}^{N}.$
(3)

From (2.2) and Assumption 2(c), it follows that the pseudogradient $\phi(x)$
is continuous. Since each player-specific problem is convex, by [34,
Proposition 1.4.2], $\mathbf{x}^{*}$ is a NE of (1) if and only if
$\mathbf{x}^{*}$ is a solution to a variational inequality problem
VI$(X,\phi)$, i.e., finding $\mathbf{x}^{*}\in X$ such that

$\displaystyle(\tilde{\mathbf{x}}-\mathbf{x}^{*})^{\top}\phi(\mathbf{x}^{*})\,\geq\,0,\quad\forall\tilde{\mathbf{x}}\in
X\triangleq\prod_{j=1}^{N}X_{j}.$ (4)

Since $\phi$ is continuous, $X$ is convex and compact, the existence of NE
follows immediately by [34, Corollary 2.2.5].

For each $i\in\mathcal{N},$ define

$\displaystyle M_{i}\triangleq\max_{x_{i}\in
X_{i}}\|x_{i}\|,~{}M_{H}\triangleq\sum_{j=1}^{N}M_{i},$ (5)
$\displaystyle\tilde{C}\triangleq\theta
M_{H}+2\theta\sum_{j=1}^{N}M_{j}/(1-\beta).$ (6)

We impose the following Lipschitz continuous conditions.

###### Assumption 3

(a) $\phi(\mathbf{x})$ is $L$-Lipschitz continuous in $x$, i.e.,

$\|\phi(\mathbf{x})-\phi(\mathbf{y})\|\leq
L\|x-y\|,\quad\forall\mathbf{x},\mathbf{y}\in X.$

(b) For each $i\in\mathcal{N}$ and any $x_{i}\in X_{i}$, $F_{i}(x_{i},z)$ is
$L_{fi}$-Lipschitz continuous in $z$ over some compact set, particularly, for
all $z_{1},z_{2}\in{\mathbb{R}^{d}}$ with $\|z_{1}\|\leq N\tilde{C}+M_{H}$ and
$\|z_{2}\|\leq N\tilde{C}+M_{H}$:

$\|F_{i}(x_{i},z_{1})-F_{i}(x_{i},z_{2})\|\leq L_{fi}\|z_{1}-z_{2}\|.$

In addition, we require each $\psi_{i}(x_{i},\sigma;\xi_{i})$ to be
differentiable, and assume there exists a stochastic oracle that returns
unbiased gradient sample with bounded variance.

###### Assumption 4

For each player $i\in\mathcal{N}$ and any ${\xi_{i}\in\mathbb{R}^{d_{i}}}$,
(a) $\psi_{i}(x_{i},\sigma;\xi_{i})$ is differentiable in $x_{i}\in X_{i}$ and
$\sigma\in\mathbb{R}^{m}$;
(b) for any $x_{i}\in X_{i}$ and $z\in\mathbb{R}^{m}$,
$q_{i}(x_{i},z;\xi_{i})=(\nabla_{x_{i}}\psi_{i}(x_{i},\sigma;\xi_{i})+\nabla_{\sigma}\psi_{i}(x_{i},\sigma;\xi_{i}))\mid_{\sigma=z}$
satisfies
$\mathbb{E}\left[q_{i}(x_{i},z_{i};\xi_{i})|x_{i}\right]=F_{i}(x_{i},z_{i})$
and
$\mathbb{E}[\|q_{i}(x_{i},z_{i};\xi_{i})-F_{i}(x_{i},z_{i})\|^{2}|x_{i}]\leq\nu_{i}^{2}$
for some constant $\nu_{i}>0.$

##  3 Algorithm Design and Main Results

In this section, we design a distributed NE seeking algorithm and prove its
optimal convergence to the Nash equilibrium in the mean-squared sense.

###  3.1 Distributed Mirror Descent Algorithm with Operator Extrapolation

We assume throughout that the paper that the regularization function $h$ is
$1$-strongly convex, i.e.,

$h(y)\geq h(x)+\langle\nabla h(x),y-x\rangle+{1\over
2}\|x-y\|^{2},\quad\forall x,y\in\mathbb{R}^{m}.$

We define the Bregman divergence associated with $h$ as follows.

$D(x,y)\triangleq h(y)-(h(x)+\langle\nabla h(x),y-x)\rangle,\quad\forall
x,y\in\mathbb{R}^{m}.$ (7)

Recall that the operator extrapolation method for VI in [32] requires a simple
recursion at each iteration given by

$x_{t+1}=\operatorname{argmin}_{x\in X}\gamma_{t}\left\langle
F\left(x_{t}\right)+\lambda_{t}\left(F\left(x_{t}\right)-F\left(x_{t-1}\right)\right),x\right\rangle+D\left(x_{t},x\right),$

which only involves one operator evaluation $F\left(x_{t}\right)$ and one
prox-mapping over the set $X$.

Suppose that each player $i$ at stage $k=1,2,\cdots$ selects a strategy
$x_{i,k}\in X_{i}$ as an estimate of its equilibrium strategy, and holds an
estimate $v_{i,k}$ for the average aggregate. At stage $k+1$, player $i$
observes its neighbors’ past information $v_{j,k},j\in\mathcal{N}_{i,k}$ and
updates an intermediate estimate by the consensus step (8), then it computes
its partial gradient based on the sample observation, and updates its strategy
$x_{i,k+1}$ by a mirror descent scheme (9) with operator extrapolation [32] by
setting
$q_{i}(x_{i,0},N\hat{v}_{i,1};\xi_{i,0})=q_{i}(x_{i,1},N\hat{v}_{i,2};\xi_{i,1})$
without loss of generality. Finally, player $i$ updates the average aggregate
with the renewed strategy $x_{i,k+1}$ by the dynamic average tracking scheme
(10). The procedures are summarized in Algorithm 1.

Algorithm 1 Distributed Mirror Descent Method with Operator Extrapolation

Initialize: Set $k=1,$ $x_{i,1}\in X_{i}$ and $v_{i,1}=x_{i,1}$ for each
$i\in\mathcal{N}$.

Iterate until convergence
Consensus. Each player computes an intermediate estimate by

$\displaystyle\hat{v}_{i,k+1}=\sum_{j\in\mathcal{N}_{i,k}}w_{ij,k}v_{j,k}.$
(8)

Strategy Update. Each player $i\in\mathcal{N}$ updates its equilibrium
strategy and its estimate of the average aggregate by

$\displaystyle x_{i,k+1}$ $\displaystyle=\mathop{\rm argmin}_{x_{i}\in
X_{i}}\Big{(}\alpha_{k}\big{\langle}(1+\lambda_{k})q_{i}(x_{i,k},N\hat{v}_{i,k+1};\xi_{i,k})$
$\displaystyle-\lambda_{k}q_{i}(x_{i,k-1},N\hat{v}_{i,k};\xi_{i,k-1}),x_{i}\big{\rangle}+D(x_{i,k},x_{i})\Big{)},$
(9) $\displaystyle v_{i,k+1}$
$\displaystyle=\hat{v}_{i,k+1}+x_{i,k+1}-x_{i,k}$ (10)

where $\alpha_{k}>0,\lambda_{k}>0$, and $\xi_{i,k}$ denotes a random
realization of $\xi_{i}$ at time $k$.

Define the gradient noise

$\displaystyle\zeta_{i,k}\triangleq
q_{i}(x_{i,k},N\hat{v}_{i,k+1};\xi_{i,k})-F_{i}(x_{i,k},N\hat{v}_{i,k+1}),$
(11)

$\mathbf{x}_{k}=(x_{1,k}^{T},\cdots,x_{N,k}^{T})^{T},$ and
$\mathcal{F}_{k}\triangleq\\{\mathbf{x}_{1},\xi_{i,l},i\in\mathcal{N},l=1,2,\cdots,k-1\\}.$
Then with Algorithm 1, $x_{i,k}$ and $\hat{v}_{i,k+1}$ are adapted to
$\mathcal{F}_{k}.$ Define
$\bm{\zeta}_{k}=(\zeta_{1,k}^{T},\cdots,\zeta_{N,k}^{T})^{T}.$ From Assumption
4, it follows that for each $i\in\mathcal{N}:$

$\begin{split}&\mathbb{E}[\bm{\zeta}_{k}|\mathcal{F}_{k}]=0,~{}\mathbb{E}[\|\zeta_{i,k}\|^{2}|\mathcal{F}_{k}]\leq\nu_{i}^{2},{\rm~{}and~{}}\\\
&\mathbb{E}[\|\bm{\zeta}_{k}\|^{2}|\mathcal{F}_{k}]\leq\sum_{i=1}^{N}\nu_{i}^{2}\triangleq\nu^{2}.\end{split}$
(12)

###  3.2 Main Results

Define

$\displaystyle\widehat{D}$
$\displaystyle(\mathbf{x},\mathbf{y})=\sum_{i=1}^{N}D(x_{i},y_{i})~{}{\rm
for~{}any~{}}\mathbf{x},\mathbf{y}\in X,$ (13)

With the definition of the Bregman’s distance, we can replace the strong
monotonicity assumption by the following assumption. This assumption taken
from [32] includes
$\langle\phi(\mathbf{x}),\mathbf{x}-\mathbf{x}^{*}\rangle\geq\mu\|\mathbf{x}-\mathbf{x}^{*}\|,\forall\mathbf{x}\in
X$ as the special case when $h(x)=\|x\|^{2}/2$.

###### Assumption 5

There exists a constat $\mu>0$ such that

$\langle\phi(\mathbf{x}),\mathbf{x}-\mathbf{x}^{*}\rangle\geq
2\mu\widehat{D}(\mathbf{x},\mathbf{x}^{*}),\quad\forall\mathbf{x}\in X.$

With this condition, we now state a convergence property of Algorithm 1, for
which the proof can be found in Section 4.2.

###### Proposition 1

Consider Algorithm 1. Let Assumptions 1-5 hold. Assume, in addition, that
there exists a positive sequence $\\{\theta_{k}\\}_{k\geq 1}$ satisfying

$\displaystyle\theta_{k+1}\alpha_{k+1}\lambda_{k+1}=\theta_{k}\alpha_{k},$
(14) $\displaystyle\theta_{k-1}\geq
16\alpha_{k}^{2}\lambda_{k}^{2}L^{2}\theta_{k},$ (15)
$\displaystyle\theta_{k}\leq\theta_{k-1}(2\mu\alpha_{k-1}+1){\rm~{}and~{}}8L^{2}\alpha_{k}^{2}\leq
1.$ (16)

Then

$\begin{split}&\theta_{t}(2\mu\alpha_{t}+1/2)\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]+\sum_{k=1}^{t-1}{\theta_{k}\over
4}\mathbb{E}[\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})]\\\
&\leq\theta_{1}\mathbb{E}[\widehat{D}(\mathbf{x}_{1},\mathbf{x}^{*})]+8\nu^{2}\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}+2\theta_{t}\alpha_{t}^{2}\nu^{2}\\\
&+\theta_{1}\alpha_{1}\lambda_{1}\mathbb{E}[\varepsilon_{1}]+2\sum_{k=1}^{t}\theta_{k}\alpha_{k}\mathbb{E}[\varepsilon_{k+1}],\end{split}$
(17)

where

$\begin{split}&\varepsilon_{k}\triangleq
2\sum_{i=1}^{N}L_{f_{i}}M_{i}\big{(}\|N\hat{v}_{i,k+1}-\sigma(\mathbf{x}_{k})\|+\|N\hat{v}_{i,k}-\sigma(\mathbf{x}_{k-1})\|\big{)}.\end{split}$
(18)

###### Remark 1

Consider the special case where the digraph $\mathcal{G}_{k}$ is a complete
graph for each time $k$ with $W_{k}={\mathbf{1}_{N}\mathbf{1}_{N}^{T}\over
N}$. Then (17) becomes

$\begin{split}&\theta_{t}(2\mu\alpha_{t}+1/2)\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]+\sum_{k=1}^{t-1}{\theta_{k}\over
4}\mathbb{E}[\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})]\\\
&\leq\theta_{1}\mathbb{E}[\widehat{D}(\mathbf{x}_{1},\mathbf{x}^{*})]+8\nu^{2}\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}+2\theta_{t}\alpha_{t}^{2}\nu^{2}.\end{split}$

This recovers the bound of [32, Theorem 3.3].

In the following, we establish a bound on the consensus error
$\|\sigma(\mathbf{x}_{k})-N\hat{v}_{i,k+1}\|$ of the aggregate, for which the
proof can be found in Section 4.3.

###### Proposition 2

Consider Algorithm 1. Let Assumptions 1, 2, and 3 hold. Then

$\begin{split}&\|\sigma(\mathbf{x}_{k})-N\hat{v}_{i,k+1}\|\leq\theta
M_{H}N\beta^{k}\\\ &+\theta
N\sum_{s=1}^{k}\beta^{k-s}\alpha_{s-1}\sum_{i=1}^{N}\Big{(}(1+2\lambda_{s-1})C_{i}+\lambda_{s-1}\|\zeta_{i,s-2}\|\\\
&\qquad\qquad\qquad\qquad\quad+(1+\lambda_{s-1})\|\zeta_{i,s-1}\|\Big{)},\end{split}$
(19)

where the constants $\theta,\beta$ are defined in (2), and

$\displaystyle C_{i}\triangleq N\tilde{C}L_{fi}+\max_{\mathbf{x}\in
X}\|\phi_{i}(\mathbf{x})\|.$ (20)

By combining Proposition 1 and Proposition 2, we can show that the proposed
method can achieve the optimal convergence rate for solving the stochastic
smooth and strongly monotone aggregative games. The proof can be found in
Section 4.4.

###### Theorem 1

Consider Algorithm 1. Suppose Assumptions 1-5 hold. Set

$\displaystyle c_{0}={4L\over\mu},~{}\alpha_{k}={1\over\mu(k+c_{0}-1)},$
$\displaystyle\theta_{k}=(k+c_{0}+1)(k+c_{0}),~{}\lambda_{k}={\theta_{k-1}\alpha_{k-1}\over\theta_{k}\alpha_{k}}.$

Then the following hold with $c_{e}\triangleq
4N\tilde{C}\sum_{i=1}^{N}L_{f_{i}}M_{i}$.

$\begin{split}&\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]\leq{2(c_{0}+2)(c_{0}+1)\widehat{D}(\mathbf{x}_{1},\mathbf{x}^{*})\over(t+c_{0}+1)(t+c_{0})}\\\
&+{8(\nu^{2}+c_{1})\over\mu^{2}(t+c_{0}+1)(t+c_{0})}+{2c_{0}(c_{0}+1)c_{e}\over\mu(c_{0}-1)(t+c_{0}+1)(t+c_{0})}\\\
&+{8(c_{2}+4\nu^{2})t\over\mu^{2}(t+c_{0}+1)(t+c_{0})},\end{split}$ (21)

where

$\displaystyle c_{1}$ $\displaystyle\triangleq 2\theta
M_{H}N(1+\beta)\mu\beta\left({(c_{0}-1)\over
1-\beta}+{1\over(1-\beta)^{2}}\right)\sum_{i=1}^{N}L_{f_{i}}M_{i},$ (22)
$\displaystyle c_{2}$ $\displaystyle\triangleq{8\theta
N(3-\beta)\sum_{i=1}^{N}(C_{i}+\nu_{i})\over(1-\beta)^{2}}\sum_{i=1}^{N}L_{f_{i}}M_{i}.$
(23)

###### Corollary 1

The number of iterations (the same as communication rounds) required by
Algorithm 1 for obtaining an approximate Nash equilibrium $\bar{\mathbf{x}}\in
X$ satisfying
$\mathbb{E}[\widehat{D}(\bar{\mathbf{x}},\mathbf{x}^{*})]\leq\epsilon$ is
bounded by

$\max\left({L\sqrt{\widehat{D}(\mathbf{x}_{1},\mathbf{x}^{*})}\over\mu\sqrt{\epsilon}},{\sqrt{\nu^{2}+c_{1}}\over\mu\sqrt{\epsilon}},{\sqrt{L}\over\mu\sqrt{\epsilon}},{c_{2}+4\nu^{2}\over\mu^{2}\epsilon}\right).$

##  4 Proof of Main Results

###  4.1 Preliminary Results

We now state a property from [35, Proposition B.3], and a well-known technical
results regarding the optimality condition of (9) (see [36, Lemma 3.1]).

###### Lemma 2

Let $h$ be a smooth and 1-strongly convex regularizer. Then

$\displaystyle D(x_{i},y_{i})\geq\frac{1}{2}\|x_{i}-y_{i}\|^{2},\quad\forall
x_{i},y_{i}\in X_{i},$ (24)

and

$\displaystyle\alpha_{k}\big{\langle}(1+\lambda_{k})q_{i}(x_{i,k},N\hat{v}_{i,k+1};\xi_{i,k}),x_{i,k+1}-x_{i}\big{\rangle}$
$\displaystyle-\alpha_{k}\big{\langle}\lambda_{k}q_{i}(x_{i,k-1},N\hat{v}_{i,k};\xi_{i,k-1}),x_{i,k+1}-x_{i}\big{\rangle}+D(x_{i,k},x_{i,k+1})$
$\displaystyle\leq D(x_{i,k},x_{i})-D(x_{i,k+1},x_{i}),\quad\forall x_{i}\in
X_{i}.$ (25)

Then from (13) and (24) it follows that

$\displaystyle\widehat{D}(\mathbf{x},\mathbf{y})\geq{1\over
2}\|\mathbf{x}-\mathbf{y}\|^{2},\quad\forall\mathbf{x},\mathbf{y}\in X.$ (26)

Furthermore, we define

$\displaystyle\Delta D_{i,k}(x_{i})\triangleq
D(x_{i,k},x_{i})-D(x_{i,k+1},x_{i}),$ (27) $\displaystyle\Delta
q_{i,k}\triangleq
q_{i}(x_{i,k},N\hat{v}_{i,k+1};\xi_{i,k})-q_{i}(x_{i,k-1},N\hat{v}_{i,k};\xi_{i,k-1}),$
(28) $\displaystyle\Delta\widehat{D}_{k}(\mathbf{x})=\sum_{i=1}^{N}\Delta
D_{i,k}(x_{i})=\widehat{D}(\mathbf{x}_{k},\mathbf{x})-\widehat{D}(\mathbf{x}_{k+1},\mathbf{x}).$
(29)

We are now ready to show the convergence properties of the proposed method.

###### Lemma 3

Let $\\{\mathbf{x}_{k}\\}$ be generated by Algorithm 1. Let Assumptions 2 and
3 hold. Suppose, in addition, that (14) and (15) hold for some positive
sequence $\\{\theta_{k}\\}_{k\geq 1}$.

Then

$\begin{split}&\sum_{k=1}^{t}\theta_{k}\left(\alpha_{k}\sum_{i=1}^{N}\big{\langle}q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}+\right.\\\
&\left.\widehat{D}(\mathbf{x}_{k+1},\mathbf{x})+{1\over
4}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})\right)\\\
&-2\alpha_{t}^{2}\theta_{t}L^{2}\|\mathbf{x}_{t+1}-\mathbf{x}\|^{2}-\alpha_{t}\theta_{t}\langle\bm{\zeta}_{t+1}-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}\\\
&\leq\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x})+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\varepsilon_{k}\\\
&+2\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}\|\bm{\zeta}_{k}-\bm{\zeta}_{k-1}\|^{2}+\alpha_{t}\theta_{t}\varepsilon_{t+1}.\end{split}$
(30)

Proof. With the definitions (27) and (28), (25) becomes

$\begin{split}\Delta
D_{i,k}(x_{i})&\geq\alpha_{k}\big{\langle}q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}\\\
&-\alpha_{k}\big{\langle}\Delta q_{i,k+1},x_{i,k+1}-x_{i}\big{\rangle}\\\
&+\alpha_{k}\lambda_{k}\langle\Delta q_{i,k},x_{i,k}-x_{i}\big{\rangle}\\\
&+\alpha_{k}\lambda_{k}\langle\Delta
q_{i,k},x_{i,k+1}-x_{i,k}\big{\rangle}+D(x_{i,k},x_{i,k+1}).\end{split}$ (31)

By multiplying both sides of the above inequality by $\theta_{k}$, and summing
up from $i=1$ to $N$, and $k=1$ to $t$, we obtain that

$\displaystyle\sum_{k=1}^{t}\theta_{k}\Delta\widehat{D}_{k}(\mathbf{x})\geq\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\sum_{i=1}^{N}\langle\Delta
q_{i,k},x_{i,k}-x_{i}\big{\rangle}$
$\displaystyle+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\sum_{i=1}^{N}\big{\langle}q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}$
$\displaystyle-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\sum_{i=1}^{N}\big{\langle}\Delta
q_{i,k+1},x_{i,k+1}-x_{i}\big{\rangle}$ (32)
$\displaystyle+\sum_{k=1}^{t}\theta_{k}\sum_{i=1}^{N}D(x_{i,k},x_{i,k+1})+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\sum_{i=1}^{N}\langle\Delta
q_{i,k},x_{i,k+1}-x_{i,k}\big{\rangle}.$

Then by recalling that $\Delta q_{i,1}=0$, using (13) and (14), we obtain

$\displaystyle\sum_{k=1}^{t}\theta_{k}\Delta\widehat{D}_{k}(\mathbf{x})\geq-\theta_{t}\alpha_{t}\sum_{i=1}^{N}\big{\langle}\Delta
q_{i,t+1},x_{i,t+1}-x_{i}\big{\rangle}$
$\displaystyle+\underbrace{\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\sum_{i=1}^{N}\langle\Delta
q_{i,k},x_{i,k+1}-x_{i,k}\big{\rangle}}_{\triangleq Q_{t}}$ (33)
$\displaystyle+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\sum_{i=1}^{N}\big{\langle}q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}.$

By the definitions of (11) and (28), we have

$\begin{split}&\sum_{i=1}^{N}\langle\Delta
q_{i,k},x_{i,k+1}-x_{i,k}\big{\rangle}=\sum_{i=1}^{N}\langle\zeta_{i,k}-\zeta_{i,k-1},x_{i,k+1}-x_{i,k}\big{\rangle}\\\
&+\sum_{i=1}^{N}\langle
F_{i}(x_{i,k},N\hat{v}_{i,k+1})-F_{i}(x_{i,k-1},N\hat{v}_{i,k}),x_{i,k+1}-x_{i,k}\big{\rangle}\\\
&=\langle\bm{\zeta}_{k}-\bm{\zeta}_{k-1},\mathbf{x}_{k+1}-\mathbf{x}_{k}\big{\rangle}\\\
&+\sum_{i=1}^{N}\langle
F_{i}(x_{i,k},N\hat{v}_{i,k+1})-F_{i}(x_{i,k-1},N\hat{v}_{i,k}),x_{i,k+1}-x_{i,k}\big{\rangle}.\end{split}$
(34)

By using (5) and Assumption 3, we derive

$\begin{split}&\sum_{i=1}^{N}\langle
F_{i}(x_{i,k},N\hat{v}_{i,k+1})-F_{i}(x_{i,k-1},N\hat{v}_{i,k}),x_{i,k+1}-x_{i,k}\rangle\\\
&\overset{\eqref{def-phi}}{=}\sum_{i=1}^{N}\langle
F_{i}(x_{i,k},N\hat{v}_{i,k+1})-F_{i}(x_{i,k},\sigma(\mathbf{x}_{k})),x_{i,k+1}-x_{i,k}\rangle\\\
&+\sum_{i=1}^{N}\langle
F_{i}(x_{i,k-1},\sigma(\mathbf{x}_{k-1}))-F_{i}(x_{i,k-1},N\hat{v}_{i,k}),x_{i,k+1}-x_{i,k}\rangle\\\
&+\sum_{i=1}^{N}\langle\phi_{i}(\mathbf{x}_{k})-\phi_{i}(\mathbf{x}_{k-1}),x_{i,k+1}-x_{i,k}\rangle\\\
&\geq-\sum_{i=1}^{N}\|F_{i}(x_{i,k},N\hat{v}_{i,k+1})-F_{i}(x_{i,k},\sigma(\mathbf{x}_{k}))\|\|x_{i,k+1}-x_{i,k}\|\\\
&-\sum_{i=1}^{N}\|F_{i}(x_{i,k-1},\sigma(\mathbf{x}_{k-1}))-F_{i}(x_{i,k-1},N\hat{v}_{i,k})\|\|x_{i,k+1}-x_{i,k}\|\\\
&+\langle\phi(\mathbf{x}_{k})-\phi(\mathbf{x}_{k-1}),\mathbf{x}_{k+1}-\mathbf{x}_{k}\rangle\\\
&\geq-2\sum_{i=1}^{N}L_{f_{i}}M_{i}\big{(}\|N\hat{v}_{i,k+1}-\sigma(\mathbf{x}_{k})\|+\|N\hat{v}_{i,k}-\sigma(\mathbf{x}_{k-1})\|\big{)}\\\
&-L\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|.\end{split}$
(35)

This together with (34), and the definitions of $Q_{t}$ and $\varepsilon_{k}$
in (33) and (18) implies that

$\begin{split}Q_{t}&\geq\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\langle\bm{\zeta}_{k}-\bm{\zeta}_{k-1},\mathbf{x}_{k+1}-\mathbf{x}_{k}\big{\rangle}\\\
&-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}L\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\varepsilon_{k}.\end{split}$
(36)

Then by (26), the following holds with defining $\theta_{0}\triangleq 0.$

$\displaystyle Q_{t}$ $\displaystyle\geq{1\over
4}\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})$
$\displaystyle+\sum_{k=1}^{t}\left(\underbrace{\theta_{k}\alpha_{k}\lambda_{k}\langle\bm{\zeta}_{k}-\bm{\zeta}_{k-1},\mathbf{x}_{k+1}-\mathbf{x}_{k}\big{\rangle}+{\theta_{k}\over
8}\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|^{2}}_{\triangleq{\rm~{}Term}~{}1}\right)$
$\displaystyle+\sum_{k=1}^{t}\left(-\theta_{k}\alpha_{k}\lambda_{k}L\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|\right)$
$\displaystyle+\sum_{k=1}^{t}\left({\theta_{k}\over
8}\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|^{2}+{\theta_{k-1}\over
8}\|\mathbf{x}_{k-1}-\mathbf{x}_{k}\|^{2}\right)$
$\displaystyle+{\theta_{t}\over
8}\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|^{2}-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\varepsilon_{k}.$
(37)

We let

$\displaystyle{\rm Term}~{}2\triangleq$
$\displaystyle-\theta_{k}\alpha_{k}\lambda_{k}L\|\mathbf{x}_{k}-\mathbf{x}_{k-1}\|\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|\
$ $\displaystyle+{\theta_{k}\over
8}\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|^{2}+{\theta_{k-1}\over
8}\|\mathbf{x}_{k-1}-\mathbf{x}_{k}\|^{2}.$

By using (15) and the Young’s inequality $a^{2}+b^{2}\geq 2\|a\|\|b\|$, we
derive ${\rm Term}~{}2\geq 0$. Also, we obtain
${\rm~{}Term}~{}1\geq-2\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}\|\bm{\zeta}_{k}-\bm{\zeta}_{k-1}\|^{2}$.
Then by substituting the bounds of Term 1 and Term 2 into (4.1), we have

$\begin{split}Q_{t}\geq&{1\over
4}\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})-2\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}\|\bm{\zeta}_{k}-\bm{\zeta}_{k-1}\|^{2}\\\
&+{\theta_{t}\over
8}\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|^{2}-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\varepsilon_{k}.\end{split}$
(38)

This together with (33) implies that

$\begin{split}&\sum_{k=1}^{t}\theta_{k}\Delta\widehat{D}_{k}(\mathbf{x})\geq{1\over
4}\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})-2\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}\|\bm{\zeta}_{k}-\bm{\zeta}_{k-1}\|^{2}\\\
&+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\sum_{i=1}^{N}\big{\langle}q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}\\\
&-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\varepsilon_{k}+{\theta_{t}\over
8}\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|^{2}-\theta_{t}\alpha_{t}\sum_{i=1}^{N}\big{\langle}\Delta
q_{i,t+1},x_{i,t+1}-x_{i}\big{\rangle}.\end{split}$ (39)

Similarly to (34), we have that

$\begin{split}&-\sum_{i=1}^{N}\langle\Delta
q_{i,t+1},x_{i,t+1}-x_{i}\big{\rangle}=-\langle\bm{\zeta}_{t+1}-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}\\\
&-\sum_{i=1}^{N}\langle
F_{i}(x_{i,t+1},N\hat{v}_{i,t+2})-F_{i}(x_{i,t},N\hat{v}_{i,t+1}),x_{i,t+1}-x_{i}\big{\rangle}\end{split}$
(40)

Similarly to (35), by using Assumption 3, (2.2), and (5), we have that

$\begin{split}&-\sum_{i=1}^{N}\langle
F_{i}(x_{i,t+1},N\hat{v}_{i,t+2})-F_{i}(x_{i,t},N\hat{v}_{i,t+1}),x_{i,t+1}-x_{i}\rangle\\\
&\geq-\varepsilon_{t+1}-L\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|\|\mathbf{x}_{t+1}-\mathbf{x}\|.\end{split}$

Therefore,

$\begin{split}&{1\over
8}\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|^{2}-\alpha_{t}\sum_{i=1}^{N}\big{\langle}\Delta
q_{i,t+1},x_{i,t+1}-x_{i}\big{\rangle}\\\ &\geq{1\over
8}\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|^{2}-\alpha_{t}\langle\bm{\zeta}_{t+1}-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}\\\
&-\alpha_{t}\varepsilon_{t+1}-\alpha_{t}L\|\mathbf{x}_{t+1}-\mathbf{x}_{t}\|\|\mathbf{x}_{t+1}-\mathbf{x}\|\\\
&\geq-\alpha_{t}\langle\bm{\zeta}_{t+1}-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}-2\alpha_{t}^{2}L^{2}\|\mathbf{x}_{t+1}-\mathbf{x}\|^{2}-\alpha_{t}\varepsilon_{t+1},\end{split}$

where the last inequality follows by the Young’ s inequality. This
incorporating with (39) produces

$\begin{split}&\sum_{k=1}^{t}\theta_{k}\Delta\widehat{D}_{k}(\mathbf{x})\geq{1\over
4}\sum_{k=1}^{t}\theta_{k}\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})-2\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}\|\bm{\zeta}_{k}-\bm{\zeta}_{k-1}\|^{2}\\\
&+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\sum_{i=1}^{N}\big{\langle}q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}\\\
&-\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\varepsilon_{k}-\alpha_{t}\theta_{t}\langle\bm{\zeta}_{k+1}-\bm{\zeta}_{k},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}\\\
&-2\alpha_{t}^{2}\theta_{t}L^{2}\|\mathbf{x}_{t+1}-\mathbf{x}\|^{2}-\alpha_{t}\theta_{t}\varepsilon_{t+1}.\end{split}$

By recalling the definition (29) and rearranging the terms, we prove the
lemma. $\Box$

###  4.2 Proof of Proposition 1

Since $\mathbf{x}_{k}$ is adapted to $\mathcal{F}_{k},$ from (12)it follows
that

$\mathbb{E}[\langle\bm{\zeta}_{k},\mathbf{x}_{k}-\mathbf{x}\big{\rangle}]=\mathbb{E}\big{[}\mathbb{E}[\langle\bm{\zeta}_{k},\mathbf{x}_{k}-\mathbf{x}\big{\rangle}|\mathcal{F}_{k}]\big{]}=\mathbb{E}\big{[}\mathbb{E}[\langle\bm{\zeta}_{k}|\mathcal{F}_{k}],\mathbf{x}_{k}-\mathbf{x}\big{\rangle}\big{]}=0$

for any $k\geq 1.$ Therefore,

$\displaystyle\mathbb{E}[\langle\bm{\zeta}_{t+1}-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}]=\mathbb{E}[\langle-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}\big{\rangle}]$
$\displaystyle=\mathbb{E}[\langle-\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}_{t}\big{\rangle}].$
(41)

By noting that
$q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1})=F_{i}(x_{i,k+1},N\hat{v}_{i,k+2})+\zeta_{i,k+1},$
we have

$\displaystyle\mathbb{E}[\langle
q_{i}(x_{i,k+1},N\hat{v}_{i,k+2};\xi_{i,k+1}),x_{i,k+1}-x_{i}\big{\rangle}]$
$\displaystyle=\mathbb{E}[\langle
F_{i}(x_{i,k+1},N\hat{v}_{i,k+2}),x_{i,k+1}-x_{i}\big{\rangle}].$ (42)

In addition, by (12) we have

$\mathbb{E}\big{[}\|\bm{\zeta}_{k}-\bm{\zeta}_{k-1}\|^{2}\big{]}\leq
2\mathbb{E}\big{[}\|\bm{\zeta}_{k}\|^{2}\big{]}+2\mathbb{E}\big{[}\|\bm{\zeta}_{k-1}\|^{2}\big{]}\leq
4\nu^{2}.$

Then by taking unconditional expectations on both sides of (30), using (4.2)
and (4.2), we obtain that

$\begin{split}&\sum_{k=1}^{t}\theta_{k}\left(\alpha_{k}\sum_{i=1}^{N}\mathbb{E}[\langle
F_{i}(x_{i,k+1},N\hat{v}_{i,k+2}),x_{i,k+1}-x_{i}\big{\rangle}]\right)\\\
&+\sum_{k=1}^{t}\theta_{k}\left(\mathbb{E}[\widehat{D}(\mathbf{x}_{k+1},\mathbf{x})]+{1\over
4}\mathbb{E}[\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})]\right)\\\
&-2\alpha_{t}^{2}\theta_{t}L^{2}\mathbb{E}[\|\mathbf{x}_{t+1}-\mathbf{x}\|^{2}]+\alpha_{t}\theta_{t}\mathbb{E}[\langle\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}_{t}\big{\rangle}]\\\
&\leq\sum_{k=1}^{t}\theta_{k}\mathbb{E}[\widehat{D}(\mathbf{x}_{k},x)]+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\mathbb{E}[\varepsilon_{k}]\\\
&+8\nu^{2}\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}+\alpha_{t}\theta_{t}\mathbb{E}[\varepsilon_{t+1}].\end{split}$
(43)

Note by Assumption 3(b), (2.2), and (5) that

$\begin{split}&\langle F_{i}(x_{i,k},N\hat{v}_{i,k+1}),x_{i,k}-x_{i}\rangle\\\
&=\langle F_{i}(x_{i,k},\sigma(\mathbf{x}_{k})),x_{i,k}-x_{i}\rangle\\\
&+\langle
F_{i}(x_{i,k},N\hat{v}_{i,k+1})-F_{i}(x_{i,k},\sigma(\mathbf{x}_{k})),x_{i,k}-x_{i}\rangle\\\
&\geq\langle\phi_{i}(\mathbf{x}_{k}),x_{i,k}-x_{i}\rangle\\\
&-\|F_{i}(x_{i,k},\sigma(\mathbf{x}_{k}))-F_{i}(x_{i,k},N\hat{v}_{i,k+1})\|\|x_{i,k}-x_{i}\|\\\
&\geq\langle\phi_{i}(\mathbf{x}_{k}),x_{i,k}-x_{i}\rangle-2L_{fi}M_{i}\|\sigma(\mathbf{x}_{k})-N\hat{v}_{i,k+1})\|.\end{split}$

Then by recalling the definition of $\varepsilon_{k}$ in (18), and using
Assumption 5, we obtain that

$\begin{split}\sum_{i=1}^{N}\langle
F_{i}(x_{i,k},N\hat{v}_{i,k+1}),x_{i,k}-x_{i}^{*}\rangle&\geq\langle\phi(\mathbf{x}_{k}),\mathbf{x}_{k}-\mathbf{x}^{*}\rangle-2\varepsilon_{k}\\\
&\geq
2\mu\widehat{D}(\mathbf{x}_{k},\mathbf{x}^{*})-\varepsilon_{k}.\end{split}$
(44)

Note by (26) and (12) that

$\begin{split}&{1\over
4}\mathbb{E}[\widehat{D}(\mathbf{x}_{t},\mathbf{x}_{t+1})]+\alpha_{t}\mathbb{E}[\langle\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}_{k}\big{\rangle}]\\\
&\geq{1\over
8}\mathbb{E}[\|\mathbf{x}_{t}-\mathbf{x}_{t+1}\|^{2}]+\alpha_{t}\mathbb{E}[\langle\bm{\zeta}_{t},\mathbf{x}_{t+1}-\mathbf{x}_{t}\big{\rangle}]\\\
&\geq{1\over
8}\mathbb{E}[\|\mathbf{x}_{t}-\mathbf{x}_{t+1}\|^{2}]-\alpha_{t}\sqrt{\mathbb{E}[\|\mathbf{x}_{t}-\mathbf{x}_{t+1}\|^{2}]}\sqrt{\mathbb{E}[\|\bm{\zeta}_{t}\|^{2}]}\\\
&\geq-2\alpha_{t}^{2}\mathbb{E}[\|\bm{\zeta}_{t}\|^{2}]\geq-2\alpha_{t}^{2}\nu^{2}.\end{split}$

This together with (43) and (44) implies that

$\begin{split}&\sum_{k=1}^{t}\theta_{k}(2\mu\alpha_{k}+1)\mathbb{E}[\widehat{D}(\mathbf{x}_{k+1},\mathbf{x}^{*})]\\\
&+\sum_{k=1}^{t-1}{\theta_{k}\over
4}\mathbb{E}[\widehat{D}(\mathbf{x}_{k},\mathbf{x}_{k+1})]-2\alpha_{t}^{2}\theta_{t}L^{2}\mathbb{E}[\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}]\\\
&\leq\sum_{k=1}^{t}\theta_{k}\mathbb{E}[\widehat{D}(\mathbf{x}_{k},\mathbf{x}^{*})]+\sum_{k=1}^{t}\theta_{k}\alpha_{k}\lambda_{k}\mathbb{E}[\varepsilon_{k}]+2\theta_{t}\alpha_{t}^{2}\nu^{2}\\\
&+8\nu^{2}\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}+\alpha_{t}\theta_{t}\mathbb{E}[\varepsilon_{t+1}]+{\sum_{k=1}^{t}\theta_{k}\alpha_{k}\mathbb{E}[\varepsilon_{k+1}]}.\end{split}$
(45)

By recalling $\theta_{k}\leq\theta_{k-1}(2\mu\alpha_{k-1}+1)$ from the
condition (16), we obtain that

$\begin{split}&\sum_{k=1}^{t}\theta_{k}(2\mu\alpha_{k}+1)\mathbb{E}[\widehat{D}(\mathbf{x}_{k+1},\mathbf{x}^{*})]-2\alpha_{t}^{2}\theta_{t}L^{2}\mathbb{E}[\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}]\\\
&\geq\theta_{t}\left((2\mu\alpha_{t}+1/2)\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]\right)+\sum_{k=1}^{t-1}\theta_{k+1}\mathbb{E}[\widehat{D}(\mathbf{x}_{k+1},\mathbf{x}^{*})]\\\
&+{\theta_{t}\over
2}\left(\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]-4\alpha_{t}^{2}L^{2}\mathbb{E}[\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}]\right)\\\
&\geq\theta_{t}(2\mu\alpha_{t}+1/2)\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]+\sum_{k=1}^{t-1}\theta_{k+1}\mathbb{E}[\widehat{D}(\mathbf{x}_{k+1},\mathbf{x}^{*})],\end{split}$
(46)

where the last inequality follows by (26) and $8L^{2}\alpha_{k}^{2}\leq 1$ of
(16). This together with (45) and (14) proves the result. $\Box$

###  4.3  Proof of Proposition 2.

Since $v_{i,1}=x_{i,1}$ and $W(k)$ is doubly stochastic, similarly to [25,
Lemma 2], we can show by induction that

$\displaystyle\sum_{i=1}^{N}v_{i,k}=\sum_{i=1}^{N}x_{i,k}=\sigma(\mathbf{x}_{k}),\quad\forall
k\geq 0.$ (47)

Akin to [25, Eqn. (16)], we obtain the following bound.

$\displaystyle\left\|{\sigma(\mathbf{x}_{k})\over
N}-\hat{v}_{i,k+1}\right\|\leq$ $\displaystyle\sum_{j=1}^{n}\left|{1\over
N}-[\Phi(k,1)]_{ij}\right|\|v_{j,1}\|$
$\displaystyle+\sum_{s=1}^{k}\sum_{j=1}^{N}\left|{1\over
N}-[\Phi(k,s)]_{ij}\right|\big{\|}x_{j,s}-x_{j,s-1}\big{\|}.$

Then by using (2) and $v_{i,1}=x_{i,1}$, we obtain that

$\begin{split}\left\|{\sigma(\mathbf{x}_{k})\over
N}-\hat{v}_{i,k+1}\right\|\leq\theta\beta^{k}\sum_{j=1}^{N}\|x_{j,1}\|+\theta\sum_{s=1}^{k}\beta^{k-s}\sum_{j=1}^{N}\big{\|}x_{j,s}-x_{j,s-1}\big{\|}.\end{split}$
(48)

This combined with (5) proves

$\displaystyle\left\|{\sigma(\mathbf{x}_{k})\over N}-\hat{v}_{i,k+1}\right\|$
$\displaystyle\leq\theta\beta^{k}M_{H}+2\theta\sum_{s=1}^{k}\beta^{k-s}\sum_{j=1}^{N}M_{j}$
$\displaystyle\leq\theta
M_{H}+2\theta\sum_{j=1}^{N}M_{j}/(1-\beta)\overset{\eqref{def-
tC}}{=}\tilde{C}.$ (49)

By (5) and (47), we have that $\|\sigma(\mathbf{x}_{k})\|\leq M_{H}$ for any
$k\geq 0$. Thus by (4.3), we obtain that for each $i\in\mathcal{N}$,

$\|N\hat{v}_{i,k+1}\|\leq N\tilde{C}+M_{H},\quad\forall k\geq 0.$

Then by using Assumption 3(b) and (2.2), we obtain that for each
$j\in\mathcal{N}$ and any $s\geq 0:$

$\displaystyle\|F_{j}(x_{j,s},N\hat{v}_{j,s})\|$
$\displaystyle\leq\|F_{j}(x_{j,s},N\hat{v}_{j,s})-F_{j}(x_{j,s},\sigma(\mathbf{x}_{s})\|$
$\displaystyle+\|F_{j}(x_{j,s},\sigma(\mathbf{x}_{s})\|$ $\displaystyle\leq
L_{fj}\|N\hat{v}_{j,s}-\sigma(\mathbf{x}_{s})\|+\|\phi_{j}(\mathbf{x}_{s})\|$
$\displaystyle\overset{\eqref{bd-con-
err}}{\leq}N\tilde{C}L_{fj}+\max_{\mathbf{x}\in
X}\|\phi_{j}(\mathbf{x})\|\overset{\eqref{def-C}}{=}C_{j}.$ (50)

By applying the optimality condition to (9), using the definitions (7) and
(11), we have that

$\begin{split}&(x_{i}-x_{i,k+1})^{T}\alpha_{k}\big{(}(1+\lambda_{k})(\zeta_{i,k}+F_{i}(x_{i,k},N\hat{v}_{i,k+1})\\\
&-\lambda_{k}(\zeta_{i,k-1}+F_{i}(x_{i,k-1},N\hat{v}_{i,k}))\big{)}\\\
&+\nabla h(x_{i,k+1})-\nabla h(x_{i,k})\big{)}\geq 0,\quad\forall x_{i}\in
X_{i}.\end{split}$ (51)

By setting $x_{i}=x_{i,k}$ in (51), rearranging the terms, and using the
assumption that $h$ is 1-strongly convex, we obtain

$\begin{array}[]{l}\alpha_{k}(x_{i,k}-x_{i,k+1})^{T}\big{(}(1+\lambda_{k})(\zeta_{i,k}+F_{i}(x_{i,k},N\hat{v}_{i,k+1})\\\
-\lambda_{k}(\zeta_{i,k-1}+F_{i}(x_{i,k-1},N\hat{v}_{i,k}))\big{)}\\\
\geq(x_{i,k}-x_{i,k+1})^{T}\big{(}\nabla h(x_{i,k})-\nabla
h(x_{i,k+1})\big{)}\geq\|x_{i,k}-x_{i,k+1}\|^{2}.\end{array}$

Then from (4.3) it follows that

$\displaystyle\|x_{i,k}-x_{i,k+1}\|$
$\displaystyle\leq\alpha_{k}\big{(}(1+\lambda_{k})(\|\zeta_{i,k}\|+\|F_{i}(x_{i,k},N\hat{v}_{i,k+1}\|)+$
$\displaystyle\lambda_{k}(\|\zeta_{i,k-1}\|+\|F_{i}(x_{i,k-1},N\hat{v}_{i,k})\|)\big{)}$
$\displaystyle{\leq}\alpha_{k}\big{(}(1+2\lambda_{k})C_{i}+(1+\lambda_{k})\|\zeta_{i,k}\|+\lambda_{k}\|\zeta_{i,k-1}\|\big{)}.$

This together with (5) and (48) produces (19). $\Box$

###  4.4 Proof of Theorem 1

By substituting (4.3) into the definition of $\varepsilon_{k}$ in (18), we
obtain that

$\begin{split}\varepsilon_{k}&=2\sum_{i=1}^{N}L_{f_{i}}M_{i}\big{(}\|N\hat{v}_{i,k+1}-\sigma(\mathbf{x}_{k})\|+\|N\hat{v}_{i,k}-\sigma(\mathbf{x}_{k-1})\|\big{)}\\\
&\leq 4N\tilde{C}\sum_{i=1}^{N}L_{f_{i}}M_{i}=c_{e},\quad\forall k\geq
1.\end{split}$ (52)

With the selection of parameters, similar to the proof of [32, Corollary 3.4],
we can verify that (14), (15), and (16) hold. Furthermore, by $c_{0}\geq 4$
and the simple calculations we obtain

$\displaystyle\theta_{t}(2\mu\alpha_{t}+1/2)>\theta_{t}/2=(t+c_{0}+1)(t+c_{0})/2,$
$\displaystyle\theta_{t}\alpha_{t}^{2}={(t+c_{0}+1)(t+c_{0})\over\mu^{2}(t+c_{0}-1)^{2}}\leq{2\over\mu^{2}},$
$\displaystyle\theta_{k}\alpha_{k}\leq{2\over\mu^{2}\alpha_{k}},~{}\theta_{k}\alpha_{k}\lambda_{k}=\theta_{k-1}\alpha_{k-1}\leq{2\over\mu^{2}\alpha_{k-1}}\leq{2\over\mu^{2}\alpha_{k}},$
$\displaystyle\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}\lambda_{k}^{2}=\sum_{k=1}^{t}\theta_{k}\alpha_{k}^{2}{\theta_{k-1}^{2}\alpha_{k-1}^{2}\over\theta_{k}^{2}\alpha_{k}^{2}}=\sum_{k=1}^{t}{\theta_{k-1}^{2}\alpha_{k-1}^{2}\over\theta_{k}}$
$\displaystyle=\sum_{k=1}^{t}{(k+c_{0})(k+c_{0}-1)^{2}\over\mu^{2}(k+c_{0}-2)^{2}(k+c_{0}+1)}\leq{2t\over\mu^{2}}.$

This incorporating with (17) and (52) implies that

$\begin{split}&{(t+c_{0}+1)(t+c_{0})\over
2}\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]\\\
&\leq\theta_{1}\mathbb{E}[\widehat{D}(\mathbf{x}_{1},\mathbf{x}^{*})]+{(16t+4)\nu^{2}\over\mu^{2}}\\\
&+{\theta_{1}\alpha_{1}\lambda_{1}c_{e}+{4\over\mu^{2}}\sum_{k=1}^{t}{\mathbb{E}[\varepsilon_{k+1}]\over\alpha_{k}}}.\end{split}$
(53)

In the following, we establish an upper bound of
$\sum_{k=1}^{t}{\mathbb{E}[\varepsilon_{k+1}]\over\alpha_{k}}.$ Note that

$\begin{split}&\sum_{k=1}^{t}k\beta^{k}=\beta\left(\sum_{k=1}^{t}k\beta^{k-1}\right)=\beta\left(\sum_{k=1}^{t}\beta^{k}\right)^{\prime}=\beta\left({\beta-\beta^{t+1}\over
1-\beta}\right)^{\prime}\\\
&={\beta\left(1-\beta^{t+1}-(t+1)\beta^{t}(1-\beta)\right)\over(1-\beta)^{2}}\leq{\beta\over(1-\beta)^{2}}.\end{split}$
(54)

Hence for any $t\geq 1,$

$\displaystyle\sum_{k=1}^{t}{\beta^{k}\over\alpha_{k}}=\sum_{k=1}^{t}\mu(k+c_{0}-1)\beta^{k}=\mu(c_{0}-1)\sum_{k=1}^{t}\beta^{k}+\mu\sum_{k=1}^{t}k\beta^{k}$
$\displaystyle\leq{\mu\beta(c_{0}-1)\over
1-\beta}+{\mu\beta\over(1-\beta)^{2}}.$ (55)

Also, note by $c_{0}\geq 4$ and (54) that

$\begin{split}&\sum_{k=1}^{t}{1\over\alpha_{k}}\sum_{s=1}^{k}\beta^{k-s}\alpha_{s-1}=\sum_{k=1}^{t}\sum_{s=1}^{k}\beta^{k-s}{\alpha_{s-1}\over\alpha_{k}}=\sum_{k=1}^{t}\sum_{s=1}^{k}\beta^{k-s}{k+c_{0}-1\over
s+c_{0}-2}\\\ &=\sum_{k=1}^{t}\sum_{s=1}^{k}\beta^{k-s}\left(1+{k-s+1\over
s+c_{0}-2}\right)\leq\sum_{k=1}^{t}\sum_{s=1}^{k}\beta^{k-s}\left(k-s+2\right)\\\
&=2\sum_{k=1}^{t}\sum_{s=0}^{k-1}\beta^{s}+\sum_{k=1}^{t}\sum_{s=0}^{k-1}s\beta^{s}\\\
&\leq\sum_{k=1}^{t}\left({2\over
1-\beta}+{\beta\over(1-\beta)^{2}}\right)={(2-\beta)t\over(1-\beta)^{2}}.\end{split}$
(56)

Similarly to (56), we have

$\begin{split}&\sum_{k=1}^{t}{1\over\alpha_{k}}\sum_{s=1}^{k+1}\beta^{k+1-s}\alpha_{s-1}=\sum_{k=1}^{t}\sum_{s=1}^{k+1}\beta^{k+1-s}\left(1+{k-s+1\over
s+c_{0}-2}\right)\\\
&\leq\sum_{k=1}^{t}\sum_{s=1}^{k+1}\beta^{k-s+1}\big{(}(k-s+1)+1\big{)}\\\
&\leq\sum_{k=1}^{t}\left(\sum_{s=0}^{k}s\beta^{s}+\sum_{s=0}^{k}\beta^{s}\right)\\\
&\leq\sum_{k=1}^{t}\left({1\over
1-\beta}+{\beta\over(1-\beta)^{2}}\right)={t\over(1-\beta)^{2}}.\end{split}$
(57)

Note by $c_{0}\geq 4$ that for any $\geq 1,$

$\displaystyle\lambda_{k}={\theta_{k-1}\alpha_{k-1}\over\theta_{k}\alpha_{k}}={(k+c_{0}-1)^{2}\over(k+c_{0}+1)(k+c_{0}-2)}\leq
3/2.$

By (12), we obtain that
$\mathbb{E}[\|\zeta_{i,k}\|]\leq\sqrt{\mathbb{E}[\|\zeta_{i,k}\|^{2}]}\leq\nu_{i},~{}\forall
k\geq 1.$ This together with (19) produces

$\begin{split}&\mathbb{E}[\|\sigma(\mathbf{x}_{k})-N\hat{v}_{i,k+1}\|]\leq\theta
M_{H}N\beta^{k}\\\ &+\theta
N\sum_{s=1}^{k}\beta^{k-s}\alpha_{s-1}\sum_{i=1}^{N}\big{(}(1+2\lambda_{s-1})C_{i}+(1+2\lambda_{s-1})\nu_{i}\big{)}\\\
&\leq\theta M_{H}N\beta^{k}+4\theta
N\sum_{i=1}^{N}\big{(}C_{i}+\nu_{i}\big{)}\sum_{s=1}^{k}\beta^{k-s}\alpha_{s-1}.\end{split}$

Then by using (4.4) and (56), we obtain that

$\begin{split}&\sum_{k=1}^{t}{1\over\alpha_{k}}\mathbb{E}[\|\sigma(\mathbf{x}_{k})-N\hat{v}_{i,k+1}\|]\\\
&\leq\theta M_{H}N\mu\beta\left({c_{0}-1\over
1-\beta}+{1\over(1-\beta)^{2}}\right)+{4\theta
N\sum_{i=1}^{N}(C_{i}+\nu_{i})(2-\beta)t\over(1-\beta)^{2}}.\end{split}$ (58)

Similarly, by using (4.4) and (57), we obtain that

$\begin{split}&\sum_{k=1}^{t}{1\over\alpha_{k}}\mathbb{E}[\|\sigma(\mathbf{x}_{k+1})-N\hat{v}_{i,k+2}\|]\\\
&\leq\theta M_{H}N\mu\beta^{2}\left({c_{0}-1\over
1-\beta}+{1\over(1-\beta)^{2}}\right)+{4\theta
N\sum_{i=1}^{N}(C_{i}+\nu_{i})t\over(1-\beta)^{2}}.\end{split}$ (59)

By recalling the definition of $\varepsilon_{k}$ in (18), we obtain that

$\begin{split}&\sum_{k=1}^{t}{\mathbb{E}[\varepsilon_{k+1}]\over\alpha_{k}}=2\sum_{i=1}^{N}L_{f_{i}}M_{i}\sum_{k=1}^{t}{\mathbb{E}[\|N\hat{v}_{i,k+1}-\sigma(\mathbf{x}_{k})\|]\over\alpha_{k}}\\\
&+2\sum_{i=1}^{N}L_{f_{i}}M_{i}\sum_{k=1}^{t}{\mathbb{E}[\|N\hat{v}_{i,k+2}-\sigma(\mathbf{x}_{k+1})\|]\over\alpha_{k}}\\\
&\leq 2\theta M_{H}N\mu\beta(1+\beta)\left({c_{0}-1\over
1-\beta}+{1\over(1-\beta)^{2}}\right)\sum_{i=1}^{N}L_{f_{i}}M_{i}\\\
&+t{8\theta
N\sum_{i=1}^{N}(C_{i}+\nu_{i})(3-\beta)\over(1-\beta)^{2}}\sum_{i=1}^{N}L_{f_{i}}M_{i}=c_{1}+c_{2}t,\end{split}$

where the last equality follows from the definitions of $c_{1}$ and $c_{2}$ in
(22) and (23). This together with (53) and
$\alpha_{1}\theta_{1}\lambda_{1}=\alpha_{0}\theta_{0}={c_{0}(c_{0}+1)\over\mu(c_{0}-1)}$
implies that

$\begin{split}&{(t+c_{0}+1)(t+c_{0})\over
2}\mathbb{E}[\widehat{D}(\mathbf{x}_{t+1},\mathbf{x}^{*})]\\\
&\leq\theta_{1}\mathbb{E}[\widehat{D}(\mathbf{x}_{1},\mathbf{x}^{*})]+{(16t+4)\nu^{2}\over\mu^{2}}\\\
&+{c_{0}(c_{0}-1)C_{e}\over\mu(c_{0}-1)}+{4(c_{1}+c_{2}t)\over\mu^{2}}.\end{split}$
(60)

Then the result follows immediately. $\Box$

##  5 Numerical Simulations

In this section, we validate the algorithm performance through numerical
simulation on the Nash-Cournot games (see e.g., [21, 14]). There is a
collection of $N$ factories denoted by $\mathcal{N}=\\{1,\dots,N\\}$ competing
over $l$ markets denoted by $\mathcal{L}\triangleq\\{1,\cdots,L\\}$, where
each factory $i$ needs to decide its production $x_{i,l}$ at markets $l$.
Then, the cost of production of factory $i$ is defined as
$c_{i}\left(x_{i};\xi_{i}\right)=\left(c_{i}+\xi_{i}\right)^{T}x_{i}$, where
$x_{i}\triangleq\operatorname{col}(x_{i,j})_{l=1}^{L}$, $c_{i}>0$ is parameter
for factory $i$, and $\xi_{i}$ is a random disturbance or noise with zero-mean
and bounded variance. The income of factory $i$ is
$p_{l}\left(S_{l};\zeta_{l}\right)x_{i}$, where $S_{l}=\sum_{i=1}^{N}x_{i,l}$
denote the aggregate products of all factories delivered to market $l$. By the
law of supply and demand, the price function $p_{l}$ can be represented by the
reverse demand function and is defined as
$p_{l}\left(S_{l};\zeta_{l}\right)=d_{l}+\zeta_{l}-b_{l}S_{l}$, where
$d_{l}>0$, $b_{l}>0$, and $\zeta_{l}$ is zero-mean random disturbance or
noise. Then, the factory $i$’s cost is
$\psi_{i}\left(x;\xi_{i},\zeta_{i}\right)=c_{i}\left(x_{i};\xi_{i}\right)-(d+\zeta-
BS)^{T}x_{i}$, with
$d=\left(d_{1},\cdots,d_{l}\right)^{T},B=\operatorname{diag}\left\\{b_{1},\cdots,b_{l}\right\\}$,
$S=\operatorname{col}\left\\{S_{1},\dots,S_{l}\right\\}$ and
$\zeta=\operatorname{col}\left\\{\zeta_{1},\cdots,\zeta_{l}\right\\}.$
Finally, factory $i$’s local optimization problem is $\min_{x_{i}\in
X_{i}}f_{i}(x)\triangleq\mathbb{E}\left[\psi_{i}\left(x;\xi_{i},\zeta_{i}\right)\right],$
while satisfying a finite capacity constraint $X_{i}$.

It is straightforward to verify that the aforementioned Nash-cournot example
satisfy the aggregative game formulation (1) with
$\sigma(x)=\sum_{i=1}^{N}x_{i}$, and
$f_{i}(x,\sigma)=c_{i}^{T}x_{i}-(d-B\sigma)^{T}x_{i}.$ Then by (2.2),
$F_{i}\left(x_{i},z\right)=c_{i}-d+B\left(z+x_{i}\right)$ and $\phi_{i}(x)=$
$c_{i}-d+B\left(\sum_{i=1}^{N}x_{i}+x_{i}\right)$. We can verify that
Assumptions 2, 3 and 4 hold when the random variables
$\xi_{i},\zeta_{l},i\in\mathcal{N},l\in\mathcal{L}$ are zero mean with bounded
variance. In addition, we set $h(x)=\frac{1}{2}\|x\|_{2}^{2}$, then the
Bregman divergence becomes $D_{h}(x,y)=\|x-y\|_{2}^{2}$.

Set $N=20,l=3$, and let the communication among the factories be described by
an undirected time-varying graph. The graph at each iteration is randomly
drawn from a set of four graphs, whose union graph is connected. Set the
adjacency matrix $W=\left[w_{ij}\right]$, where
$w_{ij}=\frac{1}{\max\left\\{\left|\mathcal{N}_{i}\right|,\left|\mathcal{N}_{j}\right|\right\\}}$
for any $i\neq j$ with $(i,j)\in\mathcal{E},w_{ii}=1-\sum_{j\neq i}w_{ij}$,
and $w_{ij}=0$, otherwise. Let $c_{i}$ is drawn from the uniform distribution
$U[3,4]$. The pricing parameters $d_{l},b_{l}$ of market $l\in\mathcal{L}$ are
derived from uniform distributions $U[10,10.5]$ and $U[0.5,1]$, respectively.
The capacity constraint for each $x_{i,l}$ is the same as
$x_{i,l}\in[2,10],\forall i\in\mathcal{N},l\in\mathcal{L}$. After fixing the
problem parameters $c_{i}$ and $d_{l}$, let the random variables
$\xi_{i},i\in\mathcal{N}$ and $\zeta_{l},l\in\mathcal{L}$ be randomly and
uniformly from $U\left[-c_{i}/8,c_{i}/8\right]$ and
$U\left[-d_{l}/8,d_{l}/8\right]$, respectively.

We implement Algorithm 1 and display the empirical results in Fig.1 by
averaging of 20 sampling trajectories with the same initial points, with
$\alpha_{k}$ and $\lambda_{k}$ choosing as in Theorem 1. Besides, we compare
our algorithm with the projected gradient method [25] and extra-gradient
method [30] when applied to the stochastic aggregative game considered in this
work, but the network aggregate value is still estimated with the dynamical
average tracking. The projected gradient method (PGA) requires a simple
projection given by

$x_{i}^{k+1}:=\operatorname{P}_{X_{i}}\left[x_{i}^{k}-\alpha_{k}q_{i}(x_{i,k},N\hat{v}_{i,k+1};\xi_{i,k})\right],$

while the extra-gradient method (Extra-G) consists of two steps, requiring two
steps of projection evaluation and gradient sample at each iteration.

$\begin{array}[]{cl}\text{ (extrapolation)
}&x_{i}^{k+1/2}=\operatorname{P}_{{X_{i}}}\left[x_{i}^{k}-\alpha_{k}q_{i}(x_{i,k},N\hat{v}_{i,k+1};\xi_{i,k})\right],\\\
\text{ (update)
}&x_{i}^{k+1}=\operatorname{P}_{{X_{i}}}\left[x_{i}^{k}-\alpha_{k}q_{i}(x_{i}^{k+1/2},N\hat{v}_{i,k+1};\xi_{i,k})\right].\end{array}$

Figure 1: Comparisons of Algorithm 1 with the projected gradient method (PGA)
and the extra-gradient method (Extra-G). The trajectories are derived by
averaging with 20 sample paths.

Fig.1 displays the convergence of the three algorithms, and it shows the
superior convergence speed of Algorithm 1 compared to projected gradient
method. Though the convergence speed of Algorithm 1 is almost the same as that
of the extra-gradient method, the advantage of our method lies in that it only
requires a single projection and one sampled gradient, greatly reducing the
computing and sampling cost.

##  6 Conclusions

This paper proposes a distributed operator extrapolation method for stochastic
aggregative game based on mirror descent, and shows that the proposed method
can achieve the optimal convergence for the class of strongly monotone games.
In addition, empirical results demonstrate that our method indeed brings
speed-ups. It is of interest to explore the algorithm convergence for monotone
games, and extend the operator extrapolation method to the other classes of
network games in distributed and stochastic settings.

## References

  * [1] J. F. Nash, Jr., “Equilibrium points in $n$-person games,” _Proc. Nat. Acad. Sci. U. S. A._ , vol. 36, pp. 48–49, 1950.
  * [2] M. K. Jensen, “Aggregative games and best-reply potentials,” _Economic theory_ , vol. 43, no. 1, pp. 45–66, 2010.
  * [3] 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_ , vol. 61, no. 11, pp. 3315–3329, 2015.
  * [4] S. Grammatico, “Dynamic control of agents playing aggregative games with coupling constraints,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 9, pp. 4537–4548, 2017.
  * [5] J. Ghaderi and R. Srikant, “Opinion dynamics in social networks with stubborn agents: Equilibrium and convergence rate,” _Automatica_ , vol. 50, no. 12, pp. 3209–3215, 2014. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0005109814004154
  * [6] F. Parise, S. Grammatico, B. Gentile, and J. Lygeros, “Distributed convergence to nash equilibria in network and average aggregative games,” _Automatica_ , vol. 117, p. 108959, 2020.
  * [7] B. G. Bakhshayesh and H. Kebriaei, “Decentralized equilibrium seeking of joint routing and destination planning of electric vehicles: A constrained aggregative game approach,” _IEEE Transactions on Intelligent Transportation Systems_ , pp. 1–10, 2021.
  * [8] G. Belgioioso, W. Ananduta, S. Grammatico, and C. Ocampo-Martinez, “Operationally-safe peer-to-peer energy trading in distribution grids: A game-theoretic market-clearing mechanism,” _arXiv preprint arXiv:2107.13444_ , 2021.
  * [9] Z. Deng, “Distributed algorithm design for aggregative games of euler-lagrange systems and its application to smart grids,” _IEEE Transactions on Cybernetics_ , pp. 1–11, 2021.
  * [10] L. S. Shapley, “Stochastic games,” _Proceedings of the national academy of sciences_ , vol. 39, no. 10, pp. 1095–1100, 1953.
  * [11] J. Lei, U. V. Shanbhag, J.-S. Pang, and S. Sen, “On synchronous, asynchronous, and randomized best-response schemes for stochastic nash games,” _Mathematics of Operations Research_ , vol. 45, no. 1, pp. 157–190, 2020\.
  * [12] J. Lei and U. V. Shanbhag, “Linearly convergent variable sample-size schemes for stochastic Nash games: Best-response schemes and distributed gradient-response schemes,” in _2018 58th IEEE Conference on Decision and Control_. IEEE, 2018, pp. 3547–3552.
  * [13] E. Meigs, F. Parise, and A. Ozdaglar, “Learning in repeated stochastic network aggregative games,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 6918–6923.
  * [14] J. Lei, P. Yi, and L. Li, “Distributed no-regret learning for stochastic aggregative games over networks,” in _2021 40th Chinese Control Conference (CCC)_. IEEE, 2021, pp. 7512–7519.
  * [15] B. Franci and S. Grammatico, “Stochastic generalized nash equilibrium seeking under partial-decision information,” _Automatica_ , vol. 137, p. 110101, 2022\.
  * [16] P. Frihauf, M. Krstic, and T. Basar, “Nash equilibrium seeking in noncooperative games,” _IEEE Transactions on Automatic Control_ , vol. 57, no. 5, pp. 1192–1207, 2012.
  * [17] G. Belgioioso and S. Grammatico, “Semi-decentralized Nash equilibrium seeking in aggregative games with separable coupling constraints and non-differentiable cost functions,” _IEEE control systems letters_ , vol. 1, no. 2, pp. 400–405, 2017.
  * [18] M. Ye and G. Hu, “Distributed Nash equilibrium seeking by a consensus based approach,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 9, pp. 4811–4818, 2017.
  * [19] U. V. Shanbhag, J.-S. Pang, and S. Sen, “Inexact best-response schemes for stochastic Nash games: Linear convergence and iteration complexity analysis,” in _Decision and Control (CDC), 2016 IEEE 55th Conference on_. IEEE, 2016, pp. 3591–3596.
  * [20] F. Salehisadaghiani and L. Pavel, “Distributed Nash equilibrium seeking in networked graphical games,” _Automatica_ , vol. 87, pp. 17–24, 2018.
  * [21] P. Yi and L. Pavel, “An operator splitting approach for distributed generalized Nash equilibria computation,” _Automatica_ , vol. 102, pp. 111–121, 2019.
  * [22] Z. Zhou, P. Mertikopoulos, A. L. Moustakas, N. Bambos, and P. Glynn, “Mirror descent learning in continuous games,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017, pp. 5776–5783.
  * [23] J. Lei, U. V. Shanbhag, J.-S. Pang, and S. Sen, “On synchronous, asynchronous, and randomized best-response schemes for stochastic Nash games,” _Mathematics of Operations Research_ , vol. 45, no. 1, pp. 157–190, 2020\.
  * [24] T. Tatarenko, W. Shi, and A. Nedić, “Geometric convergence of gradient play algorithms for distributed Nash equilibrium seeking,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 11, pp. 5342–5353, 2021.
  * [25] J. Koshal, A. Nedić, and U. V. Shanbhag, “Distributed algorithms for aggregative games on graphs,” _Operations Research_ , vol. 64, no. 3, pp. 680–704, 2016.
  * [26] J. Lei, U. V. Shanbhag, and J. Chen, “Distributed computation of Nash equilibria for monotone aggregative games via iterative regularization,” in _2020 59th IEEE Conference on Decision and Control (CDC)_ , 2020, pp. 2285–2290.
  * [27] R. Zhu, J. Zhang, K. You, and T. Başar, “Asynchronous networked aggregative games,” _Automatica_ , vol. 136, Feb. 2022.
  * [28] S. Liang, P. Yi, and Y. Hong, “Distributed nash equilibrium seeking for aggregative games with coupled constraints,” _Automatica_ , vol. 85, pp. 179–185, 2017.
  * [29] J. Lei and U. V. Shanbhag, “Linearly convergent variable sample-size schemes for stochastic Nash games: Best-response schemes and distributed gradient-response schemes,” in _2018 IEEE Conference on Decision and Control (CDC)_ , 2018, pp. 3547–3552.
  * [30] C. D. Enrich, S. Jelassi, C. Domingo-Enrich, D. Scieur, A. Mensch, and J. Bruna, “Extragradient with player sampling for faster Nash equilibrium finding,” _arXiv preprint arXiv:1905.12363_ , 2019.
  * [31] Y. Malitsky, “Projected reflected gradient methods for monotone variational inequalities,” _SIAM Journal on Optimization_ , vol. 25, no. 1, pp. 502–520, 2015.
  * [32] G. Kotsalis, G. Lan, and T. Li, “Simple and optimal methods for stochastic variational inequalities, I: operator extrapolation,” _arXiv preprint arXiv:2011.02987_ , 2020.
  * [33] A. Nedić and A. Ozdaglar, “Distributed subgradient methods for multi-agent optimization,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 1, pp. 48–61, 2009.
  * [34] F. Facchinei and J.-S. Pang, _Finite-dimensional variational inequalities and complementarity problems. Vol. I_ , ser. Springer Series in Operations Research. New York: Springer-Verlag, 2003\.
  * [35] M. Bravo, D. Leslie, and P. Mertikopoulos, “Bandit learning in concave n-person games,” _Advances in Neural Information Processing Systems_ , vol. 31, pp. 5661–5671, 2018.
  * [36] G. Lan, _First-order and Stochastic Optimization Methods for Machine Learning_. Springer Nature, 2020.