A distributed primal-dual algorithm for computation of generalized Nash equilibria with
shared affine coupling constraints via operator splitting methods
Peng Yi, Lacra Pavela,
a Department

of Electrical and Computer Engineering, University of Toronto

arXiv:1703.05388v1 [math.OC] 15 Mar 2017

Abstract
In this paper, we propose a distributed primal-dual algorithm for computation of a generalized Nash equilibrium (GNE) in
noncooperative games over network systems. In the considered game, not only each player’s local objective function depends
on other players’ decisions, but also the feasible decision sets of all the players are coupled together with a globally shared affine
inequality constraint. Adopting the variational GNE, that is the solution of a variational inequality, as a refinement of GNE, we
introduce a primal-dual algorithm that players can use to seek it in a distributed manner. Each player only needs to know its local
objective function, local feasible set, and a local block of the affine constraint. Meanwhile, each player only needs to observe the
decisions on which its local objective function explicitly depends through the interference graph and share information related to
multipliers with its neighbors through a multiplier graph. Through a primal-dual analysis and an augmentation of variables, we
reformulate the problem as finding the zeros of a sum of monotone operators. Our distributed primal-dual algorithm is based on
forward-backward operator splitting methods. We prove its convergence to the variational GNE for fixed step-sizes under some
mild assumptions. Then a distributed algorithm with inertia is also introduced and analyzed for variational GNE seeking. Finally,
numerical simulations for network Cournot competition are given to illustrate the algorithm efficiency and performance.
Keywords: Network system; generalized Nash equilibrium; multi-agent systems; distributed algorithm; operator splitting;

1. Introduction
Engineering network systems, like power grids,
communication networks, transportation networks and
sensor networks, play a foundation role in modern society.
The efficient and secure operation of various network systems
relies on efficiently solving decision and control problems
arising in those large scale network systems. In many decision
problems, the nodes can be regarded as agents that need to
make local decisions possibly limited by the shared network
resources within local feasible sets. Meanwhile, each agent
has a local cost/utility function to be optimized, which depends
on the decisions of other agents. The traditional manner
for solving such decision problems over networks is the
centralized optimization approach, which relies on a control
center to gather the data of the problem and to optimize the
social welfare (usually taking the form of the sum of local
objective functions) within the local and global constraints.
The centralized optimization approach may not be suitable for
decision problems over large scale networks, since it needs
bidirectional communication between all the network nodes
and the control center, it is not robust to the failure of the
center node, and the computational burden for the center is
unbearable. It is also not preferable because the privacy of
each agent might be compromised when the data is transferred
I This

work was supported by NSERC Discovery Grant (261764).
Email address: peng.yi@utoronto.ca,
pavel@control.toronto.edu (Peng Yi, Lacra Pavel)
Preprint submitted to Elsevier

to the center. Recently, a distributed optimization approach is
proposed as an alternative methodology for solving decision
problems in network systems (Yi, Hong & Liu (2016), Shi,
Ling, Wu and Yin (2015) and Zeng, Yi, Hong & Xie (2016)).
In the distributed optimization approach, the data is distributed
throughout the network nodes and there is no control center,
and each agent in the network can just utilize its local data
and share information with its neighbour agents to compute its
local decision that corresponds to the optimal solution of the
social welfare optimization problem. Therefore, the distributed
optimization approach overcomes the drawbacks of the
centralized optimization approach by decomposing the data,
computation and communication to each agent. Moreover,
each agent has the authority and autonomy to formulate its own
objective function without worrying about privacy leaking out.
However, both approaches adopt the same solution concept,
that is the optimal social welfare solution with local and global
constraints, as the solution criterion of decision problems in
network systems.
However, optimal solutions of social welfare may not be
proper solution concepts in many applications. In fact, with the
deregulation and liberalization of markets, there is no guarantee
that the agents will not deliberately deviate from their local
optimal solutions to increase (decrease) own local utility (cost),
possibly by deceiving to utilize more network resources. In
this paper, we consider the game theoretic approach where
each agent in the network has its own local autonomy and
rationality. In such a setup of multiple interacting rational
March 17, 2017

players making decisions in a noncooperative environment,
Nash Equilibrium (NE) is a more reasonable solution. In
an NE, no player can increase (decrease) its local utility
(cost) by unilaterally changing its local decision, therefore,
no agent has the incentive to deviate from it. In other
words, NE is a self-enforceable solution in the sense that
once NE is computed all the agents will execute that NE.
Recently, there has been increasing interest in adopting game
theory and NE as the modeling framework and solution
concept for various network decision problems, like wireless
communication systems (Scutari,Palomar, Facchinei, and Pang,
(2010) and Menache & Ozdaglar (2011)), network flow
control (Alpcan & Basar (2005)), optical networks (Pan &
Pavel (2009)) and smart grids (Ye & Hu (2016)).
Moreover, in engineering network systems, not only the
local objective function of each agent depends on other
players’ decisions, but also the feasible set of each local
decision could depend on other agents’ decisions, because
the agents may compete for the utilization of some shared
or common network resources like bandwidth, spectrum or
power. This type of network decision problems can be
modeled as noncooperative games with coupling constraints,
and generalized Nash equilibrium (GNE) can be adopted as its
solution. The study of GNE dates back to the social equilibrium
concept proposed by Debreu (1952), and flourished in the
last two decades with applications in practical problems like
environment pollution games (Krawczyk & Uryasev (2000)),
power market design (Contreras, Klusch, Krawczyk (2004)),
optical networks (Pavel (2007)), wireless communication
(Facchinei & Pang (2010)). Interested readers can refer to
Facchinei & Kanzow (2010) for a historical review of GNE,
and refer to Fischer, Herrich & Schonefeld (2014) for recent
developments, and to Facchinei & Pang (2010) for a technical
treatment.
Even though NE or GNE is a reasonable and expectable
solution as a result of multiple rational agents making decisions
in a noncooperative manner, how to arrive at an NE (GNE) is
by no way a self-evident task. Each player needs to know the
complete game information, including objective functions and
feasible decision sets of all the other agents, in order to compute
NE in an introspective manner. It gets much more complicated
for computing GNE, because the agents also have to consider
the coupling in the feasible decision sets. In fact, as of yet there
is no universal manner to efficiently compute GNE in games
with coupling constraints (Harker (1991) and Fischer, Herrich
& Schonefeld (2014)), except for games with shared coupling
constraints (Facchinei, Fischer & Piccialli (2007)). Moreover,
for games in large scale network systems, it is quite unrealistic
and undesirable to assume that each agent could have the
complete information of the whole network, because this
implies prohibitive communication and computation burden
and no privacy protection. Therefore, each player (agent)
should compute its local decision corresponding with an NE
or GNE in a distributed manner, somehow resembling the
distributed optimization approach. In other words, each agent
should only utilize its local objective function, local feasible
set and possible local data related to coupling constraints, and

should only share information with its neighbouring agents to
compute its local decision in the NE (GNE). This turns out to be
an emerging research topic and gets studied in Salehisadaghiani
& Pavel (2016a), Koshal,Nedić & Shanbhag (2016), Parise,
Gentile, Grammatico & Lygeros (2015), Ye & Hu (2016) and
Swenson, Kar & Xavier (2015), etc.
Motivated by the above, in this work we consider a
distributed algorithm for iterative computation of GNE in
noncooperative games with shared affine coupling constraints
over network systems. The considered noncooperative game
has each agent’s local objective function depending on other
agents’ decisions as specified by an interference graph, and
has also an affine constraint shared by all agents, coupling all
players’ feasible decision sets. The considered game model
covers many practical problems, like the power market model
in Contreras, Klusch, Krawczyk (2004), environment pollution
game in Krawczyk & Uryasev (2000), power allocation
game in communication systems in Yin, Shanbhag & Mehta
(2011). A (centralized) numerical algorithm was recently
studied in Schiro, Pang & Shanbhag (2013) for quadratic
objective functions and in Dreves & Sudermann (2016) for
linear objective functions. Generally speaking, the GNE of the
considered game may not be unique. In this work, we adopt
the variational GNE, that corresponds with the solution of a
variational inequality proposed in Facchinei, Fischer & Piccialli
(2007), to be a refinement GNE solution. The variational GNE
is a particular type of the normalized equilibrium proposed in
Rosen (1965), and enjoys a nice economical interpretation that
all the agents have the same shadow price for shared network
resources without any discrimination as pointed in Kulkarni &
Shanbhag (2012). Furthermore, the variational GNE enjoys
a sensitivity and stability property (Facchinei & Pang (2010)
and Facchinei & Kanzow (2010)), hence we adopt it as the
desirable solution.
We propose a new type of distributed algorithm that agents
can use to compute the variational GNE by only manipulating
their local data and communicating with neighbouring agents.
Observing that the KKT condition of the corresponding
variational inequality requires all agents to reach consensus
on the multiplier of the shared affine constraint, we introduce
a local copy of the multiplier and an auxiliary variable for
each player. To enforce the consensus of local multipliers,
we use a reformulation that incorporates the Laplacian matrix
of a connected graph. Motivated by the forward-backward
operator splitting method for finding zeros of a sum of
monotone operators (refer to Bauschke & Combettes (2011))
and the recent primal-dual algorithm proposed in Condat
(2013) for optimization problems with linear composition
terms, we propose a novel distributed algorithm for iterative
computation of GNE. The main idea is to introduce a suitable
metric matrix and to split the equivalent reformulation into
two monotone operators. An operator splitting method has
been adopted for NE computation in a centralized manner in
Briceno-Arias & Combettes (2013). A different splitting idea
is adopted here appropriate for distributed GNE computation.
Moreover, a distributed algorithm with inertia is also proposed
and investigated, motivated by the acceleration algorithms in
2

Alvarez & Attouch (2001), Attouch, Chbani, Peypouquet, and
Redont (2016), Iutzeler & Hendrickx (2016) and Lorenz &
Pock (2015), most of which only focused on optimization
problems. The convergence of the proposed algorithms is
verified under suitable fixed step-size choice and some mild
assumptions on the objective functions and communication
graphs.
The recent works of Zhu & Frazzoli (2016), Yu, van
der Schaar & Sayed (2016), Liang, Yi and Hong (2016)
and Paccagnan, Gentile, Parise, Kamgarpour & Lygeros
(2016) are closely related with this work since all of them
are concerned with the distributed algorithm for seeking GNE
of noncooperative games with coupling constraints. Zhu
& Frazzoli (2016) address the GNE seeking for the case
where each player has non-shared local coupling constraints.
Assuming that each player can observe other players’ decisions
on which its local objective function and local constraint
functions depend through the interference graph, Zhu &
Frazzoli (2016) propose a distributed primal-dual GNE seeking
algorithm based on variational inequality methods, and show
algorithm convergence under diminishing step-sizes. Yu,
van der Schaar & Sayed (2016) investigate the distributed
GNE seeking under stochastic data observations. The authors
assume that the coupling constraints have a locally shared
property that if one player has its one of local constraints
dependent on the decision of another player, then this constraint
must be shared between those two players. Their algorithm
design is based on a penalty-type gradient method. Under
the assumption that each player can observe the decisions
on which its local objective function and constraint functions
depend, Yu, van der Schaar & Sayed (2016) utilize a gradient
type algorithm to seek the pure NE of the game derived by
penalizing the coupling and local constraints. They show that
their algorithm can reach a region near the pure penalized NE
with a constant step-size, which will approach a GNE if the
penalizing parameter goes to infinity. Both Liang, Yi and
Hong (2016) and Paccagnan, Gentile, Parise, Kamgarpour &
Lygeros (2016) consider the distributed algorithm for seeking a
variational GNE of the aggregative game with globally shared
affine coupling constraints. This represents a particular type of
game where the players’ local objective functions depend on
some aggregative variables of all agents’ decisions. Liang, Yi
and Hong (2016) assume that each player has local copies of
both the aggregative variables and the multipliers, and combine
a finite-time convergent continuous-time consensus dynamics
and a projected gradient flow to derive their distributed
GNE seeking dynamics. Meanwhile, Paccagnan, Gentile,
Parise, Kamgarpour & Lygeros (2016) adopt the asymmetric
projection algorithm for variational inequalities to design their
variational GNE seeking algorithm. However Paccagnan,
Gentile, Parise, Kamgarpour & Lygeros (2016) assume that
there is an additional central node for the update of the common
multiplier, and only address quadratic objective functions.
Compared with these works, our paper has following
contributions,
(i): The considered noncooperative game model is
completely general, thus a generalization of the aggregative

game in Liang, Yi and Hong (2016) and Paccagnan, Gentile,
Parise, Kamgarpour & Lygeros (2016). We further assume
that the shared affine coupling constraint is also decomposed
such that each player only knows its local contribution to the
global constraint, that is only a sub-block matrix of the whole
constraint matrix. In this sense, no player knows exactly the
shared constraints, hence, our problem model is also different
from the ones in Yu, van der Schaar & Sayed (2016) and
Zhu & Frazzoli (2016). The decomposition of the coupling
constraints, together with the localization of player’s local
objective function and local feasibility set, is quite appealing for
iterative computation of GNE in large-scale network systems
because this reduces the data transmission and computation
burden, and protects the players’ privacies.
(ii):The proposed distributed algorithms can compute the
variational GNE iteratively under a more localized data
structure and information observing structure compared to
previous ones. Firstly, each player only utilizes the local
objective function and local feasible set, and its local sub-block
matrix of the affine constraints. Secondly, we assume the
players have two (different) information observing graphs, i.e.,
interference graph and multiplier graph. Each player only
needs to observe the decisions that its local objective function
directly depends on through the interference graph. This type
of information observation assumption has also been adopted in
Yu, van der Schaar & Sayed (2016) and Zhu & Frazzoli (2016).
Meanwhile, each player only needs to share information related
to multipliers with its neighbouring agents through another
multiplier graph. Here it is not required that each player
should know the decisions that coupling constraints depend
on, as assumed in Yu, van der Schaar & Sayed (2016) and
Zhu & Frazzoli (2016). Therefore, our information sharing
(observing) structure is more localized and sparse.
(iii): The algorithm development and convergence analysis
is motivated by the operator splitting method (Bauschke &
Combettes (2011)), different from the penalized method
adopted in Yu, van der Schaar & Sayed (2016) and the
variational inequality approach in Liang, Yi and Hong (2016),
Zhu & Frazzoli (2016) and Paccagnan, Gentile, Parise,
Kamgarpour & Lygeros (2016). Based on this operator
splitting approach, we prove the algorithm converges to the
variational GNE under fixed step-sizes. Note that neither
convergence nor non-bias estimation is achieved in Yu, van
der Schaar & Sayed (2016) under fixed step-sizes, while
Zhu & Frazzoli (2016) achieve convergence with diminishing
step-sizes. On the other hand, compared with Briceno-Arias &
Combettes (2013), this paper addresses the GNE seeking under
coupling constraints, adopts a different splitting technique,
and achieves fully distributed computations. The operator
splitting method is powerful and provides additional insights.
Moreover, a distributed algorithm with inertia is proposed
and analyzed, resembling the acceleration algorithms in
optimization (Nesterov (2013) and Iutzeler & Hendrickx
(2016)). The algorithm performance is illustrated via numerical
experiments of network Cournot competitions with bounded
market capacities.
The paper is organized as follows. Section 2 gives the
3

domRA = Rm if A is maximally monotone 1 . For a proper lower
semi-continuous convex (l.s.c.) function f , its sub-differential
m
is a set-valued operator ∂ f : dom f → 2R and

notations and preliminary background. Section 3 formulates
the noncooperative game and gives the distributed algorithm
for iterative computation of a GNE. Section 4 shows how the
operator splitting method motivates the algorithm development,
and Section 5 presents the algorithm convergence analysis.
Then a distributed GNE seeking algorithm with inertia is
proposed and analyzed in Section 6. Finally, a network Cournot
competition with bounded market capacities is formulated with
numerical studies in Section 7, while concluding remarks are
given in Section 8.

∂ f : x 7→ {g ∈ Rm | f (y) ≥ f (x) + hg, y − xi, ∀y ∈ dom f }. (1)
∂ f is a maximally monotone operator 2 . Then Prox f = R∂ f :
Rm → dom f is called the proximal operator of f 3 , i.e.
1
Prox f : x 7→ arg min f (u) + ||u − x||22 .
(2)
u∈dom f
2



0, x ∈ Ω;
Define the indicator function of set Ω as ιΩ (x) = 

