Projected Primal-Dual Gradient Flow of Augmented Lagrangian with Application to Distributed Maximization of the Algebraic Connectivity of a Network

TL;DR

This paper introduces a projected primal-dual gradient flow method for convex optimization with constraints, demonstrating its convergence and application to distributed network connectivity maximization through SDP relaxation.

Contribution

It develops a novel projected primal-dual gradient flow approach for convex problems with constraints and applies it to distributed algebraic connectivity maximization.

Findings

The proposed method converges to saddle points and optimal solutions.

Numerical examples confirm convergence and effectiveness.

The relaxation from SDP to NP is empirically validated.

Abstract

In this paper, a projected primal-dual gradient flow of augmented Lagrangian is presented to solve convex optimization problems that are not necessarily strictly convex. The optimization variables are restricted by a convex set with computable projection operation on its tangent cone as well as equality constraints. As a supplement of the analysis in \cite{niederlander2016distributed}, we show that the projected dynamical system converges to one of the saddle points and hence finding an optimal solution. Moreover, the problem of distributedly maximizing the algebraic connectivity of an undirected network by optimizing the port gains of each nodes (base stations) is considered. The original semi-definite programming (SDP) problem is relaxed into a nonlinear programming (NP) problem that will be solved by the aforementioned projected dynamical system. Numerical examples show the convergence of the aforementioned algorithm to one of the optimal solutions. The effect of the relaxation is illustrated empirically with numerical examples. A methodology is presented so that the number of iterations needed to reach the equilibrium is suppressed. Complexity per iteration of the algorithm is illustrated with numerical examples.

Figures

Equations

$$\underset{x\in K}{\mathrm{minimize}}$$

$$Ax-b=0,$$

$$\mathcal{L}(x,v)=f(x)+v^T(Ax-b),$$

$$\dot{x}$$

$$={\Pi }_K(x,-\frac{\partial {\mathcal{L}}_{\mathcal{A}}(x,v)}{\partial x}),$$

$$\dot{v}$$

$$d(x,v)=\frac{1}{2}(\Vert x-x^{\ast }{\Vert }^2+\Vert v-v^{\ast }{\Vert }^2).$$

$${\left[L_w\right]}_{ij}$$

$$\underset{{\lambda }_2,\mu ,\{w_k^{(i)}\}}{\mathrm{maximize}}$$

$${\lambda }_{2}I-\mu \text{mathbbm}1⪯\sum _{i\in \mathcal{V}}\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)}{E}_k,$$

$$\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)}=1,$$

$$w_k^{(i)}\ge 0,\forall k\in \mathcal{E},i\in \mathcal{V}.$$

$$\underset{\{{\lambda }_2^{(i)}\},\{{\mu }^{(i)}\},\{w_k^{(i)}\},\{Z^{(i)}\}}{\mathrm{maximize}}$$

$${\lambda }_2^{(i)}I-{\mu }^{(i)}\text{mathbbm}1-\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)}{E}_k$$

$$+\sum _{j\in \mathcal{N}(i)}(Z^{(i)}-Z^{(j)})⪯0,$$

$$\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)}=1,$$

$$w_k^{(i)}\ge 0,\forall k\in \mathcal{E}(i),\forall i\in \mathcal{V},$$

$$tr({\Phi }^{\ast })=1,{\Phi }^{\ast }⪰0,{\phi }_k^{(i)\ast }\ge 0,tr(\text{mathbbm}1{\Phi }^{\ast })=0,$$

$$tr(E_k{\Phi }^{\ast })-v^{(i)\ast }+{\phi }_k^{(i)\ast }=0,w_k^{(i)\ast }\ge 0,$$

$${\phi }_k^{(i)\ast }{w}_k^{(i)\ast }=0,\forall k\in \mathcal{E}(i),$$

$${\lambda }_2^{\ast }-\sum _{i\in \mathcal{V}}\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)\ast }tr(E_k{\Phi }^{\ast })=0,$$

$${\lambda }_2^{\ast }I-{\mu }^{\ast }\text{mathbbm}1-\sum _{i\in \mathcal{V}}\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)\ast }{E}_k⪯0,$$

$$\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)\ast }=1,$$

$$tr({\Phi }^{(i)\ast })=1,{\Phi }^{(i)\ast }⪰0,{\phi }_k^{(i)\ast }\ge 0,tr(\text{mathbbm}1{\Phi }^{(i)\ast })=0,$$

$$tr(E_k{\Phi }^{(i)\ast })-v^{(i)\ast }+{\phi }_k^{(i)\ast }=0,w_k^{(i)\ast }\ge 0,$$

$${\phi }_k^{(i)\ast }{w}_k^{(i)\ast }=0,\forall k\in \mathcal{E}(i),$$

$${\lambda }_2^{(i)\ast }-\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)\ast }tr(E_k{\Phi }^{(i)\ast })$$

$$+\sum _{j\in \mathcal{N}(i)}tr[(Z^{(i)}-Z^{(j)}){\Phi }^{(i)\ast }]=0,$$

$${\lambda }_2^{(i)\ast }I-{\mu }^{(i)\ast }\text{mathbbm}1-\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)\ast }{E}_k$$

$$+\sum _{j\in \mathcal{N}(i)}(Z^{(i)\ast }-Z^{(j)\ast })⪯0,$$

$$\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)\ast }=1,\sum _{j\in \mathcal{N}(i)}{\Phi }^{(i)\ast }-{\Phi }^{(j)\ast }=0,$$

$${\nabla }_X{f}_{\epsilon }(X)=[\sum _{i=1}^N{e}^{{\lambda }_i(X)/\epsilon }{]}^{-1}[\sum _{i=1}^N{e}^{{\lambda }_i(X)/\epsilon }{u}_i{u}_i^T],$$

$${\lambda }_{\max }(X)\le f_{\epsilon }(X)\le {\lambda }_{\max }(X)+\epsilon \ln N.$$

$$\underset{\{{\mu }^{(i)}\},\{w_k^{(i)}\},\{Z^{(i)}\}}{\mathrm{minimize}}$$

$$\sum _{k\in \mathcal{E}(i)}{w}_k^{(i)}=1$$

$$w_k^{(i)}\ge 0,\forall k\in \mathcal{E}(i),\forall i\in \mathcal{V}$$

$$-\sum _{i\in \mathcal{V}}{\lambda }_2^{(i)\ast }\le \sum _{i\in \mathcal{V}}{f}_{{\epsilon }_i}(X^{(i)\ast })\le -\sum _{i\in \mathcal{V}}{\lambda }_2^{(i)\ast }+\sum _{i\in \mathcal{V}}{\epsilon }_i\ln N,$$

$$-\sum _{i\in \mathcal{V}}{\lambda }_2^{(i)\ast }$$

$$\le -\sum _{i\in \mathcal{V}}{\lambda }_2^{(i)\ast }+\sum _{i\in \mathcal{V}}{\epsilon }_i\ln N,$$

$$\sum _{i\in \mathcal{V}}{\lambda }_{\max }(X^{(i)\ast })\le -\sum _{i\in \mathcal{V}}{\lambda }_2^{(i)\ast }.$$

$$-\sum _{i\in \mathcal{V}}{\lambda }_2^{(i)\ast }=-{\lambda }_2^{\ast }={\lambda }_{\max }(\sum _{i\in \mathcal{V}}{X}^{(i)\ast }).$$

$$\sum _{i\in \mathcal{V}}{\lambda }_{\max }(X^{(i)\ast })={\lambda }_{\max }(\sum _{i\in \mathcal{V}}{X}^{(i)\ast }).$$

$${\dot{\mu }}^{(i)}$$

Full text

††thanks: This work is supported by China Scholarship Council.

, , , ,

1 Introduction

