Projected Primal-Dual Gradient Flow of Augmented Lagrangian with Application to Distributed Maximization of the Algebraic Connectivity of a Network
Han Zhang Jieqiang Wei Peng Yi Xiaoming Hu



   

  [footnoteinfo]This work is supported by China Scholarship Council.

  KTH_M]Han Zhanghanzhang@kth.se,
  KTH_A]Jieqiang Weijieqiang@kth.se,
  Toronto]Peng Yipeng.yi@utoronto.ca,
  KTH_M]Xiaoming Huhu@kth.se,              

  [KTH_M]Department of Mathematics, KTH Royal Institute of Technology, SE-100 44, Stockholm, Sweden     
  [KTH_A]Department of Automatic Control, KTH Royal Institute of Technology, SE-100 44, Stockholm, Sweden      
  [Toronto]Department of Electrical and Systems Engineering, Washington University in St. Louis, USA                                    


  
	Projected Dynamical Systems, semi-definite programming, distributed optimization              
                        

  
	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 <cit.>, 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.
  





§ INTRODUCTION

Recently, due to the fact of rising interest in distributed optimization, there has been a number of advances in the field of optimization that involves using continuous time dynamical systems to solve smooth convex problems, such as <cit.>, <cit.>, <cit.>, etc, as well as nonsmooth convex problems <cit.>.
Though it seems to be quite different from the classical iterative optimization methods, it has been found that classical methods such as steepest descent and proximal point method are nothing but forward Euler and backward Euler discretization of the same saddle point dynamics <cit.>. On the other hand, continuous-time algorithms enable many powerful mathematical tools to analyze the convergence of the algorithms. Hence studying continuous-time algorithms provides us with more insights about the theory behind some iterative algorithms as well as more ideas of designing them.


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 <cit.>, <cit.> and has been revisited by <cit.>.
<cit.> 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. <cit.> 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, <cit.> proposes a smooth vector field for seeking the saddle points of strictly convex problems. <cit.> 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, <cit.> 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 <cit.>. In a recent work <cit.> and its conference version <cit.>, 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 <cit.>, 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 <cit.> that the algebraic connectivity is monotonically increasing with the algorithm, while the convergence to one optimal solution is not explicitly given. <cit.>,<cit.> and <cit.> 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 <cit.> and its conference version <cit.>, 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, <cit.>).
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.



§ PRELIMINARIES AND NOTATIONS


We denote 1=11^T as an N dimensional all-one matrix, where 1 is an N dimensional all-one vector. The element located on the ith row and and jth 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≽ A_2. We use · 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(·) denotes the trace of a matrix. ⟨·,·⟩_2 is denoted as the inner-product in Euclidean space and ⟨ A_1,A_2⟩_M=tr(A_1A_2) denotes the inner-product in 𝒮^n, which is the Hilbert space of n× n symmetric matrix. 




Assume K⊂ℝ^n is a closed and convex set, the projection of a point x to the set K is defined as P_K(x)=min_y∈ Kx-y. For x∈ K, v∈ℝ^n, the projection of the vector v at x with respect to K is defined as: ( see <cit.>,<cit.>) Π_K(x,v)=lim_δ→ 0P_K(x+δ v)-x/δ=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), ∂ K and cl(K), respectively.
The set of inward normals of K at x is defined as n(x)={γ | γ=1,⟨γ,x-y⟩_2≤ 0,∀ y∈ K }, and Π_K(x,v) fulfills the following lemma:

  If x∈ int(K), then Π_K(x,v)=v; if x∈∂ K, then Π_K(x,v)=v+β(x)n^*(x), where n^*(x)=max_n∈ n(x)⟨ v,-n⟩ and β(x)=max{0,⟨ v,-n^*(x)⟩}.





Let F be a vector field such that F:K↦ℝ^n, the projected dynamical system is given by ẋ=Π_K(x,F(x)).
Note that the right hand side of above dynamics can be discontinuous on the ∂ K. Hence given an initial value x_0∈ 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 <cit.>.
Now we introduce some basic knowledge about SDP. An SDP problem in standard form can be expressed as

    maximize ∑_i=1^m c_ix_i
       subjected to ∑_i=1^mx_iA_i≼ S,

where A_i,S∈𝒮^N.
Introducing the Lagrangian multiplier Φ≽ 0, the Lagrangian function of (<ref>) can be written as

    ℒ(x,Φ)=-∑_i=1^m c_ix_i+tr(Φ(∑_i=1^m x_iA_i-S))

Taking the gradient of the Lagrangian function, together with the constraint and the complementarity slack condition, we get its KKT conditions

    ⟨ A_i,Φ⟩_M=-c_i,  ∑_i=1^mx_iA_i≼ S,
    Φ≽ 0,  ⟨∑_i=1^mx_iA_i-S,Φ⟩_M=0.




§ PROBLEM FORMULATION AND PROJECTED SADDLE POINT DYNAMICS

In this section, we consider the following optimization problem defined on ℝ^n:

    x∈ Kminimize       f(x)
       subject to       Ax-b=0,

where f:ℝ^n↦ℝ and A∈ℝ^m× 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 (<ref>). Hence strong duality holds for (<ref>).




The Lagrangian ℒ:K×ℝ^m↦ℝ for the problem (<ref>) is given by

    ℒ(x,v)=f(x)+v^T(Ax-b),
 
where v∈ℝ^m is the Lagrangian multiplier of the constraint Ax-b=0. Since strong duality holds for (<ref>), then (x^*,v^*) is a saddle point of ℒ(x,v) if and only if x^* is an optimal solution to (<ref>) and v^* is optimal solution to its dual problem. 
The augmented Lagrangian ℒ_𝒜:K×ℝ^m↦ℝ for (<ref>) is given by ℒ_𝒜(x,v)=f(x)+v^T(Ax-b)+ρ/2(Ax-b)^T(Ax-b), where ρ>0 is the damping parameter that will help to suppress the oscillation of x during optimization algorithms. Without loss of generality, we choose ρ=1.




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

  
    ẋ   =Π_K(x,-∇ f(x)-A^Tv-A^T(Ax-b))
       =Π_K(x,-∂ℒ_𝒜(x,v)/∂ x),
    v̇   =Ax-b=∂ℒ_𝒜(x,v)/∂ v.