∞, x < Ω.
For a closed convex set Ω, ιΩ is a proper l.s.c. function. ∂ιΩ is
just the normal cone operator of set Ω, that is ∂ιΩ (x) = NΩ (x) 4
and



∅
x<Ω




(3)
NΩ (x) = 
{v|hv,
y
−
xi
≤
0,
∀y
∈
Ω}
x
∈ bd(Ω)




0
x ∈ int(Ω)

2. Notations and preliminary background
In this section, we review the notations and some preliminary
notions in monotone operators and graph theory.
Notations: In the following, Rm (Rm
+ ) denotes the
m−dimesional (nonnegative) Euclidean space. For a column
vector x ∈ Rm (matrix A ∈ Rm×n ), xT (AT ) denotes its transpose.
√
xT y = hx, yi denotes the inner product of x, y, and ||x||2 = xT x
denotes the norm induced by inner product h·, ·i. Given a
symmetric positive definite matrix G, denote the G-induced
inner product hx, yiG = hGx,
√ yi. The G-matrix induced norm,
|| · ||G , is defined as ||x||G = hGx, xi. Denote by || · || any matrix
induced norm in the Euclidean space. Denote 1m = (1, ..., 1)T ∈
Rm and 0m = (0, ..., 0)T ∈ Rm . For column vectors x, y, x ≥ (>)y
is understood componentwise. diag{A1 , ..., AN } represents the
block diagonal matrix with matrices A1 , ..., AN on its main
diagonal. Null(A) and Range(A) denote the null space and
range space of matrix A, respectively. Denote col(x1 , ...., xN )
as the column vector stacked with column vectors x1 , ..., xN . In
denotes the identity matrix in Rn×n . For a matrix A = [ai j ], ai j or
[A]i j stands for the matrix entry in the ith row and jth column
of A. We also use [x]k to denote the k−th element in column
Q
vector x. Denote ×i=1,...,n Ωi or ni=1 Ωi as the Cartesian product
of the sets Ωi , i = 1, ..., n. Denote int(Ω) as the interior of Ω, and
bd(Ω) as the boundary set of Ω. Define the projection of x onto
a set Ω by PΩ (x) = arg miny∈Ω ||x − y||2 . A set Ω is a convex set
if λx + (1 − λ)y ∈ Ω, ∀λ ∈ [0, 1], ∀x, y ∈ Ω. An extended
value proper function f : Rm → R is a convex function if
f (λx+(1−λ)y) ≤ λ f (x)+(1−λ) f (y), ∀λ ∈ [0, 1], ∀x, y ∈ dom f .

In this case, we also have 5
ProxιΩ (x) = RNΩ (x) = PΩ (x).

(4)

For a single-valued operator T : Rm → Rm , a point x ∈ Rm
is a fixed point of T if T x = x, and the set of fixed points of
T is denoted as FixT . The composition of operators A and
B, denoted by A ◦ B, is defined via its graph graA ◦ B =
{(x, z)|∃y ∈ ranB, (x, y) ∈ graB, (y, z) ∈ graA}. We also use AB
to denote the composition A ◦ B when they are single-valued.
Similarly, their sum A + B is defined as gra(A + B) = {(x, y +
z)|(x, y) ∈ graA, (x, z) ∈ graB}. Suppose operators A and B
are maximally monotone and 0 ∈ int(domA − domB), then
A + B is also maximally monotone6 . Further suppose that
A is single-valued, then zer(A + B) = FixRB ◦ (Id − A)7 ,
which helps to formulate the basic forward-backward operator
splitting algorithm for finding zeros of a sum of monotone
operators.
2.2. Graph theory
The following concepts are reviewed from Mesbahi &
Egerstedt (2010). The information sharing or exchanging
among the agents is described by graph G = (N, E). N =
{1, · · · , N} is the set of agents, and the edge set E ⊂ N × N
contains all the information interactions. If agent i can get
information from agent j, then ( j, i) ∈ E and agent j belongs
to agent i’s neighbor set Ni = { j|( j, i) ∈ E}, and i < Ni . G is
said to be undirected when (i, j) ∈ E if and only if ( j, i) ∈ E. A

2.1. Monotone operators
The following concepts are reviewed from Bauschke &
m
Combettes (2011). Let A : Rm → 2R be a set-valued operator.
Denote Id as the identity operator, i.e, Id(x) = x. The domain
of A is domA = {x ∈ Rm |Ax , ∅} where ∅ stands for the empty
set, and the range of A is ranA = {y ∈ Rm |∃x ∈ Rm , y ∈ Ax}.
The graph of A is graA = {(x, u) ∈ Rm × Rm |u ∈ Ax}, then
the inverse of A is defined through its graph as graA−1 =
{(u, x) ∈ Rm × Rm |(x, u) ∈ graA}. The zero set of operator
A is zerA = {x ∈ Rm |0 ∈ Ax}. Define the resolvent of operator
A as RA = (Id + A)−1 . An operator A is called monotone if
∀(x, u), ∀(y, v) ∈ graA, we have hx − y, u − vi ≥ 0. Moreover, it
is maximally monotone if graA is not strictly contained in the
graph of any other monotone operator. RA is single-valued and

1

Proposition 23.7 in Bauschke & Combettes (2011)
Theorem 20.40 in Bauschke & Combettes (2011)
3 Proposition 16.34 in Bauschke & Combettes (2011)
4 Example 16.12 of Bauschke & Combettes (2011)
5 Example 23.4 of Bauschke & Combettes (2011)
6 Corollary 24.4 in Bauschke & Combettes (2011)
7 Proposition 25.1 in Bauschke & Combettes (2011)
2

4

path of graph G is a sequence of distinct agents in N such that
any consecutive agents in the sequence correspond to an edge
of graph G. Agent j is said to be connected to agent i if there is
a path from j to i. G is said to be connected if any two agents
are connected.
Define the weighted adjacency matrix W = [wi j ] ∈ RN×N of
G with wi j > 0 if j ∈ Ni and wi j = 0 otherwise. Assume W =
W T for undirected graphs. Define the weighted degree matrix
P
P
Deg = diag{d1 , · · · , dN } = diag{ Nj=1 w1 j , ..., Nj=1 wN j }. Then
the weighted Laplacian of graph G is L = Deg−W. When graph
G is a connected and undirected graph, 0 is a simple eigenvalue
of Laplacian L with the eigenspace {α1N |α ∈ R}, and L1N = 0N ,
1TN L = 0TN , while all other eigenvalues are positive. Denote the
eigenvalues of L in an ascending order as 0 < s2 ≤ · · · ≤ sN ,
then by Courant-Fischer Theorem,
min xT Lx = s2 ||x||22 ,
x,0,
1T x=0

max xT Lx = sN ||x||22 .
x,0

by globally shared affine constraints. Denote
X ⊂ Rm :=

where Ωi ⊂ R is a private feasible decision set of player i, and
A = [A1 , · · · , AN ] ∈ Rm×n with Ai ∈ Rm×ni , and b ∈ Rm . Denote
QN
Ω = i=1
Ωi . Given the globally shared set X (which may not
be known by any agents), the following set-valued map gives
the feasible decision set map of agent i: Xi (x−i ) := {xi ∈ Rni :
(xi , x−i ) ∈ X}, or in other words:
X
A j x j }.
Xi (x−i ) := {xi ∈ Ωi |Ai xi ≥ b −
j,i, j∈N

Hence, each agent has a local feasible constraint xi ∈ Ωi ,
and there exists a coupling constraint shared by all agents with
sub-matrix Ai characterizing how agent i is involved in the
coupling constraint (shares the global resource). Notice that
agent i may only know its local Ai , in which case the globally
shared affine constraint couples the agents’ feasible decision
sets, but is not known by any agents.

(5)

Remark 3.1. We consider affine coupling constraints for
various reasons. Even though not as general as the nonlinear
constraints considered in Pavel (2007) and Zhu & Frazzoli
(2016), this setup does enjoy quite strong modeling flexibility.
As pointed out in page 191 of Facchinei & Kanzow (2010), “
However, it should be noted that the jointly convex assumption
on the constraints · · · practically is likely to be satisfied only
when the joint constraints gµ = g, µ = 1, ..., N are linear, i.e. of
the form Ax ≤ b for some suitable matrix A and vector b.” In
fact, many existing generalized Nash game models adopt affine
coupling constraints, as well documented in Schiro, Pang &
Shanbhag (2013) and Dreves & Sudermann (2016).

(6)

3.1. Game formulation
Consider a group of agents (players) N = {1, · · · , N}
that seek the generalized Nash equilibrium (GNE) of a
noncooperative game with coupling constraints defined as
follows. Each player i ∈ N controls its local decision
(strategy or action) xi ∈ Rni . Denote x = col(x1 , · · · , xN ) ∈
Rn as the decision profile, i.e., the stacked vector of all
PN
ni = n. Denote x−i =
the agents’ decisions where i=1
col(x1 , · · · , xi−1 , xi+1 , · · · , xN ) as the decision profile stacked
vector of the agents’ decisions except player i. Agent i aims to
optimize its local objective function within its feasible decision
set. The local objective function for agent i is fi (xi , x−i ) : Rn →
R. Notice that the local objective function of agent i is coupled
with other players’ decisions (however, may not be explicitly
coupled with all other players’ decisions). Moreover, the
feasible decision set of player i also depends on the decisions
n
of the other players with Xi (x−i ) : Rn−ni → 2R i denoting a
set-valued map that maps x−i to the feasible decision set of
agent i. The aim of agent i is to find the best-response strategy
set given the other players’ decision x−i ,
xi ∈R

i

s.t., xi ∈ Xi (x−i ).

Assumption 1. For player i, fi (xi , x−i ) is a differentiable
convex function with respect to xi given any fixed x−i , and Ωi
is a closed convex set. X has nonempty interior point (Slater’s
condition).
Suppose x∗ is a GNE of game (7), then for agent i, xi∗ is the
optimal solution to the following convex optimization problem:
X
min fi (xi , x∗−i ), s.t. xi ∈ Ωi , Ai xi ≥ b −
A j x∗j .
(9)
xi

∀i ∈ N.

j,i, j∈N

Define a local Lagrangian function for agent i with multiplier
λi ∈ Rm
+ as
Li (xi , λi ; x−i ) = fi (xi , x−i ) + λTi (b − Ax).

(7)

λ∗i

xi∗

(10)
Rm
+

When is an optimal solution to (9), there exists ∈
such
that the following optimality conditions (KKT) are satisfied:

The GNE x∗ = col(x1∗ , · · · , x∗N ) of the game in (7) is obtained
at the intersection of all the players’ best-response sets, and is
defined as:
xi∗ ∈ arg min fi (xi , x∗−i ), s.t. xi ∈ Xi (x∗−i ),

\
{x ∈ Rn |Ax ≥ b},

i=1

3. Problem formulation and distributed algorithm

minn fi (xi , x−i )

Ωi

ni

Denote d∗ = max{d1 , · · · , dN } as the maximal weighted degree
of graph G, then we have the following estimation,
d∗ ≤ sN ≤ 2d∗ .

N
Y

0 ∈ ∇ xi Li (xi∗ , λ∗i ; x∗−i ) + NΩi (xi∗ ), xi∗ ∈ Ωi
hλ∗i , b − Ax∗ i = 0, b − Ax∗ ≤ 0, λ∗i ≥ 0,

(11)

These can be equivalently written in the following form by
using (10) and the definition of the normal cone operator in (3)

(8)

0 ∈ ∇ xi fi (xi∗ , x∗−i ) − ATi λ∗i + NΩi (xi∗ )
0 ∈ (Ax∗ − b) + NRm+ (λ∗i )

Here we consider the GNE seeking in noncooperative games
where the couplings between players’ feasible sets are specified
5

(12)

In fact, since NΩ (x) = ∅ when x < Ω, it must hold that xi∗ ∈
Ωi and λ∗i ∈ Rm
=
+ when (12) is satisfied. Furthermore, NRm
+
Qm
∗
∗
N
.
If
[λ
]
=
0,
then
N
([λ
]
)
=
−R
,
and
hence
R+
+
j=1 R+
i k
i k
[b − Ax∗ ]k ≤ 0; and if [λ∗i ]k > 0, we have NR+ ([λ∗i ]k ) = 0, hence
[b − Ax∗ ]k = 0. Therefore, b − Ax∗ ≤ 0 and hλ∗i , b − Ax∗ i = 0.
Denote λ̄ = col(λ1 , · · · , λN ). Therefore, by Theorem 4.6 in
Facchinei & Kanzow (2010) when (x∗ , λ̄∗ ) satisfies KKT (12)
for i = 1, ..., N, x∗ is a GNE of the game in (7)
According to the above discussions, given x∗ as a GNE of
game in (7), its corresponding Lagrangian multipliers for the
globally shared affine coupling constraint may be different for
the agents, i.e., λ∗1 , λ∗2 ,, · · · , , λ∗N . In this work, we aim
to seek a GNE with the same Lagrangian multiplier for all the
agents, which has a nice interpretation from the viewpoint of
variational inequality.
Define
F(x) = col(∇ x1 f1 (x1 , x−1 ), · · · , ∇ xN fN (xN , x−N )),

the coupling constraint is met and a variational GNE of the
game is found.
n
m
Define two operators A and B, both from Ω × Rm
+ to R × R
as follows,
!
!
x
F(x)
A:
7→
λ !
−b
!
!
!
(17)
0 −AT
x
x
NΩ (x)
+
B:
7→
A
0
λ
λ
NRm+ (λ))
By the definition of F(x) in (13), the KKT conditions (16)
can be equivalently written as 0 ∈ (A + B)col(x∗ , λ∗ ). Notice
that B is a maximally monotone operator (similar arguments
for this can be found in Lemma 5.4). Hence, if F(x) has some
additional properties, then solving V I(F, X) can be converted to
the problem of finding zeros of a sum of monotone operators.
Assumption 2. F(x) defined in (13) is strongly monotone with
parameter η over Ω: hF(x) − F(y), x − yi ≥ η||x − y||22 , ∀x, y ∈ Ω,
and θ− Lipschitz continuous over Ω: ||F(x) − F(y)||2 ≤ θ||x −
y||2 , ∀x, y ∈ Ω.

(13)

which is usually called the pseudo-gradient. The variational
inequality (VI) approach to find a GNE of game (7) is to find
the solution of the following V I(F, X):
Find x∗ ∈ X, hF(x∗ ), x − x∗ i ≥ 0, ∀x ∈ X.

Remark 3.3. Assumption 2 has also been adopted in Yu, van
der Schaar & Sayed (2016), Paccagnan, Gentile, Parise,
Kamgarpour & Lygeros (2016) and Zhu & Frazzoli (2016).
Assumption 2 guarantees that there exists a unique solution
to V I(F, X) in (14) 9 , thus guarantees the existence of a GNE
for game in (7). However, the GNE of (7) may not be unique.
The algorithm for computing all the GNE of noncooperative
games with coupling constraints is still an opening research
topic, and interested readers can refer to Nabetani, Tseng, and
Fukushima (2011). This work aims to provide a distributed
algorithm for iterative computation of a variational GNE of the
considered game, which enjoys nice economical interpretations
and stability properties.

(14)

Let us check the KKT condition for V I(F, X) in (14). In fact, x∗
is a solution to V I(F, X) in (14) if and only if x∗ is the optimal
solution to the following optimization problem:
minn hF(x∗ ), xi,
x∈R

s.t.,

x ∈ Ω, Ax ≥ b

(15)

According to the optimization formulation of V I(F, X) in (15),
if x∗ solves V I(F, X), there exists λ∗ ∈ Rm such that the
following optimality conditions (KKT) are satisfied:
0 ∈ ∇ xi fi (xi∗ , x∗−i ) − ATi λ∗ + NΩi (xi∗ ), i = 1, · · · , N
0 ∈ (Ax∗ − b) + NRm+ (λ∗ )

(16)
3.2. Distributed algorithm

By comparing the two sets of KKT conditions in (12) and
(16) we obtain,

In practice, each player only knows its private information
in game (7), especially when the players interact over large
scale networks. It is quite natural that each player can only
know its local objective function fi (xi , x−i ) and local feasible set
Ωi , which cannot be shared with other players, because those
data contain its local private information, such as cost function,
preference and action ability. Moreover, matrix Ai specifies
how player i participates in the resource allocation or market
behavior, hence also contains private information, and b can be
PN
decomposed as b = i=1
bi . Thus matrix Ai can be regarded as a
map from decision space Rni to resource space Rm , while vector
bi can be regarded as a local contribution or observation for
the global resource. For example, if there are total m markets,
and each player i produces a kind of product with amount of
T
xi ∈ R+ , then Ai ∈ Rm
+ satisfying 1m Ai = 1 and 0 ≤ Ai ≤ 1
just specifies how each player allocates its production to each
market. In this case, if b̃i ∈ Rm
+ is a local observation of market

