Recovering network topologies via Taylor expansion and compressive sensing Guangjun Li, Xiaoqun Wu, Juan Liu, Jun-an Lu, and Chi Guo Citation: Chaos 25, 043102 (2015); doi: 10.1063/1.4916788 View online: http://dx.doi.org/10.1063/1.4916788 View Table of Contents: http://scitation.aip.org/content/aip/journal/chaos/25/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in A new protocol for finite-time consensus of detail-balanced multi-agent networks Chaos 22, 043134 (2012); 10.1063/1.4768662 Announcement: Focus Issue on “Mesoscales in Complex Networks” Chaos 20, 010202 (2010); 10.1063/1.3298887 Stochastic bimodalities in deterministically monostable reversible chemical networks due to network topology reduction J. Chem. Phys. 131, 195103 (2009); 10.1063/1.3264948 Correlations between structure and random walk dynamics in directed complex networks Appl. Phys. Lett. 91, 054107 (2007); 10.1063/1.2766683 Topology for Dominance for Network of Multi‐Agent System AIP Conf. Proc. 913, 96 (2007); 10.1063/1.2746731

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

CHAOS 25, 043102 (2015)

Recovering network topologies via Taylor expansion and compressive sensing Guangjun Li,1 Xiaoqun Wu,2,a) Juan Liu,1,a) Jun-an Lu,2 and Chi Guo3 1

Computer School, Wuhan University, Hubei 430072, China School of Mathematics and Statistics, Wuhan University, Hubei 430072, China 3 Global Navigation Satellite System Research Center, Wuhan University, Hubei 430072, China 2

(Received 2 July 2014; accepted 23 March 2015; published online 6 April 2015) Gaining knowledge of the intrinsic topology of a complex dynamical network is the precondition to understand its evolutionary mechanisms and to control its dynamical and functional behaviors. In this article, a general framework is developed to recover topologies of complex networks with completely unknown node dynamics based on Taylor expansion and compressive sensing. Numerical simulations illustrate the feasibility and effectiveness of the proposed method. Moreover, this method is found to have good robustness to weak stochastic perturbations. Finally, the impact of two major factors on the topology identification performance is evaluated. This method provides a natural and direct point to reconstruct network topologies from measurable data, which is likely to C 2015 AIP Publishing LLC. have potential applicability in a wide range of fields. V [http://dx.doi.org/10.1063/1.4916788]

Recovering network topologies from measurable data is a fundamental problem for understanding and controlling collective dynamics of complex systems. In the past few years, various methods have been developed for topology identification, among which is the famous compressive sensing method. The main feature of compressive sensing is that it only needs a small amount of measured data to infer the topological structures of complex networks. However, this method is limited to those complex networks with previously known basis functions of node dynamics, such as polynomials. Therefore, it does not work if the basis functions of network node dynamics are unknown, whereas this is a practically common situation. Inspired by the idea of polynomial approximation of nonlinear functions provided by Taylor expansion, a new method based on Taylor expansion and compressive sensing is proposed to overcome this difficulty.

I. INTRODUCTION

Complex dynamical networks widely exist, and gaining knowledge of underlying structures can help us to understand, predict, and control the behaviors of these networks. For example, for synchronized networked systems which are harmful, one can regulate the structures in order to avoid synchronization. For vibrating dynamical networked systems with a precisely known topological structure which behave in an unstable state that is unexpected, one can control some of the network nodes based on the underlying structure so that the network can arrive at some stable behavior.1–4 However, in many practical dynamical networks, the topological structures are often unknown or unobservable. Therefore, recovering the underlying topological structures a)

Authors to whom correspondence should be addressed. Electronic addresses: [email protected] and [email protected].

1054-1500/2015/25(4)/043102/9/$30.00

of complex networks is of theoretical significance and practical value. Due to the importance of inferring network topologies, there have emerged various approaches in recent years. For example, Pompe and Frenzel detected links from purely observed data by information theoretic based approaches, such as the techniques based on partial mutual information,5 conditional mutual sorting information,6 and inner composition alignment,7 Yu et al. identified network topologies based on adaptive control and outer synchronization between two networks,8–12 Jansen et al. by inferring Bayesian networks,13,14 Marwan et al. by recurrence,15,16 Wang et al. by noisy scaling17 and dynamical correlations of noisy time series,18,19 Wu et al. by Granger causality test,20,21 He et al. by optimization algorithm based on the projected conjugate gradient method,22 and Wang et al. by compressive sensing.23–28 It is worth mentioning that in Ref. 28, Shen et al. developed a framework based on compressive sensing to reconstruct complex networks on which stochastic spreading dynamics take place and established a paradigm for tracing and controlling epidemic invasion and information diffusion in complex networked systems. As a new type of sampling theory, the great success of compressive sensing relies essentially on the ability to efficiently predict an approximately sparse solution to an underdetermined linear system. Although the term compressive sensing was coined only recently with the paper by Donoho,29 followed by vast research activity,23–28,30–34 such a development did not start out of thin air. It has many potential applications in signal processing, medical imaging, computer science as well as pure mathematics. As a main feature, compressive sensing only needs a small number of measured data to estimate multiple unknown parameters. For complex networks with known basis functions of node dynamics, the compressive sensing method can not only tell topological structures but also give detailed node dynamical