Note that it is assumed that ∇ f(x) is locally Lipschitz continuous, therefore there is a unique Carathéodory solution for the dynamics (<ref>).

  Indeed, one can choose to project the negative gradient flow of the objective function onto the entire feasible set. However, this increases the computational complexity since it adds more constraints while computing the projection. Recall that we assume that calculating the projection to T_K(x) is computationally cheap, primal-dual gradient flow is used. 


§ CONVERGENCE ANALYSIS

In this section, we analyse the convergence for (<ref>) and start with the analysis of the equilibrium point of (<ref>). <cit.> 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.

  (x^*,v^*) is a saddle point to (<ref>) if and only if it is an equilibrium of (<ref>).


  Since strong duality holds for (<ref>), the optimality conditions become necessary and sufficient conditions. The optimality condition for (<ref>) is given by -∇ f(x^*)-A^Tv^*∈ N_K(x^*),   Ax^*-b=0, <cit.>, which implies -∇ f(x^*)-A^Tv^*+A^T(Ax^*-b)∈ N_K(x^*), where N_K(x^*) denotes the normal cone of K at x^*.
  This implies Π_K(x^*,-∇ f(x^*)-A^Tv^*-A^T(Ax^*-b))=0, therefore, (x^*,v^*) is an equilibrium point of (<ref>). On the other hand, if (x^*,v^*) is an equilibrium point of (<ref>), it must have -∇ f(x^*)-A^Tv^*+A^T(Ax^*-b)∈ N_K(x^*) and Ax^*-b=0, which implies the optimality condition.
  


  Assume γ is a solution to ẋ=f(x), where f:ℝ^n↦ℝ^n. The omega-limit set is defined as
  
    Ω(γ)=   {y∈ℝ^n|∃{t_k}_k=1^∞⊂[0,∞), lim_k→∞t_k=∞,
       lim_k→∞γ(t_k)=y}

  
 
  
	(LaSalle invariance principle for Carathéodory solutions) Let 𝒮∈ℝ^n be compact and invariant. Assume for each x_0∈𝒮, there exists a unique solution for ẋ=f(x) starting at x_0 and its omega-limit set is invariant. Let V:ℝ^n↦ℝ be a continuously differentiable map such that ℒ_fV(x)≤ 0, where ℒ_f=∂ V/∂ xf(x) denotes the Lie derivative along the vector field f(x). Then any solution of ẋ=f(x) starting at 𝒮 converges to the largest invariant set in cl({x∈𝒮|ℒ_fV(x)=0}).
   
	Given an initial value (x(0),v(0)), where x(0)∈ K, the trajectory of the projected dynamical system (<ref>) asymptotically converges to one of the saddle points of (<ref>).
  
  
  We use LaSalle invariance principle for Carathéodory solutions <cit.> to prove the proposition.
	Suppose (x^*,v^*) is a saddle point of the Lagrangian (<ref>), namely, x^* is the optimal solution of (<ref>) and v^* is the optimal solution of its dual problem. Construct the following Lyapunov function
	
    d(x,v)=1/2(x-x^*^2+v-v^*^2).

	Note that d(x,v) is continuously differentiable and denote the right hand side of the dynamics (<ref>) as vector field F. The Lie derivative along the vector field f(x) of a function V(x) is defined as ℒ_fV(x)=∂ V(x)/∂ xf(x).
	By the definition of saddle points, ℒ(x^*,v)≤ℒ(x^*,v^*)≤ℒ(x,v^*). The Lie derivative of d(x,v) along the vector field F is given by ℒ_Fd(x,v)=(x-x^*)^TΠ_K(x,-∇ f(x)-A^Tv-A^T(Ax-b))+(v-v^*)^T(Ax-b).
	By Lemma <ref>, it holds that Π_K(x,-∇ f(x)-A^Tv-A^T(Ax-b))=-∇ f(x)-A^Tv-A^T(Ax-b)+β(x)n^*(x), where β(x)≥ 0 and ⟨ n^*(x),x-y⟩≤ 0, ∀ y∈ K. Since Ax^*-b=0, hence it follows that ℒ_Fd(x,v)=(x-x^*)^T(-∇ f(x)-A^Tv-A^T(Ax-b)+β(x)n^*(x))+(v-v^*)^T(Ax-b)≤ (x-x^*)^T(-∇ f(x)-A^Tv-A^T(Ax-b))+(v-v^*)^T(Ax-b)=-(x-x^*)^T∂ℒ(x,v)/∂ x-(Ax-b)^T(Ax-b)+(v-v^*)^T ∂ℒ(x,v)/∂ v.
	
	

	Since ℒ(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 <cit.> that ℒ_Fd(x,v)≤ℒ(x^*,v)-ℒ(x,v^*)-(Ax-b)^T(Ax-b).
	Since ℒ(x^*,v)≤ℒ(x,v^*), we have ℒ_Fd(x,v)≤ -(Ax-b)^T(Ax-b)≤ 0. Note that d(x,v) is convex and differentiable. By Proposition <ref> and using the result in <cit.>, 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 (<ref>) 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 d/dtd(x(t),v(t))=ℒ_Fd(x,v) holds almost everywhere on t≥ 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 𝒮̂={(x,v)∈ℝ^n×ℝ^m|d(x,v)≤ d(x(0),v(0))} is a compact invariant set for the system (<ref>). 
	The invariance of the ω-limit set <cit.> can be proved by using the same methodology as Lemma 4.1 in <cit.> which is based on the continuity and the uniqueness of the solution. 
	Now by LaSalle invariance principle for Carathéodory solutions <cit.>, the trajectory of (<ref>) converges to the largest invariant set in cl(Ω), where Ω={(x,v)∈𝒮̂|ℒ_Fd(x,v)=0}.



For (x^',v^')∈Ω, we have Ax^'-b=0 and ℒ(x^*,v^')=ℒ(x^',v^*)⇔ f(x^*)=f(x^')+v^*T(Ax^'-b)⇔ f(x^*)=f(x^').
Since x^* is an optimal solution to (<ref>), x^' is also an optimal solution to (<ref>).
	Denote 𝒟={(x,v)∈𝒮̂|x } and note that Ω⊂𝒟. Hence cl(Ω)⊂ cl(𝒟).
	On the other hand, since the set of optimal solutions for a convex optimization problem is closed, 𝒟 is also closed and hence cl(Ω)⊂ cl(𝒟)=𝒟. 




	Denote ℳ as the largest invariant set in cl(Ω). 
	Recall that (<ref>) is not necessarily a strictly convex problem; for the non-strictly convex case, since there might be infinite number of optimizers and they are not isolated, the trajectory of x(t) might still evolve in the optimal solution set. To show that the trajectory converges to a point, we need to show there that if a trajectory starts in ℳ, it will remain constant for all times; since the Carathéodory solution to the projected dynamics is unique and continuously depend on the initial value, it can be concluded from the above fact that the trajectory asymptotically converges to a point.
	Assume the initial value (x̂(0),v̅(0))∈ℳ⊂ cl(Ω), then x̂(0) is also an optimal solution to (<ref>). Hence there exists some v̂ such that (x̂(0),v̂) is a saddle point of (<ref>). Then the Lyapunov function can be constructed similarly as d̂(x,v)=1/2(x-x̂(0)^2+v-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 d̂(x(t),v(t)). 
	Furthermore, since (x̂(0),v̅(0))∈ℳ and ℳ is invariant, (x̂(t),v̅(t))∈ℳ,∀ t≥ 0. This implies that x̂(t) is an optimizer for all t≥ 0 and hence we have v̇̅̇(t)=Ax̂(t)-b=0,∀ t≥ 0. And this implies v̅(t)=v̅(0),∀ t≥ 0 and hence d̂(x̂(t),v̅(t))≤d̂(x̂(0),v̅(0))=1/2v̅(0)-v̂^2. This implies x̂(t)-x̂(0)^2=0 and hence x(t)=x̂(0),∀ t≥ 0. Therefore, any trajectory that starts in ℳ remains constant for all times, i.e., any point in ℳ is an equilibrium point. By Proposition <ref>, we can conclude that they are saddle points of (<ref>).
	This implies that the trajectory starts at (x̂(0),v̅(0)) can not be time-varying, namely, x̂(t) remains constant for all t>0. This is because d̂(x̂(0),v̅(0))=1/2v̅(0)-v̂^2 and d/dtv̅(0)=0, hence x̂(t) can not deviate from its initial value x̂(0), otherwise it would contradict the fact that d̂(x(t),v(t)) is non-increasing with respect to time. Therefore, any trajectory that starts in ℳ remains constant for all times, i.e., any point in ℳ is an equilibrium point. By Proposition <ref>, we can conclude that they are saddle points of (<ref>).
	


	We have shown that given an initial value (x(0),v(0)), where x(0)∈ K, the trajectory of (<ref>) asymptotically converges to a set ℳ whose elements are saddle points of (<ref>). Now we show the trajectory asymptotically converges to a point in ℳ by contradiction. To abbreviate the notation, we denote η=(x^T,v^T)^T.
	Suppose η(t) does not converge to a point in ℳ, namely, the trajectory's ω-limit set Γ(η(t)) is not a singleton. This means we can choose η̅_1,η̅_2∈Γ(η(t))⊂ℳ, such that η̅_1-η̅_2=ζ,ζ>0. Since η̅_1,η̅_2∈ℳ, η̅_1,η̅_2 are saddle points and we have shown that all saddle points are Lyapunov stable.
	This means that there exists δ(ζ/2), such that if η(T)-η̅_1<δ(ζ/2), then η(t)-η̅_1<ζ/2,∀ t≥ T. Since η̅_1 is an ω-limit point in Γ(η(t)), there exists such T so that η(T)-η̅_1<δ(ζ/2), and hence the trajectory can never leave ζ/2-neighbourhood of η̅_1 after time instant T. But η̅_2 is also an ω-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 ℳ.
	
	
		
		
	   


§ 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.
  

 §.§ Motivation and Modeling

  The problem motivates from a physical communication network. Consider an undirected communication network 𝒢(𝒱,ℰ,𝒲) whose nodes i∈𝒱={1,2,⋯,N} are homogeneous base stations and can control their communication port gains w_k^(i)∈𝒲. (To abbreviate the notation, the edges are labelled with numbers.) The set of neighbours of node i is denoted as 𝒩(i). The set of edges (communication channel) adjacent to node i is denoted as ℰ(i) and ℰ=⋃_i∈𝒱ℰ(i). 
  As illustrated in Fig. <ref>, the communication gain (strength) on each link k∈ℰ is the sum of the port gains w_k^(i) and w_k^(j), (i,j)=k∈ℰ 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 𝒢(𝒱,ℰ,𝒲) we consider is undirected, and hence the weighted Laplacian matrix L_w is symmetric and can be expressed as
  
    [L_w]_ij   =
    	  ∑_lw_il     i= j  (i,l)∈ℰ
    
    		-w_ij     i≠ j  (i,j)∈ℰ
    
    		0

  and L_w=∑_k∈ℰw_kE_k, where k is the label of the edges, and 0≤ w_k∈𝒲, ∀ k∈ℰ are the edge-weights. If node i and j are connected via edge k, then [E_k]_ii=[E_k]_jj=1,   [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=λ_1<λ_2≤⋯≤λ_N and λ_2 is the algebraic connectivity of 𝒢(𝒱,ℰ,𝒲). Let L=∑_k∈ℰE_k be the unweighted Laplacian matrix of the graph. We suppose L has only one zero eigenvalue, namely, the unweighted graph 𝒢(𝒱,ℰ) 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 ∑_k∈ℰ(i)w_k^(i)=1. Note that this differs from the formulation in <cit.>, while they assume ∑_k∈ℰ w_k=1, and if written in the form of our model, ∑_i∈𝒱∑_k∈ℰ(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:
  
    λ_2,μ,{w_k^(i)}maximize      λ_2
       subject to      λ_2I-μ1≼∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)E_k,
             ∑_k∈ℰ(i)w_k^(i)=1,
             w_k^(i)≥ 0,  ∀ k∈ℰ,   i∈𝒱.
    	PC





  In (<ref>), the variable μ is used to shift the zero eigenvalue of L_w with its eigenvector 1. When the optimal value is reached, λ_2^* would be the smallest eigenvalue of ∑_i∈𝒱∑_k∈ℰ w_k^(i)*E_k+μ^*1. Since for any positive semi-definite matrix G, it holds that ξ I≼ G, where ξ is the smallest eigenvalue of G, we get the above constraints. 
  Moreover, since λ_2^* is the second smallest eigenvalue of ∑_i∈𝒱∑_k∈ℰ 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 (<ref>) can be attained. On the other hand, we can choose λ_2 small enough and μ large enough to make the first matrix inequality constraint strictly holds, hence strong duality holds for (<ref>).



 §.§ Problem Equivalence

  In order to solve (<ref>) distributedly, we consider the following problem
  
    {λ_2^(i)},{μ^(i)},{w_k^(i)},{Z^(i)}maximize      ∑_i∈𝒱λ_2^(i)
       subject to      λ_2^(i)I-μ^(i)1- ∑_k∈ℰ(i)w_k^(i)E_k
             +∑_j∈𝒩(i)(Z^(i)-Z^(j))≼ 0,
             ∑_k∈ℰ(i)w_k^(i)=1,
             w_k^(i)≥ 0, ∀ k∈ℰ(i),  ∀ i∈𝒱,
    	PD

  where Z^(i),∀ i∈𝒱 are symmetric matrices. They can be written as Z^(i)=∑_l=1^N(N+1)/2z^(i)_lB_l, where z^(i)_l∈ℝ is the matrix entry and B_l is the basis matrix for 𝒮^N. Both are labeled by l. To be more precise, if z_l^(i) is the entry located on pth row and qth 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 (<ref>) in the KKT conditions.
  The next proposition describes the relationship between (<ref>) and (<ref>). 
  
  
  
	
		If {λ_2^(i)*, μ^(i)*,{w_k^(i)*}, Z^(i)*}, i∈𝒱 solves (<ref>), then λ_2^*=∑_i∈𝒱λ_2^(i)*, μ^*=∑_i∈𝒱μ^(i)*, {w_k^(i)*} solves (<ref>). 
		On the other hand, if there exists an optimal solution for (<ref>), then there also exists an optimal solution for (<ref>).
  
  
	The KKT conditions of (<ref>) for all i∈𝒱 are:
	
	  
    tr(Φ^*)=1, Φ^*≽ 0, φ^(i)*_k≥ 0, tr(1Φ^*)=0,
       tr(E_kΦ^*)-v^(i)*+φ_k^(i)*=0,  w_k^(i)*≥ 0,
       φ_k^(i)*w_k^(i)*=0, ∀ k∈ℰ(i),
       λ_2^*-∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)*tr(E_kΦ^*)=0,
       λ_2^*I-μ^*1- ∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)*E_k≼ 0,
       ∑_k∈ℰ(i)w_k^(i)*=1,

	
	while the KKT conditions of (<ref>) for all i∈𝒱 reads
	
	  
    tr(Φ^(i)*)=1,Φ^(i)*≽ 0,φ^(i)*_k≥ 0,tr(1Φ^(i)*)=0,
       tr(E_kΦ^(i)*)-v^(i)*+φ_k^(i)*=0,  w_k^(i)*≥ 0,
       φ_k^(i)*w_k^(i)*=0, ∀ k∈ℰ(i),
       λ_2^(i)*-∑_k∈ℰ(i)w_k^(i)*tr(E_kΦ^(i)*)
                 +∑_j∈𝒩(i)tr[(Z^(i)-Z^(j))Φ^(i)*]=0,
       λ_2^(i)*I-μ^(i)*1- ∑_k∈ℰ(i)w_k^(i)*E_k
                       +∑_j∈𝒩(i)(Z^(i)*-Z^(j)*)≼ 0,
       ∑_k∈ℰ(i)w_k^(i)*=1,  ∑_j∈𝒩(i)Φ^(i)*-Φ^(j)*=0,

	
	where Φ^*, φ_k^(i)* and v^(i)* are the Lagrange multipliers corresponding to the matrix inequality constraint, inequality constraints and the equality constraint of (<ref>) respectively. Similarly, Φ^(i)*, φ_k^(i)* and v^(i)* are the Lagrange multipliers of (<ref>).




	Since (<ref>) is convex, the KKT conditions (<ref>) becomes necessary and sufficient conditions for optimality. Hence {λ_2^(i)*, μ^(i)*,{w_k^(i)*}, Z^(i)*}, i∈𝒱 solves (<ref>), if and only if there exist Lagrange multipliers  {Φ^(i)*}, {v^(i)*} and {φ^(i)*_k} such that (<ref>) holds.
	
	Meanwhile, recall that {B_l} are basis matrices for 𝒮^N and since Φ^(i)*∈𝒮^N,∀ i∈𝒱, Φ^(i)* can be written as ∑_lϕ_l^(i)*B_l. Denote ϕ_l^*=[ϕ_l^(1)*,⋯,ϕ_l^(N)*]^T and (<ref>) can be written as Lϕ_l^*=0,∀ l, where L=∑_k∈ℰE_k is the unweighted Laplacian matrix. Since we assume that the graph is connected, then ϕ_l^*∈ ker(L)= span(1),∀ l and hence implies Φ^(i)*=Φ^(j)* for all i,j∈𝒱.
	This means that (<ref>)-(<ref>) are the same as (<ref>)-(<ref>). 
	Further, by adding (<ref>) and (<ref>) for each node i∈𝒱, the terms tr[(Z^(i)-Z^(j))Φ^(i)] and Z^(i)*-Z^(j)* are cancelled. Denote λ_2^*=∑_i∈𝒱λ_2^(i)*, μ^*=∑_i∈𝒱μ^(i)*, we get (<ref>). Since (<ref>) 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 λ_2^*,μ^*,{w_k^(i)*} is an optimal solution to (<ref>), then there must exist Lagrange multipliers Φ^*, {v^(i)*} and {φ^(i)*_k} such that (<ref>) holds. Now choose {λ̂_2^(i)*,μ̂^(i)*,{ŵ_k^(i)*},Ẑ^(i)*} and Lagrange multipliers {Φ̂^(i)*,v̂^(i)*,{φ̂_k^(i)*}} such that ∑_i∈𝒱λ̂_2^(i)*=λ_2^*,∑_i∈𝒱μ̂^(i)*=μ^*,ŵ_k^(i)*=w_k^(i)* and Φ̂^(i)*=Φ^*,  v̂^(i)*=v^(i)*,  φ̂_k^(i)*=φ_k^(i)*,∀ k∈ℰ(i),∀ i∈𝒱. The KKT conditions (<ref>)-(<ref>) and (<ref>) is trivially satisfied by the above construction. What remains to show is that there exists such {Ẑ^(i)*} so that (<ref>) and (<ref>) are satisfied.
	


