There are a lot of complex networks in the real world, such as the power grids [1], the World Wide Web [2], gene regulatory networks [3] in systems biology, just name a few. In general, if dynamics and topological structure of a network are known, there are lots of studies on its dynamical properties [4-8] and the control of the network [9-13]. However, in many real applications, topological structure of a complex network is usually unknown and one needs to reconstruct it from the observed data.

In this article, we will consider a dynamical network consisting of $N$ nodes as follows:

$\dot {\mathbf{x}}_i(t)=\mathbf{f}_i (\mathbf{x}_i(t), t)+\sum\limits _{j=1}^N c_{ij}\mathbf{x}_j(t), \quad i =1, 2, \cdots, N,$

(1.1)

where each $\mathbf{x}_i(t)=(x_{i1}(t), x_{i2}(t), \cdots, x_{in}(t))^\top$ is an $n$-dimensional state vector corresponding to the $i$-th node, and $\mathbf{f}_i:R^n \times R^+ \mapsto R^n$ is dynamical system of the $i$-th node. The coupling matrix $C=(c_{ij})_{N \times N}$ describes the connections between nodes. If there is a connection between node $j$ and node $i (j\neq i)$, then let $c_{ij}\neq 0$; otherwise $c_{ij}=0$. Its magnitude $|c_{ij}|$ represents the strength of the connection. The coupling matrix $C$ is unknown and one needs to derive it from the observed data.

There are many different techniques to reveal the topological structure of a given network. For example, in [14-19], people use synchronization strategy to derive the structure of a given network, methods based on time series can be found in [20-24] and optimization methods were given in [25, 26]. In addition, the perturbation method and the steady-state control method were used in [27] and [28]. Very recently, methods based on compressive sensing were applied to reconstruct networks of sparse structure [29-31].

If the evolution of nodes (i.e., $\{\mathbf{x}_i(t)\}_{i=1}^N$) can be observed from experiments, to derive the network's structure (coupling matrix $C$) from evolution of the nodes is a typical parameter identification problem. Tikhonov regularization method was widely applied in inverse problem community [32-34], and this method can be very efficient if the regularization term is chosen properly. Since the real world networks often have sparse structures, we can make the sparsity of coupling matrix $C$ as the priori information. Stimulated by the priori sparsity information, we use the sparsity-promoting $\ell^1$-regularized technique in this work. In recent years, sparse regularization method [35, 36] attracted a lot of attention in inverse problem community. In cases of linear problems with $\ell^1$ regularization, there are many efficient algorithms, such as iterative soft shrinkage method [37], split Bregman iteration method [38], alternative direction method of multiplier (ADMM) [39], iteratively reweighted least squares (IRLS) method [40], semi-smooth Newton method [41, 47] and etc. In our case, estimation of matrix $C$ is a non-linear problem and our objective function includes two parts. The first data fitting part is smooth but non-convex and the second $\ell^1$ regularization term is non-smooth but convex. To overcome the non-smoothness arising from $\ell^1$ regularization term, we will apply IRLS method [42] to change the original problem into a series of smooth optimization subproblems. Then Newton method was applied to solve subproblems with a warm start technique [48].

The rest of this paper is organized as follows. In Section 2, we introduce the mathematical formula of the $\ell^1$-regularized problem. The detailed algorithm is also given in this section. Several numerical examples are given in Section 3 to show the efficiency of the method.

Recall the nodal dynamics satisfying equation (1.1). For simplicity, we rewrite the dynamical network (1.1) in a compact form as in [26]. Let

$X=(x_{11}, x_{12}, \cdots, x_{1n}, \cdots, x_{N1}, x_{N2}, \cdots, x_{Nn})^\top \in \mathbb{R}^{Nn}$

and

$F =(f_{11}, f_{12}, \cdots, f_{1n}, \cdots, f_{N1}, f_{N2}, \cdots, f_{Nn})^\top \in \mathbb{R}^{Nn}, $

then network model can be rewritten as follows:

$\dot{X} = F(X)+C\otimes I_{n \times n} \cdot X,$

(2.1)

where $\otimes$ is a Kronecker product. Our aim is to estimate the coupling matrix $C$ from the information of state vector $X(t)$.

Equipped with an initial data $X_0$, we can get the state vector $X(t)$ by equation (2.1). Let the observation data be $X_u(t)$. To deduce $C$ from $X(t)$ is a typical inverse problem and the regularization technique is usually needed. Define a map $M:R^{N \times N} \mapsto L^2(R^+, R^{Nn})$ as $M:C \mapsto X(t)$. The identification of $C$ can be formulated as an output least squares form:

$\min\limits_{C \in \mathcal{B}} J_0(C)=\dfrac{1}{2} \int_0^T \|M(C)-X_u\|^2 dt,$