25, 043102-1

C 2015 AIP Publishing LLC V

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-2

Li et al.

Chaos 25, 043102 (2015)

parameters and coupling strengths between nodes. If the basis functions of node dynamical equations are unknown, the traditional compressive sensing approach no longer works. However, in many real-world networks, the basis functions of node dynamics are rarely known. A question thus arises: can we still recover network topologies based on compressive sensing in this case? The answer is yes. In Ref. 35, Napoletani and Sauer proposed to use linearization of system dynamics around the center of a neighborhood of data points and constrained optimization techniques based on the L1 vector norm to reconstruct topologies for sparsely connected networks. Simultaneously, they added extra random vectors to the basis to absorb noise and model error in the data. Compared with network reconstruction by compressive sensing,23–28 the superiority of Napoletani and Sauer method lies in its applicability to networks with unknown node dynamics. However, it is less accurate for networks with higher connectivity rates. The need is thus raised for developing a more accurate topology identification framework for either sparsely or densely connected networks in the absence of knowledge for node dynamics. Based on above discussions, a new method to recover topological structures of complex networks by combining different orders of Taylor expansion and compressive sensing is proposed, where the node dynamics as well as the basis functions are completely unknown. Numerical simulations are provided to illustrate the effectiveness of the proposed method. Moreover, it is shown that the method is robust to weak stochastic perturbations. The rest of the paper is organized as follows. In Sec. II, the network model and the new method are presented. In Sec. III, numerical examples are provided to illustrate the effectiveness and evaluate the performance of the proposed method. Finally, some conclusions are given in Sec. IV. II. METHOD

In the real world, topologies of complex dynamical networks are often unknown, such as protein action networks and neural networks. Therefore, identifying network topologies is the most fundamental problem among many practical ones associated with complex networks. For a network, when neither the node dynamics nor the basis function is known, how to identify a network’s topological structure is proposed via the combined technique of compressive sensing and Taylor expansion. Consider a network consisting of N coupled nodes, with the dynamics of the i th node being described by x_ i ¼ f i ðt; xi Þ þ r

N X

cij Axj ;

i ¼ 1; 2; …; N;

(1)

j¼1

where xi ¼ ðxi1 ; xi2 ; …; xid Þ> is the state vector, f i ðt; xi Þ ¼ ðfi1 ; fi2 ; …; fid Þ> is a vector function governing the ith node dynamics, r is the coupling strength of the network, A 2 Rdd is the inner coupling matrix, and C ¼ ðcij ÞNN is the outer coupling matrix representing the network topology,