We first show there exists {Ẑ^(i)*} such that (<ref>) is satisfied. Denote Â^(i)=λ̂_2^(i)*I-μ̂^(i)*1-∑_k∈ℰ(i)ŵ_k^(i)*E_k. Since λ_2^*,μ^*,{w_k^(i)*} satisfy (<ref>) and λ_2^*=∑_i∈𝒱λ̂_2^(i)*, μ^*=∑_i∈𝒱μ̂^(i)*, ŵ_k^(i)*=w_k^(i)*, we know that ∑_i∈𝒱Â^(i)≼ 0. By choosing P^(i)=-Â^(i) for i=1,⋯,N-1 and P^(N)=∑_i=1^N-1Â^(i), we have Â^(i)+P^(i)= 0 for i=1,⋯ N-1 and Â^(N)+P^(N)=∑_i=1^NÂ^(i)≼ 0. Therefore, we have Â^(i)+P^(i)≼ 0 for all i∈𝒱 and ∑_i∈𝒱P^(i)=0. What remains to show is that there exists {Ẑ^(i)*}, such that P^(i)=∑_j∈𝒩(i)(Ẑ^(i)*-Ẑ^(j)*) for all i∈𝒱. Recall that P^(i)=∑_lp^(i)_lB_l, Ẑ^(i)*=∑_lẑ^(i)*_lB_l and denote p_l=[p_l^(1),⋯,p_l^(N)]^T, ẑ_l^*=[ẑ_l^(1)*,⋯,ẑ_l^(N)*]^T. Since ∑_i∈𝒱P^(i)=0, it follows that 1^Tp_l=0 for all l. Therefore, p_l∈ ker(1^T)=Im(L), where L=∑_k∈ℰE_k. This implies p_l can be expressed as Lẑ_l^* for some ẑ_l^*, namely, there exists {Ẑ^̂(̂î)̂} such that P^(i)=∑_j∈𝒩(j)(Ẑ^(i)*-Ẑ^(j)*) for all i∈𝒱 and hence (<ref>) is satisfied.
	


