Distributed generalized Nash equilibria computation of monotone games via a preconditioned proximal point algorithm

TL;DR

This paper develops distributed algorithms for computing generalized Nash equilibria in monotone games with affine coupling constraints, leveraging a preconditioned proximal point framework and allowing inexact subgame solutions.

Contribution

It introduces novel center-free distributed GNE algorithms for both equality and inequality constraints, based on a preconditioned proximal point approach, without requiring Lipschitz continuity.

Findings

Algorithms converge under monotonicity assumptions.

No need for inner-loop NE solving algorithms.

Applicable to both equality and inequality coupling constraints.

Abstract

In this paper, we investigate distributed generalized Nash equilibrium (GNE) computation of monotone games with affine coupling constraints. Each player can only utilize its local objective function, local feasible set and a local block of the coupling constraint, and can only communicate with its neighbours. We assume the game has monotone pseudo-subdifferential without Lipschitz continuity restrictions. We design novel center-free distributed GNE seeking algorithms for equality and inequality affine coupling constraints, respectively. A proximal alternating direction method of multipliers(ADMM) is proposed for the equality case, while for the inequality case, a parallel splitting type algorithm is proposed. In both algorithms, the GNE seeking task is decomposed into a sequential NE computation of regularized subgames and distributed update of multipliers and auxiliary variables, based on local data and local communication. Our two double-layer GNE algorithms need not specify the inner-loop NE seeking algorithm and moreover, only require that the strongly monotone subgames are inexactly solved. We prove their convergence by showing that the two algorithms can be seen as specific instances of preconditioned proximal point algorithms} (PPPA) for finding zeros of monotone operators. Applications and numerical simulations are given for illustration.

Equations

$$X_i({\mathbf{x}}_{-i}):=\{x_i\in {\Omega }_i|(x_i,{\mathbf{x}}_{-i})\in X\}.$$

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

$$x_i^{\ast }\in \arg \underset{x_i}{\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^e:=\prod _{i=1}^N{\Omega }_i\bigcap \{\mathbf{x}\in {\mathbf{R}}^n|\sum _{i=1}^N{A}_i{x}_i=\sum _{i=1}^N{b}_i\}.$$

$$X^i:=\prod _{i=1}^N{\Omega }_i\bigcap \{\mathbf{x}\in {\mathbf{R}}^n|\sum _{i=1}^N{A}_i{x}_i\le \sum _{i=1}^N{b}_i\}.$$

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

$$\begin{array}[]{l}\mathbf{0}\in\partial_{i}f_{i}(x_{i}^{*},\mathbf{x}^{*}_{-i})+A_{i}^{T}\lambda_{i}^{*}+N_{\Omega_{i}}(x_{i}^{*}),\;\forall i\in\mathcal{N},\\ \sum_{i=1}^{N}A_{i}x_{i}^{*}=\sum_{i=1}^{N}b_{i}.\end{array}$$

$$\begin{array}[]{l}\mathbf{0}\in\partial_{i}f_{i}(x_{i}^{*},\mathbf{x}^{*}_{-i})+A_{i}^{T}\lambda^{*}+N_{\Omega_{i}}(x_{i}^{*}),\quad\forall i\in\mathcal{N},\\ \sum_{i=1}^{N}A_{i}x_{i}^{*}=\sum_{i=1}^{N}b_{i}.\end{array}$$

$$\begin{array}[]{l}\mathbf{0}\in\partial_{i}f_{i}(x_{i}^{*},\mathbf{x}^{*}_{-i})+A_{i}^{T}\lambda^{*}+N_{\Omega_{i}}(x_{i}^{*}),\quad\forall i\in\mathcal{N},\\ \mathbf{0}\in-\sum_{i=1}^{N}(A_{i}x_{i}^{*}-b_{i})+N_{\mathbf{R}^{m}_{+}}(\lambda^{*}).\end{array}$$

$${\tilde{f}}_i(x_i,{\mathbf{x}}_{-i})=f_i(x_i,{\mathbf{x}}_{-i})+\frac{1}{2}||x_i-x_{i,k}|{|}_{R_i}^2$$

$$+[{\lambda }_{i,k}+H_i(A_i{x}_{i,k}+\sum _{l\in {\mathcal{E}}_i}{V}_{il}{z}_{l,k}-b_i){]}^T{A}_i{x}_i.$$

$$x_{i,k+1}=x_{i,k}+\rho ({\tilde{x}}_{i,k}-x_{i,k}).$$

$${\lambda }_{i,k+1}={\lambda }_{i,k}+\rho H_i(A_i{\tilde{x}}_{i,k}+\sum _{l\in {\mathcal{E}}_i}{V}_{il}{z}_{l,k}-b_i).$$

$$\begin{array}[]{l}z_{l,k+1}=z_{l,k}-\rho W_{l}(s_{j,k}-s_{i,k}).\end{array}$$

$$\mathfrak{M}^{e}:\left(\begin{array}[]{c}\mathbf{x}\\ \mathbf{Z}\\ \bar{\lambda}\\ \end{array}\right)\mapsto\left(\begin{array}[]{c}\Lambda^{T}\bar{\lambda}+(N_{\bar{\Omega}}+\partial F)\mathbf{x}\\ \bar{V}^{T}\bar{\lambda}\\ -\Lambda\mathbf{x}-\bar{V}\mathbf{Z}+\bar{b}\\ \end{array}\right)$$

$$\begin{array}[]{l}\mathbf{0}\in N_{\Omega_{i}}(\hat{x}_{i,k})+\partial_{i}f_{i}(\hat{x}_{i,k},\hat{\mathbf{x}}_{-i,k})+R_{i}(\hat{x}_{i,k}-x_{i,k})\\ \quad+A^{T}_{i}[\lambda_{i,k}+H_{i}(A_{i}x_{i,k}+\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}-b_{i})].\end{array}$$

$$\begin{array}[]{l}\mathbf{0}\in N_{\bar{\Omega}}(\hat{\mathbf{x}}_{k})+\partial F(\hat{\mathbf{x}}_{k})+R(\hat{\mathbf{x}}_{k}-\mathbf{x}_{k})\\ \qquad+\Lambda^{T}[\bar{\lambda}_{k}+H(\Lambda\mathbf{x}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})].\end{array}$$

$$\begin{array}[]{ll}\bar{\tilde{\lambda}}_{k}&=\bar{\lambda}_{k}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b}),\\ \bar{\lambda}_{k+1}&=\bar{\lambda}_{k}+\rho(\bar{\tilde{\lambda}}_{k}-\bar{\lambda}_{k}).\end{array}$$

$$\begin{array}[]{ll}\tilde{\mathbf{Z}}_{k}&=\mathbf{Z}_{k}-W\bar{V}^{T}[\tilde{\bar{\lambda}}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})],\\ \mathbf{Z}_{k+1}&=\mathbf{Z}_{k}+\rho(\tilde{\mathbf{Z}}_{k}-\mathbf{Z}_{k}).\end{array}$$

$$\begin{array}[]{l}R\mathbf{x}_{k}-\Lambda^{T}[\bar{\lambda}_{k}+H(\Lambda\mathbf{x}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})]\in(N_{\bar{\Omega}}+\partial F+R)(\hat{\mathbf{x}}_{k})\\ ||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||\leq\mu_{k}\\ \tilde{\bar{\lambda}}_{k}=\bar{\lambda}_{k}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})\\ \tilde{\mathbf{Z}}_{k}=\mathbf{Z}_{k}-W\bar{V}^{T}(\tilde{\bar{\lambda}}_{k}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})\\ \mathbf{x}_{k+1}=\mathbf{x}_{k}+\rho(\tilde{\mathbf{x}}_{k}-\mathbf{x}_{k}),\;\bar{\lambda}_{k+1}=\bar{\lambda}_{k}+\rho(\tilde{\bar{\lambda}}_{k}-\bar{\lambda}_{k})\\ \mathbf{Z}_{k+1}=\mathbf{Z}_{k}+\rho(\tilde{\mathbf{Z}}_{k}-\mathbf{Z}_{k})\end{array}$$

$$-{\Lambda }^T[{\bar{\lambda }}^{\ast }+H(\Lambda {\mathbf{x}}^{\ast }+\bar{V}{\mathbf{Z}}^{\ast }-\bar{b})]\in (N_{\bar{\Omega }}+\partial F)({\mathbf{x}}^{\ast })$$

$${\bar{\lambda }}^{\ast }={\bar{\lambda }}^{\ast }+H(\Lambda {\mathbf{x}}^{\ast }+\bar{V}{\mathbf{Z}}^{\ast }-\bar{b})$$

$${\mathbf{Z}}^{\ast }={\mathbf{Z}}^{\ast }-W{\bar{V}}^T({\bar{\lambda }}^{\ast }+H(\Lambda {\mathbf{x}}^{\ast }+\bar{V}{\mathbf{Z}}^{\ast }-\bar{b})$$

$${\tilde{f}}_i(x_i,{\mathbf{x}}_{-i})=f_i(x_i,{\mathbf{x}}_{-i})+\frac{1}{2}||x_i-x_{i,k}|{|}_{R_i}^2+{\lambda }_{i,k}^T{A}_i{x}_i,$$

$$x_{i,k+1}=x_{i,k}+\rho ({\tilde{x}}_{i,k}-x_{i,k}).$$

$$\begin{array}[]{l}\tilde{z}_{l,k}=z_{l,k}-W_{l}(\lambda_{j,k}-\lambda_{i,k}),\\ z_{l,k+1}=z_{l,k}+\rho(\tilde{z}_{l,k}-z_{l,k}).\end{array}$$

$$\begin{array}[]{l}\tilde{{\lambda}}_{i,k}=P^{H^{-1}_{i}}_{\mathbf{R}^{m}_{+}}[\lambda_{i,k}+H_{i}(A_{i}(2\tilde{x}_{i,k}-x_{i,k})\\ \qquad+\sum_{l\in\mathcal{E}_{i}}V_{il}(2\tilde{z}_{l,k}-z_{l,k})-b_{i})],\\ \lambda_{i,k+1}=\lambda_{i,k}+\rho(\tilde{{\lambda}}_{i,k}-\lambda_{i,k}).\end{array}$$

$${\mathfrak{M}}^{i}:\left(\begin{array}[]{c}\mathbf{x}\\ \mathbf{Z}\\ \bar{\lambda}\\ \end{array}\right)\mapsto\left(\begin{array}[]{c}\Lambda^{T}\bar{\lambda}+(N_{\bar{\Omega}}+\partial F)\mathbf{x}\\ \bar{V}^{T}\bar{\lambda}\\ -\Lambda\mathbf{x}-\bar{V}\mathbf{Z}+\bar{b}+N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda})\\ \end{array}\right)$$

$$\begin{array}[]{l}\mathbf{0}\in(N_{\bar{\Omega}}+\partial F)(\hat{\mathbf{x}}_{k})+R(\mathbf{x}_{k}-\hat{\mathbf{x}}_{k})+\Lambda^{T}\bar{\lambda}_{k},\\ ||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||\leq\mu_{k},\quad\mathbf{x}_{k+1}=\mathbf{x}_{k}+\rho(\tilde{\mathbf{x}}_{k}-\mathbf{x}_{k})\\ \tilde{\mathbf{Z}}_{k}=\mathbf{Z}_{k}-W\bar{V}^{T}{\bar{\lambda}}_{k},\quad\mathbf{Z}_{k+1}=\mathbf{Z}_{k}+\rho(\tilde{\mathbf{Z}}_{k}-\mathbf{Z}_{k}),\\ \mathbf{0}\in N_{\mathbf{R}^{mN}_{+}}(\tilde{\bar{\lambda}}_{k})+H^{-1}\big{\{}\tilde{\bar{\lambda}}_{k}-\bar{\lambda}_{k}\\ -H[\Lambda(2\tilde{\mathbf{x}}_{k}-\mathbf{x})+\bar{V}(2\tilde{\mathbf{Z}}_{k}-\mathbf{Z}_{k})-\bar{b})]\big{\}}\\ \bar{\lambda}_{k+1}=\bar{\lambda}_{k}+\rho(\tilde{\bar{\lambda}}_{k}-\bar{\lambda}_{k})\end{array}$$

$$\begin{array}[]{l}\Phi(\varpi_{k}-\hat{\varpi}_{k})\in{\mathfrak{M}}\hat{\varpi}_{k},\;\;||\hat{\varpi}_{k}-\tilde{\varpi}_{k}||\leq\nu_{k},\\ \varpi_{k+1}=\varpi_{k}+\rho(\tilde{\varpi}_{k}-\varpi_{k}),\end{array}$$

$${\varpi }_{k+1}=R_{\mathfrak{M}}{\varpi }_k=(\mathrm{Id}+\mathfrak{M}{)}^{-1}{\varpi }_k.$$