When solving a convex minimization problem with strong duality, it is well-known that the optimal solution is the saddle point of the Lagrangian. Hence it is natural to consider the gradient flow of Lagrangians (also known as saddle point dynamics) where the primal variable follows the negative gradient flow while the dual variable follows the gradient flow. Gradient flow of Lagrangians is first studied by (Arrow et al., 1959), (Kose, 1956) and has been revisited by (Feijer & Paganini, 2010). (Feijer & Paganini, 2010) studies the case of strictly convex problems and provides methodologies to transform non-strictly convex problems to strictly convex problems to fit the framework. The convergence is shown by employing the invariance principle for hybrid automata. (Cherukuri et al., 2016) studies the same strictly convex problem from the perspective of projected dynamical systems and is able to show the convergence by a LaSalle-like invariant principle for Carathéodory solutions. Instead of considering discontinuous dynamics, (Dürr & Ebenbauer, 2011) proposes a smooth vector field for seeking the saddle points of strictly convex problems. (Wang & Elia, 2011) considers a strictly convex problem with equality constraints and with inequality constraints respectively. Saddle point dynamics is also used therein, however, it is worth noticing that their problem is still strictly convex. When they consider the problem with inequality constraints, logarithmic barrier function is used. Though considering nonsmooth problems, (Zeng et al., 2017) uses the projected saddle point dynamics of augmented Lagrangian whose equality constraint is the variable consensus constraint, and can be viewed as a special case of our problem. Instead of using the continuous-time saddle point dynamics, an iterative distributed augmented Lagrangian method is developed in (Chatzipanagiotis et al., 2015). In a recent work (Niederländer & Cortés, 2016) and its conference version (Niederländer et al., 2016), the authors consider the nonsmooth case of projected saddle point dynamics and the dynamics are the same as the ones in the current paper when the objective function is smooth.

In this paper, we will focus on maximizing network algebraic connectivity distributedly. In (Simonetto et al., 2013), the authors maximize the algebraic connectivity of a mobile robot network distributedly. The authors use first-order Taylor expansion to approximate the original non-convex problem and get a convex problem. A more general linear dynamics are considered and a two-step algorithm is proposed to solve the problem distributedly. It is shown in (Simonetto et al., 2013) that the algebraic connectivity is monotonically increasing with the algorithm, while the convergence to one optimal solution is not explicitly given. (Schuresko & Cortés, 2008),(Yang et al., 2010) and (Zavlanos & Pappas, 2008) focus on assuring the connectivity distributedly, while the algebraic connectivity maximization is not considered.

The main contribution of this paper is as follows. As a supplement to (Niederländer & Cortés, 2016) and its conference version Niederländer et al. (2016), we propose a novel analysis line regarding the convergence of the dynamical system to reach comparable results. Moreover, the problem of distributedly maximizing the algebraic connectivity of an undirected network by adjusting the “port gains” of each nodes (base stations) is considered. It is worth noticing that the problem motivates from a physical system and the goal is to enable each base station to compute its own optimal port gains only using its neighbours’ information, the total number of nodes $N$ and the information belonging to itself; one can not “design” the communication network according to the structure of the problem or the algorithm. (For example, (Pakazad et al., 2015)). We solve the original problem, which is an SDP, by relaxing it into an NP problem. The NP problem is not strictly convex, hence we adapt the projected saddle point dynamics method proposed in this work to solve the aforementioned NP problem. Numerical examples show that the aforementioned algorithm converges to one of the optimal solutions.

2 Preliminaries and Notations

We denote $\mathbbm{1}=\bm{1}\bm{1}^{T}$ as an $N$ dimensional all-one matrix, where $\bm{1}$ is an $N$ dimensional all-one vector. The element located on the $i$th row and and $j$th column of a matrix $A$ is denoted as $[A]_{ij}$. If matrix $A_{1}-A_{2}$ is positive semi-definite, then it will be denoted as $A_{1}\succeq A_{2}$. We use $\|\cdot\|$ to denote 2-norm of vectors. $|S|$ denotes the cardinality of set $S$. And any notation with the superscript $*$ is denoted as the optimal solution to the corresponding optimization problem. $tr(\cdot)$ denotes the trace of a matrix. $\langle\cdot,\cdot\rangle_{2}$ is denoted as the inner-product in Euclidean space and $\langle A_{1},A_{2}\rangle_{M}=tr(A_{1}A_{2})$ denotes the inner-product in $\mathcal{S}^{n}$, which is the Hilbert space of $n\times n$ symmetric matrix.

Assume $K\subset\mathbb{R}^{n}$ is a closed and convex set, the projection of a point $x$ to the set $K$ is defined as $P_{K}(x)=\arg\min_{y\in K}\|x-y\|$. For $x\in K$, $v\in\mathbb{R}^{n}$, the projection of the vector $v$ at $x$ with respect to $K$ is defined as: ( see (Nagurney & Zhang, 2012),(Brogliato et al., 2006)) $\Pi_{K}(x,v)=\lim_{\delta\rightarrow 0}\frac{P_{K}(x+\delta v)-x}{\delta}=P_{T_{K}(x)}(v)$, where $T_{K}(x)$ denotes the tangent cone of $K$ at $x$. The interior, the boundary and the closure of $K$ is denoted as $int(K)$, $\partial K$ and $cl(K)$, respectively. The set of inward normals of $K$ at $x$ is defined as n(x)=\big{\{}\gamma\>|\>\|\gamma\|=1,\langle\gamma,x-y\rangle_{2}\leq 0,\forall y\in K\big{\}}, and $\Pi_{K}(x,v)$ fulfills the following lemma:

Lemma 1** (Nagurney & Zhang (2012)).**

If $x\in int(K)$, then $\Pi_{K}(x,v)=v$; if $x\in\partial K$, then $\Pi_{K}(x,v)=v+\beta(x)n^{*}(x)$, where $n^{*}(x)=\arg\max_{n\in n(x)}\langle v,-n\rangle$ and $\beta(x)=\max\{0,\langle v,-n^{*}(x)\rangle\}$.

Let $F$ be a vector field such that $F:K\mapsto\mathbb{R}^{n}$, the projected dynamical system is given by $\dot{x}=\Pi_{K}(x,F(x))$. Note that the right hand side of above dynamics can be discontinuous on the $\partial K$. Hence given an initial value $x_{0}\in K$, the system does not necessarily have a classical solution. However, if $F(x)$ is Lipschitz continuous, then it has a unique Carathéodory solution that continuously depends on the initial value (Nagurney & Zhang, 2012).

3 Problem Formulation and Projected Saddle Point Dynamics

In this section, we consider the following optimization problem defined on $\mathbb{R}^{n}$:

[TABLE]

where $f:\mathbb{R}^{n}\mapsto\mathbb{R}$ and $A\in\mathbb{R}^{m\times n}$. $K$ is a convex set such that calculating the projection on its tangent cone is computationally cheap. $f(x)$ is a convex function but not necessarily strictly convex. It is also assumed that the gradient of $f(x)$ is locally Lipschitz continuous and the Slater’s condition holds for (1). Hence strong duality holds for (1).

The Lagrangian $\mathscr{L}:K\times\mathbb{R}^{m}\mapsto\mathbb{R}$ for the problem (1) is given by

[TABLE]

where $v\in\mathbb{R}^{m}$ is the Lagrangian multiplier of the constraint $Ax-b=0$. Since strong duality holds for (1), then $(x^{*},v^{*})$ is a saddle point of $\mathscr{L}(x,v)$ if and only if $x^{*}$ is an optimal solution to (1) and $v^{*}$ is optimal solution to its dual problem. The augmented Lagrangian $\mathscr{L}_{\mathcal{A}}:K\times\mathbb{R}^{m}\mapsto\mathbb{R}$ for (1) is given by $\mathscr{L}_{\mathcal{A}}(x,v)=f(x)+v^{T}(Ax-b)+\frac{\rho}{2}(Ax-b)^{T}(Ax-b)$, where $\rho>0$ is the damping parameter that will help to suppress the oscillation of $x$ during optimization algorithms. Without loss of generality, we choose $\rho=1$.