Now we show that {λ̂_2^(i)*,μ̂^(i)*,Ẑ^(i)*} and {Φ̂^(i)*} chosen above satisfy (<ref>). Since Φ̂^(i)*=Φ^*, ŵ_k^(i)*=w_k^(i)* and P^(i)=∑_j∈𝒩(i)(Ẑ^(i)*-Ẑ^(j)*), the left hand side of (<ref>) can be written as λ̂_2^(i)*-∑_k∈ℰ(i)w_k^(i)*tr(E_kΦ^(i)*)+tr(P^(i)Φ^*).
	Recall that P^(i)=-Â^(i)=-λ̂_2^(i)*I+μ̂^(i)*1+∑_k∈ℰ(i)ŵ_k^(i)*E_k,i=1,⋯ N-1. In view of (<ref>), for i=1,⋯,N-1, we have λ̂_2^(i)*-∑_k∈ℰ(i)w_k^(i)*tr(E_kΦ^*)-λ̂_2^(i)*tr(Φ^*)+μ̂^(i)*tr(1Φ^*)+∑_k∈ℰ(i)w_k^(i)*tr(E_kΦ^*)=0.
	
	

	For i=N, P^(N)=∑_i=1^N-1Â^(i)=∑_i=1^N-1λ̂_2^(i)*I-∑_i=1^N-1μ̂^(i)*1-∑_i=1^N-1∑_k∈ℰ(i)ŵ_k^(i)*E_k. In view of (<ref>), (<ref>) and since ∑_i∈𝒱λ̂_2^(i)*=λ_2^*,∑_i∈𝒱μ̂^(i)*=μ^*, we have λ̂_2^(N)-∑_k∈ℰ(N)w_k^(N)*tr(E_kΦ^*)+∑_i=1^N-1λ̂_2^(i)*tr(Φ^*)-∑_i=1^N-1μ̂^(i)*tr(1Φ^*)-∑_i=1^N-1∑_k∈ℰ(i)w_k^(i)*tr(E_kΦ^*)=0.
	Therefore, by the variables construction above, the second part of the statement follows.Now we show the second part of the statement. Since λ̂_2, μ̂ and {ŵ_k^(i)} is a feasible solution of (<ref>), it satisfies the constraints of (<ref>). Since the last two constraints on w_k^(i) is the same for both (<ref>) and (<ref>), what remains to show is that there exists {λ̂_2^(i),μ̂^(i),{ŵ_k^(i)},Ẑ^(i)},i∈𝒱 such that λ̂_2=∑_i∈𝒱λ̂_2^(i), μ̂=∑_i∈𝒱μ̂^(i) and λ̂_2^(i)I-μ̂^(i)1-∑_k∈ℰ(i)ŵ_k^(i)E_k+∑_j∈𝒩(i)(Ẑ^(i)-Ẑ^(j))≼ 0 for all i∈𝒱. 
	We prove this by showing the following equivalent statement: for all  {λ̂_2^(i)} and {μ̂^(i)} that satisfy λ̂_2=∑_i∈𝒱λ̂_2^(i), μ̂=∑_i∈𝒱μ̂^(i), there always exists {Ẑ^(i)} such that λ̂_2^(i)I-μ̂^(i)1-∑_k∈ℰ(i)ŵ_k^(i)E_k+∑_j∈𝒩(i)(Ẑ^(i)-Ẑ^(j))≼ 0 holds for all i∈𝒱.

	Denote Â^(i)=λ̂_2^(i)I-μ̂^(i)1-∑_k∈ℰ(i)ŵ_k^(i)E_k. From the feasibility of {λ̂_2,μ̂,{ŵ_k^(i)}} and λ̂_2=∑_i∈𝒱λ̂_2^(i), μ̂=∑_i∈𝒱μ̂^(i), we know that ∑_i∈𝒱Â^(i)≼ 0. By choosing P^(i)=-Â^(i) for i=1,⋯,N-1 and P^(N)=∑_i=1^N-1Â^(i), we have Â^(i)+P^(i)= 0 for i=1,⋯ N-1 and Â^(N)+P^(N)=∑_i=1^NÂ^(i)≼ 0. Therefore, we have Â^(i)+P^(i)≼ 0 for all i∈𝒱 and ∑_i∈𝒱P^(i)=0. What remains to show is that there exists {Ẑ^(i)}, such that P^(i)=∑_j∈𝒩(i)(Ẑ^(i)-Ẑ^(j)) for all i∈𝒱. Recall that P^(i)=∑_lp^(i)_lB_l, Ẑ^(i)=∑_lẑ^(i)_lB_l and denote p_l=[p_l^(1),⋯,p_l^(N)]^T, ẑ_l=[ẑ_l^(1),⋯,ẑ_l^(N)]^T. Since ∑_i∈𝒱P^(i)=0, it follows that 1^Tp_l=0 for all l. Therefore, p_l∈ ker(1^T)=Im(L), where L=∑_k∈ℰE_k. This implies p_l can be expressed as Lẑ_l for some ẑ_l, namely, there exists {Ẑ^̂(̂î)̂} such that P^(i)=∑_j∈𝒩(j)(Ẑ^(i)-Ẑ^(j)) for all i∈𝒱. Thus the second part of statement follows.
	
  
  

 §.§ Relaxing SDP into NP
  
  (<ref>) can be solved distributedly by using a similar method as <cit.> 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 (<ref>), the problem needs to be relaxed into an NP first. This is because (<ref>) is still an SDP problem, and its inequality matrices constraints would lead to positive semidefinite matrix Lagrangian multipliers Φ^(i). This makes it hard to apply the saddle point dynamics in <cit.> to this problem since by the definition of the projection operator Π_K, it is clear that Π_K:ℝ^n×ℝ^n↦ℝ^n. The projected saddle point dynamics is not defined on the cone of positive semidefinite matrices.




  Now we introduce the convex function proposed by <cit.> which can be used to approximate the largest eigenvalue of a symmetric matrix. Given X∈𝒮^N, function f_ε:𝒮^N↦ℝ and reads f_ε(X)=εln tr(e^X/ε)=εln[∑_i=1^Ne^λ_i(X)/ε] and its derivative with respect to X reads
  
    ∇_Xf_ε(X)=[∑_i=1^Ne^λ_i(X)/ε]^-1[∑_i=1^Ne^λ_i(X)/εu_iu_i^T],

  where (λ_i(X),u_i) are eigen-pairs of X with u_i=1,∀ i.
  It has been proved in <cit.> that 
   
    λ_max(X)≤ f_ε(X)≤λ_max(X)+εln N.

  Hence when ε is sufficiently small, f_ε(X)≈λ_max(X). 




	Consider the following NP problem:
  
    {μ^(i)},{w_k^(i)},{Z^(i)}minimize      ∑_i∈𝒱f_ε_i(X^(i))
       subject to      ∑_k∈ℰ(i)w_k^(i)=1
             w_k^(i)≥ 0,  ∀ k∈ℰ(i),  ∀ i∈𝒱NP

  where X^(i)=-μ^(i)1-∑_k∈ℰ(i)w_k^(i)E_k+∑_j∈𝒩(i)(Z^(i)-Z^(j)) to abbreviate the notation.
  It is clear that (<ref>) is a convex problem, the Slater's condition also holds for (<ref>) and hence strong duality holds. Therefore, KKT conditions becomes necessary and sufficient conditions for (<ref>).