$$\left(\begin{array}[]{cccc}\mathbf{0}&\Lambda^{T}&\mathbf{0}&-\Lambda^{T}\\ -\Lambda&\mathbf{0}&-\bar{V}&\mathbf{0}\\ \mathbf{0}&\bar{V}^{T}&\mathbf{0}&-\bar{V}^{T}\\ \Lambda&\mathbf{0}&\bar{V}&\mathbf{0}\\ \end{array}\right)\varpi+\left(\begin{array}[]{c}(N_{\bar{\Omega}}+\partial F)\mathbf{x}\\ \bar{b}\\ \mathbf{0}\\ -\bar{b}\\ \end{array}\right).$$

Full text

I Introductions

Generalized Nash equilibrium and its distributed computation is an important research topic in decision making problems over large-scale multi-agent networks. Examples include power allocation over cognitive radio networks, [1, 3, 2], demand response and electric vehicle charging management in smart grids, [4, 7, 5, 6], rate control over optical networks, [8, 9], and opinion evolution over social networks, [10, 11]. Each agent (player) controls its decision, and has an objective function to be optimized, which depends on other players’ decisions. Moreover, each player’s feasible set can depend on other players’ decisions through coupling constraints, such as when they share limited network resources. Generalized Nash Equilibrium (GNE), firstly proposed in [12], is a reasonable solution, since at a GNE no player can decrease/increase its cost/utility by unilaterally changing its local decision to another feasible one. Interested readers can refer to [13] for a review on GNE.

Distributed GNE computation methods are quite appealing for noncooperative games over large-scale networks, in which the local data of each player, including own objective function and own feasible set, are kept by each player. Moreover, when the coupling constraint is a sum of separable local functions, it is also appealing to have each player only knowing its local constraint function, i.e., local contribution to the coupling constraint. Since local data is not required to be transmitted to a central node, the communication burden could be relieved, and the privacy of each player gets protected. Recently, distributed NE/GNE computation methods have received increasing research attention, see [2, 7, 5, 3, 4, 6] and [14, 24, 23, 22, 20, 21, 16, 19, 15, 18, 17]. Different information structures are considered, depending on whether or not there exists a coordination center. For example, the methods in [3][5][6] all utilize a central node to update and broadcast certain coordination/incentive signals based on all players’ decisions. Notice that [5] considers aggregative games where the agents are coupled through aggregative variables, hence, it is efficient to adopt a coordination center if permitted. Meanwhile, totally center-free distributed GNE computation algorithms have been proposed in [7, 22, 24, 23] assuming that each player is able to observe the decisions on which its local objective function or constraint function explicitly depends on. On the other hand, in the distributed NE computation algorithms of [19, 20, 21], each player is only required to have local communications with its neighbours, and each player computes an estimation of other players’ decisions or aggregative variables by resorting to consensus dynamics.

Typically, the objective function of each player is convex only with respect to its own decision. Then an NE/GNE can be computed by solving a (generalized) Variational Inequality (VI) problem constructed with the game’s pseudo-gradient/subdifferential (PG/PS) [1, 2, 3, 7, 13]. Various monotonicity and Lipschitz continuity assumptions on PG/PS play a fundamental role in the design and analysis of distributed NE/GNE seeking algorithms. [3] assumes a strongly monotone PG to get the cocoercivity of the dual operator, and show the convergence of double-layer dual gradient GNE seeking methods. [22] and [24] combine strong monotonicity and Lipschitz continuity to ensure the cocoercivity of PG, and propose primal/primal-dual gradient methods for distributed GNE computation. [6, 11] and [18] consider aggregative games with quadratic objective functions, hence also adopt a strong monotone and Lipschitz PG. [20] and [23] consider games with strictly monotone and Lipschitz PG. [20] proposes a “gradient”+“consensus” algorithm for distributed NE seeking, while [23] utilizes a continuous-time gradient flow algorithm to seek a GNE of aggregative games.

For NE seeking with only monotone PGs, [2] proposes a double-layer proximal best-response algorithm that involves solving regularized subgames at each iteration, while [16] proposes a single time-scale/layer regularized (sub)gradient algorithm with diminishing step-sizes. For GNE seeking of monotone games, [9] proposes a double-layer dual extragradient method and [15] adopts the single-layer Tikhonov regularization algorithm with diminishing step-sizes, both assuming Lipschitz continuity and using a central coordinator. [7] proposes a primal-dual gradient algorithm, and [25] proposes a payoff-based algorithm for GNE seeking with pseudo-monotone PGs, both with diminishing step-sizes.

Motivated by the above, we investigate center-free distributed algorithms for computing GNE of monotone games with affine coupling constraints. The players’ decisions are coupled together with a globally shared affine constraint, while each player only knows a local block of the constraint. We consider both equality and inequality constraints which cover many task/resource allocation games, [3, 6, 23, 24]. Compared with previous works, the key difference is that we only assume a monotone pseudo-subdifferential without Lipschitz continuity restrictions. We propose center-free GNE algorithms with fixed step-sizes where each player only utilizes its local data and has a peer-to-peer communication with its neighbours. To the best of our knowledge, this distributed GNE computation has not been discussed in literature under this general form.

We adopt the variational GNE as a refined solution and use primal-dual analysis to reformulate GNE seeking as the problem of finding zeros of monotone operators for equality and inequality cases, respectively. The monotone operators are composed of a skew-symmetric linear operator (with both the constraint matrices and a matrix related to communication graph) and an operator involving PS. In general, the proximal point algorithm can be applied for solving monotone inclusion problems without Lipschitz restriction. However, it is not directly applicable to our GNE problem because it requires to compute the inverse of a graph-related skew-symmetric matrix, which is prohibitive in distributed algorithms. To overcome these challenges, we propose novel distributed GNE seeking algorithms based on Preconditioned Proximal Point Algorithm (PPPA), for equality and inequality cases, respectively. For the equality case, we call it proximal alternating direction method of multipliers (ADMM), partially motivated by [26]. For the inequality case, we call it proximal parallel splitting algorithm, partially motivated by [27]. Both algorithms use appropriately chosen operators and preconditioning matrices, which ensure that the resolvent evaluation of monotone operators is realizable by local computation and communication. The proposed algorithms decompose the GNE computation into sequential NE computation for regularized subgames and distributed update of local multipliers and auxiliary variables. Hence, our algorithms are double-layer algorithms, similar to [2, 3, 9], but the inner-loop NE seeking algorithms need not be specified while the subgame only needs to be solved inexactly. By using proximal terms, the subgame is regularized to have strongly monotone PS, hence it can be efficiently solved by existing NE seeking distributed algorithms, such as the best-response algorithm in [2]. The inexactness in solving the subgames is also considered and relaxation steps are applied to all variables, which potentially could improve convergence speed. Moreover, proximal ADMM enjoys the feature of utilizing the most recent available information whenever possible. In both cases, the algorithms’ convergence is proved for fixed step-sizes by relating them to PPPA, and showing that they can be seen as specific instances of PPPA, while PPPA’s convergence can be shown based on averaged operator theory.

To summarize, the main contributions of this work are as follows. (i): The game model only assumes a monotone PS without Lipschitz continuity, hence it is a generalization of previous ones. Both equality and inequality affine coupling constraints are considered. (ii): Novel center-free GNE seeking algorithms with peer-to-peer communication are introduced. Since only monotonicity is imposed, the double-layer algorithms could be implemented after the NE algorithm is chosen tailored to the specific practical problem. Moreover, thanks to the proximal terms, the subgames are regularized to have strongly monotone PS/PGs, hence, could be efficiently solved. (iii): The algorithms are related to PPPA for monotone inclusion, revealing the algorithms’ intrinsic structure. Their convergence is proved for fixed step-sizes.

The paper is organized as follows. Section II gives the preliminary background. Section III formulates the noncooperative game and basic assumptions. Section IV gives distributed GNE computation algorithms for both equality and inequality constraint cases, and analyzes their limiting points. Section V presents the algorithms’ convergence analysis. Section VI gives application examples and simulation studies. Section VII draws the concluding remarks.

II notations and preliminaries

In this section, we review the notations and preliminary notions in monotone and averaged operators from [28].

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||=\sqrt{x^{T}x}$ denotes the induced norm. $||x||^{2}_{G}$ denotes $\langle x,Gx\rangle$ for a symmetric matrix $G$. Denote $\mathbf{1}_{m}=(1,...,1)^{T}\in\mathbf{R}^{m}$ and $\mathbf{0}_{m}=(0,...,0)^{T}\in\mathbf{R}^{m}$. $diag\{A_{1},...,A_{N}\}$ represents the block diagonal matrix with $A_{1},...,A_{N}$ on its main diagonal. Denote $col(x_{1},....,x_{N})$ as the stacked column vector of $x_{1}$ to $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$. Denote $int(\Omega)$ as the interior of $\Omega$ and $ri(\Omega)$ as the relative interior of $\Omega$. Denote $\times_{i=1,...,N}\Omega_{i}$ or $\prod_{i=1}^{N}\Omega_{i}$ as the Cartesian product of $\Omega_{i},i=1,...,N$.

Let $\mathfrak{A}:\mathbf{R}^{m}\rightarrow 2^{\mathbf{R}^{m}}$ be a set-valued operator. ${\rm Id}$ denotes 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,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\}$. The inverse of $\mathfrak{A}$ is defined via $gra\mathfrak{A}^{-1}=\{(u,x)|(x,u)\in gra\mathfrak{A}\}$. The zero set of $\mathfrak{A}$ is $zer\mathfrak{A}=\{x\in\mathbf{R}^{m}|\mathbf{0}\in\mathfrak{A}x\}$. The sum of $\mathfrak{A}$ and $\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}\}$. Define the resolvent of $\mathfrak{A}$ as $R_{\mathfrak{A}}=({\rm Id}+\mathfrak{A})^{-1}$.

Operator $\mathfrak{A}$ is monotone if $\forall(x,u),\forall(y,v)\in gra\mathfrak{A}$, we have $\langle x-y,u-v\rangle\geq 0.$ $\mathfrak{A}$ is maximally monotone if $gra\mathfrak{A}$ is not strictly contained in the graph of any other monotone operator. A skew-symmetric matrix $A=-A^{T}$ defines a maximally monotone operator $Ax$ ([28], p. $298$). Suppose $\mathfrak{A}$ and $\mathfrak{B}$ are maximally monotone operators and $0\in int(dom\mathfrak{A}-dom\mathfrak{B})$, then $\mathfrak{A}+\mathfrak{B}$ is also maximally monotone. For a proper lower semi-continuous convex (l.s.c.) function $f$, its subdifferential operator $\partial f:domf\rightarrow 2^{\mathbf{R}^{m}}$ is $\partial f:\;x\mapsto\{g|f(y)\geq f(x)+\langle g,y-x\rangle,\forall y\in domf\}.$ $\partial f$ is maximally monotone and $Prox_{f}=R_{\partial f}:\mathbf{R}^{m}\rightarrow domf$ is called the proximal operator of $f$, i.e., $Prox_{f}:x\mapsto\arg\min_{u\in domf}f(u)+\frac{1}{2}||u-x||_{2}^{2}.$

Define the indicator function of $\Omega$ as $\iota_{\Omega}(x)=0$ if $x\in\Omega$ and $\iota_{\Omega}(x)=\infty$ if $x\notin\Omega.$ For a closed convex set $\Omega$, $\iota_{\Omega}$ is a proper l.s.c. function. $\partial\iota_{\Omega}$ is also the normal cone operator of $\Omega$, i.e., $N_{\Omega}(x)$, where $N_{\Omega}(x)=\{v|\langle v,y-x\rangle\leq 0,\forall y\in\Omega\}$ and $domN_{\Omega}=\Omega$. Given a symmetric positive definite matrix $G$, define $P^{G}_{\Omega}(x)=\arg\min_{y}(\iota_{\Omega}(y)+\frac{1}{2}||x-y||_{G}^{2})$.

For a single-valued operator $T:\Omega\subset\mathbf{R}^{m}\rightarrow\mathbf{R}^{m}$, $x\in\Omega$ is a fixed point of $T$ if $Tx=x$. $T$ is nonexpansive if it is $1-$Lipschitzian, i.e., $||T(x)-T(y)||\leq||x-y||,\forall x,y\in\Omega$. $T$ is contractive if $\exists\gamma\in(0,1)$ s.t. $||T(x)-T(y)||\leq\gamma||x-y||,\forall x,y\in\Omega.$ Let $\alpha\in(0,1)$, then $T$ is $\alpha-$averaged, denoted as $T\in\mathcal{A}(\alpha)$, if $\exists$ a nonexpansive operator $T^{{}^{\prime}}$ such that $T=(1-\alpha){\rm Id}+\alpha T^{{}^{\prime}}$. If $T\in\mathcal{A}(\frac{1}{2})$, $T$ is called firmly nonexpansive.

III Game formulation

