Concensus-Based ALADIN Method to Faster the Decentralized Estimation of Laplacian Spectrum

Abstract

With the upcoming fifth Industrial Revolution, humans and collaborative robots will dance together in production. They themselves act as an agent in a connected world, understood as a multi-agent system, in which the Laplacian spectrum plays an important role since it can define the connection of the complex networks as well as depict the robustness. In addition, the Laplacian spectrum can locally check the controllability and observability of a dynamic controlled network, etc. This paper presents a new method, which is based on the Augmented Lagrange based Alternating Direction Inexact Newton (ALADIN) method, to faster the convergence rate of the Laplacian Spectrum Estimation via factorization of the average consensus matrices, that are expressed as Laplacian-based matrices problems. Herein, the non-zero distinct Laplacian eigenvalues are the inverse of the stepsizes $\{{\alpha }_t,t=1,2,\dots \}$ of those matrices. Therefore, the problem now is to carry out the agreement on the stepsize values for all agents in the given network while ensuring the factorization of average consensus matrices to be accomplished. Furthermore, in order to obtain the entire Laplacian spectrum, it is necessary to estimate the relevant multiplicities of these distinct eigenvalues. Consequently, a non-convex optimization problem is formed and solved using ALADIN method. The effectiveness of the proposed method is evaluated through the simulation results and the comparison with the Lagrange-based method in advance.

1. Introduction

Leaders around the world obviously prefer the present era of connectivity as the Fourth Industrial Revolution [2]. Industry 4.0 has been significantly contributing to the transformation of many industries such as transportation, manufacturing, health-care, agriculture, etc. via enabling data transmission and integration between disciplines. Today, we are in the fourth one, a generation of connection between our physical, digital, social, and biological worlds. In particular, data and information from these different areas have been made available and have been connected in complex and dense networks. For better integration and utilization, control aspect of these complex networks, researchers from different communities have brought different contributions varying from topology inference to control strategy, deputizing for interacting systems, which are modeled by graphs, whose vertices represent the components of the system while edges stand for the interactions between these components.

In the last decade there has been dramatic increasing number of publications in the cooperative control of multi-agent systems. In control of multi-agent systems, the performance of the whole system depends on both structure and the connections between individuals of the systems. Here, the total connections can be defined by the graph Laplacian matrix, and its spectrum is involved in some useful properties of the network system [3,4]. For instance, the second smallest graph Laplacian eigenvalue, i.e., the so-called algebraic connectivity of the graph, has the main role in the convergence time of various distributed algorithms as well as the performance and robustness of dynamical systems [5]. Conceptually, agents share their information with each other to achieve common objectives, relative position information, or common control algorithms. This is called consensus problem [6,7], where a group of agents’ approaches average consensus in an undirected network under simple linear iteration scheme.

It is well known that in a multi-agent system, consensus is achieved if and only if the network is connected (the algebraic connectivity) being strictly greater than zero [8]. On the other hand, the largest Laplacian eigenvalue is an important factor to decide the stability of the system. For example, minimizing the spectral radius in [9] leads to maximize the robustness of the network to time delays under a linear consensus protocol. Furthermore, to speed up consensus algorithms, the optimal Laplacian-based consensus matrix is obtained with a stepsize which is the inverse of the sum of the smallest and the largest non-zero graph Laplacian eigenvalue [10]. Moreover, the authors in [11,12] have proved that the spectrum of the Laplacian matrix can be used to design consensus matrices to obtain average consensus in finite number of steps.

In order to investigate network efficiency, structural robustness of a network which is related to its performance despite changes in the network topology [13] has been also studied. The concept of natural connectivity as a spectral measure of robustness was introduced in [14]. It is expressed in mathematical form as the average eigenvalue of the adjacency matrix of the graph representing the network topology. The Laplacian spectrum $sp\left(\mathbf{L}\right)=\{{\lambda }_1^{m_1},\dots ,{\lambda }_i,\dots \}$ can also be employed to compute the robustness indices, for instance, the number of spanning trees and the effective graph resistance (Kirchhoff index) [15]. The smaller (or greater) the Kirchhoff index (or the number of spanning trees) is, the more robust the network becomes. In addition, it has been pointed out that adding an edge strictly decreases the Kirchhoff index and hence increases the robustness. In [16], the authors have proposed a method to monitor collaboratively the robustness of the networks partitioned into sub-networks by Kirchhoff index ${\mathcal{R}}_{\mathbf{L}}=N{\sum }_{i=2}^{D+1}\frac{m_i}{{\lambda }_i}$. Here, an Alternating Direction of Multipliers Method (ADMM)-based algorithm was employed to perform the factorization of the averaging matrix and to compute the average degree of the network concurrently. However, the main point in this work was the reformulation into the convex optimization problem, which is convenient to make use of the ADMM method to solve the problem. In addition to that, the impact of the Laplacian spectrum into power systems is expressed through energy management in smart grids [17] and the determination of the grid robustness against low frequency disturbance in [18]. In this work, in the framework of spectral graph theory, the authors reveal that the decomposition of frequency signal along scaled Laplacian spectrum when the damping-inertia ratios are uniform across buses not only makes the system respond faster but also helps lower the system nadir after a disturbance. In dynamic network systems, the spectrum of Laplacian matrix can also be utilized for locally checking the controllability and the observability [19].

From a short literature survey above, it is obvious that the Laplacian spectrum plays an important role in many fields. For instance, Laplacian spectrum can be used to design consensus matrices [11,12], to compute these robustness indices [15,16,17,18], or to check the controllability and the observability [19]. Hence, it is desirable to have an efficient method for monitoring the Laplacian spectrum of a dynamics network system.

One thing to remark here is that if the global network topology is known in a-priori, the Laplacian matrix can be easily deduced. However, implementing a centralized structure is an expensive task due to the high computational cost, the heavy communication infrastructure aspect and the problem from large dimensionality. Additionally, if there is a failure problem from one point, it will affect the whole network. Therefore, our study is restricted to the assumption that the network topology (represented by the Laplacian Matrix) is unknown at the first glance. A dominant contribution of this paper is the possibility of implementing this monitoring scheme in a decentralized manner.

In this paper, we present an Augmented Lagrangian based Alternating Direction Inexact Newton (ALADIN) method to estimate the Laplacian spectrum in decentralized scheme for dynamic controlled networks. The key feature of this paper is the direct solution to non-convex optimization for Laplacian spectrum estimation using ALADIN method. To simplify, the scope of this study is restricted to networks performing noise-free as well as the number of the agents N in the network in known in-priori by using random walk algorithm [20]. The network is modeled, then Laplacian eigenvalues and average consensus are retrieved respectively. Since the Laplacian spectrum matrix is not directly computable for undetermined network topology, the decentralized estimation of the Laplacian spectrum has been introduced with three main approaches in the recent literature: Fast Fourier Transform (FFT)-based methods [21,22], local eigenvalue decomposition of given observability-based matrices [23], and distributed factorization of the averaging matrix ${\mathbf{J}}_N=\frac{1}{N}{\mathbf{11}}^T$ [24]. FFT-based methods require a specific protocol. However, they do not make use of the available measurements coming from the consensus protocol. On the other hand, the method in [23] allows using the transient of the average consensus protocol but for several consecutive initial conditions. The distributed factorization of the averaging matrix in [24] yields the inverses of non-zero Laplacian eigenvalues and can be solved as a constrained consensus problem. The Laplacian eigenvalues can be deduced as the inverse of the stepsizes in each estimating factor, where these factors are constrained to be structured as Laplacian based consensus matrices. In [1], authors have applied a gradient descent algorithm to solve this optimization problem in which only local minima was guarantees accompany with slow convergence rate. In order to solve this annoying issue, in [16,24], the authors have introduced an interesting way by reformulating a non-convex optimization problem in [1] into convex one and solved by applying an ADMM-based method. However, this is an indirect approach obtaining by an adequate re-parameterization. In this paper, we inherit the idea in [1] to form the non-convex optimization for decentralized estimation of Laplacian spectrum and then directly solve it using the ALADIN method that was proposed by [25]. The proposed approach is then evaluated with two network structures for performance evaluation in comparison with gradient descent method [1].

