A distributed primal-dual algorithm for computation of generalized Nash equilibria with shared affine coupling constraints via operator splitting methods

TL;DR

This paper introduces a distributed primal-dual algorithm based on operator splitting methods for computing generalized Nash equilibria with shared affine constraints in networked noncooperative games, ensuring convergence under mild conditions.

Contribution

It develops a novel distributed primal-dual algorithm utilizing operator splitting for variational GNE computation with shared constraints, allowing decentralized implementation and convergence guarantees.

Findings

Algorithm converges to variational GNE under fixed step-sizes.

Distributed approach requires only local information and neighbor communication.

Numerical simulations demonstrate the algorithm's efficiency in network Cournot competition.

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.

Equations

$$\partial f:x\mapsto \{g\in {\mathbf{R}}^m|f(y)\ge f(x)+⟨g,y-x⟩,\forall y\in domf\}.$$

$$Pro{x}_f:x\mapsto \arg \underset{u\in domf}{\min }f(u)+\frac{1}{2}||u-x|{|}_2^2.$$

$$N_{\Omega }(x)=\{\begin{array}{ll}\varnothing & x\notin \Omega \\ \{v|⟨v,y-x⟩\le 0,\forall y\in \Omega \}& x\in bd(\Omega )\\ 0& x\in int(\Omega )\end{array}$$

$$Pro{x}_{{\iota }_{\Omega }}(x)=R_{N_{\Omega }}(x)=P_{\Omega }(x).$$

$$\underset{\genfrac{}{}{0ex}{}{x\ne 0,}{1^{T}x=0}}{\min }{x}^{T}Lx=s_2||x|{|}_2^2,\underset{x\ne 0}{\max }{x}^{T}Lx=s_N||x|{|}_2^2.$$

$$d^{\ast }\le s_N\le 2d^{\ast }.$$

$$\underset{x_i\in {\mathbf{R}}^{n_i}}{\min }{f}_i(x_i,{\mathbf{x}}_{-i})s.t.,x_i\in X_i({\mathbf{x}}_{-i}).$$

$$x_i^{\ast }\in \arg \min f_i(x_i,{\mathbf{x}}_{-i}^{\ast }),s.t.x_i\in X_i({\mathbf{x}}_{-i}^{\ast }),\forall i\in \mathcal{N}.$$

$$X\subset {\mathbf{R}}^m:=\prod _{i=1}^N{\Omega }_i\bigcap \{\mathbf{x}\in {\mathbf{R}}^n|A\mathbf{x}\ge b\},$$

$$X_i({\mathbf{x}}_{-i}):=\{x_i\in {\Omega }_i|A_i{x}_i\ge b-\sum _{j\ne i,j\in \mathcal{N}}{A}_j{x}_j\}.$$

$$\underset{x_i}{\min }{f}_i(x_i,{\mathbf{x}}_{-i}^{\ast }),s.t.x_i\in {\Omega }_i,A_i{x}_i\ge b-\sum _{j\ne i,j\in \mathcal{N}}{A}_j{x}_j^{\ast }.$$

$$L_i(x_i,{\lambda }_i;{\mathbf{x}}_{-i})=f_i(x_i,{\mathbf{x}}_{-i})+{\lambda }_i^T(b-A\mathbf{x}).$$

$$\begin{array}[]{lll}\mathbf{0}\in\nabla_{x_{i}}L_{i}(x^{*}_{i},\lambda^{*}_{i};\mathbf{x}^{*}_{-i})+N_{\Omega_{i}}(x_{i}^{*}),x^{*}_{i}\in\Omega_{i}\\ \langle\lambda^{*}_{i},b-A\mathbf{x}^{*}\rangle=0,b-A\mathbf{x}^{*}\leq\mathbf{0},\lambda_{i}^{*}\geq\mathbf{0},\end{array}$$

$$\begin{array}[]{lll}\mathbf{0}\in\nabla_{x_{i}}f_{i}(x^{*}_{i},\mathbf{x}^{*}_{-i})-A_{i}^{T}\lambda_{i}^{*}+N_{\Omega_{i}}(x^{*}_{i})\\ \mathbf{0}\in(A\mathbf{x}^{*}-b)+N_{\mathbf{R}^{m}_{+}}(\lambda_{i}^{*})\end{array}$$

$$F(\mathbf{x})=col({\nabla }_{x_1}{f}_1(x_1,{\mathbf{x}}_{-1}),\cdots ,{\nabla }_{x_N}{f}_N(x_N,{\mathbf{x}}_{-N})),$$

$$Find{\mathbf{x}}^{\ast }\in X,⟨F({\mathbf{x}}^{\ast }),\mathbf{x}-{\mathbf{x}}^{\ast }⟩\ge 0,\forall \mathbf{x}\in X.$$

$$\underset{\mathbf{x}\in {\mathbf{R}}^n}{\min }⟨F({\mathbf{x}}^{\ast }),\mathbf{x}⟩,s.t.,\mathbf{x}\in \Omega ,A\mathbf{x}\ge b$$

$$\begin{array}[]{lll}\mathbf{0}\in\nabla_{x_{i}}f_{i}(x^{*}_{i},\mathbf{x}^{*}_{-i})-A_{i}^{T}\lambda^{*}+N_{\Omega_{i}}(x^{*}_{i}),i=1,\cdots,N\\ \mathbf{0}\in(A\mathbf{x}^{*}-b)+N_{\mathbf{R}^{m}_{+}}(\lambda^{*})\end{array}$$

$$\begin{array}[]{l}\mathfrak{A}:\;\left(\begin{array}[]{c}\mathbf{x}\\ \lambda\\ \end{array}\right)\mapsto\left(\begin{array}[]{c}F(\mathbf{x})\\ -b\\ \end{array}\right)\\ \mathfrak{B}:\;\left(\begin{array}[]{c}\mathbf{x}\\ \lambda\\ \end{array}\right)\mapsto\left(\begin{array}[]{c}N_{\Omega}(\mathbf{x})\\ N_{\mathbf{R}^{m}_{+}}(\lambda))\\ \end{array}\right)+\left(\begin{array}[]{cc}\mathbf{0}&-A^{T}\\ A&\mathbf{0}\\ \end{array}\right)\left(\begin{array}[]{c}\mathbf{x}\\ \lambda\\ \end{array}\right)\par\end{array}$$

$$\displaystyle x_{i,k+1}=P_{\Omega_{i}}\big{[}x_{i,k}-\tau_{i}(\nabla_{x_{i}}f_{i}(x_{i,k},\{x_{j,k}\}_{j\in\mathcal{N}^{f}_{i}})-A_{i}^{T}\lambda_{i,k})\big{]}$$

$$z_{i,k+1}=z_{i,k}+{\nu }_i\sum _{j\in {\mathcal{N}}_i^{\lambda }}{w}_{ij}({\lambda }_{i,k}-{\lambda }_{j,k})$$

$$\displaystyle\lambda_{i,k+1}=P_{\mathbf{R}^{m}_{+}}\Big{\{}\lambda_{i,k}-\sigma_{i}\big{[}2A_{i}x_{i,k+1}-A_{i}x_{i,k}-b_{i}+\sum_{j\in\mathcal{N}^{\lambda}_{i}}w_{ij}(\lambda_{i,k}-\lambda_{j,k})$$

$$\displaystyle\quad\qquad+2\sum_{j\in\mathcal{N}^{\lambda}_{i}}w_{ij}(z_{i,k+1}-z_{j,k+1})-\sum_{j\in\mathcal{N}^{\lambda}_{i}}w_{ij}(z_{i,k}-z_{j,k})\big{]}\Big{\}}$$

$$\displaystyle\mathbf{x}_{k+1}=P_{\Omega}\big{[}\mathbf{x}_{k}-\bar{\tau}(F(\mathbf{x}_{k})-\Lambda^{T}\bar{\lambda}_{k})\big{]}$$

$${\bar{z}}_{k+1}={\bar{z}}_k+\bar{\nu }\bar{L}{\bar{\lambda }}_k$$

$$\displaystyle\bar{\lambda}_{k+1}=P_{\mathbf{R}^{mN}_{+}}\Big{\{}\bar{\lambda}_{k}-\bar{\sigma}\big{[}\Lambda(2\mathbf{x}_{k+1}-\mathbf{x}_{k})-\bar{b}+\bar{L}\bar{\lambda}_{k}+\bar{L}(2\bar{z}_{k+1}-\bar{z}_{k})\big{]}\Big{\}}$$

$${\mathbf{x}}_k-\bar{\tau }(F({\mathbf{x}}_k)-{\Lambda }^T{\bar{\lambda }}_k)\in {\mathbf{x}}_{k+1}+N_{\Omega }({\mathbf{x}}_{k+1}).$$

$$-F({\mathbf{x}}_k)\in N_{\Omega }({\mathbf{x}}_{k+1})-{\Lambda }^T{\bar{\lambda }}_{k+1}+{\bar{\tau }}^{-1}({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)+{\Lambda }^T({\bar{\lambda }}_{k+1}-{\bar{\lambda }}_k).$$

$$\bar{\lambda}_{k}-\bar{\sigma}\Big{[}\Lambda(2\mathbf{x}_{k+1}-\mathbf{x}_{k})-\bar{b}+\bar{L}\bar{\lambda}_{k}+\bar{L}(2\bar{z}_{k+1}-\bar{z}_{k})\Big{]}\in\bar{\lambda}_{k+1}+N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda}_{k+1})$$

$$\begin{array}[]{lll}-[\bar{L}\bar{\lambda}_{k}-\bar{b}]&\in&N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda}_{k+1})+\Lambda\mathbf{x}_{k+1}+\bar{L}\bar{z}_{k+1}\\ &+&\Lambda(\mathbf{x}_{k+1}-\mathbf{x}_{k})+\bar{L}(\bar{z}_{k+1}-\bar{z}_{k})+\bar{\sigma}^{-1}(\bar{\lambda}_{k+1}-\bar{\lambda}_{k})\end{array}$$

$$\displaystyle-\left(\begin{array}[]{c}F(\mathbf{x}_{k})\\ \mathbf{0}\\ \bar{L}\bar{\lambda_{k}}-\bar{b}\\ \end{array}\right)\in\left(\begin{array}[]{c}N_{\Omega}(\mathbf{x}_{k+1})-\Lambda^{T}\bar{\lambda}_{k+1}\\ -\bar{L}\bar{\lambda}_{k+1}\\ N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda}_{k+1})+\Lambda\mathbf{x}_{k+1}+\bar{L}\bar{z}_{k+1}\\ \end{array}\right)$$

$$\displaystyle\qquad+\left(\begin{array}[]{ccc}\bar{\tau}^{-1}&\mathbf{0}&\Lambda^{T}\\ \mathbf{0}&\bar{\nu}^{-1}&\bar{L}\\ \Lambda&\bar{L}&\bar{\sigma}^{-1}\\ \end{array}\right)\left(\begin{array}[]{c}\mathbf{x}_{k+1}-\mathbf{x}_{k}\\ \bar{z}_{k+1}-\bar{z}_{k}\\ \bar{\lambda}_{k+1}-\bar{\lambda}_{k}\\ \end{array}\right)$$

$$\Phi=\left(\begin{array}[]{ccc}\bar{\tau}^{-1}&\mathbf{0}&\Lambda^{T}\\ \mathbf{0}&\bar{\nu}^{-1}&\bar{L}\\ \Lambda&\bar{L}&\bar{\sigma}^{-1}\\ \end{array}\right)$$

Full text

††journal: Elseviert1t1footnotetext: This work was supported by NSERC Discovery Grant (261764).

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 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 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,Nedić & 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 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 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.

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, $\mathbf{R}^{m}$ ($\mathbf{R}^{m}_{+}$) denotes the $m-$dimesional (nonnegative) Euclidean space. For a column vector $x\in\mathbf{R}^{m}$ (matrix $A\in\mathbf{R}^{m\times n}$), $x^{T}$ ($A^{T}$) denotes its transpose. $x^{T}y=\langle x,y\rangle$ denotes the inner product of $x,y$, and $||x||_{2}=\sqrt{x^{T}x}$ denotes the norm induced by inner product $\langle\cdot,\cdot\rangle$. Given a symmetric positive definite matrix $G$, denote the $G$-induced inner product $\langle x,y\rangle_{G}=\langle Gx,y\rangle$. The $G$-matrix induced norm, $||\cdot||_{G}$, is defined as $||x||_{G}=\sqrt{\langle Gx,x\rangle}$. Denote by $||\cdot||$ any matrix induced norm in the Euclidean space. Denote $\mathbf{1}_{m}=(1,...,1)^{T}\in\mathbf{R}^{m}$ and $\mathbf{0}_{m}=(0,...,0)^{T}\in\mathbf{R}^{m}$. For column vectors $x,y$, $x\geq(>)y$ is understood componentwise. $diag\{A_{1},...,A_{N}\}$ represents the block diagonal matrix with matrices $A_{1},...,A_{N}$ on its main diagonal. $Null(A)$ and $Range(A)$ denote the null space and range space of matrix $A$, respectively. Denote $col(x_{1},....,x_{N})$ as the column vector stacked with column vectors $x_{1},...,x_{N}$. $I_{n}$ denotes the identity matrix in $\mathbf{R}^{n\times n}$. For a matrix $A=[a_{ij}]$, $a_{ij}$ or $[A]_{ij}$ stands for the matrix entry in the $i$th row and $j$th column of $A$. We also use $[x]_{k}$ to denote the $k-$th element in column vector $x$. Denote $\times_{i=1,...,n}\Omega_{i}$ or $\prod_{i=1}^{n}\Omega_{i}$ as the Cartesian product of the sets $\Omega_{i},i=1,...,n$. Denote $int(\Omega)$ as the interior of $\Omega$, and $bd(\Omega)$ as the boundary set of $\Omega$. Define the projection of $x$ onto a set $\Omega$ by $P_{\Omega}(x)=\arg\min_{y\in\Omega}||x-y||_{2}$. A set $\Omega$ is a convex set if $\lambda x+(1-\lambda)y\in\Omega,\;\forall\lambda\in[0,1],\forall x,y\in\Omega$. An extended value proper function $f:\mathbf{R}^{m}\rightarrow\mathbf{R}$ is a convex function if $f(\lambda x+(1-\lambda)y)\leq\lambda f(x)+(1-\lambda)f(y),\forall\lambda\in[0,1],\forall x,y\in domf$.