Consider a set of players (agents) $\mathcal{N}=\{1,\cdots,N\}$ that are involved in the following noncooperative game with shared coupling constraints. Player $i\in\mathcal{N}$ controls its own decision (strategy or action) $x_{i}\in\Omega_{i}\subset\mathbf{R}^{n_{i}}$, where $\Omega_{i}$ is its private feasible set. Let $\mathbf{x}=col(x_{1},\cdots,x_{N})\in\mathbf{R}^{n}$ denote the decision profile, i.e., the stacked vector of all agents’ decisions, with $\sum_{i=1}^{N}n_{i}=n$. Let $\mathbf{x}_{-i}=col(x_{1},\cdots,x_{i-1},x_{i+1},\cdots,x_{N})$ denote the decision profile of all agents except player $i$. Player $i$ aims to optimize its own objective function within its feasible set, $f_{i}(x_{i},\mathbf{x}_{-i}):\bar{\Omega}\rightarrow\mathbf{R}$ where $\bar{\Omega}=\prod_{i=1}^{N}\Omega_{i}\subset\mathbf{R}^{n}$. Note that $f_{i}(x_{i},\mathbf{x}_{-i})$ is coupled with other players’ decisions $\mathbf{x}_{-i}$. Moreover, all the players’ decisions are coupled together through a *globally shared set * $X\subset\mathbf{R}^{n}$. Hence, player $i$ has a set-valued map $X_{i}(\mathbf{x}_{-i}):\mathbf{R}^{n-n_{i}}\rightarrow 2^{\mathbf{R}^{n_{i}}}$ that specifies its feasible set defined as

[TABLE]

Given $\mathbf{x}_{-i}$, player $i$’s best-response strategy is

[TABLE]

A generalized Nash equilibrium (GNE) $\mathbf{x}^{*}=col(x_{1}^{*},\cdots,x_{N}^{*})$ is defined at the intersection of all players’ best-response sets,

[TABLE]

We consider the set $X$ defined via two types of shared affine coupling constraints, equality and inequality constraints. For the equality constraint case, $X=X^{e}$ where we denote

[TABLE]

For the inequality constraint case, $X=X^{i}$ where

[TABLE]

In both (3) and (4), $A_{i}\in\mathbf{R}^{m\times n_{i}}$ and $b_{i}\in\mathbf{R}^{m}$ as well as $\Omega_{i}$ are private data of player $i$. Thereby, the shared set $X$ couples all players’ feasible sets, but is not known by any agent. We consider the following assumption on the game in (1).

Assumption 1

For player $i$, $f_{i}(x_{i},\mathbf{x}_{-i})$ is a proper l.s.c. function with respect to $x_{i}$ given any fixed $\mathbf{x}_{-i}$, and its subdifferential with respect to $x_{i}$ is $\partial_{i}f_{i}(x_{i},\mathbf{x}_{-i})$. The pseudo-subdifferential of the game in (1) defined as $\partial F(\mathbf{x}):\mathbf{x}\rightarrow\prod_{i=1}^{N}\partial_{i}f_{i}(x_{i},\mathbf{x}_{-i})$ is maximally monotone. $\Omega_{i}$ is a closed convex set with nonempty interior. $X^{e}$ in (3) has nonempty relative interiors, and $X^{i}$ in (4) has nonempty interiors. $X_{i}(\mathbf{x}_{-i})$ has nonempty relative interiors for $\mathbf{x}_{-i}\in\prod_{j=1,j\neq i}^{N}\Omega_{j}$ when $X=X^{e}$, and $X_{i}(\mathbf{x}_{-i})$ has nonempty interiors for $\mathbf{x}_{-i}\in\prod_{j=1,j\neq i}^{N}\Omega_{j}$ when $X=X^{i}$.

Remark 1

In many practical cases, $f_{i}(x_{i},\mathbf{x}_{-i})$ has a splitting structure such as $f_{i}(x_{i},\mathbf{x}_{-i})=g_{i}(x_{i},\mathbf{x}_{-i})+l_{i}(x_{i})$, [15], where $g_{i}(x_{i},\mathbf{x}_{-i})$ is differentiable and convex with respect to $x_{i}$, and $l_{i}(x_{i})$ is a local l.s.c. regularization/cost term. Denote $\nabla_{p}G(\mathbf{x})=col(\nabla_{1}g_{1}(x_{1},\mathbf{x}_{-1}),\cdots,\nabla_{N}g_{N}(x_{N},\mathbf{x}_{-N}))$ where $\nabla_{i}g_{i}(x_{i},\mathbf{x}_{-i})$ is the gradient of $g_{i}$ with respect to $x_{i}$ and $\partial L(\mathbf{x}):\mathbf{x}\rightarrow\prod_{i=1}^{N}\partial l_{i}(x_{i})$. Then $\partial L(\mathbf{x})$ is maximally monotone, since it is the subdifferential of $\sum_{i=1}^{N}l_{i}(x_{i})$. In this case, $\partial F(\mathbf{x})=\partial L(\mathbf{x})+\nabla_{p}G(\mathbf{x})$ is maximally monotone when $\nabla_{p}G(\mathbf{x})$ is monotone.

Define the generalized variational inequality (GVI) problem

[TABLE]

According to Proposition 12.4 in [1], any solution of (5) is a GNE of game in (1), called variational GNE.

Let us first analyze the equality constraint case, $X=X^{e}$. Under Assumption 1, $\mathbf{x}^{*}$ is a GNE of the game in (1) if and only if $\forall i\in\mathcal{N}$ there exists $\lambda_{i}^{*}\in\mathbf{R}^{m}$ such that,

[TABLE]

Meanwhile, based on the Lagrangian duality for GVI (Equation (12.4) of [1]), $\mathbf{x}^{*}$ is a solution of GVI in (5) with $X=X^{e}$ if and only if there exists a multiplier $\lambda^{*}\in\mathbf{R}^{m}$ such that

[TABLE]

By comparing the KKT conditions in (6) and (7), we have that any solution to GVI in (5) with $X=X^{e}$ is a GNE of the game in (1) with all players having the same local multiplier.

Similarly, for the inequality case $X=X^{i}$, $\mathbf{x}^{*}$ is a solution of GVI (5) with $X=X^{i}$ if and only if there exists a multiplier $\lambda^{*}\in\mathbf{R}_{+}^{m}$ such that

[TABLE]

Not every GNE of the considered game in (1) is a solution to the GVI in (5). Since the variational GNE has an economic interpretation of no price discrimination and enjoys a stability and sensitivity property (refer to [1]), we aim to propose novel distributed algorithms for computing a variational GNE of the monotone game for $X=X^{e}$ and $X=X^{i}$, respectively.

Assumption 2

The solution set of GVI in (5) is nonempty for both $X=X^{e}$ and $X=X^{i}$, or equivalently, the considered game in (1) has at least a variational GNE.

Remark 2

Some sufficient conditions for the existence of solutions to monotone GVI can be found in [1] and [15]. For example, compactness of $\Omega_{i},\forall i\in\mathcal{N}$ ensures Assumption 2.

IV Distributed GNE computation algorithms

In this section, we propose distributed algorithms that players can use to find a solution of GVI (5) for $X=X^{e}$ and $X=X^{i}$, respectively. We focus on distributed variational GNE computation because of two reasons. Firstly, player $i$ can only manipulate its local $f_{i}(x_{i},\mathbf{x}_{-i})$, $A_{i}$, $b_{i}$ and $\Omega_{i}$ for local computation, since these contain its private information. Secondly, we assume there is no central node that has bidirectional communications with all players, either because this could be inefficient from a communication point of view, or because it might be not possible to have such a central node. Thus, each player only uses its local data for local computation, and has peer-to-peer communication with its neighbours for local coordination.

We first introduce the communication graph and algorithm notations in IV-A. We give the proximal ADMM for equality constraint case in IV-B, and the distributed algorithm for inequality constraint case in IV-C.

IV-A Communication graph and algorithm variables

To facilitate the distributed coordination, players are able to communicate with their neighbours through a connected and undirected graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$. The edge set is $\mathcal{E}\subset\mathcal{N}\times\mathcal{N}$, $(i,j)\in\mathcal{E}$ if agent $i$ and agent $j$ can share information with each other, and agents $j$, $i$ are called neighbours. A path of graph $\mathcal{G}$ is a sequence of distinct agents in $\mathcal{N}$ such that any consecutive agents in the sequence are neighbours. Agent $j$ is said to be connected to agent $i$ if there is a path from $j$ to $i$. $\mathcal{G}$ is connected if any two agents are connected.

Obviously, $|\mathcal{N}|=N$, and we denote $|\mathcal{E}|=M$. The edges are labeled with $e_{l}$, $l=1,\cdots,M$. Without loss of generality, $e_{l}=(i,j)$ is arbitrarily ordered and denoted by $i\rightarrow j$. Define $\mathcal{E}^{in}_{i}$ and $\mathcal{E}^{out}_{i}$ for agent $i$ as follows: $e_{l}\in\mathcal{E}^{in}_{i}$ if agent $i$ is the targeted point of $e_{l}$; $e_{l}\in\mathcal{E}^{out}_{i}$ if agent $i$ is the starting point of $e_{l}$. Then denote $\mathcal{E}_{i}=\mathcal{E}^{in}_{i}\bigcup\mathcal{E}^{out}_{i}$ as the set of edges adjoint to agent $i$. Define the incidence matrix of $\mathcal{G}$ as ${V}\in\mathbf{R}^{N\times M}$ with ${V}_{il}=1$ if $e_{l}\in\mathcal{E}_{i}^{in}$, and ${V}_{il}=-1$ if $e_{l}\in\mathcal{E}_{i}^{out}$, otherwise ${V}_{il}=0$. We have $\mathbf{1}^{T}_{N}V=\mathbf{0}^{T}_{M}$, and $V^{T}x=\mathbf{0}_{M}$ if and only if $x\in\{\alpha\mathbf{1}_{N}|\alpha\in\mathbf{R}\}$ when $\mathcal{G}$ is connected. Denote $\mathcal{N}_{l}=\{i,j\}$ as the pair of agents connected by edge $e_{l}=(i,j)$.

We introduce the variables. Firstly, each player has a local decision $x_{i}\in\Omega_{i}$ and a local multiplier $\lambda_{i}\in\mathbf{R}^{m}$. According to KKT (7) and (8), in steady-state all players should have the same local multiplier, i.e., $\lambda_{i}=\lambda^{*},\forall i\in\mathcal{N}$. To facilitate the coordination for the consensus of local multipliers and to ensure the coupling constraint, we consider an auxiliary variable $z_{l}\in\mathbf{R}^{m}$ associated with edge $e_{l}$ of graph $\mathcal{G}$. Notice that $\mathcal{G}$ is undirected and the edges are arbitrarily ordered, therefore, we can have any agent from $\mathcal{N}_{l}$ to maintain $z_{l}$. For clarity, we let the starting agent of an edge to maintain the corresponding edge variable. That is agent $i$ will take the responsibility for maintaining $z_{l}$ if $e_{l}\in\mathcal{E}^{out}_{i}$.

Before presenting the algorithms we first make some observations. The algorithms are based on decomposing the GNE computation into sequential NE computation for regularized subgames and distributed update of local multipliers and auxiliary variables. The regularized subgames are made to have strongly monotone PS with the help of proximal terms, hence can be efficiently solved by existing distributed algorithms, such as the best-response algorithm in [2]. The update of the local multipliers has to be done so that in steady-state they are the same, and satisfy the optimality conditions (7), (8) involving the constraints, while using only local information. Towards this we use the auxiliary variables $z_{l}$, which have a double role: to help in estimating the contribution of the other players’ in the constraints and to enforce consensus.

Let $x_{i,k},\lambda_{i,k}$ and $z_{l,k}$ denote $x_{i},\lambda_{i}$ and $z_{l}$ at iteration $k$.

IV-B Proximal ADMM for $X=X^{e}$

The distributed algorithm for computing a variational GNE of game in (1) when $X=X^{e}$ is given as follows.

Algorithm 1

Step 1–update of $x_{i,k}$:*

• •

Player $i$ receives $z_{l,k},l\in\mathcal{E}^{in}_{i}$ through $\mathcal{G}$.

• •

Construct a subgame where player $i$ has a decision $x_{i}\in\Omega_{i}$ and an objective function $\tilde{f}_{i}(x_{i},\mathbf{x}_{-i})$,

[TABLE]

and denote its NE by $\hat{\mathbf{x}}_{k}=col(\hat{x}_{1,k},\cdots,\hat{x}_{N,k})$.

• •

Players compute $\tilde{\mathbf{x}}_{k}=col(\tilde{x}_{1,k},\cdots,\tilde{x}_{N,k})$ as an inexact solution to subgame (• ‣ 1) such that $||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||\leq\mu_{k}$, where $\mu_{k}$ is described below.

• •

Player $i$ updates its local decision $x_{i,k}$ with

[TABLE]

Step 2–update of $\lambda_{i,k}$:

[TABLE]

Step 3–update of $z_{l,k}$:

Let $s_{i,k}=\frac{1}{\rho}\lambda_{i,k+1}+\frac{\rho-1}{\rho}\lambda_{i,k}+H_{i}(A_{i}\tilde{x}_{i,k}+\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}-b_{i}).$ For $e_{l}\in\mathcal{E}_{i}^{out}$, player $i$ receives $s_{j,k}$, $j\in\mathcal{N}_{l}\setminus\{i\}$, and updates $z_{l,k}$ with

[TABLE]

$\{\mu_{k}\}$ is a nonnegative sequence s.t. $\sum_{k=1}^{\infty}\mu_{k}<\infty$, $\rho\in[1,2)$ is a fixed relaxation/extrapolation step-size, and $R_{i}\in\mathbf{R}^{n_{i}\times n_{i}}$, $H_{i}\in\mathbf{R}^{m\times m}$ and $W_{l}\in\mathbf{R}^{m\times m}$, $l\in\mathcal{E}_{i}^{out}$ are local parameters (step-sizes) that are symmetric positive definite matrices.*

We give next some intuition behind Algorithm 1’s design. Since $A_{i}$, $b_{i}$ are private data, the coupling constraint $X=X^{e}$ is not completely known by any player. Note that in steady-state we should have $A_{i}x_{i}^{*}-b_{i}=\sum_{j=1,j\neq i}^{N}(A_{j}x^{*}_{j}-b_{j})$ due to (7), where the right-hand side is unknown information for player $i$. The penalized cost $\tilde{f}_{i}$ in (• ‣ 1) is composed of a proximal term $f_{i}(x_{i},x_{-i})+\frac{1}{2}||x_{i}-x_{i,k}||_{R_{i}}^{2}$ (to regularize the subgames), a Lagrangian term $\langle\lambda_{i},A_{i}x_{i}\rangle$ and a penalty term. The penalty term is based on linearizing the quadratic penalty $\frac{1}{2}||A_{i}x_{i}-b_{i}+\sum_{j=1,j\neq i}^{N}(A_{j}x_{j,k}-b_{j})||_{H_{i}}^{2}$ at $x_{i,k}$, which should be zero in steady-state cf. (7). This gives $\langle H_{i}(A_{i}x_{i,k}+\sum_{j=1,j\neq i}^{N}(A_{j}x_{j,k}-b_{j})-b_{i}),A_{i}x_{i}\rangle$ after dropping all constants. To overcome the need for information about the other players $j\neq i$, in (• ‣ 1) this term is estimated as $(H_{i}[A_{i}x_{i,k}+\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}-b_{i}])^{T}A_{i}x_{i}$, via the auxiliary variables $z_{l}$. A similar term is used in the local multiplier $\lambda_{i}$’s update, (11). Player $i$ uses $A_{i}\tilde{x}_{i,k}+\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}-b_{i}$ as an estimation of $\sum_{j=1}^{N}(A_{j}x_{j,k}-b_{j})$ to update its $\lambda_{i}$. The update for $z_{l}$, (12), has an integrator dynamics form driven by the difference between $\lambda_{i}$ and $\lambda_{j}$, since $e_{l}=(i,j)$, and ensures the consensus of local multipliers. Meanwhile, (11) also utilizes $z_{l,k}$, that is the integrator for differences between multipliers, as the feedback signal to reach consensus of local multipliers.

We show in Theorem 1 that at the limit point of the algorithm $A_{i}x_{i}^{*}-b_{i}=\sum_{l\in\mathcal{E}_{i}}V_{il}z^{*}_{l},\forall i\in\mathcal{N}$, while $A_{i}x_{i}^{*}-b_{i}=\sum_{j=1,j\neq i}^{N}(A_{j}x^{*}_{j}-b_{j})$ due to (7). Hence, $\sum_{l\in\mathcal{E}_{i}}V_{il}z^{*}_{l}$, generated as an output of (12) is an estimation of $\sum_{j=1,j\neq i}^{N}(A_{j}x^{*}_{j}-b_{j})$, and $\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}$ as used by player $i$ is a dynamical estimator for $\sum_{j=1,j\neq i}^{N}(A_{j}x_{j,k}-b_{j})$. Motivated by [29], the auxiliary variable $z_{l}$ has an interpretation of network flow. In fact, if we regard $A_{i}x_{i}$ as in-flow at node $i$ and $b_{i}$ as out-flow at node $i$, and $\sum_{i=1}^{N}A_{i}x_{i}=\sum_{i=1}^{N}b_{i}$ is a conservative network flow balancing constraint. Thereby, $z_{l}$ can be regarded as flow on each edge to ensure the balancing constraint. All in all, variables $z_{l}$ estimate the other players’ contribution to coupling constraints, and ensure local multipliers reach consensus.

Algorithm 1 updates each coordinate with the most recent information in a Gauss-Seidel manner and uses proximal terms, hence is called proximal ADMM. It uses relaxation steps, $\rho\in[1,2)$, to perform extrapolations of all variables, which in practice could accelerate convergence (refer to Figure 2 of [30] and numerical studies in [26, 27]). It is a center-free distributed algorithm with peer-to-peer communications. In Step 1, player $i$ communicates with its neighbours to get $z_{l,k},l\in\mathcal{E}_{i}^{in}$. The NE of subgames can be computed in a distributed manner with existing algorithms such as best-response algorithms in [2] and gradient algorithms in [16], [20] and [21], which only involve local computations and communications. In Step 2, player $i$ uses its $\tilde{x}_{i,k}$ and locally available $z_{l,k},l\in\mathcal{E}_{i}$ to update its local multiplier $\lambda_{i}$. In Step 3, player $i$ computes $s_{i,k}$ with its local information, and receives $s_{j,k},j\in\mathcal{N}_{l}\setminus\{i\}$ to update $z_{l,k}$.

Next, we put Algorithm 1 in a compact form and show that its limiting point $\mathbf{x}^{*}$ is a variational GNE of game in (1) when $X=X^{e}$. We use the following compact notations. Denote $\bar{\lambda}=col(\lambda_{1},\cdots,\lambda_{N})$ and $\mathbf{Z}=col(z_{1},\cdots,z_{M})$. Denote $R=diag\{R_{1},\cdots,R_{N}\}$, $W=diag\{W_{1},\cdots,W_{M}\}$, $H=diag\{H_{1},\cdots,H_{N}\}$, $\bar{V}=V\otimes I_{m}$, $\Lambda=diag\{A_{1},\cdots,A_{N}\}$, and $\bar{b}=col(b_{1},\cdots,b_{N})$.

Theorem 1

Suppose that Assumption 1 and 2 hold for game (1) when $X=X^{e}$. Then any limiting point $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ of Algorithm 1 belongs to the zeros of operator $\mathfrak{M}^{e}$ defined by

[TABLE]

Meanwhile, any zero $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ of $\mathfrak{M}^{e}$ (13) has the $\mathbf{x}^{*}$ component as a variational GNE of game (1) when $X=X^{e}$.

Proof: We write Algorithm 1 in a compact form. Due to proximal terms $\frac{1}{2}||x_{i}-x_{i,k}||_{R_{i}}^{2}$ and Assumption 1, the subgame in Step 1 has a strongly monotone pseudo-subdifferential, hence its NE $\hat{\mathbf{x}}_{k}$ exists and is also unique. Therefore, $\hat{x}_{i,k}=\arg\min_{x_{i}\in\Omega_{i}}\tilde{f}_{i}(x_{i},\hat{\mathbf{x}}_{-i,k})$, and its KKT condition is

[TABLE]

Concatenating all KKT conditions together and using the compact notations defined before, yields for $\hat{\mathbf{x}}_{k}$

[TABLE]

We also have $||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||\leq\mu_{k}$ and $\mathbf{x}_{k+1}=\mathbf{x}_{k}+\rho(\tilde{\mathbf{x}}_{k}-\mathbf{x}_{k})$.

Let $\tilde{\lambda}_{i,k}=\lambda_{i,k}+H_{i}(A_{i}\tilde{x}_{i,k}+\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}-b_{i})$ and $\bar{\tilde{\lambda}}_{k}=col(\tilde{\lambda}_{1,k},\cdots,\tilde{\lambda}_{N,k})$. The compact form of Step 2 is

[TABLE]

Noticing that $\frac{1}{\rho}\lambda_{i,k+1}+\frac{\rho-1}{\rho}\lambda_{i,k}=\tilde{\lambda}_{i,k}$, we have $s_{i,k}=\tilde{\lambda}_{i,k}+H_{i}(A_{i}\tilde{x}_{i,k}+\sum_{l\in\mathcal{E}_{i}}V_{il}z_{l,k}-b_{i})$. Denote $\bar{s}_{k}=col(s_{1,k},\cdots,s_{N,k})$, then $\bar{s}_{k}=\tilde{\bar{\lambda}}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})$. Denote $\tilde{z}_{l,k}=z_{l,k}-W_{i}[s_{j,k}-s_{i,k}],e_{l}=i\rightarrow j,$ then $z_{l,k+1}=z_{l,k}+\rho(\tilde{z}_{l,k}-z_{l,k}).$ Denote $\tilde{\mathbf{Z}}_{k}=col(\tilde{z}_{1,k},\cdots,\tilde{z}_{M,k})$, then $\tilde{\mathbf{Z}}_{k}=\mathbf{Z}_{k}-W\bar{V}^{T}\bar{s}_{k}$. The updates $z_{l,k}$ are in compact form

[TABLE]

Using (14), (15), (16), Algorithm 1 is written compactly as

[TABLE]

We verify next that any limiting point of Algorithm 1, or (17), is a zero of operator $\mathfrak{M}^{e}$, (13). Since $\{\mu_{k}\}$ satisfies $\sum_{k=1}^{\infty}\mu_{k}<\infty$ and $\mu_{k}\geq 0$, we have $\mu_{k}\rightarrow 0$ as $k\rightarrow\infty$. Assume (17) has a limiting point $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$, then we have $\mathbf{x}_{k+1}=\mathbf{x}_{k}=\tilde{\mathbf{x}}_{k}=\hat{\mathbf{x}}_{k}=\mathbf{x}^{*}$, $\bar{\lambda}_{k+1}=\bar{\lambda}_{k}=\tilde{\bar{\lambda}}_{k}=\bar{\lambda}^{*}$, and $\mathbf{Z}_{k+1}=\mathbf{Z}_{k}=\tilde{\mathbf{Z}}_{k}=\mathbf{Z}^{*}$. By (17), $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ satisfies

[TABLE]

Since $H$, $R$ and $W$ are symmetric positive definite, (19) implies that $\mathbf{0}=\Lambda{\mathbf{x}^{*}}+\bar{V}\mathbf{Z}^{*}-\bar{b}$, i.e., $A_{i}x_{i}^{*}-b_{i}=\sum_{l\in\mathcal{E}_{i}}V_{il}z^{*}_{l},\forall i\in\mathcal{N}$. Then (18) and (20) imply $\mathbf{0}\in\Lambda^{T}\bar{\lambda}^{*}+(N_{\bar{\Omega}}+\partial F)({\mathbf{x}}^{*})$ and $\mathbf{0}=\bar{V}^{T}{\bar{\lambda}^{*}}$. Using (13) for operator $\mathfrak{M}^{e}$, it follows that any limit point of Algorithm 1 belongs to $zer\mathfrak{M}^{e}$.

We show that any $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})\in zer\mathfrak{M}^{e}$ has $\mathbf{x}^{*}$ as a variational GNE of game (1). Since $\mathcal{G}$ is undirected and connected, $\bar{V}^{T}\bar{\lambda}^{*}=\mathbf{0}$ implies $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$, $\lambda^{*}\in\mathbf{R}^{m}$. Using $\mathbf{1}^{T}_{N}V=\mathbf{0}^{T}_{M}$ and $-\Lambda\mathbf{x}^{*}-\bar{V}\mathbf{Z}^{*}+\bar{b}=\mathbf{0}$, $\mathbf{1}_{N}^{T}\otimes I_{m}(-\Lambda\mathbf{x}^{*}-V\otimes I_{m}\mathbf{Z}^{*}+\bar{b})=\mathbf{0}$ implies $\sum_{i=1}^{N}A_{i}x^{*}_{i}=\sum_{i=1}^{N}b_{i}.$ Moreover, $\mathbf{0}\in\Lambda^{T}\bar{\lambda}^{*}+(N_{\bar{\Omega}}+\partial F)(\mathbf{x}^{*})$ and $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$ imply $\mathbf{0}\in A_{i}^{T}\lambda^{*}+N_{\Omega_{i}}(x_{i}^{*})+\partial_{i}f_{i}(x^{*}_{i},\mathbf{x}^{*}_{-i}),\;\forall i\in\mathcal{N}.$ Therefore, $\mathbf{x}^{*}$ and $\lambda^{*}$ satisfy the KKT condition (7) for the GVI (5), hence, $\mathbf{x}^{*}$ is a variational GNE of game (1) with $X=X^{e}$. , and all players have the same local multipliers, i.e., $\lambda_{i}^{*}=\lambda^{*},\forall i\in\mathcal{N}$.