The next proposition shows that how well the approximation could be.
		
		Suppose ∑_i∈𝒱λ_2^(i)* is the optimal objective function value of (<ref>), then
		
    -∑_i∈𝒱λ_2^(i)*≤∑_i∈𝒱f_ε_i(X^(i)*)≤-∑_i∈𝒱λ_2^(i)*+∑_i∈𝒱ε_iln N,

		where X^(i)*=-μ^(i)*1-∑_k∈ℰ(i)w_k^(i)*E_k+∑_j∈𝒩(i)(Z^(i)*-Z^(j)*).
		Moreover, suppose {μ̂^(i)*,{ŵ_k^(i)*},Ẑ^(i)*} is the optimal solution to (<ref>), then
		
    -∑_i∈𝒱λ_2^(i)*   ≤∑_i∈𝒱λ_max(X̂^(i)*)≤∑_i∈𝒱f_ε_i(X̂^(i)*)
       ≤ -∑_i∈𝒱λ_2^(i)*+∑_i∈𝒱ε_iln N,

		where X̂^(i)*=-μ̂^(i)*1-∑_k∈ℰ(i)ŵ_k^(i)*E_k+∑_j∈𝒩(i)(Ẑ^(i)*-Ẑ^(j)*).
		
		
			From the KKT condition (<ref>) of (<ref>), we know that λ_max(X^(i)*)≤ -λ_2^(i)*. 			
			This implies 
			
    ∑_i∈𝒱λ_max(X^(i)*)≤-∑_i∈𝒱λ_2^(i)*.

			On the other hand, by Proposition <ref>, ∑_i∈𝒱X^(i)*=-∑_i∈𝒱μ^(i)*1-∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)*E_k=-μ^*1-∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)*E_k, where λ_2^*, μ^* and {w_k^(i)*} is the optimal solution to (<ref>). 
		We know from (<ref>) that -λ_2^*=λ_max(-μ^*1-∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)*)=λ_max(∑_i∈𝒱X^(i)*) <cit.>. Hence by Proposition <ref>, we have
		
    -∑_i∈𝒱λ_2^(i)*=-λ_2^*=λ_max(∑_i∈𝒱X^(i)*).

		Then it follows from (<ref>) that ∑_i∈𝒱λ_max(X^(i)*)≤λ_max(∑_i∈𝒱X^(i)*). On the other hand, by eigenvalue inequality, we know that ∑_i∈𝒱λ_max(X^(i)*)≥λ_max(∑_i∈𝒱X^(i)*), and hence 
		
    ∑_i∈𝒱λ_max(X^(i)*)=λ_max(∑_i∈𝒱X^(i)*).
				
		Since λ_max(X^(i)*)≤ f_ε_i(X^(i)*)≤λ_max(X^(i)*)+ε_iln N, we have ∑_i∈𝒱λ_max(X^(i)*)≤∑_i∈𝒱f_ε_i(X^(i)*)≤∑_i∈𝒱λ_max(X^(i)*)+∑_i∈𝒱ε_iln N. Hence it follows from (<ref>) and (<ref>) that (<ref>) holds.
		
		

		On the other hand, by (<ref>), it follows that λ_max(X̂^(i)*)≤ f_ε_i(X̂^(i)*) and hence ∑_i∈𝒱λ_max(X̂^(i)*)≤∑_i∈𝒱f_ε_i(X̂^(i)*).
		Note that {-λ_max(X̂^(i)*),μ̂^(i)*,{ŵ_k^(i)*},Ẑ^(i)*} is also a feasible solution to (<ref>), then this implies -∑_i∈𝒱λ_2^(i)*≤∑_i∈𝒱λ_max(X̂^(i)*). In addition, since {μ̂^(i)*,{ŵ_k^(i)*},Ẑ^(i)*} is an optimal solution to (<ref>), it follows that ∑_i∈𝒱f_ε_i(X̂^(i)*)≤∑_i∈𝒱f_ε_i(X^(i)*). Moreover, using (<ref>), we have -∑_i∈𝒱λ_2^(i)*≤∑_i∈𝒱λ_max(X̂^(i)*)≤∑_i∈𝒱f_ε_i(X̂^(i)*)≤∑_i∈𝒱f_ε_i(X^(i)*)≤ -∑_i∈𝒱λ_2^(i)*+∑_i∈𝒱ε_iln N, which proves the statement.	