In this paper, we firstly introduce the background of average consensus and state the problem in Section 2, then present the distributed estimation of Laplacian spectrum in Section 3. The structure of this section can be illustrated as in Figure 1. Before concluding the paper, the simulation results are described in Section 4 to evaluate the efficiency of the proposed method.

2. Background and Problem Statement

Consider a dynamic network, in which interconnection is represented by $G(V,E)$, an undirected graph with components’ set V and links’ set E, consisting of $N=\left|V\right|$ nodes, let us denote by ${\mathcal{N}}_i=\{j\in V:(i,j)\in E\}$ the set of neighbors of node i and $d_i=\left|{\mathcal{N}}_i\right|$ its degree. Interactions between nodes can be captured by the Laplacian matrix $\mathbf{L}\in {\Re }^N$ with entries $l_{ii}=d_i$, $l_{ij}=-1$ if $j\in {\mathcal{N}}_i$ and $l_{ij}=0$ elsewhere. Denote the Laplacian spectrum by $sp\left(\mathbf{L}\right)=\{{\lambda }_1^{m_1},{\lambda }_2^{m_2},\dots ,{\lambda }_{D+1}^{m_{D+1}}\}$, where the different Laplacian eigenvalues are in increasing order $0={\lambda }_1<{\lambda }_2<\dots <{\lambda }_{D+1}$ and superscripts stand for multiplicities $m_i=m\left({\lambda }_i\right)$, while ${\mathbf{S}}_2=\{{\lambda }_2,\dots ,{\lambda }_{D+1}\}$ stands for the set of the non-zero distinct Laplacian eigenvalues.

2.1. Average Consensus

For each node $i\in V$, let $x_i\left(t\right)$ denotes the value of node i at timestep t. Define $\mathbf{x}\left(t\right)={[x_1\left(t\right),x_2\left(t\right),\dots ,x_N\left(t\right)]}^T$, where N is the number of nodes in the network. Average consensus algorithms can be achieved by using the following linear iteration scheme as follows:

$$\mathbf{x}\left(t\right)=({\mathbf{I}}_N-\alpha \mathbf{L})\mathbf{x}(t-1),$$

(1)

where $\alpha$ is an appropriately selecting stepsize [26], by which all nodes converge asymptotically to the same value $\overline{x}$ that is the average of the initial ones $\overline{x}\mathbf{1}={lim}_{t\to \infty }\mathbf{x}\left(t\right)=\frac{1}{N}{\mathbf{11}}^T\mathbf{x}\left(0\right)$.

On the other hand, it has been shown in [11,12] that the average consensus matrix can be factored as

$$\prod _{t=D}^1{\mathbf{W}}_t=\frac{1}{N}{\mathbf{11}}^T,$$

(2)

where ${\mathbf{W}}_t={\vartheta }_t{\mathbf{I}}_N+{\alpha }_t\mathbf{L}$, ${\vartheta }_t$ and ${\alpha }_t$ being parameters to be designed. In [12], the solution was given by ${\vartheta }_t=1$ and ${\alpha }_t=-\frac{1}{{\lambda }_{t+1}}$, ${\lambda }_t$ being a non-zero Laplacian eigenvalue. Owing to the above factorization, average consensus can then be reached in D steps, D being the number of distinct non-zero Laplacian eigenvalues:

$$\overline{\mathbf{x}}=\mathbf{x}\left(D\right)=\prod _{t=D}^1{\mathbf{W}}_t\mathbf{x}\left(0\right)=\frac{1}{N}{\mathbf{11}}^T\mathbf{x}\left(0\right) \mathrm{for} \mathrm{all} \mathbf{x}\left(0\right)\in {\Re }^N.$$

(3)

2.2. Problem Statement

It can be noted that by factorizing the average consensus matrix, while constraining the factor matrices to be in the form ${\mathbf{I}}_N-{\alpha }_t\mathbf{L}$, the eigenvalues of the Laplacian matrix as the inverse of ${\alpha }_t$ can be deduced. The uniqueness has been proved in [1].

Lemma 1.

[1] Let ${\lambda }_2,\cdots ,{\lambda }_{D+1}\ne 0$ be the D distinct non-zero eigenvalues of the graph Laplacian matrix $\mathit{L}$, then, up to permutation, the sequence ${\left\{{\alpha }_i\right\}}_{i=1,\cdots ,D}$, with ${\alpha }_i=\frac{1}{{\lambda }_{i+1}}$, $i=1,2,\cdots ,D$, is the unique sequence allows getting the minimal factorization of the average consensus matrix as $\frac{1}{N}{\mathbf{11}}^T={\prod }_{i=1}^D({\mathit{I}}_N-{\alpha }_i\mathit{L})$.

Therefore, in order to implement the proposed method, the knowledge of the network should be known. Meaning that, the number of the components N of the given network should be known by adding a learning mechanism as a configuration step. Practically, in most systems where communications are involved, learning sequences are used for communication channel identification or for synchronization. In [20], the authors have proposed a method using random walks to estimate the global properties of large connected undirected graphs such as number of vertices, edges, etc. However, it is not in the scope of this paper. Indeed, assuming that the number of agents N is known in a-priori, a consensus protocol in [10] is to be uploaded to each agent to compute the average consensus value $\overline{x}$. The main task in our study is to estimate the whole Laplacian spectrum.

3. Distributed Estimation of Laplacian Spectrum

Given an initial input-output pair $\{\mathbf{x}\left(0\right),\overline{\mathbf{x}}\}$, with $\overline{\mathbf{x}}=\frac{1}{N}{\mathbf{11}}^T\mathbf{x}\left(0\right)$, the matrix factorization problem (3) is equivalent to minimize the cost function $E\left(\mathbf{W}\right)={∥\mathbf{x}\left(D\right)-\overline{\mathbf{x}}∥}^2$ that can also be rewritten as follows:

$$E\left(\mathbf{W}\right)={∥\prod _{t=D}^1{\mathbf{W}}_t\mathbf{x}\left(0\right)-\overline{\mathbf{x}}∥}^2,$$

(4)

where D is the number of steps before reaching average consensus and ${\mathbf{W}}_t={\mathbf{I}}_N-{\mathit{\alpha }}_t\mathbf{L}$.

Note that there is no need for a central node to set the initial input-output pair. Indeed, such a pair can be obtained after running a standard average consensus algorithm. Each node keeps in memory its own initial value and the consensus value.