2.1 Monotone operators

The following concepts are reviewed from Bauschke & Combettes (2011). Let $\mathfrak{A}:\mathbf{R}^{m}\rightarrow 2^{\mathbf{R}^{m}}$ be a set-valued operator. Denote ${\rm Id}$ as the identity operator, i.e, ${\rm Id}(x)=x$. The domain of $\mathfrak{A}$ is $dom\mathfrak{A}=\{x\in\mathbf{R}^{m}|\mathfrak{A}x\neq\emptyset\}$ where $\emptyset$ stands for the empty set, and the range of $\mathfrak{A}$ is $ran\mathfrak{A}=\{y\in\mathbf{R}^{m}|\exists x\in\mathbf{R}^{m},y\in\mathfrak{A}x\}$. The graph of $\mathfrak{A}$ is $gra\mathfrak{A}=\{(x,u)\in\mathbf{R}^{m}\times\mathbf{R}^{m}|u\in\mathfrak{A}x\}$, then the inverse of $\mathfrak{A}$ is defined through its graph as $gra\mathfrak{A}^{-1}=\{(u,x)\in\mathbf{R}^{m}\times\mathbf{R}^{m}|(x,u)\in gra\mathfrak{A}\}$. The zero set of operator $\mathfrak{A}$ is $zer\mathfrak{A}=\{x\in\mathbf{R}^{m}|\mathbf{0}\in\mathfrak{A}x\}$. Define the resolvent of operator $\mathfrak{A}$ as $R_{\mathfrak{A}}=({\rm Id}+\mathfrak{A})^{-1}$. An operator $\mathfrak{A}$ is called monotone if $\forall(x,u),\forall(y,v)\in gra\mathfrak{A}$, we have $\langle x-y,u-v\rangle\geq 0.$ Moreover, it is maximally monotone if $gra\mathfrak{A}$ is not strictly contained in the graph of any other monotone operator. $R_{\mathfrak{A}}$ is single-valued and $domR_{\mathfrak{A}}=\mathbf{R}^{m}$ if $\mathfrak{A}$ is maximally monotone *** Proposition 23.7 in Bauschke & Combettes (2011). For a proper lower semi-continuous convex (l.s.c.) function $f$, its sub-differential is a set-valued operator $\partial f:domf\rightarrow 2^{\mathbf{R}^{m}}$ and

[TABLE]

$\partial f$ is a maximally monotone operator *** Theorem 20.40 in Bauschke & Combettes (2011) . Then $Prox_{f}=R_{\partial f}:\mathbf{R}^{m}\rightarrow domf$ is called the proximal operator of $f$ *** Proposition 16.34 in Bauschke & Combettes (2011), i.e.

[TABLE]

Define the indicator function of set $\Omega$ as $\iota_{\Omega}(x)=\begin{cases}0,&x\in\Omega;\\ \infty,&x\notin\Omega.\end{cases}$ For a closed convex set $\Omega$, $\iota_{\Omega}$ is a proper l.s.c. function. $\partial\iota_{\Omega}$ is just the normal cone operator of set $\Omega$, that is $\partial\iota_{\Omega}(x)=N_{\Omega}(x)$ ***Example 16.12 of Bauschke & Combettes (2011) and

[TABLE]

In this case, we also have ***Example 23.4 of Bauschke & Combettes (2011)

[TABLE]

For a single-valued operator $T:\mathbf{R}^{m}\rightarrow\mathbf{R}^{m}$, a point $x\in\mathbf{R}^{m}$ is a fixed point of $T$ if $Tx=x$, and the set of fixed points of $T$ is denoted as $FixT$. The composition of operators $\mathfrak{A}$ and $\mathfrak{B}$, denoted by $\mathfrak{A}\circ\mathfrak{B}$, is defined via its graph $gra\mathfrak{A}\circ\mathfrak{B}=\{(x,z)|\exists y\in ran\mathfrak{B},(x,y)\in gra\mathfrak{B},(y,z)\in gra\mathfrak{A}\}$. We also use $\mathfrak{A}\mathfrak{B}$ to denote the composition $\mathfrak{A}\circ\mathfrak{B}$ when they are single-valued. Similarly, their sum $\mathfrak{A}+\mathfrak{B}$ is defined as $gra(\mathfrak{A}+\mathfrak{B})=\{(x,y+z)|(x,y)\in gra\mathfrak{A},(x,z)\in gra\mathfrak{B}\}$. Suppose operators $\mathfrak{A}$ and $\mathfrak{B}$ are maximally monotone and $0\in int(dom\mathfrak{A}-dom\mathfrak{B})$, then $\mathfrak{A}+\mathfrak{B}$ is also maximally monotone***Corollary 24.4 in Bauschke & Combettes (2011). Further suppose that $\mathfrak{A}$ is single-valued, then $zer(\mathfrak{A}+\mathfrak{B})=FixR_{\mathfrak{B}}\circ({\rm Id}-\mathfrak{A})$***Proposition 25.1 in Bauschke & Combettes (2011), 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 $\mathcal{G}=(\mathcal{N},\mathcal{E})$. $\mathcal{N}=\{1,\cdots,N\}$ is the set of agents, and the edge set $\mathcal{E}\subset\mathcal{N}\times\mathcal{N}$ contains all the information interactions. If agent $i$ can get information from agent $j$, then $(j,i)\in\mathcal{E}$ and agent $j$ belongs to agent $i$’s neighbor set $\mathcal{N}_{i}=\{j|(j,i)\in\mathcal{E}\}$, and $i\notin\mathcal{N}_{i}$. $\mathcal{G}$ is said to be undirected when $(i,j)\in\mathcal{E}$ if and only if $(j,i)\in\mathcal{E}$. A path of graph $\mathcal{G}$ is a sequence of distinct agents in $\mathcal{N}$ such that any consecutive agents in the sequence correspond to an edge of graph $\mathcal{G}$. Agent $j$ is said to be connected to agent $i$ if there is a path from $j$ to $i$. $\mathcal{G}$ is said to be connected if any two agents are connected.

Define the weighted adjacency matrix $W=[w_{ij}]\in\mathbf{R}^{N\times N}$ of $\mathcal{G}$ with $w_{ij}>0$ if $j\in\mathcal{N}_{i}$ and $w_{ij}=0$ otherwise. Assume $W=W^{T}$ for undirected graphs. Define the weighted degree matrix $Deg=diag\{d_{1},\cdots,d_{N}\}=diag\{\sum_{j=1}^{N}w_{1j},...,\sum_{j=1}^{N}w_{Nj}\}.$ Then the weighted Laplacian of graph $\mathcal{G}$ is $L=Deg-W.$ When graph $\mathcal{G}$ is a connected and undirected graph, 0 is a simple eigenvalue of Laplacian $L$ with the eigenspace $\{\alpha\mathbf{1}_{N}|\alpha\in\mathbf{R}\}$, and $L\mathbf{1}_{N}=\mathbf{0}_{N}$, $\mathbf{1}^{T}_{N}L=\mathbf{0}^{T}_{N}$, while all other eigenvalues are positive. Denote the eigenvalues of $L$ in an ascending order as $0<s_{2}\leq\cdots\leq s_{N}$, then by Courant-Fischer Theorem,

[TABLE]

Denote $d^{*}=\max\{d_{1},\cdots,d_{N}\}$ as the maximal weighted degree of graph $\mathcal{G}$, then we have the following estimation,

[TABLE]

3 Problem formulation and distributed algorithm

3.1 Game formulation