Proposition <ref> shows that one can have a good approximation on ∑_i∈𝒱λ_2^(i) by choosing ε_i, ∀ i∈𝒱 sufficiently small. Without losing generality, we choose ε_i=ε, ∀ i∈𝒱. Note that (<ref>) explains the reason of the introduction of μ in (<ref>) instead of writing the constraint as λ_2(I-1/N1)≼∑_i∈𝒱∑_k∈ℰ(i)w_k^(i)E_k as in <cit.>. The λ_2I term is needed for the relaxation. It is worth noticing that μ^* does not necessarily equals to λ_2^*/N in (<ref>). In fact, any (λ_2^*,μ^*,{w_k^(i)*}) such that μ^*≥λ_2^*/N is an optimal solution for (<ref>). Also note that f_ε(X) is not a strictly convex function though it is convex. Indeed, for all 0≤α≤ 1, it holds that α f_ε(I)+(1-α)f_ε(2I)=αεln(Ne^1/ε)+(1-α)εln(Ne^2/ε)=εln N+2-α=εln(Ne^2-α/ε)=f_ε((2-α)I)=f_ε(α I+(1-α)2I).
	Hence f_ε(X) is not strictly convex.
  

 §.§ Projected Dynamics and Numerical Examples

  Now we apply the projected dynmaical system to solve (<ref>). To abbreviate the notation, denote x=[x^(1)T,⋯,x^(N)T]^T, where x^(i)=[μ^(i),{w^(i)_k},{z^(i)_l}]^T and v=[v^(1),⋯,v^(N)]^T. 