which is defined as follows: if node j affects node i ðj 6¼ iÞ; cij ¼ 1, otherwise, cij ¼ 0 ðj 6¼ iÞ; the diagonal eleP ments of matrix C are defined as cii ¼  Nj¼1;j6¼i cij : In network (1), only the inner coupling matrix A is supposed to be known, while the node dynamics or the types of basis functions are supposed to be unknown. Our purpose is to recover the unknown topological structure represented by C through observed data xi ðtÞ ði ¼ 1; 2; …; NÞ at M time points tk ðk ¼ 1; 2; …; MÞ. The problem of recovering topologies can be described as the estimation of an unknown vector X 2 Rv from linear measurements Y about X in the form of Y ¼ UX, where Y 2 Rw , w is the number of measurements, and U is a w  v matrix. In general, the number of collected measurements of a discrete finite-dimensional signal should be at least as large as its length in order to ensure reconstruction. However, compressive sensing—a signal processing technique for efficiently acquiring and reconstructing a signal by finding solutions to underdetermined linear systems—provides a fundamentally new approach to reconstruct a sparse signal exactly from what was previously believed to be highly incomplete measurements by solving a convex optimization program.36 Specifically, to reconstruct a sparse signal X, the number w of measurements can be much less than the number v of components of the unknown signal.23 Here, X is sparse means that most components of X are zero. The sparse vector X can be reconstructed by solving the following convex-optimization problem: min kXk1 subject to Y ¼ UX; where kXk1 ¼ Rvi¼1 jXi j is the L1 norm of the sparse vector X. Solutions to the convex-optimization problem have been worked out by Donoho, Candes, and Baraniuk et al.29,30,32,33 and improved or extended by many others.34 Particularly, according to Fornasier and Rauhut,36 under the standard restricted isometry property (RIP) assumption on the matrix U, all ksparse vector X 2 Rv (which means that the support of a vector X is less than k) can be stably recovered with high probability provided the number w of measurements w  Ck logðv=wÞ, where C is a fixed constant. The topology inference problem of the unknown network (1) can be cast in the above form, where the topological parameters cij ði; j ¼ 1; 2; …; NÞ are contained in the sparse vector X, and dynamical observations from network (1) are in Y and U. To recover the topology of network (1), the main challenge lies in the fact that node dynamical functions f i ðt; xi Þ are completely unknown. In mathematics, a Taylor series can usually represent a function as an infinite sum of terms that are calculated from the values of the function’s derivatives at a single point. The foundation of this research is the common practice to locally approximate an unknown function by using finite terms of its Taylor series. For example, if the vector function f i ðt; xi Þ is differentiable at a point xi , one has the following first-order Taylor approximation for f i ðt; xi Þ:

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-3

Li et al.

Chaos 25, 043102 (2015)

f i ðt; xi Þ  f i ðt; xi Þ þ Df i ðt; xi Þðxi  xi Þ;

(2)

where 0

@fi1 B @x B i1 B . Df i ¼ B B .. B @ @fid @xi1

1 @fi1 @xid C C C .. C: .. . C . C @fid A … @xid …

Therefore, the i-th node dynamics (1) can be approximately written into the following linear form: x_ i  f i ðt; xi Þ þ Df i ðt; xi Þðxi  xi Þ þ

N X

cij Axj :

(3)

j¼1

Since xi is a d-dimensional vector, one has x_ i ¼ ðx_ i1 x_ i2    x_ id Þ> . Denote the dynamics of the k-th (k ¼ 1; 2; …; d) component of x_ i as x_ ik ¼ /ik vik ; where /ik ¼ ½1; xi1  xi1 ; xi2  xi2 ; …; xid  xid ; x11 ; x12 ; …; x1d ; x21 ; x22 ; …; x2d ; …; xN1 ; xN2 ; …; xNd ;

(4)

and    @fik ðt; xi Þ @fik ðt; xi Þ @fik ðt; xi Þ ; ; …; ; vik ¼ fik t; xi ; @xi1 @xi2 @xid Ak1 Ci1 ; Ak2 Ci1 ; …; Akd Ci1 ; Ak1 Ci2 ; Ak2 Ci2 ; …; > (5) Akd Ci2 ; …; Ak1 CiN ; Ak2 CiN ; …; Akd CiN : For M observations, one obtains the following equation for the kth component of node i in the form of Y ¼ UX: 1 0 1 x_ ik ðt1 Þ /ik ðt1 Þ B C B C B x_ ik ðt2 Þ C B /ik ðt2 Þ C B C B C Yik ¼ B ¼B vik ¼ Uik vik : .. C .. C B B C . A @ . C @ A /ik ðtM Þ x_ ik ðtM Þ

It is noted that the measurement vector Yik and matrix Uik for node i are related to the reference point xi , and linear approximation is employed in the first-order Taylor expansion, which is basically the method proposed in Ref. 32. Therefore, it is assumed that M observations can be collected in the neighborhood of xi . Furthermore, x_ ik ðtj Þ is calculated using the first-order forward difference quotient. However, according to Napoletani and Sauer,35 the considered network has to be sparse enough so that the L1 minimization can yield an accurate solution. Moreover, since linear approximation is employed, measurements have to be available in a very small neighborhood of each reference point. In practical applications, to obtain more accurate topology identification, one had better use the second-order or even higher order Taylor expansion. High-order Taylor expansion can increase the sparsity of the unknown vector to be recovered and approximate nonlinear dynamics much better than the first-order one. Therefore, a higher-order Taylor expansion can lead to more accurate identification of network topologies. For example, suppose the node state vector xi in a complex network has two components, let ‘i ðtj Þ ¼ xi1 ðtj Þ  xi1 ; hi ðtj Þ ¼ xi2 ðtj Þ  xi2 , where tj ðj ¼ 1; 2; …; MÞ is the j-th observation time point, then the second-order Taylor expansion results in a linear equation of compressive sensing as follows: 0 1 0 1 x_ ik ðt1 Þ uik ðt1 Þ B C B C B x_ ik ðt2 Þ C B uik ðt2 Þ C C B C B Yik ¼ B ¼B #ik ¼ Wik #ik ; (7) .. C .. C B C B . A @ . C @ A x_ ik ðtM Þ uik ðtM Þ where uik ðtj Þ ¼ ½1; ‘i ðtj Þ; hi ðtj Þ; ‘i ðtj Þ2 =2; ‘i ðtj Þhi ðtj Þ; hi ðtj Þ2 =2;