Solving this factorization problem consists in finding the sequence of stepsize ${\left\{{\mathit{\alpha }}_t\right\}}_{t=1,\cdots ,D}$. It is obvious that ${\mathit{\alpha }}_t$ are global parameters. To relax these constraints, define the factor matrices as ${\mathbf{W}}_t={\mathbf{I}}_N-{\mathbf{\Lambda }}_t\mathbf{L}$, where ${\mathbf{\Lambda }}_t=diag\left({\mathit{\alpha }}_t\right)$, ${\mathit{\alpha }}_t=[{\alpha }_{t,1},{\alpha }_{t,2},\dots ,{\alpha }_{t,N}]$, $t=1,2,\dots ,D$. The problem above can be reformulated as a constrained consensus problem, that is to compute the sequence of stepsize $\left\{{\mathit{\alpha }}_t\right\}$ so that ${\alpha }_{t,1}={\alpha }_{t,2}=\dots ={\alpha }_{t,N}$. Moreover, in Section 2.1, D is denoted as number of non-zero distinct Laplacian eigenvalues. However, in this work, Laplacian matrix is assumed to be not-known in-a-priori. Therefore, the authors have assigned D as $h=N-1$ since N can be estimated in the configuration step through the Random Walk Algorithm proposed in [20]. Furthermore, the Laplacian Spectrum estimation procedure is divided in following stages:

• Stage 1: Distributed estimation of the set of non-zero Laplacian eigenvalues ${\mathcal{S}}_1=\{{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_h\}$, composing of the set of D non-zero distinct Laplacian eigenvalues ${\mathcal{S}}_2=\{{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_D\}$.
• Stage 2: Eliminating the wrong eigenvalues in the set ${\mathcal{S}}_1$ to obtain the set ${\mathcal{S}}_2$
• Stage 3: Estimating the multiplicities $\mathbf{m}$ corresponding to each eigenvalues in the set ${\mathcal{S}}_2$ to achieve the whole Laplacian spectrum $sp\left(\mathbf{L}\right)$.

3.1. Distributed Estimation of Non-Zero Laplacian Eigenvalues

For distributively carrying out the factorization of the average consensus matrix as factors of Laplacian based consensus matrices, the idea is to minimize the disagreement between neighbors on the value of ${\mathit{\alpha }}_t$ while ensuring that the factorization of the average consensus matrix is achieved. Such a factorization is assessed by constraining the values of the nodes after h iterations of the consensus algorithm to be equal to the average of the initial values:

$$\begin{array}{cc}\underset{{\mathit{\alpha }}_t\in R^{N\times 1},t=1,2,\dots ,h}{\min }& \frac{1}{2}\sum _{t=1}^h\sum _{i\in V}\sum _{j\in {\mathcal{N}}_i}{({\alpha }_{t,j}-{\alpha }_{t,i})}^2\hfill \\ \mathrm{subject}\mathrm{to}& \mathbf{x}\left(h\right)=\overline{\mathbf{x}}\hfill \end{array}$$

(5)

or, it can be rewritten in the following form:

$$\begin{array}{cc}\underset{{\mathit{\alpha }}_t\in R^{N\times 1}}{\min }& \frac{1}{2}\sum _{t=1}^h{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t.\hfill \\ \mathrm{subject}\mathrm{to}& \mathbf{x}\left(h\right)=\overline{\mathbf{x}}\hfill \end{array}$$

(6)

This optimization has been solved by applying Augmented Lagrange Method [1]. However, the disadvantage of this method is the slow convergence rate due to the fact that it is a non-convex optimization problem. To overcome this unexpected issue, the authors have suggested an interesting variant by converting the non-convex function into the convex one. By that, the optimization can be easily and effectively solved by Alternating Direction of Multipliers Methods (ADMM) [16,27]. In this paper, we proposed a method that can solve a non-convex problem effectively by employing ALADIN method, which is described as following:

(1)

Step 1: ALADIN solves in parallel a sequence of equality-constrained non-linear problems (NLP) by introducing an augmented variables $\mathbf{y}$ as follows:

$$\begin{array}{cc}\hfill \underset{{\mathit{\alpha }}_t\in {\mathbf{R}}^{N\times 1},t=1,2,\dots ,h}{min}& \frac{1}{2}{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t+{\mathbf{\lambda }}^T({\mathit{\alpha }}_t-{\mathbf{y}}_t)+\frac{\rho }{2}{\parallel ({\mathit{\alpha }}_t-{\mathbf{y}}_t)\parallel }_{{\mathbf{\Sigma }}_t}^2\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \mathbf{x}\left(h\right)-\overline{\mathbf{x}}=0|\mathit{\beta }\hfill \\ \hfill & {\mathbf{y}}_{t,i}={\mathbf{y}}_{t,j}i=1,\dots ,N;j\in {\mathcal{N}}_i\left(\mathcal{C}\right)\hfill \end{array}$$

(7)

where $\rho ,\beta$ are penalty parameter and multiplier of the equality constraint, respectively. $\mathcal{C}=\{{\mathbf{y}}_t:y_{t,i}=y_{t,j},i=1,\dots ,N;j\in {\mathcal{N}}_i\}$.

One thing to note here is that with the given initial information $x_i\left(0\right),i=1,2,\dots ,N$, running a standard consensus algorithm can determine the average value $\overline{x}=\frac{1}{N}{\sum }_{i=1}^N{x}_i\left(0\right)$. In addition to that, the positive semi-definite scaling matrices ${\Sigma }_t$ can be randomly initialized or even be an identity matrix ${\mathbf{I}}_N$. In this optimization problem (7), augmented variables ${\mathbf{y}}_t$ are introduced with respect to the constraint $\mathcal{C}$. However, this NLP is to be solved to define the variables ${\mathit{\alpha }}_t$, hence, the constraint $\mathcal{C}$ is going to be relaxed.

The solution ${\mathit{\alpha }}_t[k+1],\mathit{\beta }[k+1]$ obtained from (7) is then used to check the stopping criteria the the next steps of the ALADIN-based algorithm procedure with k being an iteration of the optimization process. Herein, if $\frac{1}{2}{\sum }_{t=1}^h{\sum }_{i\in V}{\sum }_{j\in {\mathcal{N}}_i}{({\alpha }_{t,j}-{\alpha }_{t,i})}^2<ϵ$ and $\rho \parallel {\Sigma }_t({\mathit{\alpha }}_t-{\mathbf{y}}_t){\parallel }_1\le ϵ$, then one can get ${\mathit{\alpha }}_t^{*}$ as well as the Algorithm stops.

(2)

Step 2: The Gradients, Jacobian matrices, and Hessian matrices are estimated for the next quadratic programming (QP) subproblems the as follows:

$$\begin{array}{ccc}\hfill {\mathbf{g}}_t& =& \frac{\partial }{\partial {\mathit{\alpha }}_t}\{\frac{1}{2}{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t+{\mathit{\beta }}_t^T(\mathbf{x}\left(h\right)-\overline{\mathbf{x}})\}\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\mathbf{C}}_t& =& \frac{\partial }{\partial {\mathit{\alpha }}_t}{(\mathbf{x}\left(h\right)-\overline{\mathbf{x}})}^T\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\mathbf{B}}_t& =& \frac{{\partial }^2}{\partial {\mathit{\alpha }}_t^2}\left(\frac{1}{2}{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t\right).\hfill \end{array}$$

(10)

(3)

Step 3: Analogically to inexact SQP method, the QP problem is solved to find the $\Delta {\mathit{\alpha }}_t\left[k\right]$ and the affine multiplier ${\lambda }_{QP}\left[k\right]$ as follows:

$$\begin{array}{cc}\hfill \underset{\Delta \mathit{\alpha }\in {\mathbf{R}}^{N\times h},s\in {\mathbf{R}}^{N\times 1},t=1,2,\dots ,h}{min}& \sum _{t=1}^h\{\frac{1}{2}\Delta {\mathit{\alpha }}_t^T{\mathbf{B}}_t\Delta {\mathit{\alpha }}_t+{\mathbf{g}}_t^T\Delta {\mathit{\alpha }}_t\}+{\mathbf{\lambda }}^T\mathbf{s}+\frac{\mu }{2}{\parallel \mathbf{s}\parallel }^2\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \sum _{t=1}^h({\mathit{\alpha }}_t+\Delta {\mathit{\alpha }}_t-{\mathbf{y}}_t)=\mathbf{s}|{\mathbf{\lambda }}_{QP}\hfill \\ \hfill & {\mathbf{C}}_t\Delta {\mathit{\alpha }}_t=0,t=1,2,\dots ,h|\eta \hfill \end{array}$$

(11)

where $\mathbf{s}$ is the slack variable, introduced into the QP sub-problem to attenuate the numerical reasons when the penalty parameter $\mu$ becomes large.

(4)

Step 4: The final step is to update $\mathbf{\lambda }[k+1],{\mathbf{y}}_t[k+1],{\mathbf{B}}_t[k+1]$:

$$\begin{array}{ccc}\hfill \mathbf{\lambda }[k+1]& =& {\lambda }_{QP}\left[k\right]\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{y}}}_t\left[k\right]& =& {\mathit{\alpha }}_t\left[k\right]+\Delta {\mathit{\alpha }}_t\left[k\right]\hfill \end{array}$$

(13)

This update rule is relevant to the full-size step where $a_1=a_2=a_3=1$ for the steplength computation, which proposed in [25]. Then, projecting $\widehat{\mathbf{y}}\left[k\right]$ on the constraint $\mathcal{C}$ to derive $\mathbf{y}[k+1]$.

The steps are repeated until the stopping criteria is satisfied.

In order to derive the distributed algorithm, let us take a closer look at each step of the proposed ALADIN-based method.

3.1.1. Implementation of Decoupled Nonlinear Problems

Firstly, in order to solve the decoupled NLP (7) for $t=1,\dots ,h$, the Augmented Lagrange method is applied here. Hence, we introduce the Augmented Lagrange function with the Lagrange multiplier $\mathit{\beta }$ and penalty parameter $\mathbf{c}$ as below:

$${\mathbf{H}}_{1,t}=\frac{1}{2}{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t+{\mathit{\lambda }}^T({\mathit{\alpha }}_t-{\mathbf{y}}_t)+\frac{\rho }{2}\parallel ({\mathit{\alpha }}_t-{\mathbf{y}}_t){\parallel }_{{\Sigma }_t}^2+{\beta }^T(\mathbf{x}\left(h\right)-\overline{\mathbf{x}})+\frac{c}{2}{\parallel \mathbf{x}\left(h\right)-\overline{\mathbf{x}}\parallel }_2^2$$

(14)

The solution of this problem can be obtained by applying a gradient descent method iteratively:

$${\mathit{\alpha }}_t[k+1]={\mathit{\alpha }}_t\left[k\right]-b\frac{\partial {\mathbf{H}}_{1,t}}{\partial {\mathit{\alpha }}_t}$$

(15)

$$\mathit{\beta }[k+1]=\mathit{\beta }\left[k\right]+c(\mathbf{x}\left(h\right)-\overline{\mathbf{x}})$$

(16)

where b stands for stepsizes of the gradient descent method, which can be chosen by fix constants or be determined by deploying a line-search algorithms such as Wolfe Condition, Back-stracking, or Armijo ones in [28], while c dedicates for the penalty parameter of the Augmented Lagrange method.

Lemma 2.

The derivatives of this Lagrange function (14) is obtained as follows:

$$\frac{\partial {\mathit{H}}_{1,t}}{\partial {\mathit{\alpha }}_t}=\mathit{L}{\mathit{\alpha }}_t+\mathit{\lambda }+\rho {\mathbf{\Sigma }}_t({\mathit{\alpha }}_t-{\mathit{y}}_t)-dia{g}^{-1}\left({\mathit{\alpha }}_t\right)diag({\mathit{x}}_{t-1}-{\mathit{x}}_t){\mathit{\delta }}_t-dia{g}^{-1}\left({\mathit{\alpha }}_t\right)diag({\mathit{x}}_{t-1}-{\mathit{x}}_t){\mathit{e}}_t$$

(17)

where ${\mathit{\delta }}_h=\beta$, and ${\mathit{\delta }}_t={\mathit{W}}_{t+1}{\mathit{\delta }}_{t+1}$, while ${\mathit{e}}_h=\mathit{x}\left(h\right)-\overline{\mathit{x}}$ and ${\mathit{e}}_t={\mathit{W}}_{t+1}{\mathit{e}}_{t+1}.$

The proof is showed in the Appendix A.

3.1.2. Implementation of the Coupling Quadratic Programming (QP)

Now, we apply the Karush–Kuhn–Tucker (KKT) conditions for solving the quadratic programming (QP) (11). The Augmented Lagrange Function is described as follows:

$$\begin{array}{cc}\hfill {\mathbf{H}}_2& =\sum _{t=1}^h\left\{\frac{1}{2}\Delta {\mathit{\alpha }}_t^T{\mathbf{B}}_t\Delta {\mathit{\alpha }}_t+{\mathbf{g}}_t^T\Delta {\mathit{\alpha }}_t+{\mathit{\lambda }}_{QP}^T\Delta {\mathit{\alpha }}_t\right\}+({\mathit{\lambda }}^T-{\mathit{\lambda }}_{QP}^T)\mathbf{s}+\frac{\mu }{2}{\parallel \mathbf{s}\parallel }^2\hfill \\ & +{\mathit{\lambda }}_{QP}^T\sum _{t=1}^h({\mathit{\alpha }}_t-{\mathbf{y}}_t)+\sum _{t=1}^h{\eta }_t^T{\mathbf{C}}_t\Delta {\mathit{\alpha }}_t\hfill \end{array}$$

(18)

The KKT conditions, which are showed in the Appendix B yields a system of equations as follows:

$$\{\begin{array}{ccc}\hfill {\mathbf{B}}_t\Delta {\alpha }_t-{\mathit{\lambda }}_{QP}+{\mathbf{C}}_t^T{\eta }_t& =& -g_t,t=1,2,\dots ,h\hfill \\ \hfill \sum _{t=1}^h\Delta {\mathit{\alpha }}_t-\frac{1}{\mu }{\mathit{\lambda }}_{QP}& =& -\frac{1}{\mu }\mathit{\lambda }-\sum _{t=1}^h({\mathit{\alpha }}_t-{\mathbf{y}}_t)\hfill \\ \hfill {\mathbf{C}}_t\Delta {\mathit{\alpha }}_t& =& 0\hfill \end{array}$$

Solutions of this system of equations are $\Delta {\alpha }_t^{*},{\mathit{\lambda }}_{QP}^{*},{\eta }_t^{*}$ respectively in the equivalent matrix form as follows:   

$$\left(\begin{array}{ccccccccc}{\mathbf{B}}_1& \mathbf{0}& \dots & \mathbf{0}& \mathbf{I}& {\mathbf{C}}_1^T& \mathbf{0}& \dots & \mathbf{0}\\ \mathbf{0}& {\mathbf{B}}_2& \dots & \mathbf{0}& \mathbf{I}& \mathbf{0}& {\mathbf{C}}_2^T& \dots & \mathbf{0}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \dots & {\mathbf{B}}_h& \mathbf{I}& \mathbf{0}& \mathbf{0}& \dots & {\mathbf{C}}_h^T\\ \mathbf{I}& \mathbf{I}& \dots & \mathbf{I}& -\frac{1}{\mu }\mathbf{I}& \mathbf{0}& \mathbf{0}& \dots & \mathbf{0}\\ {\mathbf{C}}_1& \mathbf{0}& \dots & \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \dots & \mathbf{0}\\ \mathbf{0}& {\mathbf{C}}_2& \dots & \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \dots & \mathbf{0}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \dots & {\mathbf{C}}_h& \mathbf{0}& \mathbf{0}& \mathbf{0}& \dots & \mathbf{0}\end{array}\right)\left(\begin{array}{c}\Delta {\alpha }_1^{*}\\ \Delta {\alpha }_2^{*}\\ ⋮\\ \Delta {\alpha }_h^{*}\\ {\mathit{\lambda }}_{QP}^{*}\\ {\eta }_1^{*}\\ {\eta }_2^{*}\\ ⋮\\ {\eta }_h^{*}\end{array}\right)=\left(\begin{array}{c}-{\mathbf{g}}_1\\ -{\mathbf{g}}_2\\ ⋮\\ -{\mathbf{g}}_h\\ -\frac{\mathit{\lambda }}{\mu }-\sum _{t=1}^h({\mathit{\alpha }}_t-{\mathbf{y}}_t)\\ \mathbf{0}\\ \mathbf{0}\\ ⋮\\ \mathbf{0}\end{array}\right)$$

(19)

This linear system can be solved by any linear solver. One thing to note here is the adaptive parameter $\mu$. We can start with a quite small $\mu$ and adapt it during the optimization progress to relax the coupling conditions.

3.1.3. Implementation of an ADMM-Based Algorithm

Now, the update steps (12) and (13) are executed. Since the achieved ${\widehat{y}}_{t,i}$ have to agree with the constraint ($\mathcal{C}$), we solve the following optimization problem:

$$\begin{array}{cc}\hfill \underset{{\mathbf{y}}_i\in {\Re }^{h\times 1}}{min}& \frac{1}{2}\sum _{i=1}^N{\parallel {\mathbf{y}}_i-{\widehat{\mathbf{y}}}_i\parallel }^2\hfill \\ \hfill \mathrm{subject}\mathrm{to}& {\mathbf{y}}_j={\mathbf{y}}_{i}i=1,\dots ,N;j\in {\mathcal{N}}_i\left(\mathcal{C}\right)\hfill \end{array}$$

To solve this optimization problem, an Alternating Direction Method of Multipliers (ADMM) in [16,27] can be employed by introducing an augmented parameters ${\mathbf{z}}_{ij}$. Hence, the optimization can be rewritten as follows:

$$\begin{array}{cc}\hfill \underset{{\mathbf{y}}_t}{min}& \frac{1}{2}\sum _{i=1}^N{\parallel {\mathbf{y}}_i-{\widehat{\mathbf{y}}}_i\parallel }^2\hfill \\ \hfill \mathrm{subject}\mathrm{to}& {\mathbf{y}}_i={\mathbf{z}}_{ij}i=1,\dots ,N;j\in {\mathcal{N}}_i\hfill \\ \hfill & {\mathbf{z}}_{ji}={\mathbf{z}}_{ij}\hfill \end{array}$$

(20)

The Augmented Lagrange Function is defined as follows:

$${\mathbf{H}}_3(\mathbf{y},\mathbf{z},\mathbf{\tau })=\frac{1}{2}\sum _{i=1}^N\parallel {\mathbf{y}}_i-{\widehat{\mathbf{y}}}_i{\parallel }^2+\sum _{j\in {\mathcal{N}}_i}{\mathbf{\tau }}_{ij}^T({\mathbf{y}}_i-{\mathbf{z}}_{ij})+\frac{\nu }{2}\sum _{i=1}^N\parallel {\mathbf{y}}_i-{\mathbf{z}}_{ij}\parallel .$$

The ADMM solution acts in three steps repetitively until the tolerance achieves:

• Compute ${\mathbf{y}}_i$:

  $${\mathbf{y}}_i[p+1]={(1+\nu d_i)}^{-1}\{{\widehat{y}}_i\left[k\right]+\nu \sum _{j\in {\mathcal{N}}_i}{\mathbf{z}}_{ij}\left[p\right]-\sum _{j\in {\mathcal{N}}_i}{\mathbf{\tau }}_{ij}\left[p\right]\}.$$

  (21)
• Compute ${\mathbf{z}}_{ij}$:

  $${\mathbf{z}}_{ij}[p+1]=\frac{{\mathbf{y}}_i[p+1]+{\mathbf{y}}_j[p+1]}{2}+\frac{{\mathbf{\tau }}_{ij}\left[p\right]+{\mathbf{\tau }}_{ji}\left[p\right]}{2\nu }.$$

  (22)
• Lagrange multiplier update:

  $${\mathbf{\tau }}_{ij}[p+1]={\mathbf{\tau }}_{ij}\left[p\right]+\nu ({\mathbf{y}}_i[p+1]-{\mathbf{z}}_{ij}[p+1])$$

  (23)

Herein, p is a iteration of the ADMM optimization process, then this output of this process is ${\mathbf{y}}_i[k+1]={\mathbf{y}}_i^{*}[p+1]$.

The distributed algorithm is illustrated in Algorithm 1.

The convergence analysis of the ALADIN method has been studied clearly for both non-convex and convex optimization problem in [25]. Lemma 3 in [25] has been proven that with the cost function $f\left(\mathit{\alpha }\right)=\frac{1}{2}{\sum }_{t=1}^h{\sum }_{i\in V}{\sum }_{j\in {\mathcal{N}}_i}{({\alpha }_{t,j}-{\alpha }_{t,i})}^2$ being twice continuously differentiable and letting $({\mathit{\alpha }}_t^{*},{\mathit{\lambda }}^{*})$ for $t=1,\dots ,h$ of problem (5) be a regular KKT point. On the other hand, the Hessian ${\mathbf{B}}_t=\frac{{\partial }^2}{\partial {\mathit{\alpha }}_t^2}\left(\frac{1}{2}{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t\right)+\rho {\Sigma }_t\succ 0$ obviously, since ${\Sigma }_t⪰0$. There exists constants ${\chi }_1,{\chi }_2$ such that for every point ${\mathit{\alpha }}_t,\lambda$ satisfying the condition convergence of the decoupled minimization problem (7) have unique locally minimizers $\{{\mathbf{y}}_t,t=1,\dots ,h\}$ that satisfy $\parallel {\mathbf{y}}_t-{\mathit{\alpha }}_t\parallel \le {\chi }_1\parallel {\mathit{\alpha }}_t-{\mathit{\alpha }}_t^{*}\parallel +{\chi }_2\parallel \mathit{\lambda }-{\mathit{\lambda }}^{*}\parallel$.

Moreover, if $\frac{1}{\mu }<\mathbf{0}(\parallel {\mathbf{y}}_t-{\mathit{\alpha }}_t\parallel )$ when solving the QP (11), with the ($({\mathit{\alpha }}_t^{*},{\mathit{\lambda }}^{*})$) is a regular KKT point, then ${\chi }_1\parallel {\mathit{\alpha }}_t-{\mathit{\alpha }}_t^{*}\parallel +{\chi }_2\parallel \mathit{\lambda }-{\mathit{\lambda }}^{*}\parallel \le \frac{({\chi }_1+{\chi }_2)\omega }{2}({\chi }_1\parallel {\mathit{\alpha }}_t-{\mathit{\alpha }}_t^{*}\parallel +{\chi }_2\parallel \mathit{\lambda }-{\mathit{\lambda }}^{*}{\parallel )}^2$. This is sufficient to prove local quadratic convergence of the algorithm as ${\chi }_1,{\chi }_2$ are strictly positive constants. As a result, it is effectively applied to our proposed method since it is obviously an equality constrained non-convex optimization problem.

One thing to remark here is that in our study, the penalty parameters $\rho ,\mu$ can be updated using the following rules:

$$\begin{array}{cc}\rho [k+1]& =\{\begin{array}{cc}\hfill {\iota }_{\rho }\rho \left[k\right]& \mathrm{if}\rho \left[k\right]<{\rho }_{max}\\ \hfill \rho \left[k\right]& \mathrm{elsewhere}\end{array}\end{array}$$

$$\begin{array}{cc}\mu [k+1]& =\{\begin{array}{cc}\hfill {\iota }_{\mu }\mu \left[k\right]& \mathrm{if}\mu \left[k\right]<{\mu }_{max}\\ \hfill \mu \left[k\right]& \mathrm{elsewhere}\end{array}\end{array}$$

where ${\iota }_{\rho },{\iota }_{\mu }>1$ and ${\rho }_{max}=50,{\mu }_{max}=30$ obtained by experience to avoid the numerical problem. Moreover, we can use the blockwise and damped Broyden–Fletcher–Goldfarb–Shanno (BFGS) update, which ensures positive definiteness of the ${\mathbf{B}}_t\left[k\right]$ to preserve the convergence properties of ALADIN proposed in [25].

3.2. Retrieving the Non-Zero Laplacian Eigenvalues

As stated before, the set of eigenvalues deriving from the Algorithm 1, denoted ${\mathcal{S}}_1$, composing of the set of non-zero distinct Laplacian eigenvalues ${\mathcal{S}}_2$.

Let ${\widehat{\overline{x}}}_i$ be the final consensus value reconstructed by ${\widehat{\overline{x}}}_i={\widehat{x}}_i\left(h\right)=\frac{1}{N}{\sum }_{i=1}^N{x}_i\left(0\right)$ and the iteration scheme of the finite-time average consensus is implemented as in (1). Following the idea of the Proposition 3 in [27], we step-by-step assume to leave one element of the ${\mathcal{S}}_1$, then, the remaining elements of this set are used to reconstruct ${\widehat{x}}_i\left(h\right)$. If ${\widehat{\overline{x}}}_i={\widehat{x}}_i\left(h\right)$ is satisfied, then the left element is not the expected Laplacian eigenvalue. Hence, we can eliminate it out of the set ${\mathcal{S}}_1$. Otherwise, the left element is one of non-zero distinct eigenvalues. We restore it in the set ${\mathcal{S}}_1$ and marked as an element in the set ${\mathcal{S}}_2$. Now, the procedure is continued with another element to the end.

The distributed non-zero Laplacian eigenvalues are described in Algorithm 2.

Algorithm 1. ALADIN-based Laplacian eigenvalues estimation

Initialization: Number of nodes N, tolerance $ϵ$, initial input-output pairs $\{x_i\left(0\right),{\overline{x}}_i,i=1,2,\dots ,N\}$, where $\overline{x}=\frac{1}{N}{\sum }_{i=1}^N{x}_i\left(0\right)$ is retrieved from a standard average consensus algorithm. Each node $i,i=1,\dots ,N$ initializes: (a) random stepsizes ${\alpha }_{t,i}\left(0\right)$ for $t=1,\dots ,h=N-1$. (b) random Lagrange multipliers ${\beta }_{t,i},{\eta }_{t,i}$, for $t=1,\dots ,h$ and ${\lambda }_i,{\lambda }_{QP,i}$. (c) Semi-positive definite matrices ${\Sigma }_t\in {\Re }^N$, learning rates b, penalty parameters $\rho ,c,\mu$. Set $k=0$; Repeat: Set $k:=k+1$, Solving decouple NLP problem (7) for $t=1,\dots ,h=N-1$: - Propagate Lagrange multipliers ${\beta }_{t,i}\left[k\right]$ for $t=h,\dots ,2$ and $i=1,\dots ,N$: (a) Set ${\delta }_{h,i}\left[k\right]={\beta }_{t,i}\left[k\right]$. (b) ${\delta }_{t-1,i}\left[k\right]={\delta }_{t,i}\left[k\right]+{\alpha }_{t,i}\left[k\right]{\sum }_{j\in {\mathcal{N}}_i}({\delta }_{t,j}\left[k\right]-{\delta }_{t,i}\left[k\right]).$ - Finite-time average Consensus steps: $$x_{t,i}\left[k\right]=x_{t-1,i}\left[k\right]+{\alpha }_{t,i}\left[k\right]\sum _{j\in {\mathcal{N}}_i}(x_{t-1,j}\left[k\right]-x_{t-1,i}\left[k\right]).$$ - Propagate the error $e_{t,i}\left[k\right]$ by setting $e_{h,i}\left[k\right]={\mathbf{x}}_{h,i}\left[k\right]-{\overline{x}}_i\left[k\right]$: $$e_{t-1,i}\left[k\right]=e_{t,i}\left[k\right]+{\alpha }_{t,i}\sum _{j\in {\mathcal{N}}_i}(e_{t,j}\left[k\right]-e_{t,i}\left[k\right]).$$ - Update ${\alpha }_{t,i}$ for $t=1,\dots ,h$: $$\begin{array}{cc}\hfill {\alpha }_{t,i}[k+1]& ={\alpha }_{t,i}[k+1]-b\sum _{j\in {\mathcal{N}}_i}({\alpha }_{t,j}\left[k\right]-{\alpha }_{t,i}\left[k\right])\hfill \\ \hfill & -b{\lambda }_i-b\rho \left[k\right]\sum _{j=1}^N{\Sigma }_t(i,j)({\alpha }_{t,j}\left[k\right]-y_{t,j}\left[k\right])\hfill \\ \hfill & +b\sum _{j\in {\mathcal{N}}_i}(x_{t-1,j}\left[k\right]-x_{t-1,i}\left[k\right]){\delta }_{t,i}\left[k\right]\hfill \\ \hfill & +b\sum _{j\in {\mathcal{N}}_i}(x_{t-1,j}\left[k\right]-x_{t-1,i}\left[k\right])e_{t,i}\left[k\right]\hfill \end{array}$$ - Update NLP Lagrange multipliers ${\beta }_{t,i}$ for $t=1,\dots ,h$ by (16). Stopping criteria: if $\parallel {\sum }_{t=1}^h({\mathit{\alpha }}_t-{\mathbf{y}}_t)\parallel \le ϵ$ and $\parallel {\sum }_{t=1}^h{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t\parallel \le ϵ$ are simultaneously satisfied, then stop the optimization procedure. Compute Gradient, Jacobian matrices, and Hessian matrices as in (8)–(10), respectively. Solving the coupling QP problem (11) via solving the linear Equation (19) to define $\Delta {\alpha }_{t,i}^{*}\left[k\right],{\lambda }_{QP,i}^{*}\left[k\right]$. Update $\lambda ,{\widehat{y}}_{t,i}$: (a) ${\lambda }_i[k+1]={\lambda }_{QP,i}^{*}\left[k\right]$ (b) ${\widehat{y}}_{t,i}\left[k\right]={\alpha }_{t,i}[k+1]+\Delta {\alpha }_{t,i}^{*}\left[k\right]$ (c) Projecting ${\widehat{y}}_{t,i}\left[k\right]$ onto the constraint $\mathcal{C}$ by solving an ADMM-based optimization subproblem (20) to derive $y_{t,i}[k+1]$

Algorithm 2. Non-zero eigenvalues Estimation

Input: set of stepsizes, obtaining form Algorithm 1 ${\mathcal{S}}_1$, the input-output pairs $\{{\mathbf{x}}_i\left(0\right),{\overline{\mathbf{x}}}_i\},i=1,\dots ,N$, the threshold $ϵ$. Set ${\mathcal{S}}_2=\varnothing$, $\mathcal{S}={\mathcal{S}}_1$. Repeat: While $\mathcal{S}\ne \varnothing$, pick an element ${\alpha }_t$ out of $\mathcal{S}$, hence $\mathcal{S}\backslash \left\{{\alpha }_t\right\}$. Construct average consensus iteration scheme from the remain stepsizes to determine ${\widehat{x}}_i\left(h\right),i=1,\dots ,N$ as follows: $$x_{t,i}=x_{t-1,i}+{\alpha }_{t,i}\sum _{j\in {\mathcal{N}}_i}(x_{t-1,j}-x_{t-1,i}),t=1,\dots ,h$$ If ${\overline{x}}_i-{\widehat{\overline{x}}}_i\le ϵ$, then $\mathcal{S}=\mathcal{S}\backslash {\alpha }_t$ and return to (3) If ${\overline{x}}_i-{\widehat{\overline{x}}}_i>ϵ$, then ${\alpha }_t$ is included in ${\mathcal{S}}_2$ and $\mathcal{S}=\mathcal{S}\cup {\alpha }_t$ and return to (3) Output: If $\mathcal{S}$ is empty, then the non-zero distinct Laplacian eigenvalues can be derived by taking the inverse of the set ${\mathcal{S}}_2$’s elements.

3.3. Multiplicities Estimation

Now, turning to the last stage, which is the corresponding Laplacian eigenvalues multiplicities estimation. In [16], the authors have proposed a linear integer programming optimization problem to figure out the multiplicities:

Proposition 1

([16]). Consider a connected undirected graph of N vertices with degree sequence $\left\{d_i\right\}$ and Laplacian matrix $\mathit{L}$, having $D=|{\mathcal{S}}_2|$ distinct non-zero Laplacian eigenvalues ${\mathcal{S}}_2$. Let $\mathit{m}\in {\mathbb{Z}}^{{+}^{D\times 1}}$ be the vector of the corresponding multiplicities and be obtained by solving the integer programming below:

$$\begin{array}{cc}\hfill \underset{\mathit{m}\in {\mathbb{Z}}^{{+}^{D\times 1}}}{min}& {\mathcal{S}}_2^T\mathit{m}\\ \hfill \mathit{subject}\mathit{to}& {\mathcal{S}}_2^T\mathit{m}={\sum }_{i=1}^D{d}_i\\ & {\mathbf{1}}^T\mathit{m}=N-1\\ & \mathit{m}\in {{\mathbb{Z}}^{+}}^{D\times 1}\end{array}$$

(24)

The proof was given in [16]. Since all the multiplicities $\mathbf{m}$ are positive integers, then Brand-and-Bound method has been deployed to derive $\mathbf{m}$. Therefore, the problem (24) can be rewritten equivalently in linear integer programming form as follows:

$$\begin{array}{cc}\hfill \underset{\mathbf{m}\in {\mathbb{Z}}^{{+}^{D\times 1}}}{min}& {\mathcal{S}}_2^T\mathbf{m}-{\sum }_{i=1}^D{d}_i\\ \hfill \mathrm{subject}\mathrm{to}& {\mathbf{1}}^T\mathbf{m}=N-1\\ & \mathbf{m}\in {{\mathbb{Z}}^{+}}^{D\times 1}\end{array}$$

(25)

The Algorithm for this problem has been described clearly in [16]. At this step, the whole Laplacian spectrum $sp\left(\mathbf{L}\right)$ has been obtained.

In fact, the estimation problem can be converted into a convex form and can be solved effectively by using ADMM-based method proposed in [16]. However, the purpose of this study is extremely appropriate for non-convex optimization problem through deploying the promising ALADIN-based method.

4. Simulation Results

In this section, the efficiency of the proposed ALADIN-based method to estimate the Laplacian spectrum is evaluated by considering the two following case studies.

Firstly, it can be said that the Laplacian spectrum can decide the performance of the network since it reveals the connection of the network. For example, the robustness of the network can be estimated before starting the operations to avoid the interruption during these operations and, as a result, enhance the benefits technically and economically.

On the purpose of monitoring the connection of a large network $G^{*}(V^{*},E^{*})$, one may face with the numerical issue in step 3 of the ALADIN-based method to solve the linear system (19) due to the huge dimension of the obtained matrix that leads to the common ill-conditioning problem with the inverse matrix calculation.

In [16], the authors have suggested to partition the large network into M disjoint sub-networks $U_{\ell }$, $\ell =1,2,\cdots ,M$, [29]. Let us define ${\mathcal{N}}_i^{*}=\{j\in V^{*}:(i,j)\in E^{*}\}$ and its cardinality $|{\mathcal{N}}_i^{*}|$ as the set of neighbors of node i and its degree in $G^{*}$, respectively. Each sub-network is monitored by a super-node i which knows the number $N_{\ell }$ of agents in the sub-network and the associated average information state $x_{\ell }$. Node $i\in U_{\ell }$ is a super-node if it has at least one neighbor in a different subset, i.e., $\exists {\ell }^{*}\ne \ell$ s.t. ${\mathcal{N}}_i\bigcap U_{{\ell }^{*}}\ne \varnothing$. Let us consider that two sub-networks are connected if there exist edges linking at least two agents of these sub-networks. If two sub-networks are connected then their super-nodes are linked as showed in Figure 2. Here, the network can be social network, power system network, molecular network, etc.

Let $G=(V,E)$ be the undirected graph representing a network with $N=\left|V\right|$ super-nodes, which are black nodes in Figure 2. G captures the interaction between sub-networks of $G^{*}$. Therefore, the large network is to be robust if the partitions are strongly linked to each other and if the critical threshold is high enough in [16]. Then, the Laplacian spectrum of network of super-nodes G can monitor the connection of the large network $G^{*}$ via the robustness index.

4.1. Case Study 1

Let us consider a large network partitioning into 4 disjoint sub-networks. Each sub-network has only one super-node. These super-nodes interact with each other by the graph $G=(V,E)$, depicted in Figure 3.

This network does have the Laplacian eigenvalues $sp\left(\mathbf{L}\right)=\{0,4,4,4\}$. With the parameter of each super-node, denoted as $\mathbf{x}\left(0\right)=\{0.7417,0.7699,0.3216,0.5466\}$, after deploying Algorithm 1, the set of ${\mathit{\alpha }}_t,t=1,\dots ,3$ is obtained. The nodes trajectories are described as in Figure 4.

As can be seen, all ${\mathit{\alpha }}_t$ execute the consensus problem at the beginning of the procedure, and then dig into satisfying the constraint.

Figure 5 illustrate the convergence of the cost function $\frac{1}{2}{\sum }_{t=1}^{N-1}{\mathit{\alpha }}_t^T\mathbf{L}{\mathit{\alpha }}_t$ in according to the constraint $\mathbf{x}(N-1)=\overline{x}\mathbf{1}$.

Furthermore, the Algorithm 2 is applied to eliminate the unexpected eigenvalues. As a result, we receive only one eigenvalue $\lambda =\frac{1}{0.25}=4$. In order to accomplish the Laplacian spectrum estimation’s procedure, we make use the Proposition 1 to pick out $m=3$. Now, the entire Laplacian spectrum is achieved.

At this step, the robustness index such as Kirchhoff index or the number of spanning trees can be calculated.

$${\mathcal{R}}_{\mathbf{L}}=N\sum _{i=2}^{D+1}\frac{m_i}{{\lambda }_i}=4\frac{3}{4}=3.$$

Now, let us see how the proposed procedure works via the table below.

Table 1. The achievement of the proposed procedure in the sense of 4-node topology.

${\mathit{S}}_1$	$\{0.25,0.5919,0.5543\}$
${\mathit{S}}_2$	$\left\{0.25\right\}$
$\mathbf{m}$	3
Iterations	130
${\mathcal{R}}_{\mathbf{L}}$	3

is obtained after executing the proposed procedure.

As can be seen in Figure 4, at the iteration of around 90 the procedure can be stopped. In order to access the robustness of the whole large network, it is necessary to define the critical threshold, introduced in [16].

Next, we implement a comparison with the Lagrange method described in [1] by considering Case study 2.

4.2. Case Study 2

Let us consider a 6-node network described in Figure 6.

It is known that for this topology, the Laplacian spectrum is $sp\left(\mathbf{L}\right)=\{0,1,2,3,3,5\}$.

Let us define the same initial information state of each node at time $t=0$: $\mathbf{x}\left(0\right)=\{0.5832,0.74,0.2348,0.7350,0.9706,0.8669\}$ for both methods. By using a standard consensus algorithm [26], the consensus value $\overline{x}=0.6884$ can be easily inferred.

It can be seen in Figure 7 that the Lagrange-based method in [1] takes a long time to achieve the consensus first and then track the constraint to obtain the stepsizes ${\alpha }_t,t=1,\dots ,N-1$.

Now, with the same initial ${\mathit{\alpha }}_t\left(0\right),t=1,\dots ,5$, for the iterative procedure of the proposed method, the ${\mathit{\alpha }}_t$ after using the Algorithm 1 with the convergence trajectories of the ${\mathit{\alpha }}_t$ are illustrated in Figure 8.

Figure 7 and Figure 8 express the significant pros of our proposed method since the number of iterations is much less than that of Lagrange-based method. The consensus term is executed in advance from the start of the procedure and then try to reach the constraint term to figure out the expected values of the stepsizes ${\mathit{\alpha }}_t=\{1,0.5,0.1027,0.2,0.3333\}$. Obviously, since the authors operate the proposed algorithm with the number of ${\mathit{\alpha }}_t$ being $N-1=5$, it is needed to run the next stage to eliminate the residual values. As can be seen clearly, the executive time for ALADIN-based method is significantly faster than the proposed method as showed in Figure 9.

Figure 9 shows that the ALADIN-based method approaches the destined values earlier than Lagrange-based method. Furthermore, in order to get the non-zero Laplacian eigenvalues of the given network, the Algorithm 2 is carried out to obtain vector of stepsizes ${\mathit{\alpha }}_t=\{1,0.5,0.2,0.3333\}$.

Finally, by constructing the Brand-and-Bound based method to solve the Problem 1, which has been proposed in [16], we achieve the vector of multiplicities $\mathbf{m}=\{1,1,1,2\}$, hence deduce the Laplacian spectrum $sp\left(\mathbf{L}\right)=\{1,2,3,3,5\}$.

Now, let us see how the proposed procedure works via the table below.

Table 2. The achievement of two methods in the sense of 6-node topology.

	ALADIN-Based Method	Lagrange-Based Method
${\mathit{S}}_1$	$\{1,0.5,0.1027,0.2,0.3333\}$	$\{0.9992,0.2,0.5001,0.0895,0.3333\}$
${\mathit{S}}_2$	$\{1,0.5,0.2,0.3333\}$	$\{0.9992,0.2,0.5001,0.3333\}$
$\mathbf{m}$	$1,1,1,2$
Iterations	33950	372800
${\mathcal{R}}_{\mathbf{L}}$	14.2

is obtained after executing the proposed procedure and the Lagrange-based method.

Notice that the proposed method gives results much better than the Lagrange-based method. Let us see that at the iteration of 33,950, the Lagrange-based method gives the set of ${\mathcal{S}}_1$ that has still not satisfied the constraint.

Recently, besides the Laplacian spectrum estimation basing on the optimization approaches, there are also some works that approximate the Laplacian spectrum via iterative dynamics process (taking random walk for an example). However, from our point of view, the dominant contribution of our proposed method is the possibility of implementation in a decentralized manner. Moreover, another method to faster the convergence rate of the Laplacian spectrum estimation procedure is the ADMM-based method, proposed in [16]. The important step in this work is to re-parameterize adequately the non-convex formulation into convex one. It is hard to compare the efficiency between ADMM-based method and the proposed method. Since, our study focuses on the non-convex formulation.

5. Conclusions

In this paper, the authors have proposed a promising ALADIN-based method to find out the Laplacian spectrum of a given dynamic network in a distributed way. First and foremost, the study has assumed that the number of agents N in the network can be accumulated by random walk algorithm constructed in the configuration step. Briefly speaking, the proposed procedure is divided into 3 stages. The first stage is to determine the $N-1$ Laplacian eigenvalues. Then, retrieve the non-zero distinct Laplacian eigenvalues in stage 2 before estimating the corresponding multiplicities in stage 3. The ALADIN-based method is appropriate for carrying out the factorization of the average consensus matrix as factors of the Laplacian-based consensus matrices to minimize the disagreement between neighbors on the values of ${\alpha }_{t,i}$. Then, the Laplacian eigenvalues can be defined by taking the inverse of these ${\mathit{\alpha }}_t$. Herein, the authors are obviously interested in dealing with the non-convex optimization problems. From the simulation evaluation, it can be concluded that the proposed method converges much faster in comparison with the gradient descent method in [1] for the estimation of Laplacian spectrum.