Theorem 3.2. 8 Suppose Assumption 1 holds. Every solution
x∗ of V I(F, X) in (14) is a GNE of game in (7). Furthermore,
if x∗ together with λ∗ satisfies the KKT conditions for the
V I(F, X), i.e., (16), then x∗ together with λ∗1 =, · · · , = λ∗N = λ∗
satisfies the KKT conditions for the GNE, i.e., (12).
The solution x∗ of V I(F, X) in (14) is termed as a variational
GNE or normalized equilibrium of the game with coupling
constraints in (7). A variational GNE enjoys an economical
interpretations of no price discrimination and has better stability
and sensitivity properties, therefore, can be regarded as a
refinement of GNE (refer to Kulkarni & Shanbhag (2012) for
more discussions). This paper will propose a novel distributed
algorithm for the agents to find a solution of V I(F, X) in (14),
thus provides a distributed coordination mechanism such that
8 Theorem

9 Theorem

2.1 of Facchinei, Fischer & Piccialli (2007)

6

2.3.3 of Facchinei and Pang (2007)

coupling constraint: Ax ≥ b. To ensure that the globally shared
coupling constraint is satisfied and all the agents have the same
local multipliers, all players must coordinate which necessarily
requires that multiplier graph Gλ must be connected. Therefore,
G f and Gλ could be two different information observation or
information sharing graphs because they serve for different
purposes.

capacity vector, then the true market capacities can be taken as
PN
PN 1
b = i=1
bi = i=1
N b̃i .
Therefore, we assume that player i only knows its local
fi (xi , x−i ), Ωi , and matrix Ai and bi with A = [A1 , · · · , AN ] and
PN
b = i=1
bi . In other words, player i has a local first-order
oracle of fi (xi , x−i ) which returns ∇ xi fi (xi , x−i ) given (xi , x−i ),
meanwhile player i can manipulate Ωi , Ai and bi for its local
computation.

We are ready to present the main distributed algorithm
after the introduction of algorithm notations. Agent (player)
i controls its local decision variable xi ∈ Rni and local
Lagragian multiplier λi ∈ Rm . Meanwhile, we also assume
that player i has a local auxiliary variable zi ∈ Rm for the
coordinations needed to satisfy the affine coupling constraint
and to reach consensus of the local multipliers λi . As indicated
before, player i can compute ∇ xi fi (xi , x−i ) by observing the
adversary players’ decisions that its local objective function
fi (xi , x−i ) directly depends on, that is the decisions of its
interference neighbors in Nif . On the other hand, player i can
also share information related to the local multiplier λi and
local auxiliary variable zi with its multiplier neighbours in Niλ
through multiplier network Gλ .
Next we show the distributed algorithm for agent i.

The players need to find the solution to V I(F, X) in (14)
in a distributed manner by local information observation and
sharing, hence find a variational GNE of the game in (7) without
any coordinator. To facilitate the local coordination between
agents, here we specify two graphs, G f and Gλ , related to the
local information observations or exchanging between players.
Graph G f , termed as interference graph, is defined according
to the dependence relationships between the agents’ objective
functions and the other players’ decisions, which is also called
graphical model for games in computer science (refer to
Kearns, Littman, & Singh (2001)). For graph G f = (N, E f ),
( j, i) ∈ E f if the objective function of agent i, fi (xi , x−i )
explicitly depends on the decision of player j. We define
Nif = { j|( j, i) ∈ E f } as the set of interference neighbors
whose decisions directly influence the objective function of
player i. Therefore, the objective function of player i can also
be written as fi (xi , {x j } j∈N f ), and the local oracle of player i
i
returns ∇ xi fi (xi , {x j } j∈N f ), i.e., ∇ xi fi (xi , x−i ), given {x j } j∈N f . The
i
i
local oracle might compute ∇ xi fi (xi , {x j } j∈N f ) by approximating
i
∇ xi fi (xi , x−i ) with local objective function value observations
(taking the simultaneous perturbation stochastic approximation
in Spall (1992) as an example), or by utilizing the estimation
techniques developed in Salehisadaghiani & Pavel (2016b).
On the other hand, for the coordination of the feasibility of
action sets and the consensus of local multipliers λ∗1 =, · · · , =
λ∗N = λ∗ in Theorem 3.2, we also assume that the agents can
exchange certain local information through a multiplier graph
Gλ = (N, Eλ ). ( j, i) ∈ Eλ if player i can receive certain
information from player j, while the information to be shared
through Gλ will be specified later. Thereby, player i has its
multiplier neighbors Niλ = { j|( j, i) ∈ Eλ }. W = [wi j ] is the
weighted adjacency matrix associated with multiplier graph Gλ ,
and L is the corresponding weighted Laplacian matrix.

Algorithm 3.5.


xi,k+1 = PΩi xi,k − τi (∇ xi fi (xi,k , {x j,k } j∈N f ) − ATi λi,k )
i
X
zi,k+1 = zi,k + νi
wi j (λi,k − λ j,k )

(18)
(19)

j∈Niλ

X
n

λi,k+1 = PRm+ λi,k − σi 2Ai xi,k+1 − Ai xi,k − bi +
wi j (λi,k − λ j,k )
j∈Niλ

+2

X
j∈Niλ

wi j (zi,k+1 − z j,k+1 ) −

X

wi j (zi,k − z j,k )

o

(20)

j∈Niλ

τi , νi , σi > 0 are fixed constant step-sizes of player i, and
W = [wi j ] is the weighted adjacency matrix of Gλ .
Algorithm 3.5 runs sequentially as follows. At the iteration
time k, player i gets ∇ xi fi (xi,k , x−i,k ) by observing x j,k , j ∈ Nif
through interference graph G f , and updates xi,k with (18);
meanwhile, player i gets λ j,k , j ∈ Niλ through multiplier graph
Gλ , and updates zi,k by (19). Then player i gets z j,k+1 , j ∈ Niλ
through multiplier graph Gλ and updates λi,k with (20) that also
employs the most recent local information xi,k+1 .
Intuitively speaking, (18) employs the projected gradient
descent of the local Lagrangian function in (10). (19) can be
regarded as the discrete-time intergration for the consensual
errors of local copies of multipliers, which will ensure the
consensus of λi eventually. In fact, a similar dynamics has
been employed in distributed optimization in Lei, Chen & Fang
(2016). Finally, (20) updates local multiplier by a combination
of the projected gradient ascent of local Lagrangian function
(10) and a proportional-integral dynamics for consensual errors
of multipliers. Section 4 will give a detailed algorithm
development from the viewpoint of operator splitting methods.

Assumption 3. Gλ = {N, Eλ } is undirected and connected.
W = WT .
Remark 3.4. We assume that each agent can observe the
decisions which its local objective function directly depends
on through interference graph G f . Therefore, player i can
get its local gradient ∇ xi fi (xi , {x j } j∈N f ). This type of local
i
information observation model has also been adopted in Yu,
van der Schaar & Sayed (2016) and Zhu & Frazzoli (2016).
On the other hand, player i’s feasible decision set may depend
on any other player k’s decision even if fi (xi , x−i ) does not
explicitly depend on the decision of player k, i.e. k < Nif .
In fact, player i’s feasible decision set implicitly depends on
all other players’ decisions through the globally shared affine
7

Algorithm 4.1.

Algorithm 3.5 is a totally distributed algorithm and has
following key features:
i). The full data decomposition and privacy protection
is achieved since each player only needs to know its local
objective function fi (xi , x−i ) and local feasible set Ωi .
ii). The matrix A is decomposed and each block Ai is
kept by player i, hence the privacy of each player is protected
because Ai describes how player i is involved in the market or
competition.
iii). Each player only needs to observe the decisions which
its local objective function directly depends on, and only needs
to share information with its multiplier neighbours, through
G f and Gλ , respectively. Both graphs usually have sparse
edge connections, therefore, the observation or communication
burden is relieved. Furthermore, the information observation
related with decisions and the information sharing related
with multipliers is decoupled and accomplished with different
graphs G f and Gλ , respectively.
Therefore, those two
information sharing processes can work in a parallel manner,
and can be designed independently.
iv). The algorithm converges with fixed step-sizes under
some mild conditions, and works in a Gauss-Seidel manner
that utilizes the most recent information when updating λi .
Moreover, (18) and (19) can even be computed in parallel for
player i.



xk+1 = PΩ xk − τ̄(F(xk ) − ΛT λ̄k )

z̄k+1 = z̄k + ν̄L̄λ̄k
(22)
n

o
λ̄k+1 = PRmN
λ̄k − σ̄ Λ(2xk+1 − xk ) − b̄ + L̄λ̄k + L̄(2z̄k+1 − z̄k )
+
(23)
Using the fact that PΩ (x) = ProxιΩ (x) = RNΩ (x) in (4) and
the definition of resolvent operator as RNΩ (x) = (Id + NΩ )−1 ,
equation (21) be be written as xk+1 = (Id + NΩ )−1 [xk − τ̄(F(xk ) −
ΛT λ̄k )], or equivalently,
xk − τ̄(F(xk ) − ΛT λ̄k ) ∈ xk+1 + NΩ (xk+1 ).

(24)

diag{ τ11 In1 , · · ·
−1

, τ1N InN }, ν̄−1
=
1
1
diag{ σ1 Im , · · · , σN Im }.
QN
= i=1
NΩi (xi ), hence

Notice that τ̄
=
1
1
diag{ ν1 Im , · · · , νN Im } and σ̄
=
Furthermore, NΩ is a cone and NΩ (x)
τ̄−1 NΩ (x) = NΩ (x). Therefore, (24) can be written as
−1

− F(xk ) ∈ NΩ (xk+1 ) − ΛT λ̄k+1 + τ̄−1 (xk+1 − xk ) + ΛT (λ̄k+1 − λ̄k ).
(25)
QN
m (λ ), and
mN (λ̄) =
Moreover, NRmN
is
a
cone,
N
N
R
i
R+
i=1
+
+
σ̄−1 NRmN
= NRmN
. Then with similar arguments, equation (23)
+
+
can be written as:
h
i
λ̄k − σ̄ Λ(2xk+1 −xk )− b̄+ L̄λ̄k + L̄(2z̄k+1 − z̄k ) ∈ λ̄k+1 +NRmN
(λ̄k+1 )
+
(26)
or equivalently,

4. Algorithm development
In this section, we first show how Algorithm 3.5 is developed
and provide the motivations behind the algorithm’s convergence
analysis. Then we verify that the limiting point of Algorithm
3.5 solves the V I(F, X) in (14), and thus finds a variational GNE
of the game in (7).
Algorithm 3.5 is inspired by the forward-backward splitting
methods for finding zeros of the sum of monotone operators
(Bauschke & Combettes (2011)) and the primal-dual algorithm
for optimization with linear composition terms by Condat
(2013). The key difference are the specific operator splitting
form and the augmentation of variables to achieve distributed
computations. Next, we systematically show how to write
Algorithm 3.5 in the form of a forward-backward operator
splitting algorithm.
Let us define some notations to write Algorithm 3.5 in a
compact form. Denote xk = col(x1,k , · · · , xN,k ) ∈ Rn , λ̄k =
col(λ1,k , · · · , λN,k ) ∈ RmN , z̄k = col(z1,k , · · · , zN,k ) ∈ RmN , b̄ =
col(b1 , · · · , bN ) ∈ RmN , τ̄ = diag{τ1 In1 , · · · , τN InN } ∈ Rn×n , ν̄ =
diag{ν1 Im , · · · , νN Im } ∈ RmN×mN , σ̄ = diag{σ1 Im , · · · , σN Im } ∈
RmN×mN , Λ = diag{A1 , · · · , AN } ∈ RmN×n and ΛT =
diag{AT1 , · · · , ATN } ∈ Rn×mN . L̄ = L ⊗ Im where L ∈ RN×N is
the weighted Laplacian matrix of multiplier graph Gλ .
Using these notations, the definition of pseudo-gradient
N
F(x)
in
(13)
and
PQi=1
=
Ωi (col(x1 , · · · , xN ))
10
col(PΩ1 (x1 ), · · · , PΩN (xN )) , Algorithm 3.5 can be written in a
compact form as:
10 Proposition

(21)

−[L̄λ̄k − b̄] ∈
+

NRmN
(λ̄k+1 ) + Λxk+1 + L̄z̄k+1
+
Λ(xk+1 − xk ) + L̄(z̄k+1 − z̄k ) + σ̄−1 (λ̄k+1 − λ̄k )
(27)
Therefore, equation(22) together with (25) and (27) can be
written in a compact form as:

 

NΩ (xk+1 ) − ΛT λ̄k+1
 F(xk )  


 

−
L̄
λ̄
0
− 
∈
k+1



 
 ¯
NRmN
(
λ̄
)
+
Λx
+
L̄z̄
L̄λk − b̄
k+1
k+1
k+1
+
 −1


 τ̄
0
ΛT   xk+1 − xk 




L̄   z̄k+1 − z̄k 
+  0 ν̄−1
(28)




Λ
L̄ σ̄−1
λ̄k+1 − λ̄k
Denote

 −1
 τ̄

Φ =  0

Λ

0
ν̄−1
L̄

ΛT
L̄
σ̄−1






(29)

Notice that matrix Φ ∈ R(n+2mN)×(n+2mN) is symmetric due to
L = LT and L̄ = L ⊗ Im .
Denote Ω̄ = Ω × RmN × RmN
+ . Define the operators Ā : Ω̄ →
n+2mN
and B̄ : Ω̄ → 2R
as follows,
R




 x 
 F(x) 



 ,
0
A :=  z̄  7→ 




λ̄
L̄λ̄ − b̄



 


 x 
 NΩ (x)   0 0 −ΛT   x 



 +  0 0 −L̄   z̄  (30)
0
B :=  z̄  7→ 
 
 




NRmN
(λ̄)
λ̄
Λ L̄
0
λ̄
+
n+2mN

23.16 of Bauschke & Combettes (2011)

8

Remark 4.2. Operators Ā and B̄ in (30) can be regarded as an
extension of operators A and B in (17) by augmenting λ ∈ Rm
to λ̄ ∈ RmN and introducing auxiliary variables z̄ ∈ RmN .
Moreover, operators in (30) utilize L1N = 0N of Laplacian
matrix L to ensure the consensus of λi , and utilize 1TN L = 0TN
to ensure the feasibility of affine coupling constraints.

The next result shows that any limiting point of Algorithm
3.5, that is, any zero point of operator Ā + B̄, is a variational
GNE of game (7).
Theorem 4.5. Suppose that Assumptions 1-3 hold. Consider
operators Ā and B̄ defined in (30), and operators A and B
defined in (17). Then the following statements hold:
(i): Given any col(x∗ , z̄∗ , λ̄∗ ) ∈ zer(Ā + B̄), then x∗ solves the
V I(F, X) in (14), hence x∗ is a variational GNE of game in (7).
Moreover λ̄∗ = 1N ⊗ λ∗ , and the multiplier λ∗ together with x∗
satisfy the KKT condition in (16), i.e., col(x∗ , λ∗ ) ∈ zer(A + B).
(ii): zer(A + B) , ∅ and zer(Ā + B̄) , ∅.

The next result shows that Algorithm 3.5, or equivalently
Algorithm 4.1, can be regarded as a forward-backward operator
splitting method for finding zeros of a sum of operators, or
an iterative computation of fixed points of a composition of
operators.
Lemma 4.3. Suppose Φ in (29) is positive definite, and
operators Ā and B̄ in (30) are maximally monotone. Denote
T 1 := Id − Φ−1 Ā and T 2 := (Id + Φ−1 B̄)−1 . Then any limiting
point of Algorithm 3.5, i.e., col(x∗ , z̄∗ , λ̄∗ ), is a zero of Ā + B̄ and
is a fixed point of T 2 ◦ T 1 .

Proof: (i): By the definition of operators Ā, B̄ in (30), we
have,
Ā + B̄ : col(x, z̄, λ̄) 7→
 
 


 0 0 −ΛT   x   0   F(x) + NΩ (x)
 
 
 