The projected dynamics for each agent i is given by
  
	
    μ̇^(i)   =⟨∇_X^(i)f_ε(X^(i)),1⟩_M=-∂ℒ_𝒜(x,v)/∂μ^(i),
    ẇ_k^(i)   =Π_ℝ_+(w_k^(i),⟨∇_X^(i)f_ε(X^(i)),E_k⟩_M-v^(i)
       -(∑_p∈ℰ(i)w_p^(i)-1))
       =Π_ℝ_+(w_k^(i),-∂ℒ_𝒜(x,v)/∂ w_k^(i)),  ∀ k∈ℰ(i),

  
  
	
    ż^(i)_l   =-∑_j∈𝒩(i)⟨∇_X^(i)f_ε(X^(i))-∇_X^(j)f_ε(X^(j)),B_l⟩_M
       =-∂ℒ_𝒜(x,v)/∂ z_l^(i),
    v̇^(i)   =∑_k∈ℰ(i)w_k^(i)-1=∂ℒ_𝒜(x,v)/∂ v^(i),

  
  where ℒ_𝒜(x,v)=∑_i∈𝒱{f_ε(X^(i))+v^(i)(∑_k∈ℰ(i)w_k^(i)-1)+1/2(∑_k∈ℰ(i)w_k^(i)-1)^2} and ∇_X^(i)f_ε(X^(i)) is given by (<ref>).
  Note that in (<ref>), (<ref>), 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 ∇_X^(i)f_ε(X^(i)). We would like to remark that, although each time step each agent has to communicate with its neighbour a vector of size 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.
  
	The system (<ref>), (<ref>) is well-defined and the trajectory asymptotically converges to one of the saddle points of (<ref>) for all initial values μ^(i)(0),z_l^(i)(0),v^(i)(0) ∈ℝ, w_k^(i)(0)∈ℝ_+.
  
  
	f_ε(X) has a Lipschitz continuous gradient with respect to x given that X=∑_ix_iA_i, where all A_i are symmetric matrices <cit.>. By Theorem 2.5 in <cit.>, for any initial value μ^(i)(0)∈ℝ, w_k^(i)(0)∈ℝ_+, v^(i)(0)∈ℝ and z_l^(i)(0)∈ℝ, there exists a unique Carathéodory solution which continuously depends on the initial value. Therefore the system (<ref>), (<ref>) is well-defined and by Proposition <ref>, the system (<ref>), (<ref>) asymptotically converges to one of the saddle points of (<ref>).