We propose to find the saddle point of (2) via the saddle point dynamics projected on the set $K$, i.e.,

[TABLE]

Note that it is assumed that $\nabla f(x)$ is locally Lipschitz continuous, therefore there is a unique Carathéodory solution for the dynamics (3).

4 Convergence Analysis

In this section, we analyse the convergence for (3) and start with the analysis of the equilibrium point of (3). Niederländer & Cortés (2016) consider the nonsmooth case of projected saddle point dynamics and the dynamics are the same as the ones in the current paper when the objective function is smooth. As a supplement, we propose a novel analysis line regarding the stability of the dynamical system to reach comparable results.

Proposition 2**.**

$(x^{*},v^{*})$* is a saddle point to (1) if and only if it is an equilibrium of (3).*

{pf}

Since strong duality holds for (1), the optimality conditions become necessary and sufficient conditions. The optimality condition for (1) is given by $-\nabla f(x^{*})-A^{T}v^{*}\in N_{K}(x^{*}),\quad Ax^{*}-b=0$, (Eskelinen, 2007), which implies $-\nabla f(x^{*})-A^{T}v^{*}+A^{T}(Ax^{*}-b)\in N_{K}(x^{*})$, where $N_{K}(x^{*})$ denotes the normal cone of $K$ at $x^{*}$. This implies $\Pi_{K}(x^{*},-\nabla f(x^{*})-A^{T}v^{*}-A^{T}(Ax^{*}-b))=0$, therefore, $(x^{*},v^{*})$ is an equilibrium point of (3). On the other hand, if $(x^{*},v^{*})$ is an equilibrium point of (3), it must have $-\nabla f(x^{*})-A^{T}v^{*}+A^{T}(Ax^{*}-b)\in N_{K}(x^{*})$ and $Ax^{*}-b=0$, which implies the optimality condition. ∎

Proposition 3**.**

Given an initial value $(x(0),v(0))$, where $x(0)\in K$, the trajectory of the projected dynamical system (3) asymptotically converges to one of the saddle points of (1).

{pf}

We use LaSalle invariance principle for Carathéodory solutions (Cherukuri et al., 2015) to prove the proposition. Suppose $(x^{*},v^{*})$ is a saddle point of the Lagrangian (2), namely, $x^{*}$ is the optimal solution of (1) and $v^{*}$ is the optimal solution of its dual problem. Construct the following Lyapunov function

[TABLE]

Note that $d(x,v)$ is continuously differentiable and denote the right hand side of the dynamics (3) as vector field $F$. The Lie derivative along the vector field $f(x)$ of a function $V(x)$ is defined as $\mathcal{L}_{f}V(x)=\frac{\partial V(x)}{\partial x}f(x)$. By the definition of saddle points, $\mathscr{L}(x^{*},v)\leq\mathscr{L}(x^{*},v^{*})\leq\mathscr{L}(x,v^{*})$. The Lie derivative of $d(x,v)$ along the vector field $F$ is given by $\mathcal{L}_{F}d(x,v)=(x-x^{*})^{T}\Pi_{K}(x,-\nabla f(x)-A^{T}v-A^{T}(Ax-b))+(v-v^{*})^{T}(Ax-b)$. By Lemma 1, it holds that $\Pi_{K}(x,-\nabla f(x)-A^{T}v-A^{T}(Ax-b))=-\nabla f(x)-A^{T}v-A^{T}(Ax-b)+\beta(x)n^{*}(x)$, where $\beta(x)\geq 0$ and $\langle n^{*}(x),x-y\rangle\leq 0,\>\forall y\in K$. Since $Ax^{*}-b=0$, hence it follows that $\mathcal{L}_{F}d(x,v)=(x-x^{*})^{T}(-\nabla f(x)-A^{T}v-A^{T}(Ax-b)+\beta(x)n^{*}(x))+(v-v^{*})^{T}(Ax-b)\leq(x-x^{*})^{T}(-\nabla f(x)-A^{T}v-A^{T}(Ax-b))+(v-v^{*})^{T}(Ax-b)=-(x-x^{*})^{T}\frac{\partial\mathscr{L}(x,v)}{\partial x}-(Ax-b)^{T}(Ax-b)+(v-v^{*})^{T}\frac{\partial\mathscr{L}(x,v)}{\partial v}.$

Since $\mathscr{L}(x,v)$ is convex with respect to $x$ and concave with respect to $v$, it follows from the first order property of convex and concave function (Boyd & Vandenberghe, 2004) that $\mathcal{L}_{F}d(x,v)\leq\mathscr{L}(x^{*},v)-\mathscr{L}(x,v^{*})-(Ax-b)^{T}(Ax-b)$. Since $\mathscr{L}(x^{*},v)\leq\mathscr{L}(x,v^{*})$, we have $\mathcal{L}_{F}d(x,v)\leq-(Ax-b)^{T}(Ax-b)\leq 0$. Note that $d(x,v)$ is convex and differentiable. By Proposition 2 and using the result in (Bacciotti & Ceragioli, 2004), we can conclude that any saddle point is Lyapunov stable. Recall that $f(x)$ has a locally Lipschitz continuous gradient, hence the uniqueness of the Carathéodory solution to (3) can be guaranteed. Since $d(x,v)$ is differentiable and by the definition of Carathéodory solution, $x(t),v(t)$ are absolute continuous, $d(x(t),v(t))$ is differentiable almost everywhere with respect to $t$ and $\frac{d}{dt}d(x(t),v(t))=\mathcal{L}_{F}d(x,v)$ holds almost everywhere on $t\geq 0$. Therefore $d(x(t),v(t))$ is continuous and non-increasing with respect to time. Note that $d(x,v)$ is radially unbounded, hence the set $\hat{\mathcal{S}}=\{(x,v)\in\mathbb{R}^{n}\times\mathbb{R}^{m}|d(x,v)\leq d(x(0),v(0))\}$ is a compact invariant set for the system (3). The invariance of the $\omega$-limit set (Cherukuri et al., 2015) can be proved by using the same methodology as Lemma 4.1 in (Khalil & Grizzle, 2002) which is based on the continuity and the uniqueness of the solution. Now by LaSalle invariance principle for Carathéodory solutions (Cherukuri et al., 2015), the trajectory of (3) converges to the largest invariant set in $cl(\Omega)$, where $\Omega=\{(x,v)\in\hat{\mathcal{S}}|\mathcal{L}_{F}d(x,v)=0\}$.

For $(x^{\prime},v^{\prime})\in\Omega$, we have $Ax^{\prime}-b=0$ and $\mathscr{L}(x^{*},v^{\prime})=\mathscr{L}(x^{\prime},v^{*})\Leftrightarrow f(x^{*})=f(x^{\prime})+v^{*T}(Ax^{\prime}-b)\Leftrightarrow f(x^{*})=f(x^{\prime})$. Since $x^{*}$ is an optimal solution to (1), $x^{\prime}$ is also an optimal solution to (1). Denote $\mathcal{D}=\{(x,v)\in\hat{\mathcal{S}}|x\mbox{ is an optimizer}\}$ and note that $\Omega\subset\mathcal{D}$. Hence $cl(\Omega)\subset cl(\mathcal{D})$. On the other hand, since the set of optimal solutions for a convex optimization problem is closed, $\mathcal{D}$ is also closed and hence $cl(\Omega)\subset cl(\mathcal{D})=\mathcal{D}$.