$\Box$

IV-C Proximal parallel splitting algorithm for $X=X^{i}$

The distributed variational GNE computation algorithm for game (1) when $X=X^{i}$ is given as follows.

Algorithm 2

Step 1a–update of $x_{i,k}$:*

• •

Construct a subgame where player $i$ has a decision $x_{i}\in\Omega_{i}$ and an objective function $\tilde{f}_{i}(x_{i},\mathbf{x}_{-i})$,

[TABLE]

and denote its NE by $\hat{\mathbf{x}}_{k}=col(\hat{x}_{1,k},\cdots,\hat{x}_{N,k})$.

• •

Compute an inexact solution $\tilde{\mathbf{x}}_{k}=col(\tilde{x}_{1,k},\cdots,\tilde{x}_{N,k})$ to game in (21) such that $||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||\leq\mu_{k}$.

• •

Player $i$ updates its decision $x_{i,k}$ with

[TABLE]

Step 1b–update of $z_{l,k}$: If $e_{l}\in\mathcal{E}_{i}^{out}$, then player $i$ receives $\lambda_{j,k}$, $j\in\mathcal{N}_{l}\setminus\{i\}$, and updates $z_{l,k}$ with

[TABLE]

Step 2–update of $\lambda_{i,k}$: player $i$ receives $\tilde{z}_{l,k},z_{l,k}$, $e_{l}\in\mathcal{E}_{i}^{in}$.

[TABLE]

All variables have the same meaning as in Algorithm 1.*

Remark 3

Compared with Algorithm 1, Algorithm 2 has a different update order, that is \left\{\begin{array}[]{l}\mathbf{x}\\ \mathbf{Z}\end{array}\right.\rightarrow\bar{\lambda} rather than $\mathbf{x}\rightarrow\bar{\lambda}\rightarrow\mathbf{Z}$. Algorithm 2 is called a proximal parallel splitting algorithm since $\mathbf{x}$ and $\mathbf{Z}$ can be updated in parallel, and only the update of $\lambda$ utilizes the most recent information. Another difference lies in the construction of subgame at Step 1a, i.e., (21) only utilizes the proximal term and Lagrangian term without considering a (linearized) quadratic penalty term. The quadratic term in proximal ADMM is motivated by the augmented Lagrangian method for equality constrained optimization. However, the augmented Lagrangian method for inequality constrained optimization is less understood and may involve non-differentiable terms. Hence, augmented Lagrangian methods for distributed GNE computation of inequality constrained games is beyond the scope of this paper.

Using the same compact notations in subsection IV-B, such as $\bar{\lambda}$, $\mathbf{Z}$, $R$, $W$, $H$, $\bar{V}$, $\Lambda$, $\Lambda^{T}$, and $\bar{b}$, we give next the limiting point analysis of Algorithm 2.

Theorem 2

Suppose that Assumption 1, 2 hold for the game (1) when $X=X^{i}$. Then any limiting point $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ of Algorithm 2 belongs to the zeros of operator $\mathfrak{M}^{i}$ defined by

[TABLE]

Meanwhile, any zero $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ of ${\mathfrak{M}}^{i}$, (25) has the $\mathbf{x}^{*}$ component as a variational GNE of game in (1) when $X=X^{i}$.

Proof: We first write Algorithm 2 in a compact form. Since $\hat{\mathbf{x}}_{k}$ is an NE of subgame (21), $\hat{x}_{i,k}$ is an optimal solution to $\min_{x_{i}\in\Omega_{i}}\;f_{i}(x_{i},\hat{\mathbf{x}}_{-i,k})+\frac{1}{2}||x_{i}-x_{i,k}||_{R_{i}}^{2}+\lambda_{i,k}^{T}A_{i}x_{i}.$ Under Assumption 1 and 2, its optimality condition is

$\mathbf{0}\in N_{\Omega_{i}}(\hat{x}_{i,k})+\partial_{i}f_{i}(\hat{x}_{i,k},\hat{\mathbf{x}}_{-i,k})+R_{i}(\hat{x}_{i,k}-x_{i,k})+A_{i}^{T}\lambda_{i,k}$.

On the other hand, $\tilde{x}=P^{G}_{\Omega}(x)=\arg\min_{y}(\iota_{\Omega}(y)+\frac{1}{2}||x-y||_{G}^{2})$ if and only if $0\in N_{\Omega}(\tilde{x})+G(\tilde{x}-x)$. Therefore, the first line of (24) can be written as $\mathbf{0}\in N_{\mathbf{R}^{m}_{+}}(\tilde{\lambda}_{i,k})+H_{i}^{-1}\{\tilde{\lambda}_{i,k}-\lambda_{i,k}-H_{i}[A_{i}(2\tilde{x}_{i,k}-x_{i,k})+\sum_{l\in\mathcal{E}_{i}}V_{il}(2\tilde{z}_{l,k}-z_{l,k})-b_{i}]\}$. Hence, for all players we can write in compact form,

[TABLE]

Since $R,H,W$ are positive definite, with similar arguments in Theorem 1, it can be verified that any limiting point of (26) is a zero of ${\mathfrak{M}}^{i}$ in (25).

Suppose $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ is a zero of ${\mathfrak{M}}^{i}$ in (25). Then with similar arguments as in Theorem 1, we obtain $\bar{\lambda}^{*}$ to be $1_{N}\otimes\lambda^{*},\lambda^{*}\in\mathbf{R}^{m}$. And $\mathbf{x}^{*}$ together with $\lambda^{*}$ satisfies the first line of (8). Moreover, by $\mathbf{0}\in-\Lambda\mathbf{x}^{*}-\bar{V}\mathbf{Z}^{*}+\bar{b}+N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda}^{*})$ and $\bar{\lambda}^{*}=1_{N}\otimes\lambda^{*}$, there exist $v_{1},v_{2},\cdots,v_{N}\in N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$, such that $\mathbf{0}_{mN}=-\Lambda\mathbf{x}^{*}-{V}\otimes I_{m}\mathbf{Z}^{*}+\bar{b}+col(v_{1},\cdots,v_{N}).$ Multiplying both sides of above equation with $\mathbf{1}^{T}_{N}\otimes I_{m}$ and combining with $\mathbf{1}^{T}V=\mathbf{0}^{T}$, we have $\mathbf{0}_{m}=-\sum_{i=1}^{N}(A_{i}x^{*}_{i}-b_{i})+\sum_{i=1}^{N}v_{i}$ We have $\sum_{i=1}^{N}v_{i}\in N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$ due to $v_{i}\in N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$ and $N_{\mathbf{R}^{m}_{+}}(\lambda^{*})$ is a convex cone. This implies that the second line of KKT condition (8) is satisfied. The conclusion follows. $\Box$

Remark 4

Both Algorithm 1 and Algorithm 2 are double-layer algorithms since at each outer-layer iteration, players need to compute inexactly an NE of regularized subgames with a given accuracy. Since only monotonicity is assumed here, various problems could be solved with our algorithms, but the specified choice of the inner-layer algorithm should be determined according to the problem at hand. Thus, Algorithm 1 and Algorithm 2 are “prototype” algorithms. The inner-layer NE seeking algorithm is not specified for the following reasons.

• •

The NE seeking algorithm can and should be tailored according to the structure of the objective functions, such as the splitting form in Remark 1. For example, if $\partial F(\mathbf{x})=\partial L(\mathbf{x})+\nabla_{p}G(\mathbf{x})$ as in Remark 1, the possible Lipschitz continuity of $\nabla_{p}G(\mathbf{x})$ should be considered when choosing the NE algorithm.

• •

The subgames are regularized to have strongly monotone PS/PGs due to the proximal term $\frac{1}{2}||x_{i}-x_{i,k}||_{R_{i}}^{2}$, hence efficient NE seeking algorithms available in the literature can be used, e.g., **[2, 20, 18, 4, 14]**.

For example, if the objective functions satisfy the assumptions in [2], the asynchronous distributed best-response algorithm in [2] could be adopted for NE seeking. Particularly, denote $\mathcal{B}(\mathbf{x})=col(\mathcal{B}_{1}(\mathbf{x}_{-1}),\cdots,\mathcal{B}_{N}(\mathbf{x}_{-N}))$ where $\mathcal{B}_{i}(\mathbf{x}_{-i})=\arg\min_{x_{i}\in\Omega_{i}}\tilde{f}_{i}(x_{i},\mathbf{x}_{-i})$ given fixed $\mathbf{x}_{-i}$. Using Lemma 14 of [2], $R_{i}$ can be chosen such that $\mathcal{B}(\mathbf{x})$ is a contractive map, hence the best-response algorithm enjoys a geometric convergence rate. This is even more preferable if $\mathcal{B}(\mathbf{x})$ has a closed form. If it does not, since $\tilde{f}_{i}(x_{i},\mathbf{x}_{-i})$ is strongly convex in $x_{i}$ given $\mathbf{x}_{-i}$ due to $\frac{1}{2}||x_{i}-x_{i,k}||_{R_{i}}^{2}$, $\mathcal{B}_{i}(\mathbf{x}_{-i})$ can be computed locally with the proximal gradient method that also enjoys a geometric convergence rate.

The stopping criterion for the inner layer should be decided after the NE seeking algorithm is selected. For example, for the best-response algorithm in [2] a termination criterion to meet a given solution accuracy can be determined as in Remark 18 of [2].

We note that a single-layer GNE seeking algorithm has been proposed in [15], but uses diminishing step-sizes and a coordination center. Our double-layer GNE algorithm could be preferable when there is no central node and the subgames can be easily solved.

Remark 5

The challenges involved in GNE seeking of game (1) are as follows. Firstly, the game has monotone PS without Lipschitz continuity (or the Lipschitz constant is not known prior). Secondly, the players can only communicate peer-to-peer to coordinate to ensure coupling constraints, even though neither $X^{e}$ nor $X^{i}$ is available to any agent. The key idea of the proposed algorithms, i.e., Algorithm 1 and 2, is to decompose the complicated GNE seeking into sequential NE computation of regularized subgames and local coordinations. Notice that double-layer algorithms have been adopted for GNE seeking in [3] and [9], but only for strongly monotone games. The proximal terms regularize the subgame such that its NE can be much easier computed. The edge variables, motivated by network flow, [29], are introduced to assist agents to reach consensus on local multipliers and to satisfy the coupling constraints.

V Convergence analysis

In this section, we first show that both Algorithm 1 and 2 can be derived from a preconditioned proximal point algorithm (PPPA) for finding zeros of monotone operators. Then, based on this relationship we prove their convergence under a sufficient choice for the parameters $R_{i}$, $H_{i}$ and $W_{l}$, $\forall i\in\mathcal{N},l\in\mathcal{E}$.

Given a maximally monotone operator ${\mathfrak{M}}$ and a symmetric positive definite matrix $\Phi$, the inexact PPPA with relaxation steps for finding a zero of ${\mathfrak{M}}$ is given below.

Algorithm 3

**

[TABLE]

*where $\nu_{k}>0$, $\sum_{k=1}^{\infty}\nu_{k}<\infty$, and $\rho\in[1,2)$. * **

Remark 6

The proximal algorithm for solving $\mathbf{0}\in\mathfrak{M}(x)$ (referring to Theorem 23.41 of [28]) is

[TABLE]

which can be equivalently written as $\varpi_{k}-\varpi_{k+1}\in\mathfrak{M}\varpi_{k+1}$. Intuitively speaking, when $\mathfrak{M}(\varpi)$ is a linear operator $\mathfrak{M}\varpi$, each iteration of (28) involves computing an inverse of $I+\mathfrak{M}$. Hence, compared with (28), Algorithm 3 introduces a preconditioning matrix $\Phi$, considers the inexactness when evaluating the resolvent of $\mathfrak{M}$ at some specified point, and adopts an extrapolation/relexation step.

Particularly, the preconditioning matrix $\Phi$ plays a crucial role in our algorithm design:

• •

It adds a proximal term $\frac{1}{2}||x_{i}-x_{i,k}||_{R_{i}}^{2}$ to (• ‣ 1) and (21) that regularizes the subgames.

• •

It helps to compute the resolvent of the linear parts of ${\mathfrak{M}}^{e}$ and ${\mathfrak{M}}^{i}$ with just one step of local communication and local computation, without any matrix inverse.

The next result shows the convergence of Algorithm 3.

Theorem 3

Suppose ${\mathfrak{M}}$ is maximally monotone, and $\Phi$ is symmetric positive definite. Suppose $\varpi_{k}$ is generated by PPPA Algorithm 3 with $\sum_{k=1}^{\infty}\nu_{k}<\infty,$ $\rho\in[1,2)$. Then $\varpi_{k}$ converges to $\varpi^{*}$ and $\varpi^{*}\in zer{\mathfrak{M}}$.