It seems that when simulating the projected dynamics (<ref>),(<ref>), one has to do eigenvalue decomposition on X^(i) to compute ∇_X^(i) f(X^(i)) at each time step. However, since the factors e^λ_i(X^(i))/ε decrease very rapidly, the gradient numerically only depends on few largest eigenvalues and correpondant eigenvectors <cit.>. 
  
  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 <ref>, one can first choose an ε and get an “optimal” algebraic connectivity under the current choice of ε.
  If the approximation error ε Nln N compared to the “optimal” algebraic connectivity under the current choice of ε is not satisfying (for example, ε Nln N is approximately 10% of the current “optimal” algebraic connectivity), one can decrease ε until the desired relative error is achieved.
  
	We run a simple numerical example using the graph illustrated in Fig. <ref> to show that the variables do converge to the optimal solution. Using CVX, we get the optimal solution of (<ref>): w_1^(1)*=w_2^(3)*=1,  w_1^(2)*=w_2^(2)*=0.5 and λ_2^*=1.5. Forward Euler method is used to discretize (<ref>), (<ref>) and we choose ε=0.01, time step size Δ t=0.01.

	Fig. <ref> shows that the edge weights converge to the optimal solution.
	
  
	Now we consider a more complicated graph generated by ten nodes. 
	Forward Euler method is also used to discretize (<ref>), (<ref>) and we choose different ε and time step size Δ t to illustrate the effect of ε on Δ t. According to <cit.>, the choice of ε affects the Lipschitz constant of ∇_X f_ε(X) as well as the Hessian of f_ε(X). The smaller ε is, the bigger the Lipschitz constant of ∇_X f_ε(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 (<ref>) is 1.141. We do the simulation for t∈[0,50] with different ε and time step sizes. In the end, we get λ_2(L_w) equals to 1.085, 1.091 and 1.128 when choosing ε=10^-2, Δ t=10^-3, ε=10^-3, Δ t=10^-3 and ε=10^-3, Δ t=10^-4, respectively. As illustrated in Fig. <ref>, 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 ε, the limiting algebraic connectivity gets closer to the optimal value of (P_c). This illustrates the relaxation effect. In addition, the evolution of λ_2(L_w) involves a lot of oscillations when ε=10^-3 and Δ t=10^-3, while it behaves much nicer when Δ t=10^-4. This emperically shows that smaller ε requires smaller time step length.
	
  
  
  Now we do a numerical example on a graph with 50 nodes. The optimal value of (<ref>) is computed by CVX. The evolution of λ_2(L_w) is illustrated in Fig. <ref>. We choose ε=10^-4 and Δ t=10^-2. In this case, the approximation error ε Nln N≈ 0.0196 and we can see from Fig. <ref> that the gap between the limiting algebraic connectivity and the optimal value computed by CVX is below 0.0196.
  
  Fig. <ref> illustrates that a smaller ε requires a smaller step size when discretizing system (<ref>), (<ref>). This means that when the number of nodes N goes large, in order to get a good approximation of (<ref>), we need a very small ε 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 ∑_k∈ℰ(i)w_k^(i)=1 in (<ref>) as ∑_k∈ℰ(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 (<ref>), we conclude that {λ_2^(i)*,μ^(i)*,{w_k^(i)*},Z^(i)*} is the optimal solution to (<ref>) iff {aλ_2^(i)*,aμ^(i)*,a{w_k^(i)*},aZ^(i)*} is the optimal solution to (P_M) (since all the optimality conditions are linear). Using Proposition <ref>, we can conclude that -∑_i∈𝒱λ_2^(i)*≤1/a∑_i∈𝒱λ_max(X̂^(i)*)≤1/a∑_i∈𝒱f_ε_i(X̂^(i)*)≤ -∑_i∈𝒱λ_2^(i)*+1/aε Nln N provided that {μ̂^(i)*,{ŵ_k^(i)*},Ẑ^(i)*} is the optimal solution to (NP_M). Therefore, apart from choosing ε 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 ε 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.
  
  Consider a graph with 30 nodes. If we do not use the methodology above, namely, a=1, ε needs to be at least 2.1416× 10^-5 so that the relative approximation error of the optimal algebraic connectivity ε Nln N/λ_2(L_w^*) is within 5%. For comparison, if we fix ε=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. <ref>. 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.
  
  
  
  	Consider the same graph used in Example <ref>. We use forward Euler for discretization. Same initial values, time step sizes and ε are used (ε=0.5, Δ t=0.02). Fig. <ref> shows how a affects the approximation error. 
  	
  
  



  The complexity per iteration for agent i is 𝒪(|ℰ(i)|· N^2). Since the right hand side of (<ref>), (<ref>) involves only with special matrices such as 1, E_k and B_l, one does not need to do matrix multiplication and hence the complexity is greatly reduced.
  
  We test the algorithm on larger scale networks and plot the ratio between running time and |ℰ(i)|· 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. ε=0.5,Δ=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 (<ref>), (<ref>) is smaller than 10^-3. The result is shown as Fig. <ref>. It can be seen that the ratio between running time and |ℰ(i)|· N^2 is approximately constant when N changes.
  
    
  


§ 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 <cit.>, 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 ε 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. 
  
  model5-names

   