Consider a group of agents (players) $\mathcal{N}=\{1,\cdots,N\}$ that seek the generalized Nash equilibrium (GNE) of a noncooperative game with coupling constraints defined as follows. Each player $i\in\mathcal{N}$ controls its local decision (strategy or action) $x_{i}\in\mathbf{R}^{n_{i}}$. Denote $\mathbf{x}=col(x_{1},\cdots,x_{N})\in\mathbf{R}^{n}$ as the decision profile, i.e., the stacked vector of all the agents’ decisions where $\sum_{i=1}^{N}n_{i}=n$. Denote $\mathbf{x}_{-i}=col(x_{1},\cdots,x_{i-1},x_{i+1},\cdots,x_{N})$ 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 $f_{i}(x_{i},\mathbf{x}_{-i}):\mathbf{R}^{n}\rightarrow\mathbf{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 of the other players with $X_{i}(\mathbf{x}_{-i}):\mathbf{R}^{n-n_{i}}\rightarrow 2^{\mathbf{R}^{n_{i}}}$ denoting a set-valued map that maps $\mathbf{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 $\mathbf{x}_{-i}$,

[TABLE]

The GNE $\mathbf{x}^{*}=col(x_{1}^{*},\cdots,x_{N}^{*})$ of the game in (7) is obtained at the intersection of all the players’ best-response sets, and is defined as:

[TABLE]

Here we consider the GNE seeking in noncooperative games where the couplings between players’ feasible sets are specified by globally shared affine constraints. Denote

[TABLE]

where $\Omega_{i}\subset\mathbf{R}^{n_{i}}$ is a private feasible decision set of player $i$, and $A=[A_{1},\cdots,A_{N}]\in\mathbf{R}^{m\times n}$ with $A_{i}\in\mathbf{R}^{m\times n_{i}}$, and $b\in\mathbf{R}^{m}$. Denote $\Omega=\prod_{i=1}^{N}\Omega_{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$: $X_{i}(\mathbf{x}_{-i}):=\{x_{i}\in\mathbf{R}^{n_{i}}:(x_{i},\mathbf{x}_{-i})\in X\}$, or in other words:

[TABLE]

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

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 $\cdots$ practically is likely to be satisfied only when the joint constraints $g_{\mu}=g,\mu=1,...,N$ are linear, i.e. of the form $Ax\leq 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).

Assumption 1

For player $i$, $f_{i}(x_{i},\mathbf{x}_{-i})$ is a differentiable convex function with respect to $x_{i}$ given any fixed $\mathbf{x}_{-i}$, and $\Omega_{i}$ is a closed convex set. $X$ has nonempty interior point (Slater’s condition).

Suppose $\mathbf{x}^{*}$ is a GNE of game (7), then for agent $i$, $x_{i}^{*}$ is the optimal solution to the following convex optimization problem:

[TABLE]

Define a local Lagrangian function for agent $i$ with multiplier $\lambda_{i}\in\mathbf{R}^{m}_{+}$ as

[TABLE]

When $x_{i}^{*}$ is an optimal solution to (9), there exists $\lambda_{i}^{*}\in\mathbf{R}_{+}^{m}$ such that the following optimality conditions (KKT) are satisfied:

[TABLE]

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

[TABLE]

In fact, since $N_{\Omega}(x)=\emptyset$ when $x\notin\Omega$, it must hold that $x^{*}_{i}\in\Omega_{i}$ and $\lambda^{*}_{i}\in\mathbf{R}^{m}_{+}$ when (12) is satisfied. Furthermore, $N_{\mathbf{R}^{m}_{+}}=\prod_{j=1}^{m}N_{\mathbf{R}_{+}}$. If $[\lambda^{*}_{i}]_{k}=0$, then $N_{\mathbf{R}_{+}}([\lambda_{i}^{*}]_{k})=-\mathbf{R}_{+},$ and hence $[b-A\mathbf{x}^{*}]_{k}\leq 0$; and if $[\lambda^{*}_{i}]_{k}>0$, we have $N_{\mathbf{R}_{+}}([\lambda_{i}^{*}]_{k})=0,$ hence $[b-A\mathbf{x}^{*}]_{k}=0$. Therefore, $b-A\mathbf{x}^{*}\leq\mathbf{0}$ and $\langle\lambda^{*}_{i},b-A\mathbf{x}^{*}\rangle=0$. Denote $\bar{\lambda}=col(\lambda_{1},\cdots,\lambda_{N})$. Therefore, by Theorem 4.6 in Facchinei & Kanzow (2010) when $(\mathbf{x}^{*},\bar{\lambda}^{*})$ satisfies KKT (12) for $i=1,...,N$, $\mathbf{x}^{*}$ is a GNE of the game in (7)

According to the above discussions, given $\mathbf{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., $\lambda_{1}^{*}\neq\lambda_{2}^{*}\neq,\cdots,\neq\lambda_{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

[TABLE]

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 $VI(F,X)$:

[TABLE]

Let us check the KKT condition for $VI(F,X)$ in (14). In fact, $\mathbf{x}^{*}$ is a solution to $VI(F,X)$ in (14) if and only if $\mathbf{x}^{*}$ is the optimal solution to the following optimization problem:

[TABLE]

According to the optimization formulation of $VI(F,X)$ in (15), if $\mathbf{x}^{*}$ solves $VI(F,X)$, there exists $\lambda^{*}\in\mathbf{R}^{m}$ such that the following optimality conditions (KKT) are satisfied:

[TABLE]

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

Theorem 3.2

***Theorem 2.1 of Facchinei, Fischer & Piccialli (2007)

Suppose Assumption 1 holds. Every solution $\mathbf{x}^{*}$ of $VI(F,X)$ in (14) is a GNE of game in (7). Furthermore, if $\mathbf{x}^{*}$ together with $\lambda^{*}$ satisfies the KKT conditions for the $VI(F,X)$, i.e., (16), then $\mathbf{x}^{*}$ together with $\lambda_{1}^{*}=,\cdots,=\lambda^{*}_{N}=\lambda^{*}$ satisfies the KKT conditions for the GNE, i.e., (12).

The solution $\mathbf{x}^{*}$ of $VI(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 $VI(F,X)$ in (14), thus provides a distributed coordination mechanism such that the coupling constraint is met and a variational GNE of the game is found.

Define two operators $\mathfrak{A}$ and $\mathfrak{B}$, both from $\Omega\times\mathbf{R}^{m}_{+}$ to $\mathbf{R}^{n}\times\mathbf{R}^{m}$ as follows,

[TABLE]

By the definition of $F(\mathbf{x})$ in (13), the KKT conditions (16) can be equivalently written as $\mathbf{0}\in(\mathfrak{A}+\mathfrak{B})col(\mathbf{x}^{*},\lambda^{*})$. Notice that $\mathfrak{B}$ is a maximally monotone operator (similar arguments for this can be found in Lemma 5.14). Hence, if $F(\mathbf{x})$ has some additional properties, then solving $VI(F,X)$ can be converted to the problem of finding zeros of a sum of monotone operators.

Assumption 2

$F(\mathbf{x})$* defined in (13) is strongly monotone with parameter $\eta$ over $\Omega$: $\langle F(\mathbf{x})-F(\mathbf{y}),\mathbf{x}-\mathbf{y}\rangle\geq\eta||\mathbf{x}-\mathbf{y}||_{2}^{2},\forall\mathbf{x},\mathbf{y}\in\Omega$, and $\theta-$ Lipschitz continuous over $\Omega$: $||F(\mathbf{x})-F(\mathbf{y})||_{2}\leq\theta||\mathbf{x}-\mathbf{y}||_{2},\forall\mathbf{x},\mathbf{y}\in\Omega$.*

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 $VI(F,X)$ in (14) **Theorem 2.3.3 of Facchinei and Pang (2007), 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.

3.2 Distributed algorithm

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 $f_{i}(x_{i},\mathbf{x}_{-i})$ and local feasible set $\Omega_{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 $A_{i}$ specifies how player $i$ participates in the resource allocation or market behavior, hence also contains private information, and $b$ can be decomposed as $b=\sum_{i=1}^{N}b_{i}$. Thus matrix $A_{i}$ can be regarded as a map from decision space $\mathbf{R}^{n_{i}}$ to resource space $\mathbf{R}^{m}$, while vector $b_{i}$ 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 $x_{i}\in\mathbf{R}_{+}$, then $A_{i}\in\mathbf{R}^{m}_{+}$ satisfying $\mathbf{1}_{m}^{T}A_{i}=1$ and $\mathbf{0}\leq A_{i}\leq\mathbf{1}$ just specifies how each player allocates its production to each market. In this case, if $\tilde{b}_{i}\in\mathbf{R}^{m}_{+}$ is a local observation of market capacity vector, then the true market capacities can be taken as $b=\sum_{i=1}^{N}b_{i}=\sum_{i=1}^{N}\frac{1}{N}\tilde{b}_{i}$.

Therefore, we assume that player $i$ only knows its local $f_{i}(x_{i},\mathbf{x}_{-i})$, $\Omega_{i}$, and matrix $A_{i}$ and $b_{i}$ with $A=[A_{1},\cdots,A_{N}]$ and $b=\sum_{i=1}^{N}b_{i}$. In other words, player $i$ has a local first-order oracle of $f_{i}(x_{i},\mathbf{x}_{-i})$ which returns $\nabla_{x_{i}}f_{i}(x_{i},\mathbf{x}_{-i})$ given $(x_{i},\mathbf{x}_{-i})$, meanwhile player $i$ can manipulate $\Omega_{i}$, $A_{i}$ and $b_{i}$ for its local computation.

The players need to find the solution to $VI(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, $\mathcal{G}_{f}$ and $\mathcal{G}_{\lambda}$, related to the local information observations or exchanging between players.

Graph $\mathcal{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 $\mathcal{G}_{f}=(\mathcal{N},\mathcal{E}_{f})$, $(j,i)\in\mathcal{E}_{f}$ if the objective function of agent $i$, $f_{i}(x_{i},\mathbf{x}_{-i})$ explicitly depends on the decision of player $j$. We define $\mathcal{N}^{f}_{i}=\{j|(j,i)\in\mathcal{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 $f_{i}(x_{i},\{x_{j}\}_{j\in\mathcal{N}^{f}_{i}})$, and the local oracle of player $i$ returns $\nabla_{x_{i}}f_{i}(x_{i},\{x_{j}\}_{j\in\mathcal{N}^{f}_{i}})$, i.e., $\nabla_{x_{i}}f_{i}(x_{i},\mathbf{x}_{-i})$, given $\{x_{j}\}_{j\in\mathcal{N}^{f}_{i}}$. The local oracle might compute $\nabla_{x_{i}}f_{i}(x_{i},\{x_{j}\}_{j\in\mathcal{N}^{f}_{i}})$ by approximating $\nabla_{x_{i}}f_{i}(x_{i},\mathbf{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 $\lambda_{1}^{*}=,\cdots,=\lambda^{*}_{N}=\lambda^{*}$ in Theorem *, we also assume that the agents can exchange certain local information through a multiplier graph $\mathcal{G}_{\lambda}=(\mathcal{N},\mathcal{E}_{\lambda})$. $(j,i)\in\mathcal{E}_{\lambda}$ if player $i$ can receive certain information from player $j$, while the information to be shared through $\mathcal{G}_{\lambda}$ will be specified later. Thereby, player $i$ has its multiplier neighbors $\mathcal{N}^{\lambda}_{i}=\{j|(j,i)\in\mathcal{E}_{\lambda}\}$. $W=[w_{ij}]$ is the weighted adjacency matrix associated with multiplier graph $\mathcal{G}_{\lambda}$, and $L$ is the corresponding weighted Laplacian matrix.

Assumption 3

$\mathcal{G}_{\lambda}=\{\mathcal{N},\mathcal{E}_{\lambda}\}$* is undirected and connected. $W=W^{T}$.*

Remark 3.4

We assume that each agent can observe the decisions which its local objective function directly depends on through interference graph $\mathcal{G}_{f}$. Therefore, player $i$ can get its local gradient $\nabla_{x_{i}}f_{i}(x_{i},\{x_{j}\}_{j\in\mathcal{N}^{f}_{i}})$. This type of local 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 $f_{i}(x_{i},\mathbf{x}_{-i})$ does not explicitly depend on the decision of player $k$, i.e. $k\notin\mathcal{N}_{i}^{f}$. In fact, player $i$’s feasible decision set implicitly depends on all other players’ decisions through the globally shared affine coupling constraint: $A\mathbf{x}\geq 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 $\mathcal{G}_{\lambda}$ must be connected. Therefore, $\mathcal{G}_{f}$ and $\mathcal{G}_{\lambda}$ could be two different information observation or information sharing graphs because they serve for different purposes.

We are ready to present the main distributed algorithm after the introduction of algorithm notations. Agent (player) $i$ controls its local decision variable $x_{i}\in\mathbf{R}^{n_{i}}$ and local Lagragian multiplier $\lambda_{i}\in\mathbf{R}^{m}$. Meanwhile, we also assume that player $i$ has a local auxiliary variable $z_{i}\in\mathbf{R}^{m}$ for the coordinations needed to satisfy the affine coupling constraint and to reach consensus of the local multipliers $\lambda_{i}$. As indicated before, player $i$ can compute $\nabla_{x_{i}}f_{i}(x_{i},\mathbf{x}_{-i})$ by observing the adversary players’ decisions that its local objective function $f_{i}(x_{i},\mathbf{x}_{-i})$ directly depends on, that is the decisions of its interference neighbors in $\mathcal{N}^{f}_{i}$. On the other hand, player $i$ can also share information related to the local multiplier $\lambda_{i}$ and local auxiliary variable $z_{i}$ with its multiplier neighbours in $\mathcal{N}^{\lambda}_{i}$ through multiplier network $\mathcal{G}_{\lambda}$.

Next we show the distributed algorithm for agent $i$.

Algorithm 3.5

**

[TABLE]

$\tau_{i},\nu_{i},\sigma_{i}>0$ are fixed constant step-sizes of player $i$, and $W=[w_{ij}]$ is the weighted adjacency matrix of $\mathcal{G}_{\lambda}$.*

Algorithm 3.5 runs sequentially as follows. At the iteration time $k$, player $i$ gets $\nabla_{x_{i}}f_{i}(x_{i,k},\mathbf{x}_{-i,k})$ by observing $x_{j,k},j\in\mathcal{N}^{f}_{i}$ through interference graph $\mathcal{G}_{f}$, and updates $x_{i,k}$ with (18); meanwhile, player $i$ gets $\lambda_{j,k},j\in\mathcal{N}^{\lambda}_{i}$ through multiplier graph $\mathcal{G}_{\lambda}$, and updates $z_{i,k}$ by (19). Then player $i$ gets $z_{j,k+1},j\in\mathcal{N}^{\lambda}_{i}$ through multiplier graph $\mathcal{G}_{\lambda}$ and updates $\lambda_{i,k}$ with (20) that also employs the most recent local information $x_{i,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 $\lambda_{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.

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 $f_{i}(x_{i},\mathbf{x}_{-i})$ and local feasible set $\Omega_{i}$.

ii). The matrix $A$ is decomposed and each block $A_{i}$ is kept by player $i$, hence the privacy of each player is protected because $A_{i}$ 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 $\mathcal{G}_{f}$ and $\mathcal{G}_{\lambda}$, 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 $\mathcal{G}_{f}$ and $\mathcal{G}_{\lambda}$, 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 $\lambda_{i}$. Moreover, (18) and (19) can even be computed in parallel for player $i$.

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 $VI(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 $\mathbf{x}_{k}=col(x_{1,k},\cdots,x_{N,k})\in\mathbf{R}^{n}$, $\bar{\lambda}_{k}=col(\lambda_{1,k},\cdots,\lambda_{N,k})\in\mathbf{R}^{mN}$, $\bar{z}_{k}=col(z_{1,k},\cdots,z_{N,k})\in\mathbf{R}^{mN}$, $\bar{b}=col(b_{1},\cdots,b_{N})\in\mathbf{R}^{mN}$, $\bar{\tau}=diag\{\tau_{1}I_{n_{1}},\cdots,\tau_{N}I_{n_{N}}\}\in\mathbf{R}^{n\times n}$, $\bar{\nu}=diag\{\nu_{1}I_{m},\cdots,\nu_{N}I_{m}\}\in\mathbf{R}^{mN\times mN}$, $\bar{\sigma}=diag\{\sigma_{1}I_{m},\cdots,\sigma_{N}I_{m}\}\in\mathbf{R}^{mN\times mN}$, $\Lambda=diag\{A_{1},\cdots,A_{N}\}\in\mathbf{R}^{mN\times n}$ and $\Lambda^{T}=diag\{A_{1}^{T},\cdots,A_{N}^{T}\}\in\mathbf{R}^{n\times mN}$. $\bar{L}=L\otimes I_{m}$ where $L\in\mathbf{R}^{N\times N}$ is the weighted Laplacian matrix of multiplier graph $\mathcal{G}_{\lambda}$.

Using these notations, the definition of pseudo-gradient $F(\mathbf{x})$ in (13) and $P_{\prod_{i=1}^{N}\Omega_{i}}(col(x_{1},\cdots,x_{N}))=col(P_{\Omega_{1}}(x_{1}),\cdots,P_{\Omega_{N}}(x_{N}))$***Proposition 23.16 of Bauschke & Combettes (2011), Algorithm 3.5 can be written in a compact form as:

Algorithm 4.6

**

[TABLE]

**

Using the fact that $P_{\Omega}(x)=Prox_{\iota_{\Omega}}(x)=R_{N_{\Omega}}(x)$ in (4) and the definition of resolvent operator as $R_{N_{\Omega}}(x)=({\rm Id}+N_{\Omega})^{-1}$, equation (21) be be written as $\mathbf{x}_{k+1}=({\rm Id}+N_{\Omega})^{-1}[\mathbf{x}_{k}-\bar{\tau}(F(\mathbf{x}_{k})-\Lambda^{T}\bar{\lambda}_{k})]$, or equivalently,

[TABLE]

Notice that $\bar{\tau}^{-1}=diag\{\frac{1}{\tau_{1}}I_{n_{1}},\cdots,\frac{1}{\tau_{N}}I_{n_{N}}\}$, $\bar{\nu}^{-1}=diag\{\frac{1}{\nu_{1}}I_{m},\cdots,\frac{1}{\nu_{N}}I_{m}\}$ and $\bar{\sigma}^{-1}=diag\{\frac{1}{\sigma_{1}}I_{m},\cdots,\frac{1}{\sigma_{N}}I_{m}\}$. Furthermore, $N_{\Omega}$ is a cone and $N_{\Omega}(\mathbf{x})=\prod_{i=1}^{N}N_{\Omega_{i}}(x_{i})$, hence $\bar{\tau}^{-1}N_{\Omega}(\mathbf{x})=N_{\Omega}(\mathbf{x})$. Therefore, (24) can be written as

[TABLE]

Moreover, $N_{\mathbf{R}^{mN}_{+}}$ is a cone, $N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda})=\prod_{i=1}^{N}N_{\mathbf{R}^{m}_{+}}(\lambda_{i})$, and $\bar{\sigma}^{-1}N_{\mathbf{R}^{mN}_{+}}=N_{\mathbf{R}^{mN}_{+}}$. Then with similar arguments, equation (23) can be written as:

[TABLE]

or equivalently,

[TABLE]

Therefore, equation(22) together with (25) and (27) can be written in a compact form as:

[TABLE]

Denote

[TABLE]

Notice that matrix $\Phi\in\mathbf{R}^{(n+2mN)\times(n+2mN)}$ is symmetric due to $L=L^{T}$ and $\bar{L}=L\otimes I_{m}$.

Denote $\bar{\Omega}=\Omega\times\mathbf{R}^{mN}\times\mathbf{R}_{+}^{mN}$. Define the operators $\bar{\mathfrak{A}}:\bar{\Omega}\rightarrow\mathbf{R}^{n+2mN}$ and $\bar{\mathfrak{B}}:\bar{\Omega}\rightarrow 2^{\mathbf{R}^{n+2mN}}$ as follows,

[TABLE]

Remark 4.7

Operators $\bar{\mathfrak{A}}$ and $\bar{\mathfrak{B}}$ in (49) can be regarded as an extension of operators $\mathfrak{A}$ and $\mathfrak{B}$ in (17) by augmenting $\lambda\in\mathbf{R}^{m}$ to $\bar{\lambda}\in\mathbf{R}^{mN}$ and introducing auxiliary variables $\bar{z}\in\mathbf{R}^{mN}$. Moreover, operators in (49) utilize $L\mathbf{1}_{N}=\mathbf{0}_{N}$ of Laplacian matrix $L$ to ensure the consensus of $\lambda_{i}$, and utilize $\mathbf{1}^{T}_{N}L=\mathbf{0}_{N}^{T}$ to ensure the feasibility of affine coupling constraints.

The next result shows that Algorithm 3.5, or equivalently Algorithm 4.6, 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.8

Suppose $\Phi$ in (42) is positive definite, and operators $\bar{\mathfrak{A}}$ and $\bar{\mathfrak{B}}$ in (49) are maximally monotone. Denote $T_{1}:={\rm Id}-\Phi^{-1}\bar{\mathfrak{A}}$ and $T_{2}:=({\rm Id}+\Phi^{-1}\bar{\mathfrak{B}})^{-1}$. Then any limiting point of Algorithm 3.5, i.e., $col(\mathbf{x}^{*},\bar{z}^{*},\bar{\lambda}^{*})$, is a zero of $\bar{\mathfrak{A}}+\bar{\mathfrak{B}}$ and is a fixed point of $T_{2}\circ T_{1}$.

Proof: Denote $\varpi=col(\mathbf{x},\bar{z},\bar{\lambda})$, then using (34), (42) and (49), Algorithm 3.5 can written in a compact form as follows:

[TABLE]

Since $\Phi$ is symmetric and positive definite, we can write equation (63) as $\varpi_{k}-\Phi^{-1}\bar{\mathfrak{A}}(\varpi_{k})\in\varpi_{k+1}+\Phi^{-1}\bar{\mathfrak{B}}(\varpi_{k+1})$, or equivalently,

[TABLE]

Since $\Phi^{-1}\bar{\mathfrak{B}}$ is maximally monotone (refer to Lemma 5.16), $({\rm Id}+\Phi^{-1}\bar{\mathfrak{B}})^{-1}$ is single-valued. Then by the definition of the inverse of a set-valued operator, Algorithm 3.5 is written as

[TABLE]

Suppose that Algorithm 3.5, or equivalently (65), converges to a limiting point $\varpi^{*}$. 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), $\varpi^{*}=T_{2}T_{1}\varpi^{*}$. Therefore, any limiting point of Algorithm 3.5 is a fixed point of the composition $T_{2}\circ T_{1}$, and Algorithm 3.5 can be regarded as an iterative computation of fixed points of $T_{2}\circ T_{1}$.

By Theorem 25.8 of Bauschke & Combettes (2011), (65) is also the forward-backward splitting algorithm for finding zeros of a sum of monotone operators, hence $\varpi^{*}\in zer(\Phi^{-1}\bar{\mathfrak{A}}+\Phi^{-1}\bar{\mathfrak{B}})$ for any limiting point $\varpi^{*}$. Since $\Phi$ is positive definite, any limiting point $\varpi^{*}=col(\mathbf{x}^{*},\bar{z}^{*},\bar{\lambda}^{*})$, i.e., $\varpi^{*}=T_{2}T_{1}\varpi^{*}$, also belongs to $zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$. In fact,

[TABLE]

$\Box$

Remark 4.9

The iteration $\varpi_{k+1}=T_{2}T_{1}\varpi_{k}$ is also known as Picard iteration for iteratively approximating fixed points of $T_{2}T_{1}$ (refer to Berinde (2007)). Lemma 5.15 will give a sufficient condition for $\Phi$ to be positive definite. Lemma 5.14 will give the condition that ensures $\bar{\mathfrak{A}}$ and $\bar{\mathfrak{B}}$ to be maximally monotone.

The next result shows that any limiting point of Algorithm 3.5, that is, any zero point of operator $\bar{\mathfrak{A}}+\bar{\mathfrak{B}}$, is a variational GNE of game (7).

Theorem 4.10

Suppose that Assumptions 1-3 hold. Consider operators $\bar{\mathfrak{A}}$ and $\bar{\mathfrak{B}}$ defined in (49), and operators $\mathfrak{A}$ and $\mathfrak{B}$ defined in (17). Then the following statements hold:

(i): Given any $col(\mathbf{x}^{*},\bar{z}^{*},\bar{\lambda}^{*})\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, then $\mathbf{x}^{*}$ solves the $VI(F,X)$ in (14), hence $\mathbf{x}^{*}$ is a variational GNE of game in (7). Moreover $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$, and the multiplier $\lambda^{*}$ together with $\mathbf{x}^{*}$ satisfy the KKT condition in (16), i.e., $col(\mathbf{x}^{*},\lambda^{*})\in zer(\mathfrak{A}+\mathfrak{B})$.

(ii): $zer(\mathfrak{A}+\mathfrak{B})\neq\emptyset$ and $zer(\bar{{\mathfrak{A}}}+\bar{{\mathfrak{B}}})\neq\emptyset$.

Proof: (i): By the definition of operators $\bar{\mathfrak{A}}$, $\bar{\mathfrak{B}}$ in (49), we have,

[TABLE]

Suppose that $col(\mathbf{x}^{*},\bar{z}^{*},\bar{\lambda}^{*})\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$. From the second line of (66),

[TABLE]

It follows that $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*},\lambda^{*}\in\mathbf{R}^{m}$ since $L$ is the weighted Laplacian of multiplier graph $\mathcal{G}_{\lambda}$ and $\mathcal{G}_{\lambda}$ is connected due to Assumption 3.

Then by the first line of (66), combined with $\Lambda^{T}=diag\{A_{1}^{T},\cdots,A_{N}^{T}\}$ and $\lambda^{*}_{1}=\lambda^{*}_{2}=,\cdots,=\lambda_{N}^{*}=\lambda^{*}$, we have

[TABLE]

or equivalently,

[TABLE]

By the third line of (66) and using $\bar{L}\bar{\lambda}^{*}=\mathbf{0}$, it follows that

[TABLE]

This implies that there exist $v_{1},v_{2},\cdots,v_{N}\in N_{R^{m}_{+}}(\lambda^{*})$, such that

[TABLE]

Multiplying both sides of above equation with $\mathbf{1}^{T}_{N}\otimes I_{m}$ and combining with $\mathbf{1}^{T}L=\mathbf{0}^{T}$, we have

[TABLE]

Due to the fact that $N_{\bigcap_{i=1}^{N}\Omega_{i}}=\sum_{i=1}^{N}N_{\Omega_{i}}$ if $\bigcap_{i=1}^{N}int(\Omega_{i})\neq\emptyset$ ***Corollary 16.39 of Bauschke & Combettes (2011), we have

[TABLE]

By (68) and (69), for any $col(\mathbf{x}^{*},\bar{\lambda}^{*},\bar{z}^{*})\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, the KKT condition for $VI(F,X)$ in (14), i.e. (16), is satisfied for $\mathbf{x}^{*}$, $\lambda^{*}$. We conclude that $\mathbf{x}^{*}$ solves $VI(F,X)$ in (14), and is a variational GNE of game (7) by Theorem *. It also follows that $\lambda^{*}$ together with $\mathbf{x}^{*}$ satisfy the KKT condition in (16). This also implies $col(\mathbf{x}^{*},\lambda^{*})\in zer(\mathfrak{\mathfrak{A}}+\mathfrak{\mathfrak{B}})$.

(ii) Under Assumptions 1 and 2, the considered game in (7) has a unique variational GNE $\mathbf{x}^{*}$, and there exists $\lambda^{*}\in\mathbf{R}^{m}$ such that the KKT condition (16) is satisfied, i.e. $col(\mathbf{x}^{*},\lambda^{*})\in zer(\mathfrak{A}+\mathfrak{B})$. Therefore $zer(\mathfrak{A}+\mathfrak{B})\neq\emptyset$.

Then we need to show that there exists $col(\mathbf{x}^{*},\bar{\lambda}^{*},\bar{z}^{*})$ such that $col(\mathbf{x}^{*},\bar{\lambda}^{*},\bar{z}^{*})\in zer(\bar{{\mathfrak{A}}}+\bar{{\mathfrak{B}}})$.

Take $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$. Because $L\mathbf{1}_{N}=\mathbf{0}$, the second line of (66) is satisfied.

Since $col(\mathbf{x}^{*},\lambda^{*})\in zer(\mathfrak{A}+\mathfrak{B})$, $\mathbf{0}\in F(x^{*})-A^{T}\lambda^{*}+N_{\Omega}(\mathbf{x}^{*})$. Using $\lambda^{*}_{1}=\lambda^{*}_{2}=,\cdots,=\lambda_{N}^{*}=\lambda^{*}$, $\Lambda^{T}\bar{\lambda}^{*}=col(A_{1}^{T}\lambda^{*},\cdots,A_{N}^{T}\lambda^{*})=A^{T}\lambda^{*}$. Therefore, the first line of (66) is satisfied with $\mathbf{x}^{*},\mathbf{1}_{N}\otimes\lambda^{*}$.

Moreover, with $col(\mathbf{x}^{*},\lambda^{*})\in zer(\mathfrak{A}+\mathfrak{B})$, $\mathbf{0}_{m}\in A\mathbf{x}^{*}-b+N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$. Then we need to show that there exists $\bar{z}^{*}=col(z_{1}^{*},\cdots,z_{N}^{*})\in\mathbf{R}^{mN}$, such that the third line of (66) is satisfied. Take $v^{*}\in N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$ such that $\mathbf{0}=A\mathbf{x}^{*}-b+v^{*}$. Since $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$ and $N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda}^{*})=\prod_{i=1}^{N}N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$, take $v_{1}^{*}=v_{2}^{*}=,\cdots,=v_{N}^{*}=\frac{1}{N}v^{*}\in N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$. Then $(\mathbf{1}^{T}_{N}\otimes I_{m})(\Lambda\mathbf{x}^{*}-\bar{b}+\bar{L}\bar{\lambda}^{*}+col(v_{1}^{*},\cdots,v^{*}_{N}))=\sum_{i=1}^{N}A_{i}x_{i}^{*}-\sum_{i=1}^{N}b_{i}+v^{*}=A\mathbf{x}^{*}-b+v^{*}=\mathbf{0}_{m}$. That is $\Lambda\mathbf{x}^{*}-\bar{b}+\bar{L}\bar{\lambda}^{*}+col(v_{1}^{*},\cdots,v^{*}_{N})\in Null(\mathbf{1}^{T}_{N}\otimes I_{m})$. By the fundamental theorem of linear algebra***Page 405 of Meyer (2000), $Null(\mathbf{1}^{T}_{N}\otimes I_{m})=Range(\mathbf{1}_{N}\otimes I_{m})^{\bot}$ and $Range(\bar{L})=Null(\bar{L})^{\bot}$ since $\bar{L}$ is also symmetric. Notice that $Null(\bar{L})=Range(\mathbf{1}_{N}\otimes I_{m})$, hence $\Lambda\mathbf{x}^{*}-\bar{b}+\bar{L}\bar{\lambda}^{*}+col(v_{1}^{*},\cdots,v^{*}_{N})\in Range(\bar{L})$. Noticing that $col(v_{1}^{*},\cdots,v^{*}_{N})\in N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda}^{*})$, there exists $\bar{z}^{*}\in\mathbf{R}^{mN}$ such that the third line of (66) is satisfied with $\mathbf{x}^{*},\bar{z}^{*},\mathbf{1}_{N}\otimes\lambda^{*}$. Therefore, $zer(\bar{{\mathfrak{A}}}+\bar{{\mathfrak{B}}})\neq\emptyset$. $\Box$

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 (65). We will first show that all the prerequisites in Lemma 4.8 can be satisfied under suitable step-sizes. Then (65), 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 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 operators***Chapter 4 and Chapter 20 of Bauschke & Combettes (2011). All the following results are valid in Hilbert spaces, thus they hold in Euclidean spaces with any $G-$matrix induced norm $||\cdot||_{G}$, given $G$ as a symmetric positive definite matrix. Denote by $||\cdot||$ an arbitrary matrix induced norm in a finite dimensional Euclidean space.

An operator $T:\Omega\subset\mathbf{R}^{m}\rightarrow\mathbf{R}^{m}$ is nonexpansive if it is $1-$Lipschitzian, i.e., $||T(x)-T(y)||\leq||x-y||,\forall x,y\in\Omega.$ An operator $T$ is $\alpha-$averaged if there exists a nonexpansive operator $T^{{}^{\prime}}$ such that $T=(1-\alpha){\rm Id}+\alpha T^{{}^{\prime}}$. Denote the class of $\alpha-$averaged operators as $\mathcal{A}(\alpha)$. If $T\in\mathcal{A}(\frac{1}{2})$, then $T$ is also called firmly nonexpansive operator.

Lemma 5.11

***Proposition 4.25 of Bauschke & Combettes (2011)

Given an operator $T:\Omega\subset\mathbf{R}^{m}\rightarrow\mathbf{R}^{m}$ and $\alpha\in(0,1)$, the following three statements are equivalent:

(i): $T\in\mathcal{A}(\alpha)$;

(ii): $||Tx-Ty||^{2}\leq||x-y||^{2}-\frac{1-\alpha}{\alpha}||(x-y)-(Tx-Ty)||^{2},\forall x,y\in\Omega$

(iii): $||Tx-Ty||^{2}+(1-2\alpha)||x-y||^{2}\leq 2(1-\alpha)\langle x-y,Tx-Ty\rangle,\forall x,y\in\Omega$.

By (iii) of Lemma (*), $T\in\mathcal{A}(\frac{1}{2})$ if and only if

[TABLE]

The operator $T$ is called $\beta-$cocoercive (or $\beta-$inverse strongly monotone) if $\beta T$ is firmly nonexpansive, i.e.,

[TABLE]

Lemma 5.12

***Theorem 18.15 of Bauschke & Combettes (2011)

For a convex differentiable function $f$ with $\vartheta-$ Lipschitzian gradient, we have $\nabla f$ to be $\frac{1}{\vartheta}-$cocoercive, i.e.,

[TABLE]

Lemma * is known as Baillon-Haddad theorem, and one elementary proof can be found in Theorem 2.1.5 of Nesterov (2013).

Lemma 5.13

***Proposition 23.7 of Bauschke & Combettes (2011)

If operator $\Delta$ is maximally monotone, then $T=R_{\Delta}=({\rm Id}+\Delta)^{-1}$ is firmly nonexpasive and $2R_{\Delta}-{\rm Id}$ is nonexpansive.

Hence, the projection operator onto a closed convex set is firmly nonexpansive since $P_{\Omega}=Prox_{\iota_{\Omega}}=R_{\partial\iota_{\Omega}}=R_{N_{\Omega}}$ and $N_{\Omega}$ is maximally monotone***Example 20.41 and Proposition 4.8 of Bauschke & Combettes (2011).

In the following, we analyze the maximal monotonicity of operators $\bar{\mathfrak{A}}$, $\bar{\mathfrak{B}}$ in (49), the positive definite property of matrix $\Phi$, and the properties of operators $T_{1}$ and $T_{2}$ defined in Lemma 4.8 by giving sufficient step-sizes choice conditions, which are shown in Lemma 5.14, 5.15 and 5.16, respectively.

Lemma 5.14

Suppose Assumptions 1-3 hold. Given an Euclidean space with norm $||\cdot||_{2}$, then

(i): Operator $\bar{\mathfrak{A}}$ in (49) is $\beta-$cocoercive with $0<\beta\leq\min\{\frac{1}{2d^{*}},\frac{\eta}{\theta^{2}}\}$ where $d^{*}$ is the maximal weighted degree of multiplier graph $\mathcal{G}_{\lambda}$, i.e., $d^{*}=\max\{\sum_{j=1}^{N}w_{1j},\cdots,\sum_{j=1}^{N}w_{Nj}\}$, and $\eta,\theta$ are parameters in Assumption 2;

(ii): Operator $\bar{\mathfrak{B}}$ in (49) is maximally monotone.

Proof: (i): According to the definition of $\bar{\mathfrak{A}}$ in (49) and the definition of $\beta-$cocercive in (71), we need to prove that

[TABLE]

Notice that $\bar{L}\bar{\lambda}-\bar{b}$ is the gradient of function $\tilde{f}(\bar{\lambda}):=\frac{1}{2}\bar{\lambda}^{T}\bar{L}\bar{\lambda}-\bar{b}^{T}\bar{\lambda}$. Moveover, $\tilde{f}(\bar{\lambda})$ is a convex function since $\nabla^{2}\tilde{f}(\bar{\lambda})=\bar{L}$ is positive semi-definite due to Assumption 3 ***Proposition 17.10 of Bauschke & Combettes (2011). It can easily be verified that $\bar{L}\bar{\lambda}-\bar{b}$ is $||L||_{2}-$Lipschitz continuous (notice that the eigenvalues of $\bar{L}$ are just the elements in $col(0,s_{2}\cdots,s_{N})\otimes\mathbf{1}_{m}$), therefore, by Lemma *

[TABLE]

Since $||L||_{2}\leq s_{N}$ and $d^{*}\leq s_{N}\leq 2d^{*}$ by (6) where $d^{*}=\max\{d_{1},\cdots,d_{N}\}$ is the maximal weighted degree of multiplier graph $\mathcal{G}_{\lambda}$, we have $\frac{1}{||L||_{2}}\geq\frac{1}{2d^{*}}$.

Meanwhile by Assumption 2, $F(\mathbf{x})$ is $\eta-$strongly monotone and $\theta-$ Lipschitz continuous over $\Omega$. By $||F(\mathbf{x}_{1})-F(\mathbf{x}_{2})||_{2}^{2}\leq\theta^{2}||\mathbf{x}_{1}-\mathbf{x}_{2}||_{2}^{2}$, $\forall\;\mathbf{x}_{1},\mathbf{x}_{2}\in\Omega$ we have

[TABLE]

Taking $0<\beta\leq\min\{\frac{1}{2d^{*}},\frac{\eta}{\theta^{2}}\}$ and adding (74) and (75) yields (73). Thus operator $\bar{\mathfrak{A}}$ is $\beta-$cocoercive. By the definition of $\beta-$cocoercive in (71), operator $\bar{\mathfrak{A}}$ is also monotone. Since operator $\bar{\mathfrak{A}}$ is also single-valued, it is also maximally monotone.

(ii): The operator $\bar{\mathfrak{B}}$ in (49) can be written as:

[TABLE]

Since $\bar{L}$ is symmetric, $\mathfrak{B}_{1}$ is a skew-symmetric matrix, i.e., $\mathfrak{B}^{T}_{1}=-\mathfrak{B}_{1}$. Hence, $\mathfrak{B}_{1}$ is maximally monotone***Example 20.30 of Bauschke & Combettes (2011).

$\mathfrak{B}_{2}$ can be written as the direct sum of $N_{\Omega}\times\mathbf{0}_{mN}\times N_{\mathbf{R}^{mN}_{+}}$. Both $N_{\Omega}$ and $N_{\mathbf{R}^{mN}_{+}}$ are maximally monotone as normal cones of closed convex sets. Obviously, $\mathbf{0}_{mN}$ is also maximally monotone as a single-valued operator. Furthermore, the direct sum of maximally monotone operators is also maximally monotone***Proposition 20.23 of Bauschke & Combettes (2011), hence $\mathfrak{B}_{2}$ is also maximally monotone.

Obviously, $dom\mathfrak{B}_{1}=\mathbf{R}^{n+2mN}$, hence $\bar{\mathfrak{B}}=\mathfrak{B}_{1}+\mathfrak{B}_{2}$ is also maximally monotone*** Corollary 24.4 of Bauschke & Combettes (2011). $\Box$

Lemma 5.15

Given any $\delta>0$, if each player $i$ takes $\tau_{i}>0$, $\nu_{i}>0$ and $\sigma_{i}>0$ as its local fixed step-sizes in Algorithm 3.5 that satisfy:

[TABLE]

then matrix $\Phi$ defined in (42) is positive definite, and $\Phi-\delta I_{n+2mN}$ is positive semi-definite.

Proof: It is sufficient to show that $\Phi-\delta I_{n+2mN}$ is positive semi-definite.

[TABLE]

where $\bar{\tau}^{-1}-\delta I_{n}=diag\{(\frac{1}{\tau_{1}}-\delta)I_{n_{1}},\cdots,(\frac{1}{\tau_{N}}-\delta)I_{n_{N}}\}$, $\bar{\nu}^{-1}-\delta I_{mN}=diag\{(\frac{1}{\nu_{1}}-\delta)I_{m},\cdots,(\frac{1}{\nu_{N}}-\delta)I_{m}\}$, and $\bar{\sigma}^{-1}-\delta I_{mN}=diag\{(\frac{1}{\sigma_{1}}-\delta)I_{m},\cdots,(\frac{1}{\sigma_{N}}-\delta)I_{m}\}$. One sufficient condition for matrix $\Phi-\delta I_{n+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,

[TABLE]

It can easily be verified that if each agent chooses its local step-sizes satisfying (76), then (81) is satisfied. $\Box$

Given a globally known parameter $\delta$, each agent can independently choose its local step sizes $\tau_{i}$ ,$\nu_{i}$, and $\sigma_{i}$ with the rule given in (76).

Suppose that the step-sizes $\tau_{i},\mu_{i},\sigma_{i}$ in Algorithm 3.5 are chosen such that $\Phi$ in (42) is positive definite. Thus we can define a norm induced by matrix $\Phi$ for a finite Euclidean space as $||x||_{\Phi}=\sqrt{\langle x,x\rangle_{\Phi}}=\sqrt{\langle\Phi x,x\rangle}$. The next result investigates the properties of operators $\Phi^{-1}\bar{\mathfrak{A}},\Phi^{-1}\bar{\mathfrak{B}}$ and $T_{1},T_{2}$ defined in Lemma 4.8 under $\Phi-$induced norm $||\cdot||_{\Phi}$.

Lemma 5.16

Suppose Assumptions 1-3 hold. Take $0<\beta\leq\min\{\frac{1}{2d^{*}},\frac{\eta}{\theta^{2}}\}$ where $d^{*}$ is the maximal weighted degree of multiplier graph $\mathcal{G}_{\lambda}$, and $\eta,\theta$ are parameters in Assumption 2. Take $\delta>\frac{1}{2\beta}$. Suppose that the step-sizes $\tau_{i},\nu_{i},\sigma_{i}$ in Algorithm 3.5 are chosen to satisfy (76). Then the operators $\Phi^{-1}\bar{\mathfrak{A}}$ and $\Phi^{-1}\bar{\mathfrak{B}}$ with $\Phi$ in (42) and $\bar{\mathfrak{A}},\bar{\mathfrak{B}}$ in (49), and operators $T_{1}={\rm Id}-\Phi^{-1}\bar{\mathfrak{A}},T_{2}=R_{\Phi^{-1}\bar{\mathfrak{B}}}=({\rm Id}+\Phi^{-1}\bar{\mathfrak{B}})^{-1}$ defined as Lemma 4.8 satisfy the following properties under the $\Phi-$induced norm $||\cdot||_{\Phi}$:

(i). $\Phi^{-1}\bar{\mathfrak{A}}$ is $\beta\delta-$cocoercive, and $T_{1}\in\mathcal{A}(\frac{1}{2\delta\beta})$.

(ii). $\Phi^{-1}\bar{\mathfrak{B}}$ is maximally monotone, and $T_{2}\in\mathcal{A}(\frac{1}{2})$.

Proof: (i): By the definition of cocoercivity in (71), we need to prove $\langle\Phi^{-1}\bar{\mathfrak{A}}(x)-\Phi^{-1}\bar{\mathfrak{A}}(y),x-y\rangle_{\Phi}\geq{\beta}{\delta}||\Phi^{-1}\bar{\mathfrak{A}}(x)-\Phi^{-1}\bar{\mathfrak{A}}(y)||^{2}_{\Phi}$, $\forall x,y\in\bar{\Omega}$. Noticing that by the choice of parameters $\tau_{i},\nu_{i},\sigma_{i}$, we have that matrix $\Phi-\delta I_{n+2mN}$ is positive semi-definite from Lemma 5.15. Denote $s_{max}(\Phi)$ and $s_{min}(\Phi)$ as the maximal and minimal eigenvalues of matrix $\Phi$. It must hold that $s_{\max}(\Phi)\geq s_{\min}(\Phi)\geq\delta$. Furthermore, $||\Phi||_{2}=s_{max}(\Phi)\geq s_{min}(\Phi)=\frac{1}{||\Phi^{-1}||_{2}}\geq\delta$ ***Proposition 5.2.7 and 5.2.8 in Meyer (2000), therefore, we also have $||\Phi^{-1}||_{2}\leq\frac{1}{\delta}$. Notice that the operator $\bar{\mathfrak{A}}$ is single-valued and $\Phi^{-1}$ is also nonsingular, so that for any $x,y\in\bar{\Omega}$,

[TABLE]

By the $\beta-$cocoercive property of $\bar{\mathfrak{A}}$ in Lemma 5.14 and the above inequality,

[TABLE]

Therefore, the operator $\Phi^{-1}\bar{\mathfrak{A}}$ is $\beta\delta-$cocoercive under the $\Phi-$induced norm $||\cdot||_{\Phi}$.

Moreover, $\beta\delta\Phi^{-1}\bar{\mathfrak{A}}$ is firmly nonexpansive by the definition of cocoercive operator. This implies that there exists a nonexpansive operator $\breve{T}$ such that $\beta\delta\Phi^{-1}\bar{\mathfrak{A}}=\frac{1}{2}\breve{T}+\frac{1}{2}{\rm Id}$. Then

[TABLE]

since $1<2\beta\delta$ by the assumption that $\delta>\frac{1}{2\beta}$ and $-\breve{T}$ is also nonexpansive.

(ii). $\Phi$ is symmetric positive definite and nonsingular. For any $(x,u)\in gra\Phi^{-1}\bar{\mathfrak{B}}$ and $(y,v)\in gra\Phi^{-1}\bar{\mathfrak{B}}$, $\Phi u\in\Phi\Phi^{-1}\bar{\mathfrak{B}}(x)\in\bar{\mathfrak{B}}(x)$ and $\Phi v\in\Phi\Phi^{-1}\bar{\mathfrak{B}}(y)\in\bar{\mathfrak{B}}(y)$. Then $\langle x-y,u-v\rangle_{\Phi}=\langle x-y,\Phi(u-v)\rangle\geq 0,\forall x,y$ since $\bar{\mathfrak{B}}$ is monotone by Lemma 5.14. Therefore, $\Phi^{-1}\bar{\mathfrak{B}}$ is monotone under the $\Phi-$matrix induced product $\langle\cdot,\cdot\rangle_{\Phi}$.

Furthermore, take $(y,v)\in\bar{\Omega}\times\mathbf{R}^{n+2mN}$, and $\langle x-y,u-v\rangle_{\Phi}\geq 0$, for any other $(x,u)\in gra(\Phi^{-1}\bar{\mathfrak{B}})$. For any $(x,\tilde{u})\in gra\bar{\mathfrak{B}}$, we have $(x,\Phi^{-1}\tilde{u})\in gra(\Phi^{-1}\bar{\mathfrak{B}})$. $\langle x-y,\Phi(\Phi^{-1}\tilde{u}-v)\rangle\geq 0$, or equivalently, $\langle x-y,\tilde{u}-\Phi v)\rangle\geq 0$. Since $\bar{\mathfrak{B}}$ is maximally monotone, then $(y,\Phi v)\in gra\bar{\mathfrak{B}}$. We conclude that $v\in\Phi^{-1}\bar{\mathfrak{B}}(y)$ which implies that $\Phi^{-1}\bar{\mathfrak{B}}$ is maximally monotone.

Therefore, by Lemma * $T_{2}=({\rm Id}+\Phi^{-1}\bar{\mathfrak{B}})^{-1}$ is firmly nonexpansive under the $\Phi-$matrix induced norm $||\cdot||_{\Phi}$. $\Box$

Summarizing the above results, take $0<\beta\leq\min\{\frac{1}{2d^{*}},\frac{\eta}{\theta^{2}}\}$, $\delta>\frac{1}{2\beta}$, and $\tau_{i},\nu_{i},\sigma_{i}$ satisfying (76), then $T_{1}:={\rm Id}-\Phi^{-1}\bar{\mathfrak{A}}\in\mathcal{A}(\frac{1}{2\delta\beta})$, and $T_{2}:=({\rm Id}+\Phi^{-1}\bar{\mathfrak{B}})^{-1}\in\mathcal{A}(\frac{1}{2})$ in the Euclidean space with $\Phi-$matrix induced norm $||\cdot||_{\Phi}$. The next result shows the convergence of Algorithm 3.5 based on its compact reformulation (65) and the properties of $T_{1}$ and $T_{2}$.

Theorem 5.17

Suppose Assumptions 1-3 hold. Take $0<\beta\leq\min\{\frac{1}{2d^{*}},\frac{\eta}{\theta^{2}}\}$ where $d^{*}$ is the maximal weighted degree of multiplier graph $\mathcal{G}_{\lambda}$, and $\eta,\theta$ are parameters in Assumption 2. Take $\delta>\frac{1}{2\beta}$. The step-sizes $\tau_{i},\nu_{i},\sigma_{i}$ in Algorithm 3.5 are chosen to satisfy (76). Then with Algorithm 3.5, each player has its local strategy $x_{i,k}$ converging to its corresponding component in the variational GNE of game (7), and the local Lagrangian multipliers $\lambda_{i,k}$ of all the agents converge to the same Lagrangian multiplier corresponding with KKT condition (16), i.e.,

[TABLE]

Proof: With Lemma 5.14 and 5.15, Algorithm 3.5 can be written in a compact form (65) according to Lemma 4.8, i.e., $\varpi_{k+1}=T_{2}T_{1}\varpi_{k}$. The convergence analysis will be conducted via the analysis of this iterative computation of fixed points of $T_{2}\circ T_{1}$.

Firstly, by (i) and (ii) of Lemma * and the fact that $T_{1},T_{2}$ are averaged operators due to Lemma 5.16, $T_{1},T_{2}$ are also nonexpansive operators under the $\Phi-$matrix induced norm $||\cdot||_{\Phi}$. Take any $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, or equivalently any fixed point of $T_{2}\circ T_{1}$, i.e., $\varpi^{*}=T_{2}T_{1}\varpi^{*}$, and then by Lemma 4.8 and (65),

[TABLE]

Hence the sequence $\{||\varpi_{k}-\varpi^{*}||_{\Phi}\}$ is non-increasing and bounded from below. By the monotonic convergence theorem, $\{||\varpi_{k}-\varpi^{*}||\}$ is bounded and converges for every $\varpi^{*}\in zer(\mathfrak{A}+\mathfrak{B})$.

By Lemma 5.16, $T_{1}\in\mathcal{A}(\frac{1}{2\delta\beta})$ and $T_{2}\in\mathcal{A}(\frac{1}{2})$. Denote $\xi=\frac{1}{2\delta\beta}\in(0,1)$. Then with (ii) of Lemma * and (65) we have,

[TABLE]

where the first inequality follows by $T_{2}\in\mathcal{A}(\frac{1}{2})$ and the second inequality follows by $T_{1}\in\mathcal{A}(\frac{1}{2\delta\beta})$, both utilizing (ii) of Lemma *. Notice that

[TABLE]

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

[TABLE]

where the second equality follows from (86) by setting $\alpha=\xi$, $x=(T_{1}\varpi_{k}-T_{1}\varpi^{*})-(T_{2}T_{1}\varpi_{k}-T_{2}T_{1}\varpi^{*})$ and $y=(T_{1}\varpi_{k}-T_{1}\varpi^{*})-(\varpi_{k}-\varpi^{*})$.

Combining (85) and (LABEL:equ_thm_57_2) yields $\forall k\geq 0$,

[TABLE]

Using (88) from [math] to $k$ and adding all $k+1$ inequalities yields

[TABLE]

Taking limit as $k\rightarrow\infty$ we have, $(1-\xi)\sum_{k=1}^{\infty}||\varpi_{k}-\varpi_{k+1}||^{2}_{\Phi}\leq||\varpi_{0}-\varpi^{*}||^{2}_{\Phi}$. Since $1-\xi>0$, it follows that $\sum_{k=1}^{\infty}||\varpi_{k}-\varpi_{k+1}||^{2}_{\Phi}$ converges and $\lim_{k\rightarrow\infty}\varpi_{k}-\varpi_{k+1}=0$.

Since $\{||\varpi_{k}-\varpi^{*}||\}$ is bounded and converges, $\{\varpi_{k}\}$ is a bounded sequence. There exists a subsequence$\{\varpi_{n_{k}}\}$ that converges to $\tilde{\varpi}^{*}$. Notice that the composition $T_{2}\circ T_{1}$ is (Lipschitz) continuous and single-valued, because (65) is just an equivalent expression of Algorithm 3.5, and obviously the right hand side of Algorithm 3.5 is continuous. ${\varpi}_{n_{k}+1}=T_{2}T_{1}\varpi_{n_{k}}$. Since $T_{2}T_{1}$ is continuous, and $\lim_{n_{k}\rightarrow\infty}\varpi_{n_{k}}-{\varpi}_{n_{k}+1}=0$, passing to limiting point, we have $\tilde{\varpi}^{*}=T_{2}T_{1}\tilde{\varpi}^{*}$. Therefore, the limiting point $\tilde{\varpi}^{*}$ is a fixed point of $T_{2}T_{1}$, or equivalently, $\tilde{\varpi}^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$.

Setting $\varpi^{*}=\tilde{\varpi}^{*}$ in (LABEL:equ_thm_5_2), we have $\{||\varpi_{k}-\tilde{\varpi}^{*}||\}$ is bounded and converges. Since there exists a subsequence $\{\varpi_{n_{k}}\}$ that converges to $\tilde{\varpi}^{*}$, it follows that $\{||\varpi_{k}-\tilde{\varpi}^{*}||\}$ converges to zero. Therefore, $\lim_{k\rightarrow\infty}\varpi_{k}\rightarrow\tilde{\varpi}^{*}$. By Theorem 4.10, this just implies (83). $\Box$

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 the current operator output with previous iterate, leads to the well-known Krasnosel’skiĭ-Mann type of fixed point iteration***Chapter 5 of Bauschke & Combettes (2011), 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 & Vũ (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:

Algorithm 6.18

**

[TABLE]

$\alpha>0$ is a fixed step-size in the acceleration phase, and $\tau_{i}>0,\nu_{i}>0,\sigma_{i}>0$ are fixed step-sizes of player $i$, and $W=[w_{ij}]$ is the weighted adjacency matrix of multiplier graph $\mathcal{G}_{\lambda}$.*

Compared with Algorithm 3.5, Algorithm 6.18 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.18 has only an additional simple local computation burden. Obviously, Algorithm 6.18 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 $\alpha$.

In the following two subsections, we will first give some intuitive interpretation of Algorithm 6.18 from the viewpoint of a discretization of continuous-time dynamical systems, and then prove its convergence.

6.1 Interpretations from viewpoints of dynamical systems

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.18 just for illustrations of the differences behind the algorithms.

Firstly, let us show that Algorithm 3.5, or equivalently its compact reformulation (65) in Lemma 4.8, can be interpreted as the discretization of the following dynamical system:

[TABLE]

In fact, for differential inclusion (95), we have the following implicit/explicit discretization with step-size of $h$,

[TABLE]

Denote $\varpi_{k}=\varpi(kh)$ and take $h=1$, then (96) can be written as

[TABLE]

Therefore, the implicit/explicit discretization of (95) is exactly (64) that leads to (65), 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 $\mathbf{\tilde{x}}=col(\tilde{x}_{1},\cdots,\tilde{x}_{N})$, $\bar{\tilde{\lambda}}=col(\tilde{\lambda}_{1},\cdots,\tilde{\lambda}_{N})$, and $\bar{\tilde{z}}=col(\tilde{z}_{1},\cdots,\tilde{z}_{N})$. And further denote $\tilde{\varpi}=col(\mathbf{\tilde{x}},\bar{\tilde{\lambda}},\bar{\tilde{z}})$. Then by similar arguments as in Section 4 and using operators $\bar{\mathfrak{A}}$ and $\bar{\mathfrak{B}}$ defined in (49) and $\Phi$ defined in (42), Algorithm 6.18 can be written in a compact form (assume $\bar{\mathfrak{B}}$ is maximally monotone),

[TABLE]

where $\varpi_{k}$ is defined as in Section 4.

We can show that Algorithm 6.18, or equivalently (98)-(99), can be interpreted as the discretization of the following second-order continuous-time dynamical system,

[TABLE]

In fact, for differential inclusion (100) consider the following type of implicit/explicit discretization,

[TABLE]

where $\tilde{\varpi}(kh)$ is an interpolation point to be determined later. Denote $\varpi_{k}=\varpi(kh)$ and take $h=1$, then (101) can be written as

[TABLE]

Denote $\alpha=1-\tilde{\alpha}$ and take $\tilde{\varpi}_{k}=\varpi_{k}+\alpha(\varpi_{k}-\varpi_{k-1})$, then (102) can be written as

[TABLE]

(103) leads to equations (98)-(99), or equivalently Algorithm 6.18.

Remark 6.19

Compared with (95), (100) is a second order dynamical system with an additional inertial term $\alpha\dot{\varpi}$, hence (100) enjoys better convergence properties than (95). Therefore, it is expected that Algorithm 6.18, as a discretization of (100), would have better convergence properties than Algorithm 3.5.

6.2 Convergence analysis

The following result proves the convergence of Algorithm 6.18 by providing sufficient step-size choices for $\alpha$ as well as $\tau_{i},\nu_{i},\sigma_{i}$. The sufficient choice condition for $\alpha$ 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 & Vũ (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.20’s proof is also provided for completeness.

Theorem 6.20

Suppose Assumptions 1-3 hold. Take $0<\beta\leq\min\{\frac{1}{2d^{*}},\frac{\eta}{\theta^{2}}\}$ where $d^{*}$ is the maximal weighted degree of multiplier graph $\mathcal{G}_{\lambda}$, and $\eta,\theta$ are parameters in Assumption 2. Given a sufficient small $0<\epsilon<1$, take $\delta>\frac{1}{2\beta}$ and $0<\alpha<1$ in Algorithm 6.18 such that $2\beta\delta(1-3\alpha-\epsilon)\geq(1-\alpha)^{2}$. Suppose that player $i$ chooses its step-sizes $\tau_{i},\nu_{i},\sigma_{i}$ in Algorithm 6.18 satisfying (76). Then with Algorithm 6.18, players’ local strategies converge to the variational GNE of game in (7), and the local multipliers $\lambda_{i,k}$ of all the agents converge to the same multiplier corresponding with KKT condition (16), i.e.,

[TABLE]

Proof: By the choice of $\beta$, operator $\bar{\mathfrak{A}}$ is $\beta-$cocoercive and operator $\bar{\mathfrak{B}}$ is maximally monotone due to Lemma 5.14. By the choice of $\delta,\tau_{i},\nu_{i}$ and $\sigma_{i}$, $\Phi-\delta I_{n+2mN}$ is positive semi-definite due to Lemma 5.15, and $\Phi^{-1}\bar{\mathfrak{A}}$ and $\Phi^{-1}\bar{\mathfrak{B}}$ satisfy the properties in Lemma 5.16. Therefore, Algorithm 6.18 can be exactly written in the compact form of (98)-(99) with similar arguments as in Lemma 4.8.

Resorting to Theorem 4.10, we only need to show that Algorithm 6.18 converges and its limiting point belongs to $zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$. In fact, any limiting point of Algorithm 6.18 satisfies $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$ as shown next. Suppose $\lim_{k\rightarrow\infty}\varpi_{k}\rightarrow\varpi^{*}$, then $\tilde{\varpi}_{k}\rightarrow\varpi^{*}$ and $\varpi^{*}=({\rm Id}+\Phi^{-1}\bar{\mathfrak{B}})^{-1}({\rm Id}-\Phi^{-1}\bar{\mathfrak{A}})\varpi^{*}$ using the continuity of the right hand of (99). Therefore, $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$ because $\Phi$ is a positive definite matrix.

The following relationship (similar with the cosine rule) will be heavily utilized in the convergence analysis.

[TABLE]

which can be verified by directly expanding with $||a+b||^{2}_{Q}=||a||^{2}_{Q}+2\langle a,b\rangle_{Q}+||b||^{2}_{Q}$.

The proof is divided into three parts:

Part 1: Given any $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, $||\varpi_{k+1}-\varpi^{*}||^{2}_{\Phi}-||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}$ follows a recursive inequality as follows,

[TABLE]

• 2.

Part 2: Given any $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, $\sum_{k}^{\infty}||\varpi_{k+1}-\varpi_{k}||^{2}_{\Phi}<\infty$ and $\lim_{k\rightarrow\infty}\varpi_{k+1}-\varpi_{k}=0$.

• 3.

Part 3: We first show the convergence of $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}\}$ given any $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, and then show the convergence of Algorithm 6.18.

Part 1: Given any point $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, we first prove a recursive inequality (106) for $||\varpi_{k+1}-\varpi^{*}||^{2}_{\Phi}-||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}$.

Using (105) to expand the left hand of (106) yields,

[TABLE]

where the last step is derived by incorporating (98).

To tackle the second term on the right hand of (LABEL:equ_thm_6_1), we proceed as follows. By (99) we have, $\tilde{\varpi}_{k}-\Phi^{-1}\bar{\mathfrak{A}}(\tilde{\varpi}_{k})\in\varpi_{k+1}+\Phi^{-1}\bar{\mathfrak{B}}(\varpi_{k+1}),$ or equivalently,

[TABLE]

Because $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, we also have

[TABLE]

Due to the maximal monotonicity of $\bar{\mathfrak{B}}$ proved in Lemma 5.14,

[TABLE]

By incorporating (108) and (109) into (110), we have

[TABLE]

or equivalently,

[TABLE]

Using (111) for the second term on the right hand of (LABEL:equ_thm_6_1) yields

[TABLE]

By Lemma 5.14, $\bar{\mathfrak{A}}$ is $\beta-$cocoercive. For the second term on the right hand of (112), we have

[TABLE]

where the first inequality is obtained by the cocoercive property (71) and the second inequality is obtained by $2\langle a,b\rangle\geq-||a||_{2}^{2}-||b||_{2}^{2}$.

Combining (112) with (LABEL:equ_thm_6_7) we have

[TABLE]

or equivalently we have

[TABLE]

Utilizing the equality (105) we also have,

[TABLE]

Combining (114) and (115), we have,

[TABLE]

Next using (98) and (105) for the first and third terms on the right hand of (LABEL:equ_thm_6_10) yields

[TABLE]

For the second term on the right hand of (LABEL:equ_thm_6_10), $\frac{1}{2\beta}||{\varpi}_{k+1}-\tilde{\varpi}_{k}||_{2}^{2}\leq\frac{1}{2\beta\delta}||{\varpi}_{k+1}-\tilde{\varpi}_{k}||^{2}_{\Phi}$, since $\Phi>\delta I_{n+2mN}$. Incorporating this and (LABEL:equ_thm_6_11) into (LABEL:equ_thm_6_10)

[TABLE]

Since $\delta>\frac{1}{2\beta}$, we derive (106).

Part 2: In this step, we will prove $\sum_{k}^{\infty}||\varpi_{k+1}-\varpi_{k}||^{2}_{\Phi}<\infty$ and $\lim_{k\rightarrow\infty}\varpi_{k+1}-\varpi_{k}=0$.

Denote $S=\Phi-\frac{1}{2\beta}I_{n+2mN}$. Then $S$ is symmetric and positive definite since $\Phi\geq\delta I_{n+2mN}$ and $\delta>\frac{1}{2\beta}$. The first inequality of (LABEL:equ_thm_6_12) can also be written as:

[TABLE]

For the first term on the right hand of (119),

[TABLE]

where the first equality follows from (105), the second equality follows from (98), and the third inequality follows from $-2\langle x,y\rangle\leq||x||^{2}+||y||^{2}$.

Denote $Q=2\Phi-\frac{1-\alpha}{2\beta}I_{n+2mN}$. Then $Q$ is also symmetric and positive definite, since $\alpha<1$ and $\Phi\geq\delta I_{n+2mN}$. Combining (119) with (LABEL:equ_thm_6_15),

[TABLE]

Denote $\mu_{k}=||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\alpha||\varpi_{k-1}-\varpi^{*}||^{2}_{\Phi}+\alpha||\varpi_{k}-\varpi_{k-1}||^{2}_{Q}$, then

[TABLE]

where the third inequality follows by (LABEL:equ_thm_6_16).

Given a sufficient small $0<\epsilon<1$, choose $0<\alpha<1$ and $\delta>\frac{1}{2\beta}$ such that $2\beta\delta(1-3\alpha-\epsilon)\geq(1-\alpha)^{2}$, then

[TABLE]

Therefore, (LABEL:equ_thm_6_17) yields $\mu_{k+1}-\mu_{k}\leq-\epsilon||\varpi_{k+1}-\varpi_{k}||^{2}_{\Phi}$. Therefore, $\mu_{k+1}\leq\mu_{k}\leq\mu_{1}$. By the definition of $\mu_{k}$, $||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\alpha||\varpi_{k-1}-\varpi^{*}||^{2}_{\Phi}\leq||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\alpha||\varpi_{k-1}-\varpi^{*}||^{2}_{\Phi}+\alpha||\varpi_{k}-\varpi_{k-1}||^{2}_{Q}\leq\mu_{1}$. Therefore, $||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}\leq\alpha||\varpi_{k-1}-\varpi^{*}||^{2}_{\Phi}+\mu_{1}$. By Lemma 1 at Page 44 of Polyak (1987), $||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}\leq\alpha^{k}(||\varpi_{1}-\varpi^{*}||^{2}_{\Phi}-\frac{\mu_{1}}{1-a})+\frac{\mu_{1}}{1-\alpha}$, and $\{\varpi_{k}\}$ is bounded sequence.

We also have $\mu_{k+1}-\mu_{1}\leq-\epsilon\sum_{i=1}^{k}||\varpi_{i+1}-\varpi_{i}||^{2}_{\Phi}$. Then $\epsilon\sum_{i=1}^{k}||\varpi_{i+1}-\varpi_{i}||^{2}_{\Phi}\leq\mu_{1}-\mu_{k+1}\leq\mu_{1}+\alpha||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}\leq\mu_{1}+\alpha^{k+1}(||\varpi_{1}-\varpi^{*}||^{2}_{\Phi}-\frac{\mu_{1}}{1-a})+\frac{\alpha\mu_{1}}{1-\alpha}$. Let $k$ goes to infinity, we have,

[TABLE]

Therefore, $\lim_{k\rightarrow\infty}\varpi_{k+1}-\varpi_{k}=0$, and $\sum_{k=1}^{\infty}(\alpha+\alpha^{2})||\varpi_{k}-\varpi_{k-1}||^{2}_{\Phi}<\infty$.

Part 3: In this part, we first show the convergence of $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}\}$ given any $\varpi^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$, and then show the convergence of Algorithm 6.18.