0

x11 ðtj Þ; x12 ðtj Þ; x21 ðtj Þ; …; xN2 ðtj Þ; j ¼ 1; 2; …; M; (8) (6)

According to Eq. (6) and based on the compressive sensing principle, if vik is sparse, one only needs a small set of measurements (M < Nd þ d þ 1) to infer the underlying topology of a complex network, which is contained in vik ði ¼ 1; 2; …; N; k ¼ 1; 2; …; dÞ: For a complex network consisting of unknown types of node dynamics, the following steps are needed for recovering the network’s structure in our method: (i) for the i th node (i ¼ 1; 2; …; N), a reference point xi is randomly selected, and M observed data are collected in its spherical neighborhood; (ii) an optimal solution for vik in Eq. (6) is found via the compressive sensing method, with the first 1 þ d elements of vik being decided by the unknown node dynamics and the reference point, and the other elements being the topological parameters; and (iii) carry out (i) and (ii) for k ¼ 1; 2; …; d and i ¼ 1; 2; …; N:



  @fik ðt; xi Þ @fik ðt; xi Þ @ 2 fik ðt; xi Þ #ik ¼ fik t; xi ; ; ; ; @xi1 @xi2 @ 2 xi1 @ 2 fik ðt; xi Þ @ 2 fik ðt; xi Þ ; ; Ak1 Ci1 ; Ak2 Ci1 ; @xi1 xi2 @ 2 xi2 > Ak1 Ci2 ; Ak2 Ci2 ; …; Ak1 CiN ; Ak2 CiN ; k ¼ 1; 2:

(9)

A similar algorithm as that for solving the sparse vector vik in Eq. (6) based on compressive sensing can be employed to estimate #ik in Eq. (7) for i ¼ 1; 2; …; N and k ¼ 1, 2. III. NUMERICAL SIMULATIONS

In this section, complex networks with different node dynamics and topological structures are used to illustrate the feasibility and effectiveness of our method. In simulation study, the Runge-Kutta method is employed to generate time series from the complex network model (1) with an equal time step of 0.01 as observed data. Once a reference point xi for node i is selected, M measurable data of node i which lies

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-4

Li et al.

in xi ’s spherical neighborhood are collected, with the corresponding time points tk ðk ¼ 1; 2; :::; MÞ being recorded. The dynamical states of all the other nodes at the M time points are also collected as needed measurements for topology inference. Additionally, the coupling strength r of the network under study is supposed to be known and one only focuses on inferring the existence or nonexistence of edges in the network. To quantify the identification performance of the proposed method, one defines the true positive rate (TPR) as the fraction of correctly identified edges out of all existent edges, and the true negative rate (TNR) as the fraction of correctly identified non-existent edges out of all nonexistent edges. Furthermore, the prediction errors are defined separately for existent and non-existent edges in the network. The relative prediction error of an existent edge, that is, nonzero edge weight rcij ði 6¼ jÞ, is defined as the ratio to the true weight of the absolute difference between the predicted and true weights. The average error over all nonzero edge weights is the prediction error Enz. In contrast, the absolute error Ez is used for non-existent edges or zero edge weights.23 The data ratio Rm is defined as the fraction of the number of measurements to the total number of terms to be predicted. To be more representative, for each set of parameters, 30 trajectories are generated for randomly set initial values. For each trajectory, a reference point is randomly chosen and the identification procedure is carried out. Then, TPR, TNR, Enz, and Ez are averaged over the 30 realizations. A. A network composed of the Lorenz oscillators

First, consider a network consisting of 10 bidirectionally coupled Lorenz oscillators, with the underlying topological structure being shown in Fig. 1. The coupling strength r is supposed to be 1. The density of this network is already higher than 10% as required by the method of Napoletani and Sauer.35

Chaos 25, 043102 (2015)