Denote $\mathcal{M}$ as the largest invariant set in $cl(\Omega)$. Assume the initial value $(\hat{x}(0),\bar{v}(0))\in\mathcal{M}\subset cl(\Omega)$, then $\hat{x}(0)$ is also an optimal solution to (1). Hence there exists some $\hat{v}$ such that $(\hat{x}(0),\hat{v})$ is a saddle point of (2). Then the Lyapunov function can be constructed similarly as $\hat{d}(x,v)=\frac{1}{2}(\|x-\hat{x}(0)\|^{2}+\|v-\hat{v}\|^{2})$. We have shown that for any arbitrarily chosen saddle point $(x^{*},v^{*})$, $d(x(t),v(t))$ is non-increasing with respect to time; such statement also holds for $\hat{d}(x(t),v(t))$. Furthermore, since $(\hat{x}(0),\bar{v}(0))\in\mathcal{M}$ and $\mathcal{M}$ is invariant, $(\hat{x}(t),\bar{v}(t))\in\mathcal{M},\forall t\geq 0$. This implies that $\hat{x}(t)$ is an optimizer for all $t\geq 0$ and hence we have $\dot{\bar{v}}(t)=A\hat{x}(t)-b=0,\forall t\geq 0$. And this implies $\bar{v}(t)=\bar{v}(0),\forall t\geq 0$ and hence $\hat{d}(\hat{x}(t),\bar{v}(t))\leq\hat{d}(\hat{x}(0),\bar{v}(0))=\frac{1}{2}\|\bar{v}(0)-\hat{v}\|^{2}$. This implies $\|\hat{x}(t)-\hat{x}(0)\|^{2}=0$ and hence $x(t)=\hat{x}(0),\forall t\geq 0$. Therefore, any trajectory that starts in $\mathcal{M}$ remains constant for all times, i.e., any point in $\mathcal{M}$ is an equilibrium point. By Proposition 2, we can conclude that they are saddle points of (1).

We have shown that given an initial value $(x(0),v(0))$, where $x(0)\in K$, the trajectory of (3) asymptotically converges to a set $\mathcal{M}$ whose elements are saddle points of (1). Now we show the trajectory asymptotically converges to a point in $\mathcal{M}$ by contradiction. To abbreviate the notation, we denote $\eta=(x^{T},v^{T})^{T}$. Suppose $\eta(t)$ does not converge to a point in $\mathcal{M}$, namely, the trajectory’s $\omega$-limit set $\Gamma(\eta(t))$ is not a singleton. This means we can choose $\bar{\eta}_{1},\bar{\eta}_{2}\in\Gamma(\eta(t))\subset\mathcal{M}$, such that $\|\bar{\eta}_{1}-\bar{\eta}_{2}\|=\zeta,\zeta>0$. Since $\bar{\eta}_{1},\bar{\eta}_{2}\in\mathcal{M}$, $\bar{\eta}_{1},\bar{\eta}_{2}$ are saddle points and we have shown that all saddle points are Lyapunov stable. This means that there exists $\delta(\zeta/2)$, such that if $\|\eta(T)-\bar{\eta}_{1}\|<\delta(\zeta/2)$, then $\|\eta(t)-\bar{\eta}_{1}\|<\zeta/2,\forall t\geq T$. Since $\bar{\eta}_{1}$ is an $\omega$-limit point in $\Gamma(\eta(t))$, there exists such $T$ so that $\|\eta(T)-\bar{\eta}_{1}\|<\delta(\zeta/2)$, and hence the trajectory can never leave $\zeta/2$-neighbourhood of $\bar{\eta}_{1}$ after time instant $T$. But $\bar{\eta}_{2}$ is also an $\omega$-limit point, there must exists a sequence of point on the trajectory that tends to it. Hence we have a contradiction and $(x(t),v(t))$ converges to a saddle point in $\mathcal{M}$. ∎

5 Distributed Algebraic Connectivity Maximization

In this section, we apply the aforementioned algorithm to maximize the algebraic connectivity of a network in a distributed manner. The problem is first formulated as a Semi-definite Programming (SDP) problem. With an equivalent formulation of the original SDP problem, the problem is relaxed into a Nonlinear Programming (NP) problem in order to apply the aforementioned projected saddle point dynamics.

5.1 Motivation and Modeling

The problem motivates from a physical communication network. Consider an undirected communication network $\mathcal{G(V,E,W)}$ whose nodes $i\in\mathcal{V}=\{1,2,\cdots,N\}$ are homogeneous base stations and can control their communication port gains $w_{k}^{(i)}\in\mathcal{W}$. (To abbreviate the notation, the edges are labelled with numbers.) The set of neighbours of node $i$ is denoted as $\mathcal{N}(i)$. The set of edges (communication channel) adjacent to node $i$ is denoted as $\mathcal{E}(i)$ and $\mathcal{E}=\bigcup_{i\in\mathcal{V}}\mathcal{E}(i)$. As illustrated in Fig. 1, the communication gain (strength) on each link $k\in\mathcal{E}$ is the sum of the port gains $w_{k}^{(i)}$ and $w_{k}^{(j)}$, $(i,j)=k\in\mathcal{E}$ contributed by the two end nodes connected by the edge. It is assumed that each agent can only get access to the information of its neighbours as well as the information of itself. Our goal is to develop a method so that each base station can adjust its own port gains only according to its neighbours’ information, the number of nodes $N$ and the information belonging to itself, so that the algebraic connectivity of the total communication network is maximized. The graph $\mathcal{G(V,E,W)}$ we consider is undirected, and hence the weighted Laplacian matrix $L_{w}$ is symmetric and can be expressed as

[TABLE]

and $L_{w}=\sum_{k\in\mathcal{E}}w_{k}E_{k}$, where $k$ is the label of the edges, and $0\leq w_{k}\in\mathcal{W},\>\forall k\in\mathcal{E}$ are the edge-weights. If node $i$ and $j$ are connected via edge $k$, then $[E_{k}]_{ii}=[E_{k}]_{jj}=1,\quad[E_{k}]_{ij}=[E_{k}]_{ji}=-1$, and the other elements of $E_{k}$ are zero. If the graph is connected, then the eigenvalues of $L_{w}$ satisfy: $0=\lambda_{1}<\lambda_{2}\leq\cdots\leq\lambda_{N}$ and $\lambda_{2}$ is the algebraic connectivity of $\mathcal{G(V,E,W)}$. Let $L=\sum_{k\in\mathcal{E}}E_{k}$ be the unweighted Laplacian matrix of the graph. We suppose $L$ has only one zero eigenvalue, namely, the unweighted graph $\mathcal{G(V,E)}$ is connected.

It is assumed that the total amount of port gain that each base station can provide is fixed. Without loss of generality, we assume $\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}=1$. Note that this differs from the formulation in (Göring et al., 2008), while they assume $\sum_{k\in\mathcal{E}}w_{k}=1$, and if written in the form of our model, $\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}=1$. In other words, it implies that the budget of the port gains in the entire network is a constant and the power can be allocated differently to each node. Hence each node is not homogeneous any more. Therefore, the algebraic connectivity maximization of an edge-weighted Laplacian matrix can be formulated as the following SDP problem:

[TABLE]

In (PC), the variable $\mu$ is used to shift the zero eigenvalue of $L_{w}$ with its eigenvector $\bm{1}$. When the optimal value is reached, $\lambda_{2}^{*}$ would be the smallest eigenvalue of $\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}}w_{k}^{(i)*}E_{k}+\mu^{*}\mathbbm{1}$. Since for any positive semi-definite matrix $G$, it holds that $\xi I\preceq G$, where $\xi$ is the smallest eigenvalue of $G$, we get the above constraints. Moreover, since $\lambda_{2}^{*}$ is the second smallest eigenvalue of $\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}}w_{k}^{(i)*}E_{k}$, it is a continuous function with respect to $\{w_{k}^{(i)*}\}$; and $\{w_{k}^{(i)*}\}$ lives in a compact set, hence the optimal value of (PC) can be attained. On the other hand, we can choose $\lambda_{2}$ small enough and $\mu$ large enough to make the first matrix inequality constraint strictly holds, hence strong duality holds for (PC).