0
 0 0 −L̄   z̄  +  0  + 
NRmN
(λ̄)
−b̄
λ̄
Λ L̄
L̄
+

Proof: Denote $ = col(x, z̄, λ̄), then using (28), (29) and (30),
Algorithm 3.5 can written in a compact form as follows:
− Ā($k ) ∈ B̄($k+1 ) + Φ($k+1 − $k ).

(34)
Suppose that col(x∗ , z̄∗ , λ̄∗ ) ∈ zer(Ā + B̄). From the second
line of (34),
−L̄λ̄∗ = −L ⊗ Im λ̄∗ = 0mN .

(31)

Since Φ is symmetric and positive definite, we can write
equation (31) as $k − Φ−1 Ā($k ) ∈ $k+1 + Φ−1 B̄($k+1 ), or
equivalently,
(Id − Φ−1 Ā)($k ) ∈ (Id + Φ−1 B̄)($k+1 ).

It follows that λ̄∗ = 1N ⊗ λ∗ , λ∗ ∈ Rm since L is the weighted
Laplacian of multiplier graph Gλ and Gλ is connected due to
Assumption 3.
Then by the first line of (34), combined with ΛT =
diag{AT1 , · · · , ATN } and λ∗1 = λ∗2 =, · · · , = λ∗N = λ∗ , we have

(32)

Since Φ−1 B̄ is maximally monotone (refer to Lemma 5.6), (Id+
Φ−1 B̄)−1 is single-valued. Then by the definition of the inverse
of a set-valued operator, Algorithm 3.5 is written as
$k+1 = (Id + Φ B̄) (Id − Φ Ā)$k := T 2 ◦ T 1 $k .
−1

−1

−1






0n ∈ −ΛT 1N ⊗ λ∗ + F(x∗ ) + NΩ (x∗ ),

(33)

(35)

or equivalently,

Suppose that Algorithm 3.5, or equivalently (33), converges
to a limiting point $∗ . Then by the continuity of the right
hand of Algorithm 3.5 (In fact, the right hand of Algorithm
3.5 is Lipschitz continuous due to Assumption 2 and the
nonexpansive property of projection operator), $∗ = T 2 T 1 $∗ .
Therefore, any limiting point of Algorithm 3.5 is a fixed point
of the composition T 2 ◦ T 1 , and Algorithm 3.5 can be regarded
as an iterative computation of fixed points of T 2 ◦ T 1 .
By Theorem 25.8 of Bauschke & Combettes (2011), (33) is
also the forward-backward splitting algorithm for finding zeros
of a sum of monotone operators, hence $∗ ∈ zer(Φ−1 Ā+Φ−1 B̄)
for any limiting point $∗ . Since Φ is positive definite, any
limiting point $∗ = col(x∗ , z̄∗ , λ̄∗ ), i.e., $∗ = T 2 T 1 $∗ , also
belongs to zer(Ā + B̄). In fact,

0ni ∈ −ATi λ∗ + ∇ xi fi (xi∗ , x∗−i ) + NΩi (xi∗ ), i = 1, · · · , N.

(36)

By the third line of (34) and using L̄λ̄∗ = 0, it follows that
0mN ∈ Λx∗ − b̄ + L̄z̄∗ + NRmN
(λ̄∗ ).
+
This implies that there exist v1 , v2 , · · · , vN ∈ NRm+ (λ∗ ), such that
0mN = Λx∗ − b̄ + L ⊗ Im z̄∗ + col(v1 , · · · , vN ).
Multiplying both sides of above equation with 1TN ⊗ Im and
combining with 1T L = 0T , we have
0m = (1TN ⊗ Im )(Λx∗ − b̄ + L ⊗ Im z̄∗ + col(v1 , · · · , vN ))
=

$∗ = (Id + Φ−1 B̄)−1 (Id − Φ−1 Ā)$∗
⇔ (Id − Φ−1 Ā)$∗ ∈ (Id + Φ−1 B̄)$∗
⇔ 0 ∈ Φ−1 (Ā + B̄)($∗ ) ⇔ 0 ∈ (Ā + B̄)($∗ ).

N
X
i=1

Ai xi∗ −

N
X
i=1

bi +

N
X

vi

i=1

PN
TN
11
N
Due to the fact that NTi=1
,
Ωi =
i=1 NΩi if
i=1 int(Ωi ) , ∅
we have
PN
PN
PN
0m ∈
Ai xi∗ − i=1
bi + i=1
NRm+ (λ∗ )
Pi=1
P
N
N
∗
T
N
∈
Ai xi − i=1 bi + N i=1 Rm+ (λ∗ )
(37)
Pi=1
PN
N
∗
∗
m (λ ).
∈
A
x
−
b
+
N
i
i
R
i=1
i=1
+
i

2
Remark 4.4. The iteration $k+1 = T 2 T 1 $k is also known as
Picard iteration for iteratively approximating fixed points of
T 2 T 1 (refer to Berinde (2007)). Lemma 5.5 will give a sufficient
condition for Φ to be positive definite. Lemma 5.4 will give the
condition that ensures Ā and B̄ to be maximally monotone.

11 Corollary

9

16.39 of Bauschke & Combettes (2011)

By (36) and (37), for any col(x∗ , λ̄∗ , z̄∗ ) ∈ zer(Ā + B̄), the
KKT condition for V I(F, X) in (14), i.e. (16), is satisfied for
x∗ , λ∗ . We conclude that x∗ solves V I(F, X) in (14), and is a
variational GNE of game (7) by Theorem 3.2. It also follows
that λ∗ together with x∗ satisfy the KKT condition in (16). This
also implies col(x∗ , λ∗ ) ∈ zer(A + B).

analysis on the contractive property of underlying operators.
However, the contractivity assumption on operators is usually
quite restrictive. Herein we resort to the theory of averaged
operators and firmly nonexpansive operators for convergence
analysis. Firstly we give some basic definitions and properties
of averaged operators and firmly nonexpansive operators13 . All
the following results are valid in Hilbert spaces, thus they hold
in Euclidean spaces with any G−matrix induced norm || · ||G ,
given G as a symmetric positive definite matrix. Denote by
|| · || an arbitrary matrix induced norm in a finite dimensional
Euclidean space.
An operator T : Ω ⊂ Rm → Rm is nonexpansive if it
is 1−Lipschitzian, i.e., ||T (x) − T (y)|| ≤ ||x − y||, ∀x, y ∈ Ω.
An operator T is α−averaged if there exists a nonexpansive
0
0
operator T such that T = (1 − α)Id + αT . Denote the class
of α−averaged operators as A(α). If T ∈ A( 12 ), then T is also
called firmly nonexpansive operator.

(ii) Under Assumptions 1 and 2, the considered game in (7)
has a unique variational GNE x∗ , and there exists λ∗ ∈ Rm
such that the KKT condition (16) is satisfied, i.e. col(x∗ , λ∗ ) ∈
zer(A + B). Therefore zer(A + B) , ∅.
Then we need to show that there exists col(x∗ , λ̄∗ , z̄∗ ) such
that col(x∗ , λ̄∗ , z̄∗ ) ∈ zer(Ā + B̄).
Take λ̄∗ = 1N ⊗ λ∗ . Because L1N = 0, the second line of (34)
is satisfied.
Since col(x∗ , λ∗ ) ∈ zer(A + B), 0 ∈ F(x∗ ) − AT λ∗ +
NΩ (x∗ ). Using λ∗1 = λ∗2 =, · · · , = λ∗N = λ∗ , ΛT λ̄∗ =
col(AT1 λ∗ , · · · , ATN λ∗ ) = AT λ∗ . Therefore, the first line of (34) is
satisfied with x∗ , 1N ⊗ λ∗ .
Moreover, with col(x∗ , λ∗ ) ∈ zer(A + B), 0m ∈ Ax∗ −
b + NRm+ (λ∗ ). Then we need to show that there exists z̄∗ =
col(z∗1 , · · · , z∗N ) ∈ RmN , such that the third line of (34) is
satisfied. Take v∗ ∈ NRm+ (λ∗ ) such that 0 = Ax∗ − b + v∗ . Since
QN
λ̄∗ = 1N ⊗ λ∗ and NRmN
NRm+ (λ∗ ), take v∗1 = v∗2 =
(λ̄∗ ) = i=1
+
, · · · , = v∗N = N1 v∗ ∈ NRm+ (λ∗ ). Then (1TN ⊗ Im )(Λx∗ − b̄ + L̄λ̄∗ +
PN
PN
bi + v∗ = Ax∗ − b + v∗ = 0m .
Ai xi∗ − i=1
col(v∗1 , · · · , v∗N )) = i=1
∗
∗
∗
That is Λx − b̄ + L̄λ̄ + col(v1 , · · · , v∗N ) ∈ Null(1TN ⊗ Im ). By
the fundamental theorem of linear algebra12 , Null(1TN ⊗ Im ) =
Range(1N ⊗ Im )⊥ and Range(L̄) = Null(L̄)⊥ since L̄ is also
symmetric. Notice that Null(L̄) = Range(1N ⊗ Im ), hence
Λx∗ − b̄ + L̄λ̄∗ + col(v∗1 , · · · , v∗N ) ∈ Range(L̄). Noticing that
col(v∗1 , · · · , v∗N ) ∈ NRmN
(λ̄∗ ), there exists z̄∗ ∈ RmN such that the
+
third line of (34) is satisfied with x∗ , z̄∗ , 1N ⊗ λ∗ . Therefore,
zer(Ā + B̄) , ∅.
2

Lemma 5.1. 14 Given an operator T : Ω ⊂ Rm → Rm and
α ∈ (0, 1), the following three statements are equivalent:
(i): T ∈ A(α);
2
(ii): ||T x−T y||2 ≤ ||x−y||2 − 1−α
α ||(x−y)−(T x−T y)|| , ∀x, y ∈ Ω
2
2
(iii): ||T x − T y|| + (1 − 2α)||x − y|| ≤ 2(1 − α)hx − y, T x −
T yi, ∀x, y ∈ Ω.
By (iii) of Lemma (5.1), T ∈ A( 12 ) if and only if
||T (x) − T (y)||2 ≤ hx − y, T (x) − T (y)i, ∀x, y ∈ Ω.

The operator T is called β−cocoercive (or β−inverse strongly
monotone) if βT is firmly nonexpansive, i.e.,
β||T (x) − T (y)||2 ≤ hx − y, T (x) − T (y)i, ∀x, y ∈ Ω.

(39)

Lemma 5.2. 15 For a convex differentiable function f with ϑ−
Lipschitzian gradient, we have ∇ f to be ϑ1 −cocoercive, i.e.,
1
||∇ f (x) − ∇ f (y)||22 ≤ hx − y, ∇ f (x) − ∇ f (y)i.
(40)
ϑ
Lemma 5.2 is known as Baillon-Haddad theorem, and one
elementary proof can be found in Theorem 2.1.5 of Nesterov
(2013).

5. Convergence Analysis
In this section, we prove the convergence of Algorithm 3.5
by giving a sufficient step-size choice condition. The analysis
is based on the compact reformulation (33). We will first
show that all the prerequisites in Lemma 4.3 can be satisfied
under suitable step-sizes. Then (33), i.e., Algorithm 3.5,
can be regarded as a forward-backward splitting algorithm for
finding zeros of a sum of monotone operators, or equivalently,
an iterative computation of fixed points of a composition of
operators.
In fact, some existing NE (GNE) algorithms can also be
regarded as a type of iterative computation of fixed points
of operators, such as the best-response learning dynamics
(Parise, Gentile, Grammatico & Lygeros (2015)), relaxation
algorithms based on Nikaido-Isoda function (Krawczyk &
Uryasev (2000) and Contreras, Klusch, Krawczyk (2004))
and the proximal-best response algorithm in Facchinei & Pang
(2010). Most of above works built their algorithm convergence
12 Page

(38)

Lemma 5.3. 16 If operator ∆ is maximally monotone, then
T = R∆ = (Id + ∆)−1 is firmly nonexpasive and 2R∆ − Id is
nonexpansive.
Hence, the projection operator onto a closed convex set is firmly
nonexpansive since PΩ = ProxιΩ = R∂ιΩ = RNΩ and NΩ is
maximally monotone17 .
In the following, we analyze the maximal monotonicity of
operators Ā, B̄ in (30), the positive definite property of matrix
Φ, and the properties of operators T 1 and T 2 defined in Lemma
4.3 by giving sufficient step-sizes choice conditions, which are
shown in Lemma 5.4, 5.5 and 5.6, respectively.
13 Chapter

4 and Chapter 20 of Bauschke & Combettes (2011)
4.25 of Bauschke & Combettes (2011)
15 Theorem 18.15 of Bauschke & Combettes (2011)
16 Proposition 23.7 of Bauschke & Combettes (2011)
17 Example 20.41 and Proposition 4.8 of Bauschke & Combettes (2011)
14 Proposition

405 of Meyer (2000)

10

Obviously, domB1 = Rn+2mN , hence B̄ = B1 + B2 is also
maximally monotone21 .
2

Lemma 5.4. Suppose Assumptions 1-3 hold.
Given an
Euclidean space with norm || · ||2 , then
(i): Operator Ā in (30) is β−cocoercive with 0 < β ≤
min{ 2d1 ∗ , θη2 } where d∗ is the maximal weighted degree of
P
P
multiplier graph Gλ , i.e., d∗ = max{ Nj=1 w1 j , · · · , Nj=1 wN j },
and η, θ are parameters in Assumption 2;
(ii): Operator B̄ in (30) is maximally monotone.

Lemma 5.5. Given any δ > 0, if each player i takes τi > 0,
νi > 0 and σi > 0 as its local fixed step-sizes in Algorithm 3.5
that satisfy:

Proof: (i): According to the definition of Ā in (30) and the
definition of β−cocercive in (39), we need to prove that

 
 
 

 F(x1 )   F(x2 )   x1   x2 

 
 
 

0
0
h
 − 
 ,  z̄1  −  z̄2 i ≥


λ̄2
L̄λ̄1 − b̄  L̄λ̄2 − b̄  λ̄1 
 x1   x2 
 F(x1 )   F(x2 ) 

 
 2 
 

0
0
β|| 
 ||2 , ∀  z̄1  ,  z̄2  ∈ Ω̄
 − 

λ̄2
λ̄1
L̄λ̄2 − b̄
L̄λ̄1 − b̄
(41)
Notice that L̄λ̄− b̄ is the gradient of function f˜(λ̄) := 21 λ̄T L̄λ̄−
b̄T λ̄. Moveover, f˜(λ̄) is a convex function since ∇2 f˜(λ̄) = L̄ is
positive semi-definite due to Assumption 3 18 . It can easily be
verified that L̄λ̄− b̄ is ||L||2 −Lipschitz continuous (notice that the
eigenvalues of L̄ are just the elements in col(0, s2 · · · , sN ) ⊗ 1m ),
therefore, by Lemma 5.2
hL̄λ̄1 − b̄ − (L̄λ̄2 − b̄), λ̄1 − λ̄2 i ≥

1
||L̄λ̄1 − L̄λ̄2 ||22 .
||L||2

τi

≤

νi
σi

≤
≤

P1
,
T
max j=1,...ni { m
k=1 |[Ai ] jk |}+δ
1
2di +δ
1
Pn i
max j=1,...m { k=1
|[Ai ] jk |}+2di +δ

(44)

then matrix Φ defined in (29) is positive definite, and Φ−δIn+2mN
is positive semi-definite.
Proof: It is sufficient to show that Φ − δIn+2mN is positive
semi-definite.

 −1

 τ̄ − δIn
0
ΛT


−1
0
ν̄ − δImN
L̄
Φ − δIn+2mN = 
 ,

−1
Λ
L̄
σ̄ − δImN
(45)
where τ̄−1 − δIn = diag{( τ11 − δ)In1 , · · · , ( τ1N − δ)InN }, ν̄−1 −
δImN = diag{( ν11 − δ)Im , · · · , ( ν1N − δ)Im }, and σ̄−1 − δImN =
diag{( σ11 − δ)Im , · · · , ( σ1N − δ)Im }. One sufficient condition
for matrix Φ − δIn+2mN to be positive semi-definite is that it
is diagonally dominant with nonnegative diagonally elements,
that is for every row of the matrix the diagonal entry is larger
than or equal to the sum of the magnitudes of all the other
(non-diagonal) entries in that row. This is equivalent to require
that,
P
1
T
≥ max j=1,...ni { m
k=1 |[Ai ] jk |}
τi − δ
P
m
1
≥
2di
j=1 |Li j | =P
νi − δ
P
(46)
ni
1
−
δ
≥
max
{
|[Ai ] jk |} + mj=1 |Li j |
j=1,...m
σi
Pk=1
i
= max j=1,...m { nk=1
|[Ai ] jk |} + 2di