The traditional Lorenz system is37 8 > x_ ¼ 10ð y  xÞ; > > < y_ ¼ 28x  y  xz; > 8 > > : z_ ¼ xy  z: 3 Assume that the interaction between sub-variables only occurs between components x and z. That is, the inner coupling matrix is supposed to be 2 3 0 0 0 6 7 A ¼ 4 0 0 0 5: 1

0

0

In this example, the second-order Taylor expansion is used. If the estimated cij falls into the e neighbourhood of 1, that is, cij 2 ½1  e; 1 þ e, it is regarded that there exists an edge from node j to node i. Similarly, when cij 2 ½e; þe, the edge is considered to be non-existent. In this example, e is set to be 0.1. A reference point xi is randomly selected for node i ði ¼ 1; 2; …; N), and observation data are collected in its spherical neighborhood with a radius of 0.2. Then, 30 trajectories are generated from the network with 30 randomly set initial values. The identification procedure is carried out for 30 rounds and the corresponding average indexes such as the TPR, TNR, relative error, and absolute error are calculated. Identification results are displayed in Figs. 2 and 3, with standard error bars being included. From Fig. 2 it can be seen that the TPR gradually increases with the increase of data ratio Rm. When Rm reaches 80%, the TPR arrives at 100%. Meanwhile, the relative error of identification decreases with the increase of Rm and approaches 0 when the data ratio reaches 80%. From Fig. 3 it is seen that the TNR is kept at a high level. It also increases with increasing Rm. When the data ratio is increased to 80%, the TNR reaches 100%. The absolute error of identification is kept at a low level and declines with the increase of data ratio. It approaches 0 when Rm reaches 80%. Actually, since the second-order Taylor expansion is employed and node dynamics of the Lorenz system is determined by power series up to order 2, our method can also infer the node dynamical parameters. If the first-order Taylor expansion is used in this example, the neighborhood radius for data collection will affect the performance of topology identification. Further details are provided in Subsection III C. B. Networks composed of Duffing oscillators

FIG. 1. A network consisting of 10 coupled Lorenz oscillators.

To demonstrate the applicability of our method to a wide range of complex networks with unknown node dynamics, the Duffing oscillator is used as node dynamics. Since the system contains a trigonometric function as well as polynomials up to order of 3, the traditional compressive sensing based on polynomial basis cannot work well for reconstruction.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-5

Li et al.

Chaos 25, 043102 (2015)

FIG. 2. The true positive rate and relative error Enz versus the data ratio Rm for the Lorenz oscillator network, where e ¼ 0:1 and data collection neighborhood radius is 0.2.

The dynamical equation of the Duffing system is  x_ 1 ¼ x2 ; x_ 2 ¼ p2 x1  p3 x31  p1 x2 þ q cos xt:

(10)

When p1 ¼ 0:4;p2 ¼ 1:1;p3 ¼ 2; q ¼ 0:6; x ¼ 1:8, the system exhibits chaotic behavior, as shown in Fig. 4. The inner coupling matrix A is supposed to be an identity matrix, and the network structure is shown in Fig. 5. The coupling strength r ¼ 0:1. 1. Without system noise

First, consider a Duffing oscillator network without dynamical noise. In the following results, different orders of Taylor expansion are used. It is worth noting that the firstorder Taylor expansion illustrates the performance of the method proposed in Ref. 35. Time series of 30 trajectories are generated from network model (1) for randomly set initial values using the Runge-Kutta method with an equal time step of 0.01. Reference points are randomly chosen and data used for recovering the topology are sampled in the corresponding circular neighborhood of each reference point with a radius of 0.005. The first-order Taylor expansion is used for dynamical approximation. When the inferred cij 2 r½1  e; 1 þ e, the edge from node j to node i is regarded to be existent. If cij 2 r½e; þe, the edge from node j to node i is deemed as nonexistent. Here, e ¼ 0:1. The identification performance is shown in Fig. 6. It is seen from panel (1) that by using compressive sensing, one only needs 40% data to make the true positive rate to reach 74%. Furthermore, the rate roughly increases with increasing Rm. When Rm reaches 70%, the TPR is approximately 90%.