The proof of Theorem 3 is adapted from [30], and can be found in the Appendix.

In the next two subsections, we show the convergence of Algorithm 1 and 2 by relating them to PPPA Algorithm 3, for appropriately chosen monotone operators and preconditioning matrices, and by using Theorem 3.

V-A Convergence analysis for $X=X^{e}$

We introduce two auxiliary variables $\eta\in\mathbf{R}^{mN}$ and $\theta\in\mathbf{R}^{mN}$ and denote $\varpi=col(\mathbf{x},\eta,\mathbf{Z},\theta)$. Consider another operator $\bar{\mathfrak{M}}^{e}$ related to $\mathfrak{M}^{e}$ in (13), defined as $\bar{\mathfrak{M}}^{e}:\varpi\mapsto$

[TABLE]

Define a preconditioning matrix $\Phi^{e}$,

[TABLE]

where $W=diag\{W_{1},\cdots,W_{M}\}$, $H=diag\{H_{1},\cdots,H_{N}\}$.

The following result relates Algorithm 1 to the PPPA Algorithm 3 for $\mathfrak{M}=\bar{\mathfrak{M}}^{e}$ and $\Phi=\Phi^{e}$.

Theorem 4

Suppose Assumption 1 and 2 hold. Denote $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$, $\hat{\mathbf{x}}_{k}$ and $col(\tilde{\mathbf{x}}_{k},\tilde{\mathbf{Z}}_{k},\tilde{\bar{\lambda}}_{k})$ as points generated by Algorithm 1 for initial points $\mathbf{x}_{0},\mathbf{Z}_{0},\bar{\lambda}_{0}$. Denote $\varpi_{k}=col(\mathbf{x}^{{}^{\prime}}_{k},\eta_{k},\mathbf{Z}^{{}^{\prime}}_{k},\theta_{k})$, $\hat{\varpi}_{k}=col(\hat{\mathbf{x}}^{{}^{\prime}}_{k},\hat{\eta}_{k},\hat{\mathbf{Z}}^{{}^{\prime}}_{k},\hat{\theta}_{k})$, and $\tilde{\varpi}_{k}=col(\tilde{\mathbf{x}}^{{}^{\prime}}_{k},\tilde{\eta}_{k},\tilde{\mathbf{Z}}^{{}^{\prime}}_{k},\tilde{\theta}_{k})$ as the points generated by the PPPA Algorithm 3 with $\mathfrak{M}=\bar{\mathfrak{M}}^{e}$ and $\Phi=\Phi^{e}$ for initial points $\mathbf{x}^{{}^{\prime}}_{0}=\mathbf{x}_{0},\eta_{0}=\bar{\lambda}_{0}+H(\Lambda\mathbf{x}_{0}+\bar{V}^{T}\bar{\lambda}_{0}-\bar{b}),\mathbf{Z}^{{}^{\prime}}_{0}=\mathbf{Z}_{0},\theta_{0}=\mathbf{0}$. Then, any sequence $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ can be derived from some sequence $\varpi_{k}=col(\mathbf{x}^{{}^{\prime}}_{k},\eta_{k},\mathbf{Z}^{{}^{\prime}}_{k},\theta_{k})$ as follows

[TABLE]

for some nonnegative sequence $\{\nu_{k}\}$ such that $\sum_{k=1}^{\infty}{\nu_{k}}<\infty$.

The proof of Theorem 4 is based on an induction argument and is given in the Appendix.

Remark 7

The standard ADMM for optimization can be derived from the Douglas-Rachford (DS) splitting method for dual optimization problems, and analyzed as a proximal-point algorithm, see [30] and [28]. For proximal ADMM, the analysis in [26] shows that the posterior second coordinate is not available when updating the first one. That is the reason why we split $\bar{\lambda}$ into $\eta$ and $\theta$, to have a higher order dynamics. The preconditioned DS splitting method recently introduced in [31], might lead to proximal ADMM. Compared to [31], our algorithm applies relaxation steps to all coordinates and considers inexactness in solving the subproblems.

We prove the convergence of Algorithm 1, by exploiting the relationship given in Theorem 4 and using Theorem 3.

Theorem 5

Suppose Assumption 1 and 2 hold for game (1) when $X=X^{e}$, and parameters (step-sizes) $R_{i},H_{i},W_{l}$ are symmetric positive definite, chosen such that $R-\Lambda^{T}H\Lambda$ and $W^{-1}-\bar{V}^{T}H\bar{V}$ are positive definite. Then, any $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ generated by Algorithm 1 converges to $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})\in zer\mathfrak{M}^{e}$. Furthermore, $\mathbf{x}^{*}$ is a variational GNE of game in (1) when $X=X^{e}$, and $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$, $\lambda^{*}\in\mathbf{R}^{m}$.

Proof: By Theorem 4, Algorithm 1 is related to PPPA Algorithm 3 for $\bar{\mathfrak{M}}^{e}$, (29), $\Phi^{e}$, (30). Convergence follows by Theorem 3 if we show that $\bar{\mathfrak{M}}^{e}$ is maximally monotone and $\Phi^{e}$ is positive definite. Denote $\varpi=col(\mathbf{x},\eta,\mathbf{Z},\theta)$, then

[TABLE]

Since $R-\Lambda^{T}H\Lambda$ and $W^{-1}-\bar{V}^{T}H\bar{V}$ are positive definite, it follows immediately that $\Phi^{e}$ is positive definite.

Operator $\bar{\mathfrak{M}}^{e}$, (29), is written as the sum of two operators. The first is a skew-symmetric linear operator, hence, is maximally monotone with domain of whole space. $N_{\bar{\Omega}}$ is maximally monotone as a normal cone operator of a closed convex set, and $\partial F(\mathbf{x})$ is also maximally monotone by Assumption 1. Since their domains coincide, $N_{\bar{\Omega}}+\partial F$ is maximally monotone, and the 2nd term in (29) is maximally monotone as the Cartesian product of maximally monotone operators.

By Theorem 4, for any sequence $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ generated from Algorithm 1, we can find $\varpi_{k}=col(\mathbf{x}^{{}^{\prime}}_{k},\eta_{k},\mathbf{Z}^{{}^{\prime}}_{k},\theta_{k})$ generated from Algorithm 3 such that (31) holds for all $k$ and $\sum_{k=1}^{\infty}\nu_{k}<\infty$. By Theorem 3, $\varpi_{k}$ converges to $\varpi^{*}=col(\mathbf{x}^{{}^{\prime}*},\eta^{*},\mathbf{Z}^{{}^{\prime}*},\theta^{*})$ and $\varpi^{*}\in zer\bar{\mathfrak{M}}^{e}$. By (31), $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ also converges to $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})$ such that $\mathbf{x}^{*}=\mathbf{x}^{{}^{\prime}*}$, $\mathbf{Z}^{*}=\mathbf{Z}^{{}^{\prime}*}$, and $\bar{\lambda}^{*}=\eta^{*}-\theta^{*}-H(\Lambda\mathbf{x}^{{}^{\prime}*}+\bar{V}\mathbf{Z}^{{}^{\prime}*}-\bar{b})$.

Since $\varpi^{*}=col(\mathbf{x}^{{}^{\prime}*},\eta^{*},\mathbf{Z}^{{}^{\prime}*},\theta^{*})$ is a zero of $\bar{\mathfrak{M}}^{e}$, (29), we have $\Lambda\mathbf{x}^{{}^{\prime}*}+\bar{V}\mathbf{Z}^{{}^{\prime}*}-\bar{b}=\mathbf{0}$, so that $\bar{\lambda}^{*}=\eta^{*}-\theta^{*}$, and

[TABLE]

Using $\mathbf{x}^{{}^{\prime}*}=\mathbf{x}^{*}$, $\mathbf{Z}^{{}^{\prime}*}=\mathbf{Z}^{*}$ and the definition of $\mathfrak{M}^{e}$ in (13), it follows that $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ generated from Algorithm 1 converges to $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})\in zer\mathfrak{M}^{e}$. By Theorem 1, $\mathbf{x}_{k}$ converges to $\mathbf{x}^{*}$, a variational GNE of the game in (1), and players’ local multipliers converge to the same $\lambda^{*}$, which together with $\mathbf{x}^{*}$ satisfies KKT condition in (7). $\Box$

Remark 8

If $H_{i}$ is chosen to be a diagonal positive matrix, $R_{i}$ and $W_{l}$ can be chosen using diagonally dominance to ensure $R-\Lambda^{T}H\Lambda$ and $W^{-1}-\bar{V}^{T}H\bar{V}$ are positive definite. In this case, the parameters $R_{i},H_{i}$ and $W_{l}$ can be chosen independently by player $i$ with just local data and computation.

V-B Convergence analysis for $X=X^{i}$

The next result shows the convergence of Algorithm 2.

Theorem 6

Suppose Assumption 1 and 2 hold for game (1) when $X=X^{i}$, and parameters $R_{i},H_{i},W_{l}$ are symmetric positive definite, such that the matrix $\Phi^{i}$ is positive definite,

[TABLE]

Then any $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ generated by Algorithm 2 converges to $col(\mathbf{x}^{*},\mathbf{Z}^{*},\bar{\lambda}^{*})\in zer{\mathfrak{M}}^{i}$ in (25). Furthermore, $\mathbf{x}^{*}$ is a variational GNE of game in (1) and $\bar{\lambda}^{*}=\mathbf{1}_{N}\otimes\lambda^{*}$, $\lambda^{*}\in\mathbf{R}_{+}^{m}$.

Proof: Consider the PPPA Algorithm 3 with $\varpi=col(\mathbf{x},\mathbf{Z},\bar{\lambda})$, for $\Phi=\Phi^{i}$ and ${\mathfrak{M}}={\mathfrak{M}}^{i}$ and $\nu_{k}=\mu_{k}$. After manipulations, the PPPA algorithm gives (26). Hence, Algorithm 2 can be derived from Algorithm 3 via a one-to-one correspondence relation. Notice that ${\mathfrak{M}}^{i}$ in (25) can be written as the sum of a skew-symmetric linear operator and a product of $(N_{\bar{\Omega}}+\partial F)\mathbf{x}\times\mathbf{0}\times N_{\mathbf{R}^{mN}_{+}}(\bar{\lambda})$. Under Assumption 1 and 2, with similar arguments as in Theorem 5, we can show that ${\mathfrak{M}}^{i}$ in (25) is maximally monotone. Since $\Phi^{i}$ is symmetric positive definite, by Theorem 3, PPPA Algorithm 3 converges. Therefore, Algorithm 2 converges to a zero of $\mathfrak{M}^{i}$, and the conclusion follows by invoking Theorem 2. $\Box$

Remark 9

Our recent work [24] considers GNE computation for games with inequality affine constraints, but assumes a strongly monotone and Lipschitz continuous PG, with inertial steps for possible acceleration. In this paper, we only assume a monotone PS, consider the inexactness when solving subproblems, and use relaxation steps for possible acceleration. Moreover, as seen in the convergence analysis, both Algorithm 1 and 2 can be regarded as fixed-point iterations for averaged operators, hence the convergence rate for fixed-point residuals could be derived based on an analysis as in [32].

VI Application and simulation studies

VI-A Rate control game over wireless ad-hoc networks

This example is adapted from [21]. Consider a wireless ad-hoc network (WANET) with $16$ nodes and $16$ links $\{L_{1},\cdots,L_{16}\}$ as shown in Fig. 1. There are $15$ users $\{U_{1},...,U_{15}\}$ who want to transfer data through the links. $R_{i}$ is the path adopted by user $U_{i}$, and $L_{j}\in R_{i}$ if user $U_{i}$ transfers data through link $L_{j}$. User $U_{i}$ decides its data rate $x_{i}$, and should satisfy a local constraint $0\leq x_{i}\leq B_{i}$. In Fig. 1, the solid lines represent the links $\{L_{1},\cdots,L_{16}\}$, and dashed line displays each path $R_{i}$. Denote $A=[A_{1},\cdots,A_{15}]\in\mathbf{R}^{16\times 15}$ where $A_{i}\in\mathbf{R}^{16}$, and $A_{i}$ has its $j$th element to be $1$ if $U_{i}$ uses $L_{j}$ and to be [math], otherwise. Link $L_{j}$ has a maximal capacity $C_{j}>0$. Denote $C=col(C_{1},\cdots,C_{16})$, hence all users’ data rate $\mathbf{x}$ should satisfy the inequality coupling constraint $A\mathbf{x}\leq C$. The objective function of user $U_{i}$ is $f_{i}(x_{i},\mathbf{x}_{-i})=-u_{i}(x_{i})+D^{T}(\mathbf{x})A_{i}x_{i}$, where $u_{i}(x_{i})=\chi_{i}\log(x_{i}+1)$ is user $i$’s utility function, and $D(\mathbf{x})=col(d_{1}(\mathbf{x}),\cdots,d_{16}(\mathbf{x}))$ with $d_{j}(\mathbf{x})=\frac{\kappa_{j}}{C_{j}-[A\mathbf{x}]_{j}+\xi_{j}}$ maps $\mathbf{x}$ to the unit delays of each link. The parameters are randomly drawn as follows: $C_{j}\in[10,15]$, $B_{i}\in[5,10]$, $\chi_{i}\in[10,20]$, $\kappa_{j}\in[10,30]$ and $\xi_{j}\in[20,40]$, and are numerically verified to ensure Assumption 1.