5.2 Problem Equivalence

In order to solve (PC) distributedly, we consider the following problem

[TABLE]

where $Z^{(i)},\forall i\in\mathcal{V}$ are symmetric matrices. They can be written as $Z^{(i)}=\sum_{l=1}^{\frac{N(N+1)}{2}}z^{(i)}_{l}B_{l}$, where $z^{(i)}_{l}\in\mathbb{R}$ is the matrix entry and $B_{l}$ is the basis matrix for $\mathcal{S}^{N}$. Both are labeled by $l$. To be more precise, if $z_{l}^{(i)}$ is the entry located on $p$th row and $q$th column of $Z^{(i)}$, then $[B_{l}]_{pq}=[B_{l}]_{qp}=1$ and other entries of $B_{l}$ remain zeros. The purpose of the introduction of $Z^{(i)}$ is to derive the consensus condition (6f) in the KKT conditions. The next proposition describes the relationship between (PC) and (PD).

Proposition 4**.**

If $\{\lambda_{2}^{(i)*},\mu^{(i)*},\{w_{k}^{(i)*}\},Z^{(i)*}\}$, $i\in\mathcal{V}$ solves (PD), then $\lambda_{2}^{*}=\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}$, $\mu^{*}=\sum_{i\in\mathcal{V}}\mu^{(i)*}$, $\{w_{k}^{(i)*}\}$ solves (PC). On the other hand, if there exists an optimal solution for (PC), then there also exists an optimal solution for (PD).

{pf}

The KKT conditions of (PC) for all $i\in\mathcal{V}$ are:

[TABLE]

while the KKT conditions of (PD) for all $i\in\mathcal{V}$ reads

[TABLE]

where $\Phi^{*}$, $\varphi_{k}^{(i)*}$ and $v^{(i)*}$ are the Lagrange multipliers corresponding to the matrix inequality constraint, inequality constraints and the equality constraint of (PC) respectively. Similarly, $\Phi^{(i)*}$, $\varphi_{k}^{(i)*}$ and $v^{(i)*}$ are the Lagrange multipliers of (PD).

Since (PD) is convex, the KKT conditions (6) becomes necessary and sufficient conditions for optimality. Hence $\{\lambda_{2}^{(i)*},\mu^{(i)*},\{w_{k}^{(i)*}\},Z^{(i)*}\}$, $i\in\mathcal{V}$ solves (PD), if and only if there exist Lagrange multipliers $\{\Phi^{(i)*}\}$, $\{v^{(i)*}\}$ and $\{\varphi^{(i)*}_{k}\}$ such that (6) holds. Meanwhile, recall that $\{B_{l}\}$ are basis matrices for $\mathcal{S}^{N}$ and since $\Phi^{(i)*}\in\mathcal{S}^{N},\forall i\in\mathcal{V}$, $\Phi^{(i)*}$ can be written as $\sum_{l}\phi_{l}^{(i)*}B_{l}$. Denote $\phi_{l}^{*}=[\phi_{l}^{(1)*},\cdots,\phi_{l}^{(N)*}]^{T}$ and (6f) can be written as $L\phi_{l}^{*}=0,\forall l$, where $L=\sum_{k\in\mathcal{E}}E_{k}$ is the unweighted Laplacian matrix. Since we assume that the graph is connected, then $\phi_{l}^{*}\in ker(L)=span(\bm{1}),\forall l$ and hence implies $\Phi^{(i)*}=\Phi^{(j)*}$ for all $i,j\in\mathcal{V}$. This means that (5a)-(5b) are the same as (6a)-(6b). Further, by adding (6d) and (6e) for each node $i\in\mathcal{V}$, the terms $tr[(Z^{(i)}-Z^{(j)})\Phi^{(i)}]$ and $Z^{(i)*}-Z^{(j)*}$ are cancelled. Denote $\lambda_{2}^{*}=\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}$, $\mu^{*}=\sum_{i\in\mathcal{V}}\mu^{(i)*}$, we get (5). Since (PC) is convex, then the KKT conditions are necessary and sufficient conditions for optimality. Hence the first part of the statement follows.

Now we show the second part of the statement. Suppose $\lambda_{2}^{*},\mu^{*},\{w_{k}^{(i)*}\}$ is an optimal solution to (PC), then there must exist Lagrange multipliers $\Phi^{*}$, $\{v^{(i)*}\}$ and $\{\varphi^{(i)*}_{k}\}$ such that (5) holds. Now choose $\{\hat{\lambda}_{2}^{(i)*},\hat{\mu}^{(i)*},\{\hat{w}_{k}^{(i)*}\},\hat{Z}^{(i)*}\}$ and Lagrange multipliers $\{\hat{\Phi}^{(i)*},\hat{v}^{(i)*},\{\hat{\varphi}_{k}^{(i)*}\}\}$ such that $\sum_{i\in\mathcal{V}}\hat{\lambda}_{2}^{(i)*}=\lambda_{2}^{*},\sum_{i\in\mathcal{V}}\hat{\mu}^{(i)*}=\mu^{*},\hat{w}_{k}^{(i)*}=w_{k}^{(i)*}$ and $\hat{\Phi}^{(i)*}=\Phi^{*},\quad\hat{v}^{(i)*}=v^{(i)*},\quad\hat{\varphi}_{k}^{(i)*}=\varphi_{k}^{(i)*},\forall k\in\mathcal{E}(i),\forall i\in\mathcal{V}$. The KKT conditions (6a)-(6c) and (6f) is trivially satisfied by the above construction. What remains to show is that there exists such $\{\hat{Z}^{(i)*}\}$ so that (6d) and (6e) are satisfied.

We first show there exists $\{\hat{Z}^{(i)*}\}$ such that (6e) is satisfied. Denote $\hat{A}^{(i)}=\hat{\lambda}_{2}^{(i)*}I-\hat{\mu}^{(i)*}\mathbbm{1}-\sum_{k\in\mathcal{E}(i)}\hat{w}_{k}^{(i)*}E_{k}$. Since $\lambda_{2}^{*},\mu^{*},\{w_{k}^{(i)*}\}$ satisfy (5e) and $\lambda_{2}^{*}=\sum_{i\in\mathcal{V}}\hat{\lambda}_{2}^{(i)*}$, $\mu^{*}=\sum_{i\in\mathcal{V}}\hat{\mu}^{(i)*}$, $\hat{w}_{k}^{(i)*}=w_{k}^{(i)*}$, we know that $\sum_{i\in\mathcal{V}}\hat{A}^{(i)}\preceq 0$. By choosing $P^{(i)}=-\hat{A}^{(i)}$ for $i=1,\cdots,N-1$ and $P^{(N)}=\sum_{i=1}^{N-1}\hat{A}^{(i)}$, we have $\hat{A}^{(i)}+P^{(i)}=0$ for $i=1,\cdots N-1$ and $\hat{A}^{(N)}+P^{(N)}=\sum_{i=1}^{N}\hat{A}^{(i)}\preceq 0$. Therefore, we have $\hat{A}^{(i)}+P^{(i)}\preceq 0$ for all $i\in\mathcal{V}$ and $\sum_{i\in\mathcal{V}}P^{(i)}=0$. What remains to show is that there exists $\{\hat{Z}^{(i)*}\}$, such that $P^{(i)}=\sum_{j\in\mathcal{N}(i)}(\hat{Z}^{(i)*}-\hat{Z}^{(j)*})$ for all $i\in\mathcal{V}$. Recall that $P^{(i)}=\sum_{l}p^{(i)}_{l}B_{l}$, $\hat{Z}^{(i)*}=\sum_{l}\hat{z}^{(i)*}_{l}B_{l}$ and denote $p_{l}=[p_{l}^{(1)},\cdots,p_{l}^{(N)}]^{T}$, $\hat{z}_{l}^{*}=[\hat{z}_{l}^{(1)*},\cdots,\hat{z}_{l}^{(N)*}]^{T}$. Since $\sum_{i\in\mathcal{V}}P^{(i)}=0$, it follows that $\bm{1}^{T}p_{l}=0$ for all $l$. Therefore, $p_{l}\in ker(\bm{1}^{T})=Im(L)$, where $L=\sum_{k\in\mathcal{E}}E_{k}$. This implies $p_{l}$ can be expressed as $L\hat{z}_{l}^{*}$ for some $\hat{z}_{l}^{*}$, namely, there exists $\{\hat{Z^{(i)}}\}$ such that $P^{(i)}=\sum_{j\in\mathcal{N}(j)}(\hat{Z}^{(i)*}-\hat{Z}^{(j)*})$ for all $i\in\mathcal{V}$ and hence (6e) is satisfied.