From panel (2), one can see that the true negative rate is larger than 90%, even when the data ratio is as low as 30%. Panel (3) displays the relative error Enz for all existent edges, which decreases with the increase of Rm. It already declines below 2% when the data ratio reaches 40%. Panel (4) shows the absolute error Ez of all non-existent edges, which is kept below 0.001. On the whole, identification performance of the Duffing oscillator network gets better with increasing number of observations. There is an unexpected slight decline in the TPR and TNR when Rm is 90%. Moreover, the TPR and TNR cannot reach 100% as expected, even if Rm is 1. These might be due to the inaccuracy caused by linear approximation and perturbation by the trigonometric function cos xt. To obtain more accurate topology identification, the second-order Taylor expansion is used to approximate node dynamics. The identification results are shown in Fig. 7, where parameters are the same as those for the first-order Taylor expansion. From panel (1) of Fig. 7, it is seen that the TPR increases smoothly with increasing data ratio Rm. It rises to 80% when the data ratio is 40% and arrives at 100% when the data ratio is 80%. Panel (2) shows that the TNR is kept above 90% even with a very small Rm and rises with increasing data ratio until it arrives at 100% when Rm is 80%. It can be seen from panel (3) that the relative error for existent edges decreases with the increase of Rm. It is smaller than 1% when the data ratio reaches 50%. Panel (4) of Fig. 7 shows the absolute error for non-existent edges, which is kept below 104 . Compare Fig. 6 with Fig. 7, one obtains that the identification performance based on the second-order Taylor expansion is better than the first-order one. This might be due to

FIG. 3. The true negative rate and absolute error Ez varying with the data ratio Rm for the Lorenz oscillator network, where e ¼ 0:1 and data collection neighborhood radius is 0.2.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-6

Li et al.

Chaos 25, 043102 (2015)

FIG. 4. The Duffing chaotic attractor.

the more accurate dynamical approximation of the secondorder expansion. Furthermore, the second-order Taylor expansion renders a more sparse unknown vector X than the first-order expansion. However, if one does not require a very precise topology identification result in practical application, the first-order Taylor expansion is enough since 60% data ratio already renders a TPR over 85% and a TNR over 90%, although there is a cubic term and a trigonometric term in the node dynamics. That is to say, the first-order Taylor expansion together with available local data collections can also work well for topology identification of complex networks composed of completely unknown node dynamics. Furthermore, we also employ Akaike information criterion (AIC)38 to measure the relative quality of the first-order and second-order Taylor expansion models. For the data used for Figs. 6 and 7, the AIC values are calculated and it is found that the AIC value of second-order Taylor expansion

FIG. 5. A network consisting of 20 coupled Duffing oscillators.

model is smaller than that of the first-order model, which means that the former model outperforms the latter one. This is consistent with the results shown in Figs. 6 and 7. 2. With system noise

In the real world, noise is omnipresent in nature and man-made systems due to various uncertainties such as stochastic forces on physical systems and noisy environments.39 In this subsection, topology identification of unknown networks with system noise will be probed into based on the proposed method. The dynamics of a network in the presence of system noise can be described by the following stochastic differential equations:

FIG. 6. Identification performance varying with the data ratio Rm for the 20-node Duffing network based on the first-order Taylor expansion and compressive sensing, where e ¼ 0:1 and data collection neighborhood radius is 0.005.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-7

Li et al.

Chaos 25, 043102 (2015)

FIG. 7. Identification performance versus the data ratio Rm for the 20-node Duffing network based on the secondorder Taylor expansion and compressive sensing, where e ¼ 0:1 and data collection neighborhood radius is 0.005.

x_ i ¼ f i ðt; xi Þ þ r

N X

_ i ðtÞ; cij Axj þ gi ðt; xi ÞB

(11)

j¼1

where gi ðt; xi Þ : Rþ  Rn ! Rnm indicates the noise intensity function related to dynamics of the i-th node, and Bi is an m-dimensional Brownian motion defined on the probability space ðX; F ; fF t gt0 ; PÞ with EfBi ðtÞg ¼ 0; EfB2i ðtÞg ¼ 1; EfBi ðsÞBi ðtÞg ¼ 0; s 6¼ t: A typical five-node network as shown in Fig. 8 is used as the testing network, with the coupling strength being 0.1. The non-autonomous Duffing system is still taken as node dynamics. The data collection radius is 0.005 and the threshold e ¼ 0:1. The impact of system noise on topology identification is displayed in Fig. 9. As shown in the four panels, when the

FIG. 8. A typical testing network consisting of 5 nodes.

noise intensity gradually changes from 0 to 0.012, the TPR and TNR fluctuate in a certain range but keep in a stable level. For the first-order Taylor expansion, the TPR and TNR are, respectively, maintained at about 60% and 73%, while for the third-order Taylor expansion, they are at about 87% and 91%. When the noise intensity is increased to 0.013, the identification rate encounters a sharp drop and the errors rise. However, for the third-order Taylor expansion, when the system noise intensity is no larger than 0.012, the TPR and TNR are still in a high level as over 80% in the average. Therefore, our method is robust against weak stochastic perturbations, which might be due to approximation of high-order Taylor expansion and the optimization nature of compressive sensing. C. Impact of data collection radius