We use Algorithm 2. Each player has a local $C_{i}=\frac{1}{15}C$, and has local step-sizes $R_{i}=10$, $H_{i}=0.5I_{16}$, $W_{l}=0.5I_{16}$ and $\rho=1.1$. Players communicate over the graph in Fig. 1, with edges arbitrarily ordered. The initial point ${x}_{i,0}$ is randomly chosen within $[0,B_{i}]$, and initial $\lambda_{i}$, $z_{l}$ are chosen to be zero. The subgames are solved using gradient methods in [20] to get the exact NE $\hat{\mathbf{x}}_{k}$, and each $\tilde{\mathbf{x}}_{k}$ is chosen as the first point satisfying $||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||<\frac{1}{k^{2}}$. The simulation results are shown in Fig. 2-4.

VI-B Task allocation game

In this part, we consider a task allocation game with $8$ tasks $\{T_{1},\cdots,T_{8}\}$ and $14$ processors (workers) $\{w_{1},\cdots,w_{14}\}$. Each task $T_{j}$ is quantified as a load of $C_{j}>0$ that should be met by the workers. Each worker $w_{i}$ decides its working output $x_{i}=col(x^{1}_{i},x^{2}_{i},x^{3}_{i},x^{4}_{i})\in\mathbf{R}^{4}$ within its capacity $\mathbf{0}\leq x_{i}\leq B_{i},B_{i}\in\mathbf{R}_{+}^{4}$. If worker $w_{i}$ allocates a part of its output to task $T_{j}$, there is an arrow $w_{i}\rightarrow T_{j}$ in Fig. 5, either blue or red. Specifically, if $w_{i}$ allocates $x^{1}_{i},x^{2}_{i}$ to $T_{j}$, there is a dashed blue arrow in Fig. 5, and if $w_{i}$ allocates $x^{3}_{i},x^{4}_{i}$ to $T_{j}$, there is a solid red arrow in Fig. 5 Define a matrix $A=[A_{1},\cdots,A_{15}]\in\mathbf{R}^{8\times 56}$ with $A_{i}=[a^{1}_{i},a_{i}^{2},a_{i}^{3},a_{i}^{4}]\in\mathbf{R}^{8\times 4}$ quantifying how the output of worker $w_{i}$ is allocated to each task. Each column $a_{i}^{k}$ has only one element being nonzero, and the $j$th element of $a_{i}^{1}$ or $a_{i}^{2}$ is nonzero if there is a dashed blue arrow $w_{i}\rightarrow T_{j}$ on Fig. 5, and the $j$th element of $a_{i}^{3}$ or $a_{1}^{4}$ is nonzero if there is a red arrow $w_{i}\rightarrow T_{j}$ on Fig. 5. The nonzero elements in $A_{i}$ are randomly chosen from $[0.5,1]$. It is required that the tasks should be met by the working output of the players. Denote $C=col(C_{1},\cdots,C_{8})$, then the workers have an equality coupling constraint: $A\mathbf{x}=C$. The objective function of player (worker) $w_{i}$ is $f_{i}(x_{i},\mathbf{x}_{-i})=c_{i}(x_{i})-R^{T}(\mathbf{x})A_{i}x_{i}$. Here, $c_{i}(x_{i})$ is a cost function of worker $w_{i}$ and is taken as $c_{i}(x_{i})=\sum_{s=1}^{4}\max\{q^{s}_{i}{x^{s}_{i}}^{2}-\xi^{s}_{i}x^{s}_{i},l^{s}_{i}x^{s}_{i}\}+(p_{i}^{T}x_{i}-d_{i})^{2}+x_{i}^{T}S_{i}x_{i}$. $R(\mathbf{x})=col(R_{1}(\mathbf{x}),\cdots,R_{8}(\mathbf{x}))$ is a vector function that maps the workers’ output to the award price of each task, and $R_{j}(\mathbf{x})=\kappa_{j}-\chi_{j}\log([A\mathbf{x}]_{j}+1)$. Parameters of the problem are randomly drawn as follows: $C_{j}\in[1,2]$, $\chi_{j}\in[0.1,0.6]$, $\kappa_{j}\in[10,20]$, $q^{s}_{i}\in[1,2]$, $\xi^{s}_{i}\in[6,12]$, $d_{i}\in[1,2]$, and $l^{s}_{i}\in[1,3]$. $p_{i}\in\mathbf{R}^{4}$ is a randomly generated stochastic vector, $S_{i}\in\mathbf{R}^{4\times 4}$ is a randomly generated positive definite matrix, and each element of $B_{i}$ is drawn from $[1,3]$. The parameters are numerically checked to ensure Assumption 1.

We apply Algorithm 1 to this problem, over a communication graph as in Fig. 1, without node 15 and its adjacent edge, and with the remaining edges arbitrarily ordered. Each player has a local $C_{i}=\frac{1}{15}C$, and local step-sizes $R_{i}$, $H_{i}$, $W_{l}$ that are all diagonal matrices with nonzero elements uniformly drawn from $[4,8]$, $[0.2,0.4]$ and [0.2, 0.4] respectively. The relaxation step-size is taken as $\rho=1.1$. The initial ${x}_{i,0}$ is randomly chosen within $0\leq x_{i,0}\leq B_{i}$, and initial $\lambda_{i}$, $z_{l}$ are chosen to be zeros. The subgames are solved using subgradient methods in [16] to get the exact NE $\hat{\mathbf{x}}_{k}$, and each $\tilde{\mathbf{x}}_{k}$ is chosen to be the first point on the trajectory satisfying $||\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k}||<\frac{1}{k^{2}}$. The simulation results are shown in Fig. 6-7.

VII Conclusions

In this paper, we considered GNE computation of monotone games with affine coupling constraints. We proposed center-free distributed algorithms for both equality and inequality constraints, based on a preconditioned proximal point algorithm. We decomposed the GNE computation into sequential NE computation of regularized subgames and local coordination of multipliers and auxiliary variables. We considered inexactness in solving the subgames and incorporated relaxation steps. We proved their convergence by resorting to the theory of proximal algorithms and averaged operators.

There are still a lot of promising open problems. Motivated by [20] and [21], it is appealing to consider distributed GNE seeking when players cannot observe all other players’ decisions. Motivated by [23, 5] and [6], center-free GNE seeking of monotone aggregative games with discrete-time algorithms is still open. It is appealing to develop asynchronous distributed GNE computation algorithms with delayed information, and consider the problem when the players interact over switching and directed communication graphs. As important is to consider computational GNE seeking algorithms together with the mechanism design which can ensure that players faithfully report their states and auxiliary variables, possibly by providing proper incentive or punishment.

Appendix

Essentially, the proof of Theorem 3 utilizes the following facts: $\Phi^{-1}{\mathfrak{M}}$ is maximally monotone under the $\Phi-$induced norm $||\cdot||_{\Phi}$; $R_{\Phi^{-1}\mathfrak{M}}$ is a $\frac{1}{2}$-averaged operator; Proposition 4.25 of [28] for averaged operators and Robbins-Siegmund lemma for sequence convergence, given as follows.

Lemma 1** (Proposition 23.7 of [28])**

If operator $\mathfrak{A}$ is maximally monotone, then $T=R_{\mathfrak{A}}=({\rm Id}+\mathfrak{A})^{-1}$ is firmly nonexpasive, and $domR_{\mathfrak{A}}=\mathbf{R}^{m}$.

Lemma 2** (Proposition 4.25 of [28])**

Given an operator $T$ and $\alpha\in(0,1)$, then $T\in\mathcal{A}(\alpha)$ is equivalent with any following statements:

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

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

Lemma 3** (Robbins-Siegmund)**

Suppose nonnegative sequences $\{\alpha_{k}\}$, $\{\beta_{k}\}$ and $\{v_{k}\}$ satisfy the recursive relations $\alpha_{k+1}\leq\alpha_{k}-\beta_{k}+v_{k},\forall k$ and $\sum_{k=1}^{\infty}v_{k}<\infty$, then $\{\alpha_{k}\}$ converges, $\sum_{k=1}^{\infty}\beta_{k}<\infty$ and $\lim_{k\rightarrow\infty}\beta_{k}=0$.

Proof of Theorem 3: $\Phi(\varpi_{k}-\hat{\varpi}_{k})\in{\mathfrak{M}}\hat{\varpi}_{k}$ implies that $\exists u_{k}\in{\mathfrak{M}}\hat{\varpi}_{k}$ such that $\Phi(\varpi_{k}-\hat{\varpi}_{k})=u_{k}$.

Since $\Phi$ is positive definite, $\varpi_{k}-\hat{\varpi}_{k}=\Phi^{-1}u_{k}$. That is $\varpi_{k}-\hat{\varpi}_{k}\in\Phi^{-1}{\mathfrak{M}}\hat{\varpi}_{k}$. Since $\Phi$ is positive definite and ${\mathfrak{M}}$ is maximally monotone, we have $\Phi^{-1}{\mathfrak{M}}$ is maximally monotone under the $\Phi-$induced norm $||\cdot||_{\Phi}$. In fact, $\Phi$ is positive definite and nonsingular. For any $(x,u)\in gra\Phi^{-1}{\mathfrak{M}}$ and $(y,v)\in gra\Phi^{-1}{\mathfrak{M}}$, $\Phi u\in\Phi\Phi^{-1}{\mathfrak{M}}(x)\in{\mathfrak{M}}(x)$ and $\Phi v\in\Phi\Phi^{-1}{\mathfrak{M}}(y)\in{\mathfrak{M}}(y)$. Then $\langle x-y,u-v\rangle_{\Phi}=\langle x-y,\Phi(u-v)\rangle\geq 0,\forall x,y\in dom\mathfrak{M}$, since ${\mathfrak{M}}$ is monotone. Therefore, $\Phi^{-1}{\mathfrak{M}}$ is monotone under the $\Phi-$induced inner product $\langle\cdot,\cdot\rangle_{\Phi}$. Furthermore, take $(y,v)$ with $y\in dom\mathfrak{M}$, and $\langle x-y,u-v\rangle_{\Phi}\geq 0$, for any other $(x,u)\in gra(\Phi^{-1}{\mathfrak{M}})$. For any $(x,\tilde{u})\in gra{\mathfrak{M}}$, we have $(x,\Phi^{-1}\tilde{u})\in gra(\Phi^{-1}{\mathfrak{M}})$. $\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 ${\mathfrak{M}}$ is maximally monotone, then $(y,\Phi v)\in gra\bar{\mathfrak{M}}$. We conclude that $v\in\Phi^{-1}\bar{\mathfrak{B}}(y)$ which implies that $\Phi^{-1}{\mathfrak{M}}$ is maximally monotone under $||\cdot||_{\Phi}$. In the later proof, we will use $||\cdot||$ for $||\cdot||_{\Phi}$.

Therefore, $\hat{\varpi}_{k}=({\rm Id}+\Phi^{-1}{\mathfrak{M}})^{-1}\varpi_{k}$. Denote $T=({\rm Id}+\Phi^{-1}{\mathfrak{M}})^{-1}$, then $T$ a firmly nonexpansive operator by Lemma 1. In other words, their exists a nonexpansive operator $T^{{}^{\prime}}$ such that $T=\frac{1}{2}{\rm Id}+\frac{1}{2}T^{{}^{\prime}}$. Hence $\hat{\varpi}_{k}=T\varpi_{k}=\frac{1}{2}\varpi_{k}+\frac{1}{2}T^{{}^{\prime}}\varpi_{k}$. Moreover, given any $\varpi^{*}\in zer\Phi^{-1}{\mathfrak{M}}$, or equivalently, $\varpi^{*}\in zer{\mathfrak{M}}$, $\varpi^{*}$ is a fixed point of $T$ and $T^{{}^{\prime}}$, i.e., $T\varpi^{*}=\varpi^{*}$ and $T^{{}^{\prime}}\varpi^{*}=\varpi^{*}$, with the definition of resolvent.