(42)

Since ||L||2 ≤ sN and d∗ ≤ sN ≤ 2d∗ by (6) where d∗ =
max{d1 , · · · , dN } is the maximal weighted degree of multiplier
1
graph Gλ , we have ||L||
≥ 2d1 ∗ .
2
Meanwhile by Assumption 2, F(x) is η−strongly monotone
and θ− Lipschitz continuous over Ω. By ||F(x1 ) − F(x2 )||22 ≤
θ2 ||x1 − x2 ||22 , ∀ x1 , x2 ∈ Ω we have
η
||F(x1 ) − F(x2 )||22 .
θ2
(43)
Taking 0 < β ≤ min{ 2d1 ∗ , θη2 } and adding (42) and (43) yields
(41). Thus operator Ā is β−cocoercive. By the definition
of β−cocoercive in (39), operator Ā is also monotone. Since
operator Ā is also single-valued, it is also maximally monotone.
(ii): The operator B̄ in (30) can be written as:


 

 0 0 −ΛT   x   NΩ (x) 
 0 0 −L̄   z̄  +  0  := B1 + B2 .

 
 

NRmN
Λ L̄
0
λ̄
+
hF(x1 ) − F(x2 ), x1 − x2 i ≥ η||x1 − x2 ||22 ≥

It can easily be verified that if each agent chooses its local
step-sizes satisfying (44), then (46) is satisfied.
2
Given a globally known parameter δ, each agent can
independently choose its local step sizes τi ,νi , and σi with the
rule given in (44).
Suppose that the step-sizes τi , µi , σi in Algorithm 3.5 are
chosen such that Φ in (29) is positive definite. Thus we can
define a norm
√ induced by√matrix Φ for a finite Euclidean space
as ||x||Φ = hx, xiΦ = hΦx, xi. The next result investigates
the properties of operators Φ−1 Ā, Φ−1 B̄ and T 1 , T 2 defined in
Lemma 4.3 under Φ−induced norm || · ||Φ .

Since L̄ is symmetric, B1 is a skew-symmetric matrix, i.e.,
BT1 = −B1 . Hence, B1 is maximally monotone19 .
B2 can be written as the direct sum of NΩ × 0mN × NRmN
.
+
Both NΩ and NRmN
are
maximally
monotone
as
normal
cones
+
of closed convex sets. Obviously, 0mN is also maximally
monotone as a single-valued operator. Furthermore, the direct
sum of maximally monotone operators is also maximally
monotone20 , hence B2 is also maximally monotone.

Lemma 5.6. Suppose Assumptions 1-3 hold. Take 0 < β ≤
min{ 2d1 ∗ , θη2 } where d∗ is the maximal weighted degree of
multiplier graph Gλ , and η, θ are parameters in Assumption
1
2. Take δ > 2β
. Suppose that the step-sizes τi , νi , σi in
Algorithm 3.5 are chosen to satisfy (44). Then the operators
Φ−1 Ā and Φ−1 B̄ with Φ in (29) and Ā, B̄ in (30), and operators

18 Proposition

17.10 of Bauschke & Combettes (2011)
20.30 of Bauschke & Combettes (2011)
20 Proposition 20.23 of Bauschke & Combettes (2011)
19 Example

21

11

Corollary 24.4 of Bauschke & Combettes (2011)

T 1 = Id − Φ−1 Ā, T 2 = RΦ−1 B̄ = (Id + Φ−1 B̄)−1 defined as Lemma
4.3 satisfy the following properties under the Φ−induced norm
|| · ||Φ :
1
(i). Φ−1 Ā is βδ−cocoercive, and T 1 ∈ A( 2δβ
).
−1
(ii). Φ B̄ is maximally monotone, and T 2 ∈ A( 12 ).

1
A( 2δβ
), and T 2 := (Id + Φ−1 B̄)−1 ∈ A( 12 ) in the Euclidean
space with Φ−matrix induced norm || · ||Φ . The next result
shows the convergence of Algorithm 3.5 based on its compact
reformulation (33) and the properties of T 1 and T 2 .

Proof: (i): By the definition of cocoercivity in (39), we need to
prove hΦ−1 Ā(x) − Φ−1 Ā(y), x − yiΦ ≥ βδ||Φ−1 Ā(x) − Φ−1 Ā(y)||2Φ ,
∀x, y ∈ Ω̄. Noticing that by the choice of parameters τi , νi , σi ,
we have that matrix Φ − δIn+2mN is positive semi-definite from
Lemma 5.5. Denote smax (Φ) and smin (Φ) as the maximal and
minimal eigenvalues of matrix Φ. It must hold that smax (Φ) ≥
smin (Φ) ≥ δ. Furthermore, ||Φ||2 = smax (Φ) ≥ smin (Φ) =
1
≥ δ 22 , therefore, we also have ||Φ−1 ||2 ≤ 1δ . Notice that
||Φ−1 ||2
the operator Ā is single-valued and Φ−1 is also nonsingular, so
that for any x, y ∈ Ω̄,

Theorem 5.7. Suppose Assumptions 1-3 hold. Take 0 <
β ≤ min{ 2d1 ∗ , θη2 } where d∗ is the maximal weighted degree of
multiplier graph Gλ , and η, θ are parameters in Assumption
1
. The step-sizes τi , νi , σi in Algorithm 3.5
2. Take δ > 2β
are chosen to satisfy (44). Then with Algorithm 3.5, each
player has its local strategy xi,k converging to its corresponding
component in the variational GNE of game (7), and the local
Lagrangian multipliers λi,k of all the agents converge to the
same Lagrangian multiplier corresponding with KKT condition
(16), i.e.,

||Φ−1 Ā(x) − Φ−1 Ā(y)||2Φ
= hΦΦ−1 (Ā(x) − Ā(y)), Φ−1 (Ā(x) − Ā(y))i
= ||Ā(x) − Ā(y)||2Φ−1 ≤ 1δ ||Ā(x) − Ā(y)||22

lim xi,k → xi∗ ,

k→∞

By the β−cocoercive property of Ā in Lemma 5.4 and the above
inequality,

Firstly, by (i) and (ii) of Lemma 5.1 and the fact that T 1 , T 2
are averaged operators due to Lemma 5.6, T 1 , T 2 are also
nonexpansive operators under the Φ−matrix induced norm ||·||Φ .
Take any $∗ ∈ zer(Ā + B̄), or equivalently any fixed point of
T 2 ◦ T 1 , i.e., $∗ = T 2 T 1 $∗ , and then by Lemma 4.3 and (33),

1
1
1
)Id +
(−T̆ ) ∈ A(
)
2βδ
2βδ
2βδ

||$k+1 − $∗ ||Φ = ||T 2 (T 1 $k ) − T 2 (T 1 $∗ )||Φ
≤ ||T 1 $k − T 1 $∗ || ≤ ||$k − $∗ ||Φ

1
since 1 < 2βδ by the assumption that δ > 2β
and −T̆ is also
nonexpansive.
(ii). Φ is symmetric positive definite and nonsingular. For
any (x, u) ∈ graΦ−1 B̄ and (y, v) ∈ graΦ−1 B̄, Φu ∈ ΦΦ−1 B̄(x) ∈
B̄(x) and Φv ∈ ΦΦ−1 B̄(y) ∈ B̄(y). Then hx − y, u − viΦ =
hx − y, Φ(u − v)i ≥ 0, ∀x, y since B̄ is monotone by Lemma
5.4. Therefore, Φ−1 B̄ is monotone under the Φ−matrix induced
product h·, ·iΦ .
Furthermore, take (y, v) ∈ Ω̄ × Rn+2mN , and hx − y, u − viΦ ≥
0, for any other (x, u) ∈ gra(Φ−1 B̄). For any (x, ũ) ∈ graB̄,
we have (x, Φ−1 ũ) ∈ gra(Φ−1 B̄). hx − y, Φ(Φ−1 ũ − v)i ≥ 0,
or equivalently, hx − y, ũ − Φv)i ≥ 0. Since B̄ is maximally
monotone, then (y, Φv) ∈ graB̄. We conclude that v ∈ Φ−1 B̄(y)
which implies that Φ−1 B̄ is maximally monotone.
Therefore, by Lemma 5.3 T 2 = (Id + Φ−1 B̄)−1 is firmly
nonexpansive under the Φ−matrix induced norm || · ||Φ .
2

(49)