Now we show that $\{\hat{\lambda}_{2}^{(i)*},\hat{\mu}^{(i)*},\hat{Z}^{(i)*}\}$ and $\{\hat{\Phi}^{(i)*}\}$ chosen above satisfy (6d). Since $\hat{\Phi}^{(i)*}=\Phi^{*}$, $\hat{w}_{k}^{(i)*}=w_{k}^{(i)*}$ and $P^{(i)}=\sum_{j\in\mathcal{N}(i)}(\hat{Z}^{(i)*}-\hat{Z}^{(j)*})$, the left hand side of (6d) can be written as $\hat{\lambda}_{2}^{(i)*}-\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}tr(E_{k}\Phi^{(i)*})+tr(P^{(i)}\Phi^{*})$. Recall that $P^{(i)}=-\hat{A}^{(i)}=-\hat{\lambda}_{2}^{(i)*}I+\hat{\mu}^{(i)*}\mathbbm{1}+\sum_{k\in\mathcal{E}(i)}\hat{w}_{k}^{(i)*}E_{k},i=1,\cdots N-1$. In view of (5a), for $i=1,\cdots,N-1$, we have $\hat{\lambda}_{2}^{(i)*}-\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}tr(E_{k}\Phi^{*})-\hat{\lambda}_{2}^{(i)*}tr(\Phi^{*})+\hat{\mu}^{(i)*}tr(\mathbbm{1}\Phi^{*})+\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}tr(E_{k}\Phi^{*})=0$.

For $i=N$, $P^{(N)}=\sum_{i=1}^{N-1}\hat{A}^{(i)}=\sum_{i=1}^{N-1}\hat{\lambda}_{2}^{(i)*}I-\sum_{i=1}^{N-1}\hat{\mu}^{(i)*}\mathbbm{1}-\sum_{i=1}^{N-1}\sum_{k\in\mathcal{E}(i)}\hat{w}_{k}^{(i)*}E_{k}$. In view of (5a), (5d) and since $\sum_{i\in\mathcal{V}}\hat{\lambda}_{2}^{(i)*}=\lambda_{2}^{*},\sum_{i\in\mathcal{V}}\hat{\mu}^{(i)*}=\mu^{*}$, we have $\hat{\lambda}_{2}^{(N)}-\sum_{k\in\mathcal{E}(N)}w_{k}^{(N)*}tr(E_{k}\Phi^{*})+\sum_{i=1}^{N-1}\hat{\lambda}_{2}^{(i)*}tr(\Phi^{*})-\sum_{i=1}^{N-1}\hat{\mu}^{(i)*}tr(\mathbbm{1}\Phi^{*})-\sum_{i=1}^{N-1}\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}tr(E_{k}\Phi^{*})=0$. Therefore, by the variables construction above, the second part of the statement follows.∎

5.3 Relaxing SDP into NP

(PD) can be solved distributedly by using a similar method as (Zhang & Hu, 2016) when the graph is regular. However, here we would like to consider general graphs, not only regular ones. In order to apply the projected saddle point dynamics to solve (PD), the problem needs to be relaxed into an NP first. This is because (PD) is still an SDP problem, and its inequality matrices constraints would lead to positive semidefinite matrix Lagrangian multipliers $\Phi^{(i)}$. This makes it hard to apply the saddle point dynamics in (Cherukuri et al., 2016) to this problem since by the definition of the projection operator $\Pi_{K}$, it is clear that $\Pi_{K}:\mathbb{R}^{n}\times\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$. The projected saddle point dynamics is not defined on the cone of positive semidefinite matrices.

Now we introduce the convex function proposed by Nesterov (2007) which can be used to approximate the largest eigenvalue of a symmetric matrix. Given $X\in\mathcal{S}^{N}$, function $f_{\varepsilon}:\mathcal{S}^{N}\mapsto\mathbb{R}$ and reads $f_{\varepsilon}(X)=\varepsilon\ln tr(e^{X/\varepsilon})=\varepsilon\ln[\sum_{i=1}^{N}e^{\lambda_{i}(X)/\varepsilon}]$ and its derivative with respect to $X$ reads

[TABLE]

where $(\lambda_{i}(X),u_{i})$ are eigen-pairs of $X$ with $\|u_{i}\|=1,\forall i$. It has been proved in (Nesterov, 2007) that

[TABLE]

Hence when $\varepsilon$ is sufficiently small, $f_{\varepsilon}(X)\approx\lambda_{\max}(X)$.

Consider the following NP problem:

[TABLE]

where $X^{(i)}=-\mu^{(i)}\mathbbm{1}-\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}E_{k}+\sum_{j\in\mathcal{N}(i)}(Z^{(i)}-Z^{(j)})$ to abbreviate the notation. It is clear that (NP) is a convex problem, the Slater’s condition also holds for (NP) and hence strong duality holds. Therefore, KKT conditions becomes necessary and sufficient conditions for (NP). The next proposition shows that how well the approximation could be.

Proposition 5**.**

Suppose $\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}$ is the optimal objective function value of (PD), then

[TABLE]

where $X^{(i)*}=-\mu^{(i)*}\mathbbm{1}-\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}E_{k}+\sum_{j\in\mathcal{N}(i)}(Z^{(i)*}-Z^{(j)*})$. Moreover, suppose $\{\hat{\mu}^{(i)*},\{\hat{w}_{k}^{(i)*}\},\hat{Z}^{(i)*}\}$ is the optimal solution to (NP), then

[TABLE]

where $\hat{X}^{(i)*}=-\hat{\mu}^{(i)*}\mathbbm{1}-\sum_{k\in\mathcal{E}(i)}\hat{w}_{k}^{(i)*}E_{k}+\sum_{j\in\mathcal{N}(i)}(\hat{Z}^{(i)*}-\hat{Z}^{(j)*})$.

{pf}

From the KKT condition (6e) of (PD), we know that $\lambda_{\max}(X^{(i)*})\leq-\lambda_{2}^{(i)*}$. This implies

[TABLE]

On the other hand, by Proposition 4, $\sum_{i\in\mathcal{V}}X^{(i)*}=-\sum_{i\in\mathcal{V}}\mu^{(i)*}\mathbbm{1}-\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}E_{k}=-\mu^{*}\mathbbm{1}-\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*}E_{k}$, where $\lambda_{2}^{*}$, $\mu^{*}$ and $\{w_{k}^{(i)*}\}$ is the optimal solution to (PC). We know from (PC) that $-\lambda_{2}^{*}=\lambda_{\max}(-\mu^{*}\mathbbm{1}-\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)*})=\lambda_{\max}(\sum_{i\in\mathcal{V}}X^{(i)*})$ (Boyd & Vandenberghe, 2004). Hence by Proposition 4, we have