Denote $\phi_{k}=\max\{0,||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-||\varpi_{k-1}-\varpi^{*}||^{2}_{\Phi}\}$ and $\psi_{k}=(\alpha+\alpha^{2})||\varpi_{k}-\varpi_{k-1}||^{2}_{\Phi}$, and recall (106), we have $\phi_{k+1}\leq\alpha\phi_{k}+\psi_{k}.$ Apply this relationship recursively,

[TABLE]

Summing (124) from $k=1$ to $k=J$,

[TABLE]

Let $J\rightarrow\infty$, then since $0<\alpha<1$,

[TABLE]

Noticing that $\sum_{t=1}^{\infty}\psi_{t}=\sum_{k=1}^{\infty}(\alpha+\alpha^{2})||\varpi_{k}-\varpi_{k-1}||^{2}_{\Phi}<\infty$, $\sum_{k=1}^{\infty}\phi_{k}<\infty$, and hence the sequence $\{\sum_{i=1}^{k}\phi_{i}\}$, being a nonnegative and non-decreasing sequence, converges and is bounded.

Consider another sequence $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\sum_{i=1}^{k}\phi_{i}\}$. Since $||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}$ is nonnegative and $\{\sum_{i=1}^{k}\phi_{i}\}$ is bounded, $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\sum_{i=1}^{k}\phi_{i}\}$ is bounded from below. Furthermore, $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\sum_{i=1}^{k}\phi_{i}\}$ is a non-increasing sequence. In fact,