Hence the sequence {||$k − $∗ ||Φ } is non-increasing and
bounded from below. By the monotonic convergence theorem,
{||$k − $∗ ||} is bounded and converges for every $∗ ∈ zer(A +
B).
1
) and T 2 ∈ A( 12 ). Denote ξ =
By Lemma 5.6, T 1 ∈ A( 2δβ
1
2δβ ∈ (0, 1). Then with (ii) of Lemma 5.1 and (33) we have,

||$k+1 − $∗ ||2Φ = ||T 2 (T 1 $k ) − T 2 (T 1 $∗ )||2Φ
≤ ||T 1 $k − T 1 $∗ ||2Φ
−||(T 1 $k − T 1 $∗ ) − (T 2 T 1 $k − T 2 T 1 $∗ )||2Φ
≤ ||$k − $∗ ||2Φ − ||(T 1 $k − T 1 $∗ ) − (T 2 T 1 $k − T 2 T 1 $∗ )||2Φ
∗
∗ 2
− 1−ξ
ξ ||$k − $ − (T 1 $k − T 1 $ )||Φ
(50)
where the first inequality follows by T 2 ∈ A( 21 ) and the second
1
inequality follows by T 1 ∈ A( 2δβ
), both utilizing (ii) of Lemma
5.1. Notice that

Summarizing the above results, take 0 < β ≤ min{ 2d1 ∗ , θη2 },
1
δ > 2β
, and τi , νi , σi satisfying (44), then T 1 := Id − Φ−1 Ā ∈
22 Proposition

(48)

Proof: With Lemma 5.4 and 5.5, Algorithm 3.5 can be written
in a compact form (33) according to Lemma 4.3, i.e., $k+1 =
T 2 T 1 $k . The convergence analysis will be conducted via the
analysis of this iterative computation of fixed points of T 2 ◦ T 1 .

hΦ−1 Ā(x) − Φ−1 Ā(y), x − yiΦ = hĀ(x) − Ā(y), x − yi
≥ β||Ā(x) − Ā(y)||22 ≥ βδ||Φ−1 Ā(x) − Φ−1 Ā(y)||2Φ .
(47)
Therefore, the operator Φ−1 Ā is βδ−cocoercive under the
Φ−induced norm || · ||Φ .
Moreover, βδΦ−1 Ā is firmly nonexpansive by the definition
of cocoercive operator. This implies that there exists a
nonexpansive operator T̆ such that βδΦ−1 Ā = 12 T̆ + 21 Id. Then
T 1 = Id − Φ−1 Ā = (1 −

lim λi,k → λ∗ , ∀i = 1, · · · , N.

k→∞

α||x||2 + (1 − α)||y||2 = ||αx + (1 − α)y||2 + α(1 − α)||x − y||2 . (51)

5.2.7 and 5.2.8 in Meyer (2000)

12

For the second and third terms on the right hand side of (50),

the current operator output with previous iterate, leads to
the well-known Krasnosel’skiĭ-Mann type of fixed point
iteration23 , and has been utilized in (generalized) Nash
equilibrium computation in Contreras, Klusch, Krawczyk
(2004) and Krawczyk & Uryasev (2000). Meanwhile, inertial
algorithms in operator splitting methods have received attention
in recent years, such as Alvarez & Attouch (2001), Attouch,
Chbani, Peypouquet, and Redont (2016), Lorenz & Pock
(2015) and Rosasco, Villa & Vũ (2016). These efforts
are partially motivated by the heavy ball method in Polyak
(1987) and Nesterov’s acceleration algorithm in Nesterov
(2013)) for optimization problems and their recent success
in machine learning applications (refer to Wibisono, Wilson,
and Jordan, (2016)). In particular, Nesterov’s acceleration
algorithm is proved to enjoy an optimal convergence speed with
a specific step-size choice. Thereby, in this work we consider a
distributed algorithm with inertia for variational GNE seeking
given as below:

||(T 1 $k − T 1 $∗ ) − (T 2 T 1 $k − T 2 T 1 $∗ )||2Φ
∗
∗ 2
+ 1−ξ
ξ ||$k − $ − (T 1 $k − T 1 $ )||Φ
1
∗
= ξ [ξ||(T 1 $k − T 1 $ ) − (T 2 T 1 $k − T 2 T 1 $∗ )||2Φ
+(1 − ξ)||(T 1 $k − T 1 $∗ ) − ($k − $∗ )||2Φ ]
1
= ξ ||(T 1 $k − T 1 $∗ ) − ξ(T 2 T 1 $k − T 2 T 1 $∗ )
−(1 − ξ)($k − $∗ )||2Φ
1
+ ξ ξ(1 − ξ)||(T 1 $k − T 1 $∗ ) − (T 2 T 1 $k − T 2 T 1 $∗ )
−(T 1 $k − T 1 $∗ ) + ($k − $∗ )||2Φ
≥ (1 − ξ)||($k − $∗ ) − (T 2 T 1 $k − T 2 T 1 $∗ )||2Φ
= (1 − ξ)||$k − T 2 T 1 $k ||2Φ
= (1 − ξ)||$k − $k+1 ||2Φ
(52)
where the second equality follows from (51) by setting α = ξ,
x = (T 1 $k − T 1 $∗ ) − (T 2 T 1 $k − T 2 T 1 $∗ ) and y = (T 1 $k −
T 1 $∗ ) − ($k − $∗ ).
Combining (50) and (52) yields ∀k ≥ 0,
||$k+1 − $∗ ||2Φ ≤ ||$k − $∗ ||2Φ − (1 − ξ)||$k − $k+1 ||2Φ .

Algorithm 6.1.

(53)

Using (53) from 0 to k and adding all k + 1 inequalities yields
||$k+1 − $∗ ||2Φ ≤ ||$0 − $∗ ||2Φ − (1 − ξ)

k
X

Acceleration phase

||$l − $l+1 ||2Φ .

l=0

P
Taking limit as k → ∞ we have, (1 − ξ) ∞
||$ − $k+1 ||2Φ ≤
k=1
P∞ k
∗ 2
||$0 − $ ||Φ . Since 1 − ξ > 0, it follows that k=1 ||$k − $k+1 ||2Φ
converges and limk→∞ $k − $k+1 = 0.
Since {||$k − $∗ ||} is bounded and converges, {$k } is a
bounded sequence. There exists a subsequence{$nk } that
converges to $̃∗ . Notice that the composition T 2 ◦ T 1 is
(Lipschitz) continuous and single-valued, because (33) is just
an equivalent expression of Algorithm 3.5, and obviously the
right hand side of Algorithm 3.5 is continuous. $nk +1 =
T 2 T 1 $nk . Since T 2 T 1 is continuous, and limnk →∞ $nk − $nk +1 =
0, passing to limiting point, we have $̃∗ = T 2 T 1 $̃∗ . Therefore,
the limiting point $̃∗ is a fixed point of T 2 T 1 , or equivalently,
$̃∗ ∈ zer(Ā + B̄).
Setting $∗ = $̃∗ in (49), we have {||$k − $̃∗ ||} is bounded
and converges. Since there exists a subsequence {$nk } that
converges to $̃∗ , it follows that {||$k − $̃∗ ||} converges to zero.
Therefore, limk→∞ $k → $̃∗ . By Theorem 4.5, this just implies
(48).
2

x̃i,k = xi,k + α(xi,k − xi,k−1 )

(54)

z̃i,k = zi,k + α(zi,k − zi,k−1 )

(55)

λ̃i,k = λi,k + α(λi,k − λi,k−1 )

(56)

U pdate phase


xi,k+1 = PΩi x̃i,k − τi (∇ xi fi ( x̃i,k , { x̃ j,k } j∈N f ) − ATi λ̃i,k )
i
X
zi,k+1 = z̃i,k + νi
wi j (λ̃i,k − λ̃ j,k )

(57)
(58)

j∈Niλ

X
n

λi,k+1 = PRm+ λ̃i,k − σi 2Ai xi,k+1 − Ai x̃i,k − bi +
wi j (λ̃i,k − λ̃ j,k )
j∈Niλ

+2

X
j∈Niλ

wi j (zi,k+1 − z j,k+1 ) −

X

wi j (z̃i,k − z̃ j,k )

o

(59)

j∈Niλ

α > 0 is a fixed step-size in the acceleration phase, and
τi > 0, νi > 0, σi > 0 are fixed step-sizes of player i, and
W = [wi j ] is the weighted adjacency matrix of multiplier graph
Gλ .
Compared with Algorithm 3.5, Algorithm 6.1 has two
phases. In the acceleration phase, each player uses the local
state information of the last two steps to get predictive variables
by a simple linear extrapolation. In the update phase, the
players just feed the predictive variables to Algorithm 3.5 to
get the next iterates. Hence compared with Algorithm 3.5,
Algorithm 6.1 has only an additional simple local computation
burden. Obviously, Algorithm 6.1 is also totally distributed,
and shares all the features of Algorithm 3.5. However, there is
an additional need to choose a proper step-size α.
In the following two subsections, we will first give some
intuitive interpretation of Algorithm 6.1 from the viewpoint of a

6. Distributed algorithm with inertia
In this section, we propose a distributed algorithm with
inertia for variational GNE seeking, which possibly accelerates
the convergence under some mild additional computation
burden.
There are various modifications of Picard fixed point iteration
to achieve the possible acceleration of convergence speed, and
most of them fall into the domains of relaxation algorithm
and inertial algorithm (Refer to Iutzeler & Hendrickx
(2016) for reviews and numerical comparisons for optimization
problems). The relaxation algorithm that simply combines

23 Chapter

13

5 of Bauschke & Combettes (2011)

where $̃(kh) is an interpolation point to be determined later.
Denote $k = $(kh) and take h = 1, then (66) can be written as

discretization of continuous-time dynamical systems, and then
prove its convergence.

$k +(1−α̃)($k −$k−1 )−Φ−1 Ā($̃k ) ∈ $k+1 +Φ−1 B̄($k+1 ). (67)

6.1. Interpretations from viewpoints of dynamical systems

Denote α = 1 − α̃ and take $̃k = $k + α($k − $k−1 ), then (67)
can be written as

The interpretation of inertial (acceleration) algorithms from
a continuous-time dynamical system viewpoint can be found
in Polyak (1987) and most recently in Wibisono, Wilson,
and Jordan, (2016) for optimization problems and in Attouch,
Chbani, Peypouquet, and Redont (2016) for proximal point
algorithms. Here we give a comparative development of
Algorithm 3.5 and Algorithm 6.1 just for illustrations of the
differences behind the algorithms.
Firstly, let us show that Algorithm 3.5, or equivalently its
compact reformulation (33) in Lemma 4.3, can be interpreted
as the discretization of the following dynamical system:
$̇ ∈ −Φ−1 Ā($) − Φ−1 B̄($).

$̃k = $k + α($k − $k−1 ),
$̃k − Φ−1 Ā($̃k ) ∈ $k+1 + Φ−1 B̄($k+1 ).

(68) leads to equations (63)-(64), or equivalently Algorithm 6.1.
Remark 6.2. Compared with (60), (65) is a second order
dynamical system with an additional inertial term α$̇, hence
(65) enjoys better convergence properties than (60). Therefore,
it is expected that Algorithm 6.1, as a discretization of (65),
would have better convergence properties than Algorithm 3.5.

(60)
6.2. Convergence analysis
The following result proves the convergence of Algorithm
6.1 by providing sufficient step-size choices for α as well as
τi , νi , σi . The sufficient choice condition for α can be ensured
by solving a simple algebra inequality. The proof idea of the
following result is motivated by inertial algorithms works for
optimization and operator splitting such as Alvarez & Attouch
(2001), Attouch, Chbani, Peypouquet, and Redont (2016),
Rosasco, Villa & Vũ (2016), Iutzeler & Hendrickx (2016),
and especially Lorenz & Pock (2015). However, since this
work considers a noncooperative game setup and adopts a fixed
step-size in the distributed algorithm, Theorem 6.3’s proof is
also provided for completeness.

In fact, for differential inclusion (60), we have the following
implicit/explicit discretization with step-size of h,
$(kh + h) − $(kh)
∈ −Φ−1 Ā($(kh))−Φ−1 B̄($(kh+h)). (61)
h
Denote $k = $(kh) and take h = 1, then (61) can be written as
$k − Φ−1 Ā($k ) ∈ $k+1 + Φ−1 B̄($k+1 ).

(62)

Therefore, the implicit/explicit discretization of (60) is exactly
(32) that leads to (33), or equivalently Algorithm 3.5.
Moreover, the explicit discretization corresponds with the
forward step, and the implicit discretization corresponds with
the backward step. That’s the reason why Algorithm 3.5 is
called a forward-backward splitting algorithm.
Adopt similar compact notations as in Section 4, and
denote x̃ = col( x̃1 , · · · , x̃N ), λ̃¯ = col(λ̃1 , · · · , λ̃N ), and z̃¯ =
col(z̃1 , · · · , z̃N ). And further denote $̃ = col(x̃, λ̃¯ , z̃¯). Then by
similar arguments as in Section 4 and using operators Ā and
B̄ defined in (30) and Φ defined in (29), Algorithm 6.1 can be
written in a compact form (assume B̄ is maximally monotone),
$̃k = $k + α($k − $k−1 )
−1

Theorem 6.3. Suppose Assumptions 1-3 hold. Take 0 <
β ≤ min{ 2d1 ∗ , θη2 } where d∗ is the maximal weighted degree of
multiplier graph Gλ , and η, θ are parameters in Assumption 2.
1
and 0 < α < 1
Given a sufficient small 0 <  < 1, take δ > 2β
in Algorithm 6.1 such that 2βδ(1 − 3α − ) ≥ (1 − α)2 . Suppose
that player i chooses its step-sizes τi , νi , σi in Algorithm 6.1
satisfying (44). Then with Algorithm 6.1, players’ local
strategies converge to the variational GNE of game in (7), and
the local multipliers λi,k of all the agents converge to the same
multiplier corresponding with KKT condition (16), i.e.,

(63)

$k+1 = (Id + Φ B̄) (Id − Φ Ā)$̃k
−1

−1

(64)

lim xi,k → xi∗ ,

k→∞

where $k is defined as in Section 4.
We can show that Algorithm 6.1, or equivalently (63)-(64),
can be interpreted as the discretization of the following
second-order continuous-time dynamical system,
$̈ + α̃$̇ ∈ −Φ−1 Ā($) − Φ−1 B̄($).

(68)

lim λi,k → λ∗ , ∀i = 1, · · · , N.

k→∞

(69)

Proof: By the choice of β, operator Ā is β−cocoercive and
operator B̄ is maximally monotone due to Lemma 5.4. By the
choice of δ, τi , νi and σi , Φ − δIn+2mN is positive semi-definite
due to Lemma 5.5, and Φ−1 Ā and Φ−1 B̄ satisfy the properties
in Lemma 5.6. Therefore, Algorithm 6.1 can be exactly written
in the compact form of (63)-(64) with similar arguments as in
Lemma 4.3.
Resorting to Theorem 4.5, we only need to show that
Algorithm 6.1 converges and its limiting point belongs to
zer(Ā + B̄). In fact, any limiting point of Algorithm 6.1 satisfies
$∗ ∈ zer(Ā + B̄) as shown next. Suppose limk→∞ $k → $∗ ,
then $̃k → $∗ and $∗ = (Id+Φ−1 B̄)−1 (Id−Φ−1 Ā)$∗ using the

(65)

In fact, for differential inclusion (65) consider the following
type of implicit/explicit discretization,
$(kh) − $(kh − h)
$(kh + h) − 2$(kh) + $(kh − h)
+ α̃
∈
2
h
h
− Φ−1 Ā($̃(kh)) − Φ−1 B̄($(kh + h)),
(66)
14

continuity of the right hand of (64). Therefore, $∗ ∈ zer(Ā+ B̄)
because Φ is a positive definite matrix.
The following relationship (similar with the cosine rule) will
be heavily utilized in the convergence analysis.
||a −

c||2Q

− ||b −

c||2Q

= 2ha − b, a − ciQ − ||a −

b||2Q ,

Using (76) for the second term on the right hand of (72)
yields
||$k − $∗ ||2Φ − ||$k+1 − $∗ ||2Φ
≥ ||$k − $k+1 ||2Φ + 2hĀ($̃k ) − Ā($∗ ), $k+1 − $∗ i
−2αh$k − $k−1 , $k+1 − $∗ iΦ

(70)

By Lemma 5.4, Ā is β−cocoercive. For the second term on
the right hand of (77), we have

which can be verified by directly expanding with ||a + b||2Q =
||a||2Q + 2ha, biQ + ||b||2Q .
The proof is divided into three parts:
• Part 1: Given any $∗ ∈ zer(Ā + B̄), ||$k+1 − $∗ ||2Φ − ||$k −
$∗ ||2Φ follows a recursive inequality as follows,
||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ
≤ α(||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ ) + (α + α2 )||$k − $k−1 ||2Φ
(71)
P
2
• Part 2: Given any $∗ ∈ zer(Ā+ B̄), ∞
k ||$k+1 −$k ||Φ < ∞
and limk→∞ $k+1 − $k = 0.

or equivalently we have
||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ ≤ −||$k − $k+1 ||2Φ
1
||$k+1 − $̃k ||22 + 2αh$k − $k−1 , $k+1 − $∗ iΦ .
+ 2β

Part 1: Given any point $∗ ∈ zer(Ā + B̄), we first prove a
recursive inequality (71) for ||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ .
Using (70) to expand the left hand of (71) yields,

α(||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ )
= α(2h$k − $k−1 , $k − $∗ iΦ − ||$k − $k−1 ||2Φ )

(80)

Combining (79) and (80), we have,
||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ
−α(||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ )
1
≤ −||$k − $k+1 ||2Φ + 2β
||$k+1 − $̃k ||22
+2αh$k − $k−1 , $k+1 − $∗ iΦ
−α(2h$k − $k−1 , $k − $∗ iΦ − ||$k − $k−1 ||2Φ )
1
≤ −||$k − $k+1 ||2Φ + 2β
||$k+1 − $̃k ||22
+2αh$k − $k−1 , $k+1 − $k iΦ + α||$k − $k−1 ||2Φ

(81)

Next using (63) and (70) for the first and third terms on the
right hand of (81) yields

(73)

Because $∗ ∈ zer(Ā + B̄), we also have

2αh$k − $k−1 , $k+1 − $k iΦ − ||$k − $k+1 ||2Φ
= 2αh α1 ($̃k − $k ), $k+1 − $k iΦ − ||$k − $k+1 ||2Φ
= 2h$̃k − $k , $k+1 − $k iΦ − ||$k − $k+1 ||2Φ
= ||$k − $̃k ||2Φ − ||$k+1 − $̃k ||Φ
= α2 ||$k − $k−1 ||2Φ − ||$k+1 − $̃k ||2Φ

(74)

Due to the maximal monotonicity of B̄ proved in Lemma 5.4,
hu − v, $k+1 − $∗ i ≥ 0, ∀ u ∈ B̄($k+1 ), v ∈ B̄($∗ )

(79)

Utilizing the equality (70) we also have,

||$k − $∗ ||2Φ − ||$k+1 − $∗ ||2Φ
= 2h$k − $k+1 , $k − $∗ iΦ − ||$k − $k+1 ||2Φ
= 2h$k − $k+1 , $k − $k+1 + $k+1 − $∗ iΦ − ||$k − $k+1 ||2Φ
= ||$k − $k+1 ||2Φ + 2h$k − $k+1 , $k+1 − $∗ iΦ
= ||$k − $k+1 ||2Φ + 2h$̃k − $k+1 , $k+1 − $∗ iΦ
−2αh$k − $k−1 , $k+1 − $∗ iΦ
(72)
where the last step is derived by incorporating (63).
To tackle the second term on the right hand of (72), we
proceed as follows. By (64) we have, $̃k − Φ−1 Ā($̃k ) ∈
$k+1 + Φ−1 B̄($k+1 ), or equivalently,

− Ā($∗ ) ∈ B̄($∗ ).

hĀ($̃k ) − Ā($∗ ), $k+1 − $∗ i
= hĀ($̃k ) − Ā($∗ ), $k+1 − √
$̃k + $̃k − $∗ i
≥ β||Ā($̃k ) − Ā($∗ )||22 + 2h β(Ā($̃k ) − Ā($∗ )), 2 √1 β ($k+1 − $̃k )i
1
≥ β||Ā($̃k ) − Ā($∗ )||22 − β||Ā($̃k ) − Ā($∗ )||22 − 4β
||$k+1 − $̃k ||22
1
2
= − 4β ||$k+1 − $̃k ||2
(78)
where the first inequality is obtained by the cocoercive property
(39) and the second inequality is obtained by 2ha, bi ≥ −||a||22 −
||b||22 .
Combining (77) with (78) we have
||$k − $∗ ||2Φ − ||$k+1 − $∗ ||2Φ ≥ ||$k − $k+1 ||2Φ
1
||$k+1 − $̃k ||22 − 2αh$k − $k−1 , $k+1 − $∗ iΦ ,
− 2β

• Part 3: We first show the convergence of {||$k − $∗ ||2Φ }
given any $∗ ∈ zer(Ā+ B̄), and then show the convergence
of Algorithm 6.1.

Φ($̃k − $k+1 ) − Ā($̃k ) ∈ B̄($k+1 ).

(77)

(75)

(82)

1
For the second term on the right hand of (81), 2β
||$k+1 −
1
2
2
$̃k ||2 ≤ 2βδ ||$k+1 − $̃k ||Φ , since Φ > δIn+2mN . Incorporating
this and (82) into (81)

By incorporating (73) and (74) into (75), we have
hΦ($̃k − $k+1 ) − Ā($̃k ) + Ā($∗ ), $k+1 − $∗ i ≥ 0

||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ
−α(||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ )
1
≤ −||$k+1 − $̃k ||2Φ + 2β
||$k+1 − $̃k ||22 + (α + α2 )||$k − $k−1 ||Φ
1
≤ −(1 − 2δβ
)||$k+1 − $̃k ||2Φ + (α + α2 )||$k − $k−1 ||2Φ
(83)

or equivalently,
h$̃k − $k+1 , $k+1 − $∗ iΦ − hĀ($̃k ) − Ā($∗ ), $k+1 − $∗ i ≥ 0.
(76)
15

Since δ >

1
2β ,

we derive (71).

P
2
Part 2: In this step, we will prove ∞
k ||$k+1 − $k ||Φ < ∞ and
limk→∞ $k+1 − $k = 0.
1
Denote S = Φ − 2β
In+2mN . Then S is symmetric and positive
1
definite since Φ ≥ δIn+2mN and δ > 2β
. The first inequality of
(83) can also be written as:
||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ − α(||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ )
≤ −||$k+1 − $̃k ||2S + (α + α2 )||$k − $k−1 ||2Φ .
(84)
For the first term on the right hand of (84),
−||$k+1 − $̃k ||2S
= 2h$k − $k+1 , $k − $̃k iS − ||$k − $k+1 ||2S − ||$k − $̃k ||2S
= −2h$k − $k+1 , α($k − $k−1 )iS
−||$k − $k+1 ||2S − α2 ||$k − $k−1 ||2S
≤ α||$k − $k+1 ||2S + α||$k − $k−1 ||2S
−||$k − $k+1 ||2S − α2 ||$k − $k−1 ||2S
= (α − 1)||$k − $k+1 ||2S + (α − α2 )||$k − $k−1 ||2S .
(85)
where the first equality follows from (70), the second equality
follows from (63), and the third inequality follows from
−2hx, yi ≤ ||x||2 + ||y||2 .
Denote Q = 2Φ − 1−α
2β In+2mN . Then Q is also symmetric and
positive definite, since α < 1 and Φ ≥ δIn+2mN . Combining (84)
with (85),

k

Therefore, limk→∞ $k+1 − $k = 0, and
$k−1 ||2Φ < ∞.

P∞

k=1 (α

+ α2 )||$k −

Part 3: In this part, we first show the convergence of
{||$k − $∗ ||2Φ } given any $∗ ∈ zer(Ā + B̄), and then show the
convergence of Algorithm 6.1.
Denote φk = max{0, ||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ } and ψk =
(α + α2 )||$k − $k−1 ||2Φ , and recall (71), we have φk+1 ≤ αφk + ψk .
Apply this relationship recursively,
φk+1 ≤ αk φ1 +

k−1
X

αi ψk−i .

(89)

i=0

Summing (89) from k = 1 to k = J,
P J Pk−1 i
PJ
1−α J
k=1 φk+1 ≤ φ1 1−α + k=1 i=0 α ψk−i
P J−1 i P J
J
≤ φ1 1−α
k=i+1 ψk−i
1−α + i=0 α

||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ
−α(||$k − $∗ ||2Φ − ||$k−1 − $∗ ||2Φ )
≤ (α − 1)||$k − $k+1 ||2S + (α − α2 )||$k − $k−1 ||2S
+(α + α2 )||$k − $k−1 ||2Φ
= (α − 1)||$k − $k+1 ||2S
1
+αh[(1 − α)(Φ − 2β
I) + (1 + α)Φ]$k − $k−1 , $k − $k−1 i
= (α − 1)||$k − $k+1 ||2S + α||$k − $k−1 ||2Q .
(86)
Denote µk = ||$k − $∗ ||2Φ − α||$k−1 − $∗ ||2Φ + α||$k − $k−1 ||2Q ,
then
µk+1 − µk = ||$k+1 − $∗ ||2Φ − α||$k − $∗ ||2Φ + α||$k+1 − $k ||2Q
−||$k − $∗ ||2Φ + α||$k−1 − $∗ ||2Φ − α||$k − $k−1 ||2Q
= ||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ − α(||$k − $∗ ||2Φ
−||$k−1 − $∗ ||2Φ ) + α||$k+1 − $k ||2Q − α||$k − $k−1 ||2Q
≤ (α − 1)||$k − $k+1 ||2S + α||$k − $k−1 ||2Q
+α||$k+1 − $k ||2Q − α||$k − $k−1 ||2Q
= (α − 1)||$k − $k+1 ||2S + α||$k+1 − $k ||2Q
1
= h[(α − 1)(Φ − 2β
I) + α(2Φ − 1−α
2β I)]$k+1 − $k ,
$k+1 − $k i
(87)
where the third inequality follows by (86).
Given a sufficient small 0 <  < 1, choose 0 < α < 1 and
1
δ > 2β
such that 2βδ(1 − 3α − ) ≥ (1 − α)2 , then
1
1−α
2β In+2mN ) − α(2Φ − 2β In+2mN )
2
− (1−α)
2β In+2mN ≥ Φ