Denote $\breve{\varpi}_{k+1}=\varpi_{k}+\rho(\hat{\varpi}_{k}-\varpi_{k})$. We have $\breve{\varpi}_{k+1}=\varpi_{k}+\rho(\frac{1}{2}\varpi_{k}+\frac{1}{2}T^{{}^{\prime}}\varpi_{k}-\varpi_{k})=(1-\frac{\rho}{2})\varpi_{k}+\frac{\rho}{2}T^{{}^{\prime}}\varpi_{k}$. Denote $\tilde{T}=(1-\frac{\rho}{2}){\rm Id}+\frac{\rho}{2}T^{{}^{\prime}}$, then $\tilde{T}\in\mathcal{A}(\frac{\rho}{2})$ since $\rho\in[1,2)$. Moreover, given any $\varpi^{*}\in zer{\mathfrak{M}}$ we have $\tilde{T}\varpi^{*}=(1-\frac{\rho}{2})\varpi^{*}+\frac{\rho}{2}T^{{}^{\prime}}\varpi^{*}=\varpi^{*}$ since $\varpi^{*}$ is a fixed point of $T^{{}^{\prime}}$.

Given any $\varpi^{*}\in zer{\mathfrak{M}}$, with (i) of Lemma (2) we have,

[TABLE]

Therefore, $||\breve{\varpi}_{k+1}-\varpi^{*}||\leq||\varpi_{k}-\varpi^{*}||$. We also have $||\breve{\varpi}_{k+1}-\varpi_{k+1}||\leq\rho\nu_{k}$ since $\breve{\varpi}_{k+1}-\varpi_{k+1}=\rho(\hat{\varpi}_{k}-\tilde{\varpi}_{k})$ and $||\hat{\varpi}_{k}-\tilde{\varpi}_{k}||\leq\nu_{k}$ due to Algorithm 3. Then by the triangle inequality

[TABLE]

Since $\sum_{k=1}^{\infty}\rho\nu_{k}<\infty$, we conclude that $\{||{\varpi}_{k}-\varpi^{*}||\}$ converges for any given $\varpi^{*}\in zer{\mathfrak{M}}$ with Lemma 3. Hence, $\{||\breve{\varpi}_{k}-\varpi^{*}||\}$ and $\{||\varpi_{k}-\varpi^{*}||\}$ are both bounded sequences, and we denote $c_{4}=\sup_{k}||\breve{\varpi}_{k}-\varpi^{*}||$.

Since $T=\frac{1}{2}{\rm Id}+\frac{1}{2}T^{{}^{\prime}}$ is firmly nonexpansive, ${\rm Id}-T=\frac{1}{2}{\rm Id}+\frac{1}{2}(-T^{{}^{\prime}})$ is also firmly nonexpansive. By (ii) of Lemma 2, ${\rm Id}-T\in\mathcal{A}(\frac{1}{2})$ if and only if $\forall\varpi_{1},\varpi_{2}\in domT,$

[TABLE]

Hence, we have

[TABLE]

where the third equality follows from $({\rm Id}-T)\varpi^{*}=\mathbf{0}$ and the last inequality follows from (34). Denote $c_{6}=(2\rho-\rho^{2})$, then we also have $||\breve{\varpi}_{k+1}-\varpi_{k+1}||^{2}\leq\rho^{2}\nu^{2}_{k}$ and

[TABLE]

We have $\sum_{k=1}^{\infty}(\rho^{2}\nu_{k}+2\rho c_{4})\nu_{k}<\infty$ due to $\sum_{k=1}^{\infty}\nu_{k}<\infty$. By Lemma 3, we conclude that $\sum_{k=1}^{\infty}||\varpi_{k}-T\varpi_{k}||^{2}<\infty$, and $\lim_{k\rightarrow\infty}\varpi_{k}-T\varpi_{k}=\mathbf{0}$.

Since $\{||\varpi_{k}-\varpi^{*}||\}$ converges, $\{\varpi_{k}\}$ is a bounded sequence. There exists a subsequence$\{\varpi_{n_{k}}\}$ that converges to $\acute{\varpi}^{*}$. Passing to limiting point of Algorithm 3, we have $T\acute{\varpi}^{*}=\acute{\varpi}^{*}$ by $\lim_{n_{k}\rightarrow\infty}T{\varpi}_{n_{k}}-\varpi_{n_{k}}=\mathbf{0}$ and (Lipschitz) continuity of $T$. Therefore, the limiting point $\acute{\varpi}^{*}$ is a fixed point of $T$ and is a zero of $\mathfrak{M}$ in (29). Setting $\varpi^{*}=\acute{\varpi}^{*}$ in (33), we have $\{||\varpi_{k}-\acute{\varpi}^{*}||\}$ is bounded and converges. Since there exists a subsequence $\{\varpi_{n_{k}}\}$ that converges to $\acute{\varpi}^{*}$, it follows that $\{||\varpi_{k}-\acute{\varpi}^{*}||\}$ converges to zero. Therefore, the whole sequence $\{\varpi_{k}\}$ generated from Algorithm 3 with any initial point converges to $\varpi^{*}$, and $\varpi^{*}\in zer{\mathfrak{M}}$. $\Box$

Proof of Theorem 4:

We first give some useful relations derived from Algorithm 3 when $\mathfrak{M}=\bar{\mathfrak{M}}^{e}$ and $\Phi=\Phi^{e}$.

Write $\Phi^{e}(\varpi_{k}-\hat{\varpi}_{k})\in\bar{\mathfrak{M}}^{e}\hat{\varpi}_{k}$ in its componentwise form

[TABLE]

Since $R,H$ and $W$ are positive definite, (35) gives

[TABLE]

Denote

[TABLE]

then by (36)

[TABLE]

Denote $\hat{\bar{\lambda}}^{{}^{\prime}}_{k}=\hat{\eta}_{k}-\hat{\theta}_{k}-H(\Lambda\hat{\mathbf{x}}^{{}^{\prime}}_{k}+\bar{V}\hat{\mathbf{Z}}^{{}^{\prime}}_{k}-\bar{b}).$ By (37) and (39), we have

[TABLE]

From (40), (42), $\hat{\bar{\lambda}}^{{}^{\prime}}_{k}=\eta_{k}-\theta_{k}+H[-\Lambda(\mathbf{x}^{{}^{\prime}}_{k}-\hat{\mathbf{x}}^{{}^{\prime}}_{k})]$.

Then by (38) and (37)

[TABLE]

Then we prove (31) by induction. Firstly, we choose $\mathbf{x}^{{}^{\prime}}_{0}=\mathbf{x}_{0},\eta_{0}=\bar{\lambda}_{0}+H(\Lambda\mathbf{x}_{0}+\bar{V}^{T}\bar{\lambda}_{0}-\bar{b}),\mathbf{Z}^{{}^{\prime}}_{0}=\mathbf{Z}_{0},\theta_{0}=\mathbf{0}$. Hence (31) holds at $k=0$.

Suppose (31) is true at time $k$, then $\bar{\lambda}^{{}^{\prime}}_{k}$ defined in (40) has $\bar{\lambda}^{{}^{\prime}}_{k}=\bar{\lambda}_{k}$. We can choose $\hat{\mathbf{x}}^{{}^{\prime}}_{k}=\hat{\mathbf{x}}_{k}$ to satisfy equation (41) due to $\mathbf{x}^{{}^{\prime}}_{k}=\mathbf{x}_{k}$, $\mathbf{Z}_{k}^{{}^{\prime}}=\mathbf{Z}_{k}$ and (17). Then we choose $\tilde{\mathbf{x}}^{{}^{\prime}}_{k}=\tilde{\mathbf{x}}_{k}$ such that $||\hat{\mathbf{x}}^{{}^{\prime}}_{k}-\tilde{\mathbf{x}}^{{}^{\prime}}_{k}||=||\hat{\mathbf{x}}_{k}-\tilde{\mathbf{x}}_{k}||\leq\mu_{k}$. Thereby, we have $\mathbf{x}^{{}^{\prime}}_{k+1}=\mathbf{x}^{{}^{\prime}}_{k}+\rho(\tilde{\mathbf{x}}^{{}^{\prime}}_{k}-\mathbf{x}^{{}^{\prime}}_{k})=\mathbf{x}_{k+1}$.

Recall that $\hat{\eta}_{k}$ and $\hat{\theta}_{k}$ are generated by (37) and (39) from $\eta_{k}$ and $\theta_{k}$. Due to (42), (31) and $\hat{\mathbf{x}}^{{}^{\prime}}_{k}=\hat{\mathbf{x}}_{k}$ we have:

[TABLE]

By (17) we also have $\tilde{\bar{\lambda}}_{k}=\bar{\lambda}_{k}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})$. Hence $\tilde{\bar{\lambda}}_{k}-\hat{\bar{\lambda}}^{{}^{\prime}}_{k}=H\Lambda(\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k})$. By (43)

[TABLE]

We choose $\tilde{\mathbf{Z}}^{{}^{\prime}}_{k}=\tilde{\mathbf{Z}}_{k}$ where $\tilde{\mathbf{Z}}_{k}=\mathbf{Z}_{k}-W\bar{V}^{T}(\tilde{\bar{\lambda}}_{k}+H(\Lambda\tilde{\mathbf{x}}_{k}+\bar{V}\mathbf{Z}_{k}-\bar{b})$ due to (17), so that

[TABLE]

Therefore, $\exists c_{1}>0$, such that $||\tilde{\mathbf{Z}}^{{}^{\prime}}_{k}-\hat{\mathbf{Z}}^{{}^{\prime}}_{k}||\leq c_{1}\mu_{k}$. Since $\tilde{\mathbf{Z}}^{{}^{\prime}}_{k}=\tilde{\mathbf{Z}}_{k}$ and $\mathbf{Z}^{{}^{\prime}}_{k}=\mathbf{Z}_{k}$, we have $\mathbf{Z}^{{}^{\prime}}_{k+1}=\mathbf{Z}^{{}^{\prime}}_{k}+\rho(\tilde{\mathbf{Z}}^{{}^{\prime}}_{k}-\mathbf{Z}^{{}^{\prime}}_{k})=\mathbf{Z}^{k+1}$.

Denote $\tilde{\bar{\lambda}}^{{}^{\prime}}_{k}=\tilde{\eta}_{k}-\tilde{\theta}_{k}-H(\Lambda\tilde{\mathbf{x}}^{{}^{\prime}}_{k}+\bar{V}\tilde{\mathbf{Z}}^{{}^{\prime}}_{k}-\bar{b})$. Then we want to find $\tilde{\eta}_{k},\tilde{\theta}_{k}$ with $||\tilde{\eta}_{k}-\hat{\eta}_{k}||\leq c_{2}\mu_{k}$, $||\tilde{\theta}_{k}-\hat{\theta}_{k}||\leq c_{3}\mu_{k}$ such that $\tilde{\bar{\lambda}}_{k}=\tilde{\bar{\lambda}}^{{}^{\prime}}_{k}$. Suppose $\tilde{\theta}_{k}$ and $\tilde{\eta}_{k}$ are chosen to ensure $\tilde{\eta}_{k}-\tilde{\theta}_{k}-H(\Lambda\tilde{\mathbf{x}}^{{}^{\prime}}_{k}+\bar{V}\tilde{\mathbf{Z}}^{{}^{\prime}}_{k}-\bar{b})=\hat{\eta}_{k}-\hat{\theta}_{k}-H(\Lambda\hat{\mathbf{x}}^{{}^{\prime}}_{k}+\bar{V}\hat{\mathbf{Z}}^{{}^{\prime}}_{k}-\bar{b})+H\Lambda(\tilde{\mathbf{x}}_{k}-\hat{\mathbf{x}}_{k})$. Then due to (17) and (44),

[TABLE]

Hence, let $\tilde{\theta}_{k}$ and $\tilde{\eta}_{k}$ be chosen as,

[TABLE]

Obviously, in this case $\exists c_{2}>0,c_{3}>0$ such that $||\tilde{\eta}_{k}-\hat{\eta}_{k}||\leq c_{2}\mu_{k}$, $||\tilde{\theta}_{k}-\hat{\theta}_{k}||\leq c_{3}\mu_{k}$. Moreover, from (47) we have $\tilde{\bar{\lambda}}^{{}^{\prime}}_{k}=\tilde{\bar{\lambda}}_{k}$. Hence, we obtain that

[TABLE]

Therefore, when (31) holds at time $k$, it also holds at time $k+1$. Thus, we have shown by induction that given sequences $col(\mathbf{x}_{k},\mathbf{Z}_{k},\bar{\lambda}_{k})$ generated from Algorithm 1 with initial points $\mathbf{x}_{0},\mathbf{Z}_{0},\bar{\lambda}_{0}$ and $\{\mu_{k}\}$, we can find sequences $\varpi_{k}=col(\mathbf{x}^{{}^{\prime}}_{k},\eta_{k},\mathbf{Z}^{{}^{\prime}}_{k},\theta_{k})$ generated from Algorithm 3 with $\nu_{k}\leq\sqrt{1+c_{1}^{2}+c_{2}^{2}+c_{3}^{2}}\mu_{k}$ such that (31) holds. Since $\sum_{k=1}^{\infty}\mu_{k}<\infty$, we have $\sum_{k=1}^{\infty}\nu_{k}<\infty$, and the conclusion follows. $\Box$