(2.2)

where $\|.\|$ is the Euclidean norm and $\mathcal{B}$ is the set of all possible coupling matrixes,

$\mathcal{B}=\Big\{C \in \mathbb{R}^{N \times N} \Big|\sum\limits_{j=1}^N c_{ij} = 0, \, \, i=1, 2, \cdots, N \Big\}.$

(2.3)

In order to identify different kinds of coupling matrixes, choosing a proper regularization term is very important and it plays a crucial role for the successful identification. In this work, we focus on the sparse networks. Compared $\ell^1$-regularized term with $\ell^2$-regularized term, in terms of promoting the sparsity of solution, $\ell^1$-regularized technique should be better than $\ell^2$-regularized technique [43]. Inspired by this information, we apply the $\ell^1$-regularized technique as follows:

$\min\limits_{C \in \mathcal{B}} J_\alpha(C)=\dfrac{1}{2} \int_0^T \|M(C)-X_u\|^2 dt +\alpha\sum\limits_{i, j=1}^N|c_{ij}|.$

(2.4)

Unlike the standard Tikhonov $\ell^2$-regularization, $J_\alpha(C)$ is a non-smooth function and this causes a lot of difficulties to find the optimal solution. It is observed that the objective function $J_\alpha(C)$ includes two parts: the first term $\dfrac{1}{2} \int_0^T\|M(C)-X_u\|^2 dt$ is smooth but non-convex and the second term $\alpha \sum\limits_{i, j=1}^N|c_{ij}|$ is convex but non-smooth. To overcome the non-smoothness of the second term, we will apply IRLS algorithm. Firstly, the non-smooth term is approximated by

$\large |c_{ij}|=\frac{c^2_{ij}}{|c_{ij}|} \approx \frac{c^2_{ij}}{\sqrt{|c_{ij}|^2+\epsilon^2}}.$

(2.5)

Then we apply a fixed point iteration to solve a series of subproblems as follows:

$\min\limits_{C \in \mathcal{B}} J^k_\alpha(C) = J_0(C)+ \alpha \sum\limits_{i, j=1}^N\frac{c^2_{ij}}{\sqrt{|c^{k-1}_{ij}|^2+\epsilon^2_k}},$

(2.6)

where $c^{k-1}_{ij}$ is the solution from the $(k-1)^{th}$ iteration, and $\epsilon_k$ is the cutting weight factor.

To solve subproblem (2.6), one needs to find its first order optimality system. By applying the similar derivation as in [26], we can easily get the first order necessary optimality (KKT) system of problem (2.6) as follows:

Theorem 2.1 The first-order optimality system (KKT condition) of problem (2.6)