[TABLE]

where the second inequality follows from the definition of $\phi_{k}$, $\phi_{k+1}\geq||\varpi_{k+1}-\varpi^{*}||^{2}_{\Phi}-||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}$. As $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\sum_{i=1}^{k}\phi_{i}\}$ is a non-increasing sequence and bounded from below, $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\sum_{i=1}^{k}\phi_{i}\}$ converges.

Therefore, $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}\}$, being the sum of two convergent sequences $\{||\varpi_{k}-\varpi^{*}||^{2}_{\Phi}-\sum_{i=1}^{k}\phi_{i}\}$ and $\{\sum_{i=1}^{k}\phi_{i}\}$, also converges.

We are ready to show the convergence of Algorithm 6.18 using the results in Part 1 and Part 2. Since $\{\varpi_{k}\}$ is a bounded sequence, it has a convergent subsequence $\{\varpi_{n_{k}}\}$ that converges to $\tilde{\varpi}^{*}$. Because $\lim_{k\rightarrow\infty}\varpi_{k+1}-\varpi_{k}=0$ by Part 2, we have $\lim_{k\rightarrow\infty}\varpi_{n_{k}-1}-\varpi_{n_{k}}=0$ and $\lim_{k\rightarrow\infty}\varpi_{n_{k}+1}-\varpi_{n_{k}}=0$. Pass to limiting point of $\{\varpi_{n_{k}}\}$, then we have $\tilde{\varpi}^{*}=T_{1}T_{2}\tilde{\varpi}^{*}$ because the righthand side of (98)-(99) is continuous. Hence, $\tilde{\varpi}^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$. Taking $\varpi^{*}=\tilde{\varpi}^{*}$ in (106) of Part 1, we also have $\{||\varpi_{k}-\tilde{\varpi}^{*}||^{2}_{\Phi}\}$ converges by Part 3. Because a subsequence $\{||\varpi_{n_{k}}-\tilde{\varpi}^{*}||^{2}_{\Phi}\}$ converges to zero, the whole sequence $\{||\varpi_{k}-\tilde{\varpi}^{*}||^{2}_{\Phi}\}$ converges to zero. Therefore, the whole sequence of $\{\varpi_{k}\}$ converges to $\tilde{\varpi}^{*}\in zer(\bar{\mathfrak{A}}+\bar{\mathfrak{B}})$. Resorting to Theorem 4.10 gives the desired result. $\Box$