[TABLE]

Then it follows from (11) that $\sum_{i\in\mathcal{V}}\lambda_{\max}(X^{(i)*})\leq\lambda_{\max}(\sum_{i\in\mathcal{V}}X^{(i)*})$. On the other hand, by eigenvalue inequality, we know that $\sum_{i\in\mathcal{V}}\lambda_{\max}(X^{(i)*})\geq\lambda_{\max}(\sum_{i\in\mathcal{V}}X^{(i)*})$, and hence

[TABLE]

Since $\lambda_{\max}(X^{(i)*})\leq f_{\varepsilon_{i}}(X^{(i)*})\leq\lambda_{\max}(X^{(i)*})+\varepsilon_{i}\ln N$, we have $\sum_{i\in\mathcal{V}}\lambda_{\max}(X^{(i)*})\leq\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(X^{(i)*})\leq\sum_{i\in\mathcal{V}}\lambda_{\max}(X^{(i)*})+\sum_{i\in\mathcal{V}}\varepsilon_{i}\ln N$. Hence it follows from (12) and (13) that (9) holds.

On the other hand, by (8), it follows that $\lambda_{\max}(\hat{X}^{(i)*})\leq f_{\varepsilon_{i}}(\hat{X}^{(i)*})$ and hence $\sum_{i\in\mathcal{V}}\lambda_{\max}(\hat{X}^{(i)*})\leq\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(\hat{X}^{(i)*})$. Note that $\{-\lambda_{\max}(\hat{X}^{(i)*}),\hat{\mu}^{(i)*},\{\hat{w}_{k}^{(i)*}\},\hat{Z}^{(i)*}\}$ is also a feasible solution to (PD), then this implies $-\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}\leq\sum_{i\in\mathcal{V}}\lambda_{\max}(\hat{X}^{(i)*})$. In addition, since $\{\hat{\mu}^{(i)*},\{\hat{w}_{k}^{(i)*}\},\hat{Z}^{(i)*}\}$ is an optimal solution to (NP), it follows that $\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(\hat{X}^{(i)*})\leq\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(X^{(i)*})$. Moreover, using (9), we have $-\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}\leq\sum_{i\in\mathcal{V}}\lambda_{\max}(\hat{X}^{(i)*})\leq\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(\hat{X}^{(i)*})\leq\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(X^{(i)*})\leq-\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}+\sum_{i\in\mathcal{V}}\varepsilon_{i}\ln N$, which proves the statement. ∎ Proposition 5 shows that one can have a good approximation on $\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)}$ by choosing $\varepsilon_{i}$, $\forall i\in\mathcal{V}$ sufficiently small. Without losing generality, we choose $\varepsilon_{i}=\varepsilon$, $\forall i\in\mathcal{V}$. Note that (NP) explains the reason of the introduction of $\mu$ in (PC) instead of writing the constraint as $\lambda_{2}(I-\frac{1}{N}\mathbbm{1})\preceq\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}E_{k}$ as in (Ghosh & Boyd, 2006). The $\lambda_{2}I$ term is needed for the relaxation. It is worth noticing that $\mu^{*}$ does not necessarily equals to $\lambda_{2}^{*}/N$ in (PC). In fact, any (\lambda_{2}^{*},\mu^{*},\big{\{}w_{k}^{(i)*}\big{\}}) such that $\mu^{*}\geq\lambda_{2}^{*}/N$ is an optimal solution for (PC). Also note that $f_{\varepsilon}(X)$ is not a strictly convex function though it is convex. Indeed, for all $0\leq\alpha\leq 1$, it holds that $\alpha f_{\varepsilon}(I)+(1-\alpha)f_{\varepsilon}(2I)=\alpha\varepsilon\ln(Ne^{\frac{1}{\varepsilon}})+(1-\alpha)\varepsilon\ln(Ne^{\frac{2}{\varepsilon}})=\varepsilon\ln N+2-\alpha=\varepsilon\ln(Ne^{\frac{2-\alpha}{\varepsilon}})=f_{\varepsilon}((2-\alpha)I)=f_{\varepsilon}(\alpha I+(1-\alpha)2I)$. Hence $f_{\varepsilon}(X)$ is not strictly convex.

5.4 Projected Dynamics and Numerical Examples

Now we apply the projected dynmaical system to solve (NP). To abbreviate the notation, denote x=\big{[}x^{(1)T},\cdots,x^{(N)T}\big{]}^{T}, where x^{(i)}=\big{[}\mu^{(i)},\{w^{(i)}_{k}\},\{z^{(i)}_{l}\}\big{]}^{T} and v=\big{[}v^{(1)},\cdots,v^{(N)}\big{]}^{T}.

The projected dynamics for each agent $i$ is given by

[TABLE]

[TABLE]

where $\mathscr{L}_{\mathcal{A}}(x,v)=\sum_{i\in\mathcal{V}}\{f_{\varepsilon}(X^{(i)})+v^{(i)}(\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}-1)+\frac{1}{2}(\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}-1)^{2}\}$ and $\nabla_{X^{(i)}}f_{\varepsilon}(X^{(i)})$ is given by (7). Note that in (14), (15), every agent only uses the information that belongs to its neighbours as well as to itself. The exchanging information for each agent is the gradient $\nabla_{X^{(i)}}f_{\varepsilon}(X^{(i)})$. We would like to remark that, although each time step each agent has to communicate with its neighbour a vector of size $\frac{N(N+1)}{2}$, this is the price to pay in order to solve the problem distributedly. This is because of the “dense” structure of the problem (since we do not make special assumptions on the graph topology) and the constraint that the communication network is the physical network itself. The reasons above have made it hard to decompose the problem into small scales.

Theorem 6**.**

The system (14), (15) is well-defined and the trajectory asymptotically converges to one of the saddle points of (NP) for all initial values $\mu^{(i)}(0),z_{l}^{(i)}(0),v^{(i)}(0)$ $\in\mathbb{R}$, $w_{k}^{(i)}(0)\in\mathbb{R}_{+}$.

{pf}

$f_{\varepsilon}(X)$ has a Lipschitz continuous gradient with respect to $x$ given that $X=\sum_{i}x_{i}A_{i}$, where all $A_{i}$ are symmetric matrices (Nesterov, 2007). By Theorem 2.5 in (Nagurney & Zhang, 2012), for any initial value $\mu^{(i)}(0)\in\mathbb{R}$, $w_{k}^{(i)}(0)\in\mathbb{R}_{+}$, $v^{(i)}(0)\in\mathbb{R}$ and $z_{l}^{(i)}(0)\in\mathbb{R}$, there exists a unique Carathéodory solution which continuously depends on the initial value. Therefore the system (14), (15) is well-defined and by Proposition 3, the system (14), (15) asymptotically converges to one of the saddle points of (NP).∎ It seems that when simulating the projected dynamics (14),(15), one has to do eigenvalue decomposition on $X^{(i)}$ to compute $\nabla_{X^{(i)}}f(X^{(i)})$ at each time step. However, since the factors $e^{\lambda_{i}(X^{(i)})/\varepsilon}$ decrease very rapidly, the gradient numerically only depends on few largest eigenvalues and correpondant eigenvectors (Nesterov, 2007). Extreme eigenvalues will converge first in numerical methods such as Arnoldi scheme, hence one does not have to do the entire eigenvalue decomposition and the numerical complexity is reduced.