$\left\{\begin{array}{ll} \mbox{primal equation}&\dot{X}=F(X)+C\otimes I_{n \times n} \cdot X, \, \, X(0)=X_0, \\[1.5ex] \mbox{adjoint equation}&-\dot{P}=(\mathrm{diag}(\nabla\mathbf{f}_1, \cdots, \nabla \mathbf{f}_N)+C^\top \otimes I_{n \times n})P-(X-X_u), \, \, P(T)=0, \\[1.5ex] \mbox{optimality condition}&\alpha\frac{c^2_{ij}}{D_{ij}}-\displaystyle\int^T_0 \mathbf{p}_i\cdot \mathbf{x}_j dt+\lambda_i=0, \, \, i, j=1, 2, \cdots, N, \\[1.5ex] \textrm{admissible condition}&\sum\limits^N_{j=1}c_{ij} = 0, \, \, i=1, 2, \cdots, N, \end{array}\right.$

(2.7)

where $D_{ij}$ is $\sqrt{|c^{k-1}_{ij}|^2+\epsilon^2_k}$.

Because optimality condition and admissible condition are not ordinary differential equations, we couldn't solve the KKT system above using ordinary differential solver directly. Here, we discretize the primal equation and the adjoint equation firstly. Then we solve the KKT system from the perspective of nonlinear equations. We will apply the Newton iteration method, since it has the locally quadratical convergence property.

For simplicity, we discretize optimality condition and adjoint equation by forward Euler scheme. Let $h$ be the time step, $T$ be the total time and $M=\tfrac{T}{h}$. With the time grids $t_i=ih$ $(i=0, 1, 2, \cdots, M), $ the discrete KKT system reads as follows $(m=0, 1, 2, \cdots, M-1)$:

$\left\{\begin{array}{ll} X^{m+1}-X^{m}-h(F(X^{m})+C\otimes I_{n \times n}\cdot X^{m})=0, \quad X^0 = X_0, \\[1ex] P^{m}-P^{m+1}-h((\mathrm{diag}(\nabla\mathbf{f}_1(x_1^m), \cdots, \nabla \mathbf{f}_N(x_N^m))+C^\top \otimes I_{n \times n})P^{m}\\-(X^{m}-X_u^{m}))=0, \quad P^M = 0, \\[1ex] \alpha\frac{c_{ij}^2}{D_{ij}}-\sum\limits_{m=0}^{M-1} \mathbf{p}_i^{m}\cdot \mathbf{x}_j^{m}h+\lambda_i=0, \, \, i, j=1, 2, \cdots, N, \\[1ex] \sum^N_{j=1}c_{ij} = 0, \, \, i=1, 2, \cdots, N. \end{array}\right.$

(2.8)

Since the coupling matrix $C$ is a sparse matrix, only a little information from $X(t)$ is needed to recover $C$, i.e., $M$ is not necessary to be large. On the other hand, in order to apply Newton method, a good initial guess is needed. We have a nature choice for the initial guess, i.e., the solution of the $(k-1)^{th}$ subproblem is a fairly good initial guess for the $k^{th}$ subproblem.

In this subsection, we give a complete algorithm to estimate the coupling matrix. The details of the algorithm can be found in Algorithm 1.

Since the regularization parameter $\alpha$ is important for the estimation, we choose a continuation technique for regularization parameter, and a discrepancy principle is applied to choose a proper $\alpha$. Then an iteratively reweighted least squares method is proposed to solve (2.4) with a given parameter $\alpha$.

The choices of $\alpha_0$ and $\epsilon_k$ are not sensitive. We can choose any $\alpha_0 = O(1)$ and $\epsilon_k = O(1)$. Since Newton method has a locally quadratic convergence, it is very efficient when a good initial guess is equipped. The inner loop in line 5 just requires a few Newton steps (typically two or three). In general, larger $\rho$ and larger $R$ can obtain more accurate and stable solutions. But, this will cost more CPU time. The IRLS algorithm is a fixed point type algorithm and has slow convergence, but we don't need to find the very accurate solution for $\alpha_k$-subproblem. A few IRLS steps will give a fairly good initial guess to the next $\alpha_{k+1}$-subproblem. In the numerical test, we choose $\rho = 0.7$ and $R=5$.

In each IRLS step, we need to solve the discrete non-linear equations (2.8). If we choose a small time step $h$, i.e., more observed information, one needs to solve a non-linear system of large size. Newton method will cost much due to a large number of the unknown variables. However, because of the sparse structure of the coupling matrix $C$, the degree of freedom of coupling matrix $C$ is not large, then a few observed information can produce a good estimation effect. The choice of $M$ depends on the degree of freedom of the sparse coupling matrix which we will inverse.

In this section we present several numerical examples. Assume the nodal dynamics to be Lorenz system [44] which is described by $f_i(x_i)=[-ax_{i1}+ax_{i2};cx_{i1}-x_{i2}-x_{i1}x_{i3};-bx_{i3}+x_{i1}x_{i2}]$ with parameters $a=10, b=8/3, c=28$ and $i=1, 2, \cdots, N$. During the numerical experiments, let $T=0.1$ and $M=10$.

Example 1 In this experiment we estimate a small-world symmetric network [45] with $30$ nodes. Figure 2 displays the network's topological structure with weighted connections and Figure 1 displays the colormap corresponding to the coupling matrix whose true structure is C^†. Firstly, let the observed data be free of noise. Figure 3 shows the estimated structure $C$ and Figure 4 shows error of $C$ and C^†. Then, add $1\%$ noise. The estimated coupling matrix and the error one are displayed in Figure 5 and Figure 6. The reconstruction is quite good.

Example 2 In the second example we consider a BA scale-free symmetric network [46] with 30 nodes. Figure 8 and Figure 7 respectively show its topological connections with different weights and the colormap of the corresponding matrix whose true structure is C^†. Figure 9 and Figure 10 show the noise free case. Figure 11 and Figure 12 show the noisy case.

Example 3 The last example is a $80$-nodes small-world network. In coupling matrix C^†, if there is a connection from node $i$ to node $j$($i\neq j$), then let C_ij^†=1; otherwise C_ij^†=0. The topological structure and the coupling matrix with the true structure C^† are showed in Figure 13 and Figure 14. The noise free case is shown in Figure 15 and Figure 16.

In this paper we propose a sparsity-promoting $\ell^1$-regularized output least squares method to identify the large-scale sparse complex networks. The optimization problem is solved by the iteratively reweighted least squares (IRLS) method. The Newton algorithm with continuation strategy is given to solve each subproblem of the inner loop of IRLS. For the sparsity of coupling matrix, a satisfactory estimation can be obtained by only a few observed data. Several numerical examples are presented to verify the efficiency of the method.