Remark 6.21

A sufficient and simple choice of parameters to ensure the conditions in Theorem (6.20) is $\delta=\frac{1}{\beta}$, $\epsilon=\alpha$ and $0<\alpha<\sqrt{10}-3$. In fact, $\delta>\frac{1}{2\beta}$ implies that $2\beta\delta$ could be be any real number $\varrho>1$. If we take $\epsilon=\varsigma\alpha,\varsigma>0$, then the quadratic inequality becomes $\alpha^{2}-(2-3\varrho-\varsigma\varrho)\alpha+1-\varrho<0$. Since $1-\varrho<0$, $\alpha^{2}-(2-3\varrho-\varsigma\varrho)\alpha+1-\varrho$ takes value strictly less than zero when $\alpha$ takes [math]. By the continuity of quadratic equation, there always exists $0<\alpha<1$ that ensures the above quadratic inequality given any $\varrho>1$ and $\varsigma>0$.

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 $F_{1},\cdots,F_{N}$ and $m$ markets with labels $M_{1},\cdots,M_{m}$. Company $F_{i}$ decides its strategy to participate in the competition in $n_{i}$ markets by producing and delivering $x_{i}\in\mathbf{R}^{n_{i}}$ amounts of products to the markets it connects with. The production limitation of company $F_{i}$ is $x_{i}\in\Omega_{i}\subset\mathbf{R}^{n_{i}}$. Company $F_{i}$ has a local matrix $A_{i}\in\mathbf{R}^{m\times n_{i}}$ that specifies which market it will participate in. The $j$-th column of $A_{i}$, that is $[A_{i}]_{:j}$, has only one element being $1$ and all other elements being [math], and $[A_{i}]_{:j}$ has its $k$-th element being $1$ if and only if player $F_{i}$ delivers $[x_{i}]_{j}$ amount of production to the market $M_{k}$. Therefore, matrices $A_{1},\cdots,A_{N}$ can be used to specify a bipartite graph that represents the connections between the companies and the markets. Denote $n=\sum_{i=1}^{N}n_{i}$, $\mathbf{x}=col(x_{1},\cdots,x_{N})\in\mathbf{R}^{n}$, and $A=[A_{1},\cdots,A_{N}]\in\mathbf{R}^{m\times n}$. Then $A\mathbf{x}\in\mathbf{R}^{m}=\sum_{i=1}^{N}A_{i}x_{i}$ is just the total product supply to all the markets given the action profile $\mathbf{x}$ of all the companies. Market $M_{j}$ has a maximal capacity of $r_{j}>0$, therefore, it should be satisfied that $A\mathbf{x}\leq r$ where $r=col(r_{1},\cdots,r_{m})\in\mathbf{R}^{m}$. Suppose that $P:\mathbf{R}^{m}\rightarrow\mathbf{R}^{m}$ 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 $c_{i}(x_{i}):\Omega_{i}\rightarrow\mathbf{R}$. Then the local objective function of company (player) $F_{i}$ is $f_{i}(x_{i},\mathbf{x}_{-i}):c_{i}(x_{i})-P^{T}(A\mathbf{x})A_{i}x_{i}$.