−(α − 1)(Φ −
= (1 − 3α)Φ

α||$k−1 −$∗ ||2Φ ≤ ||$k −$∗ ||2Φ −α||$k−1 −$∗ ||2Φ +α||$k −$k−1 ||2Q ≤
µ1 . Therefore, ||$k − $∗ ||2Φ ≤ α||$k−1 − $∗ ||2Φ + µ1 . By Lemma
1 at Page 44 of Polyak (1987), ||$k − $∗ ||2Φ ≤ αk (||$1 − $∗ ||2Φ −
µ1
µ1
1−a ) + 1−α , and {$k } is bounded sequence.
P
We also have µk+1 − µ1 ≤ − ki=1 ||$i+1 − $i ||2Φ . Then
P
 ki=1 ||$i+1 − $i ||2Φ ≤ µ1 − µk+1 ≤ µ1 + α||$k − $∗ ||2Φ ≤
αµ1
µ1
) + 1−α
. Let k goes to infinity,
µ1 + αk+1 (||$1 − $∗ ||2Φ − 1−a
we have,
∞
X
||$k+1 − $k ||2Φ < ∞
(88)

Therefore, (87) yields µk+1 − µk ≤ −||$k+1 − $k ||2Φ .
Therefore, µk+1 ≤ µk ≤ µ1 . By the definition of µk , ||$k −$∗ ||2Φ −

(90)

Let J → ∞, then since 0 < α < 1,
P∞
P
φ1
i P∞
+ ∞
≤ 1−α
i=0 α Pk=i+1 ψk−i
k=1 φk
P
φ1
i ∞
(91)
+ ∞
≤ 1−α
t=1 ψt
i=0 α
φ1
1 P∞
≤ 1−α + 1−α t=1 ψt
P
P∞
2
2
Noticing that ∞
t=1 ψt =
k=1 (α + α )||$k − $k−1 ||Φ < ∞,
P∞
Pk
k=1 φk < ∞, and hence the sequence { i=1 φi }, being a
nonnegative and non-decreasing sequence, converges and is
bounded.
P
Consider another sequence {||$k − $∗ ||2Φ − ki=1 φi }. Since
P
||$k − $∗ ||2Φ is nonnegative and { ki=1 φi } is bounded, {||$k −
P
$∗ ||2Φ − ki=1 φi } is bounded from below. Furthermore, {||$k −
P
$∗ ||2Φ − ki=1 φi } is a non-increasing sequence. In fact,
P
||$k+1 − $∗ ||2Φ − k+1
i=1 φiP
= ||$k+1 − $∗ ||2Φ − φk+1 − ki=1 φi
P
≤ ||$k+1 − $∗ ||2Φ − ||$k+1 − $∗ ||2Φ + ||$k − $∗ ||2Φ − ki=1 φi
P
= ||$k − $∗ ||2Φ − ki=1 φi
where the second inequality follows from the definition of φk ,
P
φk+1 ≥ ||$k+1 − $∗ ||2Φ − ||$k − $∗ ||2Φ . As {||$k − $∗ ||2Φ − ki=1 φi }
is a non-increasing sequence and bounded from below, {||$k −
P
$∗ ||2Φ − ki=1 φi } converges.
Therefore, {||$k − $∗ ||2Φ }, being the sum of two convergent
P
P
sequences {||$k − $∗ ||2Φ − ki=1 φi } and { ki=1 φi }, also converges.
We are ready to show the convergence of Algorithm 6.1 using
the results in Part 1 and Part 2. Since {$k } is a bounded

16

sequence, it has a convergent subsequence {$nk } that converges
to $̃∗ . Because limk→∞ $k+1 − $k = 0 by Part 2, we have
limk→∞ $nk −1 − $nk = 0 and limk→∞ $nk +1 − $nk = 0. Pass
to limiting point of {$nk }, then we have $̃∗ = T 1 T 2 $̃∗ because
the righthand side of (63)-(64) is continuous. Hence, $̃∗ ∈
zer(Ā + B̄). Taking $∗ = $̃∗ in (71) of Part 1, we also have
{||$k − $̃∗ ||2Φ } converges by Part 3. Because a subsequence
{||$nk − $̃∗ ||2Φ } converges to zero, the whole sequence {||$k −
$̃∗ ||2Φ } converges to zero. Therefore, the whole sequence of
{$k } converges to $̃∗ ∈ zer(Ā + B̄). Resorting to Theorem
4.5 gives the desired result.
2

total product supply to all the markets given the action profile x
of all the companies. Market M j has a maximal capacity of
r j > 0, therefore, it should be satisfied that Ax ≤ r where
r = col(r1 , · · · , rm ) ∈ Rm . Suppose that P : Rm → Rm
is a price vector function that maps the total supply of each
market to the corresponding market’s price. Each company
has also a local production cost function ci (xi ) : Ωi → R.
Then the local objective function of company (player) Fi is
fi (xi , x−i ) : ci (xi ) − PT (Ax)Ai xi .
Overall, in this network Cournot game, each company needs
to solve the following optimization problem given the other
companies’ profile x−i ,

Remark 6.4. A sufficient and simple choice of parameters to
ensure the conditions in Theorem (6.3) is δ = β1 ,  = α and
√
1
0 < α < 10 − 3. In fact, δ > 2β
implies that 2βδ could be
be any real number % > 1. If we take  = ςα, ς > 0, then the
quadratic inequality becomes α2 − (2 − 3% − ς%)α + 1 − % < 0.
Since 1 − % < 0, α2 − (2 − 3% − ς%)α + 1 − % takes value strictly
less than zero when α takes 0. By the continuity of quadratic
equation, there always exists 0 < α < 1 that ensures the above
quadratic inequality given any % > 1 and ς > 0.

min xi ∈Ωi
s.t.

ci (xi ) − PT (Ax)Ai xi
P
Ai xi ≤ r − Nj=1, j,i A j x j

(92)

Obviously, the above network Cournot game in (92) is a
particular problem of game in (7). Some practical decision
problems in engineering networks can be well described by the
network Cournot game in (92), such as the rate control game in
communication network (Yin, Shanbhag & Mehta (2011)) and
the demand response game in smart grids (Ye & Hu (2016)).
Example 7.1 (Rate control game). Consider a group of
source-destination pairs (nodes) in a communication
network, that is {S 1 , ..., S N }, to decide their data rates in
a non-cooperative setting. The data is transferred through a
group of communication links (channels), that is {L1 , ..., Lm },
and each link L j has a maximal data rate capacity of c j > 0.
Assume that an additional layer has decided the routine table
for each source-destination pair S i , which is encoded by
Ai ∈ Rm×ni . Each column of Ai has only one element being 1
and all the other elements being zero, and the k-th element of
column j is 1 if S i utilizes the link Lk and transfers data rate
[xi ] j on link Lk . The local decision variable of S i is the data
rate on each link that it utilizes, denoted by xi ∈ Rni . xi also
PN
has a local feasibility constraint xi ∈ Ωi . Denote n = i=1
ni ,
x = col(x1 , · · · , xN ) ∈ Rn , A = [A1 , · · · , AN ] ∈ Rm×n and
c = col(c1 , · · · , cm ) ∈ Rm . The total data rate on each
link should be less than the capacity of that link: Ax ≤ c.
Given the data rate profile of all the nodes x, the payoff
function of S i , Ji (xi , x−i ) : Rn → R, takes the form as
Ji (xi , x−i ) = −ui (xi ) + DT (Ax)Ai xi , where ui (xi ) : Ωi → R is the
utility of source S i , and D : Rm → Rm is a delay function that
maps the total data rate on each link to the unit delay of that
link. Thereby, the data rate control game can be well described
by the network Cournot game in (92).

7. Network Cournot game and simulation studies
There are various practical problems that can be well
modeled by the game in (7), such as the river basin pollution
game in Krawczyk & Uryasev (2000), the power market
competition in Contreras, Klusch, Krawczyk (2004), plug-in
electric vehicles charging management in Paccagnan, Gentile,
Parise, Kamgarpour & Lygeros (2016), and communication
network congestion game in Yin, Shanbhag & Mehta (2011).
All above examples can be regarded as the type of the network
Cournot game described below, which is a generalization of the
network Cournot competition in Bimpikis, Ehsani, & Ilkilic
(2014) by introducing additional market capacity constraints or
equivalently globally shared coupling affine constraints. This
type of network Cournot game with affine coupling constraints
also appeared in the numerical studies of Yu, van der Schaar &
Sayed (2016).
7.1. Network Cournot game
Suppose that there are N companies (players) with labels
F1 , · · · , F N and m markets with labels M1 , · · · , Mm . Company
Fi decides its strategy to participate in the competition in ni
markets by producing and delivering xi ∈ Rni amounts of
products to the markets it connects with. The production
limitation of company Fi is xi ∈ Ωi ⊂ Rni . Company Fi has
a local matrix Ai ∈ Rm×ni that specifies which market it will
participate in. The j-th column of Ai , that is [Ai ]: j , has only
one element being 1 and all other elements being 0, and [Ai ]: j
has its k-th element being 1 if and only if player Fi delivers
[xi ] j amount of production to the market Mk . Therefore,
matrices A1 , · · · , AN can be used to specify a bipartite graph
that represents the connections between the companies and the
PN
markets. Denote n = i=1
ni , x = col(x1 , · · · , xN ) ∈ Rn , and
PN
m×n
A = [A1 , · · · , AN ] ∈ R . Then Ax ∈ Rm = i=1
Ai xi is just the

Example 7.2 (Demand response game). Given a distribution
network in power grids, suppose that there are T time periods,
and each period has a desirable minimal total load shedding
di > 0. Suppose that there are N load managers (energy
management units or players) in the network, and each load
manager i can decide a local vector xi ∈ Rti as its local
load shedding vector in some specific time periods. Each load
manager i also has a local matrix Ai ∈ RT ×ti that specifies
which time period player i will participate. For j−th column
of Ai , it has one element being 1 while all other elements
17

F5

being zero. The k−th element of the j−th column of Ai is 1
if load manager i decides to decrease its load by [xi ] j at time k.
PN
Denote x = col(x1 , · · · , xN ), and A = [A1 , · · · , AN ] ∈ RT × i=1 ti ,
d = col(d1 , · · · , dT ). Naturally, it is required that the total load
shedding of all the load managers should meet the minimal
value, Ax ≥ d. Each player has a local feasible constraint
xi ∈ Ωi ⊂ Rti , and a cost function ci (xi ) : Ωi → R due
to local load shedding. P : RT → RT is the payment price
vector function that maps total load shedding of each period to
the payment price vector, therefore, PT (Ax)Ai xi is the payment
awards of player i for its load shedding. The disutility function
of player i is Ji (xi , x−i ) = ci (xi ) − PT (Ax)Ai xi given all the
players’ action profile x. All in all, the demand response
management game is well described by the network Cournot
game model in (92).

ci (xi ) − (P̄ − DAx)T Ai xi
P
ci (xi ) − P̄T Ai xi + ( Nj=1 DA j x j )T Ai xi
ci (xi ) + xiT ATi DAi xi − ATi P̄T xi
P
+ Nj=1, j,i xTj ATj DAi xi

AT1 DA2
2AT2 DA2
···
ATN DA2

···
···
···
···

AT1 DAN
AT2 DAN
···
2ATN DAN

F20

M5 F12

F8

M3

F13 F16

F9

F4

Figure 1: Network cournot game: An edge from Fi to M j on this graph implies
that Fi participates in the competition in M j .

F5

F14

F10

F1

F17

F18

F2
F15

F7
F6
F3
F4

F19

F12
F8
F9

F11

F20
F16

F13

Figure 2: Interference graph G f : An undirected edge from Fi to F j on this
graph means that the objective function of Fi depends on x j , and vice versa.
Company Fi is assumed to be able to observe the decisions of its neighbors on
this graph .

(93)

F(x) is strongly monotone24 and Lipschitz continuous, and
Assumptions 1 and 2 are satisfied.
Remark 7.3. Notice that in the network Cournot game of (92),
each player i’s local objective function ci (xi ) − PT (Ax)Ai xi
only depends on the decisions of the players that participate
the same markets as player i. Denote Z[a, b] as the set of
integers from a to b. Mathematically, j ∈ Nif if and only if
∃k ∈ Z[1, m], p ∈ Z[1, ni ] and q ∈ Z[1, n j ] such that [Ai ]kp = 1
and [A j ]kq = 1. Since A = [A1 , · · · , AN ] is usually a sparse
matrix, the interference graph G f of network Cournot game
also has sparse edge connections.




 ,


(94)

7.2. Simulation studies
In the studies, we adopt a similar simulation setting as
Yu, van der Schaar & Sayed (2016) without considering the
stochastic factors. Consider 20 companies and 7 markets, and
the connection relationship between the companies and the
markets is depicted in Figure 1. If there is an edge from Fi to
M j in Figure 1, then company Fi participates the competition in

Notice that Q can be written as
Q = diag{AT1 DA1 , · · · , ATN DAN } + S T S

F19
M7

M4

F3

F18

M6
F15

F7
F6

then
F(x) = ∇c(x) + Qx − AT P̄.

F11

M2

Denote F(x) = col(∇ x1 J1 (x1 , x−1 ), · · · , ∇ xN JN (xN , x−N )),
∇c(x) = col(∇c1 (x1 ), ∇c2 (x2 ), · · · , ∇cN (xN )) and Q ∈ Rn×n , and

 2AT1 DA1
 AT DA
1
Q =  2
···

ATN DA1

F14 F
17

M1

F2

Moreover, the Assumptions 1 and 2 for Algorithm 3.5 and
6.1 can easily be satisfied for many practical cost functions and
price functions. For example, take company Fi ’s production
cost function to be a strongly convex function with Lipschitz
continuous gradients (A quadratic function ci (xi ) = xiT Qi xi +
bTi xi with Qi ∈ Rni ×ni being a symmetric and positive definite
matrix and bi ∈ Rni is one possible choice). The price of
market M j is taken as the linear function of the total supplying
p j (x) = P̄ j − d j [Ax] j (known as a linear inverse demand
function in economics) with P̄ j > 0, d j > 0. Denote P =
col(p1 , · · · , pm ) : Rn → Rm , P̄ = col(P̄1 , · · · , P̄m ) ∈ Rm ,
D = diag{d1 , · · · , dm } ∈ Rm×m . Then P = P̄ − DAx is the
vector price function. The payments of company Fi by selling
product xi to the markets that it connects with is just PT Ai xi .
Therefore, the objective function of company Fi is,
Ji (xi , x−i ) =
=
=

F10

F1

(95)

where
S is
√ a block matrix
√ defined√ as S √ =
√
[ DA1 , · · · , DAN ] ∈ Rm×n and D = diag{ d1 , · · · , dm }.
Therefore,
Q
is
positive
semi-definite
matrix.
Hence, the Jacobian matric of F(x), JF(x)
=
diag{∇2 c1 (x1 ), ∇2 c2 (x2 ), · · · , ∇2 cN (xN )} + Q is postive definite
since the cost functions ci (xi ) are strongly convex. Therefore,

24 Proposition

18

2.3.2 of Facchinei and Pang (2007)

F1

F2

F3

F4

F5

F6

F7

F20

F8

F19

F9

F18

F10

F17

F16

F15

F14

F13

F12

F11

Figure 3: Multiplier graph Gλ : Company Fi and F j are able to exchange their
local {λi , zi } and {λ j , z j } if there exists an edge between them on this graph.
||x −x∗ ||

||$ −$∗ ||

2
2
k
k
× 100% and ||$
× 100%
Figure 5:
The trajectories of ||x
∗ ||
∗ ||
2
2
generated with Algorithm 3.5 and 6.1: This also shows the convergence of
both algorithms, and the superior convergence speed of Algorithm 6.1.

Figure 4: The trajectories of ||xk+1 − xk ||2 and ||$k+1 − $k ||2 generated with
Algorithm 3.5 and 6.1: This shows the convergence of both algorithms, and the
superior convergence speed of the algorithm with inertia 6.1.

market M j by producing and delivering products to market M j .
Each company Fi has a local constraint as 0 < xi < Θi where
each component of Θi is randomly drawn from (10, 25). Each
market M j has a maximal capacity of r j > 0, j = 1, · · · , 7 and r j
is randomly drawn from (20, 80). The local objective function
is taken as (93). The local cost function of company i is ci (xi ) =
Pi
πi ( nj=1
[xi ] j )2 + bTi xi which is a strongly convex function with
Lipschitz continuous gradient. Here πi is randomly drawn from
(1, 8), and each component of bi is randomly drawn from (1, 4).
The price function is taken as the linear function P = P̄ − DAx,
and P̄ j and d j are randomly drawn from (250, 500) and (1, 5),
respectively.
With Figure 1 and the definition of objective function in (93),
the interference graph G f can be easily obtained and is depicted
in Figure 2. Meanwhile, we adopte the multiplier graph Gλ
shown in Figure 3. The weighted adjacency matrix W = [wi j ]
of multiplier graph Gλ has all its nonzero elements to be 1.
With Figure 1 and the definition of objective function in (93),
the interference graph G f can be easily obtained and is depicted
in Figure 2. Meanwhile, we adopte the multiplier graph Gλ
shown in Figure 3. The weighted adjacency matrix W = [wi j ]
of multiplier graph Gλ has all its nonzero elements to be 1.
Set the step-sizes in Algorithm 3.5 as τi = 0.03, νi = 0.2, σ =
0.02 for all companies, and for Algorithm 6.1 set α = 0.12
while other step-sizes are the same as Algorithm 3.5. The initial
starting points xi,0 , λi,0 and zi,0 of both algorithms are set to be
zeros.

Figure 6: The trajectory of ||(L ⊗ I7 )λ̄k ||2 and [(1T20 ⊗ I7 )λ̄k ]T (Axk − r): The
left sub-figure shows that λ̄k converges to the null space of L ⊗ I7 , hence the
local multipliers of all companies reach consensus asymptotically. (1T20 ⊗ I7 )λ̄k
P
is the averaged vector 20
i=1 λi,k , therefore, the right sub-figure shows that the
complementary condition hλ∗ , Ax∗ − ri = 0 is asymptotically satisfied.

Figure 7: The trajectories of local decisions xi,k of F1 , F6 , F10 and F16 . The
solid red lines are generated with Algorithm 3.5 and dashed lines are generated
with Algorithm 6.1

19

extended for stochastic GNE seeking with noisy gradient
observations and noisy information sharing by resorting to the
stochastic forward-backward splitting algorithm in Rosasco,
Villa & Vũ (2016). The strong monotonicity assumption on the
pseudo-gradient might be relaxed to monotonicity assumption
by utilizing the forward-backward-forward splitting method in
Briceno-Arias & Combettes (2013).
Reference
References
Figure 8: The trajectories of some components of local multiplier [λi,k ] j where
solid red lines are generated with Algorithm 3.5 and dashed lines are generated
with Algorithm 6.1. [λi,k ] j is the j−th component of local multiplier λi at
iteration time k. Here we select local multipliers of F1 , F6 , F10 and F16

Alpcan, T. and Basar, T., 2005. Distributed algorithms for Nash equilibria
of flow control games. In Advances in dynamic games (pp. 473-498).
Birkhauser Boston.
Alvarez, F. and Attouch, H., 2001. An inertial proximal method for maximal
monotone operators via discretization of a nonlinear oscillator with
damping. Set-Valued Analysis, 9(1-2), pp.3-11.
Attouch, H., Chbani, Z., Peypouquet, J. and Redont, P., 2016. Fast convergence
of inertial dynamics and algorithms with asymptotic vanishing viscosity.
Mathematical Programming, pp.1-53.
Bauschke, H.H. and Combettes, P.L., 2011. Convex analysis and monotone
operator theory in Hilbert spaces. Springer Science & Business Media.
Berinde, V., 2007. Iterative approximation of fixed points. Berlin, Germany:
Springer.
Bimpikis, K., Ehsani, S. and Ilkilic, R., 2014. Cournot competition in
networked markets. In EC (p. 733).
Briceno-Arias, L.M. and Combettes, P.L., 2013. Monotone operator methods
for Nash equilibria in non-potential games. In Computational and Analytical
Mathematics (pp. 143-159). Springer New York.
Combettes, P.L. and Vũ, B.C., 2014. Variable metric forward-backward
splitting with applications to monotone inclusions in duality. Optimization,
63(9), pp.1289-1318.
Contreras, J., Klusch, M. and Krawczyk, J.B., 2004. Numerical solutions to
Nash-Cournot equilibria in coupled constraint electricity markets. IEEE
Transactions on Power Systems, 19(1), pp.195-206.
Condat, L., 2013. A primal-dual splitting method for convex optimization
involving Lipschitzian, proximable and linear composite terms. Journal of
Optimization Theory and Applications, 158(2), pp.460-479.
Debreu, G., 1952. A social equilibrium existence theorem. Proceedings of the
National Academy of Sciences, 38(10), pp.886-893.
Dreves, A. and Sudermann-Merx, N., 2016. Solving linear generalized Nash
equilibrium problems numerically. Optimization Methods and Software,
pp.1-28.
Facchinei, F. and Pang, J.S., 2007. Finite-dimensional variational inequalities
and complementarity problems. Springer Science & Business Media.
Facchinei, F. and Pang, J.S., 2010. Nash equilibria: the variational approach.
Convex optimization in signal processing and communications, p.443.
Facchinei, F. and Kanzow, C., 2010. Generalized Nash equilibrium problems.
Annals of Operations Research, 175(1), pp.177-211.
Facchinei, F., Fischer, A. and Piccialli, V., 2007. On generalized Nash games
and variational inequalities. Operations Research Letters, 35(2), pp.159-164.
Fischer, A., Herrich, M. and Schonefeld, K., 2014. Generalized Nash
equilibrium problems-recent advances and challenges. Pesquisa
Operacional, 34(3), pp.521-558.
Harker, P.T., 1991. Generalized Nash games and quasi-variational inequalities.
European journal of Operational research, 54(1), pp.81-94.
Iutzeler, F. and Hendrickx, J.M., 2016. A Generic Linear Rate Acceleration
of Optimization algorithms via Relaxation and Inertia. arXiv preprint
arXiv:1603.05398.
Kearns, M., Littman, M.L. and Singh, S., 2001. Graphical models for game
theory. In Proceedings of the Seventeenth conference on Uncertainty in
artificial intelligence (pp. 253-260). Morgan Kaufmann Publishers Inc..
Koshal, J., Nedić, A. and Shanbhag, U.V., 2016. Distributed Algorithms for
Aggregative Games on Graphs. Operations Research, 64(3), pp.680-704.
Krawczyk, J.B. and Uryasev, S., 2000. Relaxation algorithms to find
Nash equilibria with economic applications. Environmental Modeling &
Assessment, 5(1), pp.63-73.

The trajectories of selected algorithm performance
indexes,
∗
∗
k −$ ||
k −x ||
× 100%, ||$||$
×
including ||xk+1 − xk ||2 , ||$k+1 − $k ||2 , ||x||x
∗ ||
∗ ||
100%, ||L ⊗ I7 λ̄k ||2 and [(1T20 ⊗ I7 )λ̄k ]T (Axk − r), are shown in
Figures 4, 5 and 6. The trajectories of the local decisions xi,k of
some companies are shown in Figure 7, and the trajectories of
[λi,k ] j , which stands for the j-th component of local Lagrangian
multiplier λi,k , are shown in Figure 8.
8. Conclusions
In this paper, we proposed a primal-dual distributed
algorithm based on operator splitting methods for iterative
computation of a variational GNE in noncooperative games
with globally shared affine coupling constraints. The algorithm
is motivated by the forward-backward operator splitting method
for finding zeros of a sum of monotone operators. Each player
only needs to knows its local information, especially a block
of the affine coupling constraints. The proposed algorithm
is proved to converge with fixed step-sizes under some
mild assumptions by exploiting properties of composition of
averaged operators. Furthermore, a distributed algorithm with
inertia is also proposed and analyzed for possible acceleration
of convergence speed. Numerical simulation studies for
a network Cournot game demonstrate the efficiency of the
proposed algorithms and the superior convergence speed of the
inertial algorithm.
Many challenging and exciting topics are still open for
distributed NE/GNE seeking. Here we only list some problems
with probable solution hints. Finding all the generalized
Nash equilibria has its only interests, and this could be
partially solved by combining the design in this paper with
the parameterized variational inequality method in Nabetani,
Tseng, and Fukushima (2011). The algorithm requires that
each player is able to observe all its neighbors’ decisions
through the interference graph G f . This assumption could
be relaxed by adopting the local consensus dynamics in
Salehisadaghiani & Pavel (2016b), and then it could only
be required that the players were able to observe parts of its
neighbors’ decisions through a maximal triangle-free spanning
subgraph of G f . The methodology of this paper could be
20

Kulkarni, A.A. and Shanbhag, U.V., 2012. On the variational equilibrium
as a refinement of the generalized Nash equilibrium. Automatica, 48(1),
pp.45-55.
Lei, J., Chen, H.F. and Fang, H.T., 2016. Primal-Dual Algorithm for Distributed
Constrained Optimization. Systems & Control Letters, 96, 110-117.
Liang, S., Yi, P. and Hong, Y., 2016. Distributed Nash equilibrium
seeking for aggregative games with coupled constraints. arXiv preprint
arXiv:1609.02253.
Lorenz, D.A. and Pock, T., 2015. An inertial forward-backward algorithm for
monotone inclusions. Journal of Mathematical Imaging and Vision, 51(2),
pp.311-325.
Menache, I. and Ozdaglar, A., 2011. Network games: Theory, models, and
dynamics. Synthesis Lectures on Communication Networks, 4(1), pp.1-159.
Mesbahi, M. and Egerstedt, M. (2010). Graph Theoretic Methods for
Multiagent Networks, Princeton: Princeton University Press.
Meyer, C.D., 2000. Matrix analysis and applied linear algebra (Vol. 2). SIAM.
Nabetani, K., Tseng, P. and Fukushima, M., 2011. Parametrized variational
inequality approaches to generalized Nash equilibrium problems with
shared constraints. Computational Optimization and Applications, 48(3),
pp.423-452.
Nesterov, Y., 2013. Introductory lectures on convex optimization: A basic
course (Vol. 87). Springer Science & Business Media.
Paccagnan D., Gentile B., Parise F.,Kamgarpour M. and Lygeros J. 2016.
Distributed computation of generalized Nash equilibria in quadratic
aggregative games with affine coupling constraints. In 2016 55th IEEE
Conference on Decision and Control (CDC), (pp.6123-6128), IEEE.
Pan, Y. and Pavel, L., 2009. Games with coupled propagated constraints in
optical networks with multi-link topologies. Automatica, 45(4), pp.871-880.
Parise, F., Gentile, B., Grammatico, S. and Lygeros, J., 2015, December.
Network aggregative games: Distributed convergence to Nash equilibria.
In 2015 54th IEEE Conference on Decision and Control (CDC) (pp.
2295-2300). IEEE.
Pavel, L., 2007. An extension of duality to a game-theoretic framework.
Automatica, 43(2), pp.226-237.
Poljak B. T., 1987. Introduction to optimization. Optimization Software.
Rosasco, L., Villa, S. and Vũ, B.C., 2016. A stochastic inertial
forward-backward splitting algorithm for multivariate monotone inclusions.
Optimization, 65(6), pp.1293-1314.
Rosen, J.B., 1965. Existence and uniqueness of equilibrium points for concave
n-person games. Econometrica: Journal of the Econometric Society,
pp.520-534.
Salehisadaghiani, F. and Pavel, L., 2016a. Distributed Nash equilibrium
seeking: A gossip-based algorithm. Automatica, 72, pp.209-216.
Salehisadaghiani, F. and Pavel, L., 2016b. Distributed Nash Equilibrium
Seeking By Gossip in Games on Graphs. arXiv preprint arXiv:1610.01896.
Schiro, D.A., Pang, J.S. and Shanbhag, U.V., 2013. On the solution of affine
generalized Nash equilibrium problems with shared constraints by Lemkes
method. Mathematical Programming, 142(1-2), pp.1-46.
Scutari, G., Palomar, D.P., Facchinei, F. and Pang, J.S., 2010. Convex
optimization, game theory, and variational inequality theory. IEEE Signal
Processing Magazine, 27(3), pp.35-49.
Shi, W., Ling, Q., Wu, G. and Yin, W., 2015. A proximal gradient algorithm
for decentralized composite optimization. IEEE Transactions on Signal
Processing, 63(22), pp.6013-6023.
Spall, J.C., 1992. Multivariate stochastic approximation using a simultaneous
perturbation gradient approximation. IEEE transactions on Automatic
Control, 37(3), pp.332-341.
Swenson, B., Kar, S. and Xavier, J., 2015. Empirical centroid fictitious play: An
approach for distributed learning in multi-agent games. IEEE Transactions
on Signal Processing, 63(15), pp.3888-3901.
Wibisono, A., Wilson, A.C. and Jordan, M.I., 2016. A Variational Perspective
on Accelerated Methods in Optimization. Proceedings of the National
Academy of Sciences, 133, doi:10.1073/pnas.1614734113
Ye, M. and Hu, G., 2016. Game design and analysis for price-based demand
response: An aggregate game approach. IEEE Transactions on Cybernetics,
to appear.
Yi P, Hong Y, Liu F. 2016. Initialization-free Distributed Algorithms
for Optimal Resource Allocation with Feasibility Constraints and its
Application to Economic Dispatch of Power Systems. Automatica 74,
259-269.
Yin, H., Shanbhag, U.V. and Mehta, P.G., 2011. Nash equilibrium problems

with scaled congestion costs and shared constraints. IEEE Transactions on
Automatic Control, 56(7), pp.1702-1708.
Yu, C. K., van der Schaar, M. and Sayed, A. H., 2016. Distributed Learning
for Stochastic Generalized Nash Equilibrium Problems. arXiv preprint
arXiv:1608.00039.
Zeng X, Yi P, Hong Y, Xie L. 2016. Continuous-Time Distributed
Algorithms for Extended Monotropic Optimization Problems, arXiv
preprint arXiv:1608.01167, 2016
Zhu, M. and Frazzoli, E., 2016. Distributed robust adaptive equilibrium
computation for generalized convex games. Automatica, 63, pp.82-91.

21