By Proposition 5, one can first choose an $\varepsilon$ and get an “optimal” algebraic connectivity under the current choice of $\varepsilon$. If the approximation error $\varepsilon N\ln N$ compared to the “optimal” algebraic connectivity under the current choice of $\varepsilon$ is not satisfying (for example, $\varepsilon N\ln N$ is approximately 10% of the current “optimal” algebraic connectivity), one can decrease $\varepsilon$ until the desired relative error is achieved.

Example 7**.**

We run a simple numerical example using the graph illustrated in Fig. 1 to show that the variables do converge to the optimal solution. Using CVX, we get the optimal solution of (PC): $w_{1}^{(1)*}=w_{2}^{(3)*}=1$, $w_{1}^{(2)*}=w_{2}^{(2)*}=0.5$ and $\lambda_{2}^{*}=1.5$. Forward Euler method is used to discretize (14), (15) and we choose $\varepsilon=0.01$, time step size $\Delta t=0.01$.

Fig. 2 shows that the edge weights converge to the optimal solution. ∎

Example 8**.**

Now we consider a more complicated graph generated by ten nodes. Forward Euler method is also used to discretize (14), (15) and we choose different $\varepsilon$ and time step size $\Delta t$ to illustrate the effect of $\varepsilon$ on $\Delta t$. According to (Nesterov, 2007), the choice of $\varepsilon$ affects the Lipschitz constant of $\nabla_{X}f_{\varepsilon}(X)$ as well as the Hessian of $f_{\varepsilon}(X)$. The smaller $\varepsilon$ is, the bigger the Lipschitz constant of $\nabla_{X}f_{\varepsilon}(X)$ will be. Hence intuitively, bigger Lipschitz constant of the gradient implies a smaller step size to avoid the case of moving around in the neighbourhood of optimum without converging. Using CVX, we know the optimal value of (PC) is 1.141. We do the simulation for $t\in[0,50]$ with different $\varepsilon$ and time step sizes. In the end, we get $\lambda_{2}(L_{w})$ equals to 1.085, 1.091 and 1.128 when choosing $\varepsilon=10^{-2},\Delta t=10^{-3}$, $\varepsilon=10^{-3},\Delta t=10^{-3}$ and $\varepsilon=10^{-3},\Delta t=10^{-4}$, respectively. As illustrated in Fig. 3, the algebraic connectivity in the network does not converge to the optimal value of the unrelaxed and “centralized” problem $(P_{c})$. However, as we decrease $\varepsilon$, the limiting algebraic connectivity gets closer to the optimal value of $(P_{c})$. This illustrates the relaxation effect. In addition, the evolution of $\lambda_{2}(L_{w})$ involves a lot of oscillations when $\varepsilon=10^{-3}$ and $\Delta t=10^{-3}$, while it behaves much nicer when $\Delta t=10^{-4}$. This emperically shows that smaller $\varepsilon$ requires smaller time step length.∎

Fig. 3 illustrates that a smaller $\varepsilon$ requires a smaller step size when discretizing system (14), (15). This means that when the number of nodes $N$ goes large, in order to get a good approximation of (PD), we need a very small $\varepsilon$ and hence it leads to a very small step size. This would result in the slow evolution of the system states per iteration and hence requires a large number of iterations to reach the equilibrium.

One practical solution to the issue above is presented as follows. We can solve the problem above by modifying $\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}=1$ in (PD) as $\sum_{k\in\mathcal{E}(i)}w_{k}^{(i)}=a$, where $a>0$. We call the modified optimization problem and its relaxed nonlinear programming problem as $(P_{M})$ and $(NP_{M})$ respectively. By checking the optimality conditions (6), we conclude that $\{\lambda_{2}^{(i)*},\mu^{(i)*},\{w_{k}^{(i)*}\},Z^{(i)*}\}$ is the optimal solution to (PD) iff $\{a\lambda_{2}^{(i)*},a\mu^{(i)*},a\{w_{k}^{(i)*}\},aZ^{(i)*}\}$ is the optimal solution to $(P_{M})$ (since all the optimality conditions are linear). Using Proposition 5, we can conclude that $-\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}\leq\frac{1}{a}\sum_{i\in\mathcal{V}}\lambda_{\max}(\hat{X}^{(i)*})\leq\frac{1}{a}\sum_{i\in\mathcal{V}}f_{\varepsilon_{i}}(\hat{X}^{(i)*})\leq-\sum_{i\in\mathcal{V}}\lambda_{2}^{(i)*}+\frac{1}{a}\varepsilon N\ln N$ provided that $\{\hat{\mu}^{(i)*},\{\hat{w}_{k}^{(i)*}\},\hat{Z}^{(i)*}\}$ is the optimal solution to $(NP_{M})$. Therefore, apart from choosing $\varepsilon$ to be small, we can choose $a$ sufficiently large, solve $(NP_{M})$ and divide the optimal weight realization obtained from $(NP_{M})$ by $a$ to suppress the approximation error. Namely, we do not need to choose $\varepsilon$ to be too small so that the time step size does not need to be too small. Therefore, the number of iterations needed to reach the equilibrium is suppressed when $N$ goes large.

Example 9**.**

Consider a graph with 30 nodes. If we do not use the methodology above, namely, $a=1$, $\varepsilon$ needs to be at least $2.1416\times 10^{-5}$ so that the relative approximation error of the optimal algebraic connectivity $\varepsilon N\ln N/\lambda_{2}(L_{w}^{*})$ is within 5%. For comparison, if we fix $\varepsilon=0.5$ first, and choose $a$ such that the relative approximation error is within 5%. Multiple time steps have been tried and the largest ones such that the discretized systems converge are illustrated in Fig. 4. Same initial values and forward Euler discretization are used. The algebraic connectivity that uses the methodology mentioned above uses much fewer iterations to converge to the optimal value.

Example 10**.**

Consider the same graph used in Example 8. We use forward Euler for discretization. Same initial values, time step sizes and $\varepsilon$ are used ($\varepsilon=0.5$, $\Delta t=0.02$). Fig. 5 shows how $a$ affects the approximation error.

The complexity per iteration for agent $i$ is $\mathcal{O}(|\mathcal{E}(i)|\cdot N^{2})$. Since the right hand side of (14), (15) involves only with special matrices such as $\mathbbm{1}$, $E_{k}$ and $B_{l}$, one does not need to do matrix multiplication and hence the complexity is greatly reduced.

Example 11**.**

We test the algorithm on larger scale networks and plot the ratio between running time and $|\mathcal{E}(i)|\cdot N^{2}$ versus $N$. Forward Euler is used for discretization. To eliminate the influence of the network topologies and number of edges on the convergence, we choose the same families of graphs and let $N$ varies. We consider the family of ring graphs and complete graphs. $\varepsilon=0.5,\Delta=0.01$ and $a$ is chosen such that the relative approximation error of the objective function is within 5%. The iterations are terminated when the infinity norm of the right hand side of (14), (15) is smaller than $10^{-3}$. The result is shown as Fig. 6. It can be seen that the ratio between running time and $|\mathcal{E}(i)|\cdot N^{2}$ is approximately constant when $N$ changes.

6 Conclusion

In this paper, a projected saddle point dynamics of augmented Lagrangian is presented to solve, not necessary strictly, convex optimization problems. As a supplement to the analysis in (Niederländer & Cortés, 2016), we show that the projected saddle point dynamics converges to one of the saddle points. Moreover, the problem of distributedly maximizing the algebraic connectivity of an undirected communication network by optimizing the port gains of each nodes (base stations) is considered. The original SDP problem is relaxed into an NP problem and then the aforementioned projected dynamical system is applied to solve the NP. Numerical examples are used to illustrate: 1. the convergence of the edge weights to one of the optimal solutions; 2. the effect of $\varepsilon$ on the choice of time step size; 3. complexity per iteration of the algorithm. A methodology is presented so that the number of iterations needed to reach the equilibrium is suppressed.