Overall, in this network Cournot game, each company needs to solve the following optimization problem given the other companies’ profile $\mathbf{x}_{-i}$,

[TABLE]

Obviously, the above network Cournot game in (127) 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 (127), 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.22** (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 $\{L_{1},...,L_{m}\}$, 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 $A_{i}\in\mathbf{R}^{m\times n_{i}}$. Each column of $A_{i}$ 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 $L_{k}$ and transfers data rate $[x_{i}]_{j}$ on link $L_{k}$. The local decision variable of $S_{i}$ is the data rate on each link that it utilizes, denoted by $x_{i}\in\mathbf{R}^{n_{i}}$. $x_{i}$ also has a local feasibility constraint $x_{i}\in\Omega_{i}$. Denote $n=\sum_{i=1}^{N}n_{i}$, $\mathbf{x}=col(x_{1},\cdots,x_{N})\in\mathbf{R}^{n}$, $A=[A_{1},\cdots,A_{N}]\in\mathbf{R}^{m\times n}$ and $c=col(c_{1},\cdots,c_{m})\in\mathbf{R}^{m}$. The total data rate on each link should be less than the capacity of that link: $A\mathbf{x}\leq c.$ Given the data rate profile of all the nodes $\mathbf{x}$, the payoff function of $S_{i}$, $J_{i}(x_{i},\mathbf{x}_{-i}):\mathbf{R}^{n}\rightarrow\mathbf{R}$, takes the form as $J_{i}(x_{i},\mathbf{x}_{-i})=-u_{i}(x_{i})+\mathcal{D}^{T}(A\mathbf{x})A_{i}x_{i}$, where $u_{i}(x_{i}):\Omega_{i}\rightarrow\mathbf{R}$ is the utility of source $S_{i}$, and $\mathcal{D}:R^{m}\rightarrow R^{m}$ 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 (127).

Example 7.23** (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 $d_{i}>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 $x_{i}\in\mathbf{R}^{t_{i}}$ as its local load shedding vector in some specific time periods. Each load manager $i$ also has a local matrix $A_{i}\in\mathbf{R}^{T\times t_{i}}$ that specifies which time period player $i$ will participate. For $j-$th column of $A_{i}$, it has one element being $1$ while all other elements being zero. The $k-$th element of the $j-$th column of $A_{i}$ is $1$ if load manager $i$ decides to decrease its load by $[x_{i}]_{j}$ at time $k$. Denote $\mathbf{x}=col(x_{1},\cdots,x_{N})$, and $A=[A_{1},\cdots,A_{N}]\in\mathbf{R}^{T\times\sum_{i=1}^{N}t_{i}}$, $d=col(d_{1},\cdots,d_{T})$. Naturally, it is required that the total load shedding of all the load managers should meet the minimal value, $A\mathbf{x}\geq d$. Each player has a local feasible constraint $x_{i}\in\Omega_{i}\subset\mathbf{R}^{t_{i}}$, and a cost function $c_{i}(x_{i}):\Omega_{i}\rightarrow\mathbf{R}$ due to local load shedding. $P:\mathbf{R}^{T}\rightarrow\mathbf{R}^{T}$ is the payment price vector function that maps total load shedding of each period to the payment price vector, therefore, $P^{T}(Ax)A_{i}x_{i}$ is the payment awards of player $i$ for its load shedding. The disutility function of player $i$ is $J_{i}(x_{i},\mathbf{x}_{-i})=c_{i}(x_{i})-P^{T}(A\mathbf{x})A_{i}x_{i}$ given all the players’ action profile $\mathbf{x}$. All in all, the demand response management game is well described by the network Cournot game model in (127).

Moreover, the Assumptions 1 and 2 for Algorithm 3.5 and 6.18 can easily be satisfied for many practical cost functions and price functions. For example, take company $F_{i}$’s production cost function to be a strongly convex function with Lipschitz continuous gradients (A quadratic function $c_{i}(x_{i})=x_{i}^{T}Q_{i}x_{i}+b_{i}^{T}x_{i}$ with $Q_{i}\in\mathbf{R}^{n_{i}\times n_{i}}$ being a symmetric and positive definite matrix and $b_{i}\in\mathbf{R}^{n_{i}}$ is one possible choice). The price of market $M_{j}$ is taken as the linear function of the total supplying $p_{j}(\mathbf{x})=\bar{P}_{j}-d_{j}[A\mathbf{x}]_{j}$ (known as a linear inverse demand function in economics) with $\bar{P}_{j}>0,d_{j}>0$. Denote $P=col(p_{1},\cdots,p_{m}):\mathbf{R}^{n}\rightarrow\mathbf{R}^{m}$, $\bar{P}=col(\bar{P}_{1},\cdots,\bar{P}_{m})\in\mathbf{R}^{m}$, $D=diag\{d_{1},\cdots,d_{m}\}\in\mathbf{R}^{m\times m}$. Then $P=\bar{P}-DA\mathbf{x}$ is the vector price function. The payments of company $F_{i}$ by selling product $x_{i}$ to the markets that it connects with is just $P^{T}A_{i}x_{i}$. Therefore, the objective function of company $F_{i}$ is,

[TABLE]

Denote $F(\mathbf{x})=col(\nabla_{x_{1}}J_{1}(x_{1},\mathbf{x}_{-1}),\cdots,\nabla_{x_{N}}J_{N}(x_{N},\mathbf{x}_{-N}))$, $\nabla c(\mathbf{x})=col(\nabla c_{1}(x_{1}),\nabla c_{2}(x_{2}),\cdots,\nabla c_{N}(x_{N}))$ and $Q\in\mathbf{R}^{n\times n}$, and

[TABLE]

then

[TABLE]

Notice that $Q$ can be written as

[TABLE]

where $S$ is a block matrix defined as $S=[\sqrt{D}A_{1},\cdots,\sqrt{D}A_{N}]\in\mathbf{R}^{m\times n}$ and $\sqrt{D}=diag\{\sqrt{d}_{1},\cdots,\sqrt{d}_{m}\}$. Therefore, $Q$ is positive semi-definite matrix. Hence, the Jacobian matric of $F(\mathbf{x})$, $JF(\mathbf{x})=diag\{\nabla^{2}c_{1}(x_{1}),\nabla^{2}c_{2}(x_{2}),\cdots,\nabla^{2}c_{N}(x_{N})\}+Q$ is postive definite since the cost functions $c_{i}(x_{i})$ are strongly convex. Therefore, $F(x)$ is strongly monotone***Proposition 2.3.2 of Facchinei and Pang (2007) and Lipschitz continuous, and Assumptions 1 and 2 are satisfied.

Remark 7.24

Notice that in the network Cournot game of (127), each player $i$’s local objective function $c_{i}(x_{i})-P^{T}(A\mathbf{x})A_{i}x_{i}$ only depends on the decisions of the players that participate the same markets as player $i$. Denote $\mathbb{Z}[a,b]$ as the set of integers from $a$ to $b$. Mathematically, $j\in\mathcal{N}_{i}^{f}$ if and only if $\exists k\in\mathbb{Z}[1,m]$, $p\in\mathbb{Z}[1,n_{i}]$ and $q\in\mathbb{Z}[1,n_{j}]$ such that $[A_{i}]_{kp}=1$ and $[A_{j}]_{kq}=1$. Since $A=[A_{1},\cdots,A_{N}]$ is usually a sparse matrix, the interference graph $\mathcal{G}_{f}$ of network Cournot game also has sparse edge connections.

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 $F_{i}$ to $M_{j}$ in Figure 1, then company $F_{i}$ participates the competition in market $M_{j}$ by producing and delivering products to market $M_{j}$. Each company $F_{i}$ has a local constraint as $\mathbf{0}<x_{i}<\Theta_{i}$ where each component of $\Theta_{i}$ is randomly drawn from $(10,25)$. Each market $M_{j}$ has a maximal capacity of $r_{j}>0,j=1,\cdots,7$ and $r_{j}$ is randomly drawn from $(20,80)$. The local objective function is taken as (128). The local cost function of company $i$ is $c_{i}(x_{i})=\pi_{i}(\sum_{j=1}^{n_{i}}[x_{i}]_{j})^{2}+b_{i}^{T}x_{i}$ which is a strongly convex function with Lipschitz continuous gradient. Here $\pi_{i}$ is randomly drawn from $(1,8)$, and each component of $b_{i}$ is randomly drawn from $(1,4)$. The price function is taken as the linear function $P=\bar{P}-DA\mathbf{x}$, and $\bar{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 (128), the interference graph $\mathcal{G}_{f}$ can be easily obtained and is depicted in Figure 2. Meanwhile, we adopte the multiplier graph $\mathcal{G}_{\lambda}$ shown in Figure 3. The weighted adjacency matrix $W=[w_{ij}]$ of multiplier graph $\mathcal{G}_{\lambda}$ has all its nonzero elements to be $1$.

With Figure 1 and the definition of objective function in (128), the interference graph $\mathcal{G}_{f}$ can be easily obtained and is depicted in Figure 2. Meanwhile, we adopte the multiplier graph $\mathcal{G}_{\lambda}$ shown in Figure 3. The weighted adjacency matrix $W=[w_{ij}]$ of multiplier graph $\mathcal{G}_{\lambda}$ has all its nonzero elements to be $1$.

Set the step-sizes in Algorithm 3.5 as $\tau_{i}=0.03,\nu_{i}=0.2,\sigma=0.02$ for all companies, and for Algorithm 6.18 set $\alpha=0.12$ while other step-sizes are the same as Algorithm 3.5. The initial starting points $x_{i,0},\lambda_{i,0}$ and $z_{i,0}$ of both algorithms are set to be zeros.

The trajectories of selected algorithm performance indexes, including $||\mathbf{x}_{k+1}-\mathbf{x}_{k}||_{2}$, $||\varpi_{k+1}-\varpi_{k}||_{2}$, $\frac{||\mathbf{x}_{k}-\mathbf{x}^{*}||}{||\mathbf{x}^{*}||}\times 100\%$, $\frac{||\varpi_{k}-\varpi^{*}||}{||\varpi^{*}||}\times 100\%$, $||L\otimes I_{7}\bar{\lambda}_{k}||_{2}$ and $[(\mathbf{1}^{T}_{20}\otimes I_{7})\bar{\lambda}_{k}]^{T}(A\mathbf{x}_{k}-r)$, are shown in Figures 4, 5 and 6. The trajectories of the local decisions $x_{i,k}$ of some companies are shown in Figure 7, and the trajectories of $[\lambda_{i,k}]_{j}$, which stands for the $j$-th component of local Lagrangian multiplier $\lambda_{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 $\mathcal{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 $\mathcal{G}_{f}$. The methodology of this paper could be 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 & Vũ (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