It is obvious that the neighborhood radius of data collection will affect identification performance, since Taylor series of finite terms cannot accurately represent the original dynamics. What is the relationship between this neighborhood radius and identification performance? Is a smaller data collection neighborhood radius better than a larger one? In this subsection, more numerical simulations are carried out in order to assess the impact of data collection radii on identification performance. The 10-node Lorenz oscillator network as shown in Fig. 1 and the 20-node Duffing network as shown in Fig. 5 are used to illustrate the impact of neighborhood radius on identification performance. The TPR and TNR based on the first- and second-order Taylor expansion are shown in Figs. 10 and 11, where the data ratio is 80% and e ¼ 0:1. It is clearly observed that the smaller the neighborhood radius is, the better the identification performance is.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-8

Li et al.

Chaos 25, 043102 (2015)

FIG. 9. Identification performance versus system noise for the 5-node Duffing network based on Taylor expansion and compressive sensing, where e ¼ 0:1; Rm ¼ 70% and data collection neighborhood radius is 0.005.

FIG. 10. Identification performance varying with neighborhood radius of data collection for the 10-node Lorenz network based on the first- and secondorder Taylor expansion, where Rm ¼ 80% and e ¼ 0:1.

When the data collection radius is small enough, for the sparse Duffing network, both the first-order and second-order Taylor expansion yield accurate identification, with 100% TPR and 100% TNR. This is consistent with the nature of Taylor expansion, that is, the polynomials can more precisely approximate the original node dynamics in a smaller neighborhood. Another reason for the accurate identification is that the unknown vector is sparse enough. For the Lorenz

network, when the data collection radius is small enough, the second-order Taylor expansion can yield an accurate identification result, while the first-order Taylor expansion cannot. When the neighborhood radius is getting larger, the identification performance becomes worse. Basically, both the TPR and TNR decline with the increase of data sampling radius. From Figs. 10 and 11, one can also see that identification performance based on the second-order Taylor expansion and

FIG. 11. Identification performance varying with neighborhood radius of data collection for the 20-node Duffing network based on the first- and second-order Taylor expansion, where Rm ¼ 80% and e ¼ 0:1.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

043102-9

Li et al.

compressive sensing outperforms that of the first-order one. This is due to the fact that higher-order Taylor expansion can more accurately approximate node dynamics. For the Lorenz network, the identification results for the two cases have marked difference. The TPR and TNR based on the secondorder Taylor expansion are significantly larger than that of the first-order one. However, the difference for the Duffing network is not so significant as that for the Lorenz network. This might be due to the reason that the second-order Taylor expansion can exactly represent the Lorenz system, while the first-order cannot. However, for the Duffing system, Taylor expansion of any finite order cannot accurately replace the node dynamics due to the trigonometric term. In a word, Figs. 10 and 11 indicate that higher-order Taylor approximation and a smaller neighborhood for data collection are more favorable for recovering topologies. IV. CONCLUSIONS

A new technique for topology identification of complex networks has been proposed. Based on Taylor expansion and the traditional compressive sensing algorithm, topological structures of complex networks with completely unknown node dynamics can be recovered with a relatively small amount of observations. Numerical simulations have been presented to illustrate the effectiveness of the method. Furthermore, this method works well for complex networks contaminated with weak system noise. Additionally, identification performance varying with different parameters has been probed into. It has been found that a small neighborhood radius for data collection facilitates topology identification. Furthermore, when observations are only available in a large neighborhood, then high-order Taylor expansion is needed to render satisfying results. Nevertheless, when observations are available in a sufficiently small neighborhood, the simple first-order Taylor expansion together with compressive sensing is sufficient to yield satisfying identification results. ACKNOWLEDGMENTS

This work was supported by the National Natural Science and Technology Major Project of China (Grant No. 2014ZX1004-001-014) and the National Natural Science Foundation of China (Grant Nos. 61174028, 11172215, and 61203159).

Chaos 25, 043102 (2015) 1

P. DeLellis, M. D. Bernardo, and F. Garofalo, IEEE Trans. Circuits Syst. I 60, 3033 (2013). J. Zhou, X. Wu, W. Yu, M. Small, and J. Lu, Chaos 18, 043111 (2008). 3 X. Wu, W. Zheng, and J. Zhou, Chaos 19, 013109 (2009). 4 W. Yu, G. Chen, J. L€ u, and J. Kurths, SIAM J. Control Optim. 51, 1395 (2013). 5 S. Frenzel and B. Pompe, Phys. Rev. Lett. 99, 204101 (2007). 6 B. Pompe and J. Runge, Phys. Rev. E 83, 051122 (2011). 7 S. Hempel, A. Koseska, J. Kurths, and Z. Nikoloski, Phys. Rev. Lett. 107, 054101 (2011). 8 D. Yu, M. Righero, and L. Kocarev, Phys. Rev. Lett. 97, 188701 (2006). 9 J. Zhou and J. Lu, Physica A 386, 481 (2007). 10 X. Wu, Physica A 387, 997 (2008). 11 H. Liu, J. Lu, J. L€ u, and D. J. Hill, Automatica 45, 1799 (2009). 12 J. Zhao, L. Qin, J. Lu, and Z. P. Jiang, Chaos 20, 023119 (2010). 13 R. Jansen, H. Yu, D. Greenbaum, Y. Kluger, N. J. Krogan, S. Chung, A. Emili, M. Snyder, J. F. Greenblatt, and M. Gerstein, Science 302, 449 (2003). 14 K. Sachs, O. Perez, D. Pe’er, D. A. Lauffenburger, and G. P. Nolan, Science 308, 523 (2005). 15 N. Marwan, M. C. Romano, M. Thiel, and J. Kurths, Phys. Rep. 438, 237 (2007). 16 M. C. Romano, M. Thiel, J. Kurths, and C. Grebogi, Phys. Rev. E 76, 036211 (2007). 17 W.-X. Wang, Q. Chen, L. Huang, Y. C. Lai, and M. A. F. Harrison, Phys. Rev. E 80, 016116 (2009). 18 J. Ren, W.-X. Wang, B. Li, and Y. C. Lai, Phys. Rev. Lett. 104, 058701 (2010). 19 W.-X. Wang, J. Ren, Y. C. Lai, and B. Li, Chaos 22, 033131 (2012). 20 X. Wu, W. Wang, and W. Zheng, Phys. Rev. E 86, 046106 (2012). 21 X. Wu, C. Zhou, G. Chen, and J. Lu, Chaos 21, 043129 (2011). 22 T. He, X. Lu, X. Wu, J. Lu, and W. Zheng, Physica A 392, 1038 (2013). 23 W.-X. Wang, R. Yang, Y.-C. Lai, V. Kovanis, and C. Grebogi, Phys. Rev. Lett. 106, 154101 (2011). 24 W.-X. Wang, R. Yang, Y.-C. Lai, V. Kovanis, and M. A. F. Harrison, Europhys. Lett. 94, 48006 (2011). 25 W.-X. Wang, Y.-C. Lai, C. Grebogi, and J. Ye, Phys. Rev. X 1, 021021 (2011). 26 R.-Q. Su, W.-X. Wang, and Y.-C. Lai, Phys. Rev. E 85, 065201 (2012). 27 R.-Q. Su, X. Ni, W.-X. Wang, and Y.-C. Lai, Phys. Rev. E 85, 056220 (2012). 28 Z. Shen, W.-X. Wang, Y. Fan, Z. Di, and Y.-C. Lai, Nat. Commun. 5, 4323 (2014). 29 D. Donoho, IEEE Trans. Inf. Theory 52, 1289 (2006). 30 E. Candes, J. Romberg, and T. Tao, IEEE Trans. Inf. Theory 52, 489 (2006). 31 E. Candes and T. Tao, IEEE Trans. Inf. Theory 52, 5406 (2006). 32 R. G. Baraniuk, IEEE Signal Process. Mag. 24, 118 (2007). 33 E. Candes and M. Wakin, IEEE Signal Process. Mag. 25, 21 (2008). 34 Q. Fan, Y. Jiao, and X. Lu, IEEE Trans. Signal Process. 62, 6276 (2014). 35 D. Napoletani and T. D. Sauer, Phys. Rev. E 77, 026103 (2008). 36 M. Fornasier and H. Rauhut, Handbook of Mathematical Methods in Imaging (Springer, 2011), pp. 187–228. 37 E. N. Lorenz, J. Atmos. Sci. 20, 130 (1963). 38 See http://en.wikipedia.org/wiki/Akaike_information_criterion for Akaike information criterion (AIC), which is a measure of the relative quality of a statistical model for a given set of data and provides a means for model selection. 39 P. S. Swain, Chaos 16, 026101 (2006). 2

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Wed, 01 Jul 2015 10:23:16

Recovering network topologies via Taylor expansion and compressive sensing.

Gaining knowledge of the intrinsic topology of a complex dynamical network is the precondition to understand its evolutionary mechanisms and to contro...
2MB Sizes 2 Downloads 8 Views